非模式物种的一种简单版的go，kegg富集

首先推荐Kobashttp://kobas.cbi.pku.edu.cn/kobas3/genelist/
上传你的蛋白序列 genelist，背景基因选相似物种就可以做kegg，kegg的结果画图参考
https://mp.weixin.qq.com/s/fwZ4rZOg8PdVb2ZE9wskoQ

data01<-read.table("你物种的ko结果.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

我同时也去尝试kobas的go富集 但一直显示这个画面 不过kegg倒是挺快的


image.png

接下来可以用eggnog-mapper先注释然后构建自己的org.db 再用clusterProfiler包画，可以参考别人的做法 我做完感觉图还是不好看 ，构建库也是很慢
接下来做Go我使用eggnog-mapper先对我要的物种注释获得Go号
sim.emapper.annotations这是我的注释结果命名为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

然后借用下某平台简单做一下
https://www.omicshare.com/tools/home/report/goenrich.html
把你的genelist和注释的结果一放就好了 结果会有类似这种图 和用R语言clusterProfiler包画的差不多吧 感觉应该可以这样做吧 请大家试试

image.png

6.7号更新
用eggnog-mapper画go富集分析代码 不用构建物种库 快速运行 节约时间

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号由一对多变为多对多 然后换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库比对注释图
首先diamond比对上nr库然后利用 UniProt2GO_annotate.py 脚本得到GO号
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()#根据所有列删除重复的行（完全一样的观测值）
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注释图的绘制
首先上传蛋白质到KOBAS网站获取注释文件

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号 和注释通路

## KAAS在线注释

首先介绍怎样使用网站的在线注释功能。

进入KAAS：[https://www.genome.jp/kaas-bin/kaas_main](https://www.genome.jp/kaas-bin/kaas_main)
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 202,980评论 5赞 476
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,178评论 2赞 380
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 149,868评论 0赞 336
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,498评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,492评论 5赞 364
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,521评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,910评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,569评论 0赞 256
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,793评论 1赞 296
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,559评论 2赞 319
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,639评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,342评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,931评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,904评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,144评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 42,833评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,350评论 2赞 342