TCGA ID 转化的小插曲

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

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

如果对生信技能中ID转化没什么概念的话，可以先看看我们生信技能树中曾经的一道编程题练练手生信编程直播第8题-几个ID转换咯 http://www.biotrainee.com/thread-941-1-1.html

NCBI数据库脚本转化

名字有点绕，其实是依靠NCBI提供一些ID对应关系的文件，其中有一个gene2ensembl文件就提供了gene id和ensembl id的对应关系。当然我们是想将ensembl id转化为gene symbol，所以还需要gene_info文件，其提供了gene id和gene symbol的对应关系，这样思路清楚了：整理上述2个文件，获得ensembl id与gene symbol的对应关系文件。

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

接着我们要从这两个文件中分别提取人类（taxid：9606）的gene id，ensembl id，gene symbol等信息：

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

2. 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

1. #!/bin/perl -w

2. use strict;

3. my $gene2ensembl = shift @ARGV;

4. my $gene2symbol = shift @ARGV;

5. my %hash;

6. open my $fh1, $gene2ensembl or die;

7. while (<$fh1>) {

8.    chomp;

9.    my @array = split /\t/, $_;

10.    $hash{$array[0]} = $array[1];

11. }

12. close $fh1;

13. open my $fh2, $gene2symbol or die;

14. while (<$fh2>) {

15.    chomp;

16.    my @array = split /\t/, $_;

17.    if ($hash{$array[0]}) {

18.        print "$hash{$array[0]}\t$array[1]\n";

19.    }

20. }

21. 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转化，比如TCGA肯定是用org.Hs.eg.db包了，然后利用 org.Hs.egENSEMBL2EG和 org.Hs.egSYMBOL中的数据；从命名上应该很好理解，前者是ensembl id和gene id的对应关系，后者是gene id和gene symbol的对应关系。最后整理下获得跟上述一样的ensembl id和gene symbol的对应关系。

1. library(org.Hs.eg.db)

2. ensembl2gene <- toTable(org.Hs.egENSEMBL2EG)

3. gene2symbol <- toTable(org.Hs.egSYMBOL)

4. ensemble2symbol <- merge(ensembl2gene, gene2symbol, by = "gene_id")[2:3]

5. 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也只有一个gene 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/hgnccompleteset.txt

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

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

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

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

结果还真有9个，再把9个查了下ncbi和ensembl网站，都是存在，只是有几个NCBI的symbol跟emsembl的不一致，好尴尬。

那为何这种的ensembl id与symbol的对应关系相比前两种方法多这么多呢。后来我查了下，发现NCBI的gene2ensembl文件里的gene id和ensembl id的对应关系一般都是有对应RNA和protein的（官网解释：This file reports matches between NCBI and Ensembl annotation based on comparison of rna and protein features；而且R包中数据也是基于NCBI的那个data数据库，所以可以粗略前2种方法是没有假基因存在的。那我们也剔除掉假基因，再看看对应关系，可以看到假基因在标签在hgnccompleteset.txt文件的第4列，那么代码如下

1. 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

这个剔除掉假基因的hgncensembl2symbol2.list文件只有26034个对应关系了，但是跟方法1还是只有90%的相似度。好吧，还是以这次的为准。毕竟我还是挑了好几个（有差异的）查了下的，都跟ENSEMBL数据库中收录的保持一致，毕竟NCBI收录的跟ensembl收录的还是有点差别的。

HGNC Symbol与生信人工具的比较

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

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

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

Summary

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

如果是ID转换知识，推荐我们技能树以前的总结帖：