eggNOG-mapper功能注釋
http://eggnog-mapper.embl.de/
eggNOG ( Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups) 公共數(shù)據(jù)庫(kù)V5.0版本搜集了5090個(gè)生物（477真核生物、4445個(gè)代表性細(xì)菌和168古菌）和2502個(gè)病毒的全基因組蛋白序列落竹。將這些物種分成了379類（taxonomic levels）冈绊，每類的編號(hào)以NCBI的分類編號(hào)表示。對(duì)各個(gè)類別物種的全基因組蛋白序列進(jìn)行同源基因分類晚树，總共得到4.4M個(gè)同源基因類。eggNOG數(shù)據(jù)庫(kù)是NCBI的COG數(shù)據(jù)庫(kù)的擴(kuò)展，它收集了更全面的物種和更大量的蛋白序列數(shù)據(jù)裤纹。
總結(jié)一下就是，信息更全更大丧没。并且提供了本地化軟件和網(wǎng)頁(yè)工具進(jìn)行eggNOG注釋工具鹰椒。
1.一般是輸入蛋白序列
2.選擇本地輸入文件，即蛋白序列集合
3.給定一個(gè)郵箱地址
4.點(diǎn)擊提交
5.進(jìn)入郵箱呕童，點(diǎn)擊Start
6.等待結(jié)果漆际，下載推薦csv格式
結(jié)果如下，共有21列夺饲，后面就可以根據(jù)此文件進(jìn)行功能富集分析奸汇，將#號(hào)注釋的行和query前的#刪掉即可：
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的構(gòu)建
ClusterProfiler是一個(gè)功能強(qiáng)大的R包沸呐，同時(shí)支持GO和KEGG的富集分析，而且可視化功能非常的優(yōu)秀呢燥，那么怎么利用這個(gè)R包來進(jìn)行富集分析崭添。
進(jìn)行分析時(shí)，首先就要求注釋信息叛氨。Bioconductor上提供了以下19個(gè)物種的Org類型的包呼渣，包含了這些模式物種的GO及KEGG注釋信息
對(duì)于以上19個(gè)物種，只需要安裝對(duì)應(yīng)的org包寞埠，clusterProfile就會(huì)自動(dòng)從中獲取GO注釋信息屁置，我們只需要差異基因的列表就可以了，使用起來非常方便仁连。但是沒有Orgdb的物種蓝角，就需要自己構(gòu)建或者利用他人上傳的數(shù)據(jù)，這里建議自己構(gòu)建饭冬，不存在ID的錯(cuò)誤問題使鹅。
使用eggNOG-mapper結(jié)果進(jìn)行構(gòu)建：
####注釋Orgdb
#install.packages("")
library(tidyverse)
library(stringr)
library(KEGGREST)
library(AnnotationForge)
library(clusterProfiler)
library(dplyr)
library(jsonlite)
library(purrr)
library(RCurl)
#順手設(shè)置一下options
options(stringsAsFactors = F)
#順手設(shè)置一下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")}
# 調(diào)用函數(shù)后在本地創(chuàng)建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分析
將每個(gè)簇的所有marker名稱提取出來
library(org.your_genus.eg.db)
#查看數(shù)據(jù)庫(kù)總共有多少基因
length(keys(org.your_genus.eg.db))
gene <- read.csv("14cluster",header=FALSE) #讀取一個(gè)基因列表
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()
繪制點(diǎn)圖
pdf(file="bubblecluster14.pdf",width = 15,height = 12)
dotplot(de_ego,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
