非模式生物GO、KEGG富集分析

GO、KEGG富集分析是我们做生信分析较为常用的部分,它可以将基因与功能相联系起来。
GO指的是Gene Ontology,是基因功能国际标准分类体系。目的在于建立一个适用于各种物种的,对基因和蛋白质功能进行限定和描述的,并能随着研究不断深入而更新的语言词汇标准。GO分为分子功能(Molecular Function)(MF)、生物过程(Biological Process)(BP)、和细胞组成(Cellular Component)(CC)三个部分。
KEGG指的是京都基因与基因组百科全书,通常我们使用KEGG中的pathway模块,将基因映射到某些通路上,了解基因参与生物体中的代谢过程等。
对于模式生物,GO和KEGG富集分析实现起来比较容易,对于非模式生物来说还是需要花点时间和精力。对于模式生物的GO和KEGG富集分析,网上教程案例挺多的。对于非模式生物,以小麦为例,进行下面一些基本的富集分析。

一、基本概念

做富集分析,我们需要了解一下几个概念。
1、前景基因:指的是我们所要进行富集的基因,一般是基因的ID
2、背景基因:指的是前景基因在某个基因集合进行富集,这个基因集合就是背景基因


v2-50cbf96d15658a93c6d680a1c9436eb7_720w.png

3、描述信息:每个GO的Term的属性,或者是每个KO号或者map号的属性。


image.png

image.png

我们具备前景基因,背景基因以及描述信息我们就可以做富集分析啦。

二、文件准备

1、前景基因:这是必须的啦。有时候需要进行ID转换,但是个人觉得ID转换根据需要来就行。如果前景基因里面的基因ID是包括在背景基因里面,那就需要进行转换。如果前景基因在是新的基因或者在背景基因没有被注释到的,就不用进行ID转换。下面这个就是融合基因,在背景基因里面没有注释到的,那么我就不要转换。


image.png

2、背景基因:一个基因可能具备多个GO term,一个基因也可能参与多个通路,与之相对应的有多个map号
这个案例中背景基因文件构建思路如下图


image.png

Knum是与组成成分有关的,就好像pathway是表示通路的一样。由于我们前景基因是融合基因的,那么我们就需要将前景基因与对应GO号、map号,k num加入到对应的背景基因文件中去,构成一个新的背景文件。
image.png

image.png

image.png

3、描述文件


image.png

image.png

image.png

上面的文件是有固定格式的。
三、富集分析
准备好这些文件之后进开始进行富集分析啦,使用下面R代码进行富集分析,主要使用clusterProfiler包进行富集分析,绘制一些简单的图。需要绘制高级图话,根据结果进行绘制。下面脚本里面参数根据自己需要可进行修改。
library("clusterProfiler")
library("ggplot2")
#导入数据
setwd("<path>")#设置工作目录
GO_enrichment_gene<-read.delim("<GO前景基因路径文件>",header = F)#导入GO前景文件
GO_enrichment_gene<-as.factor(GO_enrichment_gene$V1)#转化为因子
GO_background_genes<-read.table(<GO背景基因路径文件>,header = F)#导入GO背景文件
GO_Description<-read.table(<GO描述文件>,header = T)#导入GO描述文件
GO_description<-GO_Description[,-1]
#GO富集分析
GO_enrich_analysis <- enricher(GO_enrichment_gene,TERM2GENE=GO_background_genes,TERM2NAME=GO_description,pvalueCutoff = 0.01, qvalueCutoff = 0.05)
GO_enrich_analysis_data<-as.data.frame(GO_enrich_analysis)
for (id in GO_enrich_analysis_data$ID){
  if (is.na(GO_enrich_analysis_data[id,"Description"])){
    GO_enrich_analysis_data[id,"Class"]<-"NA"
  }
  else{
    GO_enrich_analysis_data[id,"Class"]<-GO_Description$Class[GO_Description$GO_ID==id]
  }
}
#绘制气泡图
dotplot(GO_enrich_analysis,showCategory = 10)
pdf("GOenrich_dotplot.pdf",width = 10,height = 10)
dotplot(GO_enrich_analysis,showCategory = 10)
dev.off()
#绘制条形图
barplot(GO_enrich_analysis,showCategory = 10)
pdf("GOenrich_barplot.pdf",width = 10,height = 10)
barplot(GO_enrich_analysis,showCategory = 10)
dev.off()
##绘制GO二级分类图
#选取出每一级P.adjust最小的前10个
#去除Description列中NA行
GO_2_classification<-GO_enrich_analysis_data[!is.na(GO_enrich_analysis_data$Description),]
#提取出各二级分类数据
GO_2_classification1<-GO_2_classification[GO_2_classification$Class=="BP",][order(GO_2_classification[GO_2_classification$Class=="BP",]$p.adjust)[1:10],]
GO_2_classification2<-GO_2_classification[GO_2_classification$Class=="CC",][order(GO_2_classification[GO_2_classification$Class=="CC",]$p.adjust)[1:10],]
GO_2_classification3<-GO_2_classification[GO_2_classification$Class=="MF",][order(GO_2_classification[GO_2_classification$Class=="MF",]$p.adjust)[1:10],]
GO_2_classification_1_2_3<-rbind(GO_2_classification1,GO_2_classification2)
#合并上面三个数据集
GO_2_classification_1_2_3<-rbind(GO_2_classification_1_2_3,GO_2_classification3)
#绘制二级分类图
ggplot(GO_2_classification_1_2_3,aes(x=reorder(Description,Count),y=Count,fill=-log10(p.adjust)))+
  geom_bar(stat="identity")+#由于是Description为变量,分类变量所以选择stat="identity"
  coord_flip()+#颠倒x和y轴
  scale_fill_gradient(low = "blue",high = "red")+#设置数值变量填充颜色,若要设置颜色可用scale_color_gradient函数
  facet_grid(Class~.,scale="free")+#分出刻面按class列分
  labs(x="Term",y="Counts",fill="FDR(-log10(P.adjust))")+#改变x,y轴标签,图例名字
  theme(axis.title=element_text(size = 15),axis.text.y = element_text(size = 13))#设置x,y轴标签大小以及y轴刻度名字
