eggNOG网站注释蛋白序列得到文件query_seqs.fa.emapper.annotations
python环境
#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数据库的obo文件
输入数据准备
首先需要去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 (必不可少的一步，不转化parse_eggNOG.py会报错, go.tb为parse_eggNOG.py用临时文件，富集时还是用go.tbold文件)
$python parse_eggNOG.py -i query_seqs.fa.emapper.annotations \
-g go.tb \
-O ath,osa \
-o ./ #生成GO和KO数据库
--------------------------参数说明----------------------------------------------------------
「-i」 eggNOG的注释结果
「-g」 上一步根据obo解析出来的文件
「-O」 参考物种(只用于KEGG注释，使用KEGG三字母物种缩写表示).设置这个参数的原因是我做KEGG富集的时候发现有的基因会出现在非常荒唐的通路上，比如某个植物基因富集到了癌症的相关通路，后来发现原因是有的比较基础的KO可能与癌症通路有关，如果不使用参考物种，直接用KO去寻找map的话就会出现上述的情况。这里使用参考物种可以把没有出现在参考物种中的通路给过滤掉。植物我选择拟南芥和水稻作为参考，同样的如果做非模式动物的话，可以考虑设置一些动物物种来排除富集到植物的通路上
「-o」 输出结果文件夹。会在该文件夹生成GOannotation.tsv和KOannotation.tsv两个文件
--------------------------------------------------------------------------------------------------
#(在界面直接跑python parse_eggNOG.py -i query_seqs.fa.emapper.annotations -g go.tb -O ath,osa -o ./) 处理的结果文件有两个：「GOannotation.tsv」和「KOannotation.tsv」 分别对应GO注释和KO注释,使用这两个文件就可以进行富集分析了
富集分析
#R环境
library(clusterProfiler)
KOannotation <- read.delim("KOannotation.tsv", stringsAsFactors=FALSE)
GOannotation <- read.delim("GOannotation.tsv", stringsAsFactors=FALSE)
GOinfo <- read.delim("go.tbold", stringsAsFactors=FALSE)
#根据具体比较组来，上面步骤不需要重做。
gene_list <- read.csv("test.txt",header=FALSE) ###提供差异基因数据集，不要表头
gene_list <- gene_list$V1 ####提取基因所在列，如gene_list <- gene_list$Row.names，根据实际情况提取；此处为类型转换，因为enricher识别的gene_list为facter，charater类型，而不识别data.drame；转化为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富集结果
dotplot(df_GO)
####柱状图用barplot(df_GO, showCategory = 10)，10表示选择前10个term###
###做有向无环图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)
优缺点
优点
理论上针对所有有基因组蛋白序列的物种都可以注释，甚至没有基因组但是有参考转录本也可以注释
相对来说没有用到特别多的依赖工具，两个python脚本也只使用了最基本的包。
缺点
中间需要解析eggNOG的结果文件，我个人用的python写的脚本，一是与R衔接不连贯，二是速度较慢，主要是我考虑做一个背景数据集能用很久，所以偷懒没有对脚本的性能进行优化。
无法一张表完成基因、转录本等多层次的富集信息。主要是因为我做的植物很多时候基因组质量不高，根本没有多个转录本的注释，所以就没考虑这一问题。
引用：https://mp.weixin.qq.com/s/Mr3YLoc_-Y1WeLKJku1TzQ
网友评论