Weighted 2D-kernel density estimations provide a new probabilistic measure for epigenetic age

Abstract

Epigenetic aging signatures provide insights into human aging, but traditional clocks rely on linear regression of DNA methylation levels, assuming linear trajectories. This study explores a non-parametric approach using 2D-kernel density estimation to determine epigenetic age. Our weighted model achieves similar predictive accuracy as conventional clocks and provides a variation score reflecting the inherent variability of age-related epigenetic changes within samples. This score is significantly increased in various diseases and associated with mortality risk in the Lothian Birth Cohort 1921. Thus, weighted 2D-kernel density estimation facilitates accurate epigenetic age predictions and offers an additional variable for biological age estimation.

Background

Aging is reflected by gains and losses of DNA methylation (DNAm) at specific sites in our genome and this can be used to determine donor age. Epigenetic clocks have become a central component of aging research [1]. Since they were first described [2, 3], epigenetic age-predictors have become more and more sophisticated and tailored toward specific applications [4]. First generation clocks are trained to correlate as highly as possible with chronological age, e.g., for applications in forensics [5,6,7]. Even in these clocks, the deviation of chronological and epigenetic age (delta age) is often indicative for all-cause mortality and affected by various diseases [8]. However, this phenomenon has been proven to decrease when the training sample size is large enough or when a correction for the white blood cells counts is applied [9]. To better reflect the biological age, multifactorial second-generation clocks have been described [10, 11]. In addition to age, they also integrate epigenetic indicators of other clinical parameters, such as blood counts, glucose levels, or blood pressure. Furthermore, third-generation clocks have been trained on large cohort studies and implement additional clinical measurements to better quantify the individual pace of aging [12].

Despite the wide range of epigenetic clocks that have already been established, most of them rely on the same approach: signatures of age-associated CG dinucleotides (CpGs) are preselected in a training dataset and then used for multivariable regressions under the assumption that there are linear changes of methylation with age. Particularly in childhood, age-associated DNAm changes follow a logarithmic pattern [13], which can be effectively addressed through non-linear corrections [5]—but the predictors still follow the assumption that there is a continuous trajectory—either linear or logarithmic. Besides, among older individuals it has been proven that the linear age predictors proposed by Horvath and Hannum regularly misestimate biological age, leading to predictions that are younger than the actual age [14]. This might partly be attributed to a survivor bias, wherein a reduced biological age in the elderly has facilitated the longevity of these individuals [15, 16]. However, this phenomenon may be also related to the nonlinear saturation effects observed in some aging-related CpG sites, where they converge to a stable value, i.e., 0 or 100% DNAm [17]. Recently, alternative non-linear epigenetic clocks have been described that are based on deep learning, e.g., by neuronal network frameworks in multidimensional space [18,19,20]. Furthermore, machine learning algorithms can be used to derive non-parametric epigenetic clocks based on Gaussian process regression models [21]. While these approaches theoretically use the full range of DNAm information of the training dataset, they are relatively complex, as they usually need a huge number of samples and computational power to work. For instance, Aliferi et al. developed a support vector machine-based clock by selecting markers from 51 different studies and more than 4000 patients [22], and the machine learning-driven clock created by de Lima Camillo and colleagues was trained and validated using 142 datasets and needs the full 20,318 CpGs that are shared among the Illumina Infinium HumanMethylation27, HumanMethylation450, and EPIC arrays as input [19]. Furthermore, all of these epigenetic clocks only provide one specific age-prediction for a given sample as a single output.

In this study, we introduce a novel methodology to build epigenetic clocks based on 2D kernel density estimation (KDE). For each individual CpG of the aging signature, the KDE creates a matrix that relates age and DNA methylation (DNAm) levels to corresponding density values, which are then translated into probabilities. Integration of these probabilities can then be used to determine the most likely age-prediction. Our KDE approach is non-parametric and does not require linear or logarithmic assumptions. We demonstrate that, particularly with a weighted KDE model, we can facilitate robust epigenetic age-predictions. Furthermore, the probability distribution provides insight into how consistent the age-associated DNAm is for a specific age-prediction in different sites of the genome. This variation score is a new complementary measure that may be useful to estimate biological age.

Results

Two-dimensional kernel distribution provides probabilistic epigenetic age estimates

To develop epigenetic clocks based on KDE, we used DNAm datasets of 13 different studies on human peripheral blood that were separated into a training and a validation set. We arbitrarily preselected CpGs with the highest Pearson squared correlation with chronological age in the training set (R² > 0.7: 27 CpGs, or R² > 0.6: 491 CpGs). For each preselected CpG, the probability distribution of the age for a given DNAm level was determined by KDE and then combined into a joined probability estimation, where the age corresponding to the maximum probability was considered the predicted age of the sample. However, due to the heterogeneous distribution of donor ages, the predictions would be enriched at ages that are overrepresented in the training set (Fig. 1A–E).

Fig. 1

