Implementation overview

Functional enrichment analysis is a process of data summarisation from genes to gene sets (pathways). This is made more complicated for Infinium methylation array data due to the presence of multiple probes per gene, meaning the data needs to be summarised from probes to genes and then to gene sets. Here, we outline eight potential approaches for FCS of Infinium methylation data implemented in R (v4.3.2), where methods 1-5 are type I and 6-8 are type II:

1. Limma Average t-test (LAT). Differential methylation analysis is conducted at the probe level with limma (v3.58.1), and the limma t-statistics for each gene are summarised (arithmetic mean). If a gene has an equal measure of positive and negative t-statistics, the net result will be zero. The mean t-statistics are used in a downstream two-sample two-way t-test of gene set enrichment.

2. Limma Top t-test (LTT). As above, except instead of calculating the mean t-statistic, the probe with the largest magnitude is selected to represent the gene.

3. Limma Average Wilcox (LAW). Similar to LAT, except instead of the two-sided t-test of gene set enrichment, a non-parametric alternative, the Wilcoxon signed-rank test is used.

4. Limma Average Mitch (LAM). Similar to LAW, however, instead of the Wilcoxon test, the mitch package¹⁴ (v1.15.0) is used to execute a two-way, two-sample ANOVA-on-ranks test as described previously¹⁵. Mitch was selected as it performs favourably in terms of accuracy, is extensible to multiple contrasts and has advanced visualisation options.¹⁴

5. Limma Rank Mitch (LRM). Probes are ranked by t-statistic. The mean rank values for probes belonging to each gene are computed and used as input for a two-way, two-sample ANOVA-on-ranks test.

6. Aggregate Limma t-test (ALT). In this method, all probe measurements (M-values) belonging to each gene are averaged before conducting differential methylation analysis with limma. The limma t-statistics are then used downstream for enrichment analysis with the two-sample two-way t-test.

7. Aggregate Limma Wilcox (ALW). Similar to ALT, however, it uses the Wilcoxon non-parametric test for gene set enrichment.

8. Aggregate Limma Mitch (ALM). Similar to ALW, however, uses the Mitch package to execute a two-way, two-sample ANOVA-on-ranks test.

The above were compared to a standard ORA-based approach (GSA), which involves limma on probes, selection of statistically significant probes (FDR<0.05) for ORA with the gsameth() function of missMethyl (v1.36.0) which conducts a modified hypergeometric test that accounts for multiple probe biases. In the case that fewer than 250 significant probes were identified, the 250 probes with the smallest p-values were selected. We conducted separate tests for increased and decreased probes, and we specified the background as all probes that passed quality control filtering.

We didn’t benchmark our methods against ebGSEA and methylGSA methods because they are non-directional tests, and so have a different null hypothesis as compared to LAM.

Method validation using simulations

To assess these methods, we adopted a simulated data approach based on the selected modification of real methylation data. We downloaded raw intensity EPIC IDAT files from NCBI Gene Expression Omnibus (GEO) for study GSE158422, which consists of lung tumour and normal adjacent tissues from 37 patients.¹⁶ For the simulations, only non-cancerous datasets were used. Probe annotations were obtained from the “IlluminaHumanMethylationEPICanno.ilm10b4.hg19” Bioconductor package. We randomly sampled datasets to serve as control and cases (no replacement), with group sizes varying between 3 and 12. One thousand random gene sets were created with sizes varying between 20 and 100 members, with member genes drawn from the annotation set. Random gene sets were used to avoid problems caused by the large overlap among real gene sets. Throughout our evaluations, 50 gene sets were selected to be differentially methylated, 25 with increased and 25 with decreased methylation. From these gene sets, half of the member genes were selected. For those selected genes, half of the annotated probes were selected. We adjusted the M-values in the case group by a specified amount, which we call the “delta”, which we varied between 0.1 and 0.5. Following the incorporation of methylation changes to selected probes in the case group, the data underwent limma differential methylation analysis, followed by enrichment analysis using gsameth, or one of the eight FCS methods outlined above. Enrichment results were then filtered for FDR<0.05 and the same direction of methylation change as the ground truth. Differentially methylated gene sets observed were compared to gene sets selected to be the ground truth, so that the number of true positives (TP), false positives (FP) and false negatives (FN) could be determined. At each setting for group size, delta and gene set size, 100 replications were conducted, each with a set seed to ensure reproducibility. Simulations were conducted on a SLURM high-performance computing cluster.¹⁷ Mean TP, FP and FN values were determined, and then precision, recall and F1 scores were calculated.

