Example of using mitch package for Infinium methylation analysis
The GMEA team
2024-04-23
Introduction
Here we are looking at putting together the easiest possible workflow for the package.
Source code: https://github.com/markziemann/gmea/blob/main/example_workflow/example_workflow.Rmd
Requirements
Load packages.
Important: ensure that the mitch version used is 1.15.1 or higher.
suppressPackageStartupMessages({
library("limma")
library("eulerr")
library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
library("IlluminaHumanMethylationEPICanno.ilm10b4.hg19")
library("HGNChelper")
library("tictoc")
library("mitch")
library("kableExtra")
library("beeswarm")
library("missMethyl")
library("gridExtra")
library("png")
})
if( packageVersion("mitch") < "1.15.1") {
warning("This workflow requires mitch version 1.15.1 or higher")
}Load methylation data
The limma1 data table is read in. This dataset is available from GEO (GSE49667)2. The limma data was obtained using a previously prepared workflow.3
dm <- read.csv("DMPs.csv.gz",header=TRUE)
head(dm)## Name logFC AveExpr t P.Value adj.P.Val B
## 1 cg07499259 3.654104 2.46652171 18.73082 7.258963e-08 0.005062817 7.454265
## 2 cg26992245 4.450696 -0.09180715 18.32680 8.603867e-08 0.005062817 7.360267
## 3 cg09747445 -3.337299 -0.25201484 -18.24369 8.913981e-08 0.005062817 7.340431
## 4 cg18808929 -2.990263 0.77522878 -17.90079 1.033329e-07 0.005062817 7.256713
## 5 cg25015733 -3.054336 0.83280190 -17.32537 1.332317e-07 0.005062817 7.109143
## 6 cg21179654 2.859016 1.32460816 17.27702 1.361568e-07 0.005062817 7.096322Note that this object has probes with their own column, not simply as row names.
Load pathways
Gene ontologies were downloaded in GMT format from MSigDB on 15th Jan 20244,5. The GMT file is read into R using the mitch function gmt_import().
gene_sets <- gmt_import("human_12-04-2023.gmt.gz")Curate the annotation
One of the critical parts of this workflow is the establishment of probe-gene relationships. This controls how the probe data is aggregated to make the gene level scores.
In this demonstration we are using all probes on the HM450K chip. Alternative approaches could use only promoters CpGs for example.
As these annotations are several years old, many of the annotated gene names are no longer current. To remedy this, the gene names are screened with the HGNChelper package and any defunct symbols get updated to the newer gene name, so they will be recognised properly in the gene sets.
The process for making a gene table for EPIC data is included down below.
tic()
anno <- getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt1 <- stack(gp2)
colnames(gt1) <- c("gene","probe")
gt1$probe <- as.character(gt1$probe)
dim(gt1)## [1] 407090 2str(gt1)## 'data.frame': 407090 obs. of 2 variables:
## $ gene : chr "TSPY4" "FAM197Y2" "TTTY14" "TMSB4Y" ...
## $ probe: chr "cg00050873" "cg00050873" "cg00212031" "cg00214611" ...length(unique(gt1$gene))## [1] 21231toc() #6.0s## 3.306 sec elapsedtic()
#new.hgnc.table <- getCurrentHumanMap()
new.hgnc.table <- readRDS("../figs/new.hgnc.table.rds")
fix <- checkGeneSymbols(gt1$gene,map=new.hgnc.table)## Warning in checkGeneSymbols(gt1$gene, map = new.hgnc.table): Human gene symbols
## should be all upper-case except for the 'orf' in open reading frames. The case
## of some letters was corrected.## Warning in checkGeneSymbols(gt1$gene, map = new.hgnc.table): x contains
## non-approved gene symbolsfix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
length(unique(fix2$x))## [1] 2788gt1$gene <- fix$Suggested.Symbol
toc() # 36 sec (2788 gene names updated)## 32.316 sec elapsedhead(gt1)## gene probe
## 1 TSPY4 cg00050873
## 2 FAM197Y2 cg00050873
## 3 TTTY14 cg00212031
## 4 TMSB4Y cg00214611
## 5 TBL1Y cg01707559
## 6 TMSB4Y cg02004872Mitch pipeline
The first part is to import the data into mitch.
This only works properly when mitch’s mapGeneIds function uses mean to aggregate (line 69 of mitch.R).
tic()
m <- mitch_import(x=dm,DEtype="limma",geneTable=gt1,geneIDcol="Name")## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 439918## Note: no. genes in output = 20092## Warning in mitch_import(x = dm, DEtype = "limma", geneTable = gt1, geneIDcol = "Name"): Warning: less than half of the input genes are also in the
## outputtoc() #2.5 sec## 2.207 sec elapsedhead(m) %>%
kbl(caption = "Example of differential gene methylation scores used for pathway analysis") %>%
kable_paper("hover", full_width = F)| x | |
|---|---|
| A1BG | 0.7093783 |
| A1BG-AS1 | 0.5498603 |
| A1CF | 1.0133746 |
| A2M | 0.2846473 |
| A2ML1 | 0.1595816 |
| A4GALT | 0.0441000 |
Now run the enrichment analysis.
Note that the results are not sorted by p-value, rather S.distance, an enrichment score. I think this works better for interpretation.
tic()
mres <- mitch_calc(x=m,genesets=gene_sets,minsetsize=5, priority="effect",cores=16)## Note: Enrichments with large effect sizes may not be
## statistically significant.toc() #7 secs## 358.377 sec elapsedmtable <- mres$enrichment_result
up <- subset(mtable,s.dist>0 & p.adjustANOVA<0.05)
dn <- subset(mtable,s.dist<0 & p.adjustANOVA<0.05)
nrow(up)## [1] 2396nrow(dn)## [1] 42059head(up,10) %>%
kbl(caption = "Top significant pathways with higher methylation") %>%
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 114514 | GSE162353-1-vs-2-human up | 5 | 0.0040316 | 0.7425599 | 0.0155825 |
| 118141 | GSE178398-0-vs-3-human up | 5 | 0.0042588 | 0.7380794 | 0.0162838 |
| 68529 | GSE166437,GSE166439-3-vs-4-human up | 5 | 0.0052811 | 0.7202768 | 0.0193669 |
| 108043 | GSE146358-0-vs-4-human dn | 6 | 0.0031908 | 0.6950944 | 0.0128963 |
| 70654 | GSE130614-3-vs-1-human dn | 5 | 0.0072981 | 0.6927963 | 0.0252028 |
| 14000 | GSE182381-2-vs-3-human dn | 6 | 0.0036982 | 0.6842743 | 0.0145423 |
| 101285 | GSE136891-6-vs-1-human up | 5 | 0.0082650 | 0.6819834 | 0.0278761 |
| 107595 | GSE129811-1-vs-3-human up | 5 | 0.0087041 | 0.6774431 | 0.0290616 |
| 121418 | GSE140208-0-vs-3-human dn | 6 | 0.0040946 | 0.6767234 | 0.0157758 |
| 71601 | GSE162179,GSE162180-4-vs-10-human up | 5 | 0.0107308 | 0.6588241 | 0.0343618 |
head(dn,10) %>%
kbl(caption = "Top significant pathways with lower methylation") %>%
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 80802 | GSE54308-1-vs-3-human up | 5 | 0.0013866 | -0.8256186 | 0.0065214 |
| 98727 | GSE173377-0-vs-5-human dn | 5 | 0.0017852 | -0.8066212 | 0.0080255 |
| 91255 | GSE113585,GSE113590-0-vs-3-human up | 5 | 0.0028025 | -0.7717130 | 0.0115973 |
| 110979 | GSE157012-3-vs-6-human up | 5 | 0.0051333 | -0.7226465 | 0.0189353 |
| 57575 | GSE102481-6-vs-7-human up | 7 | 0.0009524 | -0.7211138 | 0.0047933 |
| 29028 | GSE79735,GSE79736-0-vs-3-human dn | 5 | 0.0056134 | -0.7151591 | 0.0203459 |
| 100362 | GSE116256-4-vs-11-human up | 5 | 0.0061945 | -0.7068353 | 0.0220586 |
| 110173 | GSE117192-8-vs-11-human up | 5 | 0.0063844 | -0.7042664 | 0.0226081 |
| 43563 | GSE194125-5-vs-4-human dn | 10 | 0.0001545 | -0.6908874 | 0.0010898 |
| 52987 | GSE124761-4-vs-6-human dn | 6 | 0.0042940 | -0.6731720 | 0.0163914 |
Now make a barplot of these top findings.
top <- rbind(up[1:10,],dn[1:10,])
top <- top[order(top$s.dist),]
barnames <- gsub("_"," ",top$set)
cols <- as.character(((sign(top$s.dist)+1)/2)+1)
cols <- gsub("1","blue",gsub("2","red",cols))
par(mar = c(5.1, 29.1, 4.1, 2.1))
barplot(abs(top$s.dist),horiz=TRUE,las=1,names.arg=barnames, cex.names=0.8,
cex.axis=0.8,col=cols, xlab="Enrichment score",main="Gene Ontologies")
mtext("whole gene")
grid()
par( mar = c(5.1, 4.1, 4.1, 2.1) )There are some interesting results there that give hints to the differences between naive T and Treg cells.
Make a html report and some charts.
mitch_report(res=mres,outfile="rummageo_mitchreport_rummageo.html",overwrite=TRUE)## Dataset saved as " /tmp/RtmpuA2ogm/rummageo_mitchreport_rummageo.rds ".##
##
## processing file: mitch.Rmd## 1/34
## 2/34 [checklibraries]
## 3/34
## 4/34 [peek]
## 5/34
## 6/34 [metrics]
## 7/34
## 8/34 [scatterplot]## 9/34
## 10/34 [contourplot]
## 11/34
## 12/34 [input_geneset_metrics1]
## 13/34
## 14/34 [input_geneset_metrics2]## 15/34
## 16/34 [input_geneset_metrics3]## 17/34
## 18/34 [echart1d]
## 19/34 [echart2d]
## 20/34
## 21/34 [heatmap]
## 22/34
## 23/34 [effectsize]
## 24/34
## 25/34 [results_table]
## 26/34
## 27/34 [results_table_complete]
## 28/34
## 29/34 [detailed_geneset_reports1d]## 30/34
## 31/34 [detailed_geneset_reports2d]
## 32/34
## 33/34 [session_info]
## 34/34## output file: /home/mark.ziemann@domain.internal.burnet.edu.au/projects/gmea/example_workflow/mitch.knit.md## /usr/bin/pandoc +RTS -K512m -RTS /home/mark.ziemann@domain.internal.burnet.edu.au/projects/gmea/example_workflow/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/RtmpuA2ogm/mitch_report.html --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /usr/local/lib/R/site-library/rmarkdown/rmarkdown/lua/latex-div.lua --self-contained --variable bs3=TRUE --section-divs --template /usr/local/lib/R/site-library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /tmp/RtmpuA2ogm/rmarkdown-str2275313c7ccaf8.html##
## Output created: /tmp/RtmpuA2ogm/mitch_report.html## [1] TRUEmitch_plots(res=mres,outfile="rummageo_mitchcharts.pdf")## png
## 2Make gene table for EPIC array data
In case you are working with EPIC array data, the probe-gene table can be make with this snippet of code.
tic()
anno <- getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
myann <- data.frame(anno[,c("UCSC_RefGene_Name","UCSC_RefGene_Group","Islands_Name","Relation_to_Island")])
gp <- myann[,"UCSC_RefGene_Name",drop=FALSE]
gp2 <- strsplit(gp$UCSC_RefGene_Name,";")
names(gp2) <- rownames(gp)
gp2 <- lapply(gp2,unique)
gt <- stack(gp2)
colnames(gt) <- c("gene","probe")
gt$probe <- as.character(gt$probe)
dim(gt)## [1] 684970 2str(gt)## 'data.frame': 684970 obs. of 2 variables:
## $ gene : chr "YTHDF1" "EIF2S3" "PKN3" "CCDC57" ...
## $ probe: chr "cg18478105" "cg09835024" "cg14361672" "cg01763666" ...toc() #9.0s## 6.589 sec elapsedtic()
#new.hgnc.table <- getCurrentHumanMap()
new.hgnc.table <- readRDS("../figs/new.hgnc.table.rds")
fix <- checkGeneSymbols(gt$gene,map=new.hgnc.table)## Warning in checkGeneSymbols(gt$gene, map = new.hgnc.table): Human gene symbols
## should be all upper-case except for the 'orf' in open reading frames. The case
## of some letters was corrected.## Warning in checkGeneSymbols(gt$gene, map = new.hgnc.table): x contains
## non-approved gene symbolsfix2 <- fix[which(fix$x != fix$Suggested.Symbol),]
length(unique(fix2$x))## [1] 3253gt$gene <- fix$Suggested.Symbol
toc()## 44.696 sec elapsedSession information
sessionInfo()## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## 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_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] pkgload_1.3.4
## [2] GGally_2.2.1
## [3] ggplot2_3.5.0
## [4] reshape2_1.4.4
## [5] gplots_3.1.3.1
## [6] gtools_3.9.5
## [7] tibble_3.2.1
## [8] dplyr_1.1.4
## [9] echarts4r_0.4.5
## [10] png_0.1-8
## [11] gridExtra_2.3
## [12] missMethyl_1.36.0
## [13] beeswarm_0.4.0
## [14] kableExtra_1.4.0
## [15] mitch_1.15.6
## [16] tictoc_1.2
## [17] HGNChelper_0.8.1
## [18] IlluminaHumanMethylationEPICanno.ilm10b4.hg19_0.6.0
## [19] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
## [20] minfi_1.48.0
## [21] bumphunter_1.44.0
## [22] locfit_1.5-9.8
## [23] iterators_1.0.14
## [24] foreach_1.5.2
## [25] Biostrings_2.70.2
## [26] XVector_0.42.0
## [27] SummarizedExperiment_1.32.0
## [28] Biobase_2.62.0
## [29] MatrixGenerics_1.14.0
## [30] matrixStats_1.2.0
## [31] GenomicRanges_1.54.1
## [32] GenomeInfoDb_1.38.6
## [33] IRanges_2.36.0
## [34] S4Vectors_0.40.2
## [35] BiocGenerics_0.48.1
## [36] eulerr_7.0.0
## [37] limma_3.58.1
##
## loaded via a namespace (and not attached):
## [1] splines_4.3.3 later_1.3.2
## [3] BiocIO_1.12.0 bitops_1.0-7
## [5] filelock_1.0.3 preprocessCore_1.64.0
## [7] XML_3.99-0.16.1 lifecycle_1.0.4
## [9] lattice_0.22-5 MASS_7.3-60.0.1
## [11] base64_2.0.1 scrime_1.3.5
## [13] magrittr_2.0.3 sass_0.4.8
## [15] rmarkdown_2.25 jquerylib_0.1.4
## [17] yaml_2.3.8 httpuv_1.6.14
## [19] doRNG_1.8.6 askpass_1.2.0
## [21] DBI_1.2.2 RColorBrewer_1.1-3
## [23] abind_1.4-5 zlibbioc_1.48.0
## [25] quadprog_1.5-8 purrr_1.0.2
## [27] RCurl_1.98-1.14 rappdirs_0.3.3
## [29] GenomeInfoDbData_1.2.11 genefilter_1.84.0
## [31] annotate_1.80.0 svglite_2.1.3
## [33] DelayedMatrixStats_1.24.0 codetools_0.2-19
## [35] DelayedArray_0.28.0 xml2_1.3.6
## [37] tidyselect_1.2.0 beanplot_1.3.1
## [39] BiocFileCache_2.10.1 illuminaio_0.44.0
## [41] GenomicAlignments_1.38.2 jsonlite_1.8.8
## [43] multtest_2.58.0 ellipsis_0.3.2
## [45] survival_3.5-8 systemfonts_1.0.5
## [47] tools_4.3.3 progress_1.2.3
## [49] Rcpp_1.0.12 glue_1.7.0
## [51] SparseArray_1.2.4 xfun_0.42
## [53] HDF5Array_1.30.1 withr_3.0.0
## [55] fastmap_1.1.1 rhdf5filters_1.14.1
## [57] fansi_1.0.6 openssl_2.1.1
## [59] caTools_1.18.2 digest_0.6.34
## [61] R6_2.5.1 mime_0.12
## [63] colorspace_2.1-0 biomaRt_2.58.2
## [65] RSQLite_2.3.5 utf8_1.2.4
## [67] tidyr_1.3.1 generics_0.1.3
## [69] data.table_1.15.2 rtracklayer_1.62.0
## [71] prettyunits_1.2.0 httr_1.4.7
## [73] htmlwidgets_1.6.4 S4Arrays_1.2.0
## [75] ggstats_0.5.1 pkgconfig_2.0.3
## [77] gtable_0.3.4 blob_1.2.4
## [79] siggenes_1.76.0 htmltools_0.5.7
## [81] scales_1.3.0 rstudioapi_0.15.0
## [83] knitr_1.45 tzdb_0.4.0
## [85] rjson_0.2.21 nlme_3.1-164
## [87] curl_5.2.1 org.Hs.eg.db_3.18.0
## [89] cachem_1.0.8 rhdf5_2.46.1
## [91] stringr_1.5.1 KernSmooth_2.23-22
## [93] AnnotationDbi_1.64.1 restfulr_0.0.15
## [95] GEOquery_2.70.0 pillar_1.9.0
## [97] grid_4.3.3 reshape_0.8.9
## [99] vctrs_0.6.5 promises_1.2.1
## [101] dbplyr_2.4.0 xtable_1.8-4
## [103] evaluate_0.23 readr_2.1.5
## [105] GenomicFeatures_1.54.3 cli_3.6.2
## [107] compiler_4.3.3 Rsamtools_2.18.0
## [109] rlang_1.1.3 crayon_1.5.2
## [111] rngtools_1.5.2 nor1mix_1.3-2
## [113] mclust_6.0.1 plyr_1.8.9
## [115] stringi_1.8.3 viridisLite_0.4.2
## [117] BiocParallel_1.36.0 munsell_0.5.0
## [119] Matrix_1.6-5 hms_1.1.3
## [121] sparseMatrixStats_1.14.0 bit64_4.0.5
## [123] Rhdf5lib_1.24.2 KEGGREST_1.42.0
## [125] statmod_1.5.0 shiny_1.8.0
## [127] highr_0.10 memoise_2.0.1.9000
## [129] bslib_0.6.1 bit_4.0.5Bibliography
1. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43:e47–7.
2. Zhang Y, Maksimovic J, Naselli G, Qian J, Chopin M, Blewitt ME, Oshlack A, Harrison LC. Genome-wide DNA methylation analysis identifies hypomethylated genes regulated by FOXP3 in human regulatory T cells. Blood 2013; 122:2823–36.
3. Maksimovic J, Phipson B, Oshlack A. A cross-package bioconductor workflow for analysing methylation array data. F1000Res 2017; 5:1281.
4. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database hallmark gene set collection. Cell Syst 2015; 1:417–25.
5. The Gene Ontology Consortium, Aleksander SA, Balhoff J, Carbon S, Cherry JM, Drabkin HJ, Ebert D, Feuermann M, Gaudet P, Harris NL, et al. The gene ontology knowledgebase in 2023. Genetics 2023; 224.