Construction of the probabilistic 2D kernel age predictor. A The association of DNA methylation (DNAm) with chronological age is exemplarily depicted for the CpG sites cg00329615 for all the samples in the training set. B 2D density kernel of the scatter plot depicted in A. The horizontal line exemplarily depicts the DNAm level of a given sample. C Density distribution as estimated for the probability of age-predictions for the sample depicted in B. D Joint probability for 27 age-associated CpGs for the same sample. E Performance of the methodology shown A–D for all the samples in the training set (color code indicates different studies of the training set). F Histogram of the ages in the training set. G Normalization of 2D kernels by age (divided by the corresponding value in the histogram). H Density distribution for the exemplary sample (line) in G. I Joint probability for 27 CpGs after normalization by age. J Performance of the methodology shown F–I for all the samples in the training set. K Alternatively, the training set was split in 5-year bins from 0 to 90 years (plus one bin of samples older than 90 years) and 15 samples from each bin were taken. The Scatter plot again depicts this selection for cg00329615. L 2D density kernel of K. M Density distribution for the sample in panel L. N Joint probability for 27 CpGs of the same sample. O Performance of the methodology shown K–N for all the samples in the training set

Therefore, we normalized the KDE maps by the frequency of samples at each specific age. We excluded samples older than 85 years (N = 39), because their low frequency then resulted in high density in the 2D kernel maps. Adjustment of KDE by the histogram of donor ages provided more accurate predictions, but there were still clusters, indicating that the normalization now rather skewed for age-ranges with the highest density or with low number of samples (Fig. 1F–J).

To further reduce the effects of age distribution in the training set we subsequently divided the samples into 18 bins of 5-year intervals, ranging from 0 to 90 years, plus one 19^th bin with the samples older than 90 years. From each bin, 15 samples were selected to construct the 2D density kernels (Additional file 1: Figure S1). This resulted in 2D density kernels with a more even distribution of densities and epigenetic age predictions that correlated with chronological age (Fig. 1K–O). Thus, KDE can be used for a probabilistic approach for epigenetic age predictions.

Improvement by weighted 2D kernel age-predictions

When we initially tested our probabilistic kernel models based on either 27 or 491 CpGs, Pearson squared correlation with chronological age in the training set was R² = 0.85 or 0.79, respectively, with a median absolute error (MAE) of 5 or 6 years. However, in the independent validation set, the correlations and precisions of age-predictions were much lower (27 CpG model: R² = 0.37 and MAE = 9 years, Fig. 2A; 491 CpG model: R² = 0.22 and MAE = 11 years, Additional file 1: Figure S2A), suggesting that there may be off-sets between different studies at individual CpGs, e.g., due to batch variation.

Fig. 2

A weighted approach improves 2D kernel age predictions. A 2D kernel age prediction model was generated for 27 CpGs (R² > 0.7 in the training set). The model was trained on a subset of samples from the training set with uniform age distribution. When we applied the model to the validation dataset, the predictions were not reliable. Pearson squared correlation R² and median absolute error (MAE) are indicated. B The same CpGs were used to generate multivariate models based on the entire training set. C Alternatively, for each of the age-associated CpGs, a linear regression model was established to facilitate single CpG predictions. When we averaged these predictions, the performance was much lower than for the multivariable model. D The 2D kernel age-prediction model was further optimized by optimized weights for individual CpGs, which were determined by genetic algorithm optimization. E Heatmap of Pearson squared correlation (R²) between chronological and predicted age for WKDE and a set of widely used clocks, split in all datasets used for training and validation. F Heatmap of MAE for WKDE and the same set of clocks as in E, split in all datasets used for training and validation

For comparison, we used the same set of age associated CpGs to generate multivariate regression models. This approach provided higher correlation with chronological age in the training set (27 CpG model: R² = 0.94, Fig. 2B; 491 CpG model: R² = 0.98, Additional file 1: Figure S2B). In fact, when creating multivariate models with the training set by varying the number of CpG sites, the highest correlation in the validation set was achieved with the top 125 CpGs (R² = 0.92, Additional file 1: Figure S3A), while the highest precision in the validation set would have been achieved with the top 81 CpGs (MAE = 2.87 years, Additional file 1: Figure S3B), indicating that larger signatures are not always beneficial. In the independent validation set, the multivariable models with 27 or 491 CpGs revealed high correlation with age, too (27 CpG model R² = 0.84, Fig. 2B; 491 CpG model R² = 0.83, Additional file 1: Figure S2B). In contrast, when we just used the average of individual CpG predictions, the age-estimates were much less precise (27 CpG model: training: R² = 0.86, validation: R² = 0.58, Fig. 2C; 491 CpG model: training: R² = 0.77, validation: R² = 0.30, Additional file 1: Figure S2C). Furthermore, the ages of very young samples were overestimated due to exponential age-associated epigenetic changes in childhood [13]. These findings supported the notion that using conventional regression approaches as multivariable weighted approach was clearly advantageous.

Consequently, we have also adjusted our probabilistic kernel models to better weight the impact of individual CpGs. To this end, we used genetic algorithm optimization to identify the best weights for each CpG according to the cumulative error curves. This weighted kernel density estimate model (WKDE) improved predictions; also in the validation set (Fig. 2D, Additional file 1: Figure S2D)—particularly with the 27 CpG model, the predictions (R² = 0.94 and MAE = 5 years in the training set; R² = 0.81 and MAE = 4 years in the validation set) were now in a similar range as observed for the multivariable model. However, the 491 CpG model showed a very poor correlation in the validation set, which might be attributed to over-fitting of the weights and the increased amount features with lower age association. Therefore, we did not consider the 491 CpG WKDE model for further analysis.

Benchmarking with other commonly used epigenetic clocks

