Bioconductor開發(fā)的物種注釋包系列集合了一個物種不同來源的注釋信息，能夠根據(jù)基因ID對其進行多種來源的注釋，比如說基因的別名，基因的功能等。

我之前也寫過一篇文章用Bioconductor對基因組注釋介紹如何使用AnnotationHub下載注釋數(shù)據(jù)庫, 使用select(), mapIds等函數(shù)進行注釋操作。我自己寫一個流程也用到了它給基因ID, 如AT1G14185, 注釋別名和功能描述。 注釋結(jié)果中會出現(xiàn)一些基因無法被注釋, 比如說下面這些情況, 我一直認為只是這些基因沒有得到比較好的研究, 即便這些基因能夠在TAIR搜到Araport的注釋, 我也認為那些注釋只是同源注釋沒有多大意義。

AT1G13970   NA  NA
AT1G14120   NA  NA
AT1G14240   NA  NA
AT1G14600   NA  NA

一開始得到的結(jié)果里沒有多少個基因，所以缺少幾個注釋，通過手工去查找也行，但是目前差異表達分析動不動就給別人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”，進行測試

mapIds(org, "AT1G14185", "GENENAME", "TAIR")

img

這下就非常有趣了，在原始文件中能搜索到的基因用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)鍵字，其他列合并成一列才行，然而作者沒有那么細致。

于是我就放棄了用org.At.tair.db注釋基因功能描述和基因別名了，還是自己寫一個Python腳本進行注釋吧。

下面這個腳本只適用于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 sysfrom 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)