KeepNotes blog

Stay hungry, Stay Foolish.

0%

TCGA ID转化的小插曲

Posted on Edited on In ,
Symbols count in article: 4k Reading time ≈ 4 mins.

在下载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.egENSEMBL2EG和org.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转载请注明出处

Related Posts