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

GO、KEGG富集分析是我們做生信分析較為常用的部分，它可以將基因與功能相聯(lián)系起來。
GO指的是Gene Ontology，是基因功能國際標準分類體系。目的在于建立一個適用于各種物種的，對基因和蛋白質功能進行限定和描述的，并能隨著研究不斷深入而更新的語言詞匯標準。GO分為分子功能（Molecular Function)(MF)、生物過程（Biological Process）(BP)、和細胞組成（Cellular Component）(CC)三個部分。
KEGG指的是京都基因與基因組百科全書，通常我們使用KEGG中的pathway模塊，將基因映射到某些通路上，了解基因參與生物體中的代謝過程等。
對于模式生物，GO和KEGG富集分析實現(xiàn)起來比較容易，對于非模式生物來說還是需要花點時間和精力。對于模式生物的GO和KEGG富集分析，網(wǎng)上教程案例挺多的。對于非模式生物，以小麥為例，進行下面一些基本的富集分析。

一、基本概念

做富集分析，我們需要了解一下幾個概念。
1、前景基因：指的是我們所要進行富集的基因，一般是基因的ID
2、背景基因：指的是前景基因在某個基因集合進行富集，這個基因集合就是背景基因


v2-50cbf96d15658a93c6d680a1c9436eb7_720w.png

3、描述信息：每個GO的Term的屬性，或者是每個KO號或者map號的屬性。


image.png

image.png

我們具備前景基因，背景基因以及描述信息我們就可以做富集分析啦。

二、文件準備

1、前景基因:這是必須的啦。有時候需要進行ID轉換，但是個人覺得ID轉換根據(jù)需要來就行。如果前景基因里面的基因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包進行富集分析，繪制一些簡單的圖。需要繪制高級圖話，根據(jù)結果進行繪制。下面腳本里面參數(shù)根據(jù)自己需要可進行修改。 
library("clusterProfiler")
library("ggplot2")
#導入數(shù)據(jù)
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),]
#提取出各二級分類數(shù)據(jù)
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)
#合并上面三個數(shù)據(jù)集
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")+#設置數(shù)值變量填充顏色，若要設置顏色可用scale_color_gradient函數(shù)
  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")+#設置數(shù)值變量填充顏色，若要設置顏色可用scale_color_gradient函數(shù)
  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富集分析
#導入數(shù)據(jù)
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
最后編輯于 ：
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容