Sensitivity analysis

The full GSE158422 dataset including normal and cancer samples of all 37 patients underwent differential analysis with limma correcting for patient-of-origin effects. Downstream pathway enrichment analysis was conducted using the GSA and LAM methods with Reactome pathways downloaded from MSigDB website (version 2023.1)¹⁸. LAM relies on probe-gene associations from the “UCSC_RefGene_Name” column of the “IlluminaHumanMethylationEPICanno.ilm10b4.hg19” annotation. As gene annotations are several years old, we used HGNChelper (v0.8.1) to update them (database current as of 21-Dec-2023)¹⁹ (3,253 gene symbols were updated). Euler diagrams were created with the eulerr R package (v7.0.0). To test sensitivity, a random subset of n patients was selected followed by pathway enrichment analysis with the respective method, with n varying between 2 and 30. Significantly enriched pathways were defined as those with FDR values < 0.05. The subset significant results were compared to the full group (n=37) results to test the sensitivity of these enrichment methods, and whether findings at a smaller sample size were consistent with the full group. This process was repeated 50 times for each sample size.

Association of methylation with gene expression

RNA-seq gene counts for the same set of lung cancer patients were downloaded from GEO (accession: GSE158420). HGNChelper was used to fix gene names converted to dates, followed by DESeq2 differential analysis. This underwent filtering to remove genes with expression below 10 reads per sample on average across the dataset. Then DESeq2²⁰ was used to compare normal and cancerous tissue gene expression taking into consideration sample pairing. The differential expression results underwent enrichment analysis using the mitch package¹⁴ with default settings for DESeq2 data tables. CpG sites annotated as promoters (±1.5 kb around transcription start sites) were considered separately from those located at gene bodies (transcription start site to transcription termination site) with LAM. For comparison, GSA was used to separately examine pathways using separate analyses of promoter and gene body CpGs. For this analysis, Reactome gene sets were used. Subsequently, to generate the pathway heatmap, we used mitch in multivariate mode by analysing differential RNA expression, promoter and gene body methylation together, prioritising the results by S distance, the enrichment score provided by mitch.

Pathway-level differential methylation in aging

In order to demonstrate the utility of the LAM method for large EWAS studies, we examined EPIC methylation profile associations with age in two independent cohorts involving a total of 7,036 participants.²¹ Limma summary statistics for chronological age including discovery and replication groups were read into R and enrichment analysis was conducted with the LAM method. Gene-probe associations were obtained and gene symbols updated as above. Discovery and replication data were imported separately with mitch and multivariate enrichment analysis was conducted with gene sets from Reactome.

Pathway-level differential methylation in assisted reproductive technology

