{"paper_id":"2c3cedcd-1dcf-428c-8a42-2ea651d83c25","body_text":"Direction-aware functional class scoring\nenrichment analysis of In finium DNA\nmethylation data\nMark Ziemann1,2*, Mandhri Abeysooriya 2,3, Anusuiya Bora1,2, Séverine Lamon 3, Mary\nSravya Kasu2, Mitchell W. Norris 2, Yen Ting Wong 4,5, Jeffrey M. Craig 4,5\n1. Burnet Institute, Melbourne, Australia.\n2. Deakin University, School of Life and Environmental Sciences, Geelong, Australia.\n3. Deakin University, School of Exercise and Nutrition Sciences, Institute for Physical\nActivity and Nutrition, Geelong, Australia.\n4. Deakin University, School of Medicine, Geelong, Australia.\n5. Murdoch Children’s Research Institute, Melbourne, Australia.\n(*) Correspondence: mark.ziemann@burnet.edu.au\nAuthor ORCID\nMark Ziemann 0000-0002-7688-6974\nMandhri Abeysooriya 0000-0003-2163-6203\nAnusuiya Bora 0009-0006-2908-1352\nSéverine Lamon 0000-0002-3271-6551\nMary Sravya Kasu 0000-0002-7891-836X\nMitchell W. Norris 0009-0001-4338-6252\nYen Ting Wong 0000-0003-3120-0964\nJeffrey M Craig 0000-0003-3979-7849\n1\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nAbstract\nIn\nfi\nnium Methylation BeadChip arrays remain one of the most popular platforms for\nepigenome-wide association studies, but tools for downstream pathway analysis have\ntheir limitations. Functional class scoring (FCS) is a group of pathway enrichment\ntechniques that involve the ranking of genes and evaluation of their collective regulation in\nbiological systems, but the implementations described for In\nfi\nnium methylation array data\ndo not retain direction information, which is important for mechanistic understanding of\ngenomic regulation. Here, we evaluate several candidate FCS methods that retain\ndirectional information. According to simulation results, the best-performing method\ninvolves the mean aggregation of probe limma t-statistics by gene followed by a\nrank-ANOVA enrichment test using the mitch package. This method, which we call “LAM”,\noutperformed an existing over-representation analysis method in simulations, and showed\nhigher sensitivity and robustness in an analysis of real lung tumour-normal paired\ndatasets. Using matched RNA-seq data we examine the relationship of methylation\ndifferences at promoters and gene bodies with RNA expression at the level of pathways in\nlung cancer. To demonstrate the utility of our approach, we apply it to three other contexts\nwhere public data were available. Firstly, we examine differential pathway methylation\nassociated with chronological age. Secondly, we investigate pathway methylation\ndifferences in infants conceived with in vitro fertilisation. Lastly, we analyse differential\npathway methylation in 19 disease states, identifying hundreds of novel associations.\nThese results show LAM is a powerful method for the detection of differential pathway\nmethylation as compared to existing methods. A reproducible vignette is provided to\nillustrate how to implement this method.\nKeywords: Pathway analysis, functional enrichment analysis, In\nfi\nnium Array, DNA\nmethylation, epigenetics, epigenome-wide association study.\n2\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nIntroduction\nDNA methylation is the most widely studied epigenetic mark, for its role in development\nand disease.1 Hundreds of epigenome-wide association studies (EWASs) are conducted\neach year to understand DNA methylation patterns in disease and other contexts. 2 In\nfi\nnium\narrays remain the preferred platform for EWASs due to low cost and analytical simplicity\nas compared to high throughput methylation sequencing. 3 In\nfi\nnium arrays typically include\nmultiple probes per gene in different locations including promoters, CpG islands, gene\nbodies and enhancers, 4 which complicates downstream functional enrichment analysis,\notherwise known as pathway analysis or gene set enrichment analysis. Enrichment\nanalysis comes in two popular types: over-representation analysis (ORA) and functional\nclass scoring (FCS). 5 A third type of enrichment analysis called pathway topology\nimproves upon other methods with more sophisticated modelling of gene network\nactivations, but these are yet to be adopted widely. 6 In ORA, genes with probes that meet\nan arbitrary signi\nfi\ncance threshold are selected and compared to a background list of all\ngenes measured in the assay. The test seeks to identify sets of genes (e.g., ontologies,\npathways) that are over-represented in the gene list of interest relative to the background. 7\nFCS takes a different approach by ranking all detected genes by a differential regulation\nscore (e.g., fold change, con\nfi\ndent effect size, t-statistic) followed by a test to assess\nwhether each set of genes has a distribution of scores that is different from the null. 8 The\ngsameth() function of the missMethyl package is the state-of-the-art method for ORA\nof In\nfi\nnium methylation array data as it addresses issues related to probes belonging to\nmore than one gene and the fact that one gene can have multiple probes. 9 As gsameth is\nan ORA method, results strongly depend on the signi\nfi\ncance threshold used, 10 and also on\nthe proportion of probes that meet this threshold. FCS methods ebGSEA and methylGSA\nhave been developed for EWAS data and are suggested to have better sensitivity to\ndetermine subtle associations between pathways and differential methylation. 10,11 The\nebGSEA tool uses an empirical Bayes modeling approach followed by a one-tailed\nnonparametric test for gene set enrichment. 11 The methylGSA package implements two\ndifferent approaches. Firstly, the methylRRA function, which uses robust rank aggregation\nfor detecting differentially methylated genes taking into consideration the variable number\nof probes followed by an FCS test based on z-scores. 10 Secondly, methylglm uses a\nlogistic regression approach to model differential methylation for genes inside and outside\neach gene set. 10 These methods are indeed more sensitive than ORA approaches, but as\nthey are based on one-tailed tests irrespective of the direction of methylation change, it\nlimits the utility for downstream interpretation of the importance of methylation changes\non genomic regulation. Moreover, as pathways tend to be either upregulated or\ndownregulated in omics assays, combining the two tends to dilute their signal and reduce\nsensitivity.12 Neither gsameth nor methylGSA report computed enrichment scores, which\nhinders downstream interpretation, as the enrichment score is a useful surrogate for effect\nsize. Without an enrichment score, users are limited to using statistical signi\nfi\ncance values\n3\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nto prioritise results, which biases towards large gene sets that have non-speci\nfi\nc biological\nfunctions.\nHere we aim to develop and evaluate methods for two-tailed FCS of In\nfi\nnium methylation\narray data that addresses these limitations. Since genes have a variable number of probes,\nthere are many conceivable ways that FCS could be applied to In\nfi\nnium methylation array\ndata. The fi\nrst set of methods (type I) uses differential probe methylation data from\nlimma13 as an input. The differential methylation values of probes belonging to each gene\nare scored and these gene scores undergo an enrichment test. The second set of methods\n(type II) fi\nrst aggregate probe methylation values for each gene before differential\nmethylation analysis with limma. Then, differential gene methylation results are used for a\ndownstream enrichment test. With the above methods, there are several methodological\nchoices, such as the type of aggregation, and the type of downstream test. To determine\nthe best approach to applying FCS to In\nfi\nnium methylation data, in this study, we use\nsimulated data and compare results to existing ORA methods. We then examine the\nsensitivity of selected methods on real cancer data, investigate the association with gene\nexpression and end with additional examples related to human aging, assisted\nreproductive technologies and a large-scale EWAS of 19 diseases.\nMethods\nImplementation overview\nFunctional enrichment analysis is a process of data summarisation from genes to gene\nsets (pathways). This is made more complicated for In\nfi\nnium methylation array data due to\nthe presence of multiple probes per gene, meaning the data needs to be summarised from\nprobes to genes and then to gene sets. Here, we outline eight potential approaches for\nFCS of In\nfi\nnium methylation data implemented in R (v4.3.2), where methods 1-5 are type I\nand 6-8 are type II:\n1. Limma Average t-test (LAT). Differential methylation analysis is conducted at the\nprobe level with limma (v3.58.1), and the limma t-statistics for each gene are\nsummarised (arithmetic mean). The mean t-statistics are used in a downstream\ntwo-sample two-way t-test of gene set enrichment.\n2. Limma Top t-test (L TT). As above, except instead of calculating the mean t-statistic,\nthe probe with the largest magnitude is selected to represent the gene.\n3. Limma Average Wilcox (LAW). Similar to LAT , except instead of the two-sided t-test\nof gene set enrichment, a non-parametric alternative, the Wilcoxon signed-rank test\nis used.\n4\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\n4. Limma Average Mitch (LAM). Similar to LAW, however, instead of the Wilcoxon test,\nthe mitch package 14 (v1.15.0) is used to execute a two-way, two-sample\nANOVA-on-ranks test as described previously 15.\n5. Limma Rank Mitch (LRM). Probes are ranked by t-statistic. The mean rank values\nfor probes belonging to each gene are computed and used as input for a two-way,\ntwo-sample ANOVA-on-ranks test.\n6. Aggregate Limma t-test (AL T). In this method, all probe measurements belonging to\neach gene are averaged before conducting differential methylation analysis with\nlimma. The limma t-statistics are then used downstream for enrichment analysis\nwith the two-sample two-way t-test.\n7. Aggregate Limma Wilcox (ALW). Similar to AL T , however, it uses the Wilcoxon\nnon-parametric test for gene set enrichment.\n8. Aggregate Limma Mitch (ALM). Similar to ALW, however, it uses the Mitch package\nto execute a two-way, two-sample ANOVA-on-ranks test.\nThe above were compared to a standard ORA-based approach (GSA), which involves\nlimma on probes, selection of statistically signi\nfi\ncant probes (FDR<0.05) for ORA with the\ngsameth() function of missMethyl (v1.36.0) which conducts a modi\nfi\ned hypergeometric\ntest that accounts for multiple probe biases. In the case that fewer than 250 signi\nfi\ncant\nprobes were identi\nfi\ned, the 250 probes with the smallest p-values were selected. We\nconducted separate tests for increased and decreased probes, and we speci\nfi\ned the\nbackground as all probes that passed quality control fi\nltering.\nMethod validation using simulations\nTo assess these methods, we adopted a simulated data approach based on the selected\nmodi\nfi\ncation of real methylation data. We downloaded raw intensity EPIC IDAT fi\nles from\nNCBI Gene Expression Omnibus (GEO) for study GSE158422, which consists of lung\ntumour and normal adjacent tissues from 37 patients. 16 For the simulations, only\nnon-cancerous datasets were used. Probe annotations were obtained from the\n“IlluminaHumanMethylationEPICanno.ilm10b4.hg19” Bioconductor package. We randomly\nsampled datasets to serve as control and cases (no replacement), with group sizes varying\nbetween 3 and 12. One thousand random gene sets were created with sizes varying\nbetween 20 and 100 members, with member genes drawn from the annotation set.\nRandom gene sets were used to avoid problems caused by the large overlap among real\ngene sets. Throughout our evaluations, 50 gene sets were selected to be differentially\nmethylated, 25 with increased and 25 with decreased methylation. From these gene sets,\nhalf of the member genes were selected. For those selected genes, half of the annotated\nprobes were selected. We adjusted the M-values in the case group by a speci\nfi\ned amount,\nwhich we call the “delta”, which we varied between 0.1 and 0.5. Following the incorporation\n5\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nof methylation changes to selected probes in the case group, the data underwent limma\ndifferential methylation analysis, followed by enrichment analysis using gsameth, or one\nof the eight FCS methods outlined above. Enrichment results were then fi\nltered for\nFDR<0.05 and the correct direction of methylation change. Differentially methylated gene\nsets observed were compared to gene sets selected to be the ground truth so that the\nnumber of true positives (TP), false positives (FP) and false negatives (FN) could be\ndetermined. At each setting for group size, delta and gene set size, 100 replications were\nconducted, each with a set seed to ensure reproducibility. Simulations were conducted on\na SLURM high-performance computing cluster. 17 Mean TP , FP and FN values were\ndetermined, and then precision, recall and F1 scores were calculated.\nSensitivity analysis\nThe full GSE158422 dataset including normal and cancer samples of all 37 patients\nunderwent differential analysis with limma correcting for patient-of-origin effects.\nDownstream pathway enrichment analysis was conducted using the GSA and LAM\nmethods with Reactome pathways downloaded from MSigDB website (version 2023.1) 18.\nLAM relies on probe-gene associations from the “UCSC_RefGene_Name” column of the\n“IlluminaHumanMethylationEPICanno.ilm10b4.hg19” annotation. As gene annotations are\nseveral years old, we used HGNChelper to update them (database current as of\n21-Dec-2023)19 (3,253 gene symbols were updated). Euler diagrams were created with the\neulerr R package (v7.0.0). To test sensitivity, a random subset of n patients was selected\nfollowed by pathway enrichment analysis with the respective method, with n varying\nbetween 2 and 30. Signi\nfi\ncantly enriched pathways were de\nfi\nned as those with FDR values\n< 0.05. The subset signi\nfi\ncant results were compared to the full group (n=37) results to test\nthe sensitivity of these enrichment methods, and whether fi\nndings at a smaller sample size\nwere consistent with the full group. This process was repeated 50 times for each sample\nsize.\nAssociation of methylation with gene expression\nRNA-seq gene counts for the same set of lung cancer patients were downloaded from GEO\n(accession: GSE158420). HGNChelper was used to fi\nx gene names converted to dates,\nfollowed by DESeq2 differential analysis. This underwent fi\nltering to remove genes with\nexpression below 10 reads per sample on average across the dataset. Then DESeq2 20 was\nused to compare normal and cancerous tissue gene expression taking into consideration\nsample pairing. The differential expression results underwent enrichment analysis using\nthe mitch package 14 with default settings for DESeq2 data tables. CpG sites annotated as\npromoters were considered separately from those located at gene bodies in mitch\nanalysis. For comparison, GSA was used to separately examine pathways using separate\nanalyses of promoter and gene body CpGs. For this analysis, Reactome gene sets were\n6\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nused. Subsequently, to generate the pathway heatmap, we used mitch in multivariate mode\nby analysing differential RNA expression, promoter and gene body methylation together,\nprioritising the results by S distance, a surrogate measure of effect size.\nPathway-level differential methylation in aging\nIn order to demonstrate the utility of the LAM method for large EWAS studies, we\nexamined EPIC methylation pro\nfi\nle associations with age in two independent cohorts\ninvolving a total of 7,036 participants. 21 Limma summary statistics for chronological age\nincluding discovery and replication groups were read into R and enrichment analysis was\nconducted with the LAM method. Brie\nfl\ny, gene-probe associations were fi\nrst obtained from\nthe “IlluminaHumanMethylationEPICanno.ilm10b4.hg19” dataset, followed by an update of\ngene symbols with HGNChelper. Discovery and replication data were imported separately\nwith mitch and multivariate enrichment analysis was conducted with gene sets from\nReactome.\nPathway-level differential methylation in assisted\nreproductive technology\nIn order to demonstrate that methylation patterns across different array platforms can be\ncompared using the LAM method, we examined methylation differences in infants\nassociated with assisted reproduction as described by two independent studies. The Estill\nstudy was published in 2016 and is based on the HM450K array and the data is available\nfrom GEO Accession GSE79257. 22 The Novakovic study in was published 2019, and used\nthe EPIC methylation array and the data is available from GEO under accession\nGSE13143323. Various assisted reproductive technology conception groups are described\nin these studies, but we only examined fresh in vitro fertilisation (IVF) conceived infants\ncompared to naturally conceived infants. Raw intensity fi\nles (IDAT format) were read into R\nwith min\nfi (v1.48.0),24 normalisation was conducted with the SWAN method, 25 and probe\nfi\nltering was conducted to remove probes with a detection p-value > 0.01 in addition to\nprobes located on X or Y chromosomes. M values were computed and underwent\ndifferential methylation analysis using limma. 13 The Estill study included 43 naturally\nconceived infants and 38 with fresh embryo transfer. The Novakovic study included 58\nnaturally conceived infants and 75 with fresh embryo transfer. Limma differential analysis\nwas conducted accounting for sex, then methylation tables underwent LAM analysis with\nReactome gene sets. Probe gene associations for HM450K array were established using\nthe “IlluminaHumanMethylation450kanno.ilmn12.hg19” annotation set and the gene\nsymbols were updated with HGNChelper. As the Novakovic study involved EPIC array data,\ngene-probe associations were obtained from the\n“IlluminaHumanMethylationEPICanno.ilm10b4.hg19” annotation set, and gene symbols\nwere updated with HGNChelper.\n7\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nDifferential pathway methylation in 19 common disease\nstates\nTo examine whether LAM can identify differentially methylated pathways associated with\ndisease prevalence and incidence, we downloaded blood-based EWAS summary statistics\nfrom a recent publication, 26 made available via Zenodo\n(https:/ /zenodo.org/records/8021411). The study used a longitudinal design, with sample\ncollection baseline and diagnosis status established via electronic health record linkage up\nto 14 years after baseline. The study included 18,413 volunteers of European ancestry. We\nused summary statistics from the “full model”, which corrected for relevant potential\nconfounders. Disease status at baseline was used to classify participants into prevalence\ngroups. Disease status at follow-up was used to classify participants into incidence\ngroups. The study used the EPIC In\nfi\nnium array and 752,722 probes were described in the\nsummary statistic dataset. As limma t-statistics were not available, we used probe\ndelta-beta values as input for LAM. LAM was conducted individually for each disease\ncondition with Reactome gene sets. Finally, as a test of type I errors, we randomised each\nof the 19 incident pro\nfi\nles using the base R sample command with a set seed, followed by\nLAM enrichment and counting the number of signi\nfi\ncant pathways. This process was\nrepeated 100 times with unique seeds.\nResults\nEvaluation of FCS methods for In\nfi\nnium methylation array\ndata\nSimulations were conducted using the eight FCS methods (LAT , L TT , LAW, LAM, LRM, AL T ,\nALW, ALM) and GSA, with group sizes of 3-12 and delta values between 0.1 and 0.5, and\nrandom gene sets with sizes 20, 50 and 100. At each parameter setting, the simulations\nwere repeated 100 times with a different seed value. True positives, false positives and\nfalse negatives were used to calculate overall precision and recall at these three gene set\nsizes (Figure 1A). The LAM and ALM methods recorded an overall precision of 0.94, while\nGSA scored lowest with 0.848. Recall was highest for LAM with 0.22, while ALM scored\n0.21 and GSA scored 0.15 F1 scores were highest for LAM with 0.36 followed by ALM with\n0.35, while GSA scored 0.26.\n8\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nFigure 1. Precision, recall and F1 for enrichment tests with simulated data. (A) Overall\nresults including 20, 50 and 100 genes per set. (B-D) Results for simulations using 20, 50\nand 100 genes per set respectively.\nWe noticed that the size of gene sets used in the simulations strongly in\nfl\nuenced the\nresults, with overall recall increasing with gene set size (Figure 1B-D). The GSA method\nperformed better than FCS methods with gene sets of 20 (Figure 1B), while at 50 genes\nper set GSA performance was on par with LAM (Figure 1C). At 100 genes per set, however,\nLAM showed superior precision and recall compared to GSA (Figure 1D).\nFocusing on results from 100 genes per set, LAM precision was more consistent across\nthe parameter ranges as compared to GSA which showed lower precision when group size\nand delta were lower (Figure 2). Recall was strongly dependent on group size and delta\nparameters. LAM recorded relatively higher recall when group size and delta were lower.\nF1 performance scores at 100 genes per set were better for LAM as compared to GSA.\n9\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nFigure 2. Precision, recall and F1 for GSA and LAM enrichment tests with simulated\ndata over a range of group sizes and delta values. One hundred genes per set.\nSensitivity of GSA and LAM methods with real cancer data\nNext, GSA and LAM methods were applied to compare paired normal and cancer\nmethylation data from 37 lung squamous cell carcinoma patients. From the 839,473\nprobes passing quality control, 473,572 were considered signi\nfi\ncant at 5% FDR. Of these,\n119,785 probes were signi\nfi\ncantly increased in cancer, while 353,787 were signi\nfi\ncantly\ndecreased. These signi\nfi\ncantly increased probes mapped to 16,520 unique genes, while\nthe decreased probes mapped to 24,114 genes. These numbers are large, considering that\nthere are 26,219 genes represented on this array. LAM resulted in 406 signi\nfi\ncant\nReactome pathways at FDR<0.05. Of these, 370 involved higher methylation, while 36\ninvolved lower methylation. LAM execution took 51 seconds using eight parallel threads.\nWith GSA, there were 75 Reactome pathways with higher methylation and 32 with lower\nmethylation. 12 pathways were signi\nfi\ncant in both directions. Thirty pathways were\ncommon between LAM and GSA methods. This includes 26 and 4 with higher and lower\n10\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nmethylation respectively (Figure 3A). GSA execution for both directions in series took 22\nseconds on one thread.\nFigure 3. Pathway enrichment analysis comparing normal and cancer samples using\nLAM and GSA methods. (A) Euler diagrams showing the overlap in statistically\nsigni\nfi\ncant pathways (FDR<0.05) from LAM and GSA results. (B) Sizes of signi\nfi\ncant\ngene sets found with LAM and GSA methods. (C) Signi\nfi\ncant LAM pathways with the\nlargest S distance, an enrichment score generated by mitch and effect size proxy. (D)\nTop-ranked signi\nfi\ncant GSA pathways ranked by fold enrichment.\nAs the simulation results indicated LAM had better recall with larger gene sets, we were\ncurious about whether there were differences in the sizes of gene sets found by LAM and\nGSA. GSA signi\nfi\ncant pathways had a median size of 135, while for LAM the median was\n54.5, indicating pathways identi\nfi\ned as signi\nfi\ncant with LAM were collectively smaller than\nGSA (Figure 3B).\nAlthough algorithmic accuracy cannot be inferred from the types of pathways identi\nfi\ned in\nreal data,27 top-ranked results from LAM appeared to be more related to cell differentiation,\nidentity and development (Figure 3C), as compared to GSA (Figure 3D).\n11\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nTo compare the sensitivity of these methods with real data, sample sizes were randomly\ndownsampled, analysed as above and the signi\nfi\ncant pathways were compared to the\nresults from the full group of 37 patients. This process was repeated 50 times at each\nsample size and the results are shown in (Figure 4). Pathways that were identi\nfi\ned as\nstatistically signi\nfi\ncant (FDR<0.05) in the smaller and full group were termed “consistent”,\nwhile pathways that were identi\nfi\ned in the smaller but not in the full group were termed\n“inconsistent”. At a sample size of 10, LAM was able to detect 354 out of 406 (87%)\nconsistent pathways, while GSA detected just 13 of 107 at this sample size (12%). LAM\nidenti\nfi\ned more inconsistent pathways than GSA, but the proportion of inconsistent\nfi\nndings was lower across all sample sizes. At a sample size of 10, the proportion of\ninconsistent pathways was 10% for LAM and 28% for GSA. This result suggests that LAM\nhas superior sensitivity to detect differentially methylated pathways in real cancer data.\n12\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nFigure 4. Sensitivity of LAM and GSA methods with real cancer array data. Upper panels\nshow the number of signi\nfi\ncant pathways identi\nfi\ned in the smaller sample that were\nconsistent with the full group for LAM (left) and GSA (right). Middle panels show the number\nof signi\nfi\ncant pathways identi\nfi\ned in the smaller sample that were inconsistent with the full\ngroup. Lower panels show the proportion of signi\nfi\ncant pathways identi\nfi\ned that were\ninconsistent with the full group. Left panels correspond to LAM, and right panels to GSA\nmethods.\nIntegrating methylation and RNA expression pathways\nAssociations of epigenetic marks with gene expression are of great interest for understanding\ndisease processes. Matching tumour-normal RNA-seq datasets were analysed and we\nconducted enrichment analysis of gene expression together with DNA methylation. From\n18,704 detected genes, there were 12,380 differentially expressed genes (FDR<0.05), with\n7,415 and 4,965 up- and down-regulated in the tumour group respectively. At the pathway level,\nthere were 304 and 135 up- and down-regulated Reactome pathways with altered gene\nexpression (FDR<0.05). As the context of gene methylation is important in in\nfl\nuencing gene\nexpression, promoter and gene body methylation were considered separately. In promoters,\nGSA identi\nfi\ned 15 and 7 pathways with higher and lower methylation respectively, while at gene\nbodies, there were 18 and 75 pathways with higher and lower methylation respectively\n13\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\n(FDR<0.05). Using LAM at promoters, there were 190 and 20 pathways identi\nfi\ned with higher\nand lower methylation respectively. LAM at gene bodies identi\nfi\ned 251 and 31 pathways with\nhigher and lower methylation respectively (FDR<0.05). The relative overlap between signi\nfi\ncant\npromoter and gene body pathways was larger for LAM pathway sets (22%) as compared to\nGSA (3%) (Figure 5A,B). There was no observed overlap between GSA methylation-based\npathways and gene expression pathways (Figure 5A), but there were some overlaps between\npathways identi\nfi\ned with LAM-based promoter and gene body methylation and gene expression\n(Figure 5B).\nFigure 5. Integration of gene expression data with methylation pathways. (A,B) Euler\ndiagrams showing the overlap in statistically signi\nfi\ncant pathways (FDR<0.05) from gene\nexpression and GSA results (A) and LAM (B) results. (C) Multi-dimensional enrichment\nanalysis of promoter and gene body methylation with gene expression using mitch. Top 20\ngene sets shown with the largest S distance after FDR fi\nltering at 0.05.\nUsing LAM, 136/304 upregulated RNA expression pathways were associated with increasing\ngene body methylation. Likewise, 21 pathways with downregulated RNA expression were\nassociated with increasing promoter methylation. Some of these pathways with large\nenrichment scores across the three contrasts are depicted in heatmap form (Figure 5C). There\nwere a small number of pathways that exhibited an inverse relationship between promoter\nmethylation and RNA expression such as “unwinding of DNA”, “mucopolysaccharidoses”,\n“interaction with cumulus cells and the zona pellucida”, “glucocorticoid biosynthesis” and\n14\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\n“olfactory signaling”. These were, however, less common than pathways that showed a positive\nrelationship between RNA expression and promoter methylation. These results demonstrate\nthat LAM can illuminate the complicated relationship between DNA methylation and gene\nexpression at the pathway level in cancer.\nExploration of methylation pathways in aging\nWe sought to explore the utility of LAM for exploring large EWAS datasets. One of the\nbest-studied associations with DNA methylation is with chronological aging, although\npathway-level differential methylation is not well de\nfi\nned. Using limma summary statistics from\nindependent discovery and replication studies using the EPIC array platform described by\nMcCartney et al (2020), 21 we applied the LAM method to explore the pathway level differential\nmethylation. In the discovery cohort, we observed 304 differential pathways (FDR<0.05), while\nin the replication group, we observed 107. There was a high degree of agreement, with 43\npathways signi\nfi\ncant in both groups. This is not surprising, as the gene-level rank differential\nmethylation scores show a strong positive association between discovery and replication\nstudies (Figure 6A). Using the rank-MANOVA test of the mitch package, we identi\nfi\ned 390\npathways with altered methylation. Visualised as a scatterplot, the overall pattern is\nconcordant, although a cluster of points in the lower right of the chart indicates 268 pathways\nwith higher methylation in the discovery sample and lower methylation in the replication\nsample (Figure 6B). When prioritising results by S distance, the enrichment score reported by\nmitch and a surrogate measure of effect size, 29 of the top 30 pathways were concordant in\ntheir direction of regulation (Figure 6C). Some of the observed methylation changes make\nsense with reports in the literature. For example, the complement pathway is associated with\naging28 and in this analysis we observe a reduction in methylation of genes belonging to this\nset including C1R, MBL2, C1S and C1QC (Figure 6D).\n15\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nFigure 6. Pathway-level DNA methylation alterations with chronological age. (A) Contour\nheatmap showing the similarity in gene methylation score ranks in the discovery and\nreplication studies. (B) Mitch pathway enrichment scores (S distances) in discovery and\nreplication studies. Pathways with MANOVA FDR<0.05 are shown in red while others are\nshaded grey. (C) Heatmap of 30 pathways with most extreme S distances (a surrogate of\neffect size). Red indicates increasing methylation and blue shows lower methylation. (D) An\nexample of a pathway identi\nfi\ned with this method, “creation of C4 and C2 activators” shows\nlower methylation of member genes in both discovery and replication studies.\n16\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nExploration of methylation pathway differences in infants\nconceived with in vitro fertilisation\nAnother application for this method is the joint pathway enrichment analysis of studies\nconducted with different array systems. Epigenetic differences between infants conceived\nnaturally and with assisted reproductive technologies have been the subject of studies since\n2009.29 Two of the highest-powered studies with publicly available data are Estill 2016 which\nused the HM450K array and Novakovic 2019 study which used the EPIC array. 22,23 Joint\nenrichment analysis can uncover pathways that are common or different between these\nstudies and illuminate the biological differences between conception groups. For this example,\nwe focused on comparing naturally conceived infants to those conceived with in vitro\nfertilisation (IVF) with fresh embryo transfer. The Estill dataset on the HM450K platform had\n418,833 probes that met the quality criteria (Natural n=43; IVF n=38). Of these, 149,589 probes\nshowed statistically signi\nfi\ncant differential methylation (FDR<0.05), of which 28,767 and\n120,822 exhibited higher or lower methylation in the IVF group respectively. The Novakovic\ndataset on EPIC array had 793,844 probes that met the fi\nltering criteria (n=133). Of these, 5,562\nwere statistically signi\nfi\ncant (FDR<0.05) with 1,348 and 4,214 exhibiting higher or lower\nmethylation in the IVF group respectively. After summarisation to the gene level, the Estill\ndataset on HM450K array describes 19,240 genes, while the Novakovic dataset on EPIC array\ndescribes 22,588 genes. There were 19,234 genes common to both platforms. The rank-rank\nplot of differential gene methylation scores indicates a high degree of similarity overall, and\ninterestingly shows a trend of more genes having lower methylation levels in the IVF group\n(Figure 7A). A scatterplot of pathway enrichment scores indicates a moderate degree of\nagreement between these studies (r=0.35, p=2.2e-16) despite these studies being conducted\non independent cohorts years apart and analysed with different array systems (Figure 7B). A\nheatmap shows a high degree of agreement of pathways with higher methylation in the IVF\ngroups, while pathways with overall lower methylation in the IVF group showed some variability\nin enrichment scores between studies (Figure 7C). We were interested in which pathways had\nconsistently lower methylation in the IVF group and identi\nfi\ned the Reactome pathway\n“Adrenoreceptors” which had enrichment scores of -0.42 in the Estill study and -0.56 in the\nNovakovic study, recording a MANOVA FDR value of 0.04 (Figure 7D). Top-ranked genes from\nthis pathway include ADRA2C, ADRB2, ADRA2B and ADRA1B. A loss of methylation of these\ngenes could contribute to reported elevated blood pressure in people conceived by IVF . 30 Taken\ntogether, this analysis shows that LAM method can be used to compare pathway enrichment\nacross independent studies conducted with HM450K and EPIC arrays.\n17\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nFigure 7. Pathway level DNA methylation differences in natural and IVF conceived infants.\n(A) Contour heatmap showing the similarity in gene methylation score ranks in the Estill\n(HM450K) and Novakovic (EPIC) studies. (B) Mitch pathway enrichment scores (S distances)\nin Estil and Novakovic studies. Pathways with MANOVA FDR<0.05 are shown in red while\nothers are shaded grey. (C) Heatmap of 30 pathways with most extreme S distances (a\nsurrogate of effect size). Red indicates higher methylation and blue shows lower\nmethylation. (D) An example of a pathway identi\nfi\ned with this joint enrichment analysis\nmethod, “Adrenoreceptors” shows lower methylation of member genes in both Estill and\nNovakovic studies.\n18\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nDifferential pathway methylation in 19 common disease\nstates\nTo examine the potential for LAM to reveal differential pathway methylation associations with\ndiseases, we obtained and analysed blood-based EWAS summary statistics from a recent\npublication.26 This study examined methylation in 14 prevalent disease states and the\nincidence of 19 disease states in a group of 18,413 participants, using the EPIC array.\nDelta-beta values were used for the purpose of scoring probe differential methylation, and\nthese values underwent LAM analysis with Reactome gene sets.\nThere were 899 differentially methylated pathways in the prevalence arm of the study\n(FDR<0.05) (Figure 8A). On average, there were 64 pathways with differential methylation in\neach prevalent condition, with Chronic obstructive pulmonary disease (COPD) having the most\n(263) and Alzheimer’s disease having the fewest (3). We selected up to fi\nve statistically\nsigni\nfi\ncant pathways with an absolute S distance of >0.4 in either direction in each prevalent\ncondition to include in a heatmap (Figure 8B). We observed reduced methylation in the\nangiotensinogen metabolism pathway in patients with diabetes. Patients with colorectal\ncancer had lower methylation of the interleukin 10 signaling pathway. Reduced methylation of\nthe salty taste perception pathway was observed in patients with CKD.\n19\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nFigure 8. Differential pathway methylation associated with prevalence of 14 common\ndisease states. (A) A bar plot showing the number of statistically signi\nfi\ncant pathways with\nhigher and lower methylation identi\nfi\ned in each prevalent condition (FDR<0.05). (B) A\nheatmap of S scores for selected pathways across 14 common prevalent disease states.\nStars indicate that the pathway was identi\nfi\ned as being among the top fi\nve differentially\nmethylated pathways in each direction for each condition.\nWe also wanted to know whether the LAM method could identify signatures that appear before\npatients are diagnosed with a condition. Therefore, we examined differential pathway\nmethylation associated with the incidence of 19 disease states. Across all 19 conditions, there\nwere 1570 signi\nfi\ncant pathways (FDR<0.05) (Figure 9A). Liver cirrhosis had the most, with 464,\nfollowed by COPD with 400, while Alzheimer’s had just one and COVID-19 hospitalisation had\nnone. A heatmap shows the enrichment scores for selected pathways across these incident\nconditions (Figure 9B). Breast cancer incidence was associated with reduced methylation to\nthe metal sequestration pathway. Incident prostate cancer was associated with reduced\nmethylation of the FGFR1 pathway. Incident Parkinson’s disease was associated with reduced\nmethylation of aquaporins. In COPD incidence, higher methylation of carbohydrate-binding\nfi\ncolins was observed together with reduced methylation of the GDP mannose synthesis\n20\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\npathway. Incident stroke was associated with higher methylation of the lipoxin synthesis\npathway.\nFigure 9. Differential pathway methylation associated with incidence of 19 common\ndisease states. (A) A bar plot showing the number of statistically signi\nfi\ncant pathways with\nhigher and lower methylation identi\nfi\ned in each incident condition (FDR<0.05). (B) A heatmap\nof S scores for selected pathways across 19 common incident disease states. Stars indicate\nthat the pathway was identi\nfi\ned as being among the top three differentially methylated\npathways in each direction for each condition.\nTo show that these fi\nndings are not the result of false positives from the LAM method, we\nrandomised the incidence pro\nfi\nles for all 19 conditions prior to enrichment analysis. This was\nrepeated 100 times, and these analyses of randomised data yielded few signi\nfi\ncant pathways.\nSpeci\nfi\ncally, from 100 runs 64 yielded no signi\nfi\ncant pathways across 19 conditions.\nTwenty-eight repeats had fewer than fi\nve false positives and eight runs had more than fi\nve\nfalse positives. Across the 100 repeats, the mean number of false positives was 1.95 and the\nmedian was 0. These results indicate that LAM can identify differentially methylated pathways\n21\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nin disease groups and detect altered methylation before a diagnosis has been made, with few\nfalse positives.\nDiscussion\nOverall, we found that the FCS methods performed better across most simulation conditions.\nOf the different FCS methods tested, LAM was superior in terms of precision and recall. LAM\nalso fi\nts relatively easily into existing work\nfl\nows, as most EWAS studies seem to use limma in\ntheir analyses of differential probe methylation. The ALM method had only slightly worse\nperformance as compared to LAM, but ALM requires users to aggregate the probe methylation\nvalues for each sample for each gene before limma. This extra aggregation step is\ncomputationally intensive and not a standard part of existing methylation analysis work\nfl\nows.\nOur analyses also consistently showed that non-parametric methods showed better recall with\nonly a small decrease in precision. Summarisation of probe level t-statistics to gene level only\ntakes a few seconds with the e\nffi\ncient base R aggregate command, and once aggregated,\nsignatures can be analysed with existing packages like mitch. Mitch does not require a large\namount of system memory and due to its parallel architecture can make use of multi-threaded\nprocessors, so it is not much slower than conducting ORA with existing methods that are\nsingle-threaded.\nAn interesting observation is that ORA performed relatively better when gene sets were smaller\nin our simulations. Given that popular pathway databases like KEGG and Reactome pathways\nconsist of both large (>100) and small (<20) gene sets, a hybrid FCS and ORA approach might\nincrease sensitivity and recall.\nAt 100 genes per set, LAM had superior sensitivity compared to GSA, which agrees with a\nprevious report of mitch’s performance with simulated RNA-seq gene expression data 14. The\nimproved sensitivity of the LAM method provides a new opportunity for researchers to\nreanalyse previously conducted EWASs with contemporary pathway databases to better\nunderstand subtle signatures, as we have demonstrated with the chronological age, assisted\nreproduction examples and 19 disease EWAS analyses. The analyses conducted here used\nReactome pathways, but LAM is general-purpose, so different gene sets can be used to\nexamine various hypotheses. Transcription factor target gene sets can be used to identify\ntranscription factors associated with changes in pathway methylation. MicroRNA target gene\nsets can identify potential associations between microRNA targets and DNA methylation\npatterns. The MSigDB resource contains these and several other types of gene sets for\nexploration.30\nThe differences in pathway results obtained with LAM and GSA from lung cancer data were\nstriking, not only in the number of pathways identi\nfi\ned, but also in their relevance to the\ndisease. Interestingly, the relationship between promoter methylation and gene expression was\nnot always inverse as we had anticipated. Joint enrichment analysis showed some pathways\nwith increasing promoter methylation together with increasing RNA expression (e.g., “Negative\n22\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nregulation of AP2 transcription factors”). Likewise, a few pathways showed decreases in\nmethylation with decreases in RNA expression such as “Dectin family”, “TLR regulation” and\n“glycan biosynthesis termination”. These results are in line with previous reports of the\nwidespread existence of both negative and positive relationships between DNA methylation\nand gene expression 31–34, and support the idea that there is a complicated relationship between\nDNA methylation and gene expression in lung cancer.\nWe have demonstrated that this method works with both HM450K and EPIC array data and can\neven facilitate the comparison of EWAS across different platforms. We have provided an\nexample work\nfl\now that researchers can use as a template. It shows how to generate a\nprobe-gene table from Bioconductor annotation datasets, how to update the defunct gene\nnames and how to run the enrichment analysis. Using this example, generating similar gene\ntables for the new In\nfi\nnium Mouse and EPIC 2.0 arrays and any future arrays will be\nstraightforward, once their annotation sets are available through Bioconductor.\nUpdating defunct gene symbols is an important step as 2,788 and 3,253 gene symbols were\nupdated on the HM450K and EPIC array annotations, which represent approximately 10-14% of\ngenes on these platforms. We note that other packages for enrichment of In\nfi\nnium array data\ndo not by default address this, which could lead to a loss of these genes from pathway\nenrichment and diminished sensitivity.\nThe analysis of the large-scale disease EWAS shows that LAM can be readily applied to\nidentify differential pathway methylation associated with prevalent and incident disease. The\nauthors of the study used GOmeth for their pathway analysis and from the 33 models\nexamined, signi\nfi\ncant fi\nndings were obtained for only two, prevalent diabetes and heart\ndisease. In contrast, LAM identi\nfi\ned 2,469 pathway associations (74.8 per model). To\ndemonstrate these are not the results of false positives, permuted pro\nfi\nles yielded on average\n0.10 signi\nfi\ncant fi\nndings per model, supporting the idea that LAM is identifying real pathway\nsignatures.\nThis work could be extended to include not just proximal CpG sites, but also enhancer-based\nprobes and their target genes based on resources such as GeneHancer. 35\nThere is also potential to generate new libraries of differentially methylated probes and gene\nsets, to contribute to the pool of molecular signatures in public resources like MSigDB, which\nwill assist in understanding similarities between methylation pro\nfi\nles.\nData availability\nData sets reanalysed here are publicly available from NCBI GEO, Zenodo or from the\nsupplementary tables of journal articles as described in the Methods section.\n23\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\nCode availability\nCode to generate the results shown here are included in a GitHub repository\n(https:/ /github.com/markziemann/gmea/). To enable reproduction of the fi\nndings, a docker\nimage is made available at DockerHub (https:/ /hub.docker.com/r/mziemann/gmea). A\nreproducible example work\nfl\now report is hosted on our website\n(https:/ /ziemann-lab.net/public/gmea/example_work\nfl\now.html). These materials have been\narchived at Zenodo (https:/ /doi.org/10.5281/zenodo.10685538).\nFunding\nThis work was supported by a grant from the Australian National Health and Medical Research\nCouncil (grant number 1146333 to JMC). Severine Lamon is supported by an Australian\nResearch Council (ARC) Future Fellowship (FT210100278). This research was supported by\nuse of the Nectar Research Cloud, a collaborative Australian research platform supported by\nthe NCRIS-funded Australian Research Data Commons (ARDC). Authors gratefully\nacknowledge the contribution to this work of the Victorian Operational Infrastructure Support\nProgram received by the Burnet Institute.\nAcknowledgements\nWe thank the researchers who publicly shared data to enable this work.\nDisclosure statement\nThe authors report there are no competing interests to declare.\nBibliography\n1. Robertson KD. DNA methylation and human disease. Nat Rev Genet 2005; 6:597–610.\n2. Li M, Zou D, Li Z, Gao R, Sang J, Zhang Y , Li R, Xia L, Zhang T , Niu G, et al. EWAS atlas: A\ncurated knowledgebase of epigenome-wide association studies. Nucleic Acids Res 2019;\n47:D983–8.\n3. Wei S, Tao J, Xu J, Chen X, Wang Z, Zhang N, Zuo L, Jia Z, Chen H, Sun H, et al. Ten years of\nEWAS. Adv Sci (Weinh) 2021; 8.\n4. Pidsley R, Zotenko E, Peters T J, Lawrence MG, Risbridger GP , Molloy P , Van Djik S,\nMuhlhausler B, Stirzaker C, Clark SJ. Critical evaluation of the illumina MethylationEPIC\nBeadChip microarray for whole-genome DNA methylation pro\nfi\nling. Genome Biol 2016; 17.\n24\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\n5. Khatri P , Sirota M, Butte AJ. Ten years of pathway analysis: Current approaches and\noutstanding challenges. PLoS Comput Biol 2012; 8:e1002375.\n6. Xie C, Jauhari S, Mora A. Popularity and performance of bioinformatics software: The case\nof gene set analysis. BMC Bioinformatics 2021; 22.\n7. Tavazoie S, Hughes JD, Campbell MJ, Cho RJ, Church GM. Systematic determination of\ngenetic network architecture. Nat Genet 1999; 22:281–5.\n8. Subramanian A, Tamayo P , Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A,\nPomeroy SL, Golub TR, Lander ES, et al. Gene set enrichment analysis: A knowledge-based\napproach for interpreting genome-wide expression pro\nfi\nles. Proc Natl Acad Sci U S A 2005;\n102:15545–50.\n9. Maksimovic J, Oshlack A, Phipson B. Gene set enrichment analysis for genome-wide DNA\nmethylation data. Genome Biol 2021; 22.\n10. Ren X, Kuan PF . methylGSA: A bioconductor package and shiny app for DNA methylation\ndata length bias adjustment in gene set testing. Bioinformatics 2019; 35:1958–9.\n11. Dong D, Tian Y , Zheng SC, Teschendorff AE. ebGSEA: An improved gene set enrichment\nanalysis method for Epigenome-Wide-Association studies. Bioinformatics 2019; 35:3514–6.\n12. Hong G, Zhang W, Li H, Shen X, Guo Z. Separate enrichment analysis of pathways for up-\nand downregulated genes. J R Soc Interface 2014; 11:20130950.\n13. Ritchie ME, Phipson B, Wu D, Hu Y , Law CW, Shi W, Smyth GK. Limma powers differential\nexpression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;\n43:e47–7.\n14. Kaspi A, Ziemann M. Mitch: Multi-contrast pathway enrichment for multi-omics and\nsingle-cell pro\nfi\nling data. BMC Genomics 2020; 21.\n15. Cox J, Mann M. 1D and 2D annotation enrichment: A statistical method integrating\nquantitative proteomics with complementary high-throughput data. BMC Bioinformatics 2012;\n13.\n16. Cho J-W, Shim HS, Lee CY , Park SY , Hong MH, Lee I, Kim HR. The importance of enhancer\nmethylation for epigenetic regulation of tumorigenesis in squamous lung cancer. Exp Mol Med\n2022; 54:12–22.\n17. Yoo AB, Jette MA, Grondona M. SLURM: Simple linux utility for resource management. In:\nJob scheduling strategies for parallel processing. Berlin, Heidelberg: Springer Berlin\nHeidelberg; 2003. pages 44–60.\n18. Milacic M, Beavers D, Conley P , Gong C, Gillespie M, Griss J, Haw R, Jassal B, Matthews L,\nMay B, et al. The reactome pathway knowledgebase 2024. Nucleic Acids Res 2024;\n52:D672–8.\n25\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\n19. Oh S, Abdelnabi J, Al-Dulaimi R, Aggarwal A, Ramos M, Davis S, Riester M, Waldron L.\nHGNChelper: Identi\nfi\ncation and correction of invalid gene symbols for human and mouse.\nF1000Res 2022; 9:1493.\n20. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for\nRNA-seq data with DESeq2. Genome Biol 2014; 15.\n21. McCartney DL, Zhang F , Hillary RF , Zhang Q, Stevenson AJ, Walker RM, Bermingham ML,\nBoutin T , Morris SW, Campbell A, et al. An epigenome-wide association study of sex-speci\nfi\nc\nchronological ageing. Genome Med 2020; 12.\n22. Estill MS, Bolnick JM, Waterland RA, Bolnick AD, Diamond MP , Krawetz SA. Assisted\nreproductive technology alters deoxyribonucleic acid methylation pro\nfi\nles in bloodspots of\nnewborn infants. Fertil Steril 2016; 106:629–639.e10.\n23. Novakovic B, Lewis S, Halliday J, Kennedy J, Burgner DP , Czajko A, Kim B, Sexton-Oates A,\nJuonala M, Hammarberg K, et al. Assisted reproductive technologies are associated with\nlimited epigenetic variation at birth that largely resolves by adulthood. Nat Commun 2019; 10.\n24. Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP , Hansen KD, Irizarry RA.\nMin\nfi\n: A fl\nexible and comprehensive bioconductor package for the analysis of in\nfi\nnium DNA\nmethylation microarrays. Bioinformatics 2014; 30:1363–9.\n25. Maksimovic J, Gordon L, Oshlack A. SWAN: Subset-quantile within array normalization for\nillumina in\nfi\nnium HumanMethylation450 BeadChips. Genome Biol 2012; 13:R44.\n26. Hillary RF , McCartney DL, Smith HM, Bernabeu E, Gadd DA, Chybowska AD, Cheng Y ,\nMurphy L, Wrobel N, Campbell A, et al. Blood-based epigenome-wide analyses of 19 common\ndisease states: A longitudinal, population-based linked cohort study of 18,413 scottish\nindividuals. PLoS Med 2023; 20:e1004247.\n27. Timmons JA, Szkop KJ, Gallagher IJ. Multiple sources of bias confound functional\nenrichment analysis of global -omics data. Genome Biol 2015; 16.\n28. Katari S, Turan N, Bibikova M, Erinle O, Chalian R, Foster M, Gaughan JP , Coutifaris C,\nSapienza C. DNA methylation and gene expression differences in children conceived in vitro or\nin vivo. Hum Mol Genet 2009; 18:3769–78.\n29. Meister TA, Rimoldi SF , Soria R, Arx R von, Messerli FH, Sartori C, Scherrer U, Rexhaj E.\nAssociation of assisted reproductive technologies with arterial hypertension during\nadolescence. J Am Coll Cardiol 2018; 72:1267–74.\n30. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP , Tamayo P . The molecular\nsignatures database hallmark gene set collection. Cell Syst 2015; 1:417–25.\n26\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint \n\n31. Wagner JR, Busche S, Ge B, Kwan T , Pastinen T , Blanchette M. The relationship between\nDNA methylation, genetic and expression inter-individual variation in untransformed human\nfi\nbroblasts. Genome Biol 2014; 15:R37.\n32. Xu W, Xu M, Wang L, Zhou W, Xiang R, Shi Y , Zhang Y , Piao Y . Integrative analysis of DNA\nmethylation and gene expression identi\nfi\ned cervical cancer-speci\nfi\nc diagnostic biomarkers.\nSignal Transduct Target Ther 2019; 4.\n33. Spainhour JCG, Lim HS, Yi SV , Qiu P . Correlation patterns between DNA methylation and\ngene expression in the cancer genome atlas. Cancer Inform 2019; 18:117693511982877.\n34. Wang K, Dai R, Xia Y , Tian J, Jiao C, Mikhailova T , Zhang C, Chen C, Liu C. Spatiotemporal\nspeci\nfi\ncity of correlated DNA methylation and gene expression pairs across different human\ntissues and stages of brain development. Epigenetics 2022; 17:1110–27.\n35. Fishilevich S, Nudel R, Rappaport N, Hadar R, Plaschkes I, Iny Stein T , Rosen N, Kohn A, Twik\nM, Safran M, et al. GeneHancer: Genome-wide integration of enhancers and target genes in\nGeneCards. Database (Oxford) 2017; 2017.\n27\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted February 24, 2024. ; https://doi.org/10.1101/2024.02.22.581670doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}