ClusterProfiler通路富集教程

单个基因水平上能反映的生物学信息有限,很多时候要进行通路富集分析,来从系统水平上反映出一组基因与哪些生物学通路相关。

假如你有一组感兴趣的基因,如何进行通路富集?

假如你感兴趣的这组基因是来源于特定芯片panel(而非全基因组),如何设置合适的背景基因集合进行通路富集分析?

假如你想比较两组感兴趣基因集合的通路异同又应如何操作呢?

……

今天我们就来谈谈Y叔开发的ClusterProfiler包做通路富集时的应用场景和详细步骤(以Homo sapiens为例)。

Step1.文件准备(如图1)。

1. 将所有文件存放于本地路径中(如D:/Bioming)。

2. 感兴趣的基因集合文件,每个基因为一行。(感兴趣的基因集合可以是差异表达基因、差异甲基化基因、突变基因集合等)。文件格式如图2。

3. 若感兴趣的基因集合是基于特定的panel芯片得到的,而并非全基因组数据,则需要准备背景基因集合文件,文件中包含了panel中所覆盖的所有基因,格式同图2。

4. 若需要将两组感兴趣的基因集合进行通路比较(如原发灶突变基因与转移灶突变基因所富集通路比较等),需要准备另一组文件感兴趣的基因集合2。

图1

图2

 

Step2.R语言执行通路富集分析。

 

进入到R界面,若未安装过ClusterProfiler包,用如下命令安装:

if (!requireNamespace(“BiocManager”, quietly = TRUE))

install.packages(“BiocManager”)

BiocManager::install(“clusterProfiler”, version = “3.8”)

加载ClusterProfiler包:

library(org.Hs.eg.db)

library(clusterProfiler)

分析开始:

setwd(“D:\ Bioming”)  #设置默认路径,即存放文件的本地路径,需要注意的是在windows系统中路径层级的符号为“”,输入到R中时需要改成双反斜杠”\”或正斜杠”/”。

CaseGene=read.table(“Case.txt”,header=F)[,1]  #读取感兴趣基因集合文件,第一列是基因名。

CaseGeneSet=bitr(CaseGene, ‘SYMBOL’, “ENTREZID”, “org.Hs.eg.db”)[, “ENTREZID”]  #将Gene名字转化为基因ID,若输入的是基因SYMBOL需要用该命令转化成基因ID,若输入的是基因ID,则忽略该命令。

Case_KEGG <- enrichKEGG(CaseGeneSet, keyType = “kegg”,pvalueCutoff=0.05,qvalueCutoff=0.05,pAdjustMethod = “BH”, minGSSize = 5, maxGSSize = 500,organism = “hsa”, use_internal_data=FALSE)  #KEGG富集分析

###Parameter解析

#CaseGeneSet: 感兴趣基因集合

#keyType: 设置类型,可选择的有”kegg”, ‘ncbi-geneid’, ‘ncib-proteinid’ and ‘uniprot’

#pvalueCutoff:设置p值的阈值

#qvalueCutoff:设置q值阈值(既校正后的阈值)

# pAdjustMethod:对p值进行校正的方法,可选方法有”holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”

#minGSSize: 设置最小基因个数

#maxGSSize:设置最大基因个数

#orgnaism:物种选择,可参考:http://www.genome.jp/kegg/catalog/org_list.html

# use_internal_data:设置注释参考的数据库,TRUE—参考最新的线上KEGG数据库,FALSE—参考KEGG.db数据库。

运行完上述命令后得到的Case_KEGG结果如下图。

 

接下来进行富集结果存储。

运行命令:

write.csv(summary(Case_KEGG),”KEGG-enrich.csv”,row.names =F)#导出Case_KEGG结果至默认路径下。

 

Step3.可视化通路富集结果。

运行命令:

barplot(Case_KEGG, showCategory=15,title=”Barplot of Case KEGG Enrichment”)#绘制条形图

 

运行命令:

dotplot(Case_KEGG,title=” Dotplot of Case KEGG Enrichment “)#绘制气泡图

 

运行命令:

enrichMap(Case_KEGG, layout=igraph::layout.kamada.kawai)#网络图,每个节点是一个富集到的pathway,若pathways之间有重叠的感兴趣基因,则自动将这两个通路用线连接。默认画top30个富集到的pathways, 节点大小对应该pathway下富集到的感兴趣基因个数,节点的颜色对应p.adjust的值,从小到大,对应蓝色到红色。

 

 

到此,已经可以完整呈现出感兴趣基因的通路富集分析了。那么如果这些感兴趣的基因并非在全基因组范围内找到的,而且来源于特定的芯片,如何科学地做出通路富集分析呢?

很简单,只需要将芯片所覆盖的所有基因作为背景基因,存放于文件中(背景基因集合)。

运行命令:

BackGround=read.table(“BackGroundGene.txt”,header=F)[,1]#读取背景基因,第一列是基因名。

BackGroundSet=bitr(BackGround,  ‘SYMBOL’, “ENTREZID”, “org.Hs.eg.db”)[, “ENTREZID”]#将Gene名字转化为基因ID,若输入的是基因ID,则忽略该命令。

Case_KEGG <- enrichKEGG(CaseGeneSet, keyType = “kegg”,pvalueCutoff=0.05,qvalueCutoff=0.05,pAdjustMethod = “BH”, minGSSize = 5, maxGSSize = 500,organism = “hsa”, use_internal_data=FALSE, universe=BackGroundSet)# 只需要在Step2的enrichKEGG这一步设置universe参数。

假如你想比较两组感兴趣基因集合的通路异同又应如何操作呢?需准备另一组感兴趣的基因集合(如图1中的感兴趣基因集合2,格式如图2)。导入第二组感兴趣的基因集合,并进行通路富集分析。

运行命令:

ControlGene=read.table(“Control.txt”,header=F)[,1]

ControlGeneSet=bitr(ControlGene, ‘SYMBOL’, “ENTREZID”, “org.Hs.eg.db”)[, “ENTREZID”]

Control_KEGG <- enrichKEGG(ControlGeneSet, keyType = “kegg”,pvalueCutoff=0.05,qvalueCutoff=0.05,pAdjustMethod = “BH”, minGSSize = 5, maxGSSize = 500,organism = “hsa”)

cp = list(Case=CaseGeneSet, Control=ControlGeneSet)  # 合并两个数据集,并转换为列表

compare <- compareCluster(cp, fun=”enrichKEGG”, organism=”hsa”, pvalueCutoff=0.05,pAdjustMethod = “BH”,qvalueCutoff = 0.05)#比较两组通路富集结果。

dotplot(compare, showCategory=15,font.size = 9,includeAll=TRUE)#画两组集合的通路比较图

 

今天的基因通路富集就到此结束了,不知道你学会了么?欢迎留言讨论!

生物信息学

单细胞测序-如何选择适合的scRNA测序方法

2020-8-28 1:10:47

生物信息学

免疫浸润分析方法

2020-8-28 1:17:18

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索