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、背景基因:指的是前景基因在某个基因集合进行富集,这个基因集合就是背景基因
3、描述信息:每个GO的Term的属性,或者是每个KO号或者map号的属性。
我们具备前景基因,背景基因以及描述信息我们就可以做富集分析啦。
二、文件准备
1、前景基因:这是必须的啦。有时候需要进行ID转换,但是个人觉得ID转换根据需要来就行。如果前景基因里面的基因ID是包括在背景基因里面,那就需要进行转换。如果前景基因在是新的基因或者在背景基因没有被注释到的,就不用进行ID转换。下面这个就是融合基因,在背景基因里面没有注释到的,那么我就不要转换。
2、背景基因:一个基因可能具备多个GO term,一个基因也可能参与多个通路,与之相对应的有多个map号
这个案例中背景基因文件构建思路如下图
Knum是与组成成分有关的,就好像pathway是表示通路的一样。由于我们前景基因是融合基因的,那么我们就需要将前景基因与对应GO号、map号,k num加入到对应的背景基因文件中去,构成一个新的背景文件。
3、描述文件
上面的文件是有固定格式的。
三、富集分析
准备好这些文件之后进开始进行富集分析啦,使用下面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")
跑完之后就会得到一些结果:
生成一些简单的气泡图,条形图,GO二级分类图