Subsequently, we compared the performance of the 27 CpG WKDE model with other widely used epigenetic clocks in the field. To better test the applicability to independent datasets, we have done this analysis for each study of the training and validation sets. The WKDE predictor revealed better or similar correlations to many other clocks, while particularly the Zhang clock and Horvath Skin & Blood clock showed highest correlations (Fig. 2E). Furthermore, the WKDE model provided small median absolute errors for all datasets—apart from Yang clock and both Horvath clocks, other clocks revealed high MAEs of more than 15 years for at least one dataset (Fig. 2F). Thus, despite the small epigenetic signature of only 27 CpGs the WKDE model facilitated robust estimation of chronological age.

Performance of the 27 CpG clock in purified cell types

To gain better insight into how the cellular composition affects age-predictions with our 27 CpG WKDE clock, we used a dataset with 6 purified cell types and provided chronological ages (GSE110554). Particularly for purified CD4 T cells, monocytes, neutrophils, and NK cells, the delta age was overall very low, whereas B cells were overestimated and CD8 T cells underestimated in age (Additional file 1: Figure S4A). These results were further validated with another dataset where the cell types were isolated form the same donors, but donor age was not provided (GSE224807). Again monocytes, granulocytes, and NK cells revealed similar age-predictions as for whole blood, whereas B cells were over-estimated, and T cells were predicted younger (Additional file 1: Figure S4B). These results demonstrate that while the 27 CpG WKDE clock can generally also be applied to purified cell types, the predictions may be influenced particularly by the composition of lymphocytes.

Applicability of the WKDE method with targeted assays

To test if our new WKDE method was applicable to targeted epigenetic clocks, we used our previous published pyrosequencing data for nine age-associated CpGs in FHL2, IGSF11, CCDC102B, MEIS1-AS3, ELOVL2, COL1 A1, PDE4 C, ASPA, and ITGA2B [23]. After the 2-dimentional kernel calculation step and the genetic algorithm to assign the weight for each kernel, the 9 CpG model notably achieved a good correlation between chronological and predicted age and low prediction error in the validation set (R² = 0.82, MAE = 5.33 weeks, Additional file 1: Figure S5), indicating that the approach is also suitable for targeted DNAm analysis, e.g. by pyrosequencing.

Heterogeneity of age-associated DNAm within a given sample

Our probabilistic approach allows to calculate the probability of the age of a given sample for every year from 0 to 100 (Fig. 3A). To estimate the disparity among the selected methylation sites, we calculate a variation score. This measure reflects the heterogeneity in age-associated DNAm within our 27 CpG signature, rather than the probability that the estimates of chronological age are correct. However, for most of the samples in training set (99.22% of total samples; Fig. 3B) and validation set (99.07% of total samples; Fig. 3C), the chronological age falls into this range (predicted age ± variation score). The variation score shows small correlations with chronological age (training: R² = 0.36, validation: R² = 0.33), epigenetic age (training: R² = 0.37, validation: R² = 0.24), and delta age (training: R² = 0.07, validation: R² = 0.02) and was highest at an age-range from 25 to 75 years (Figs. 3D–F). When we compared male and female samples, we did not observe any significant sex bias with regard to delta age (P = 0.27, Fig. 3G), or variation scores (P = 0.89; Fig. 3H).

Fig. 3

WKDE allows to measure heterogeneity in age-associated DNAm. A Possible age (y-axis) vs. probability (x-axis) for eight random samples from the training set. The thickness of the plot represents how probable is that a sample belongs to that age. B Chronological vs. predicted age in the training set measured with WKDE. Vertical lines depict predicted age ± variation score. C Same as in panel B for validation set. D Chronological age vs. variation score. E Predicted age vs. variation score. F Delta age vs. variation score. Variation score reflects almost no correlation with delta age measured with WKDE. G Boxplot of delta age vs. sex. Delta age distributes evenly across males and females (Wilcoxon rank sum test P = 0.27). H Boxplot of variation score vs. sex. Variation score distributes evenly across males and females (Wilcoxon rank sum test P = 0.89)

The variation score is significantly increased in several diseases

Epigenetic age-predictions can be impacted by diseases, particularly by hematopoietic malignancies [24]. Hence, we anticipated that the heterogeneity in age-associated DNAm changes, reflected by the variation score, may be affected in such diseases. To test this hypothesis, we used six publicly available DNAm datasets of healthy and diseased individuals for acute myeloid leukemia (AML). In fact, we observed pronounced differences in the distribution of variation score with age in all datasets, and when we compared diseased with available healthy individuals the difference was significant (GSE124413 P = 2e − 16, GSE58477, P = 2.93e − 5; GSE152710, P = 0.0126; Fig. 4A). Furthermore, the variation score declined upon therapy with azacytidine (GSE152710, P = 2.15e − 7) [25]. We have then exemplarily tested datasets of myelofibrosis, Down syndrome, HIV, progeroid syndromes, Parkinson’s disease, and schizophrenia (Additional file 1: Table S1), and calculated the corresponding variation score with our 27 CpG WKDE model. When we compared healthy and diseased individuals in these datasets, the variation score was significantly higher in myelofibrosis (P = 0.0131, Fig. 4B), Down syndrome (P = 2.13e − 10, Fig. 4C), and HIV (P = 1.33e − 12, Fig. 4D), whereas no significant difference was found in progeroid syndrome (P = 0.6943, Fig. 4E), Parkinson’s disease (P = 0.265, Fig. 4F), and schizophrenia (P = 0.105, Fig. 4G). Thus, high variation scores may be indicative for underlying diseases.

Fig. 4

