Machine learning model to incorporate high-dimensional DNA methylation CpGs

MOSES aims to identify trait- or disease-associated genes through DNA methylation. Following a structure similar to TWAS, MOSES encompasses three key steps: training, imputation, and association (Fig. 2) (see Methods). In the initial step, MOSES employs elastic-net penalized regression models on a training dataset to estimate weights for CpGs within 1 Mb of promoter regions, both SNPs and CpGs within 1 Mb of promoter regions, or CpGs within 10 Mb from promoter regions for the expression of a gene (g) (Fig. 2). The elastic-net has been used to identify associated variables and predict outcomes in high-dimensional correlated data (see supplementary methods). The choice of 10 Mb is rooted in our breast cancer research, where we observed that prediction accuracy plateaus after 10 Mb, while computational time increases significantly beyond this threshold [15]. During this step, MOSES strategically eliminates many CpGs/SNPs not associated with the gene from the model due to the L1 penalty in elastic-net regression. This results in greater weights for more impactful CpGs/SNPs that either regulate or are associated with the expression of each gene. Moving to step 2, MOSES utilizes the coefficients (weights) estimated in the first step and the selected CpGs/SNPs to predict eGReX/GReX (when SNPs are involved) in a test DNA methylation dataset. Finally, in step 3, utilizing the imputed gene expression from the previous step, MOSES performs gene-level association tests for a given trait or condition.

Fig. 2

MOSES framework. MOSES (1) estimates weights for DNA methylation markers (and SNPs for MOSES-SNP + DNAm 1 M) from the elastic-net penalized regression models using the training data set, (2) predicts gene expression using a test data set, and (3) conducts the gene-level association tests using the imputed gene expression levels to identify trait-associated genes

To facilitate accessibility, we have made the MOSES R package publicly available at https://github.com/YidiQin/MOSES, and it will soon be available on Bioconductor as well The package provides two modes for users. In the first mode, users can train the model and conduct association tests using any tissue. This involves inputting training (reference) data comprising DNA methylation, gene expression data, and optional SNP data for the same subjects. The model then estimates tissue-specific effects of DNA methylation on expression for each gene. Subsequently, users input testing data for DNA methylation, optional SNPs, and a phenotype to impute gene expression and conduct gene-level association tests. In the second mode, users can streamline the process by skipping step 1 and directly perform gene-level association tests using a test dataset (DNA methylation) of nasal (airway) epithelium. By inputting DNA methylation and a phenotype, MOSES leverages precalculated coefficients from our reference data (EVA-PR) uploaded to GitHub to impute gene expression and conduct gene-level association tests, eliminating the need for additional reference data.

MOSES demonstrates superior performance in predicting gene expression in nasal epithelium

Utilizing genome-wide SNPs and CpGs along with RNA-seq data from nasal epithelial samples from EVA-PR participants, we conducted a comprehensive comparison of prediction accuracy for gene expression across different methods: cis-methylation (MOSES-DNAm 1 M), methylation within 5 Mb from promoter region (MOSES-DNAm 5 M), methylation within 10 Mb (MOSES-DNAm 10 M), both cis-DNAm and cis-SNPs (MOSES-SNP + DNAm 1 M) (see Methods), PrediXcan [8], which uses cis-SNPs (see Supplementary Materials), Biomethyl [14] and MethylXcan [13], both of which use DNA methylation in gene regions (see Methods). We measured CV-r between predicted and observed gene expression for 14,260 genes.