In order to demonstrate that methylation patterns across different array platforms can be compared using the LAM method, we examined methylation differences in infants associated with assisted reproduction as described by two independent studies. The Estill study was published in 2016 and is based on the HM450K array and the data is available from GEO accession number GSE79257.²² The Novakovic study in was published 2019, and used the EPIC methylation array and the data is available from GEO under accession number GSE131433²³. Various assisted reproductive technology conception groups are described in these studies, but we only examined fresh in vitro fertilisation (IVF) conceived infants compared to naturally conceived infants. Raw intensity files (IDAT format) were read into R with minfi²⁴ (v1.48.0), normalisation was conducted with the SWAN method,²⁵ and probe filtering was conducted to remove proves with a detection p-value > 0.01 in addition to probes located on X or Y chromosomes. M values were computed and underwent differential methylation analysis using limma.¹³ The Estill study included 43 naturally conceived infants and 38 with fresh embryo transfer. The Novakovic study included 58 naturally conceived infants and 75 with fresh embryo transfer. Limma differential analysis was conducted accounting for sex, then methylation tables underwent LAM analysis with Reactome gene sets. Probe gene associations were obtained from “IlluminaHumanMethylation450kanno.ilmn12.hg19” and “IlluminaHumanMethylationEPICanno.ilm10b4.hg19” Bioconductor packages for HM450K and EPIC array data respectively and gene symbols updated as above.

Differential pathway methylation in 19 common disease states

