Example of using mitch package for Infinium methylation analysis
The GMEA team
2024-01-15
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.
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")
})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 workflow3.
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("c5.go.v2023.2.Hs.symbols.gmt")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## 6.164 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)## 29.633 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.457 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")## Note: Enrichments with large effect sizes may not be
## statistically significant.toc() #7 secs## 27.94 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] 8nrow(dn)## [1] 10head(up,10) %>%
kbl(caption = "Top significant pathways with higher methylation") %>%
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 10188 | GOMF_STRUCTURAL_CONSTITUENT_OF_SKIN_EPIDERMIS | 35 | 1.53e-05 | 0.4223663 | 0.0155125 |
| 2100 | GOBP_KERATINIZATION | 78 | 4.00e-07 | 0.3329977 | 0.0019071 |
| 2016 | GOBP_INTERMEDIATE_FILAMENT_ORGANIZATION | 73 | 1.20e-06 | 0.3286556 | 0.0041505 |
| 7143 | GOBP_STRIATED_MUSCLE_ADAPTATION | 55 | 3.76e-05 | 0.3212865 | 0.0222716 |
| 8073 | GOCC_KERATIN_FILAMENT | 85 | 3.30e-06 | 0.2919890 | 0.0052689 |
| 2014 | GOBP_INTERMEDIATE_FILAMENT_BASED_PROCESS | 95 | 2.13e-05 | 0.2523268 | 0.0159130 |
| 1964 | GOBP_INORGANIC_ION_IMPORT_ACROSS_PLASMA_MEMBRANE | 136 | 2.14e-05 | 0.2110562 | 0.0159130 |
| 10191 | GOMF_STRUCTURAL_MOLECULE_ACTIVITY | 769 | 1.76e-05 | 0.0911287 | 0.0155125 |
head(dn,10) %>%
kbl(caption = "Top significant pathways with lower methylation") %>%
kable_paper("hover", full_width = F)| set | setSize | pANOVA | s.dist | p.adjustANOVA | |
|---|---|---|---|---|---|
| 9210 | GOMF_GLUTATHIONE_PEROXIDASE_ACTIVITY | 22 | 3.73e-05 | -0.5077864 | 0.0222716 |
| 6144 | GOBP_REGULATION_OF_PLATELET_AGGREGATION | 27 | 3.86e-05 | -0.4575555 | 0.0222716 |
| 7750 | GOCC_CENTROSOME | 630 | 1.79e-05 | -0.1002339 | 0.0155125 |
| 8135 | GOCC_MICROTUBULE_ORGANIZING_CENTER | 829 | 2.30e-06 | -0.0967837 | 0.0052689 |
| 8236 | GOCC_NUCLEAR_PROTEIN_CONTAINING_COMPLEX | 1184 | 0.00e+00 | -0.0944773 | 0.0004838 |
| 4779 | GOBP_POSITIVE_REGULATION_OF_TRANSCRIPTION_BY_RNA_POLYMERASE_II | 1226 | 3.50e-06 | -0.0789663 | 0.0052689 |
| 1182 | GOBP_DNA_METABOLIC_PROCESS | 1016 | 2.82e-05 | -0.0778273 | 0.0195531 |
| 8132 | GOCC_MICROTUBULE_CYTOSKELETON | 1339 | 3.50e-06 | -0.0757001 | 0.0052689 |
| 743 | GOBP_CELL_CYCLE_PROCESS | 1287 | 1.43e-05 | -0.0721735 | 0.0155125 |
| 736 | GOBP_CELL_CYCLE | 1791 | 4.20e-06 | -0.0657249 | 0.0055179 |
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="example_mitchreport.html",overwrite=TRUE)## Note: overwriting existing report## Dataset saved as " /tmp/Rtmp9fpDNY/example_mitchreport.rds ".##
##
## processing file: mitch.Rmd##
|
| | 0%
|
|. | 3%
|
|.. | 6% [checklibraries]
|
|.... | 9%
|
|..... | 12% [peek]
|
|...... | 15%
|
|....... | 18% [metrics]
|
|........ | 21%
|
|.......... | 24% [scatterplot]
|
|........... | 26%
|
|............ | 29% [contourplot]
|
|............. | 32%
|
|.............. | 35% [input_geneset_metrics1]
|
|................ | 38%
|
|................. | 41% [input_geneset_metrics2]
|
|.................. | 44%
|
|................... | 47% [input_geneset_metrics3]
|
|.................... | 50%
|
|...................... | 53% [echart1d]
|
|....................... | 56% [echart2d]
|
|........................ | 59%
|
|......................... | 62% [heatmap]
|
|........................... | 65%
|
|............................ | 68% [effectsize]
|
|............................. | 71%
|
|.............................. | 74% [results_table]
|
|............................... | 76%
|
|................................. | 79% [results_table_complete]
|
|.................................. | 82%
|
|................................... | 85% [detailed_geneset_reports1d]
|
|.................................... | 88%
|
|..................................... | 91% [detailed_geneset_reports2d]
|
|....................................... | 94%
|
|........................................ | 97% [session_info]
|
|.........................................| 100% ## output file: /home/rstudio/gmea/example_workflow/mitch.knit.md## /usr/lib/rstudio-server/bin/quarto/bin/tools/pandoc +RTS -K512m -RTS /home/rstudio/gmea/example_workflow/mitch.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /tmp/Rtmp9fpDNY/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 --embed-resources --standalone --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/Rtmp9fpDNY/rmarkdown-str1ae1d1e4150.html##
## Output created: /tmp/Rtmp9fpDNY/mitch_report.html## [1] TRUEmitch_plots(res=mres,outfile="example_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.497 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()## 42.326 sec elapsedSession information
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 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: Etc/UTC
## 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.3
## [2] GGally_2.2.0
## [3] ggplot2_3.4.4
## [4] reshape2_1.4.4
## [5] gplots_3.1.3
## [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.3.4
## [15] mitch_1.15.0
## [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.1
## [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.2
## [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.2 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 lifecycle_1.0.4
## [9] lattice_0.22-5 MASS_7.3-60
## [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.13
## [19] doRNG_1.8.6 askpass_1.2.0
## [21] DBI_1.1.3 RColorBrewer_1.1-3
## [23] abind_1.4-5 zlibbioc_1.48.0
## [25] rvest_1.0.3 quadprog_1.5-8
## [27] purrr_1.0.2 RCurl_1.98-1.13
## [29] rappdirs_0.3.3 GenomeInfoDbData_1.2.11
## [31] genefilter_1.84.0 annotate_1.80.0
## [33] svglite_2.1.3 DelayedMatrixStats_1.24.0
## [35] codetools_0.2-19 DelayedArray_0.28.0
## [37] xml2_1.3.6 tidyselect_1.2.0
## [39] beanplot_1.3.1 BiocFileCache_2.10.1
## [41] webshot_0.5.5 illuminaio_0.44.0
## [43] GenomicAlignments_1.38.0 jsonlite_1.8.8
## [45] multtest_2.58.0 ellipsis_0.3.2
## [47] survival_3.5-7 systemfonts_1.0.5
## [49] tools_4.3.2 progress_1.2.3
## [51] Rcpp_1.0.11 glue_1.6.2
## [53] SparseArray_1.2.2 xfun_0.41
## [55] HDF5Array_1.30.0 withr_2.5.2
## [57] fastmap_1.1.1 rhdf5filters_1.14.1
## [59] fansi_1.0.6 openssl_2.1.1
## [61] caTools_1.18.2 digest_0.6.33
## [63] R6_2.5.1 mime_0.12
## [65] colorspace_2.1-0 biomaRt_2.58.0
## [67] RSQLite_2.3.4 utf8_1.2.4
## [69] tidyr_1.3.0 generics_0.1.3
## [71] data.table_1.14.10 rtracklayer_1.62.0
## [73] prettyunits_1.2.0 httr_1.4.7
## [75] htmlwidgets_1.6.4 S4Arrays_1.2.0
## [77] ggstats_0.5.1 pkgconfig_2.0.3
## [79] gtable_0.3.4 blob_1.2.4
## [81] siggenes_1.76.0 htmltools_0.5.7
## [83] scales_1.3.0 knitr_1.45
## [85] rstudioapi_0.15.0 tzdb_0.4.0
## [87] rjson_0.2.21 nlme_3.1-163
## [89] curl_5.2.0 org.Hs.eg.db_3.18.0
## [91] cachem_1.0.8 rhdf5_2.46.1
## [93] stringr_1.5.1 KernSmooth_2.23-22
## [95] AnnotationDbi_1.64.1 restfulr_0.0.15
## [97] GEOquery_2.70.0 pillar_1.9.0
## [99] grid_4.3.2 reshape_0.8.9
## [101] vctrs_0.6.5 promises_1.2.1
## [103] dbplyr_2.4.0 xtable_1.8-4
## [105] evaluate_0.23 readr_2.1.4
## [107] GenomicFeatures_1.54.1 cli_3.6.2
## [109] compiler_4.3.2 Rsamtools_2.18.0
## [111] rlang_1.1.2 crayon_1.5.2
## [113] rngtools_1.5.2 nor1mix_1.3-2
## [115] mclust_6.0.1 plyr_1.8.9
## [117] stringi_1.8.3 viridisLite_0.4.2
## [119] BiocParallel_1.36.0 munsell_0.5.0
## [121] Matrix_1.6-1.1 hms_1.1.3
## [123] sparseMatrixStats_1.14.0 bit64_4.0.5
## [125] Rhdf5lib_1.24.1 KEGGREST_1.42.0
## [127] statmod_1.5.0 shiny_1.8.0
## [129] highr_0.10 memoise_2.0.1
## [131] 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.