Variation score is significantly different in certain diseases. A The variation score was analyzed in six different datasets of acute myeloid leukemia (AML) and plotted in relation to donor age if provided or otherwise in relation to predicted age (with the 27 CpG WKDE model). Untreated AML samples (green), treated AML samples with azacitidine (purple), healthy controls from the same datasets if provided (red), and for comparison the training set of healthy controls are depicted. B–G The same analysis was performed for datasets of B myelofibrosis, C Down syndrome, D HIV, E progeroid syndromes, F Parkinson’s disease, and G schizophrenia. The indicated P values are related to the comparisons indicated by the corresponding color code

Association of sample-intrinsic heterogeneity in epigenetic aging with all-cause mortality

It has been reported that all-cause mortality increases with accelerated epigenetic age in various clocks [8, 26]. To estimate if this also applies for our WKDE model, we analyzed methylation from the first waves of the LBC1921 and LBC1936 with subsequent mortality risk. Cox proportional hazards regression models for delta age vs. survival, adjusting for age and sex, revealed no significant association between delta age and mortality in the LBC1921 (P = 0.978, Fig. 5A) and in the LBC1936 (P = 0.415, Additional file 1: Figure S6A).

Fig. 5

Variation score shows significant association with mortality. A The association of epigenetic age-predictions with all-cause mortality was analyzed in the Lothian Birth Cohort of 1921 (LBC1921). Kaplan-Meier survival curves are depicted for highest and lowest 10% delta ages. Cox regression model for all donors, adjusted for chronological age and sex, showed no significant effect of delta age in mortality risk (HR = 0.9998, 95% CI (0.988, 1.012), P = 0.978). B Individual Cox regressions for all the 27 CpGs in WKDE (adjusting for age and sex) revealed 4 significantly mortality-associated CpGs in the LBC1921. The cg11436113 was also significantly associated with mortality in the Lothian Birth Cohort of 1936 (LBC1936). C Kaplan-Meier survival curves for donors with highest and lowest 10% DNAm at cg11436113 in the LBC1921. Cox regression model, adjusted for chronological age and sex, shows that increase of 1% in the DNAm of cg11436113 is associated with a 1.94% decrease in mortality risk (95% CI (0.9665, 0.9949), P = 0.0080). D Kaplan-Meier survival curves for highest and lowest 10% variation score in the LBC1921. Cox regression model, adjusted for chronological age and sex, shows that increase of 1 unit in the variation score is associated with a 9.2% decrease in mortality risk (95% CI (0.8387, 0.9872), P = 0.0160)

Subsequently, we tested if any of the individual CpG sites in the 27 CpG or 491 CpG models were associated with all-cause mortality, after adjustment for age and sex. Within the 491 CpGs, we found 25 mortality associated sites that overlap between LBC1921 and LBC1936 (Additional file 1: Figure S6B and Tables S2 and S3), with almost all of them being hypomethylated with aging. Out of the 27 CpG sites, four were significantly associated with mortality in the LBC1921 and one in the LBC1936 (Fig. 5B). Particularly cg11436113, which is located in an intergenic region, was highly associated with mortality in both cohorts (LBC1921: P < 0.008, Fig. 5C; LBC1936: P < 0.00003, Additional file 1: Figure S6C).

When testing for associations between variation score and mortality in the LBC1921, we observed that an increase of 1 unit in the variation score is associated with a 9.2% decrease in mortality risk (95% CI (0.8387, 0.9872), P = 0.0160, Figure 5D), after adjusting for sex and chronological age. The results remained significant when we adjusted for sex, chronological age, and delta age (P = 0.0159), indicating that the variation score might be an independent measure for biological aging. However, a significant association of the variation score with all-cause mortality was not observed in the LBC1936 (HR = 0.9894, 95% CI (0.8870, 1.1037), P = 0.849, Additional file 1: Figure S6D), which might partly be attributed to a lower number of deaths in this cohort.

Subsequently, we analyzed the association between variation score and mortality in a further adjusted model controlling for age, sex, delta age, and 22 additional parameters (Additional file 1: Table S4). Even this adjusted model indicated that a one-unit increase in the variation score is significantly associated with a 13.65% reduction in mortality risk within the LBC1921 (HR = 0.8635, 95% CI (0.7754,0.9617), P = 0.00754), which is still not reflected in the LBC1936 (HR = 1.0052, 95% CI (0.8765,1.1529), P = 0.9407).

Discussion

In this study, we demonstrate that 2D kernel density estimation can be used for robust epigenetic age-predictions. Even with a very simplistic selection of 27 age-associated CpGs based on Pearson correlation with chronological age, we could demonstrate that WKDE models can provide similar precision as conventional epigenetic clocks. However, there are several ways how feature selections might be optimized in the future. Age-associated DNAm does not necessarily follow a strictly linear pattern [27]. This is also reflected by the 2D kernel density plots of the 27 CpGs, where almost all of them followed a logarithmic pattern. Thus, it might be advantageous to select the age-associated CpGs by Spearman correlation, albeit there is a very high overlap in age-associated CpGs selected with Pearson correlation [23]. Another approach might be to select CpGs based on the slope in single-CpG linear regressions to exclude CpGs that may reveal high correlation but only small absolute changes with age [28]. The selection of CpGs would be further improved by a larger collection of datasets, including newborns and centenarians. Our training set had a predominant emphasis on individuals of European ancestry and should better reflect ethnic diversity [29, 30]. To this end, we included datasets from USA (both European and Hispanic-Mexican ethnicity), the Netherlands, Denmark, Germany, the UK, and China in the training and/or validation sets. CpGs that show higher variation between cell types might be excluded from the predictors to reduce the impact of the cellular composition [23]. This might reduce the variation in age-predictions that we observed in B cells and T cells. Theoretically, automatic CpG selection approaches, such as ElasticNet [31], could be used to further optimize the selection of CpGs for WKDE—in analogy to the ElasticNet based conventional epigenetic clocks [5, 7]—but it would require an innovative automatic feature selection algorithm for matrices (taking also different possible weights into account) and this would require very high computational resources. Thus, in this study we introduce a new and very promising concept for constructing epigenetic clocks, rather than a fully optimized clock.

