|
### RNA表达数据 在TCGA改版之前，从TCGA中下载并整理好的RNASeqV2数据，或者改版后从GDC Legacy Archive中下载的RNA数据，其格式如下： 第一列Symbol，共计20502个基因，其中包含了mRNA和lncRNA，基于文件中的Gene Symbol虽然可以提取lncRNA，但是数目较少（可能就几十或者几百个！） 而从GDC Data Portal上下载（即官网直接下载）的RNA表达数据（经重新比对定量，且基因名称全部换成了Ensembl ID），内容如下： 第一列Ensembl ID，共计60483个基因（接近GDC Legacy Archive上的3倍），其中也包含了mRNA和lncRNA，而lncRNA的数目在10000+，所以如果想基于TCGA进行lncRNA数据的分析，使用GDC Data Portal上的数据是个不错的选择，现在问题就变成了如何从一列Ensembl ID中识别lncRNA。 ### 注释数据 1、我们在之前的文章TCGA-RNA数据下载全攻略中已经讲到过，GDC Data Portal中RNA的注释数据库是GENCODE： https://docs.gdc./Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/ 而通过GENCODE数据库注释文件中的 biotype 可以得知每个基因是否属于lncRNA，而具体哪些注释算作lncRNA呢？ 先来几篇文献参考一下： A、6类：lincRNA、processed_transcript、sense_intronic、sense_overlapping、antisense、3prime_overlapping_ncrna： B、6类：lincRNA、antisense、sense_intronic、sense_overlapping、processed_transcript、processed_pseudogene： C、8类：non-coding、3-prime overlapping ncrna、antisense、bidirectional promoter lncRNA、lincRNA、macro lncRNA、sense intronic、sense overlapping 随机检索了几篇可见，不仅用的GENCODE数据库版本不一致(v17~v27基本上均有人使用)，针对lncRNA的提取标准也不一样(有的纳入的biotype种类多，有的少)... 2、GENCODE中也提供了整理好的lncRNA专属信息： ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_27/gencode.v27.long_noncoding_RNAs.gtf.gz 其中，收录了属于lncRNA的10种基因： 4、LNCipedia https:/// … 所以，提取lncRNA的方法五花八门，因为理论上来讲，只要有lncRNA注释信息的数据库都可以用于提取lncRNA数据！ ### GENCODE lncRNA 我们此处仅以上示第二种提取方法为例，讲解如何从TCGA中筛选得到lncRNA表达数据： gencode.v27.long_noncoding_RNAs.gtf文件格式如下： 比较关键的gene_id、gene_type、gene_name，对于这三列信息的提取，如果在命令行中使用awk可参考： https://www./p/77347/ 下面重点讲解基于R语言的解决方案： # 安装rtracklayer包 # https:///packages/release/bioc/html/rtracklayer.html source('http:///biocLite.R') biocLite('rtracklayer') # 使用 library(rtracklayer) AnnoData = import('gencode.v27.long_noncoding_RNAs.gtf') 将3列信息单独提出来： index = which(AnnoData$type == 'gene') Target = data.frame(Ensembl_ID = AnnoData$gene_id[index], Symbol = AnnoData$gene_name[index], Biotype = AnnoData$gene_type[index]) 共计15778个lncRNA，在与RNA数据取lncRNA交集之前需要先把Ensembl_ID中的版本号去除掉： Target$Ensembl_ID = gsub('\\..*', '', Target$Ensembl_ID) 最终格式如下： 然后与TCGA-RNA数据中的Ensembl ID（同样要去除版本号）取交集: 结合之前文章 GDC：我们不一样！中的表达数据 common = intersect(Target$Ensembl_ID, rownames(D_coad)) D_coad_Hiseq_exp = D_coad[common, intersect(D_coad_samples, LH_coad_samples)] 最终共得到15329个lncRNA在323个COAD(Hiseq)样本中的表达值！ 至此，关于如何从TCGA中下载并整理lncRNA数据已经讲解完毕，如果小伙伴们有什么问题、意见和建议请指教👇👇👇 |
|
|