pdf("GOenrich_second_class.pdf",width = 15,height = 10)
ggplot(GO_2_classification_1_2_3,aes(x=reorder(Description,Count),y=Count,fill=-log10(p.adjust)))+
  geom_bar(stat="identity")+#由于是Description为变量,分类变量所以选择stat="identity"
  coord_flip()+#颠倒x和y轴
  scale_fill_gradient(low = "blue",high = "red")+#设置数值变量填充颜色,若要设置颜色可用scale_color_gradient函数
  facet_grid(Class~.,scale="free")+#分出刻面按class列分
  labs(x="Term",y="Counts",fill="FDR(-log10(P.adjust))")+#改变x,y轴标签,图例名字
  theme(axis.title=element_text(size = 15),axis.text.y = element_text(size = 13))#设置x,y轴标签大小以及y轴刻度名字
dev.off()
#导出GO富集结果
write.table(GO_enrich_analysis_data ,"GO_enrichment_results.txt",row.names = F,col.names = T,sep="\t")
##KEGG Psthway富集分析
#导入数据
KEGG_enrichment_gene<-read.table(<KEGGpathway前景文件>,header = F)#导入pathway前景文件
KEGG_background_genes<-read.table(<KEGGpathway背景文件>,header = F)#导入pathway背景文件
KEGG_description<-read.delim(<KEGGpathway描述文件>,header = F,sep = "\t")#导入pathway描述文件
KEGG_enrichment_gene<-as.factor(KEGG_enrichment_gene$V1)
KEGG_enrich_analysis<-enricher(KEGG_enrichment_gene,TERM2GENE = KEGG_background_genes,TERM2NAME = KEGG_description,pvalueCutoff = 0.01, qvalueCutoff = 0.05)
KEGG_enrich_analysis_data<-as.data.frame(KEGG_enrich_analysis)
#绘制气泡图
dotplot(KEGG_enrich_analysis,showCategory = 10)
pdf("KEGG_pathway_enrich_dotplot.pdf",width = 10,height = 10)
dotplot(KEGG_enrich_analysis,showCategory = 10)
dev.off()
#绘制条形图
barplot(KEGG_enrich_analysis,showCategory = 10)
pdf("KEGG_pathway_enrich_barplot.pdf",width = 10,height = 10)
barplot(KEGG_enrich_analysis,showCategory = 10)
dev.off()
#导出KEGG pathway富集结果
write.table(KEGG_enrich_analysis_data,"KEGG_pathway_enrichment_results.txt",col.names = T,row.names = F,sep = "\t")
####KEGG ko富集分析
KEGG_KO_enrichment_gene<-read.table(<KEGG k num 前景文件>,header = F)
KEGG_KO_background_genes<-read.table(<KEGG k num 背景文件>,header = F)
KEGG_KO_description<-read.delim(<KEGG k num 描述文件>,header = F,sep = "\t")
KEGG_KO_enrichment_gene<-as.factor(KEGG_KO_enrichment_gene$V1)
KEGG_KO_enrich_analysis<-enricher(KEGG_KO_enrichment_gene,TERM2GENE =KEGG_KO_background_genes,TERM2NAME = KEGG_KO_description,pvalueCutoff = 0.01, qvalueCutoff = 0.05)
KEGG_KO_enrich_analysis_data<-as.data.frame(KEGG_KO_enrich_analysis)
#绘制气泡图
dotplot(KEGG_KO_enrich_analysis,showCategory = 10)
pdf("KEGG_KO_enrich_dotplot.pdf",width = 10,height = 10)
dotplot(KEGG_KO_enrich_analysis,showCategory = 10)
dev.off()
#绘制条形图
barplot(KEGG_KO_enrich_analysis,showCategory = 10)
pdf("KEGG_KO_enrich_barplot.pdf",width = 10,height = 10)
barplot(KEGG_KO_enrich_analysis,showCategory = 10)
dev.off()
#导出K num富集结果
write.table(KEGG_KO_enrich_analysis_data,"KEGG_KO_enrichment_results.txt",col.names = T,row.names = F,sep = "\t")

跑完之后就会得到一些结果:


image.png

生成一些简单的气泡图,条形图,GO二级分类图


image.png

image.png

image.png
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 200,302评论 5 470
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,232评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 147,337评论 0 332
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,977评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,920评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,194评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,638评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,319评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,455评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,379评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,426评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,106评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,696评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,786评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,996评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,467评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,043评论 2 341

推荐阅读更多精彩内容