在非模式生物中，我们想要获取相应物种的GO或KEGG注释信息，或者描述信息的渠道很少。数据分布不集中，需要我们花费很长时间去收集，整合，才能得到相应的信息去做功能富集分析。下面以小麦为例，我们如何去收集功能富集前的准备工作。

一、功能注释

对一些未知功能的序列做初步的功能注释，通常的方法就是先下载一些功能数据库(NR、NT、uniprot、GO、KEGG)，然后将序列进行比对，从而预测功能。但是这些数据库通常比较庞大，下载慢，处理起来也比较复杂。因此，我们一般采用在线做功能注释。有些网站例如eggNOG-MAPPER、PANNZER、g:Profiler等。我认为在线功能注释网站中比较快的是eggNOG-MAPPER，PANNZER专注做功能注释约有10年了，也是比较专业的，大家也可以多在几个网站中做，然后比较一下结果。取自己想要的，合适的结果。

下面就以eggNOG-MAPPER网站为例，做以下序列的功能注释：

示例数据取自小麦中国春v2.1版本的高可信蛋白序列。一般我们以蛋白序列作为功能注释的输入文件，示例文件名：Example_pep.fa

image.png

image.png

上面为填写后的网页，我一般采用默认参数，然后提交。在邮箱中会看到下面界面。

image.png

image.png

image.png

当运行完成之后会给你邮箱发送一个任务完成邮件。
任务完成啦，然后下载结果。

image.png

image.png

image.png

image.png

image.png

image.png

out.emapper.annotations.GO文件里，有GO与蛋白质ID对应关系，这样我们就初步得到我们所需要的文件啦。

image.png

image.png

image.png

grep map out.emapper.annotations.KEGG_Pathway.txt > out.emapper.annotations.KEGG_map.tx
grep ko out.emapper.annotations.KEGG_Pathway.txt > out.emapper.annotations.KEGG_ko.txt

这样我们就能提取出只有map号和KO号的文件啦。

image.png

image.png

到此我们前景基因的功能注释就完成啦。

二、获取该物种所对应的GO、KEGG注释

获取物种与GO号所对应的关系的方法有以下几种方式：

1、通过GO官网进行下载http://current.geneontology.org/annotations/index.html或者http://current.geneontology.org/products/pages/downloads.html下载对应物种的注释信息。很可惜，官网里面并没有小麦的注释信息。

2、 从COA项目中下载ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/，但是也没有小麦的注释信息。

image.png

image.png

image.png

三、GO、KEGG描述信息

GO描述信息可以在http://current.geneontology.org/ontology/go-basic.obo得到，是所有物种的描述信息，不影响做GO富集分析。进入下面页面后右击点击另存为，保存文件。然后提取文件中信息即可得到简化的描述信息。

image.png

with open("go-basic.obo","r") as file:
    lib={}
    for line in file:
        line=line.strip()
        col_name=line.split(":")[0]
        if col_name == "id":
            id=line.split(" ",maxsplit=1)[1]
            lib[id]=""
        if col_name == "name":
            name=line.split(" ",maxsplit=1)[1]
            lib[id]=lib[id]+"@"+name
        if col_name == "namespace":
            namespace=line.split(" ",maxsplit=1)[1]
            lib[id]=lib[id]+"@"+namespace
out=open("GO_basic_Description.txt","a+")
out.write("Class"+"\t"+"GO_IDs"+"\t"+"Description"+"\n")
for key in lib.keys():
    go_id=key
    go_name=lib[key].split("@")[1]
    go_namespace=lib[key].split("@")[2]
    if go_namespace == "molecular_function":
        go_namespace="MF"
        out.write(go_namespace+"\t"+go_id+"\t"+go_name+"\n")
    if go_namespace == "biological_process":
        go_namespace="BP"
        out.write(go_namespace+"\t"+go_id+"\t"+go_name+"\n")
    if go_namespace == "cellular_component":
        go_namespace="CC"
        out.write(go_namespace+"\t"+go_id+"\t"+go_name+"\n")

得到下面的文件 GO_basic_Description.txt：

image.png

image.png

image.png

image.png

image.png

image.png

image.png

image.png

#获取KEGG数据库信息
#加载包
library(KEGGREST)
##查看KEGG数据库包含的数据
listDatabases()#
##获取pathway(所有物种)数据集中的数据
pathway<- keggList("pathway")
head(pathway)
#转换数据集，导出数据集
pathway_data<-as.data.frame(pathway)
write.table(pathway_data,"<path>/KEGG_pathway_allspacies_database.txt",row.names = T,col.names = F,sep = "\t")
##对单个(小麦)数据库进行物种的选择
taes_pathway <-keggList("pathway","taes")#taes是小麦的缩写
taes_pathway_data<-as.data.frame(taes_pathway)
write.table(taes_pathway_data,"<path>/KEGG_pathway_wheat_database.txt",row.names = T,col.names = F,sep = "\t")
##获取KO(所有基因)数据集中的数据
Ko<-keggList("ko")
Ko_data<-as.data.frame(Ko)
write.table(Ko_data,"<path>/KEGG_KO_allspacies_database.txt",row.names = T,col.names = F,sep = "\t")
##获取KO(小麦)数据集中的数据
taes_Ko<-keggList("ko","taes")
taes_Ko_data<-as.data.frame(taes_Ko)
write.table(taes_Ko_data,"<path>/KEGG_KO_wheat_database.txt",row.names = T,col.names = F,sep = "\t")

得到文件如下：
KEGG_pathway_allspacies_database.txt
KEGG_pathway_wheat_database.txt
KEGG_KO_allspacies_database.txt
KEGG_KO_wheat_database.txt
这个包不仅可以提取map号、ko号还可以提取其他的东西。下图就是所有可提的部分，只要修改以下keggList("pathway","taes")参数即可。

image.png

得到KEGG的描述信息后需要用一个python脚本简单处理一下，输入文件为KEGG_pathway_allspacies_database.txt、KEGG_KO_allspacies_database.txt

#将得到的KO、pathway文件进行处理，得到可以富集的描述信息文件
with open("KEGG_KO_allspacies_database.txt","r") as file:
    out=open("KEGG_KO_allspacies_description.txt","a+")
    for line in file:
        line=line.strip()
        if "[EC:" in line.split("\t",maxsplit=1)[1] :
            desc=line.split("\t",maxsplit=1)[1].rsplit("[EC:",maxsplit=1)[0].split(";")[1].strip('"')
            id=line.split("\t",maxsplit=1)[0].split(":")[1][:-1]
            out.write(id+"\t"+desc+"\n")
        else:
            if ";" in line.split("\t",maxsplit=1)[1]:
                out.write(line.split("\t",maxsplit=1)[0][:-1].split(":")[1]+"\t"+line.split("\t",maxsplit=1)[1].split(";")[1].strip('"')+"\n")
            else:
                out.write(line.split("\t",maxsplit=1)[0][:-1].split(":")[1]+"\t"+line.split("\t",maxsplit=1)[1].strip('"')+"\n")
with open("KEGG_pathway_allspacies_database.txt") as file1:
    out1=open("KEGG_pathway__allspacies_description.txt","a+")
    for line1 in file1:
        line1=line1.strip()
        map_id=line1.split("\t") [0].strip('"').split(":")[1]
        descript=line1.split("\t") [1].strip('"')
        out1.write(map_id+"\t"+descript+"\n")

输出文件为KEGG_KO_allspacies_description.txt、KEGG_pathway__allspacies_description.txt

image.png

image.png

这样我们GO、KEGG富集前准备文件就完成啦!