TCGA中的RNA表达数据整理之Ensembl【数据挖掘专题】

TCGA中的RNA表达数据整理之Ensembl【数据挖掘专题】推文来源：生信控
the RNA-Seq count matrix had 60,483
Ensembl IDs. We used the GTF file
(Homo_sapiens.GRCh38.86.chr) to
obtain the gene symbol for each
Ensembl ID, and 57,288 genes were kept
for further study.
这里用到的注释数据库是Ensembl。

而TCGA官方RNA分析流程：
/Data/BioinformaticsPipelines/ExpressionmRNA_Pipeline/
是使用GENCODE v22作为基因组注释文件：
/download-gdc-reference-files
gencode.v22.annotation.gtf.gz
除了数据挖掘专题| ID转换的N种方式，网上随处可见丰(luan)富(qi)多(ba)彩(zao)的用法，但主要还是围绕 Ensembl 和 GENCODE ，将数据库中的GTF注释文件下载到本地，并在R中完成处理。

不过，数据库会有不同的版本、同一版本中又提供了不同的下载文件，所以...该如何选择合适的数据库或许是个值得关注的问题！
本期开始，我们将陆续解决其中的一些问题。

首先，尝试重复这个：
the RNA-Seq count matrix had 60,483
Ensembl IDs. We used the GTF file
(Homo_sapiens.GRCh38.86.chr) to
obtain the gene symbol for each
Ensembl ID, and 57,288 genes were kept
for further study.
由上示文件名可见，用到的Ensembl数据库版本为v86，注释数
据（GTF文件）下载链接如下：
ftp:///pub/release-86/gtf/homo_sapiens/
ftp:///pub/release-
86/gtf/homo_sapiens/Homo_sapiens.GRCh38.86.chr.gtf.gz 常规对于这种压缩文件，首先会手动下载到本地、然后解压读取。

不过，也可以参考 R语言| 读取压缩文件，其中提供了一种无需下载，可直接读取网页压缩文件的方法！
而针对于GTF文件的读取，推荐使用rtracklayer 包中的 import 函数，这个函数可以直接读取网页上的压缩的GTF文件，其内部实现是首先将网页上的文件下载到临时目录下，详见 R语言| 看到数据被下载但找不到？
1.p_load(rtracklayer)
2.AnnoData = import('ftp:///pub/release-86/gtf/homo_sapiens/Homo_sapiens.GRCh38.86.chr.gtf.gz') 不过，如果网速不好且本地不缺存储空间的话，还是建议先下载到本地然后读取！
报错的看 R语言 | 如果你再问我怎么安装R包
现在可以对于从数据挖掘专题 | TCGA-RNA数据下载全攻略中下载得到的表达数据（dataAssy）进行基因名的转换了：
1.AnnoData = AnnoData[which(AnnoData$type == 'gene'),]
2.tmp = data.frame(Ensembl_ID = AnnoData$gene_id, Symbol = AnnoData$gene_name, Biotype = AnnoData$gene_biotype)
3.tmp_1 = tmp[tmp$Ensembl_ID %in% rownames(dataAssy),]
4.nrow(tmp_1) # 57288
没错，最终如文中所述，得到57288个基因！
1.head(tmp_1)
那么，这意味着可以使用得到的Symbol作为dataAssy的行名了吗？
1.tmp_1[which(tmp_1$Symbol == 'HAR1A'),]
Ensembl_ID Symbol Biotype
52278 ENSG00000274915
HAR1A misc_RNA
52279 ENSG00000225978
HAR1A lincRNA
可以看到，同一个Symbol对应了两个Ensembl_ID，进一步的：
s(which(table(tmp_1$Symbol) > 1))
2.length(which(table(tmp_1$Symbol) > 1)) # 218
共有218个基因存在同一个Symbol对应多个Ensembl_ID的情况！（没见过有网上教程提到这一点...）
所以，为什么会存在这种情况：
There is no one-to-one mapping of
gene ids from one database to the other.
Ensembl (who maintain Ensembl ENSG
IDs), ncbi (who maintain EntrezGene IDs
and RefSeq transcript ids) and HUGO (who
maintain gene symbols/names) all have
different ideas about what a gene is,
which part of the genome belongs to
which gene etc.
应该如何处理这种情况？
答案是没有定论，不过主流的建议是基于Ensembl ID进行后续分析！比如在获取差异基因、预后模型之后再对其进行Symbol的转换。

但是无法避免的是，可能仍然会存在重复...
参考：
/p/352492/
/p/119540/
/p/90185/
/p/105430/
/questions/2959/duplicate-genes-with-rsem-counts-which-one-to-choos
你也可以将对应同一Symbol的多个Ensembl表达求和，方法可以参考 R语言 | 统计不同组患者的生存时间
或者干脆将这些基因排除：
1.tmp_clean = tmp_1[!tmp_1$Symbol %in% names(which(table(tmp_1$Symbol) > 1)),]
2.nrow(tmp_clean)
[1] 55414
可见，最终保留55414个基因！。