Nonetheless, WKDE method has some limitations. Predictions are only possible within the age-range that is covered by the training set and the estimated kernels—in our case, we can therefore not reliably predict ages above 100 years. The imbalance of different donor-ages in the training set was a major problem, which we could neither entirely corrected by normalization with the age-distribution, nor by a conditional density-based approach for the creation of the 2D density kernels [32]. We have therefore randomly reduced the training set for a uniform age-distribution—with the downside that not the entire training set was considered. In the future, alternative approaches based on resampling or machine learning might be further exploited to solve this limitation [33, 34]. Furthermore, our WKDE model was optimized with a genetic algorithm, which has shown several advantages over other non-linear optimization techniques, such as high efficiency, parallelization, and no need of derivative information [35]. However, other tuning parameters for the genetic algorithm or even non-linear optimization algorithms can be tested in subsequent studies.

Larger epigenetic signatures do not always increase the precision. When we only considered the top 27 CpGs in our training set to build the WKDE model, we observed a higher correlation with chronological age and lower median absolute error than using a less stringent pre-selection of 491 CpGs. This might be due to the more stringent selection of highly correlating CpGs and the greater number of CpGs increases the risk of overfitting of the genetic optimization algorithm to find optimal weighs [36]. In fact, the 491 WKDE model performed very poorly on the validation set, indicating that there was some overfitting.

Our model was specifically trained for human blood samples, while similar models might also be trained for other cell types or tissues in the future. There is a growing perception that even the heterogeneity of cells in blood can have major impact on epigenetic age-predictions with many conventional epigenetic aging clocks, and this may contribute to associations between clock-derived measures and age-related health outcomes [37]. For example, the association of epigenetic age acceleration and mortality found with Horvath and Hannum clocks seems to disappear when a correction for the white blood cell counts is performed [9]. Since we observed that purified B cells are predicted older whereas T cells are predicted younger, we anticipated that our signature will be impacted by cellular composition. Depending on the specific purpose of a study, it might be crucial to exclude CpGs with high variation between leukocyte subsets [23]. On the other hand, the 2D kernel density approach might also be applied for other epigenetic biomarkers, such as deconvolution of leukocytes in blood [38], or of cell types within a tissue [39, 40].

A big advantage of smaller signatures is that they can also be addressed by other methods for targeted DNAm analysis, such as pyrosequencing, digital PCR, or barcoded bisulfite amplicon sequencing [23]. In comparison to Illumina BeadChip technology, these approaches may facilitate faster and more cost-effective measurement. Furthermore, the measurements are independent of continuously changing microarray versions, which may anyway not be approved for clinical diagnostics [41]. In our previous work, we have measured 9 age-associated CpGs by pyrosequencing to derive an epigenetic aging clock with high accuracy based on pyrosequencing [23]. Notably, the same data can also be used to derive and validate a 9 CpG WKDE model for pyrosequencing data. This exemplifies that our approach is also applicable for targeted methods, particularly if the available training data is large enough.

Probabilistic approaches for epigenetic clocks have been described before. For example, we have previously described probabilistic epigenetic age predictions for individual sequencing reads, using the binary sequel of methylated and non-methylated CpGs in barcoded bisulfite amplicon sequencing data [23]. This method has been further developed for genome wide single-cell DNAm datasets [42]. In contrast, WKDE is not applicable to shallow sequencing data, but it provides density estimates for each CpG, a joint probability range for the sample and a variability score associated to it.

An important advantage of our approach is that the variation score provides a new measure for the intrinsic heterogeneity of age-associated DNAm changes. It is conceptionally somewhat related to a recently described “noise barometer” that is based on sums of standard deviations of DNAm at individual CpGs [43]. Notably, our variation score was clearly increased in malignancies, such as AML and myelofibrosis. Thus, even without knowing the chronological age, the combination of epigenetic age prediction and variation score can provide a “red flag” to initiate further diagnostic procedures. It was also striking that HIV, which has been shown to be associated with accelerated epigenetic clocks [44], might rather reflect increase heterogeneity of DNAm at age-associated CpGs.

Although association between accelerated age and mortality has been shown before [8, 26], we did not observe this with WKDE in the LBC1921, which might be due to the sample selection and the nature of the clock itself. We have intentionally only chosen samples from the initial waves of the LBC1921 and LBC1936 to isolate the effects of the variables under examination on mortality from the confounding factor of aging (which differs from the previously mentioned studies). It has to be noted, that our WKDE model was trained to correlate as closely as possible with chronological age and it is therefore conceivable that mortality-associated CpGs sites—such as cg11436113—were considered with a negative weight for age-predictions. Furthermore, WKDE employs a method of prediction that focuses on determining the age that maximizes the probability of a joint probability function for the age, rather than relying on a singular value. Hence, the biological age might still fall into the probability range of the age.

