eggNOG-mapper功能注释

http://eggnog-mapper.embl.de/

eggNOG ( Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups) 公共数据库V5.0版本搜集了5090个生物（477真核生物、4445个代表性细菌和168古菌）和2502个病毒的全基因组蛋白序列。将这些物种分成了379类（taxonomic levels），每类的编号以NCBI的分类编号表示。对各个类别物种的全基因组蛋白序列进行同源基因分类，总共得到4.4M个同源基因类。eggNOG数据库是NCBI的COG数据库的扩展，它收集了更全面的物种和更大量的蛋白序列数据。
总结一下就是，信息更全更大。并且提供了本地化软件和网页工具进行eggNOG注释工具。

图片.png

1.一般是输入蛋白序列
2.选择本地输入文件，即蛋白序列集合
3.给定一个邮箱地址
4.点击提交
5.进入邮箱，点击Start
6.等待结果，下载推荐csv格式

图片.png

图片.png

图片.png

query，seed_ortholog，evalue，score，eggNOG_OGs，max_annot_lvl，COG_category，Description
Preferred_name，GOs，EC，KEGG_ko，KEGG_Pathway，KEGG_Module，KEGG_Reaction，KEGG_rclass
BRITE，KEGG_TC，CAZy，BiGG_Reaction，PFAMs

非模式物种OrgDB的构建

ClusterProfiler是一个功能强大的R包，同时支持GO和KEGG的富集分析，而且可视化功能非常的优秀，那么怎么利用这个R包来进行富集分析。
进行分析时，首先就要求注释信息。Bioconductor上提供了以下19个物种的Org类型的包，包含了这些模式物种的GO及KEGG注释信息

图片.png

对于以上19个物种，只需要安装对应的org包，clusterProfile就会自动从中获取GO注释信息，我们只需要差异基因的列表就可以了，使用起来非常方便。但是没有Orgdb的物种，就需要自己构建或者利用他人上传的数据，这里建议自己构建，不存在ID的错误问题。
使用eggNOG-mapper结果进行构建：

####注释Orgdb
#install.packages("")
library(tidyverse)
library(stringr)
library(KEGGREST)
library(AnnotationForge)
library(clusterProfiler)
library(dplyr)
library(jsonlite)
library(purrr)
library(RCurl)
#顺手设置一下options
options(stringsAsFactors = F)
#顺手设置一下options
emapper <- read.table("out.emapper.annotations.txt",header = TRUE, sep = "\t")
#替换空格值
emapper[emapper==""]<-NA

gene_info <- emapper %>% dplyr::select(GID = query, GENENAME = Preferred_name) %>% na.omit()
gos <- emapper %>% dplyr::select(query, GOs) %>% na.omit()

gene2go = data.frame(GID = character(), GO = character(), EVIDENCE = character())

for (row in 1:nrow(gos)) {
  the_gid <- gos[row,"query"][[1]]
  the_gos <- str_split(gos[row,"GOs"], ",", simplify = FALSE)[[1]]
  df_temp <- data_frame(GID = rep(the_gid, length(the_gos)),
                        GO = the_gos,
                        EVIDENCE = rep("IDA", length(the_gos)))
  gene2go <- rbind(gene2go, df_temp)}

gene2go$GO[gene2go$GO=="-"]<-NA
gene2go<-na.omit(gene2go)

#Step3:提取KEGG信息

gene2ko <- emapper %>% dplyr::select(GID = query, Ko = KEGG_ko) %>% na.omit()

gene2ko$Ko <- gsub("ko:","",gene2ko$Ko)

getwd()


update_kegg <- function(json = "ko00001.json") {
  pathway2name <- tibble(Pathway = character(), Name = character())
  ko2pathway <- tibble(Ko = character(), Pathway = character())
  kegg <- fromJSON(json)
  for (a in seq_along(kegg[["children"]][["children"]])) {
    A <- kegg[["children"]][["name"]][[a]]
    for (b in seq_along(kegg[["children"]][["children"]][[a]][["children"]])) {
      B <- kegg[["children"]][["children"]][[a]][["name"]][[b]] 
      for (c in seq_along(kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]])) {
        pathway_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["name"]][[c]]
        pathway_id <- str_match(pathway_info, "ko[0-9]{5}")[1]
        pathway_name <- str_replace(pathway_info, " \\[PATH:ko[0-9]{5}\\]", "") %>% str_replace("[0-9]{5} ", "")
        pathway2name <- rbind(pathway2name, tibble(Pathway = pathway_id, Name = pathway_name))
        kos_info <- kegg[["children"]][["children"]][[a]][["children"]][[b]][["children"]][[c]][["name"]]
        kos <- str_match(kos_info, "K[0-9]*")[,1]
        ko2pathway <- rbind(ko2pathway, tibble(Ko = kos, Pathway = rep(pathway_id, length(kos))))}}}
  save(pathway2name, ko2pathway, file = "kegg_info.RData")}
# 调用函数后在本地创建kegg_info.RData文件
update_kegg()
# 载入kegg_info.RData文件
load(file = "kegg_info.RData")

gene2pathway <- gene2ko %>% left_join(ko2pathway, by = "Ko") %>% dplyr::select(GID, Pathway) %>% na.omit()

tax_id = "NCBI查询"
genus = "your genus "
species = "your species"

makeOrgPackage(gene_info=gene_info,
               go=gene2go,
               ko=gene2ko,
               pathway=gene2pathway,
               version="0.0.2",  #版本
               maintainer = "name <邮箱>",  #修改为你的名字和邮箱
               author = "name <邮箱>",  #修改为你的名字和邮箱
               outputDir = ".",  #输出文件位置
               tax_id=tax_id,  #你在NCBI上查并定义的id
               genus=genus,  #你在NCBI上查并定义的属名
               species=species,  #你在NCBI上查并定义的种名
               goTable="go")
#安装org.your_genus.eg.db
install.packages("./org.your_genus.eg.db/", repos = NULL, type = "sources")
#加载org.your_genus.eg.db
library(org.your_genus.eg.db)
#查看信息
columns(org.your_genus.eg.db)

go分析

将每个簇的所有marker名称提取出来

library(org.your_genus.eg.db)
#查看数据库总共有多少基因
length(keys(org.your_genus.eg.db))
gene <- read.csv("14cluster",header=FALSE) #读取一个基因列表
gene_list <- gene[,1]
de_ego <- enrichGO(gene = gene_list,
                   OrgDb = org.Zjujuba.eg.db,
                   keyType = 'GID',
                   ont = 'ALL',
                   qvalueCutoff = 0.05,
                   pvalueCutoff = 0.05)
pdf(file="dotplotcluster14.pdf",width = 15,height = 12)
dotplot(de_ego)
dev.off()
de_ego_df <- as.data.frame(de_ego)
head(de_ego_df)
write.table(de_ego_df,file = "./cluster14.results.txt",quote=F)

绘制条形图

pdf(file="barplotcluster14.pdf",width = 15,height = 12)
barplot(de_ego, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()

绘制点图

pdf(file="bubblecluster14.pdf",width = 15,height = 12)
dotplot(de_ego,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()