利用Perl或者R语言对TCGA肿瘤数据lincRNA提取

Posted 2021-04-17 医知圈

医知圈

        众所周知lncRNA属于RNA中的非编码RNA，在转录调控中扮演者重要角色。并且现在听说lncRNA的研究也很火热，使用TCGA的数据对lncRNA的研究也是非常常见的需求。而我们如果想对TCGA的lncRNA进行定量则必须从TCGA RNA-Seq的表达量数据中提取出lncRNA的那部分数据。

我之前也没做过lncRNA相关的研究，所以就从网上搜了下，整理了下思路

首先下载TCGA中的RNA-Seq数据，因为新版的TCGA数据库是用Ensembl id作为基因的命名的，所以这些RNA-Seq基因表达量文件中必然也包含了lncRNA的表达量数据

接着要获取TCGA RNA表达量的Ensembl id中属于lncRNA的id

最后从总的RNA表达量文件中提取lncRNA的表达量数据

1. 第一步已经在前面里做了，是以乳腺癌的TCGA数据为例

2. 第二步获取ensembl id与lncRNA的对应关系

   这问题的解决思路第一个想起来应该去GENCODE数据库中寻找，因为逛过ENCODE人类相关数据库的话，会发现有个Long non-coding RNA gene annotation选项，果断下载了其gtf文件，GENCODE对这个文件的表述：evidence-based annotation of the human genome (GRCh38), version 27 (Ensembl 90) – long non-coding RNAs，从中也能看出该文件的版本跟TCGA使用的版本也保持一致，因此只要对该gtf文件进行处理就能获得属于lncRNA的ensembl id

   粗略看几行gtf格式

   chr1    HAVANA  gene    29554   31109   .   +   .   gene_id "ENSG00000243485.5"; gene_type "lincRNA"; gene_name "MIR1302-2HG"; level 2; tag "ncRNA_host"; havana_gene "OTTHUMG00000000959.2";chr1    HAVANA  transcript  29554   31097   .   +   .   gene_id "ENSG00000243485.5"; transcript_id "ENST00000473358.1"; gene_type "lincRNA"; gene_name "MIR1302-2HG"; transcript_type "lincRNA"; transcript_name "MIR1302-2HG-202"; level 2; transcript_support_level "5"; tag "not_best_in_genome_evidence"; tag "dotter_confirmed"; tag "basic"; havana_gene "OTTHUMG00000000959.2"; havana_transcript "OTTHUMT00000002840.1";chr1    HAVANA  exon    29554   30039   .   +   .   gene_id "ENSG00000243485.5"; transcript_id "ENST00000473358.1"; gene_type "lincRNA"; gene_name "MIR1302-2HG"; transcript_type "lincRNA"; transcript_name "MIR1302-2HG-202"; exon_number 1; exon_id "ENSE00001947070.1"; level 2; transcript_support_level "5"; tag "not_best_in_genome_evidence"; tag "dotter_confirmed"; tag "basic"; havana_gene "OTTHUMG00000000959.2"; havana_transcript "OTTHUMT00000002840.1";chr1    HAVANA  exon    30564   30667   .   +   .   gene_id "ENSG00000243485.5"; transcript_id "ENST00000473358.1"; gene_type "lincRNA"; gene_name "MIR1302-2HG"; transcript_type "lincRNA"; transcript_name "MIR1302-2HG-202"; exon_number 2; exon_id "ENSE00001922571.1"; level 2; transcript_support_level "5"; tag "not_best_in_genome_evidence"; tag "dotter_confirmed"; tag "basic"; havana_gene "OTTHUMG00000000959.2"; havana_transcript "OTTHUMT00000002840.1";chr1    HAVANA  exon    30976   31097   .   +   .   gene_id "ENSG00000243485.5"; transcript_id "ENST00000473358.1"; gene_type "lincRNA"; gene_name "MIR1302-2HG"; transcript_type "lincRNA"; transcript_name "MIR1302-2HG-202"; exon_number 3; exon_id "ENSE00001827679.1"; level 2; transcript_support_level "5"; tag "not_best_in_genome_evidence"; tag "dotter_confirmed"; tag "basic"; havana_gene "OTTHUMG00000000959.2"; havana_transcript "OTTHUMT00000002840.1";

   因为TCGA是基于基因水平进行定量分析的，那么我需要将第3列是gene的行中的ensembl id，gene name提取出来，ensembl id需要去掉版本号；而至于为什么提gene name，因为我想最后以gene name来表示各个gene。

3. lncRNA_ensembl2name.list文件总共有15778行，其中12行是重复的数据，那么去重后是15766行。

4. 如果仔细研究下，会发现lncRNA应该有其特定的几个gene_type，从上述的gtf文件就可看出，那么哪些gene_type属于lncRNA呢，首先我们从GENCODE的gtf文件中先提取下看看

lncRNA_ensembl2name.list文件总共有15778行，其中12行是重复的数据，那么去重后是15766行

如果仔细研究下，会发现lncRNA应该有其特定的几个gene_type，从上述的gtf文件就可看出，那么哪些gene_type属于lncRNA呢，首先我们从GENCODE的gtf文件中先提取下看看

输出结果有3列，第一列是ensembl id，第二列是gene name，第三列是biotype，总共有14700行，然后跟GENCODE的结果比较下，后者包含了前者，说明结果是正常的，因为后者还多一个TEC的biotype。

最后我例行跟生信人工具再比较下，因为其也有一个ENSG_ID_LNC.txt文件，里面记录了ensembl id跟lncRNA的对应关系，当然也有biotype信息，先粗略看下biotype有几类，结果如下：

3prime_overlapping_ncRNA antisense bidirectional_promoter_lncRNA lincRNA macro_lncRNA processed_transcript sense_intronic sense_overlapping

跟我上面总结的结果比少了一个non_coding的biotype，那么再看一下相同的结果有多少，结果显示只有4000不到的对应关系是相同的，其实就是说至于4000不到的ensembl id对应的lncRNA的gene name是相同的，应该是lncRNA的命名的差别，只是还没查到其ensembl对应的lncRNA是什么来历；有些是lncRNA命名的版本号不同。

1. 最后检查了好久，才晓得原来生信人用的人类ENSEMBL版本是GRCh37.p13，而我用的是GRCh38.p10，TCGA用的版本也是GRCh38，具体可参看网址GDC TCGA的历史版本，其在2016年就使用GRCh38了，比如mRNA-Seq数据则是BAM files aligned to GRCh38 using STAR 2-pass strategy，所以还是纯手工的比较靠谱，软件毕竟是别人开发的。。。

2. 第三步，只要基于ensembl id和lncRNA id的对应关系，从TCGA的表达矩阵中属于lncRNA的那些基因的表达量提取出来即可，然后再按照需要再考虑要不要转化为ensembl的gene name。

   本模式代码小伍全已编写完毕。