We have previously identified individual CpG sites from our 99-CpG clock and other epigenetic clocks that showed correlations with survival in both the LBC1921 and LBC1936 [45]—and most of these were hypomethylated with age [46]. Here, we showed an overlap of 25 methylation sites that are significantly associated with mortality in both cohorts, and again most of these CpGs become hypomethylated with age. Notably, these CpGs included cg06285727, cg06647068, cg08301612, and cg20813374, which were previously reported as all-cause mortality associated in a large-scale meta-EWAS study of six large cohorts [47]. Furthermore, 7 of them are located in genes BST2, FKBP5, PRDX5, NWD1, FAM38 A, and NOD2, which have previously shown high association between DNAm and aging [48,49,50,51], and we have already identified cg16363586 (located in gene BST2) as life expectancy associated site [45]. In addition, the CpG site that was identified as mortality-associated in both cohorts with the 27 CpG signature, cg11436113, has been reported as smoking and cancer related [52, 53].

While previous studies have discussed the link between age-related DNA methylation (DNAm) and mortality [54, 55], the connection between age-related DNAm heterogeneity and all-cause mortality was so far not addressed. Our results suggest that heterogeneity among age-associated methylation sites in the genome, measured as variability score, can provide additional insight into biological age. It may appear counterintuitive that here a lower variation score was associated with higher mortality, albeit very high variation scores can be indicative for other diseases, as mentioned above. This might be due to different pathophysiological mechanisms in aging and disease, or with the observation that particularly in the elderly the variation score declines again. Nevertheless, further validation is necessary to examine the potential application of the variability score to estimate life-expectancy and to explore possible association with specific diseases.

Conclusions

Our study introduces weighted 2D-kernel density estimation (WKDE) to perform accurate epigenetic age predictions. Furthermore, we describe a variation score that may act as an additional parameter for evaluating biological age.

Methods

DNA methylation data

We compiled DNAm datasets of human peripheral blood samples of 13 different studies from Gene Expression Omnibus (GEO; Additional file 1: Table S5). We considered only samples that were classified as healthy or controls. Datasets were separated into a training set (7 Illumina HumanMethylation450 K array studies, 1029 samples, age range 1–101 years, 50.1% female) and an independent validation set (6 Illumina HumanMethylation450 K array studies, 980 samples, age range 2–79 years, 53.5% female). The data was processed in R 4.3.0. with the geoGEO function of the GEOquery package. To simplify the analysis and make it easily applicable to other groups, the preprocessed and normalized beta value matrices for each study were downloaded directly from the Gene Expression Omnibus (no additional background correction, normalization or quality control steps were performed, as they were already performed by the corresponding studies). CpGs in sex chromosomes and 6749 single-nucleotide polymorphism (SNP) probes (according to the HumanMethylation450 K v1.2 annotation) were filtered out.

Epigenetic clock based on kernel density estimation

As a proof of concept that employing two-dimensional (2D) kernel estimations is feasible for accurate age predictions, we opted for one of the simplest feature selection approaches for age predictions: we have focused on the CpGs with highest linear correlation with chronological age in the training set, either 27 CpGs with R² > 0.7 (Additional file 1: Table S6), or 491 CpGs with R² > 0.6 (Additional file 1: Table S7). To generate KDE clocks based on these CpGs, we tested various alternative approaches. Initially, we used all 1029 samples of the training set to generate 2D kernel density plots of chronological age versus DNAm. Since this approach resulted in offsets of epigenetic age-predictions due to the heterogeneous age-distribution in the training set, we subsequently adjusted densities of the original kernel by the frequency of donor ages. To this end, we only considered samples ≤ 85 years to avoid a distortion for age-categories with very few samples. Alternatively, we tested conditional density resampled estimate of mutual information (DREMI) [32], which computes a two-dimensional kernel density estimate in cases where one variable, e.g., age, stochastically influences the other variable, e.g., DNAm. Even though this approach was meant to adjust the densities for the variations in DNA methylation linked to aging, it did not yield any enhancements in age prediction (data not shown). Finally, we randomly selected 15 samples of the training set for each 5-year bin from 0 to 90 years, and for one additional bin with samples older than 90 years. This resulted in a subset of 285 samples with uniform age distribution.

For each CpG, the kernels were generated with the function kde2d of the MASS R package, with 101 grid points, resulting in kernels with age on the x-axis (0–100 years), beta values on the y-axis (0 to 100% DNAm), and the density on the z-axis. Thus, each CpG kernel corresponds to a 3D matrix K, where ${K}_{i,j,k}$ represents the density value of the CpG $k$ for a sample with age $i$ and DNAm percentage $j$. An R object containing the 27 kernels, one for each of the 27 CpGs, can be found in Additional file 2 (.Rds). At this point, a preliminary probabilistic kernel age prediction for a sample can be done by:

$$predicted\ ag{e}_{preliminary}=\underset{i\in \left[\text{0,100}\right]}{\text{argmax}}\left(\sum_{k}{K}_{i,DNAm\left(k\right),k}\right)$$

where $k$ iterates through the selected CpGs and DNAm(k) is the DNAm percentage of the CpG k in the considered sample. For each CpG k, the row corresponding to the DNAm(k) is a vector of 101 densities, each of them corresponding to the years 0 to 100, respectively. After performing the summation of all these density vectors, the algorithm returns the age $i$, the one with the maximum cumulative density of the summation, as the final age.