To compare all the methods, we defined 7885 “valid genes” in which at least one CpG site is in gene region, and which are predictable with r > 0.1 at least one of the four methouds-MOSES-DNAm10M, PrediXcan, MethylXcan, and Biomethyl as various TWAS methods suggested the threshold r > 0.1 (corresponding to R2 > 0.01)[9, 12, 23]. MOSES-DNAm 10 M (mean CV-r = 0.26) significantly outperforms PrediXcan, MethylXcan, and Biomethyl (mean CV-r = 0.03, 0.06, and 0.14, respectively) (Fig. 3A). We also simulated gene expression using DNA methylation of EVA-PR in various simulation scenarios (see Method in Supplemental Material). Consistent with our real data analysis, MOSES-DNAm 10 M also outperforms MethylXcan and Biomethyl in various simulation scenarios (Fig. S1). Therefore, for the rest of the analysis, to predict gene expression better and identify more genes, we excluded MethylXcan and Biomethyl so that we can include genes that do not include CpG site in gene region but include the CpG site within 1 M from the promoter region.

Fig. 3

Comparison of prediction accuracy in EVA-PR data. A Pearson’s correlation (r) between observed and predicted gene expression using tenfold cross-validation (CV) of EVA-PR data using MOSES-DNAm 10 M (CpGs \(\pm 10\) Mb from the promoter region), PrediXcan (SNPs \(\pm 1\) Mb from the promoter region), MethylXcan, and Biomethyl. We included 7885 genes for which at least one method predicted gene expression with CV-r > 0.1. Mean and a standard deviation bar. B Comparison of prediction accuracy between PrediXcan and MOSES-DNAm 10 M. C Comparison of prediction accuracy between MOSES-DNAm 1 M and MOSES-SNP + DNAm 1 M. D Comparison of prediction accuracy between PrediXcan and MOSES-SNP + DNAm 1 M. E Comparison of input variables (the number (#) of SNPs and/or # of DNA methylation CpGs) for the 4 methods F Comparison of the number of selected variables by the method (# of SNPs and/or # of DNA methylation CpGs) for the 4 methods. G GO Pathway analysis using 561 genes significantly better predicted by MOSES-DNAm 10 M (CV-r > 0.5 & FDR < 0.05 by Fisher’s z test). H GO Pathway analysis using 239 genes significantly better predicted by PrediXcan (CV-r > 0.5 & FDR < 0.05 by Fisher’s z test)

We predicted 14,260 genes using four versions of MOSES and PrediXcan, and of them 12,745 genes were “valid” genes using at least one method. When we compare three other versions of MOSES, MOSES-DNAm 10 M performs the best while it works very similarly to MOSES-DNAm 5 M in terms of Pearson’s correlation coefficients (Fig. S1A) and residual mean squared error (RMSE) (Fig. S1B). MOSES-DNAm 10 M predicted approximately five times more expression levels with higher accuracy than PrediXcan (10,578 vs. 2,167) (Fig. 3B). Interestingly, adding SNPs to cis-methylation data (MOSES-SNP + DNAm 1 M) does not significantly improve prediction accuracy compared to using methylation data alone (MOSES-DNAm 1 M) (Fig. 3C). In fact, MOSES-DNAm 1 M better predicts the expression of more genes than MOSES-SNP + DNAm 1 M (7,938 vs. 6,333), emphasizing the limited additional effects of adding selected SNPs to CpGs in imputing gene expression. On the other hand, adding CpGs to SNPs (MOSES-SNP + DNAm 1 M) does significantly improve prediction accuracy compared with using SNPs only (PrediXcan) (10,996 vs. 1749) (Fig. 3D). MOSES-DNAm 10 M demonstrates superior performance over MOSES-DNAm 1 M (Fig. S2), underscoring the greater impact of long-range CpGs on gene expression compared to cis-CpGs in nasal epithelial tissue.

We additionally assessed the number of input SNPs or CpGs and associated variables selected by elastic-net in the final model. Despite MOSES-DNAm 10 M incorporating fewer than half the input variables required by PrediXcan on average (p = 3,622 vs. p = 7,721) (Fig. 3E), the number of selected variables using MOSES-DNAm 10 M is approximately three times higher (s = 45.6), on average, than that using PrediXcan (s = 16.4) (Fig. 3F). This suggests that, even though there are more SNPs than CpGs, fewer SNPs exert an influence on gene expression compared to CpGs. While the addition of cis-SNPs to cis-CpGs does not significantly impact prediction accuracy, it is nearly additive concerning the average number of selected variables: using cis-DNA methylation (MOSES-DNAm 1 M) = 26.8, cis-SNPs (PrediXcan) = 16.4, and using both (MOSES-SNP + DNAm 1 M) = 40.0. Thus, the inclusion of selected (associated) SNPs alongside CpGs has minimal additional effects on imputing gene expression.

We also compared the running times of all the methods, and as expected, the running time scaled proportionally with the number of input variables (Table S1). Notably, MOSES-DNAm 10 M required significantly less time to predict gene expression than PrediXcan (4.49 mins vs 23.10 mins) and MOSES-SNP + DNAm 1 M (45.45 mins), while it required a little more time than MOSES-DNAm 1 M (4.49 mins vs. 4.10 mins).

For genes with expression predicted with high accuracy (r > 0.5) using either PrediXcan or MOSES-DNAm 10 M, we compared the prediction accuracy of PrediXcan with that of MOSES-DNAm 10 M to test whether there are statistically significant differences for each gene (see Supplemental Methods). The expression of 239 genes is significantly better predicted by PrediXcan (Table S2A), while that of 561 genes is significantly better predicted by MOSES-DNAm 10 M (Table S2B). Intriguingly, genes with superior prediction by PrediXcan are enriched in metabolic, catabolic, and biosynthetic pathways, whereas those with better prediction by MOSES-DNAm 10 M are enriched in immune pathways, including cytokine, antigen receptor, T cell, and interferon-gamma (Fig. 3G, H). The top 30 genes with significantly better prediction by MOSES-DNAm 10 M (mostly r > = 0.6) exhibit poor prediction by PrediXcan, demonstrating all negative correlations (Table 1). Notably, the top three genes—NTRK2, NLRC3, and CST1. all relevant to asthma—showcase the superior predictive capacity of MOSES-DNAm 10 M. NTRK2 (TrkB) expression is increased in the airways of subjects with asthma and may play a role in airway hyper-responsiveness [24], NLRC3 influences T cell regulation [25], and CST1 (the most significant DEG for atopic asthma in EVA-PR [20]) promotes eosinophilic airway inflammation in asthma [26]. In contrast, the top 3 genes significantly better predicted by PrediXcan—SNHG5, LINC01207, and STEAP1 (Table S1A)—are not found to be relevant to asthma in published studies. However, these genes are relevant to other lung diseases such as COPD (SNHG5) [27] and lung cancer (LINC01207 and STEAP1) [28, 29].

Table 1 The 30 most significantly improved predictions by MOSES-DNAm 10 M compared to PrediXcanMOSES-DNAm 10 M especially exhibits enhanced predictive capabilities for immunity-related genes and differentially expressed genes (DEGs) in atopic asthma

To identify which types of genes are better predicted by MOSES-DNAm 10 M, we investigated immunity-related genes, as they are relevant to asthma. First, we compared the predictive accuracy of the four methods for the expression of 753 immunity-related genes sourced from IMMPORT (Bioinformatics For the Future of Immunology)[19]. MOSES-DNAm 10 M predicts the expression of these genes significantly better than that of the same number of randomly selected genes (P < 2.2 × 10–16) (Fig. 4A). MOSES-DNAm 1 M and MOSES-SNP + DNAm 1 M also demonstrated superior prediction accuracy for these immunity-related genes over randomly selected genes (P < 2.2 × 10–16 and P = 9.6 × 10–11, respectively), while PrediXcan did not exhibit significant improvement (P = 0.99). Importantly, the prediction accuracy of immunity-related genes (r) for MOSES-DNAm 10 M (0.321) surpassed those of PrediXcan (0.012), MOSES-DNAm 1 M (0.278), and MOSES-SNP + DNAm 1 M (0.271) on average, indicating the pronounced collective impact of CpGs on the expression of immunity-related genes compared to SNPs.

Fig. 4

Prediction accuracy (CV-Pearson’s correlation, r) of expression levels for various categories of genes using four different methods. A Accuracy of predicting expression of 753 immunity-related genes versus that of the same number of randomly selected genes using the four methods. The list of immunity-related genes was obtained from the IMMPORT database [19]. B Accuracy of predicting expression of immunity-related genes in 8 pathways using MOSES-DNAm 10 M. The pathways/gene lists were obtained from the IMMPORT database [19]. C Accuracy of predicting differentially expressed genes (DEGs) in atopic asthma vs that of the same number of randomly selected genes using the four methods. The list of DEGs was obtained from Forno et al. [20]

In addition to individual gene analysis, we explored the prediction accuracy of MOSES-DNAm 10 M for eight immune gene pathways available in IMMPORT (Fig. 4B). The chemokine pathway, crucial for regulating airway inflammation in asthma [30], exhibited the highest accuracy in gene expression prediction. The antigen presentation pathway relevant to dendritic cells [31] influenced by innate immunity [32] and natural killer cells modulating immune responses in the airways of subjects with asthma [33] were the second and third most accurately predicted pathways, respectively.

The immunity-related genes whose expression was best predicted were those in the human leukocyte antigen (HLA) family, such as HLA-DQB1 (r = 0.88), HLA-DQA1 (r = 0.77), and HLA-DRB5 (r = 0.66), known to be associated with asthma and total IgE (Table 2). Other genes with well-predicted expression included the CXCR6/CXCL16 signaling axis (r = 0.69), which controls the localization of resident memory T (TRM) cells to different compartments of the lung and maintains airway TRM cells [34], and CD 247 (r = 0.66) which is part of the T cell antigen receptor (TCR) complex and has been associated with moderate to severe asthma [35].

Table 2 Prediction accuracies (CV-R) of immune genes using four different methods and the pathways the genes belong to

Given the inclusion of subjects with and without asthma in EVA-PR, we compared the performance of the four methods in imputing expression of DEGs in atopic asthma, previously identified in a differential expression analysis (DEA) within the same cohort [20]. Among 6,058 DEGs, 4,796 were considered valid genes (r > 0.1 for at least one of the four methods). MOSES-DNAm 10 M significantly outperformed PrediXcan, predicting the expression of these DEGs more accurately than that of randomly selected genes (P < 2.2 × 10–16 and 1, respectively) (Fig. 4C). This underscores the substantial combined effects of cis- and trans-CpGs on the expression of DEGs in atopic asthma, surpassing the impact of SNPs.

MOSES-DNAm 10 M exhibits superior performance over Biomethyl and MethlyXcan in conducting gene-level association tests for atopic asthma in nasal epithelium

We next conducted gene-level association tests of atopic asthma using imputed gene expression in nasal epithelium, performed through fivefold nested cross-validation using EVA-PR data (see Methods). While Biomethyl and MethylXcan are not developed to conduct gene-level association tests but to conduct pathway enrichment analysis and predict gene expression, respectively, we conducted DEA using limma for a comparison using predicted expression using each method. In this experiment, MOSES-DNAm 10 M again outperforms MethylXcan, Biomethyl, and other MOSES methods (Table S3). Especially, MOSES-DNAm 10 M successfully identified 4030 previously reported DEGs using observed gene expressions (84.9% of the 4,749 DEGs identified in the prior DEA [20] at FDR-P < 0.01), while MethylXcan identified 220 DEGs (4.6% of all the DEGs) and Biomethyl identified 237 DEGs (5.0% of all the DEGs).Table 3 showcases the top 30 DEA genes based on P values using observed gene expression, emphasizing the robust performance of MOSES-DNAm 10 M. The average cross-validated correlation coefficient (CV-r) for these 30 genes using MOSES-DNAm 10 M was 0.618, while MethylXcan and Biomethyl yielded an average CV-r of 0.111 and 0.274, respectively.

Table 3 Comparison of gene association test results across MOSES-DNAm 10 M, MethylXcan, and Biomethyl utilizing EVA-PR data

The highlighted genes encompass crucial candidate genes for asthma, such as CST1[20] and ALOX15[6]. Multiple ALOX15 SNPs have established associations with airway diseases, and the gene's DNA methylation and histone modifications are linked to airway inflammation [36,37,38,39]. The notable success of MOSES-DNAm 10 M in identifying these genes underscores its effectiveness, providing valuable insights for further research and clinical implications of long-range DNA methylation.

Gene-level association tests of atopic asthma in nasal epithelial samples from an external cohort

After building the MOSES-DNAm 10 M models trained using data from 455 participants in EVA-PRs, we tested this model using public data from 69 children in an independent cohort (GSE65205, including 36 children with atopic asthma and 33 non-atopic controls)[22]). Assessing the prediction accuracy of gene expression, the average correlation coefficient (r) was 0.229 for 10,477 genes (Figure S3), demonstrating a comparable range of prediction accuracy to that observed in EVA-PR (average r = 0.256 for these genes).

Next, we conducted a gene-level association test for atopic asthma using imputed gene expression. Leveraging observed gene expression, we identified 42 DEGs among valid genes (FDR-P < 0.05) (Fig. 5A). By employing imputed expression, MOSES-DNAm 10 M successfully identified 38 (90%) of these 42 DEGs (FDR-P < 0.05), all exhibiting concordant direction of associations (Table 4). The average test correlation coefficient (test-r) between observed and imputed gene expressions for these 38 DEGs was 0.614.

Fig. 5

MOSES-DNAm 10 M performance of the gene-level association tests of an external cohort, GSE65205. The models were trained using EVA-PR (N = 455); and tested on GSE65205 data (N = 69). a Manhattan plot of DEA results using expression levels imputed using MOSES. The red line represents the significance threshold after Bonferroni correction (P < 0.05). Green dots indicate whether the genes are among the top 200 differentially expressed genes for atopic asthma in a DEA conducted in EVA-PR [20]. b Scatter plot of CST1 expression: observed gene expression of GSE65205 vs. gene expression predicted by MOSES-DNAm 10 M. Red dots indicate atopic asthma subjects, and blue dots indicate non-atopic controls. r is Pearson’s correlation between the predicted and observed gene expression levels in the test data set. c Box plot of CST1 expression in participants with atopic asthma vs. that in non-atopic controls using observed microarray gene expression (left) and predicted RNA-seq gene expression (right). The P values are from DEA results, empirical Bayes moderated t-statistics using package limma [17]. D. Scatter plot of ALOX15 expression: observed gene expression vs. gene expression predicted by MOSES. E. box plot of ALOX15 expression in atopic asthma vs. that in non-atopic controls using observed (left) and predicted (right) gene expression levels in GSE65205

Table 4 Differentially expressed genes (DEGs) in atopic asthma as identified by MOSES-DNAm 10 M

The most significant DEG identified in the gene-level association test of atopic asthma among the 38 genes was CCK (cholecystokinin), linked to airflow limitation in obesity-related asthma [40]. Other notable genes included CTSC, which influences cell-mediated immune responses and immune cell trafficking in the airways [41], as well as other genes discussed earlier. Examples of these genes are visually presented in Fig. 5B–E. CST1 is highly predicted by MOSES-DNAm 10 M (r = 0.89), suggesting that this gene's expression is primarily influenced by DNA methylation (Fig. 5B). Utilizing imputed gene expression, MOSES-DNAm 10 M accurately identified CST1 as upregulated in atopic asthma (Fig. 5C). A similar pattern was also observed for imputed expression of ALOX15 (Fig. 5D–E). These findings underscore the robust performance of MOSES-DNAm 10 M in uncovering relevant genes associated with atopic asthma in an independent cohort.

Application of MOSES-DNAm 10 M to data from white blood cells (WBCs) in EVA-PR

We employed MOSES-DNAm 10 M to impute gene expression in white blood cells (WBCs) using data from 464 samples in EVA-PR, featuring both DNA methylation and RNA-seq data, and PrediXcan using 466 samples in the same cohort, featuring both SNPs and RNA-seq data (see Supplemental Methods). Similar to the results observed in nasal epithelium, MOSES-DNAm 10 M exhibited superior performance compared to PrediXcan (P < 2.2 × 10–16; with average Cross-Validation correlation coefficients (CV-r) of 0.212 and 0.076 for MOSES-DNAm 10 M and PrediXcan, respectively). This analysis focused on 8,298 overlapping valid genes shared between nasal tissue and WBCs (with r > 0.1 for at least one of the four methods). Notably, the prediction accuracy of gene expression in WBCs using MOSES-DNAm 10 M was comparable to that observed in nasal epithelium (average CV-R was 0.212 in WBCs and 0.246 in nasal tissue for the shared 8,298 valid genes) (Fig. 6A).

Fig. 6

Application of MOSES-DNAm 10 M to other data sets. A Comparison of prediction accuracy (cross-validated Pearson’s correlation) between MOSES-DNAm 10 M and PrediXcan in nasal epithelium (n = 455) and white blood cells (n = 469) of EVA-PR subjects for 8,298 overlapping valid genes between the two types of tissues (CV-r > 0.1). A two-sided t test was used to test the significance of differences. B Comparison of p values of the gene-level association tests between MOSES and PrediXcan in white blood cells. 9.267 genes were tested. The red color indicates significant differentially expressed genes (DEGs) for atopic asthma detected by the DEA analysis (FDR < 0.05). The labeled genes were identified from previous DEA using observed gene expression data [42]. C Comparison of CV-r using MOSES -DNAm 10 M in predicting gene expression in nasal epithelium of asthma and control (n = 455), that in WBC of asthma and control (n = 469), and that in lung tissue of idiopathic pulmonary fibrosis (IPF) (n = 382). D a Manhattan plot of DEA results using expression levels imputed by the MOSES in the IPF data set. The red line represents the significance threshold after Bonferroni correction (P < 0.05). Yellow dots indicate DEGs using observed gene expression

We also conducted gene-level association tests for atopic asthma in WBCs, testing 9,267 valid genes. MOSES identified 60 DEGs in atopic asthma (FDR-P < 0.05), whereas PrediXcan could not identify any gene in this context (Fig. 6B). Notably, the DEGs identified by MOSES included nine genes previously reported in the DEA conducted using observed expression [42]. Among these genes, SLC29A1 plays a crucial role in adenosine-regulated suppression of IgE-dependent mast cells [43], IL5RA and ALOX15 were associated with eosinophil count in a transcriptomic analysis of blood samples from subjects with asthma [44], and PIK3R6 and SMPD3 expression in whole blood was associated with total IgE in subjects with asthma [45]. Other DEGs detected by MOSES-DNAm 10 M have been previously implicated in asthma or related traits [45,46,47,48,49]. These results underscore MOSES's effectiveness in imputing gene expression in WBCs and identifying genes associated with atopic asthma.

Application of MOSES-DNAm 10 M to lung samples from persons with idiopathic pulmonary fibrosis (IPF)

To evaluate the performance of MOSES-DNAm 10 M in a disease other than asthma, we applied the model to lung samples from persons with idiopathic pulmonary fibrosis (IPF, n = 226) and healthy controls (n = 156) (see Methods). MOSES performed even better in this dataset than in those for nasal epithelium and WBCs from asthma datasets with regard to gene expression prediction (mean r = 0.52) (Fig. 6C). Using 16,365 predicted gene expression levels, MOSES identified 7,108 DEGs out of 7,110 identified using observed expression (Bonferroni-p value < 0.05), showing excellent performance in identifying DEGs and confirming that MOSES can be effectively applied to other diseases or tissues.