Bioconductor開發(fā)的物種注釋包系列集合了一個(gè)物種不同來(lái)源的注釋信息，能夠根據(jù)基因ID對(duì)其進(jìn)行多種來(lái)源的注釋，比如說(shuō)基因的別名，基因的功能等。
我之前也寫過(guò)一篇文章用Bioconductor對(duì)基因組注釋介紹如何使用AnnotationHub下載注釋數(shù)據(jù)庫(kù), 使用select(), mapIds等函數(shù)進(jìn)行注釋操作。我自己寫一個(gè)流程也用到了它給基因ID, 如AT1G14185, 注釋別名和功能描述。 注釋結(jié)果中會(huì)出現(xiàn)一些基因無(wú)法被注釋, 比如說(shuō)下面這些情況, 我一直認(rèn)為只是這些基因沒(méi)有得到比較好的研究, 即便這些基因能夠在TAIR搜到Araport的注釋, 我也認(rèn)為那些注釋只是同源注釋沒(méi)有多大意義。
AT1G13970 NA NA
AT1G14120 NA NA
AT1G14240 NA NA
AT1G14600 NA NA
一開始得到的結(jié)果里沒(méi)有多少個(gè)基因，所以缺少幾個(gè)注釋，通過(guò)手工去查找也行，但是目前差異表達(dá)分析動(dòng)不動(dòng)就給別人500多個(gè)基因，于是就有幾十個(gè)甚至上百個(gè)未注釋的基因，所以我想著要不自己更新擬南芥的物種包。
library(AnnotationHub)
ah <- AnnotationHub()
org <- ah[["AH57965"]]
org
#...
# TAIRGENEURL: ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_functional_descriptions
#...
通過(guò)上面的代碼，我找到了基因功能描述的數(shù)據(jù)庫(kù)來(lái)源文件，我下載了這個(gè)文件，并且拿用AnnotationHub注釋不到的功能的一個(gè)基因，"AT1G14185"，進(jìn)行測(cè)試
mapIds(org, "AT1G14185", "GENENAME", "TAIR")
這下就非常有趣了，在原始文件中能搜索到的基因用Bioconductor的物種注釋包時(shí)卻沒(méi)有注釋信息！為了搞清楚這個(gè)原因，我花了快半個(gè)下午的時(shí)間去折騰，終于被我找到了原因。 我分別用一個(gè)能被org.At.tair.db注釋和一個(gè)不能被org.At.tair.db的注釋去搜索原始文本.
# 無(wú)注釋
1 AT1G14185.1
2 protein_coding
3 Glucose-methanol-choline (GMC) oxidoreductase family protein
4
5 Glucose-methanol-choline (GMC) oxidoreductase family protein; FUNCTIONS IN: ...
# 有注釋
1 AT1G19610.1
2 protein_coding
3 Arabidopsis defensin-like protein
4 Predicted to encode a PR (pathogenesis-related) protein. ...
5 PDF1.4; FUNCTIONS IN: molecular_function unknown;
簡(jiǎn)單的比較之后，你差不多就知道了org.At.tair.db的在功能描述這一部分其實(shí)只用第一列和第四列(為了方便展示我轉(zhuǎn)置了原始數(shù)據(jù))。這就是非常讓人意外了，為啥它不用第一列和第三列呢？ 我于是又去看了其他幾個(gè)基因，就差不多明白了，原始的文本特別的混亂，你除了能保證第一列和第二列有信息外，其他列你根本無(wú)法保證，因此最好的策略以第一列作為檢索的關(guān)鍵字，其他列合并成一列才行，然而作者沒(méi)有那么細(xì)致。
于是我就放棄了用org.At.tair.db注釋基因功能描述和基因別名了，還是自己寫一個(gè)Python腳本進(jìn)行注釋吧。
下面這個(gè)腳本只適用于bed格式的輸入，且第四列為轉(zhuǎn)錄本ID，另外兩個(gè)輸入文件分別為"gene_aliases_20140331.txt"和"TAIR10_functional_descriptions_20140331.txt"， 用法為
python bed_anno.py to_anno.bed gene_aliases_20140331.txt TAIR10_functional_descriptions_20140331.txt > anno.xls
import sys
from collections import defaultdict
bed_file = sys.argv[1]
alias_file = sys.argv[2]
func_file = sys.argv[3]
alias_dict = defaultdict(list)
func_dict = defaultdict(list)
# read alias file
for line in open(alias_file, 'r'):
items = line.strip().split('\t')
alias_dict[items[0]] = items[1:]
# read function description file
for line in open(func_file, 'r'):
items = line.strip().split('\t')
func_dict[items[0]] = items[1:]
# annotation and output
for line in open(bed_file, 'r'):
transcript_id = line.strip().split("\t")[3]
gene_id = transcript_id.split(".")[0]
gene_alias = alias_dict[gene_id] if len(alias_dict[gene_id]) > 0 else ['']
gene_func = func_dict[transcript_id] if len(func_dict[transcript_id]) > 0 else ['']
gene_anno = '{}\t{}\t{}'.format(line.strip(), gene_alias[0], '\t'.join(gene_func))
print(gene_anno)