To examine whether LAM can identify differentially methylated pathways associated with disease prevalence and incidence, we downloaded blood-based EWAS summary statistics from a recent publication,²⁶ made available via Zenodo (https://zenodo.org/records/8021411). The study used a longitudinal design, with sample collection baseline and diagnosis status established via electronic health record linkage up to 14 years after baseline. The study included 18,413 volunteers of European ancestry. We used summary statistics from the “full model”, which corrected for relevant potential confounders. Disease status at baseline was used to classify participants into prevalence groups. Disease status at follow-up was used to classify participants into incidence groups. The study used the EPIC Infinium array and 752,722 probes were described in the summary statistic dataset. As limma t-statistics were not available, we used probe delta-beta values as input for LAM. LAM was conducted individually for each disease condition with Reactome gene sets. Finally, as a test of type I errors, we randomised each of the 19 incident profiles using the base R sample command with a set seed, followed by LAM enrichment and counting the number of significant pathways. This process was repeated 100 times with unique seeds.

Evaluation of FCS methods for Infinium methylation array data

Simulations were conducted using the eight FCS methods (LAT, LTT, LAW, LAM, LRM, ALT, ALW, ALM) and GSA, with group sizes of 3-12 and delta values between 0.1 and 0.5, and random gene sets with sizes 20, 50 and 100. At each parameter setting, the simulations were repeated 100 times with a different seed value. True positives, false positives and false negatives were used to calculate overall precision and recall at these three gene set sizes (Figure 1A). The LAM and ALM methods recorded an overall precision of 0.94, while GSA scored lowest with 0.848. Recall was highest for LAM with 0.22, while ALM scored 0.21 and GSA scored 0.15 F1 scores were highest for LAM with 0.36 followed by ALM with 0.35, while GSA scored 0.26.

We noticed that the size of gene sets used in the simulations strongly influenced the results, with overall recall increasing with gene set size (Figure 1B-D). The GSA method performed better than FCS methods with gene sets of 20 (Figure 1B), while at 50 genes per set GSA performance was on par with LAM (Figure 1C). At 100 genes per set, however, LAM showed superior precision and recall compared to GSA (Figure 1D).

Figure 1. Precision, recall and F1 for enrichment tests with simulated data. (A) Overall results including 20, 50 and 100 genes per set. (B-D) Results for simulations using 20, 50 and 100 genes per set respectively.

Focusing on results from 100 genes per set, LAM precision was more consistent across the parameter ranges as compared to GSA which showed lower precision when group size and delta were lower (Figure 2). Recall was strongly dependent on group size and delta parameters. LAM recorded relatively higher recall when group size and delta were lower. F1 performance scores at 100 genes per set were better for LAM as compared to GSA.

Figure 2. Precision, recall and F1 for GSA and LAM enrichment tests with simulated data over a range of group sizes and delta values. One hundred genes per set.

Sensitivity of GSA and LAM methods with real cancer data

Next, GSA and LAM methods were applied to compare paired normal and cancer methylation data from 37 lung squamous cell carcinoma patients. From the 839,473 probes passing quality control, 473,572 were considered significant at 5% FDR. Of these, 119,785 probes were significantly increased in cancer, while 353,787 were significantly decreased. These significantly increased probes mapped to 16,520 unique genes, while the decreased probes mapped to 24,114 genes. These numbers are large, considering that there are 26,219 genes represented on this array. LAM resulted in 406 significant Reactome pathways at FDR<0.05. Of these, 370 involved higher methylation, while 36 involved lower methylation. LAM execution took 51 seconds using eight parallel threads. With GSA, there were 75 Reactome pathways with higher methylation and 32 with lower methylation. 12 pathways were significant in both directions. Thirty pathways were common between LAM and GSA methods. This includes 26 and 4 with higher and lower methylation respectively (Figure 3A). GSA execution for both directions in series took 22 seconds on one thread.

As the simulation results indicated LAM had better recall with larger gene sets, we were curious about whether there were differences in the sizes of gene sets found by LAM and GSA. GSA significant pathways had a median size of 135, while for LAM the median was 54.5, indicating pathways identified as significant with LAM were collectively smaller than GSA (Figure 3B).

Although algorithmic accuracy cannot be inferred from the types of pathways identified in real data,²⁷ top-ranked results from LAM appeared to be more related to cell differentiation, identity and development (Figure 3C), as compared to GSA (Figure 3D).

As the number of probes per gene is thought to bias downstream enrichment analysis, we investigated the average probe number for gene sets that were identified as statistically significant or not. For LAM, the mean probe number for significant sets was 6% smaller than non-significant sets (32.9 versus 35.2) (Figure 3E). For GSA the probe number for significant sets was 24% larger than non-significant sets (45.9 versus 36.9) (Figure 3F). This result indicates that LAM is not adversely biased by the number of probes per gene.

Figure 3. Pathway enrichment analysis comparing normal and cancer samples using LAM and GSA methods. (A) Euler diagrams showing the overlap in statistically significant pathways (FDR<0.05) from LAM and GSA results. (B) Sizes of significant gene sets found with LAM and GSA methods. (C) Significant LAM pathways with the largest absolute enrichment score. (D) Top-ranked significant GSA pathways ranked by fold enrichment. (E) Mean number of probes per gene for statistically significant (FDR<0.05) and non-significant (FDR>0.05) gene sets identified with LAM. (F) Mean number of probes per gene for statistically significant (FDR<0.05) and non-significant (FDR>0.05) gene sets identified with GSA.

To compare the sensitivity of these methods with real data, sample sizes were randomly downsampled, analysed as above and the significant pathways were compared to the results from the full group of 37 patients. This process was repeated 50 times at each sample size and the results are shown in (Figure 4). Pathways that were identified as statistically significant (FDR<0.05) in the smaller and full group were termed “consistent”, while pathways that were identified in the smaller but not in the full group were termed “inconsistent”. At a sample size of 10, LAM was able to detect 354 out of 406 (87%) consistent pathways, while GSA detected just 13 of 107 at this sample size (12%). LAM identified more inconsistent pathways than GSA, but the proportion of inconsistent findings was lower across all sample sizes. At a sample size of 10, the proportion of inconsistent pathways was 10% for LAM and 28% for GSA. This result suggests that LAM has superior sensitivity to detect differentially methylated pathways in real cancer data.

Figure 4. Sensitivity of LAM and GSA methods with real cancer array data. Upper panels show the number of significant pathways identified in the smaller sample that were consistent with the full group for LAM (left) and GSA (right). Middle panels show the number of significant pathways identified in the smaller sample that were inconsistent with the full group. Lower panels show the proportion of significant pathways identified that were inconsistent with the full group. Left panels correspond to LAM, and right panels to GSA methods.

Integrating methylation and RNA expression pathways

Associations of epigenetic marks with gene expression are of great interest for understanding disease processes. Matching tumour-normal RNA-seq datasets were analysed and we conducted enrichment analysis of gene expression together with DNA methylation. From 18,704 detected genes, there were 12,380 differentially expressed genes (FDR<0.05), with 7,415 and 4,965 up- and down-regulated in the tumour group respectively. At the pathway level, there were 304 and 135 up- and down-regulated Reactome pathways with altered gene expression (FDR<0.05). As the context of gene methylation is important in influencing gene expression, promoter and gene body methylation were considered separately. In promoters, GSA identified 15 and 7 pathways with higher and lower methylation respectively, while at gene bodies, there were 18 and 75 pathways with higher and lower methylation respectively (FDR<0.05). Using LAM at promoters, there were 190 and 20 pathways identified with higher and lower methylation respectively. LAM at gene bodies identified 251 and 31 pathways with higher and lower methylation respectively (FDR<0.05). The relative overlap between significant promoter and gene body pathways was larger for LAM pathway sets (22%) as compared to GSA (3%) (Figure 5A,B). There was no observed overlap between GSA methylation-based pathways and gene expression pathways (Figure 5A), but there were some overlaps between pathways identified with LAM-based promoter and gene body methylation and gene expression (Figure 5B).

Using LAM, 136/304 upregulated RNA expression pathways were associated with increasing gene body methylation. Likewise, 21 pathways with downregulated RNA expression were associated with increasing promoter methylation. Some of these pathways with large enrichment scores across the three contrasts are depicted in heatmap form (Figure 5C). There were a small number of pathways that exhibited an inverse relationship between promoter methylation and RNA expression such as “unwinding of DNA”, “mucopolysaccharidoses”, “interaction with cumulus cells and the zona pellucida”, “glucocorticoid biosynthesis” and “olfactory signaling”. These were, however, less common than pathways that showed a positive relationship between RNA expression and promoter methylation. These results demonstrate that LAM can illuminate the complicated relationship between DNA methylation and gene expression at the pathway level in cancer.

Figure 5. Integration of gene expression data with methylation pathways. (A,B) Euler diagrams showing the overlap in statistically significant pathways (FDR<0.05) from gene expression and GSA results (A) and LAM (B) results. (C) Multi-dimensional enrichment analysis of promoter and gene body methylation with gene expression using mitch. Top 20 gene sets shown with the largest absolute enrichment scores after FDR filtering at 0.05.

Exploration of methylation pathways in aging

We sought to explore the utility of LAM for exploring large EWAS datasets. One of the best-studied associations with DNA methylation is with chronological aging, although pathway-level differential methylation is not well defined. Using limma summary statistics from independent discovery and replication studies using the EPIC array platform described by McCartney et al (2020),²¹ we applied the LAM method to explore the pathway-level differential methylation. In the discovery cohort, we observed 304 differential pathways (FDR<0.05), while in the replication group, we observed 107. There was a high degree of agreement, with 43 pathways significant in both groups. This is not surprising, as the gene-level rank differential methylation scores show a strong positive association between discovery and replication studies (Figure 6A). Using the rank-MANOVA test of the mitch package, we identified 390 pathways with altered methylation. Visualised as a scatterplot, the overall pattern is concordant, although a cluster of points in the lower right of the chart indicates 268 pathways with higher methylation in the discovery sample and lower methylation in the replication sample (Figure 6B). When prioritising results by enrichment score, 29 of the top 30 pathways were concordant in their direction of regulation (Figure 6C). Some of the observed methylation changes make sense with reports in the literature. For example, the complement pathway is associated with aging²⁸ and in this analysis we observe a reduction in methylation of genes belonging to this set including C1R, MBL2, C1S and C1QC (Figure 6D).

Figure 6. Pathway-level DNA methylation alterations with chronological age. (A) Contour heatmap showing the similarity in gene methylation score ranks in the discovery and replication studies. (B) Mitch pathway enrichment scores in discovery and replication studies. Pathways with MANOVA FDR<0.05 are shown in red while others are shaded grey. (C) Heatmap of 30 pathways with largest absolute enrichment scores. Red indicates increasing methylation and blue shows lower methylation. (D) An example of a pathway identified with this method, “creation of C4 and C2 activators” shows lower methylation of member genes in both discovery and replication studies.

Exploration of methylation pathway differences in infants conceived with in vitro fertilisation

Another application for this method is the joint pathway enrichment analysis of studies conducted with different array systems. Epigenetic differences between infants conceived naturally and with assisted reproductive technologies have been the subject of studies since 2009.²⁹ Two of the highest-powered studies with publicly available data are Estill 2016 which used the HM450K array and Novakovic 2019 study which used the EPIC array.^22,23 Joint enrichment analysis can uncover pathways that are common or different between these studies and illuminate the biological differences between conception groups. For this example, we focused on comparing naturally conceived infants to those conceived with in vitro fertilisation (IVF) with fresh embryo transfer. The Estill dataset on the HM450K platform had 418,833 probes that met the quality criteria (Natural n=43; IVF n=38). Of these, 149,589 probes showed statistically significant differential methylation (FDR<0.05), of which 28,767 and 120,822 exhibited higher or lower methylation in the IVF group respectively. The Novakovic dataset on EPIC array had 793,844 probes that met the filtering criteria (n=133). Of these, 5,562 were statistically significant (FDR<0.05) with 1,348 and 4,214 exhibiting higher or lower methylation in the IVF group respectively. After summarisation to the gene level, the Estill dataset on HM450K array describes 19,240 genes, while the Novakovic dataset on EPIC array describes 22,588 genes. There were 19,234 genes common to both platforms. The rank-rank plot of differential gene methylation scores indicates a high degree of similarity overall, and interestingly shows a trend of more genes having lower methylation levels in the IVF group (Figure 7A). A scatterplot of pathway enrichment scores indicates a moderate degree of agreement between these studies (r=0.35, p=2.2e-16) despite these studies being conducted on independent cohorts years apart and analysed with different array systems (Figure 7B). A heatmap shows a high degree of agreement of pathways with higher methylation in the IVF groups, while pathways with overall lower methylation in the IVF group showed some variability in enrichment scores between studies (Figure 7C). We were interested in which pathways had consistently lower methylation in the IVF group and identified the Reactome pathway “Adrenoreceptors” which had enrichment scores of -0.42 in the Estill study and -0.56 in the Novakovic study, recording a MANOVA FDR value of 0.04 (Figure 7D). Top-ranked genes from this pathway include ADRA2C, ADRB2, ADRA2B and ADRA1B. A loss of methylation of these genes could contribute to reported elevated blood pressure in people conceived by IVF.³⁰ Taken together, this analysis shows that LAM method can be used to compare pathway enrichment across independent studies conducted with HM450K and EPIC arrays.

Figure 7. Pathway-level DNA methylation differences in natural and IVF conceived infants. (A) Contour heatmap showing the similarity in gene methylation score ranks in the Estill (HM450K) and Novakovic (EPIC) studies. (B) Mitch pathway enrichment scores in Estil and Novakovic studies. Pathways with MANOVA FDR<0.05 are shown in red while others are shaded grey. (C) Heatmap of 30 pathways with largest absolute enrichment score. Red indicates higher methylation and blue shows lower methylation. (D) An example of a pathway identified with this joint enrichment analysis method, “Adrenoreceptors” shows lower methylation of member genes in both Estill and Novakovic studies.

Differential pathway methylation in 19 common disease states

To examine the potential for LAM to reveal differential pathway methylation associations with diseases, we obtained and analysed blood-based EWAS summary statistics from a recent publication.²⁶ This study examined methylation in 14 prevalent disease states and the incidence of 19 disease states in a group of 18,413 participants, using the EPIC array. Delta-beta values were used for the purpose of scoring probe differential methylation, and these values underwent LAM analysis with Reactome gene sets.

There were 899 differentially methylated pathways in the prevalence arm of the study (FDR<0.05) (Figure 8A). On average, there were 64 pathways with differential methylation in each prevalent condition, with Chronic obstructive pulmonary disease (COPD) having the most (263) and Alzheimer’s disease having the fewest (3). We selected up to five statistically significant pathways with an absolute enrichment score of >0.4 in either direction in each prevalent condition to include in a heatmap (Figure 8B). We observed reduced methylation in the angiotensinogen metabolism pathway in patients with diabetes. Patients with colorectal cancer had lower methylation of the interleukin 10 signaling pathway. Reduced methylation of the salty taste perception pathway was observed in patients with CKD.

Figure 8. Differential pathway methylation associated with prevalence of 14 common disease states. (A) A bar plot showing the number of statistically significant pathways with higher and lower methylation identified in each prevalent condition (FDR<0.05). (B) A heatmap of enrichment scores for selected pathways across 14 common prevalent disease states. Stars indicate that the pathway was identified as being among the top five differentially methylated pathways in each direction for each condition.

We also wanted to know whether the LAM method could identify signatures that appear before patients are diagnosed with a condition. Therefore, we examined differential pathway methylation associated with the incidence of 19 disease states. Across all 19 conditions, there were 1570 significant pathways (FDR<0.05) (Figure 9A). Liver cirrhosis had the most, with 464, followed by COPD with 400, while Alzheimer’s had just one and COVID-19 hospitalisation had none. A heatmap shows the enrichment scores for selected pathways across these incident conditions (Figure 9B). Breast cancer incidence was associated with reduced methylation to the metal sequestration pathway. Incident prostate cancer was associated with reduced methylation of the FGFR1 pathway. Incident Parkinson’s disease was associated with reduced methylation of aquaporins. In COPD incidence, higher methylation of carbohydrate-binding ficolins was observed together with reduced methylation of the GDP mannose synthesis pathway. Incident stroke was associated with higher methylation of the lipoxin synthesis pathway.

Figure 9. Differential pathway methylation associated with incidence of 19 common disease states. (A) A bar plot showing the number of statistically significant pathways with higher and lower methylation identified in each incident condition (FDR<0.05). (B) A heatmap of enrichment scores for selected pathways across 19 common incident disease states. Stars indicate that the pathway was identified as being among the top three differentially methylated pathways in each direction for each condition.

To show that these findings are not the result of false positives from the LAM method, we randomised the incidence profiles for all 19 conditions prior to enrichment analysis. This was repeated 100 times, and these analyses of randomised data yielded few significant pathways. Specifically, from 100 runs 64 yielded no significant pathways across 19 conditions. Twenty-eight repeats had fewer than five false positives and eight runs had more than five false positives. Across the 100 repeats, the mean number of false positives was 1.95 and the median was 0. These results indicate that LAM can identify differentially methylated pathways in disease groups and detect altered methylation before a diagnosis has been made, with few false positives.

Data availability

Data sets reanalysed here are publicly available from NCBI GEO, Zenodo or from the supplementary tables of journal articles as described in the Methods section.

Code availability

Code to generate the results shown here are included in a GitHub repository (https://github.com/markziemann/gmea/). To enable reproduction of the findings, a docker image is made available at DockerHub (https://hub.docker.com/repository/docker/mziemann/gmea). A reproducible example workflow report is hosted on our website (https://ziemann-lab.net/public/gmea/example_workflow.html).

Funding

This work was supported by a grant from the Australian National Health and Medical Research Council (grant number 1146333 to JMC). Severine Lamon is supported by an Australian Research Council (ARC) Future Fellowship (FT210100278). This research was supported by use of the Nectar Research Cloud, a collaborative Australian research platform supported by the NCRIS-funded Australian Research Data Commons (ARDC). Authors gratefully acknowledge the contribution to this work of the Victorian Operational Infrastructure Support Program received by the Burnet Institute.

Acknowledgements

We thank the researchers who publicly shared data to enable this work.