To improve the performance of the predictor, a genetic algorithm was implemented to include a weight coefficient ${w}_{k}$ for each CpG $k$, with possible values in the range [− 10, 10]. For this, the function minimize from the R package EmiR was used, with a number of 150 initial conditions (parameter population_size = 150, meaning that 150 different vectors of 27 or 491 randomly assigned weights between − 10 and 10 were set as initial conditions at the start of the algorithm), 100 iterations (parameter iterations = 100, indicating that the algorithm will stop only after 100 iterations have been completed, and will return the best vector of 27 or 491 weights among the 100^th generation of descendants), a keep fraction of 0.4 (parameter keep_fraction = 0.4, so that 40% of the vectors survive for the next mating step), and a mutation rate of 0.1 (parameter mutation_rate = 0.1, implying a 10% probability of random modifications in the vectors of the offspring). The minimization function was set as:

$$\text{minimize}\left(\sum_{l}\left|ag{e}_{l}-\underset{i\in \left[\text{0,100}\right]}{\mathit{argmax}}\left(\sum_{k}{w}_{k}\cdot {K}_{i,DNA{m}_{l}\left(k\right),k}\right)\right|\right)$$

where $l$ iterates for all the samples in the training set, $ag{e}_{l}$ is the real age of the sample $l$ and $DNA{m}_{l}\left(k\right)$ is the DNAm percentage at site $k$ of the sample $l$. The genetic algorithm adjusts the weights of each CpG to minimize the sum of absolute differences between the predictions and the real ages, therefore finding the coefficients that minimize the total absolute error in the training set (Additional file 1: Tables S6 and S7). An R object containing the 27 weights for each CpG can be found in the Additional file 3 (.Rds). After finding the optimal weights, the final predictions can be calculated as:

$$predicted\ age=\underset{i\in \left[\text{0,100}\right]}{\mathit{argmax}}\left(\sum_{k}{w}_{k}\cdot {K}_{i,DNAm\left(k\right),k}\right)$$

In a similar way, a probability vector can be calculated. Every vector has one entry for each year (0 to 100), which quantifies how probable it is to observe the measured methylation value assuming that a person is of the respective age. Since each CpG is associated with a weight coefficient, we have to add the individual weighed vectors for all the CpGs, which ends up in a single merged vector. The vectors can be subsequently normalized to interpret the weighted sum as a probability, i.e., to ensure that the sum of probabilities is equal to 1. This results in the following formula for weighted 2D KDE (WKDE model):

$$probability\left(age\right)=\frac{\sum_{k=1}^{n}{w}_{k}\cdot {K}_{age,DNAm\left(k\right),k}}{\sum_{i=0}^{100}\left(\sum_{k=1}^{n}{w}_{k}\cdot {K}_{i,DNAm\left(k\right),k}\right)}$$

where $n$ is the number of the CpGs that the clock uses, ${w}_{k}$ is the weight of the CpG $k$, $DNAm(k)$ is the measured DNAm percentage of the CpG $k$ in the considered sample and ${K}_{age,DNA\left(k\right),k}$ is the value in the row $DNAm(k)$ and column $age$ of the density map of the CpG $k$.

Epigenetic age predictions based on linear regression

To benchmark epigenetic age-predictions of our WKDE model with conventional approaches using the same subsets of CpGs (R² > 0.6 or 0.7 in the training set, as indicated above), we established two alternative epigenetic clocks: (1) by multivariable linear regression using the lm function from stats package in R and the predict function to calculate age predictions for both training and validation sets (Additional file 1: Tables S6 and S7) and (2) by calculating the average of individual CpG age-predictions, whereby for each selected age-associated CpGs the slopes and intercepts between DNAm (beta values) and age were calculated with the lm and apply functions in R and the mean of the CpG-specific predictions was then used to estimate donor age (Additional file 1: Tables S6 and S7). Pearson squared correlation (R²) and median absolute error (MAE) between chronological and predicted age was calculated. Furthermore, we used the R package methylCIPHER [56] to compare the performance across the individual datasets used in training and validation sets with other commonly used epigenetic clocks: the Horvath clock [5], the Hannum clock [7], PhenoAge [10], retrained principal component PhenoAge [57], Horvath Skin & Blood epigenetic clock [58], Lin clock [45], Vidal-Bralo clock [59], Han clock [23], and Zhang clock [9].

Age predictions in purified cell types

To estimate epigenetic age in purified leukocyte subsets, we used normalized matrices of beta values from two datasets: GSE110554 (30 samples after removing smokers and samples with missing information on chronological age) [60] and GSE224807 (415 samples after removing smokers, age-information is not provided) [61]. The data was downloaded from the GEO repository and processed in R 4.3.0. with the getGEO function of the GEOquery package, and no additional preprocessing or quality control steps were performed.

Pyrosequencing-based WKDE clock

Pyrosequencing data of 80 healthy human whole blood samples (age range 19–73 years) was retrieved from a previous study [23]. The methylation values of 9 age-associated CpGs were used, and the samples were redistributed into training and validation sets so that the training set had a uniform age distribution (Additional file 1: Table S8). WKDE method was applied to build a clock with the 9 CpGs (both KDE calculation and genetic algorithm were performed as described before).

Variation score

The variation score provides an estimate for how homogeneous the age-predictions for individual CpGs are. Thus, we focus on the probability vectors themselves, without considering the weights. For this, we used the R function approxfun with the non-weighted probability function to obtain probability values for all ages from 0 to 100 in steps of 0.1 years. The function approxfun performs a linear interpolation between given data points:

