TCGA ID转化的小插曲

下载TCGA数据后,我们会发现新版的TCGA表达数据有6万个ENSG(ENSEMBL ID),听说旧版时是用gene symbol的,而且我也觉得有时gene symbol更能说明问题。如果用过Firehose下载的TCGA数据的话,会发现其是用gene id和gene symbol共同来表示表达数据的。如果是使用生信人的那款工具的话,其也是提供了ID转化功能,即把ENSG转化为gene symbol。

作为一个生信工作者,当然有自己的一套ID转化方法,所以我想把从官网下载方式下载下来的TCGA数据,以不同的方法将ENSG转化为gene symbol,然后跟Firehose和生信人工具进行比较下。。。

NCBI数据库脚本转化

名字有点绕,其实就是以我之前的一篇博文的方法(也是学习别人来的)NCBI/Ensembl ID的转换,我也比较习惯用这个来做ID转化。

NCBI提供一些基因相关文件,其中有一个gene2ensembl文件就提供了gene id和ensembl id的对应关系。当然我们是想将ensembl id转化为gene symbol,所以还需要一个gene_info文件,其提供了gene id和gene symbol的对应关系,这样思路的清楚了。

两个文件的下载地址:ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/

接着我们要从这两个脚本中分别提取人类(taxid:9606)的gene id,ensembl id,gene symbol等信息:

perl -alne 'print $_ if ($F[0] == 9606)' gene2ensembl |cut -f 2,3 >9606.gene2ensembl.list
perl -alne 'print $_ if ($F[0] == 9606)' gene_info |cut -f 2,3 >9606.gene2symbol.list

获得9606.gene2ensembl.list和9606.gene2symbol.list文件后,我们还需要将两个文件合并,整合成一个9606.ensembl2symbol.list文件,写一个简单的脚本1.pl

#!/bin/perl -w
use strict;

my $gene2ensembl = shift @ARGV;
my $gene2symbol = shift @ARGV;
my %hash;

open my $fh1, $gene2ensembl or die;
while (<$fh1>) {
    chomp;
    my @array = split /\t/, $_;
    $hash{$array[0]} = $array[1];
}
close $fh1;

open my $fh2, $gene2symbol or die;
while (<$fh2>) {
    chomp;
    my @array = split /\t/, $_;
    if ($hash{$array[0]}) {
        print "$hash{$array[0]}\t$array[1]\n";
    }
}
close $fh2;

运行脚本perl 1.pl 9606.gene2ensembl.list 9606.gene2symbol.list> 9606.ensembl2symbol.list,其中ensembl id与gene symbol有对应的ID数目有25385个,其中有3个gene symbol是重复的,应该是不同的几个ensembl id对应到同一个gene symbol的原因。

R包转换

这个方法就比较好理解了,就是利用注释R包中的数据进行ID转化,比如这个肯定是用org.Hs.eg.db包了,然后利用org.Hs.egENSEMBL2EGorg.Hs.egSYMBOL中的数据;从命名上应该很好理解,前者是将ensembl id和gene id的对应关系,后者是gene id和gene symbol的对应关系。最后获得跟上述一样的ensembl id和gene symbol的对应关系

library(org.Hs.eg.db)

ensembl2gene <- toTable(org.Hs.egENSEMBL2EG)
gene2symbol <- toTable(org.Hs.egSYMBOL)
ensemble2symbol <- merge(ensembl2gene, gene2symbol, by = "gene_id")[2:3]
write.table(ensemble2symbol, file = "ensembl2symbol.txt", sep = "\t", quote = F, row.names = F)

ensembl2symbol.txt文件中有28945个ID对应关系,这比第一种方法获得的结果还多3000多个。。。然后我粗略的检查了下,发现是有多个gene id对应到同一个ensembl id上的情况,然后我也有理由相信第一种方法也会有这种情况发生(但是检查了下,这种情况比较少,大约只有39个。。)。但是我在ENSEMBL官网查到一般一个ensemble id也只有一个HGNC Symbol。

HGNC Symbol

既然从ensembl id转换到gene id,再转化到gene symbol会有这么曲折,那我干脆直接从ensembl id转化到HGNC symbol吧(现在很多官方symbol都用HGNC symbol了?),第一个想到是用EMSEMBL的biomart,然后想了想,还是从HGNC官网上下载信息然后用代码来提取symbol与ensembl id的对应关系吧

下载地址:ftp://ftp.ebi.ac.uk/pub/databases/genenames/new/tsv/hgnc_complete_set.txt

然后看了下第2列是symbol,第20列是ensembl id,所以用perl提取下

perl -F'\t' -alne 'print "$F[19]\t$F[1]" if ($F[19])' hgnc_complete_set.txt >hgnc_ensembl2symbol.list

总共有37358个对应关系,然后再查看下有无同一个ensembl id对应到多个symbol上的

perl -alne '$hash{$F[0]}++;print $_ if ($hash{$F[0]}>1)' hgnc_ensembl2symbol.list

结果还真有9个,再把9个查了下ncbi和ensembl网站,都是存在,只是有几个NCBI的symbol跟emsembl id没对应上,好尴尬。。。

那为何这种的ensembl id与symbol相比前两种多这么多呢,后来我查了下,发现gene2ensembl文件里的gene id和ensembl id的对应关系一般都是有对应RNA和protein的,而且R包中数据也是基于NCBI那个data数据库的,所以可以粗略判断是没有假基因存在的。。。那如果我们在这方法也剔除掉假基因,再看看对应关系呢。可以看到假基因在标签在hgnc_complete_set.txt文件的第4列,那么代码如下

 perl -F'\t' -alne 'next if ($F[3] eq "pseudogene"); print "$F[19]\t$F[1]" if ($F[19])' hgnc_complete_set.txt >hgnc_ensembl2symbol_2.list

这个剔除掉假基因的hgnc_ensembl2symbol_2.list文件只有26034个对应关系了,但是跟方法1还是只有90%的相似度。。。好吧,还是以这次的为准。。(毕竟我还是挑了好几个查了下的,而且NCBI收录的毕竟还是跟ensembl有点差别的)

HGNC Symbol与生信人工具的比较

在我查看生信人这款工具的时候,我猜想其ID转化功能必然有一个ID对应关系的文件,果然在软件的首目录下就有个ENSG_ID.txt文件,打开就能发现其就是ensembl id与symbol的对应关系,总共有25527条信息。

这跟我通过HGNC Symbol找到的对应关系也有初入,简单做个venn图比较下,ENSG_ID.txt文件中有2414条信息没有在hgnc_ensembl2symbol_2.list文件中找到。然后我找了几个ensembl id在ENSEMBL官网上搜下了,发现没找到的那2414个ID中有些是假基因。。。

所以我拿未剔除假基因的对应关系hgnc_ensembl2symbol.list文件与ENSG_ID.txt文件再做个venn图 看看,结果ENSG_ID.txt文件中有2010条信息没在hgnc_ensembl2symbol.list文件中找到。然后也挑了几个在ENSEMBL官网搜了下,这次结果还是可以接受,hgnc_ensembl2symbol.list结果跟官网还是保持一致的,而生信人工具一些symbol是HGNC symbol的别名。。当然也不是说别名就不能用啦。。。

Summary

这篇文章主要还是练习了下ID转化的方法,至于TCGA到底用哪个ID转化方法的还是看情况了吧,如果只是转化为gene id,而且也不考虑假基因什么的不常见的情况,那么方法1和方法2都是可以的;如果需要转化为symbol,我觉得如果是我自己来做的话,一般会选择方法3。

本文出自于http://www.bioinfo-scrounger.com转载请注明出处