首先推薦Kobashttp://kobas.cbi.pku.edu.cn/kobas3/genelist/
上傳你的蛋白序列 genelist，背景基因選相似物種就可以做kegg纠俭，kegg的結(jié)果畫圖參考
https://mp.weixin.qq.com/s/fwZ4rZOg8PdVb2ZE9wskoQ

data01<-read.table("你物種的ko結(jié)果.txt",header = T,sep="\t")
pathway=data.frame(data01[,c(1,4,6,7)])
pathway=head(pathway,n=20)
names(pathway)<-c("PATHWAY","Gene","Pvalue","Qvalue")
library(ggplot2)
ggplot(pathway,aes(Pvalue,PATHWAY))+
  geom_point(aes(size=Gene,color=-1*log10(Qvalue)))+
  scale_color_gradient(low="green",high = "red")+ 
  labs(color=expression(-log[10](Qvalue)),size="Gene",  
       x="Pvalue",     
       y="Pathway name",
       title="Top20 of Pathway enrichment")+ 
  theme_bw() 
write.table(pathway,"simpathway.txt",row.names=F,col.names=T,sep="\t")

image.png

image.png

我同時(shí)也去嘗試kobas的go富集 但一直顯示這個(gè)畫面 不過kegg倒是挺快的

image.png

接下來可以用eggnog-mapper先注釋然后構(gòu)建自己的org.db 再用clusterProfiler包畫冰木，可以參考別人的做法 我做完感覺圖還是不好看 穷劈，構(gòu)建庫也是很慢
接下來做Go我使用eggnog-mapper先對(duì)我要的物種注釋獲得Go號(hào)
sim.emapper.annotations這是我的注釋結(jié)果命名為sim

vi chang_go.py
#!/usr/bin/python
# -*- coding: utf-8 -*-
import sys, os
x = sys.argv[1]
file = open(x, "r")
lines = file.readlines()
for line in lines:
    line=line.strip()  
    if line.startswith("#"):
        continue
    else:
        tmp=line.split("\t")
        if tmp[6] == "":
            continue
        else :
             mystr=tmp[0]+"\t"+tmp[6]
             print (mystr)


python chang_go.py sim >simout
sed 's/,/\t/g' simout >simout.txt

然后借用下某平臺(tái)簡單做一下
https://www.omicshare.com/tools/home/report/goenrich.html
把你的genelist和注釋的結(jié)果一放就好了 結(jié)果會(huì)有類似這種圖 和用R語言clusterProfiler包畫的差不多吧 感覺應(yīng)該可以這樣做吧 請(qǐng)大家試試

image.png

6.7號(hào)更新
用eggnog-mapper畫go富集分析代碼 不用構(gòu)建物種庫 快速運(yùn)行 節(jié)約時(shí)間

rm(list=ls())
setwd("D:")
BiocManager::install("stringr")
BiocManager::install("dplyr")
BiocManager::install("clusterProfiler")

library(stringr)
library(dplyr)
egg<-read.csv(".emapper.annotations",sep="\t",header=T)
egg[egg==""]<-NA  
gene_ids <- egg$query_name
eggnog_lines_with_go <- egg$GOs!= ""
eggnog_lines_with_go
eggnog_annoations_go <- str_split(egg[eggnog_lines_with_go,]$GOs, ",")
gene_to_go <- data.frame(gene = rep(gene_ids[eggnog_lines_with_go],
                                    times = sapply(eggnog_annoations_go, length)),
                         term = unlist(eggnog_annoations_go))
head(gene_to_go)
library(clusterProfiler)
gene_list <- read.table("ID")
gene_list <-as.vector(gene_list$V1)
term2gene<-gene_to_go[,c(2,1)]
term2gene=na.omit(term2gene)
df<-enricher(gene=gene_list,
             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             TERM2GENE = term2gene)
head(df)
barplot(df)
dotplot(df)
df<-as.data.frame(df)
dim(df)
df1<-go2term(df$ID)
dim(df1)
head(df1)
df$term<-df1$Term
df2<-go2ont(df$ID)
dim(df2)
head(df2)
df$Ont<-df2$Ontology
head(df)
df3<-df%>%
  select(c("term","Ont","pvalue"))
#df3=subset(df,select=c("term","Ont","pvalue"))
head(df3)
library(ggplot2)
ggplot(df3,aes(x=term,y=-log10(pvalue)))+
  geom_col(aes(fill=Ont))+
  coord_flip()+labs(x="")+
  theme_bw()
write.table(df3,"go_out",row.names=T,col.names=T,sep="\t")

Go號(hào)由一對(duì)多變?yōu)槎鄬?duì)多 然后換agriGo做

library(tidyverse)
setwd("D:/kobas注釋/go")
data01<-read.table("simout",header = F,sep="\t")
names(data01)<-c("gene","GO_ID")
data02<-data01 %>%
  tidyr::separate_rows(GO_ID, sep = ",")
write.table(data02,"simgo.out",quote = FALSE,row.names=F,col.names=F,sep="\t")

