医学和生信笔记公众号主要分享4大块内容：
前期主要是以医学统计和临床预测模型为主，关于生信挖掘和机器学习的内容偏少，所以后面会逐渐增加这方面的内容，除了常见的生信分析外，还会涉及一些SCI图表学习等内容。
富集分析作为了解基因功能的常用方法，也是各种生信相关SCI中的常客，基本上是必会的内容。
今天为大家带来常见的富集分析类型，其实主要就是两类：
ORA(Over-Representation Analysis)，通过超几何分布检验实现- 基因集富集分析GSEA(gene set enrichment analysis)
除此之外，还有ssGSEA、GSVA等，我们等下次再介绍。富集分析的可视化，我们也是放到下次再介绍，因为1条推文根本说不完，太多了！
目前做ORA和GSEA，肯定是首选clusterProfiler了，不管是使用的便捷性，还是对各种富集分析的完美支持，还是对各种最新数据库的支持，clusterProfiler都是最好的选择！
关于各种富集分析的底层原理，我们这里就不说了，大家感兴趣的可以自己学习下。
本文目录：
GO富集分析
KEGG富集分析
KEGG pathway ORA
KEGG pathway GSEA
KEGG module ORA
KEGG module GSEA
KEGG通路可视化
WikiPathways富集分析
Reactome富集分析
疾病富集分析
MeSH富集分析
自定义富集分析
多个基因集富集结果比较
对富集结果取子集
参考资料
library(clusterProfiler)
##
## clusterProfiler v4.6.2 For help: https:///biomedical-knowledge-mining-book/
##
## If you use clusterProfiler in published research, please cite:
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021, 2(3):100141
##
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:stats':
##
## filter
GO富集分析
clusterProfiler为我们提供了3种GO富集分析：
GO分类
基于特定GO分布对基因进行分类，使用起来超级简单，下面我们用示例数据为大家演示。
data(geneList, package="DOSE")
gene <- names(geneList)[abs(geneList) > 2]
# Entrez gene ID
head(gene)
## [1] "4312" "8318" "10874" "55143" "55388" "991"
geneList是一个有序的数值型向量，其名字是基因的ID（这里是ENTREZID）。
除了做GSEA分析需要geneList这种格式外，其他都是直接提供基因ID的向量即可：
ggo <- groupGO(gene = gene,
keyType = "ENTREZID",#任意类型都可，不一定非要ENTREZID，这里指定即可
OrgDb = "org.Hs.eg.db",
ont = "CC",# “BP”, “CC”, “MF”
level = 3,
readable = TRUE
)
##
head(ggo)
## ID Description Count GeneRatio geneID
## GO:0000133 GO:0000133 polarisome 0 0/207
## GO:0000408 GO:0000408 EKC/KEOPS complex 0 0/207
## GO:0000417 GO:0000417 HIR complex 0 0/207
## GO:0000444 GO:0000444 MIS12/MIND type complex 0 0/207
## GO:0000808 GO:0000808 origin recognition complex 0 0/207
## GO:0000930 GO:0000930 gamma-tubulin complex 0 0/207
GO ORA
这个就是大家最常见的GO富集分析了，使用enrichGO函数实现，各种参数也是和上面差不多：
ego <- enrichGO(gene = gene,
universe = names(geneList),#背景基因
OrgDb = "org.Hs.eg.db",
ont = "ALL",#"BP", "MF", "CC"
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
ego
## #
## # over-representation test
## #
## #...@organism Homo sapiens
## #...@ontology GOALL
## #...@keytype ENTREZID
## #...@gene chr [1:207] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
## #...pvalues adjusted by 'BH' with cutoff <0.01
## #...153 enriched terms found
## 'data.frame': 153 obs. of 10 variables:
## $ ONTOLOGY : chr "BP" "BP" "BP" "BP" ...
## $ ID : chr "GO:0140014" "GO:0000280" "GO:0000070" "GO:0048285" ...
## $ Description: chr "mitotic nuclear division" "nuclear division" "mitotic sister chromatid segregation" "organelle fission" ...
## $ GeneRatio : chr "32/195" "36/195" "26/195" "37/195" ...
## $ BgRatio : chr "263/11590" "357/11590" "165/11590" "391/11590" ...
## $ pvalue : num 7.10e-19 1.96e-18 2.60e-18 5.24e-18 6.54e-18 ...
## $ p.adjust : num 2.16e-15 2.65e-15 2.65e-15 3.99e-15 3.99e-15 ...
## $ qvalue : num 2.00e-15 2.44e-15 2.44e-15 3.68e-15 3.68e-15 ...
## $ geneID : chr 省略
head(ego)
## ONTOLOGY ID Description GeneRatio
## GO:0140014 BP GO:0140014 mitotic nuclear division 32/195
## GO:0000280 BP GO:0000280 nuclear division 36/195
## GO:0000070 BP GO:0000070 mitotic sister chromatid segregation 26/195
## GO:0048285 BP GO:0048285 organelle fission 37/195
## GO:0000819 BP GO:0000819 sister chromatid segregation 27/195
## GO:0007059 BP GO:0007059 chromosome segregation 30/195
## BgRatio pvalue p.adjust qvalue
## GO:0140014 263/11590 7.101540e-19 2.163839e-15 1.998149e-15
## GO:0000280 357/11590 1.956617e-18 2.645070e-15 2.442530e-15
## GO:0000070 165/11590 2.604269e-18 2.645070e-15 2.442530e-15
## GO:0048285 391/11590 5.235148e-18 3.987643e-15 3.682300e-15
## GO:0000819 188/11590 6.543556e-18 3.987643e-15 3.682300e-15
## GO:0007059 278/11590 2.850765e-16 1.447713e-13 1.336859e-13
## geneID
## GO:0140014 CDCA8/CDC20/KIF23/CENPE/MYBL2/NDC80/NCAPH/DLGAP5/UBE2C/NUSAP1/TPX2/TACC3/NEK2/UBE2S/MAD2L1/KIF18A/CDT1/BIRC5/KIF11/TTK/NCAPG/AURKB/CHEK1/TRIP13/PRC1/KIFC1/KIF18B/AURKA/CCNB1/KIF4A/PTTG1/BMP4
#省略
## Count
## GO:0140014 32
## GO:0000280 36
## GO:0000070 26
## GO:0048285 37
## GO:0000819 27
## GO:0007059 30
GO的数据库条目可能有很多是在解释同一个东西，所clusterProfielr提供了simplify能够去除冗余的GO：
simple_ego <- simplify(ego)
目前这个函数支持来自enrichGO、gseGO、compareCluster的结果，非常实用！
GO GSEA
GO版本的GSEA富集分析， gene set enrichment analysis，通过gseGO函数实现：
ego3 <- gseGO(geneList = geneList,
OrgDb = "org.Hs.eg.db",
ont = "CC",
minGSSize = 100,
maxGSSize = 500,# 只有基因集数量在[minGSSize, maxGSSize]之间的才会被检验
pvalueCutoff = 0.05,
verbose = FALSE)
head(ego3)
## ID Description setSize
## GO:0000775 GO:0000775 chromosome, centromeric region 188
## GO:0000779 GO:0000779 condensed chromosome, centromeric region 138
## GO:0000776 GO:0000776 kinetochore 130
## GO:0000228 GO:0000228 nuclear chromosome 175
## GO:0098687 GO:0098687 chromosomal region 305
## GO:0000793 GO:0000793 condensed chromosome 210
## enrichmentScore NES pvalue p.adjust qvalue rank
## GO:0000775 0.5970689 2.652931 1e-10 1.62e-09 1.073684e-09 530
## GO:0000779 0.6216009 2.648130 1e-10 1.62e-09 1.073684e-09 798
## GO:0000776 0.6230073 2.609292 1e-10 1.62e-09 1.073684e-09 449
## GO:0000228 0.5875327 2.580497 1e-10 1.62e-09 1.073684e-09 1905
## GO:0098687 0.5419949 2.545830 1e-10 1.62e-09 1.073684e-09 1721
## GO:0000793 0.5507360 2.468350 1e-10 1.62e-09 1.073684e-09 1721
## leading_edge
## GO:0000775 tags=20%, list=4%, signal=19%
## GO:0000779 tags=24%, list=6%, signal=23%
## GO:0000776 tags=21%, list=4%, signal=20%
## GO:0000228 tags=34%, list=15%, signal=29%
## GO:0098687 tags=27%, list=14%, signal=24%
## GO:0000793 tags=28%, list=14%, signal=24%
## 省略
前面说过，做GSEA富集分析对数据格式是有要求的：
KEGG富集分析
KEGG自2012年起不再免费提供FTP下载服务，但是clusterProfiler包支持使用KEGG网站下载最新的在线版本的KEGG数据，KEGG pathway和module都支持哦。
KEGG作为一个常用的数据库是有官方网站的，大家可以自己去熟悉一下：https://www./
KEGG pathway ORA
KEGG通路ORA分析，KEGG pathway over-representation analysis，通过enrichKEGG函数实现，这些函数的名字都是有规律的，方便你记忆和使用：
kk <- enrichKEGG(gene = gene,
organism = 'hsa',
pvalueCutoff = 0.05)
## Reading KEGG annotation online: "https://rest./link/hsa/pathway"...
## Reading KEGG annotation online: "https://rest./list/pathway/hsa"...
head(kk)
## ID Description
## hsa04110 hsa04110 Cell cycle
## hsa04114 hsa04114 Oocyte meiosis
## hsa04218 hsa04218 Cellular senescence
## hsa04061 hsa04061 Viral protein interaction with cytokine and cytokine receptor
## hsa03320 hsa03320 PPAR signaling pathway
## hsa04814 hsa04814 Motor proteins
## GeneRatio BgRatio pvalue p.adjust qvalue
## hsa04110 15/106 157/8465 8.177242e-10 1.717221e-07 1.695702e-07
## hsa04114 10/106 131/8465 5.049610e-06 5.302091e-04 5.235648e-04
## hsa04218 10/106 156/8465 2.366003e-05 1.639157e-03 1.618617e-03
## hsa04061 8/106 100/8465 3.326461e-05 1.639157e-03 1.618617e-03
## hsa03320 7/106 75/8465 3.902756e-05 1.639157e-03 1.618617e-03
## hsa04814 10/106 193/8465 1.433387e-04 5.016856e-03 4.953988e-03
## geneID
## hsa04110 8318/991/9133/10403/890/983/4085/81620/7272/9212/1111/9319/891/4174/9232
## hsa04114 991/9133/983/4085/51806/6790/891/9232/3708/5241
## hsa04218 2305/4605/9133/890/983/51806/1111/891/776/3708
## hsa04061 3627/10563/6373/4283/6362/6355/9547/1524
## hsa03320 4312/9415/9370/5105/2167/3158/5346
## hsa04814 9493/1062/81930/3832/3833/146909/10112/24137/4629/7802
## Count
## hsa04110 15
## hsa04114 10
## hsa04218 10
## hsa04061 8
## hsa03320 7
## hsa04814 10
KEGG pathway GSEA
KEGG通路基因集富集分析，KEGG pathway gene set enrichment analysis，通过gseKEGG实现：
kk2 <- gseKEGG(geneList = geneList,
organism = 'hsa',
minGSSize = 120,
pvalueCutoff = 0.05,
verbose = FALSE)
head(kk2)
## ID Description setSize
## hsa04110 hsa04110 Cell cycle 139
## hsa05169 hsa05169 Epstein-Barr virus infection 193
## hsa04613 hsa04613 Neutrophil extracellular trap formation 130
## hsa05166 hsa05166 Human T-cell leukemia virus 1 infection 202
## hsa04510 hsa04510 Focal adhesion 192
## hsa04218 hsa04218 Cellular senescence 141
## enrichmentScore NES pvalue p.adjust qvalue rank
## hsa04110 0.6637551 2.876742 1.000000e-10 7.900000e-09 5.263158e-09 1155
## hsa05169 0.4335010 1.956648 2.185378e-07 8.632244e-06 5.750996e-06 2820
## hsa04613 0.4496569 1.910465 7.157902e-06 1.884914e-04 1.255772e-04 2575
## hsa05166 0.3893613 1.765888 1.482811e-05 2.928551e-04 1.951067e-04 1955
## hsa04510 -0.4233234 -1.741010 3.210562e-05 5.072688e-04 3.379539e-04 2183
## hsa04218 0.4115945 1.788002 8.805195e-05 1.159351e-03 7.723855e-04 1155
## leading_edge
## hsa04110 tags=36%, list=9%, signal=33%
## hsa05169 tags=39%, list=23%, signal=31%
## hsa04613 tags=37%, list=21%, signal=30%
## hsa05166 tags=26%, list=16%, signal=22%
## hsa04510 tags=28%, list=17%, signal=23%
## hsa04218 tags=17%, list=9%, signal=16%
## core_enrichment
## hsa04110 8318/991/9133/10403/890/983/4085/81620/7272/9212/1111/9319/891/4174/9232/4171/993/990/5347/701/9700/898/23594/4998/9134/4175/4173/10926/6502/994/699/4609/5111/26271/1869/1029/8317/4176/2810/3066/1871/1031/9088/995/1019/4172/5885/11200/7027/1875
#省略
KEGG module ORA
KEGG模块ORA分析，KEGG module over-representation analysis，通过enrichMKEGG实现：
mkk <- enrichMKEGG(gene = gene,
organism = 'hsa',
pvalueCutoff = 1,
qvalueCutoff = 1)
## Reading KEGG annotation online: "https://rest./link/hsa/module"...
## Reading KEGG annotation online: "https://rest./list/module"...
head(mkk)
## ID
## M00912 M00912
## M00095 M00095
## M00053 M00053
## M00938 M00938
## M00003 M00003
## M00049 M00049
## Description
## M00912 NAD biosynthesis, tryptophan => quinolinate => NAD
## M00095 C5 isoprenoid biosynthesis, mevalonate pathway
## M00053 Deoxyribonucleotide biosynthesis, ADP/GDP/CDP/UDP => dATP/dGTP/dCTP/dUTP
## M00938 Pyrimidine deoxyribonucleotide biosynthesis, UDP => dTTP
## M00003 Gluconeogenesis, oxaloacetate => fructose-6P
## M00049 Adenine ribonucleotide biosynthesis, IMP => ADP,ATP
## GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
## M00912 2/9 12/852 0.00620266 0.03721596 0.03264558 23475/3620 2
## M00095 1/9 10/852 0.10126209 0.18449380 0.16183667 3158 1
## M00053 1/9 11/852 0.11086855 0.18449380 0.16183667 6241 1
## M00938 1/9 14/852 0.13914269 0.18449380 0.16183667 6241 1
## M00003 1/9 18/852 0.17559737 0.18449380 0.16183667 5105 1
## M00049 1/9 19/852 0.18449380 0.18449380 0.16183667 26289 1
KEGG module GSEA
KEGG模块GSEA分析，KEGG module gene set enrichment analysis
mkk2 <- gseMKEGG(geneList = geneList,
organism = 'hsa',
pvalueCutoff = 1)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
head(mkk2)
## ID Description
## M00001 M00001 Glycolysis (Embden-Meyerhof pathway), glucose => pyruvate
## M00938 M00938 Pyrimidine deoxyribonucleotide biosynthesis, UDP => dTTP
## M00002 M00002 Glycolysis, core module involving three-carbon compounds
## M00035 M00035 Methionine degradation
## M00104 M00104 Bile acid biosynthesis, cholesterol => cholate/chenodeoxycholate
## M00009 M00009 Citrate cycle (TCA cycle, Krebs cycle)
## setSize enrichmentScore NES pvalue p.adjust qvalue rank
## M00001 24 0.5739036 1.746110 0.006607415 0.2180447 0.2086552 2886
## M00938 10 0.6648004 1.550423 0.041882396 0.4607064 0.4408673 648
## M00002 11 0.6421781 1.548820 0.039734634 0.4607064 0.4408673 1381
## M00035 11 0.5824213 1.404697 0.090476190 0.6366300 0.6092153 1555
## M00104 10 -0.5876900 -1.370515 0.108247423 0.6366300 0.6092153 961
## M00009 22 0.4504911 1.340658 0.125858124 0.6366300 0.6092153 3514
## leading_edge
## M00001 tags=54%, list=23%, signal=42%
## M00938 tags=40%, list=5%, signal=38%
## M00002 tags=55%, list=11%, signal=49%
## M00035 tags=45%, list=12%, signal=40%
## M00104 tags=50%, list=8%, signal=46%
## M00009 tags=50%, list=28%, signal=36%
## core_enrichment
## M00001 5214/3101/2821/7167/2597/5230/2023/5223/5315/3099/5232/2027/5211
## M00938 6241/7298/4830/1841
## M00002 7167/2597/5230/2023/5223/5315
## M00035 875/1789/191/1788/1786
## M00104 6342/10998/1581/3295/8309
## M00009 3418/50/4190/3419/2271/3421/55753/3417/1431/6389/4191
KEGG通路可视化
执行以下代码会直接打开网页查看：
browseKEGG(kk, 'hsa04110')
或者也可以下载到本地，通过pathview这个包实现：
library("pathview")
hsa04110 <- pathview(gene.data = geneList,
pathway.id = "hsa04110",
species = "hsa",
limit = list(gene=max(abs(geneList)), cpd=1))
WikiPathways富集分析
WikiPathways(https://www./)和KEGG类似，也是一个通路数据库。
富集分析本质上是对基因注释，以了解这些基因的功能，不同的数据库有不同的侧重点，可能今天说的这些注释数据库有几个你没听过，但是不要紧，你可以自己去了解下即可，如果它们都不符合你的要求，你可以自己定义一个注释数据！
WikiPathways ORA
通过enrichWP函数实现，不过需要良好的网络环境，不然就会报错~
enrichWP(gene, organism = "Homo sapiens")
## #
## # over-representation test
## #
## #...@organism Homo sapiens
## #...@ontology WikiPathways
## #...@keytype ENTREZID
## #...@gene chr [1:207] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...8 enriched terms found
## 'data.frame': 8 obs. of 9 variables:
## $ ID : chr "WP2446" "WP2361" "WP179" "WP3942" ...
## $ Description: chr "Retinoblastoma gene in cancer" "Gastric cancer network 1" "Cell cycle" "PPAR signaling pathway" ...
## $ GeneRatio : chr "11/113" "6/113" "10/113" "7/113" ...
## $ BgRatio : chr "90/8444" "28/8444" "120/8444" "68/8444" ...
## $ pvalue : num 2.59e-08 1.49e-06 4.18e-06 3.15e-05 3.70e-05 ...
## $ p.adjust : num 7.14e-06 2.06e-04 3.84e-04 2.04e-03 2.04e-03 ...
## $ qvalue : num 6.42e-06 1.85e-04 3.46e-04 1.84e-03 1.84e-03 ...
## $ geneID : chr "8318/9133/7153/6241/890/983/81620/7272/1111/891/24137" "4605/7153/11065/22974/6286/6790" "8318/991/9133/890/983/7272/1111/891/4174/9232" "4312/9415/9370/5105/2167/3158/5346" ...
## $ Count : int 11 6 10 7 4 7 8 10
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141
WikiPathways GSEA
通过gseWP实现，也是需要良好的网络环境：
gseWP(geneList, organism = "Homo sapiens")
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
## #
## # Gene Set Enrichment Analysis
## #
## #...@organism Homo sapiens
## #...@setType WikiPathways
## #...@keytype ENTREZID
## #...@geneList Named num [1:12495] 4.57 4.51 4.42 4.14 3.88 ...
## - attr(*, "names")= chr [1:12495] "4312" "8318" "10874" "55143" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...66 enriched terms found
## 'data.frame': 66 obs. of 11 variables:
## $ ID : chr "WP2446" "WP179" "WP466" "WP2361" ...
## $ Description : chr "Retinoblastoma gene in cancer" "Cell cycle" "DNA replication" "Gastric cancer network 1" ...
## $ setSize : int 84 111 42 23 62 201 77 102 41 61 ...
## $ enrichmentScore: num 0.731 0.663 0.792 0.837 0.665 ...
## $ NES : num 2.86 2.75 2.74 2.55 2.47 ...
## $ pvalue : num 1.00e-10 1.00e-10 1.00e-10 1.41e-09 2.66e-09 ...
## $ p.adjust : num 2.12e-08 2.12e-08 2.12e-08 2.24e-07 2.83e-07 ...
## $ qvalue : num 1.77e-08 1.77e-08 1.77e-08 1.87e-07 2.35e-07 ...
## $ rank : num 1333 1234 1002 854 1111 ...
## $ leading_edge : chr "tags=51%, list=11%, signal=46%" "tags=40%, list=10%, signal=36%" "tags=55%, list=8%, signal=51%" "tags=52%, list=7%, signal=49%" ...
## $ core_enrichment: chr "8318/9133/7153/6241/890/983/81620/7272/1111/891/24137/993/898/4998/10733/9134/4175/4173/6502/5984/994/7298/3015"| __truncated__ "8318/991/9133/890/983/7272/1111/891/4174/9232/4171/993/990/5347/9700/898/23594/4998/9134/4175/4173/10926/6502/9"| __truncated__ "8318/55388/81620/4174/4171/990/23594/4998/4175/4173/10926/5984/5111/51053/8317/5427/23649/4176/5982/5557/5558/4172/5424" "4605/7153/11065/22974/6286/6790/1894/56992/4173/1063/9585/8607" ...
## #...Citation
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
## clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
## The Innovation. 2021, 2(3):100141
Reactome富集分析
Reactome(https:///)也是一个通路数据库~这个网站自带的可视化功能其实很强大，大家可以去试试。
Reactome富集分析是通过ReactomePA包实现的，注意要先加载R包。
Reactome pathway ORA
Reactome ORA分析，通过函数enrichPathway实现：
library(ReactomePA)
## ReactomePA v1.42.0 For help: https:///biomedical-knowledge-mining-book/
##
## If you use ReactomePA in published research, please cite:
## Guangchuang Yu, Qing-Yu He. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Molecular BioSystems 2016, 12(2):477-479
x <- enrichPathway(gene=gene,
pvalueCutoff = 0.05,
readable=TRUE)
head(x)
## ID Description GeneRatio
## R-HSA-69620 R-HSA-69620 Cell Cycle Checkpoints 22/143
## R-HSA-2500257 R-HSA-2500257 Resolution of Sister Chromatid Cohesion 15/143
## R-HSA-68877 R-HSA-68877 Mitotic Prometaphase 18/143
## R-HSA-69618 R-HSA-69618 Mitotic Spindle Checkpoint 14/143
## R-HSA-68882 R-HSA-68882 Mitotic Anaphase 18/143
## R-HSA-2555396 R-HSA-2555396 Mitotic Metaphase and Anaphase 18/143
## BgRatio pvalue p.adjust qvalue
## R-HSA-69620 293/10899 2.906024e-11 1.330959e-08 1.095112e-08
## R-HSA-2500257 126/10899 8.316509e-11 1.904481e-08 1.567005e-08
## R-HSA-68877 204/10899 1.539628e-10 2.350499e-08 1.933989e-08
## R-HSA-69618 113/10899 2.130715e-10 2.439669e-08 2.007358e-08
## R-HSA-68882 236/10899 1.657466e-09 1.354116e-07 1.114166e-07
## R-HSA-2555396 237/10899 1.773951e-09 1.354116e-07 1.114166e-07
## 省略
Reactome pathway GSEA
Reactome GSEA分析，通过函数gsePathway实现：
y <- gsePathway(geneList,
pvalueCutoff = 0.2,
pAdjustMethod = "BH",
verbose = FALSE)
head(y)
## ID Description setSize
## R-HSA-69620 R-HSA-69620 Cell Cycle Checkpoints 237
## R-HSA-453279 R-HSA-453279 Mitotic G1 phase and G1/S transition 142
## R-HSA-69278 R-HSA-69278 Cell Cycle, Mitotic 457
## R-HSA-69306 R-HSA-69306 DNA Replication 137
## R-HSA-69206 R-HSA-69206 G1/S Transition 124
## R-HSA-69481 R-HSA-69481 G2/M Checkpoints 134
## enrichmentScore NES pvalue p.adjust qvalue rank
## R-HSA-69620 0.6678105 3.064037 1e-10 4.105556e-09 3.064327e-09 1905
## R-HSA-453279 0.7016040 3.001137 1e-10 4.105556e-09 3.064327e-09 1155
## R-HSA-69278 0.6077642 2.974757 1e-10 4.105556e-09 3.064327e-09 1769
## R-HSA-69306 0.6968511 2.961013 1e-10 4.105556e-09 3.064327e-09 1769
## R-HSA-69206 0.7046866 2.949856 1e-10 4.105556e-09 3.064327e-09 1155
## R-HSA-69481 0.6876905 2.925269 1e-10 4.105556e-09 3.064327e-09 1905
## leading_edge
## R-HSA-69620 tags=43%, list=15%, signal=37%
## R-HSA-453279 tags=41%, list=9%, signal=37%
## R-HSA-69278 tags=36%, list=14%, signal=32%
## R-HSA-69306 tags=49%, list=14%, signal=42%
## R-HSA-69206 tags=42%, list=9%, signal=38%
## R-HSA-69481 tags=49%, list=15%, signal=42% core_enrichment
## R-HSA-69620 8318/55143/55388/991/1062/9133/10403/11065/220134/79019/55839/27338/890/983/54821/4085/81930/332/9212/1111/891/4174/4171/11004/993/990/5347/701/51512/9156/11130/79682/57405/2491/898/23594/4998/9134/4175/4173/10926/5984/1058/699/1063/85236/5688/5709/55055/641/1029/5698/1763/8970/5693/8317/4176/5713/79980/5982/9735/5721/2810/5691/9088/995/5685/7468/4172/7336/5690/5684/83990/5686/5695/11200/10213/8345/7534/80010/23198/5983/7979/54908/6396/4683/63967/3018/5699/5714/5702/3014/5905/3619/5708/55166/5692/10393/8290/5704/580/6119
## R-HSA-453279 省略
可视化
这里简单说下Reactome富集分析的可视化，因为viewPathway函数是它专用的：
viewPathway("E2F mediated regulation of DNA replication",
readable = TRUE,
foldChange = geneList)
## Loading required package: graphite
plot of chunk unnamed-chunk-18疾病富集分析
大家比较常用的是对基因做富集分析，但是疾病也是可以做富集分析的，有基因本体论（GO），也有疾病本体论（DO）
DOSE包是专门用来做疾病富集分析的。DOSE提供了五种方法来测量DO-term和基因产物之间的语义相似性，也是通过ORA和GSEA将疾病与基因相关联。
library(DOSE)
## DOSE v3.24.2 For help: https:///biomedical-knowledge-mining-book/
##
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609
Disease ORA
DOSE支持疾病本体论（DO）、癌症基因网络和疾病基因网络（DisGeNET）的富集分析。
- Over-representation analysis for disease ontology
疾病本体论（DO）的ORA分析，和GO的ORA分析一模一样：
x <- enrichDO(gene = gene,
ont = "DO",
pvalueCutoff = 0.05,
pAdjustMethod = "BH",
universe = names(geneList),
minGSSize = 10,
maxGSSize = 500,
qvalueCutoff = 0.05,
readable = FALSE)
head(x)
## ID Description GeneRatio BgRatio
## DOID:1107 DOID:1107 esophageal carcinoma 15/157 195/7744
## DOID:5041 DOID:5041 esophageal cancer 15/157 229/7744
## DOID:3459 DOID:3459 breast carcinoma 22/157 461/7744
## DOID:3748 DOID:3748 esophagus squamous cell carcinoma 10/157 131/7744
## DOID:8719 DOID:8719 in situ carcinoma 10/157 132/7744
## DOID:0060071 DOID:0060071 pre-malignant neoplasm 11/157 159/7744
## pvalue p.adjust qvalue
## DOID:1107 8.724913e-06 0.006343012 0.005675786
## DOID:5041 5.832237e-05 0.021200181 0.018970118
## DOID:3459 1.412419e-04 0.033578155 0.030046043
## DOID:3748 3.100853e-04 0.033578155 0.030046043
## DOID:8719 3.297454e-04 0.033578155 0.030046043
## DOID:0060071 3.664508e-04 0.033578155 0.030046043
## geneID
## DOID:1107 4312/6280/3868/6279/8140/890/5918/332/2146/4321/891/2066/3169/8614/11122
## DOID:5041 4312/6280/3868/6279/8140/890/5918/332/2146/4321/891/2066/3169/8614/11122
## DOID:3459 4312/6280/6279/7153/6241/4751/890/4085/332/6286/2146/6790/891/9232/771/4250/2066/3169/10647/5304/5241/10551
## DOID:3748 6280/3868/6279/8140/5918/2146/891/2066/3169/11122
## DOID:8719 6280/7153/6278/10232/332/9212/6790/891/2066/5241
## DOID:0060071 6280/7153/6278/10232/332/9212/4321/6790/891/2066/5241
## Count
## DOID:1107 15
## DOID:5041 15
## DOID:3459 22
## DOID:3748 10
## DOID:8719 10
## DOID:0060071 11
注意：ont参数可以是"DO"或"DOLite"，DOLite是一个简化后的DO版本，合并了一些冗余的DO term，但是数据不更新了，建议使用ont="DO"
- Over-representation analysis for the network of cancer gene
癌症基因网络的ORA分析：
gene2 <- names(geneList)[abs(geneList) < 3]
head(gene2)
## [1] "8685" "597" "7153" "23397" "6278" "79733"
ncg <- enrichNCG(gene2)
head(ncg)
## ID
## pan-cancer_paediatric pan-cancer_paediatric
## triple_negative_breast_cancer triple_negative_breast_cancer
## bladder_cancer bladder_cancer
## pancreatic_cancer_(all_histologies) pancreatic_cancer_(all_histologies)
## soft_tissue_sarcoma soft_tissue_sarcoma
## paediatric_high-grade_glioma paediatric_high-grade_glioma
## Description
## pan-cancer_paediatric pan-cancer_paediatric
## triple_negative_breast_cancer triple_negative_breast_cancer
## bladder_cancer bladder_cancer
## pancreatic_cancer_(all_histologies) pancreatic_cancer_(all_histologies)
## soft_tissue_sarcoma soft_tissue_sarcoma
## paediatric_high-grade_glioma paediatric_high-grade_glioma
## GeneRatio BgRatio pvalue
## pan-cancer_paediatric 162/2281 183/3177 1.833773e-08
## triple_negative_breast_cancer 71/2281 75/3177 4.290660e-07
## bladder_cancer 97/2281 112/3177 1.253690e-04
## pancreatic_cancer_(all_histologies) 40/2281 42/3177 1.262162e-04
## soft_tissue_sarcoma 26/2281 26/3177 1.742793e-04
## paediatric_high-grade_glioma 25/2281 25/3177 2.434966e-04
## p.adjust qvalue
## pan-cancer_paediatric 1.613721e-06 7.721152e-07
## triple_negative_breast_cancer 1.887890e-05 9.032967e-06
## bladder_cancer 2.776757e-03 1.328592e-03
## pancreatic_cancer_(all_histologies) 2.776757e-03 1.328592e-03
## soft_tissue_sarcoma 3.067315e-03 1.467615e-03
## paediatric_high-grade_glioma 3.073768e-03 1.470702e-03
## geneID
## pan-cancer_paediatric 2146/55353/4609/1029/3575/22806/3418/3066/2120/30012/867/7468/7545/3195/865/64109/4613/613/11177/7490/238/10736/10054/5771/4893/140885/1785/9760/3417/6597/6476/9126/4869/10320/7307/80204/1050/10992/8028/2312/6608/896/894/2196/4849/7023/5093/5079/5293/5727/55181/171017/51322/5781/3718/55294/60/673/8085/5897/4851/1665/51176/1108/7764/10664/6098/2332/2201/6495/3845/7015/1441/2782/64919/4298/23512/8239/29102/6929/8021/6134/6598/4209/5290/22941/8726/207/3717/2033/10716/4928/6932/694/5156/10019/6886/9968/7080/2623/7874/1654/4149/3020/23219/55252/55729/10735/5728/4853/23451/51341/387/3206/6146/79718/2624/63035/3815/171023/23269/25/23592/5896/7403/2260/54880/3716/9203/57178/6777/5789/4297/29072/90/546/120/25836/8289/4345/9611/5925/4763/1997/1499/7157/3399/5295/1387/4602/51564/1027/4005/2322/2078/678/6403/55709/1277/7494/64061/2625
## triple_negative_breast_cancer 省略
## Count
## pan-cancer_paediatric 162
## triple_negative_breast_cancer 71
## bladder_cancer 97
## pancreatic_cancer_(all_histologies) 40
## soft_tissue_sarcoma 26
## paediatric_high-grade_glioma 25
- Over-representation analysis for the disease gene network
疾病基因网络的ORA分析，DisGeNET(https://www./)是一个综合的、全面的基因-疾病数据库，包含基因-疾病关联和snp-基因-疾病关联。
我做网络药理学时经常使用的一个数据库：
疾病-基因关联的富集分析由enrichDGN函数实现，snp-基因-疾病关联的分析由enrichDGNv函数实现。
dgn <- enrichDGN(gene)
head(dgn)
## ID Description GeneRatio BgRatio
## C0007124 C0007124 Noninfiltrating Intraductal Carcinoma 24/204 486/21671
## C3642347 C3642347 Basal-Like Breast Carcinoma 16/204 245/21671
## C1609538 C1609538 Latent Tuberculosis 14/204 183/21671
## C1827293 C1827293 Carcinoma of urinary bladder, invasive 16/204 287/21671
## C0278601 C0278601 Inflammatory Breast Carcinoma 13/204 183/21671
## C4733092 C4733092 estrogen receptor-negative breast cancer 17/204 356/21671
## pvalue p.adjust qvalue
## C0007124 3.806436e-11 1.568251e-07 1.271350e-07
## C3642347 1.555102e-09 3.038993e-06 2.463650e-06
## C1609538 2.212859e-09 3.038993e-06 2.463650e-06
## C1827293 1.491211e-08 1.535947e-05 1.245161e-05
## C0278601 2.050738e-08 1.689808e-05 1.369893e-05
## C4733092 5.070020e-08 3.481414e-05 2.822311e-05
## geneID
## C0007124 4312/6280/6279/7153/6278/9582/55872/3669/4751/5080/2146/1111/3887/6790/9370/125/8839/2066/3169/9547/10647/2625/7021/5241
## C3642347 2305/1062/4605/9833/7368/11065/10232/55765/5163/2146/2568/3620/6790/3169/2625/5241
## C1609538 8685/3627/3669/6373/820/3002/50852/3902/3620/6355/10578/9370/3117/23158
## C1827293 2305/6279/55165/220134/2146/9212/4321/6790/891/8544/9122/23158/652/2625/5241/10551
## C0278601 4312/6278/11065/9582/51203/5918/2146/4321/3887/2066/9547/10647/5241
## C4733092 2305/6278/79733/6241/81930/81620/2146/3620/80129/4137/8839/3169/1408/5304/2625/5241/10551
## Count
## C0007124 24
## C3642347 16
## C1609538 14
## C1827293 16
## C0278601 13
## C4733092 17
snp <- c("rs1401296", "rs9315050", "rs5498", "rs1524668", "rs147377392",
"rs841", "rs909253", "rs7193343", "rs3918232", "rs3760396",
"rs2231137", "rs10947803", "rs17222919", "rs386602276", "rs11053646",
"rs1805192", "rs139564723", "rs2230806", "rs20417", "rs966221")
dgnv <- enrichDGNv(snp)
head(dgnv)
## ID Description GeneRatio BgRatio pvalue
## C0010054 C0010054 Coronary Arteriosclerosis 6/17 440/194515 1.568917e-12
## C0151744 C0151744 Myocardial Ischemia 4/17 103/194515 1.754840e-10
## C0031099 C0031099 Periodontitis 4/17 116/194515 2.839985e-10
## C0007785 C0007785 Cerebral Infarction 4/17 123/194515 3.599531e-10
## C0003850 C0003850 Arteriosclerosis 4/17 267/194515 8.145389e-09
## C0004153 C0004153 Atherosclerosis 4/17 281/194515 9.996713e-09
## p.adjust qvalue
## C0010054 2.761295e-10 NA
## C0151744 1.544259e-08 NA
## C0031099 1.583793e-08 NA
## C0007785 1.583793e-08 NA
## C0003850 2.867177e-07 NA
## C0004153 2.932369e-07 NA
## geneID Count
## C0010054 rs5498/rs147377392/rs11053646/rs1805192/rs2230806/rs20417 6
## C0151744 rs5498/rs147377392/rs11053646/rs1805192 4
## C0031099 rs5498/rs909253/rs1805192/rs20417 4
## C0007785 rs147377392/rs11053646/rs1805192/rs2230806 4
## C0003850 rs5498/rs1805192/rs2230806/rs20417 4
## C0004153 rs5498/rs1805192/rs2230806/rs20417 4
Disease GSEA
DO GSEA分析：
y <- gseDO(geneList,
minGSSize = 120,
pvalueCutoff = 0.2,
pAdjustMethod = "BH",
verbose = FALSE)
head(y, 3)
## ID Description setSize
## DOID:399 DOID:399 tuberculosis 181
## DOID:612 DOID:612 primary immunodeficiency disease 305
## DOID:0050338 DOID:0050338 primary bacterial infectious disease 264
## enrichmentScore NES pvalue p.adjust qvalue
## DOID:399 0.5030583 2.252875 1.00000e-10 2.180000e-08 9.789474e-09
## DOID:612 0.4023288 1.898732 4.87282e-09 5.311374e-07 2.385117e-07
## DOID:0050338 0.4150405 1.945371 1.17702e-08 6.414760e-07 2.880602e-07
## rank leading_edge
## DOID:399 2082 tags=40%, list=17%, signal=34%
## DOID:612 2099 tags=33%, list=17%, signal=29%
## DOID:0050338 1864 tags=33%, list=15%, signal=29%
## 省略
癌症基因网络的GSEA分析：
ncg <- gseNCG(geneList,
pvalueCutoff = 0.5,
pAdjustMethod = "BH",
verbose = FALSE)
ncg <- setReadable(ncg, 'org.Hs.eg.db')
head(ncg, 3)
## ID
## pan-gynecological and breast pan-gynecological and breast
## pan-gastric pan-gastric
## esophageal_adenocarcinoma esophageal_adenocarcinoma
## Description setSize
## pan-gynecological and breast pan-gynecological and breast 43
## pan-gastric pan-gastric 49
## esophageal_adenocarcinoma esophageal_adenocarcinoma 92
## enrichmentScore NES pvalue p.adjust
## pan-gynecological and breast -0.5263429 -1.734359 0.004764925 0.1718695
## pan-gastric -0.4993803 -1.683334 0.006365539 0.1718695
## esophageal_adenocarcinoma -0.4235434 -1.550378 0.004980207 0.1718695
## qvalue rank leading_edge
## pan-gynecological and breast 0.1585801 2464 tags=44%, list=20%, signal=36%
## pan-gastric 0.1585801 3280 tags=49%, list=26%, signal=36%
## esophageal_adenocarcinoma 0.1585801 2464 tags=34%, list=20%, signal=27%
## core_enrichment
## pan-gynecological and breast ATM/ZC3H13/NIPBL/SPOP/ARID1A/RASA1/RB1/RNF43/MAP2K4/NF1/CTNNB1/TP53/PIK3R1/CDKN1B/CCND1/ARID5B/MAP3K1/TBX3/GATA3
## pan-gastric BCOR/SOX9/TCF7L2/ATM/CALD1/SEMG2/HTR7/ARID1A/RASA1/RB1/TTBK2/RNF43/CTNNB1/TP53/BCL9/SMAD3/APC/ZFP36L2/TGFBR2/MUC6/MAP3K1/CACNA1C/ATP8B1/CYP4B1
## esophageal_adenocarcinoma NIPBL/ARID1A/RASA1/MSH3/RNF43/CTNNB1/DNAH7/TP53/GATAD1/PIK3R1/STK11/CDKN1B/SPART/APC/SLIT3/LIN7A/SLIT2/TGFBR2/EPHA3/MUC6/ABCB1/NAV3/SYNE1/MAP3K1/SEMA5A/FTO/NUAK1/PCDH9/CHL1/ABCG2/TRIM58
疾病基因网络的GSEA分析：
dgn <- gseDGN(geneList,
pvalueCutoff = 0.2,
pAdjustMethod = "BH",
verbose = FALSE)
dgn <- setReadable(dgn, 'org.Hs.eg.db')
head(dgn, 3)
## ID Description setSize enrichmentScore
## C0024266 C0024266 Lymphocytic Choriomeningitis 120 0.5712593
## C4721414 C4721414 Mantle cell lymphoma 368 0.4107437
## C0205682 C0205682 Waist-Hip Ratio 401 -0.4425633
## NES pvalue p.adjust qvalue rank
## C0024266 2.437904 1e-10 2.05275e-07 1.757895e-07 2579
## C4721414 1.951356 1e-10 2.05275e-07 1.757895e-07 1745
## C0205682 -1.946241 1e-10 2.05275e-07 1.757895e-07 2011
## leading_edge
## C0024266 tags=48%, list=21%, signal=38%
## C4721414 tags=26%, list=14%, signal=23%
## C0205682 tags=28%, list=16%, signal=24%
## 省略
MeSH富集分析
meshes包支持使用MeSH注释（医学主题词表）对基因列表或整个表达谱进行富集分析（ORA和GSEA）。gendoo/gene2pubmed/RBBH的数据源都支持。
但是进行MeSH富集分析需要通过AnnotationHub下载数据，这一步由于网络问题经常失败...
AnnotationHub在bioconductor 3.16版本（对应R4.2.0版本）之后得到了重用，专门用于存放各种注释数据的R包，通过它可以非常方便的进行查询、下载所有的注释数据！在此之前的版本都是各种注释数据都有自己的专属R包，不管R包多大，我们都可以通过更换镜像实现快速下载！
现在统一之后确实方便大家查找，但是我觉得这样做让大陆用户很为难......基本上都会卡在下载数据这一步（因为这个东西没有镜像可以用了...），即使通过其他方式下载了数据，但是不同的注释数据还有不同的处理方法。
MeSH ORA
首先自己构建下数据库，AnnotationHub下载数据太慢了...（但是如果网络环境好也是嗖嗖嗖地下载）
# 需要bioc版本在3.14以上
library(AnnotationHub)
library(MeSHDbi)
ah <- AnnotationHub() #第一次运行会问你要不要建立本地缓存，选yes
hsa <- query(ah, c("MeSHDb", "Homo sapiens"))
file_hsa <- hsa[[1]] # 这一步是下载真正的注释，还好不用每次都下载
#file_hsa <- "C:/Users/liyue/AppData/Local/R/cache/R/AnnotationHub/mesh_hs_db_98354"#这个是我自己通过迅雷下载的，可以读取，但是不知道为啥富集分析不能用
db <- MeSHDbi::MeSHDb(file_hsa)
来自官方的提示：
MeSH-related packages (MeSH.XXX.eg.db, MeSH.db, MeSH.AOR.db, and MeSH.PCR.db) are deprecated since Bioconductor 3.14. Use AnnotationHub instead. For details, check the vignette of MeSHDbi
以前是可以安装各种注释R包的，换个镜像即可嗖嗖的下载，好日子一去不复返了！
## BioC 2.14 to BioC 3.13 ##
library(MeSH.Hsa.eg.db)
db <- MeSH.Hsa.eg.db
接下来就是进行MeSH的ORA分析：
library(meshes)
## meshes v1.24.0
##
## If you use meshes in published research, please cite the most appropriate paper(s):
##
## Guangchuang Yu. Using meshes for MeSH term enrichment and semantic analyses. Bioinformatics 2018, 34(21):3766-3767, doi:10.1093/bioinformatics/bty410
##
## Attaching package: 'meshes'
## The following object is masked from 'package:DOSE':
##
## geneSim
data(geneList, package="DOSE")
de <- names(geneList)[1:100]
x <- enrichMeSH(de, MeSHDb = db, database='gendoo', category = 'C')
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'select': object 'db' not found
head(x)
## ID Description GeneRatio BgRatio
## DOID:1107 DOID:1107 esophageal carcinoma 15/157 195/7744
## DOID:5041 DOID:5041 esophageal cancer 15/157 229/7744
## DOID:3459 DOID:3459 breast carcinoma 22/157 461/7744
## DOID:3748 DOID:3748 esophagus squamous cell carcinoma 10/157 131/7744
## DOID:8719 DOID:8719 in situ carcinoma 10/157 132/7744
## DOID:0060071 DOID:0060071 pre-malignant neoplasm 11/157 159/7744
## pvalue p.adjust qvalue
## DOID:1107 8.724913e-06 0.006343012 0.005675786
## DOID:5041 5.832237e-05 0.021200181 0.018970118
## DOID:3459 1.412419e-04 0.033578155 0.030046043
## DOID:3748 3.100853e-04 0.033578155 0.030046043
## DOID:8719 3.297454e-04 0.033578155 0.030046043
## DOID:0060071 3.664508e-04 0.033578155 0.030046043
## geneID
## DOID:1107 4312/6280/3868/6279/8140/890/5918/332/2146/4321/891/2066/3169/8614/11122
## DOID:5041 4312/6280/3868/6279/8140/890/5918/332/2146/4321/891/2066/3169/8614/11122
## DOID:3459 4312/6280/6279/7153/6241/4751/890/4085/332/6286/2146/6790/891/9232/771/4250/2066/3169/10647/5304/5241/10551
## DOID:3748 6280/3868/6279/8140/5918/2146/891/2066/3169/11122
## DOID:8719 6280/7153/6278/10232/332/9212/6790/891/2066/5241
## DOID:0060071 6280/7153/6278/10232/332/9212/4321/6790/891/2066/5241
## Count
## DOID:1107 15
## DOID:5041 15
## DOID:3459 22
## DOID:3748 10
## DOID:8719 10
## DOID:0060071 11
MeSH GSEA
然后是GSEA分析：
y <- gseMeSH(geneList, MeSHDb = db, database = 'gene2pubmed', category = "G")
## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'select': object 'db' not found
head(y)
## ID Description
## DOID:399 DOID:399 tuberculosis
## DOID:612 DOID:612 primary immunodeficiency disease
## DOID:0050338 DOID:0050338 primary bacterial infectious disease
## DOID:104 DOID:104 bacterial infectious disease
## DOID:526 DOID:526 human immunodeficiency virus infectious disease
## DOID:2237 DOID:2237 hepatitis
## setSize enrichmentScore NES pvalue p.adjust
## DOID:399 181 0.5030583 2.252875 1.000000e-10 2.180000e-08
## DOID:612 305 0.4023288 1.898732 4.872820e-09 5.311374e-07
## DOID:0050338 264 0.4150405 1.945371 1.177020e-08 6.414760e-07
## DOID:104 321 0.3899673 1.856122 9.801772e-09 6.414760e-07
## DOID:526 201 0.4147782 1.874250 4.756033e-07 2.073630e-05
## DOID:2237 406 0.3328749 1.613797 1.567985e-06 5.697013e-05
## qvalue rank leading_edge
## DOID:399 9.789474e-09 2082 tags=40%, list=17%, signal=34%
## DOID:612 2.385117e-07 2099 tags=33%, list=17%, signal=29%
## DOID:0050338 2.880602e-07 1864 tags=33%, list=15%, signal=29%
## DOID:104 2.880602e-07 2330 tags=36%, list=19%, signal=30%
## DOID:526 9.311812e-06 2099 tags=36%, list=17%, signal=31%
## DOID:2237 2.558292e-05 2082 tags=27%, list=17%, signal=23%
## 省略
自定义富集分析
有些用户需要自定义自己的注释数据库或者没有支持的物种，也可以自己DIY。
clusterProfiler提供了用于超几何检验的enricher()函数和用于基因集富集分析的GSEA()函数，可以让用户根据自己的注释信息进行富集分析。
需要两个额外的参数TERM2GENE和TERM2NAME。TERM2GENE是一个数据框，第一列是term-ID，第二列是相应的映射基因；TERM2NAME也是一个数据框，第一列是term-ID，第二列是相应的术语名称。其中TERM2NAME是可选的。
MiSigdb
下面是一个用MiSigdb数据库进行富集分析的例子。MiSigdb可以非常方便的进行查询和下载MiSigdb数据库的基因集。
查看支持的物种：
library(msigdbr)
msigdbr_show_species()
## [1] "Anolis carolinensis" "Bos taurus"
## [3] "Caenorhabditis elegans" "Canis lupus familiaris"
## [5] "Danio rerio" "Drosophila melanogaster"
## [7] "Equus caballus" "Felis catus"
## [9] "Gallus gallus" "Homo sapiens"
## [11] "Macaca mulatta" "Monodelphis domestica"
## [13] "Mus musculus" "Ornithorhynchus anatinus"
## [15] "Pan troglodytes" "Rattus norvegicus"
## [17] "Saccharomyces cerevisiae" "Schizosaccharomyces pombe 972h-"
## [19] "Sus scrofa" "Xenopus tropicalis"
选择我们想要的数据下载：
m_df <- msigdbr(species = "Homo sapiens")#下载数据
head(m_df, 2) %>% as.data.frame
## gs_cat gs_subcat gs_name gene_symbol entrez_gene ensembl_gene
## 1 C3 MIR:MIR_Legacy AAACCAC_MIR140 ABCC4 10257 ENSG00000125257
## 2 C3 MIR:MIR_Legacy AAACCAC_MIR140 ABRAXAS2 23172 ENSG00000165660
## human_gene_symbol human_entrez_gene human_ensembl_gene gs_id gs_pmid
## 1 ABCC4 10257 ENSG00000125257 M12609
## 2 ABRAXAS2 23172 ENSG00000165660 M12609
## gs_geoid gs_exact_source gs_url
## 1
## 2
## gs_description
## 1 Genes having at least one occurence of the motif AAACCAC in their 3' untranslated region. The motif represents putative target (that is, seed match) of human mature miRNA hsa-miR-140 (v7.1 miRBase).
## 2 Genes having at least one occurence of the motif AAACCAC in their 3' untranslated region. The motif represents putative target (that is, seed match) of human mature miRNA hsa-miR-140 (v7.1 miRBase).
我们使用人类的C6注释集：
m_t2g <- msigdbr(species = "Homo sapiens", category = "C6") %>%
dplyr::select(gs_name, entrez_gene)
head(m_t2g)
## # A tibble: 6 × 2
## gs_name entrez_gene
## <chr> <int>
## 1 AKT_UP.V1_DN 57007
## 2 AKT_UP.V1_DN 22859
## 3 AKT_UP.V1_DN 22859
## 4 AKT_UP.V1_DN 137872
## 5 AKT_UP.V1_DN 249
## 6 AKT_UP.V1_DN 271
有了这个数据之后，我们就可以根据这个自定义的注释信息做富集分析了：
首先是ORA分析：
em <- enricher(gene, TERM2GENE=m_t2g)
head(em)
## ID Description GeneRatio
## RPS14_DN.V1_DN RPS14_DN.V1_DN RPS14_DN.V1_DN 22/183
## GCNP_SHH_UP_LATE.V1_UP GCNP_SHH_UP_LATE.V1_UP GCNP_SHH_UP_LATE.V1_UP 16/183
## PRC2_EZH2_UP.V1_DN PRC2_EZH2_UP.V1_DN PRC2_EZH2_UP.V1_DN 15/183
## VEGF_A_UP.V1_DN VEGF_A_UP.V1_DN VEGF_A_UP.V1_DN 15/183
## RB_P107_DN.V1_UP RB_P107_DN.V1_UP RB_P107_DN.V1_UP 10/183
## E2F1_UP.V1_UP E2F1_UP.V1_UP E2F1_UP.V1_UP 12/183
## BgRatio pvalue p.adjust qvalue
## RPS14_DN.V1_DN 186/10922 4.657342e-13 7.870908e-11 6.422230e-11
## GCNP_SHH_UP_LATE.V1_UP 181/10922 5.765032e-08 4.871452e-06 3.974838e-06
## PRC2_EZH2_UP.V1_DN 192/10922 7.574546e-07 3.419660e-05 2.790255e-05
## VEGF_A_UP.V1_DN 193/10922 8.093869e-07 3.419660e-05 2.790255e-05
## RB_P107_DN.V1_UP 130/10922 6.394609e-05 2.013355e-03 1.642787e-03
## E2F1_UP.V1_UP 188/10922 7.477287e-05 2.013355e-03 1.642787e-03
## geneID
## RPS14_DN.V1_DN 10874/55388/991/9493/1062/4605/9133/23397/79733/9787/55872/83461/54821/51659/9319/9055/10112/4174/5105/2532/7021/79901
## GCNP_SHH_UP_LATE.V1_UP 55388/7153/79733/6241/9787/51203/983/9212/1111/9319/9055/3833/6790/4174/3169/1580
## PRC2_EZH2_UP.V1_DN 8318/55388/4605/23397/9787/55355/10460/6362/81620/2146/7272/9212/11182/3887/24137
## VEGF_A_UP.V1_DN 8318/9493/1062/9133/10403/6241/9787/4085/332/3832/7272/891/23362/2167/10234
## RB_P107_DN.V1_UP 8318/23397/79733/6241/4085/8208/9055/24137/4174/1307
## E2F1_UP.V1_UP 55388/7153/23397/79733/9787/2146/2842/9212/8208/1111/9055/3833
## Count
## RPS14_DN.V1_DN 22
## GCNP_SHH_UP_LATE.V1_UP 16
## PRC2_EZH2_UP.V1_DN 15
## VEGF_A_UP.V1_DN 15
## RB_P107_DN.V1_UP 10
## E2F1_UP.V1_UP 12
再试一下C3基因集：
C3_t2g <- msigdbr(species = "Homo sapiens", category = "C3") %>%
dplyr::select(gs_name, entrez_gene)
head(C3_t2g)
## # A tibble: 6 × 2
## gs_name entrez_gene
## <chr> <int>
## 1 AAACCAC_MIR140 10257
## 2 AAACCAC_MIR140 23172
## 3 AAACCAC_MIR140 81
## 4 AAACCAC_MIR140 81
## 5 AAACCAC_MIR140 90
## 6 AAACCAC_MIR140 8754
然后是GSEA分析：
em2 <- GSEA(geneList, TERM2GENE = C3_t2g)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
head(em2)
## ID Description setSize
## HSD17B8_TARGET_GENES HSD17B8_TARGET_GENES HSD17B8_TARGET_GENES 401
## E2F1_Q3 E2F1_Q3 E2F1_Q3 187
## E2F1_Q6 E2F1_Q6 E2F1_Q6 185
## SGCGSSAAA_E2F1DP2_01 SGCGSSAAA_E2F1DP2_01 SGCGSSAAA_E2F1DP2_01 129
## E2F_Q4 E2F_Q4 E2F_Q4 193
## E2F1DP1_01 E2F1DP1_01 E2F1DP1_01 184
## enrichmentScore NES pvalue p.adjust
## HSD17B8_TARGET_GENES 0.6368099 3.107103 1.000000e-10 1.655000e-07
## E2F1_Q3 0.4983913 2.220028 1.000000e-10 1.655000e-07
## E2F1_Q6 0.4906657 2.184070 2.011165e-10 2.218985e-07
## SGCGSSAAA_E2F1DP2_01 0.5408182 2.289484 1.032436e-09 8.543410e-07
## E2F_Q4 0.4685394 2.099303 2.297950e-09 9.507767e-07
## E2F1DP1_01 0.4675989 2.083101 2.123251e-09 9.507767e-07
## qvalue rank leading_edge
## HSD17B8_TARGET_GENES 1.310526e-07 1042 tags=36%, list=8%, signal=34%
## E2F1_Q3 1.310526e-07 2045 tags=36%, list=16%, signal=30%
## E2F1_Q6 1.757123e-07 1819 tags=34%, list=15%, signal=30%
## SGCGSSAAA_E2F1DP2_01 6.765174e-07 1797 tags=38%, list=14%, signal=33%
## E2F_Q4 7.528809e-07 1663 tags=29%, list=13%, signal=26%
## E2F1DP1_01 7.528809e-07 1797 tags=33%, list=14%, signal=28%
#省略
cellmarker
除了msigdb数据库，还有一个cellmarker数据库，做单细胞的一定见过，我们先下载它的注释数据：
cell_marker_data <- vroom::vroom('http://xteam./CellMarker/download/Human_cell_markers.txt')
## Rows: 2868 Columns: 15
## ── Column specification ─────────────────────────────
## Delimiter: "\t"
## chr (15): speciesType, tissueType, UberonOntologyID, cancerType, cellType, c...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## instead of `cellName`, users can use other features (e.g. `cancerType`)
cells <- cell_marker_data %>%
dplyr::select(cellName, geneID) %>%
dplyr::mutate(geneID = strsplit(geneID, ', ')) %>%
tidyr::unnest()
有了数据就可以使用通用的enricher函数进行ORA富集分析了：
x <- enricher(gene, TERM2GENE = cells)
head(x)
## ID Description GeneRatio
## Neural progenitor cell Neural progenitor cell Neural progenitor cell 33/153
## MKI67+ progenitor cell MKI67+ progenitor cell MKI67+ progenitor cell 20/153
## DCLK1+ progenitor cell DCLK1+ progenitor cell DCLK1+ progenitor cell 7/153
## BgRatio pvalue p.adjust qvalue
## Neural progenitor cell 176/11582 3.395466e-29 3.259648e-27 2.716373e-27
## MKI67+ progenitor cell 111/11582 1.426062e-17 6.845100e-16 5.704250e-16
## DCLK1+ progenitor cell 99/11582 3.298755e-04 1.055601e-02 8.796679e-03
## geneID
## Neural progenitor cell 991/9493/9833/9133/10403/7153/259266/6241/9787/11065/55355/55872/51203/83461/22974/10460/79019/55839/890/983/4085/5080/332/3832/9212/51659/9055/891/4174/9232/1602/2018/4137
## MKI67+ progenitor cell 991/9133/10403/7153/259266/6241/9787/11065/51203/22974/10460/890/983/332/9212/9055/891/9232/2167/51705
## DCLK1+ progenitor cell 1307/730/54829/10234/79901/57758/4969
## Count
## Neural progenitor cell 33
## MKI67+ progenitor cell 20
## DCLK1+ progenitor cell 7
使用GSEA进行基因集富集分析：
y <- GSEA(geneList, TERM2GENE = cells)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...
head(y)
## ID
## Neural progenitor cell Neural progenitor cell
## FOXN4+ cell FOXN4+ cell
## Paneth cell Paneth cell
## Leydig precursor cell Leydig precursor cell
## FGFR1HighNME5- epithelial cell FGFR1HighNME5- epithelial cell
## Ciliated epithelial cell Ciliated epithelial cell
## Description setSize
## Neural progenitor cell Neural progenitor cell 149
## FOXN4+ cell FOXN4+ cell 63
## Paneth cell Paneth cell 261
## Leydig precursor cell Leydig precursor cell 184
## FGFR1HighNME5- epithelial cell FGFR1HighNME5- epithelial cell 131
## Ciliated epithelial cell Ciliated epithelial cell 217
## enrichmentScore NES pvalue
## Neural progenitor cell 0.7160571 3.085727 1.000000e-10
## FOXN4+ cell 0.7136090 2.685258 1.000000e-10
## Paneth cell 0.5381280 2.489150 1.000000e-10
## Leydig precursor cell -0.5710583 -2.308894 1.000000e-10
## FGFR1HighNME5- epithelial cell -0.5739472 -2.229562 1.476993e-10
## Ciliated epithelial cell -0.5012031 -2.064339 3.336445e-10
## p.adjust qvalue rank
## Neural progenitor cell 4.125000e-09 2.315789e-09 684
## FOXN4+ cell 4.125000e-09 2.315789e-09 1419
## Paneth cell 4.125000e-09 2.315789e-09 2187
## Leydig precursor cell 4.125000e-09 2.315789e-09 2609
## FGFR1HighNME5- epithelial cell 4.874076e-09 2.736323e-09 1404
## Ciliated epithelial cell 9.175223e-09 5.151002e-09 1440
## leading_edge
## Neural progenitor cell tags=47%, list=5%, signal=45%
## FOXN4+ cell tags=51%, list=11%, signal=45%
## Paneth cell tags=45%, list=18%, signal=38%
## Leydig precursor cell tags=45%, list=21%, signal=36%
## FGFR1HighNME5- epithelial cell tags=39%, list=11%, signal=35%
## Ciliated epithelial cell tags=29%, list=12%, signal=27%
#省略
DO, GO, Reactome的分析是有readable参数的，其实是根据OrgDb进行转换的，这和setReadable函数的原理一样，所以用到OrgDb的都有这个参数。
多个基因集富集结果比较
同时对多个基因集进行富集分析，这个功能在单细胞里面简直是绝配！
data(gcSample)
str(gcSample)
## List of 8
## $ X1: chr [1:216] "4597" "7111" "5266" "2175" ...
## $ X2: chr [1:805] "23450" "5160" "7126" "26118" ...
## $ X3: chr [1:392] "894" "7057" "22906" "3339" ...
## $ X4: chr [1:838] "5573" "7453" "5245" "23450" ...
## $ X5: chr [1:929] "5982" "7318" "6352" "2101" ...
## $ X6: chr [1:585] "5337" "9295" "4035" "811" ...
## $ X7: chr [1:582] "2621" "2665" "5690" "3608" ...
## $ X8: chr [1:237] "2665" "4735" "1327" "3192" ...
同时对这8个基因集分别进行富集分析:
ck <- compareCluster(geneCluster = gcSample, fun = enrichKEGG)
ck <- setReadable(ck, OrgDb = org.Hs.eg.db, keyType="ENTREZID")
head(ck)
## Cluster ID Description GeneRatio BgRatio
## 1 X2 hsa05169 Epstein-Barr virus infection 23/401 202/8465
## 2 X2 hsa05340 Primary immunodeficiency 8/401 38/8465
## 3 X2 hsa04110 Cell cycle 18/401 157/8465
## 4 X3 hsa04512 ECM-receptor interaction 9/190 89/8465
## 5 X3 hsa04060 Cytokine-cytokine receptor interaction 17/190 297/8465
## 6 X3 hsa04151 PI3K-Akt signaling pathway 19/190 354/8465
## pvalue p.adjust qvalue
## 1 7.981195e-05 0.02466189 0.02310346
## 2 3.300857e-04 0.04463551 0.04181490
## 3 4.333545e-04 0.04463551 0.04181490
## 4 1.617209e-04 0.03272190 0.03037690
## 5 3.565676e-04 0.03272190 0.03037690
## 6 3.622350e-04 0.03272190 0.03037690
## geneID
## 1 LYN/ICAM1/TNFAIP3/E2F1/CCNA2/E2F3/BAK1/RUNX3/BID/IKBKE/TAP2/FAS/CCNE2/RELB/CD3E/ENTPD3/SYK/TRAF3/IKBKB/CD247/NEDD4/CD40/STAT1
## 2 ADA/TAP2/LCK/CD79A/CD3E/CD8A/CD40/DCLRE1C
## 3 CDC20/E2F1/CCNA2/E2F3/BUB1B/CDC6/DBF4/PKMYT1/CDC7/ESPL1/CCNE2/CDKN2A/SFN/BUB1/CHEK2/ORC6/CDC14B/MCM4
## 4 THBS1/HSPG2/COL9A3/ITGB7/CHAD/ITGA7/LAMA4/ITGB8/ITGB5
## 5 CXCL1/TNFRSF11B/LIFR/XCL1/GDF5/TNFRSF17/TNFRSF11A/IL5RA/EPOR/CSF2RA/TNFRSF25/BMP7/BMP4/GDF11/IL21R/IL23A/TGFB2
## 6 CCND2/THBS1/STK11/FGF2/COL9A3/ITGB7/FGF7/ERBB4/CHAD/FGF18/CDK6/PCK1/ITGA7/HGF/EPOR/LAMA4/IKBKB/ITGB8/ITGB5
## Count
## 1 23
## 2 8
## 3 18
## 4 9
## 5 17
## 6 19
支持使用R语言中的formula形式：
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"
str(mydf)
## 'data.frame': 1324 obs. of 4 variables:
## $ Entrez : chr "4312" "8318" "10874" "55143" ...
## $ FC : num 4.57 4.51 4.42 4.14 3.88 ...
## $ group : chr "upregulated" "upregulated" "upregulated" "upregulated" ...
## $ othergroup: chr "B" "B" "B" "B" ...
formula_res <- compareCluster(Entrez~group+othergroup, data=mydf, fun="enrichKEGG")
head(formula_res)
## Cluster group othergroup ID
## 1 downregulated.A downregulated A hsa04974
## 2 downregulated.A downregulated A hsa04510
## 3 downregulated.A downregulated A hsa04151
## 4 downregulated.A downregulated A hsa04512
## 5 downregulated.B downregulated B hsa03320
## 6 upregulated.A upregulated A hsa04110
## Description GeneRatio BgRatio pvalue p.adjust
## 1 Protein digestion and absorption 16/326 103/8465 1.731006e-06 4.864127e-04
## 2 Focal adhesion 20/326 203/8465 1.030192e-04 1.447420e-02
## 3 PI3K-Akt signaling pathway 28/326 354/8465 2.240338e-04 2.098450e-02
## 4 ECM-receptor interaction 11/326 89/8465 5.802099e-04 4.075974e-02
## 5 PPAR signaling pathway 5/43 75/8465 3.532307e-05 5.581045e-03
## 6 Cell cycle 21/224 157/8465 7.699196e-10 1.963295e-07
## qvalue
## 1 4.555279e-04
## 2 1.355516e-02
## 3 1.965209e-02
## 4 3.817170e-02
## 5 5.131141e-03
## 6 1.750554e-07
## 省略
可视化也是一样的支持：
dotplot(ck)
plot of chunk unnamed-chunk-43dotplot(formula_res)
plot of chunk unnamed-chunk-43cnetplot(ck)
plot of chunk unnamed-chunk-44除了以上的富集分析外，还有几个其他类似方法没介绍到，单基因富集分析、基因集变异分析GSVA、ssGSEA，这些等以后再介绍。
由于输出内容非常多，本文中很多输出都省略了，大家可以自己敲一敲代码看看详细的结果。
对富集结果取子集
富集分析得到的结果都是一个对象，但是可以直接用as.data.frame转化为data.frame，然后你就可以直接保存为CSV格式了。
#以ego为例，其他也是一样
class(ego)
## [1] "enrichResult"
## attr(,"package")
## [1] "DOSE"
go_res <- as.data.frame(ego)
class(go_res)
## [1] "data.frame"
write.csv(go_res,file = "gores.csv",quote = F)
所有的富集结果都是可以像这样转化的，方便你进行个性化的操作。
当然，如果你只是单纯的想要对这个结果取子集，也是很简单的，因为：
基本上data.frame取子集那些方法，包括[, [[, >$等操作符，都被我重新定义，也就是说，对于富集分析的结果，你完全可以像data.frame一样的操作。参考：https://mp.weixin.qq.com/s/k-Lbpz8M616Wo0St0TQ-2Q
不过在此之前，我们先探索下富集分析的结果。
dim(ego)
## [1] 153 10
#head(ego)
这个富集结果共有153行，10列。其实这个数据就在ego@result这个数据框中，所有的操作也是在对这个数据框进行操作。
如果你不知道这10个列都是哪些信息（即你不知道富集分析可以得到哪些信息），可以直接查看这个结果的列名：
names(ego@result)
## [1] "ONTOLOGY" "ID" "Description" "GeneRatio" "BgRatio" "pvalue"
## [7] "p.adjust" "qvalue" "geneID" "Count"
可以看到，就是常见的富集到的条目，P值等信息，其实各种数据操作和画图等，都是针对这个数据的，只不过clusterprofiler和enrichplot对这个过程进行了简化，不用你提取数据也能直接使用。
取子集非常简单：
#asis=T让子集也是富集分析结果对象，可以直接用
y = ego[ego$pvalue < 0.001, asis=T]
dim(y)
## [1] 153 10
class(y)
## [1] "enrichResult"
## attr(,"package")
## [1] "DOSE"
如何查看某个条目下的所有基因名字，很简单，也是不断的取子集操作：
# 第一个条目的所有基因
gsea_res[[gsea_res$ID[[1]]]]
## [1] "1673" "5068" "2919" "5967" "5968" "1670" "1671" "2920" "3553"
## [10] "6373" "4283" "10563" "6374" "6372" "6283" "2921" "6280" "1755"
## [19] "3627" "5266" "718" "725" "4057" "931" "629" "3426" "6278"
## [28] "28461" "1604" "6347" "54210" "3495" "1380" "1378" "1191" "722"
## [37] "1880" "5199" "5648" "5788" "1235" "9308" "717" "3569" "3500"
## [46] "710" "11005" "5196" "716" "715" "1236" "64127" "5079" "940"
## [55] "3507" "28639" "10417" "2219" "5452" "124976" "2213" "5450" "10578"
## [64] "28912" "3514" "6590" "3119" "6036" "3853" "7124" "4068" "5919"
## [73] "4049" "6406" "8547" "3605" "5473" "3458" "729230" "6480" "240"
## [82] "3123" "4239"
除此之外，dplyr中的很多函数，包括filetr/arrange/select/mutate/slice/summarise，也被重新定义了，支持直接对这些富集结果进行操作，简直不要太方便。
y <- ego %>%
filter(p.adjust < 10^-8) %>%
arrange(Count)
dim(y)
## [1] 27 10
#head(y)
#结果依然是对象！
class(y)
## [1] "enrichResult"
## attr(,"package")
## [1] "DOSE"
其他的就不多演示了，大家发挥自己的想象力即可，无非就是对数据框进行操作~
clusterprofiler的富集结果都是一样使用，所以其他的富集分析结果就不再一一演示了。
参考资料
关于clusterprofiler包的用法，y叔专门写了一本书：https:///biomedical-knowledge-mining-book/
除此之外，在y叔的公众号也有大量关于clusterprofiler包的推文，我这里直接把合集链接放这里：clusterprofiler合集，其中很多推文是对使用者问题的回答，很多推文中说过问题/bug等，现在基本都被修复了。
