Bioconductor開發(fā)的物種注釋包系列集合了一個物種不同來源的注釋信息，能夠根據(jù)基因ID對其進(jìn)行多種來源的注釋蹋宦，比如說基因的別名修档，基因的功能等。
我之前也寫過一篇文章用Bioconductor對基因組注釋介紹如何使用AnnotationHub下載注釋數(shù)據(jù)庫, 使用select(), mapIds等函數(shù)進(jìn)行注釋操作祝钢。我自己寫一個流程也用到了它給基因ID, 如AT1G14185, 注釋別名和功能描述比规。 注釋結(jié)果中會出現(xiàn)一些基因無法被注釋, 比如說下面這些情況, 我一直認(rèn)為只是這些基因沒有得到比較好的研究, 即便這些基因能夠在TAIR搜到Araport的注釋, 我也認(rèn)為那些注釋只是同源注釋沒有多大意義。
AT1G13970 NA NA
AT1G14120 NA NA
AT1G14240 NA NA
AT1G14600 NA NA
一開始得到的結(jié)果里沒有多少個基因太颤，所以缺少幾個注釋苞俘，通過手工去查找也行，但是目前差異表達(dá)分析動不動就給別人500多個基因龄章，于是就有幾十個甚至上百個未注釋的基因吃谣，所以我想著要不自己更新擬南芥的物種包。
library(AnnotationHub)
ah <- AnnotationHub()
org <- ah[["AH57965"]]
org
#...
# TAIRGENEURL: ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_functional_descriptions
#...
通過上面的代碼做裙，我找到了基因功能描述的數(shù)據(jù)庫來源文件岗憋，我下載了這個文件，并且拿用AnnotationHub注釋不到的功能的一個基因锚贱，"AT1G14185"仔戈，進(jìn)行測試
mapIds(org, "AT1G14185", "GENENAME", "TAIR")
這下就非常有趣了，在原始文件中能搜索到的基因用Bioconductor的物種注釋包時卻沒有注釋信息拧廊！為了搞清楚這個原因监徘，我花了快半個下午的時間去折騰，終于被我找到了原因吧碾。 我分別用一個能被org.At.tair.db注釋和一個不能被org.At.tair.db的注釋去搜索原始文本.
# 無注釋
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;
簡單的比較之后凰盔，你差不多就知道了org.At.tair.db的在功能描述這一部分其實只用第一列和第四列(為了方便展示我轉(zhuǎn)置了原始數(shù)據(jù))。這就是非常讓人意外了倦春，為啥它不用第一列和第三列呢户敬？ 我于是又去看了其他幾個基因落剪，就差不多明白了，原始的文本特別的混亂尿庐，你除了能保證第一列和第二列有信息外忠怖，其他列你根本無法保證，因此最好的策略以第一列作為檢索的關(guān)鍵字抄瑟，其他列合并成一列才行凡泣，然而作者沒有那么細(xì)致。
于是我就放棄了用org.At.tair.db注釋基因功能描述和基因別名了锐借，還是自己寫一個Python腳本進(jìn)行注釋吧问麸。
下面這個腳本只適用于bed格式的輸入，且第四列為轉(zhuǎn)錄本ID钞翔，另外兩個輸入文件分別為"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)