library(ggplot2)
CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")
p <- ggplot(data=data, aes(x=data$Term, y=LG10, fill=KEGG_A_class)) +
  geom_bar(stat="identity", width=0.8) + coord_flip() + 
  scale_fill_manual(values = CPCOLS) + theme_bw() + 
  scale_x_discrete(labels=data) +
  xlab("KGGG term") + 
  theme(axis.text=element_text(face = "bold", color="gray50")) +
  labs(title = "The Most Enriched KEGG Terms")

###go富集 注釋的分類圖
rm(list=ls())
setwd("D:/")
library(openxlsx)
data<- read.xlsx("Functions enriched by Blast2GO.xlsx", sheet = 1)
###處理為ID   Description   GeneNumber     type
data1=subset(data,select=c("GO.ID","GO.Name","Nr.test","GO.Category"))
library(ggplot2)
data1=data1[,c(2,3,4)]
colnames(data1)=c("name","number","type")
a=data1[data1$number>=10,]
a=a[order(a$type),]
ggplot(a,aes(x=factor(name,levels=unique(name)),y=number,fill=type))+geom_bar(stat="identity")+coord_flip()+
  theme_bw() + ylab("Number") + xlab("Name") 

###nr庫比對(duì)注釋圖
首先diamond比對(duì)上nr庫然后利用 UniProt2GO_annotate.py 腳本得到GO號(hào)
nohup diamond blastp -d nr.dmnd -q my.pep -o my.out -f 6 -p 10 --more-sensitive -e 1e-5  &
nohup python UniProt2GO_annotate.py ./idmapping/idmapping.tb.gz my.out my.go.out &
然后利用R語言
library(ggplot2)
options(stringsAsFactors = F)
library("dplyr")
library("stringr")
library(tidyverse)
library(GO.db)
male=read.table("my.go.out",header = F)
names(male)<-c("Gene","GOID")
male02<-male %>%
  tidyr::separate_rows(GOID, sep = ",")
godb <- select(GO.db, keys(GO.db), columns(GO.db))
go_male02 <- male02 %>% left_join(godb)
go_annot=go_male02
go_annot=go_annot[,-1]
go_annot=go_annot[,-2]
as=as.data.frame(table(go_annot$GOID))
colnames(as)=c("GOID","GeneNumber")
go_annot2 <- go_annot %>% left_join(as)
go_annot=go_annot2 %>% distinct()#根據(jù)所有列刪除重復(fù)的行（完全一樣的觀測值）
go_annot=go_annot[,c(1,3,2,4)]
colnames(go_annot)=c("ID","Description","GeneNumber","type")
go_enrich_df=go_annot
dat=go_enrich_df[,c(2,3,4)]
colnames(dat)=c("name","type","number")
a=dat[dat$number>=50,]
a=a[order(a$type),]
library(ggplot2)
ggplot(a,aes(x=factor(name,levels=unique(name)),y=number,fill=type))+geom_bar(stat="identity")+coord_flip()+
  theme_bw() + ylab("Number") + xlab("Name")   
###

KEGG注釋圖的繪制
首先上傳蛋白質(zhì)到KOBAS網(wǎng)站獲取注釋文件

import sys

input = sys.argv[1]
input_file = open(input, "r")

input_reads = input_file.read()
title = "query_name" + "\t" + "gene_name" + "\t" + "KO_number" + "\t" + "pathway"
print(title)
for lines in input_reads.split("Query"):
    if "Pathway" in lines:
        for line in lines.split("\n"):
            if line.startswith(":"):
                query_name = line.split("\t")[1]
            if line.startswith("Gene:"):
                gene_name = line.split("\t")[1]
            if "KEGG" in line:
                KO_number = str(line.split("\t")[3]).replace("pop", "")
                pathway = line.split("\t")[1]
                output_line = query_name + "\t" + gene_name + "\t" + KO_number + "\t" + pathway
                print(output_line)

input_file.close()
···
得到
query_name      gene_name       KO_number       pathway
Simonii00034466-RA      pop:7475026     01100   Metabolic pathways
Simonii00034466-RA      pop:7475026     01110   Biosynthesis of secondary metabolites
Simonii00034466-RA      pop:7475026     00564   Glycerophospholipid metabolism
Simonii00034466-RA      pop:7475026     00565   Ether lipid metabolism
Simonii00034466-RA      pop:7475026     00562   Inositol phosphate metabolism
Simonii00036281-RA      pop:7480371     03018   RNA degradation
Simonii00038638-RA      pop:7472530     04075   Plant hormone signal transduction
Simonii00038638-RA      pop:7472530     04016   MAPK signaling pathway - plant
Simonii00010134-RA      pop:7478365     01100   Metabolic pathways
如上ko號(hào) 和注釋通路

## KAAS在線注釋

首先介紹怎樣使用網(wǎng)站的在線注釋功能。

進(jìn)入KAAS：[https://www.genome.jp/kaas-bin/kaas_main](https://www.genome.jp/kaas-bin/kaas_main)