R|clusterProfiler-富集分析

简单总结clusterProfiler包进行GO、KEGG的富集分析方法，结果输出及内置的图形展示。

一 bioconductor包安装

#PS：安装是否顺利，看缘分了（用本站weinformatics则不用担心这些问题，包都装好了，解决了你入门的最大问题）

#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 ENTREZID
1    AASDH   132949
2   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 ENTREZID
1    AASDH ENSG00000157426   132949
2   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 SYMBOL
1     4312 ENSG00000196611   MMP1
2     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                                                Description
GO:0002887       BP GO:0002887 negative regulation of myeloid leukocyte mediated immunity
GO:0033004       BP GO:0033004                negative regulation of mast cell activation
           GeneRatio  BgRatio      pvalue  p.adjust    qvalue                   geneID Count
GO:0002887     2/121 10/11461 0.004706555 0.7796682 0.7796682              CD300A/CD84     2
GO: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-Value
qvalue：对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  BgRatio
hsa04750 hsa04750 Inflammatory mediator regulation of TRP channels      5/53  97/7387
hsa04020 hsa04020                        Calcium signaling pathway      6/53 182/7387
               pvalue   p.adjust    qvalue                              geneID Count
hsa04750 0.0006135305 0.08589427 0.0807277             40/3556/3708/5608/79054     5
hsa04020 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.gz
gda <- read.delim("all_gene_disease_associations.tsv")
head(gda,1)

 geneId geneSymbol diseaseId    diseaseName       score NofPmids NofSnps source
1      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.adjust
C3642345 C3642345 Luminal A Breast Carcinoma    11/198 110/17074 6.057328e-08 0.0001488891
               qvalue                                              geneID Count
C3642345 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等相关分析。

请关注“恒诺新知”微信公众号，感谢“R语言“，”数据那些事儿“，”老俊俊的生信笔记“，”冷🈚️思“，“珞珈R”，“生信星球”的支持！