【生信分析】clusterProfiler: universal enrichment tool for functional and comparative study（3）
clusterProfiler: universal enrichment tool for functional and comparative study
- ClusterProfiler: 用于功能和比较研究的通用富集工具
- 7. MSigDb 分析
- 8. Reactome通路分析
- 9. MeSH富集分析
- 10. 基因组协作的功能富集分析
- 11. 生物主题对比
- 11.1 compareCluster 的公式接口（Formula interface of compareCluster）
- 11.2 数据谱比较的可视化（Visualization of profile comparison）
ClusterProfiler: 用于功能和比较研究的通用富集工具
7. MSigDb 分析
MSigDB 是一个带注释的基因集的集合，它包括 8 个主要集合：
- H: hallmark gene sets （癌症特征基因集合）
- C1: positional gene sets （染色体位置基因集合）
- C2: curated gene sets （包含了已知数据框、文件和专家支持的基因集信息）
- C3: motif gene sets （调控靶基因集合，包括miRNA-target以及转录因子-target调控模式）
- C4: computational gene sets (计算机软件预测出来的基因集合，主要是和癌症相关的基因)
- C5: GO gene sets (Gene Ontology对应的基因集合)
- C6: oncogenic signatures （致癌基因集合）
- C7: immunologic signatures （免疫相关基因集）
用户可以使用enricher和 GSEA 函数来分析从分子特征数据库 (MSigDb) 下载的基因集集合。 clusterProfiler 提供了一个函数 read.gmt，用于将gmt文件解析为一个 TERM2GENE data.frame，该框架为enricher和 GSEA 函数做好了准备。
> data(geneList, package="DOSE")
> gene <- names(geneList)[abs(geneList) > 2]
> gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
> c5 <- read.gmt(gmtfile)
> egmt <- enricher(gene, TERM2GENE = c5)
> head(egmt)ID Description GeneRatio BgRatio pvalue p.adjust qvalue
SPINDLE SPINDLE SPINDLE 11/82 39/5270 7.667674e-12 5.214018e-10 4.197043e-10
MICROTUBULE_CYTOSKELETON MICROTUBULE_CYTOSKELETON MICROTUBULE_CYTOSKELETON 16/82 152/5270 8.449298e-10 2.872761e-08 2.312439e-08
CYTOSKELETAL_PART CYTOSKELETAL_PART CYTOSKELETAL_PART 15/82 235/5270 2.414879e-06 5.237096e-05 4.215619e-05
SPINDLE_MICROTUBULE SPINDLE_MICROTUBULE SPINDLE_MICROTUBULE 5/82 16/5270 3.080645e-06 5.237096e-05 4.215619e-05
MICROTUBULE MICROTUBULE MICROTUBULE 6/82 32/5270 7.740446e-06 1.052701e-04 8.473751e-05
CYTOSKELETON CYTOSKELETON CYTOSKELETON 16/82 367/5270 1.308357e-04 1.482805e-03 1.193589e-03geneID Count
SPINDLE 991/9493/9787/22974/983/332/3832/7272/9055/6790/24137 11
MICROTUBULE_CYTOSKELETON 991/9493/9133/7153/9787/22974/4751/983/332/3832/7272/9055/6790/24137/4137/7802 16
CYTOSKELETAL_PART 991/9493/7153/9787/22974/4751/983/332/3832/7272/9055/6790/24137/4137/7802 15
SPINDLE_MICROTUBULE 983/332/3832/9055/24137 5
MICROTUBULE 983/332/3832/9055/24137/4137 6
CYTOSKELETON 991/9493/9133/7153/9787/22974/4751/983/332/3832/7272/9055/6790/24137/4137/7802 16> egmt2 <- GSEA(geneList, TERM2GENE=c5, verbose=FALSE)
> head(egmt2)ID Description setSize enrichmentScore NES pvalue
EXTRACELLULAR_REGION EXTRACELLULAR_REGION EXTRACELLULAR_REGION 401 -0.3860230 -1.705466 0.001253133
EXTRACELLULAR_REGION_PART EXTRACELLULAR_REGION_PART EXTRACELLULAR_REGION_PART 310 -0.4101043 -1.778083 0.001302083
EXTRACELLULAR_MATRIX EXTRACELLULAR_MATRIX EXTRACELLULAR_MATRIX 95 -0.6229461 -2.306284 0.001503759
CELL_PROJECTION CELL_PROJECTION CELL_PROJECTION 87 -0.4729701 -1.722175 0.001526718
PROTEINACEOUS_EXTRACELLULAR_MATRIX PROTEINACEOUS_EXTRACELLULAR_MATRIX PROTEINACEOUS_EXTRACELLULAR_MATRIX 93 -0.6355317 -2.338288 0.001538462
EXTRACELLULAR_MATRIX_PART EXTRACELLULAR_MATRIX_PART EXTRACELLULAR_MATRIX_PART 54 -0.5908035 -1.969987 0.001587302p.adjust qvalues rank leading_edge
EXTRACELLULAR_REGION 0.027746 0.02167891 1797 tags=29%, list=14%, signal=26%
EXTRACELLULAR_REGION_PART 0.027746 0.02167891 1897 tags=32%, list=15%, signal=28%
EXTRACELLULAR_MATRIX 0.027746 0.02167891 1473 tags=48%, list=12%, signal=43%
CELL_PROJECTION 0.027746 0.02167891 2280 tags=28%, list=18%, signal=23%
PROTEINACEOUS_EXTRACELLULAR_MATRIX 0.027746 0.02167891 1473 tags=49%, list=12%, signal=44%
EXTRACELLULAR_MATRIX_PART 0.027746 0.02167891 1794 tags=59%, list=14%, signal=51%
8. Reactome通路分析
ReactomePA（Yu 和 He 2016）使用 Reactome 作为通路数据的来源。 ReactomePA 中enrichPathway 和gsePathway 的函数调用与enrichKEGG 和gseKEGG 一致。
9. MeSH富集分析
meshes 支持使用 MeSH 注释的基因列表或整个表达谱的富集分析（过表达分析和基因集富集分析）。 支持来自 gendoo、gene2pubmed和 RBBH 的数据源。 用户可以选择感兴趣的类别进行测试。 支持所有 16 个类别。 该分析支持在 [MeSHDb BiocView](https://bioconductor.org/packages/release/BiocViews.html#___MeSHDb) 中列出的超过 70 个物种。
有关算法的详细信息，请参阅 DOSE(Yu et al. 2015) 包。
library(meshes)
data(geneList, package="DOSE")
de <- names(geneList)[1:100]
x <- enrichMeSH(de, MeSHDb = "MeSH.Hsa.eg.db", database='gendoo', category = 'C') # 报错 object 'MeSH.Hsa.eg.db' not found
head(x)
在过表达分析中，我们使用来自 gendoo 和 C（疾病）类别的数据源。
在以下示例中，我们使用来自gene2pubmed 的数据源并使用GSEA测试类别G（现象和过程）。
y <- gseMeSH(geneList, MeSHDb = "MeSH.Hsa.eg.db", database = 'gene2pubmed', category = "G") # 报错，同上
head(y)
用户可以使用在enrichplot中实现的可视化方法（即barplot、dotplot、cnetplot、emapplot 和 gseaplot）来可视化这些富集结果。 使用这些可视化方法，可以更轻松地解释丰富的结果。
dotplot(x)
gseaplot(y, y[1,1], title=y[1,2])
10. 基因组协作的功能富集分析
使用 NGS 数据（例如，RNA-Seq 和 ChIP-Seq）的功能分析可以通过 ChIPseeker（Yu、Wang 和 He 2015）包将编码和非编码区与编码基因连接起来进行，该包可以将基因组区域注释到它们的 最近基因、宿主基因和侧翼基因。 此外，它提供了一个功能，seq2gene，它同时考虑宿主基因、启动子区域和来自基因间区域的侧翼基因，这些基因可能通过顺式调节受到控制。 该功能以多对多的方式将基因组区域映射到基因并促进功能分析。 有关更多详细信息，请参阅 ChIPseeker（Yu、Wang 和 He 2015）。
11. 生物主题对比
clusterProfiler 是为生物主题比较而开发的（Yu et al. 2012），它提供了一个函数compareCluster来自动计算每个基因簇已富集的功能类别。
> data(gcSample)
> lapply(gcSample, head)
$X1
[1] "4597" "7111" "5266" "2175" "755" "23046"$X2
[1] "23450" "5160" "7126" "26118" "8452" "3675" $X3
[1] "894" "7057" "22906" "3339" "10449" "6566" $X4
[1] "5573" "7453" "5245" "23450" "6500" "4926" $X5
[1] "5982" "7318" "6352" "2101" "8882" "7803"$X6
[1] "5337" "9295" "4035" "811" "23365" "4629" $X7
[1] "2621" "2665" "5690" "3608" "3550" "533" $X8
[1] "2665" "4735" "1327" "3192" "5573" "9528"
geneCluster 参数的输入应该是基因 ID 的命名列表。 为了加快本文档的编译速度，我们设置 use_internal_data = TRUE。
> ck <- compareCluster(geneClusters = gcSample, fun = "enrichKEGG")
> head(as.data.frame(ck))Cluster ID Description GeneRatio BgRatio pvalue p.adjust qvalue
1 X2 hsa04110 Cell cycle 18/384 124/8096 1.983449e-05 0.00604952 0.00565805
2 X2 hsa05169 Epstein-Barr virus infection 23/384 202/8096 8.044708e-05 0.01226818 0.01147429
3 X2 hsa05340 Primary immunodeficiency 8/384 38/8096 3.322186e-04 0.03377555 0.03158991
4 X3 hsa04512 ECM-receptor interaction 9/187 88/8096 1.830703e-04 0.04554147 0.04255205
5 X3 hsa04060 Cytokine-cytokine receptor interaction 17/187 295/8096 4.547531e-04 0.04554147 0.04255205
6 X3 hsa04151 PI3K-Akt signaling pathway 19/187 354/8096 5.117019e-04 0.04554147 0.04255205geneID Count
1 991/1869/890/1871/701/990/10926/9088/8317/9700/9134/1029/2810/699/11200/23594/8555/4173 18
2 4067/3383/7128/1869/890/1871/578/864/637/9641/6891/355/9134/5971/916/956/6850/7187/3551/919/4734/958/6772 23
3 100/6891/3932/973/916/925/958/64421 8
4 7057/3339/1299/3695/1101/3679/3910/3696/3693 9
5 2919/4982/3977/6375/8200/608/8792/3568/2057/1438/8718/655/652/10220/50615/51561/7042 17
6 894/7057/6794/2247/1299/3695/2252/2066/1101/8817/1021/5105/3679/3082/2057/3910/3551/3696/3693 19
11.1 compareCluster 的公式接口（Formula interface of compareCluster）
compareCluster支持以Entrez~group或者Entrez~group+othergroup的类型的公式（支持公式的代码由 Giovanni Dall’Olio 贡献）。
> mydf <- data.frame(Entrez=names(geneList), FC=geneList)
> mydf <- mydf[abs(mydf$FC) > 1,]
> mydf$group <- "upregulated"
> mydf$group[mydf$FC < 0] <- "downregulated"
> mydf$othergroup <- "A"
> mydf$othergroup[abs(mydf$FC) > 2] <- "B"
> formula_res <- compareCluster(Entrez~group+othergroup, data=mydf, fun="enrichKEGG")
> head(as.data.frame(formula_res))Cluster group othergroup ID Description GeneRatio BgRatio pvalue p.adjust qvalue
1 downregulated.A downregulated A hsa04974 Protein digestion and absorption 16/319 103/8096 2.311850e-06 6.450063e-04 6.132487e-04
2 downregulated.A downregulated A hsa04510 Focal adhesion 20/319 201/8096 1.214789e-04 1.694631e-02 1.611194e-02
3 downregulated.A downregulated A hsa04151 PI3K-Akt signaling pathway 28/319 354/8096 3.218931e-04 2.993606e-02 2.846212e-02
4 downregulated.A downregulated A hsa04512 ECM-receptor interaction 11/319 88/8096 6.348341e-04 4.427968e-02 4.209952e-02
5 downregulated.B downregulated B hsa03320 PPAR signaling pathway 5/43 76/8096 4.650631e-05 7.254985e-03 6.804608e-03
6 upregulated.A upregulated A hsa04110 Cell cycle 20/220 124/8096 1.009671e-10 2.524176e-08 2.263788e-08geneID Count
1 1281/50509/1290/477/1294/1360/1289/1292/1296/23428/1359/1300/1287/6505/2006/7373 16
2 55742/2317/7058/25759/56034/3693/3480/5159/857/1292/3908/3909/63923/3913/1287/3679/7060/3479/10451/80310 20
3 55970/5618/7058/10161/56034/3693/4254/3480/4908/5159/1292/3908/2690/3909/8817/9223/4915/3551/2791/63923/3913/9863/3667/1287/3679/7060/3479/80310 28
4 7058/3693/3339/1292/3908/3909/63923/3913/1287/3679/7060 11
5 9370/5105/2167/3158/5346 5
6 4171/993/990/5347/701/9700/898/23594/4998/9134/4175/4173/10926/6502/994/699/4609/5111/1869/1029 20
11.2 数据谱比较的可视化（Visualization of profile comparison）
我们可以使用 dotplot 方法将结果可视化。
> dotplot(ck)
> dotplot(formula_res)
> dotplot(formula_res, x=~group) + ggplot2::facet_grid(~othergroup)
默认情况下，只绘制每个集群的前 5 个（最重要的）类别。 用户可以更改参数 showCategory 以指定要绘制的每个集群的类别数，如果将 showCategory 设置为 NULL，则将绘制整个结果。
plot 函数接受一个参数 by 用于设置点大小的比例。 默认参数by设置为*“geneRatio”，对应输出的“GeneRatio”列。 如果设置为count*，则比较将基于基因计数，而如果设置为 rowPercentage，则点大小将按count/(sumofeachrow)count/(sum of each row)count/(sumofeachrow)进行归一化。
为了提供完整的信息，我们还提供了当 by 设置为 rowPercentage 时每个类别中已识别基因的数量（括号中的数字）和每个聚类标签中的基因簇数（括号中的数字）当 by 设置为geneRatio 时，如图所示 在图 3 中。如果点大小基于count，则不会显示行号。
p 值表明哪些类别更有可能具有生物学意义。 图中的点根据其相应的 p 值进行颜色编码。 从红色到蓝色的颜色渐变对应于 p 值递增的顺序。 也就是说，红色表示低 p 值（高富集），蓝色表示高 p 值（低富集）。 通过参数 pvalueCutoff 给出的阈值过滤掉 P 值和调整后的 p 值，并且可以通过 qvalue 估计 FDR。
用户可以参考 Yu (2012)(Yu et al. 2012) 中的例子； 我们分析了来自 200 名患者的乳腺肿瘤组织的公开表达数据集（GSE11121，Gene Expression Omnibus）（Schmidt 等人，2008 年）。 我们从差异表达的基因中识别出 8 个基因簇，并使用 compareCluster 来比较这些基因簇的富集生物过程。
比较函数被设计为一个框架，用于比较任何种类ontology关联的基因簇，不仅是本包中提供的 groupGO、enrichGO、enrichKEGG 和enricher，还包括其他biological 和biomedical ontologies，例如，来自 DOSE（Yu et al. 2015)，从网格中富集 MeSH 和 ReactomePA 中的enrichPathway与 compareCluster 一起工作，用于比较疾病和反应组通路方面的生物学主题。 更多细节可以在 DOSE（Yu et al. 2015）和 ReactomePA 的资料中找到。
【生信分析】clusterProfiler: universal enrichment tool for functional and comparative study（3）相关推荐
- 生信分析流程构建的几大流派
导言 构建生信分析流程是生物信息学从业人员必备的技能之一,对该项能力的评估常常是各大公司招录人员的参考项目之一. 在进行 ngsjs 项目时,我做了一张示意图来表示一些高通量测序数据分析项目重现性的要 ...
- docker卸载命令_使用docker完成生信分析环境搭建
生信开发人员最头疼的问题,可能就是平台搭建和软件安装了.部署和迁移上要费很大力气.本文讲述使用docker制作一个镜像,后续通过导入自己定制的镜像,复制文件完成分析流程的部署和迁移. 如何使用dock ...
- 在B站学习大名鼎鼎的StatQuest 系列统计和生信分析视频（中文字幕）- 也见证助理教授到创业者的华丽转身...
生物信息学习的正确姿势 NGS系列文章包括NGS基础.在线绘图.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞 ...
- 生信分析R语言助力作图----单基因批量相关性分析
单基因批量相关性分析 嘻嘻嘻~~~,晚上秒变生信分析小白,一个游走在生物学和计算机变成之间的小白,享受着里面的快乐和痛苦.不停的挣扎,不停的成长,多学习,多尝试,一定会有意想不到的收获.加油!!! 首 ...
- 生信分析（1）：单变量+多变量COX分析
从TCGA上下载数据库和临床数据之后,往往需要进行COX分析,一般的分析思路是先进行单变量,在进行多变量的分析.然而,当关注的基因比较多是,手动输入就会比较麻烦.接下来介绍一种利用循环的方法,快速的对 ...
- 肿瘤/非肿瘤/单基因/单细胞/非编码：史上最全生信分析攻略！！！
解读生信之美,探讨每篇文献背后的逻辑 非肿瘤专栏:条条大路通罗马 1.4+非肿瘤生信分析+铁死亡/焦亡/自噬/代谢/免疫的万能钥匙 短评:适合一些热门机制如铁死亡/焦亡/自噬等在非肿瘤疾病中的研究 2 ...
- 面向生信分析的高性 RStudio 服务器
因需要超大内存的拼接/比对/表达量计算发愁? 为了使用组里的服务器而被困在实验室? 浪费大量的时间龟速下载 NCBI 的数据? 快来看看云筏 HPC 吧! https://my.cloudraft.c ...
- 福建农林大学朱方捷组招聘讲师/副教授/助理——生信分析方向
福建农林大学转录系统生物学课题组招聘(生信分析.组培转化) 工作地点: 福建农林大学海峡联合研究院 薪金: 18-30万 招聘岗位: 讲师/副教授/助理 实验室网址: http://hbmcsysbi ...
- 生信分析平台方案推介，助力科研
生信分析平台方案推介,助力科研 专注 专业 共赢 目前生信分析对计算性能和存储高并发性能都提出来新的要求,例如在基因测序分析中,基因序列数目庞大,对基因进行同源性搜寻.比对.分析.系统发育分析等需要对 ...
最新文章
- Uber自动驾驶汽车被赶出了亚利桑那，近300人被裁
- 宝塔如何备份网站_学习织梦网站必需会的一件事：织梦网站数据备份
- java中min用法,java11教程--类MinguoDate用法
- 威尔士柯基犬，计算机视觉，以及深度学习的力量
- Android APK方式换肤实现原理
- mysql5.7安装教程centos_CentOS7下MySQL5.7安装配置方法图文教程（YUM）
- 关于python字符编码_关于python文件的字符编码
- JBPM工作流框架应用
- SAP傻瓜式安装教程
- 诺基亚：丑小鸭的重生
- OpenCV——修改图像像素（随心所欲）
- 支付宝网页端支付接口实现案例流程
- 怎么让照片里的人嘴巴动起来_让照片动起来软件下载-让照片动起来制作软件下载-西西软件下载...
- 0089-【生物软件】-ANNOVAR基因变异注释
- Mocha.js官方文档翻译 —— 简单、灵活、有趣
- 80%中国男人不敢主动和女人搭讪
- 自学-CAD零基础视频教程网站
- java对接天猫精灵语音助手实现对公司其下的智能设备进行控制（附上源码）
- python羊车门问题_羊车门问题简析
- 这才是2019年最新资料！