R|clusterProfiler-富集分析

04-26

簡單總結clusterProfiler包進行GO、KEGG的富集分析方法，結果輸出及內置的圖形展示。

一 bioconductor包安裝

#PS：安裝是否順利，看緣分了

#source("https://bioconductor.org/biocLite.R")#biocLite("DOSE") #biocLite("topGO")#biocLite("clusterProfiler")#biocLite("pathview")

##載入需要的R包

library(DOSE)library(org.Hs.eg.db)library(topGO)library(clusterProfiler)library(pathview)

二、數據載入及轉化

需要將輸入的基因格式轉為clusterProfiler可分析的格式，通過功能函數bitr實現各種ID之間的轉化。通過keytypes函數可查看所有支持及可轉化類型

keytypes(org.Hs.eg.db)

1.向量方式讀入#適合少量基因

x <- c("AASDH","ABCB11","ADAM12","ADAMTS16","ADAMTS18")test = bitr(x, #數據集fromType="SYMBOL", #輸入為SYMBOL格式toType="ENTREZID", # 轉為ENTERZID格式OrgDb="org.Hs.eg.db") #人類資料庫head(test,2) SYMBOL ENTREZID1 AASDH 1329492 ABCB11 8647

2.基因list文件讀入

data <- read.table("gene",header=FALSE) #單列基因名文件data$V1 <- as.character(data$V1) #需要character格式，然後進行ID轉化#將SYMBOL格式轉為ENSEMBL和ENTERZID格式 test1 = bitr(data$V1, fromType="SYMBOL", toType=c("ENSEMBL", "ENTREZID"), OrgDb="org.Hs.eg.db")head(test1,2) SYMBOL ENSEMBL ENTREZID1 AASDH ENSG00000157426 1329492 ABCB11 ENSG00000073734 8647

3. 內置示例數據

data(geneList, package="DOSE") #富集分析的背景基因集gene <- names(geneList)[abs(geneList) > 2]gene.df <- bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db)head(gene.df,2) ENTREZID ENSEMBL SYMBOL1 4312 ENSG00000196611 MMP12 8318 ENSG00000093009 CDC45

三、 GO分析

3.1 groupGO 富集分析

ggo <- groupGO(gene = test1$ENTREZID, OrgDb = org.Hs.eg.db, ont = "CC",level = 3,readable = TRUE)

3.2 enrichGO 富集分析

ego_ALL <- enrichGO(gene = test1$ENTREZID, universe = names(geneList), #背景基因集 OrgDb = org.Hs.eg.db, #沒有organism="human"，改為OrgDb=org.Hs.eg.db #keytype = ENSEMBL, ont = "ALL", #也可以是 CC BP MF中的一種 pAdjustMethod = "BH", #矯正方式 holm」, 「hochberg」, 「hommel」, 「bonferroni」, 「BH」, 「BY」, 「fdr」, 「none」中的一種 pvalueCutoff = 1, #P值會過濾掉很多，可以全部輸出 qvalueCutoff = 1, readable = TRUE) #Gene ID 轉成gene Symbol ，易讀head(ego_ALL,2) ONTOLOGY ID DescriptionGO:0002887 BP GO:0002887 negative regulation of myeloid leukocyte mediated immunityGO:0033004 BP GO:0033004 negative regulation of mast cell activation GeneRatio BgRatio pvalue p.adjust qvalue geneID CountGO:0002887 2/121 10/11461 0.004706555 0.7796682 0.7796682 CD300A/CD84 2GO:0033004 2/121 10/11461 0.004706555 0.7796682 0.7796682 CD300A/CD84 2其中：ONTOLOGY：CC BP MF GO ID: Gene Ontology資料庫中唯一的標號信息Description ：Gene Ontology功能的描述信息GeneRatio：差異基因中與該Term相關的基因數與整個差異基因總數的比值BgRation：所有（ bg）基因中與該Term相關的基因數與所有（ bg）基因的比值pvalue: 富集分析統計學顯著水平，一般情況下， P-value < 0.05 該功能為富集項p.adjust 矯正後的P-Valueqvalue：對p值進行統計學檢驗的q值geneID：與該Term相關的基因Count：與該Term相關的基因數

3.2.1 setReadable函數進行轉化

ego_MF <- enrichGO(gene = test1$ENTREZID, universe = names(geneList),OrgDb = org.Hs.eg.db,ont = "MF", pAdjustMethod = "BH",pvalueCutoff = 1,qvalueCutoff = 1,readable = FALSE)ego_MF1 <- setReadable(ego_MF, OrgDb = org.Hs.eg.db)3.3 GSEA分析（暫略）

3.4 輸出結果及圖像

3.4.1 輸出結果

write.csv(summary(ego_ALL),"ALL-enrich.csv",row.names =FALSE)

3.4.2 繪製圖形

##可視化--點圖dotplot(ego_MF,title="EnrichmentGO_MF_dot")#點圖，按富集的數從大到小的##可視化--條形圖barplot(ego_MF, showCategory=20,title="EnrichmentGO_MF")#條狀圖，按p從小到大排，繪製前20個Term

##可視化--plotGOgraph(ego_MF)

3.5 過濾

3.5.1 利用內置cutoff設置閾值，by設定指標

#可能會有問題，還不太清楚ego_MF.fil <- simplify(ego_MF, cutoff=0.05, by="pvalue", select_fun=min)

3.5.2 存儲結果然後過濾

ego_ALL.sig <- ego_ALL[ego_ALL$pvalue <= 0.01]

過濾後為數據框，不能用自帶的參數直接繪製，可以使用ggplot2進行繪製。（暫略）

四、KEGG分析

4.1 候選基因進行通路分析

kk <- enrichKEGG(gene = test1$ENTREZID, organism = hsa, #KEGG可以用organism = hsa pvalueCutoff = 1)head(kk,2) ID Description GeneRatio BgRatiohsa04750 hsa04750 Inflammatory mediator regulation of TRP channels 5/53 97/7387hsa04020 hsa04020 Calcium signaling pathway 6/53 182/7387 pvalue p.adjust qvalue geneID Counthsa04750 0.0006135305 0.08589427 0.0807277 40/3556/3708/5608/79054 5hsa04020 0.0018078040 0.12654628 0.1189345 493/1129/2066/3707/3708/4842 6

4.2 富集結果及圖形展示

4.2.1 結果輸出文件

write.csv(summary(kk),"KEGG-enrich.csv",row.names =FALSE)

4.2.1 KEGG 氣泡圖

dotplot(kk,title="Enrichment KEGG_dot")

4.2.2 查看特定通路圖

hsa04750 <- pathview(gene.data = geneList,

pathway.id = "hsa04750", #上述結果中的hsa04750通路

species = "hsa",

limit = list(gene=max(abs(geneList)), cpd=1))

五、注釋文件、注釋庫

如果clusterProfiler包沒有所需要物種的內置資料庫，可以通過自定義注釋文件或者自建注釋庫的方法進行富集分析。

5.1 自定義注釋文件

5.1.1 待富集的基因list

data(geneList, package="DOSE")deg <- names(geneList)[abs(geneList)>2]

5.1.2 讀入注釋文件：

## downloaded from http://www.disgenet.org/ds/DisGeNET/results/all_gene_disease_associations.tar.gzgda <- read.delim("all_gene_disease_associations.tsv")head(gda,1) geneId geneSymbol diseaseId diseaseName score NofPmids NofSnps source1 1 A1BG C0001418 Adenocarcinoma 0.002732912 1 0 LHGDN

自定義的兩個注釋文件（重要）

disease2gene=gda[, c("diseaseId", "geneId")]disease2name=gda[, c("diseaseId", "diseaseName")]

5.1.3 enricher函數進行富集分析

x = enricher(deg, TERM2GENE=disease2gene, TERM2NAME=disease2name)head(summary(x),1) ID Description GeneRatio BgRatio pvalue p.adjustC3642345 C3642345 Luminal A Breast Carcinoma 11/198 110/17074 6.057328e-08 0.0001488891 qvalue geneID CountC3642345 0.0001287342 9833/55355/3620/891/9122/2167/3169/5304/2625/9/5241 11

5.2 資料庫

5.2.1 查看當前支持的物種

Bioconductor - BiocViews

5.2.2 通過AnnotationHub自建OrgDb注釋庫並富集

library(AnnotationHub)hub <- AnnotationHub() #較慢query(hub, "Cricetulus")Cgriseus <- hub[["AH48061"]]sample_gene <- sample(keys(Cgriseus), 100)str(sample_gene)sample_test <- enrichGO(sample_gene, OrgDb=Cgriseus, pvalueCutoff=1, qvalueCutoff=1)head(summary(sample_test))

六、參考資料

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

以上，clusterProfiler包進行GO、KEGG的富集分析方法，結果輸出及內置的圖形展示。待補充富集結果ggplot2繪製，GSEA等相關分析。