$$f\left(x\right)=approxfun\left(\frac{{\sum }_{k=1}^{n}{K}_{age,DNAm\left(k\right),k}}{{\sum }_{i=0}^{100}\left({\sum }_{k=1}^{n}{K}_{i,DNAm\left(k\right),k}\right)}\right)\left(x\right),where\ x\in \left\{\text{0,1},\dots ,\text{999,1000}\right\}$$

The call returns a vector with 1000 components, where the component x (which is between 0 and 1000) corresponds to the probability that the sample comes from an individual of age x/10. Subsequently, we were able to calculate mean, variance, and standard deviation of the probability function:

$$mean= 0.1\cdot \sum_{x=0}^{1000}f\left(x\right)\cdot x$$

$$variance=0.1\cdot \sum_{x=0}^{1000}f\left(x\right)\cdot {\left(0.1 x-mean\right)}^{2}$$

$$variation\ score=\sqrt{variance}$$

The R script with the functions to calculate predicted age, probability age function, and variation score with WKDE can be found in the additional file 4 (.R).

Finally, the Mann–Whitney-Wilcoxon test was used in R with the function wilcox.test to verify if there were differences in delta age (predicted − chronological age) by sex.

Association between variation score and several diseases

The already preprocessed and normalized matrixes of beta values for the different studies (Additional file 1: Table S1) were directly downloaded from the GEO repositories. The data was processed in R 4.3.0. with the getGEO function of the GEOquery package, and no additional preprocessing or quality control steps were performed. When available, plots and statistics were performed with the chronological age of the samples; otherwise, the predicted age calculated with the 27 CpG WKDE model was used. The following repositories were used: GSE124413 [62], GSE133986 [63, 64], GSE153347 [65, 66], GSE62298 [67, 68], GSE58477 [69, 70], and GSE152710 [25, 71] for AML; GSE118241 [72, 73] for myelofibrosis, GSE52588 [74, 75] for Down syndrome, GSE67751 [76, 77] for HIV, GSE131752 [78, 79] for progeroid syndromes, GSE165081 [80, 81] for Parkinson’s disease, and GSE41169 [82, 83] for schizophrenia.

To assess whether there were significant differences in variation scores between healthy and diseased samples, ANOVA was utilized for each study individually. The variation score was designated as the response variable, with age (chronological or predicted), disease state, and the interaction of age and disease state included as factors. This methodology aimed to control for age, as the average age of the diseased and healthy samples differed in the majority of studies. The reported p-values for all the diseases correspond to the p-value of the factor “disease state” in the ANOVA.

Mortality association in LBC1921 and LBC1936

The Lothian Birth Cohorts of 1921 (LBC1921) and 1936 (LBC1936) are follow-up studies of the Scottish Mental Surveys of participants born in 1921 and 1936, respectively. The study was initially set up to study determinants of non-pathological cognitive aging [84, 85] and peripheral blood samples were analyzed by 450 k Illumina BeadChips. Only samples from the first waves of the LBC1921 and LBC1936 were considered, to keep all the samples around the same chronological age (79.11 ± 0.59 years for LBC1921 and 69.56 ± 0.84 years for LBC1936)—thus, in these cohorts chronological age has little impact on the association with mortality analyses. All the raw data in the form of IDAT files were processed in R with the ENmix package. Quality control was performed, probes with detection p-values greater than 0.01 were filtered out, and probes and samples with more than 10% of missing values were excluded. ENmix_oob background correction [86], RELIC dye-bias correction [87], and RCP probe-type bias adjustment [88] were performed. To minimize the potential influence of fatal acute illnesses on the methylation measurements, deaths that occurred during the first 2 years of follow-up were excluded from the analysis. After the preprocessing filter for the first wave, a total of 374 samples from the LBC1921 and 721 from the LBC1936 remained.

Mortality status was ascertained from data linkage using dates of death (which were converted to age in days at death by the LBC research team) from the National Health Service Central Resister, provided by National Records of Scotland. Time of the events for the Cox models was defined as the interval between [age in days at census wave1; age in days at death]/365.25 in case of a death or [age in days at census wave1; age in days at last census]/365.25 in case of census. At time of last census, 367/374 (98.13%) and 285/721 (39.53%) had died. Delta age was calculated as predicted age (with the WKDE model) − chronological age. Cox proportional hazards regression models were performed (adjusted for age and sex), to assess the association between delta age and mortality and between variation score and mortality. Additional regression models were tested for delta age and variation score that also considered the following parameters as covariates: smoking status (never, former or current), alcohol consumption (units per week), counts in blood for basophils, eosinophils, neutrophils, leukocytes, monocytes, white cells and platelets, triglycerides, total serum cholesterol and vitamin B12 levels in blood, heart rate (information only available for LBC1921), sitting diastolic pressure, body mass index, grip strength, forced expiratory volume, telomere length, activity lifestyle score (LBC1921), physical activity level (LBC1936), self-perceived life quality, self-perceived health status, and self-perceived life enjoyment. Additional details regarding the cohorts and the different variables that were measured can be found in the Cohort Profile Update [85]. The statistical software R was used to conduct all analyses, employing the 'survival package for the Cox models. For fitting the Cox proportional hazards regression models, the function coxph was used. Significance of the relationship between Schoenfeld residuals and time for the basic adjusted (age and sex) and fully adjusted (all additional covariates mentioned above) models was calculated with the function cox.zph.