使用clusterProfiler對非模式生物進行富集分析

來自專欄小王子的生信筆記

最近，小編有很多同學問我，非模式生物如何做富集分析？

小編本身是做小麥的，也屬於非模式生物的範疇。以前的話，非模式生物要用blast2go跑電子注釋，而blast2go又需要使用MySQL，沒有root許可權的話非常的麻煩。所以非模式生物如何做富集分析也困擾了小編很久，直到一天，小編看了公眾號「小麥研究聯盟」的一篇推送，發現了Y叔的神包「clusterProfiler 」！可以輕鬆做富集分析！

非模式生物的話，分為兩種，一種是可以在AnnotationHub上在線抓取Org.Db的非模式生物，一種是在AnnotationHub上沒有Org.Db的生物。

下面我們先來講講可以在AnnotationHub上抓取到Org.Db的非模式生物如何做富集分析：

# 載入包library("AnnotationHub")library("biomaRt")library("clusterProfiler")library("topGO")library("Rgraphviz")library("pathview")# 導入文件data <- read.table(file,header = F,sep=" ")# 搜尋 OrgDb（比如我想查找物種 Theobroma cacao）hub <- AnnotationHub::AnnotationHub()query(hub,"Theobroma cacao")

我們會搜索到Org.Db「AH59191」。

# 抓取 OrgDbThecacao.OrgDb <- hub[["AH59191"]]# 查看 OrgDb 里的Gene ID類型columns(Thecacao.OrgDb)# 查看Gene ID 類型columns(Thecacao.OrgDb)# ID轉換（將我們的ID轉換成與Org.Db一致的ID類型）data <- as.character(data$V1)data_id <- mapIds(x = Thecacao.OrgDb,keys = data,keytype = "SYMBOL",column = "ENTREZID")# 去除NA值na.omit(data_id)# 富集（其他富集函數參見R包文檔）erich.go.BP <- enrichGO(gene=data,OrgDb = Thecacao.OrgDb,keyType = "SYMBOL",ont = "BP",pvalueCutoff = 0.01,qvalueCutoff = 0.05,readable = T)# 繪製條形圖barplot(erich.go.BP)# 繪製氣泡圖dotplot(erich.go.BP)

以上部分，就是可以抓取到Org.Db的物種進行富集分析的步驟。

但是，對於抓取不到的物種，我們也有方式做富集分析。以小麥GO富集分析為例：

首先在IWGSC的官網上下載最新的注釋文件：https://wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations

得到注釋文件後，我們就可以提取到每個基因對應的GO號，這時，我們可以準備三個文件。

第一個文件：所要富集基因的列表（格式如下）

第二個文件，基因與GO號的對應關係（格式如下）.

小編之前整理的是csv文件，所以這裡是逗號分隔，大家也可以用TAB分割，但是兩列的順序一定不要調換！

第三個文件，GO號與其對應注釋的文件（格式如下）

有了這三個文件，我們就可以開始富集分析啦！

library("clusterProfiler")# 導入基因列表gene <- read.csv(file,header = T,sep=",")gene <- as.factor(gene$V1)# 導入注釋文件term2gene <- read.csv(file,header=F,sep=",")term2name <- read.csv(file,header=F,sep=",")# 富集分析x <- enricher(gene,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.01, pAdjustMethod = "BH", qvalueCutoff = 0.05)# 設置結果輸出文件ouf <- paste(out_file,sep =" ")# 輸出結果write.csv(x,ouf)# 繪製條形圖barplot(x)# 繪製氣泡圖dotplot(x)

有時候，我們想做兩列基因的富集結果比較時，可以用以下方式展現：

首先，我們要生成需要進行比較的兩個基因列表（格式如下）

# 導入文件gene <- read.csv(file,header = T,sep=",")# 富集分析x <- compareCluster(gene, fun=enricher,TERM2GENE=term2gene,TERM2NAME=term2name,pvalueCutoff = 0.01, pAdjustMethod = "BH", qvalueCutoff = 0.05)# 繪製氣泡圖dotplot(x, showCategory=10,includeAll=TRUE)

這樣，富集分析就簡單完成啦！

當然，還有很多函數和細節，小編沒有細講，有興趣的朋友可以查看「clusterProfiler」的文檔：

Statistical analysis and visualization of functional profiles for genes and gene clusters

參考資料：

Statistical analysis and visualization of functional profiles for genes and gene clusters

非模式基因GO富集分析：以玉米為例+使用OrgDb

基因的富集分析 - CSDN博客

2018-01-07 生物信息學教程-GO, KEGG, DO富集分析_嗶哩嗶哩 (゜-゜)つロ 乾杯~-bilibili

另外！推薦微信公眾號「小麥研究聯盟」，做小麥的同學們可以關注！裡面的內容都是精華！

http://weixin.qq.com/r/lj-r8_bEELtvraKQ92qE (二維碼自動識別)

歡迎關注小編的微信公眾號「生信小王子」！