2021-03-21

功能富集分析（R语言）

一、clusterProfiler包做功能分析

1、包的安装与下载

（1）一般：install.packages("clusterProfiler")

library("clusterProfiler")

（2）package_list = c("clusterProfiler")
for(p in package_list){
  if (!requireNamespace(p, quietly = TRUE))
    BiocManager::install(p)
}

library("clusterProfiler")

（3）通过BioLite或者BioManager进行安装；或者说直接去github下载安装包到本地；甚至如果是因为R版本的问题，就得去R官网下载一个合适的R版本，然后再在Rstudio里将版本调整为新下的R版本。点击tools——global options进去就是了。

2、GO功能富集

（1）gene名的转换： geneList<-bitr(x, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db")

——对于geneID, GO功能富集分析默认是用的ENTREZID，所以一般转了ID会比较好，但是实际上不转也行。不过后续的enrichGO函数的参数keyType就要改一下。

（2）GO<-enrichGO(gene, OrgDb,  keyType = "ENTREZID", ont = "MF", pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2, minGSSize = 10, maxGSSize = 500,readable = FALSE,pool = FALSE)

——gene ：就是需要基因名的一个向量。默认是ENTREZID，但是也可以是SYMBOl。

——OrgDb ：org.Hs.eg.db，这个需要安装org.Hs.eg.db这个包。

——keyType：这个就是刚才说的那个gene名方式的选择。没记错的话应该是有ENTREZID、SYMBOL、UniProt这三种。

——ont ：BP、MF、CC、ALL选一个，ALL即同时做GO中的BP、MF、CC这三个。

——universe：是背景基因集，一般都可以不用管这个参数。

——minGSSize &maxGSSize：是富集的最小/大的基因集的大小（基因数目），default一般min为10，max为500。

（3）对于GO富集分析结果的可视化：

A、一般正常的barplot就可以了，如果你需要限制显示的富集通路的条数，可以设置。如：

 得到的结果图如上，这是柱状图。

B、 ggplot也可以，如 ggplot(data=GO)+  geom_bar(aes(x=Description,y=-log10(pvalue), fill=ID), stat='identity') + coord_flip() + scale_x_discrete(limits=GO$Description) 。

并且用ggplot函数画图可以实现不同富集类型的一个区别，即同时画BP、MF、CC，只需要设置不同的填充颜色即可。不过可能显示效果不是很好，因为一起进行功能分析，可能效果不是很好。

C、还有dotplot、cnetplot、enrichMap、plotGOgraph等函数可以实现对GO功能富集分析的一个可视化。具体参数可以通过？functionname进行查看。

D、另外，对于GO，还有一个专门的R包：GOplot包。它主要有

——GOBar（subset(circ,category=="BP")）,另外这个函数还可以结合ggplot2进行可视化 。下图是从网上复制的。

GOBar(subset(circ,category=="BP"))+
  theme(axis.text.x = element_text(size=5))

——GOBubble：气泡图分布，主要有GO_result，title，color，display，labels这几个参数。可以同时展示BP、MF、CC的分布，用不同颜色可以区分开。

——GOChord(chord,space=  , gene.oder ='logFC',gene.space =   , gene.size =   )；这些参数可以自己设置，直至弦图是你想要的效果。

——GOHeat : 画GO功能分析结果的热图。GOHeat ( chord[,-8], nlfc=0 ）

——GOCluster : 画功能富集分析的聚类图，一般有EC$process, clust.by='logFC' , term.width这些参数。

——GOVeen ：画韦恩图，不过这个函数一般需要自己选定富集的通路进行韦恩图的绘制。

3、KEGG通路富集分析

（1）首先也是gene列表/向量，KEGG通路富集需要用的是ENTREZID。如果不是的话，需要用bitr函数进行ID的转换。有两种方式。

——bitr_kegg(names(geneList), fromType = 'ncbi-geneid', toType="kegg", organism="hsa")

——bitr(x, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db")

（2）enrichKEGG(gene,organism = "hsa",keyType = "kegg", pvalueCutoff = 0.05,pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500,qvalueCutoff = 0.2, use_internal_data = FALSE)

——每个参数与enrichGO函数都差不多，universe默认也不需要设置。就是合理keyType需要设置为kegg，物种是hsa。别的都没有什么需要注意的。

（3）KEGG通路富集的可视化。

——barplot()

—— dotplot() ：可以看富集到的通路的基因的情况。

——如果想直接看富集到的通路图，可用browseKEGG( kegg, 'hsa04974'），打开浏览器调到KEGG数据库的这条通路上，并且富集的基因会以不同的颜色标出。

——也可以通过pathview来看。

——另外，cnetplot()、heatplot（） 这些函数都可以对KEGG富集结果进行可视化。

4、DO（disease ontology）分析

（1）enrichDO(gene, ont = "DO", pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500, qvalueCutoff = 0.2, readable = FALSE)

具体的和GO & KEGG都差不多。

5、DAVID分析（这个其实可以通过在线的DAVID软件实现）

（1）enrichDAVID(gene, idType = "ENTREZ_GENE_ID", universe, minGSSize = 10, maxGSSize = 500, annotation = "GOTERM_BP_FAT", pvalueCutoff = 0.05,pAdjustMethod = "BH"，  qvalueCutoff = 0.2,species = NA,  david.user)

 

二、GSEA（基因集富集分析） ——用fgsea包做

1、fgsea包做功能富集分析主要是通过基因的FC值来做的。

2、fgsea ( pathways, stats, nperm, minSize = 1, maxSize = Inf, nproc = 0,  gseaParam = 1, BPPARAM = NULL)

3、gsea结果的可视化：gseaplot(  x,geneSetID,  by = "all",  title = "",  color = "black",  color.line = "green",  color.vline = "#FA5860",  ...)

——x ：为fgsea得到的结果。

——geneSetID ：是需要做基因集富集分析的基因集的ID。

——by : one of "runningScore" or "position"，或者都要。