eggNOG網(wǎng)站注釋蛋白序列得到文件query_seqs.fa.emapper.annotations
python環(huán)境
#git clone https://github.com/Hua-CM/HuaSmallTools.git下載注釋包(parse_go_obofile.py 和 parse_eggNOG.py)
#wget http://purl.obolibrary.org/obo/go/go-basic.obo下載Go數(shù)據(jù)庫的obo文件
輸入數(shù)據(jù)準(zhǔn)備
首先需要去GO下載GO的obo文件,這里我使用go-basic.obo然后?使用parse_go_obofile.py可以把obo文件解析為如下格式：
$python parse_go_obofile.py -i go-basic.obo -o go.tbold
$less go.tbold |awk 'BEGIN{FS="\t";OFS="\t"}{print $2,$1,$3}' > tmp && mv tmp go.tb (必不可少的一步柄错，不轉(zhuǎn)化parse_eggNOG.py會(huì)報(bào)錯(cuò), go.tb為parse_eggNOG.py用臨時(shí)文件呜达，富集時(shí)還是用go.tbold文件)
$python parse_eggNOG.py -i query_seqs.fa.emapper.annotations \
? ? ? ? ? ? ? ? ? ? ? -g go.tb \
? ? ? ? ? ? ? ? ? ? ? -O ath,osa \
? ? ? ? ? ? ? ? ? ? ? -o ./? #生成GO和KO數(shù)據(jù)庫
--------------------------參數(shù)說明----------------------------------------------------------
「-i」 eggNOG的注釋結(jié)果
「-g」 上一步根據(jù)obo解析出來的文件
「-O」 參考物種(只用于KEGG注釋饱亮，使用KEGG三字母物種縮寫表示).設(shè)置這個(gè)參數(shù)的原因是我做KEGG富集的時(shí)候發(fā)現(xiàn)有的基因會(huì)出現(xiàn)在非程睿荒唐的通路上，比如某個(gè)植物基因富集到了癌癥的相關(guān)通路，后來發(fā)現(xiàn)原因是有的比較基礎(chǔ)的KO可能與癌癥通路有關(guān)，如果不使用參考物種椎工，直接用KO去尋找map的話就會(huì)出現(xiàn)上述的情況。這里使用參考物種可以把沒有出現(xiàn)在參考物種中的通路給過濾掉。植物我選擇擬南芥和水稻作為參考晋渺，同樣的如果做非模式動(dòng)物的話镰绎，可以考慮設(shè)置一些動(dòng)物物種來排除富集到植物的通路上
「-o」 輸出結(jié)果文件夾。會(huì)在該文件夾生成GOannotation.tsv和KOannotation.tsv兩個(gè)文件
--------------------------------------------------------------------------------------------------
#(在界面直接跑python parse_eggNOG.py -i query_seqs.fa.emapper.annotations -g go.tb -O ath,osa -o ./)?處理的結(jié)果文件有兩個(gè)：「GOannotation.tsv」和「KOannotation.tsv」 分別對應(yīng)GO注釋和KO注釋,使用這兩個(gè)文件就可以進(jìn)行富集分析了
富集分析
#R環(huán)境
library(clusterProfiler)
KOannotation <- read.delim("KOannotation.tsv", stringsAsFactors=FALSE)
GOannotation <- read.delim("GOannotation.tsv", stringsAsFactors=FALSE)
GOinfo <- read.delim("go.tbold", stringsAsFactors=FALSE)
#根據(jù)具體比較組來木西，上面步驟不需要重做畴栖。
gene_list <- read.csv("test.txt",header=FALSE) ###提供差異基因數(shù)據(jù)集，不要表頭
gene_list <- gene_list$V1 ####提取基因所在列八千，如gene_list <- gene_list$Row.names吗讶，根據(jù)實(shí)際情況提取恋捆；此處為類型轉(zhuǎn)換照皆，因?yàn)閑nricher識別的gene_list為facter，charater類型沸停，而不識別data.drame膜毁；轉(zhuǎn)化為character類型的話#gene_list <- as.character(gene_list$V
Go富集
GOannotation = split(GOannotation, with(GOannotation, level))
df_GO <- enricher(gene_list,TERM2GENE=GOannotation[['MF']][c(2,1)],TERM2NAME=GOinfo[1:2])
df_GO_df <- as.data.frame(df_GO) ###易讀模式，查看GO富集結(jié)果
dotplot(df_GO)
####柱狀圖用barplot(df_GO, showCategory = 10)愤钾，10表示選擇前10個(gè)term###
###做有向無環(huán)圖plotGOgraph(df_GO)###要先安裝topGo####
####輸出表格write.csv(df_GO, "0_up1.csv")###
df_GO <- enricher(gene_list,TERM2GENE=GOannotation[['CC']][c(2,1)],TERM2NAME=GOinfo[1:2])
dotplot(df_GO)
df_GO <- enricher(gene_list,TERM2GENE=GOannotation[['BP']][c(2,1)],TERM2NAME=GOinfo[1:2])
dotplot(df_GO)
KEGG富集
enricher(gene_list,TERM2GENE=KOannotation[c(3,1)],TERM2NAME=KOannotation[c(3,4)])
df_KO <- enricher(gene_list,TERM2GENE=KOannotation[c(3,1)],TERM2NAME=KOannotation[c(3,4)])
dotplot(df_KO)
優(yōu)缺點(diǎn)
優(yōu)點(diǎn)
理論上針對所有有基因組蛋白序列的物種都可以注釋瘟滨，甚至沒有基因組但是有參考轉(zhuǎn)錄本也可以注釋
相對來說沒有用到特別多的依賴工具，兩個(gè)python腳本也只使用了最基本的包能颁。
缺點(diǎn)
中間需要解析eggNOG的結(jié)果文件杂瘸，我個(gè)人用的python寫的腳本，一是與R銜接不連貫伙菊，二是速度較慢败玉，主要是我考慮做一個(gè)背景數(shù)據(jù)集能用很久，所以偷懶沒有對腳本的性能進(jìn)行優(yōu)化镜硕。
無法一張表完成基因运翼、轉(zhuǎn)錄本等多層次的富集信息。主要是因?yàn)槲易龅闹参锖芏鄷r(shí)候基因組質(zhì)量不高谦疾，根本沒有多個(gè)轉(zhuǎn)錄本的注釋南蹂，所以就沒考慮這一問題。
引用：https://mp.weixin.qq.com/s/Mr3YLoc_-Y1WeLKJku1TzQ
