在非模式生物中，我們想要獲取相應(yīng)物種的GO或KEGG注釋信息，或者描述信息的渠道很少。數(shù)據(jù)分布不集中，需要我們花費(fèi)很長(zhǎng)時(shí)間去收集，整合，才能得到相應(yīng)的信息去做功能富集分析。下面以小麥為例，我們?nèi)绾稳ナ占δ芨患暗臏?zhǔn)備工作。

一、功能注釋

對(duì)一些未知功能的序列做初步的功能注釋，通常的方法就是先下載一些功能數(shù)據(jù)庫(kù)(NR、NT、uniprot、GO、KEGG)，然后將序列進(jìn)行比對(duì)，從而預(yù)測(cè)功能。但是這些數(shù)據(jù)庫(kù)通常比較龐大，下載慢，處理起來(lái)也比較復(fù)雜。因此，我們一般采用在線做功能注釋。有些網(wǎng)站例如eggNOG-MAPPER、PANNZER、g:Profiler等。我認(rèn)為在線功能注釋網(wǎng)站中比較快的是eggNOG-MAPPER，PANNZER專注做功能注釋約有10年了，也是比較專業(yè)的，大家也可以多在幾個(gè)網(wǎng)站中做，然后比較一下結(jié)果。取自己想要的，合適的結(jié)果。

下面就以eggNOG-MAPPER網(wǎng)站為例，做以下序列的功能注釋：

示例數(shù)據(jù)取自小麥中國(guó)春v2.1版本的高可信蛋白序列。一般我們以蛋白序列作為功能注釋的輸入文件，示例文件名：Example_pep.fa

image.png

image.png

上面為填寫(xiě)后的網(wǎng)頁(yè)，我一般采用默認(rèn)參數(shù)，然后提交。在郵箱中會(huì)看到下面界面。

image.png

image.png

image.png

當(dāng)運(yùn)行完成之后會(huì)給你郵箱發(fā)送一個(gè)任務(wù)完成郵件。
任務(wù)完成啦，然后下載結(jié)果。

image.png

image.png

image.png

image.png

image.png

image.png

out.emapper.annotations.GO文件里，有GO與蛋白質(zhì)ID對(duì)應(yīng)關(guān)系，這樣我們就初步得到我們所需要的文件啦。

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號(hào)和KO號(hào)的文件啦。

image.png

image.png

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

二、獲取該物種所對(duì)應(yīng)的GO、KEGG注釋

獲取物種與GO號(hào)所對(duì)應(yīng)的關(guān)系的方法有以下幾種方式：

1、通過(guò)GO官網(wǎng)進(jìn)行下載http://current.geneontology.org/annotations/index.html或者http://current.geneontology.org/products/pages/downloads.html下載對(duì)應(yīng)物種的注釋信息。很可惜，官網(wǎng)里面并沒(méi)有小麥的注釋信息。

2、 從COA項(xiàng)目中下載ftp://ftp.ebi.ac.uk/pub/databases/GO/goa/，但是也沒(méi)有小麥的注釋信息。

image.png

image.png

image.png

三、GO、KEGG描述信息

GO描述信息可以在http://current.geneontology.org/ontology/go-basic.obo得到，是所有物種的描述信息，不影響做GO富集分析。進(jìn)入下面頁(yè)面后右擊點(diǎn)擊另存為，保存文件。然后提取文件中信息即可得到簡(jiǎn)化的描述信息。

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數(shù)據(jù)庫(kù)信息
#加載包
library(KEGGREST)
##查看KEGG數(shù)據(jù)庫(kù)包含的數(shù)據(jù)
listDatabases()#
##獲取pathway(所有物種)數(shù)據(jù)集中的數(shù)據(jù)
pathway<- keggList("pathway")
head(pathway)
#轉(zhuǎn)換數(shù)據(jù)集，導(dǎo)出數(shù)據(jù)集
pathway_data<-as.data.frame(pathway)
write.table(pathway_data,"<path>/KEGG_pathway_allspacies_database.txt",row.names = T,col.names = F,sep = "\t")
##對(duì)單個(gè)(小麥)數(shù)據(jù)庫(kù)進(jìn)行物種的選擇
taes_pathway <-keggList("pathway","taes")#taes是小麥的縮寫(xiě)
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(所有基因)數(shù)據(jù)集中的數(shù)據(jù)
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(小麥)數(shù)據(jù)集中的數(shù)據(jù)
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
這個(gè)包不僅可以提取map號(hào)、ko號(hào)還可以提取其他的東西。下圖就是所有可提的部分，只要修改以下keggList("pathway","taes")參數(shù)即可。

image.png

得到KEGG的描述信息后需要用一個(gè)python腳本簡(jiǎn)單處理一下，輸入文件為KEGG_pathway_allspacies_database.txt、KEGG_KO_allspacies_database.txt

#將得到的KO、pathway文件進(jìn)行處理，得到可以富集的描述信息文件
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富集前準(zhǔn)備文件就完成啦!