MAGeCKFlute - Functional enrichment analysis in MAGeCKFlute
24 October 2023
Abstract
MAGeCKFlute provides three different methods for pathway enrichment analysis, including hypergeometric test (HGT), over-representation test (ORT), and gene set enrichment analysis (GSE). MAGeCKFlute also includes several functions for visualizing enrichment results. This vignette introduces the functions related to pathway enrichment analysis.
Contents
- 0.1 Install and load the packages
- 1 Input data
- 2 Enrichment analysis methods
- 3 Session info
- References
knitr::opts_chunk$set(tidy=FALSE, cache=TRUE,
dev="png", message=FALSE, error=FALSE, warning=TRUE)Citation: if you use MAGeCKFlute in published research, please cite: Binbin Wang, Mei Wang, Wubing Zhang. “Integrative analysis of pooled CRISPR genetic screens using MAGeCKFlute.” Nature Protocols (2019), doi: 10.1038/s41596-018-0113-7.
0.1 Install and load the packages
if(!"MAGeCKFlute" %in% installed.packages()) BiocManager::install("MAGeCKFlute")
if(!"ggplot2" %in% installed.packages()) BiocManager::install("ggplot2")
library(MAGeCKFlute)
library(ggplot2)1 Input data
MAGeCKFlute requires a weighted gene list for enrichment analysis.
file1 = file.path(system.file("extdata", package = "MAGeCKFlute"),
"testdata/rra.gene_summary.txt")
gdata = ReadRRA(file1)
genelist = gdata$Score
names(genelist) = gdata$id
genelist = sort(genelist, decreasing = TRUE)
head(genelist)## Jak1 Stat1 Ifngr2 Ifngr1 B2m H2-D1
## 10.7440 9.0330 5.6097 5.4719 4.3433 4.21102 Enrichment analysis methods
2.1 Hypergeometric test
# Alternative functions EnrichAnalyzer and enrich.HGT.
hgtRes1 = EnrichAnalyzer(genelist[1:200], method = "HGT",
type = "Pathway", organism = "mmu")
head(hgtRes1@result)## ID Description
## BIOCARTA_IFNG_PATHWAY BIOCARTA_IFNG_PATHWAY IFNG PATHWAY
## REACTOME_877312 REACTOME_877312 Regulation of IFNG signaling
## BIOCARTA_IL22BP_PATHWAY BIOCARTA_IL22BP_PATHWAY IL22BP PATHWAY
## BIOCARTA_IFNA_PATHWAY BIOCARTA_IFNA_PATHWAY IFNA PATHWAY
## BIOCARTA_TID_PATHWAY BIOCARTA_TID_PATHWAY TID PATHWAY
## REACTOME_877300 REACTOME_877300 Interferon gamma signaling
## NES pvalue p.adjust GeneRatio BgRatio
## BIOCARTA_IFNG_PATHWAY 18.307949 9.401935e-11 5.048839e-08 5/107 6/9264
## REACTOME_877312 9.714240 2.440253e-05 4.081607e-03 3/107 9/9264
## BIOCARTA_IL22BP_PATHWAY 15.317503 2.440253e-05 4.081607e-03 3/107 9/9264
## BIOCARTA_IFNA_PATHWAY 13.883283 3.040303e-05 4.081607e-03 4/107 18/9264
## BIOCARTA_TID_PATHWAY 9.455281 4.054116e-05 4.354120e-03 4/107 19/9264
## REACTOME_877300 9.714240 9.109707e-05 6.444888e-03 3/107 12/9264
## geneID
## BIOCARTA_IFNG_PATHWAY 16451/16452/20846/15979/15980
## REACTOME_877312 16452/15979/15980
## BIOCARTA_IL22BP_PATHWAY 16451/16452/20846
## BIOCARTA_IFNA_PATHWAY 16451/15975/15976/20846
## BIOCARTA_TID_PATHWAY 16452/22059/15979/15980
## REACTOME_877300 16452/15979/15980
## geneName Count
## BIOCARTA_IFNG_PATHWAY Jak1/Jak2/Stat1/Ifngr1/Ifngr2 5
## REACTOME_877312 Jak2/Ifngr1/Ifngr2 3
## BIOCARTA_IL22BP_PATHWAY Jak1/Jak2/Stat1 3
## BIOCARTA_IFNA_PATHWAY Jak1/Ifnar1/Ifnar2/Stat1 4
## BIOCARTA_TID_PATHWAY Jak2/Trp53/Ifngr1/Ifngr2 4
## REACTOME_877300 Jak2/Ifngr1/Ifngr2 3# hgtRes2 = enrich.HGT(genelist[1:200])
# head(hgtRes2@result)2.2 Over-representation test
The ORT and GSEA are implemented in the R package clusterProfiler (Yu et al. 2012). So please install it prior to the analysis.
if(!"clusterProfiler" %in% installed.packages()){
BiocManager::install("clusterProfiler")
}
library(clusterProfiler)
# Alternative functions EnrichAnalyzer and enrich.ORT.
ortRes1 = EnrichAnalyzer(genelist[1:200], method = "ORT",
type = "KEGG", organism = "mmu")
head(ortRes1@result)## ID
## KEGG_mmu05140 KEGG_mmu05140
## KEGG_mmu05235 KEGG_mmu05235
## KEGG_mmu04612 KEGG_mmu04612
## KEGG_mmu05340 KEGG_mmu05340
## KEGG_mmu05212 KEGG_mmu05212
## KEGG_mmu04658 KEGG_mmu04658
## Description NES
## KEGG_mmu05140 Leishmaniasis 17.491523
## KEGG_mmu05235 PD-L1 expression and PD-1 checkpoint pathway in cancer 17.795521
## KEGG_mmu04612 Antigen processing and presentation 9.038750
## KEGG_mmu05340 Primary immunodeficiency 5.717129
## KEGG_mmu05212 Pancreatic cancer 13.069760
## KEGG_mmu04658 Th1 and Th2 cell differentiation 18.307949
## pvalue p.adjust GeneRatio BgRatio
## KEGG_mmu05140 1.530265e-05 0.001254818 7/54 70/4536
## KEGG_mmu05235 5.495082e-04 0.016338254 6/54 88/4536
## KEGG_mmu04612 6.199206e-04 0.016338254 6/54 90/4536
## KEGG_mmu05340 7.969880e-04 0.016338254 4/54 36/4536
## KEGG_mmu05212 1.927944e-03 0.031618287 5/54 76/4536
## KEGG_mmu04658 3.670334e-03 0.044654889 5/54 88/4536
## geneID
## KEGG_mmu05140 16451/20846/15980/15979/16452/17357/16412
## KEGG_mmu05235 16451/20846/15980/15979/16452/13645
## KEGG_mmu04612 12010/14964/21354/21355/19727/53970
## KEGG_mmu05340 21354/21355/19727/53970
## KEGG_mmu05212 16451/20846/12575/13645/22059
## KEGG_mmu04658 16451/20846/15980/15979/16452
## geneName Count
## KEGG_mmu05140 Jak1/Stat1/Ifngr2/Ifngr1/Jak2/Marcksl1/Itgb1 7
## KEGG_mmu05235 Jak1/Stat1/Ifngr2/Ifngr1/Jak2/Egf 6
## KEGG_mmu04612 B2m/H2-D1/Tap1/Tap2/Rfxank/Rfx5 6
## KEGG_mmu05340 Tap1/Tap2/Rfxank/Rfx5 4
## KEGG_mmu05212 Jak1/Stat1/Cdkn1a/Egf/Trp53 5
## KEGG_mmu04658 Jak1/Stat1/Ifngr2/Ifngr1/Jak2 5# ortRes2 = enrich.ORT(genelist[genelist< -1])
# head(ortRes2@result)2.3 Gene set enrichment analysis
library(clusterProfiler)
# Alternative functions EnrichAnalyzer and enrich.GSE.
gseRes1 = EnrichAnalyzer(genelist, method = "GSEA", type = "Pathway", organism = "mmu")## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (8.16% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.2.4 Visualize enrichment results.
2.4.1 Barplot
require(ggplot2)
df = hgtRes1@result
df$logFDR = -log10(df$p.adjust)
p = BarView(df[1:5,], "Description", 'logFDR')
p = p + labs(x = NULL) + coord_flip()
p
# Or use function barplot from enrichplot package
barplot(hgtRes1, showCategory = 5)
2.4.2 Dot plot
## top: up-regulated pathways;
## bottom: down-regulated pathways
EnrichedView(hgtRes1, top = 5, bottom = 0, mode = 1)
EnrichedView(hgtRes1, top = 5, bottom = 0, mode = 2)
dotplot(hgtRes1, showCategory = 5)
2.4.3 Other visualization functions from enrichplot (Yu 2018).
library(enrichplot)
hgtRes1@result$geneID = hgtRes1@result$geneName
cnetplot(hgtRes1, 5)
heatplot(hgtRes1, showCategory = 5, foldChange = genelist)
tmp <- pairwise_termsim(hgtRes1)
emapplot(tmp, showCategory = 5, layout = "kk")## Warning in emapplot.enrichResult(x, showCategory = showCategory, ...): Use 'layout.params = list(layout = your_value)' instead of 'layout'.
## The layout parameter will be removed in the next version.
2.4.4 Visulization for GSEA enriched categories
# show GSEA results of one pathway
idx = which(gseRes1$NES>0)[1]
gseaplot(gseRes1, geneSetID = idx, title = gseRes1$Description[idx])
# show GSEA results of multiple pathways
gseaplot2(gseRes1, geneSetID = which(gseRes1$NES>0)[1:3])
2.5 Type of gene sets for enrichment analysis
For enrichment analysis, MAGeCKFlute signifies the public available gene sets, including Pathways (PID, KEGG, REACTOME, BIOCARTA, C2CP), GO terms (GOBP, GOCC, GOMF), Complexes (CORUM) and molecular signature from MsigDB (c1, c2, c3, c4, c5, c6, c7, HALLMARK).
2.5.1 Functional enrichment analysis on GO terms and pathways
Analysis of high-throughput data increasingly relies on pathway annotation and functional information derived from Gene Ontology, which is also useful in the analysis of CRISPR screens.
## combination of the gene sets
enrichComb = EnrichAnalyzer(genelist[1:200], type = "KEGG+REACTOME", organism = "mmu")
EnrichedView(enrichComb, top = 5)
## All pathways
enrich = EnrichAnalyzer(geneList = genelist[1:200], type = "REACTOME", organism = "mmu")
EnrichedView(enrich, top = 5)
## All gene ontology
enrichGo = EnrichAnalyzer(genelist[1:200], type = "GOBP", organism = "mmu")2.5.2 Protein complex analysis
Functional annotations from the pathways and GO are powerful in the context of network dynamics. However, the approach has limitations in particular for the analysis of CRISPR screenings, in which elements within a protein complex rather than complete pathways might have a strong selection. So we incorporate protein complex resource from CORUM database, which enable the identification of essential protein complexes from the CRISPR screens.
enrichPro = EnrichAnalyzer(genelist[1:200], type = "Complex", organism = "mmu")
EnrichedView(enrichPro, top = 5)
2.5.3 Enrichment analysis on the combination of the gene sets
2.6 Limit the size of gene sets for testing
Usually, only the core genes in a pathway could be selected in a CRISPR screen, so we recommend to perform enrichment analysis on small pathways instead of pathways involving hundreds or even more genes.
enrich = EnrichAnalyzer(genelist[1:200], type = "GOBP", limit = c(2, 80), organism = "mmu")
EnrichedView(enrich, top = 5)
2.7 Remove redundant results using EnrichedFilter.
The EnrichedFilter function tries to eliminate redundant pathways from the enrichment results by calculating the Jaccard index between pathways.
enrich1 = EnrichAnalyzer(genelist[1:200], type = "Pathway", organism = "mmu")
enrich2 = EnrichedFilter(enrich1)
EnrichedView(enrich1, top = 15)
EnrichedView(enrich2, top = 15)
3 Session info
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] enrichplot_1.22.0 ggplot2_3.4.4 clusterProfiler_4.10.0
## [4] MAGeCKFlute_2.6.0 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.1 later_1.3.1
## [3] bitops_1.0-7 ggplotify_0.1.2
## [5] filelock_1.0.2 tibble_3.2.1
## [7] polyclip_1.10-6 graph_1.80.0
## [9] XML_3.99-0.14 lifecycle_1.0.3
## [11] edgeR_4.0.0 lattice_0.22-5
## [13] MASS_7.3-60 magrittr_2.0.3
## [15] limma_3.58.0 sass_0.4.7
## [17] rmarkdown_2.25 jquerylib_0.1.4
## [19] yaml_2.3.7 httpuv_1.6.12
## [21] depmap_1.15.0 cowplot_1.1.1
## [23] DBI_1.1.3 RColorBrewer_1.1-3
## [25] zlibbioc_1.48.0 purrr_1.0.2
## [27] ggraph_2.1.0 BiocGenerics_0.48.0
## [29] msigdbr_7.5.1 RCurl_1.98-1.12
## [31] yulab.utils_0.1.0 tweenr_2.0.2
## [33] rappdirs_0.3.3 sva_3.50.0
## [35] GenomeInfoDbData_1.2.11 IRanges_2.36.0
## [37] S4Vectors_0.40.0 ggrepel_0.9.4
## [39] tidytree_0.4.5 genefilter_1.84.0
## [41] pheatmap_1.0.12 annotate_1.80.0
## [43] codetools_0.2-19 DOSE_3.28.0
## [45] ggforce_0.4.1 tidyselect_1.2.0
## [47] aplot_0.2.2 farver_2.1.1
## [49] viridis_0.6.4 matrixStats_1.0.0
## [51] stats4_4.3.1 BiocFileCache_2.10.0
## [53] pathview_1.42.0 jsonlite_1.8.7
## [55] ellipsis_0.3.2 tidygraph_1.2.3
## [57] survival_3.5-7 ggnewscale_0.4.9
## [59] tools_4.3.1 treeio_1.26.0
## [61] HPO.db_0.99.2 Rcpp_1.0.11
## [63] glue_1.6.2 gridExtra_2.3
## [65] xfun_0.40 mgcv_1.9-0
## [67] qvalue_2.34.0 MatrixGenerics_1.14.0
## [69] GenomeInfoDb_1.38.0 dplyr_1.1.3
## [71] withr_2.5.1 BiocManager_1.30.22
## [73] fastmap_1.1.1 fansi_1.0.5
## [75] digest_0.6.33 R6_2.5.1
## [77] mime_0.12 gridGraphics_0.5-1
## [79] colorspace_2.1-0 GO.db_3.18.0
## [81] RSQLite_2.3.1 utf8_1.2.4
## [83] tidyr_1.3.0 generics_0.1.3
## [85] data.table_1.14.8 graphlayouts_1.0.1
## [87] httr_1.4.7 scatterpie_0.2.1
## [89] pkgconfig_2.0.3 gtable_0.3.4
## [91] blob_1.2.4 XVector_0.42.0
## [93] shadowtext_0.1.2 htmltools_0.5.6.1
## [95] bookdown_0.36 fgsea_1.28.0
## [97] scales_1.2.1 Biobase_2.62.0
## [99] png_0.1-8 ggfun_0.1.3
## [101] knitr_1.44 reshape2_1.4.4
## [103] nlme_3.1-163 curl_5.1.0
## [105] org.Hs.eg.db_3.18.0 cachem_1.0.8
## [107] stringr_1.5.0 BiocVersion_3.18.0
## [109] parallel_4.3.1 HDO.db_0.99.1
## [111] AnnotationDbi_1.64.0 pillar_1.9.0
## [113] grid_4.3.1 vctrs_0.6.4
## [115] promises_1.2.1 dbplyr_2.3.4
## [117] xtable_1.8-4 Rgraphviz_2.46.0
## [119] evaluate_0.22 KEGGgraph_1.62.0
## [121] magick_2.8.1 cli_3.6.1
## [123] locfit_1.5-9.8 compiler_4.3.1
## [125] rlang_1.1.1 crayon_1.5.2
## [127] labeling_0.4.3 plyr_1.8.9
## [129] fs_1.6.3 stringi_1.7.12
## [131] viridisLite_0.4.2 BiocParallel_1.36.0
## [133] babelgene_22.9 MPO.db_0.99.7
## [135] munsell_0.5.0 Biostrings_2.70.0
## [137] lazyeval_0.2.2 GOSemSim_2.28.0
## [139] Matrix_1.6-1.1 ExperimentHub_2.10.0
## [141] patchwork_1.1.3 bit64_4.0.5
## [143] KEGGREST_1.42.0 statmod_1.5.0
## [145] shiny_1.7.5.1 interactiveDisplayBase_1.40.0
## [147] AnnotationHub_3.10.0 igraph_1.5.1
## [149] memoise_2.0.1 bslib_0.5.1
## [151] ggtree_3.10.0 fastmatch_1.1-4
## [153] bit_4.0.5 ape_5.7-1
## [155] gson_0.1.0References
Yu, Guangchuang. 2018. Enrichplot: Visualization of Functional Enrichment Result. https://github.com/GuangchuangYu/enrichplot.
Yu, Guangchuang, Li-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “ClusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters.” OMICS: A Journal of Integrative Biology 16 (5): 284–87. https://doi.org/10.1089/omi.2011.0118.