Early research has unveiled a profound link between DNA methylation and aging, highlighting its significant role in the aging process [18, 19]. This groundbreaking insight spurred researchers to develop the first epigenetic clock, a remarkable tool created using saliva samples in 2011 [20]. This innovation marked a pivotal moment in epigenetics. Subsequently, Hannum et al. [21] published their 71 site methylation clock trained on blood followed by Horvath’s seminal work on a pan-tissue epigenetic clock on 353 CpG sites, demonstrating that the chronological age of an organism can be estimated based on epigenetic modifications, particularly DNA methylation patterns. Unlike chronological age, which simply counts the number of years a person has lived, epigenetic age aims to reflect the biological condition of the body and its systems, potentially providing a more accurate indicator of aging and health status. Many subsequent studies have aimed to refine and enhance these predictions (e.g., [22, 23, 24]). Horvath’s model, using linear regression on DNA methylation data, provided a robust framework for estimating chronological ages with impressive accuracy [25].

In 2021, Galkin et. al. [22] introduced DeepMAge, a deep neural network model that improved prediction performance over Horvath’s original model, particularly in blood samples. They employed gradient-based feature selection followed by a sequential neural network with four hidden layers, each containing 512 neurons. They used pathway analysis on the 1000 selected CpG sites by feature selection. They also emphasized the need for multimodal dataset, a lacking that still prevails as of today. As will be clear later, despite the current lack of such comprehensive datasets, our proposed model shows promise as a multimodal model.

Similarly, Levy et. al. [23] achieved encouraging results in age prediction using a multi-layer perceptrons. Another significant development, AltumAge [26], employed multi-layer perceptron layers and experimented their model performance across multiple tissues. They further utilized SHAP (Shapley Additive Explanations) [27] to interpret the contributions of different CpG sites towards age prediction. AltumAge also used interaction scores obtained from SHAP to demonstrate relationships between hyper-and hypomethylating CpG sites. Note that, hypomethylation and hypermethylation of DNA refer to relatively less or more methylation than in some standard DNA. Additionally, they identified pathways involving their top CpG sites. Thus the authors attempted to interpret and explain the results instead of just using AltmAge as a black-box. As will be shown later we also take a big stride along this dimension and hence we use AltumAge as the current state of the art for benchmarking.

Ying et. al. [24] proposed using epigenome-wide Mendelian Randomization [28] (EWMR) on 420,509 CpG sites to identify the sites that are causal to twelve aging-related traits, and then running the selected sites through an elastic net regression [29] model. A notable concern raised by Ying et al. was that correlation, rather than causality, dominated the results due to the exclusive use of methylation data. Therefore, they focused on using other biological data about the CpG sites for feature selection to achieve a more biologically significant outcome. This brings up the question of whether we could have a generalized formulation that would incorporate any additional feature (e.g., the ones mentioned by Ying et. al.) in the model under consideration, assuming of course, that the data of all relevant sites are available. As will be clear later, our model, GraphAge, naturally lends itself to such a formulation.

Before concluding this section, a brief discussion on a related interesting concept, called co-methylation is in order. Affinito et al. [30] showed that co-methylation depends on spatial and structural information and is related to the nucleotide distance among CpG sites. The co-methylation of closer CpG sites have more influence than that of those that are further apart. Also, research findings indicate that after DNA replication, De Novo genes, depending on specific sites, may or may not exhibit the same methylation rates as their parents, suggesting that methylation is influenced by structural and positional information [31]. This has motivated the researchers to model co-methylation networks along with constructing models to predict aging. For example, Wu et al. [32] created a co-methylation network based on MeRIP-Seq data by taking correlations among all sites and visualizing it using Cytoscape [33].

All previous clock models had fairly good accuracy but they did not adequately take into account crucial structural information and the dynamic interactions between CpG sites through co-methylation. This limits our understanding of interaction among these sites thereby preventing us from achieving a comprehensive interpretation. On the other hand, works that tried to model the co-methylation network did not take into account how CpG sites influence the aging process and which co-methylation edges/CpG sites are important to give a meaningful interpretation.

We propose to solve both the problems by utilizing all possible given information we currently have available regarding the CpG sites to formulate a graph where nodes encapsulate the information about CpG sites and edges represent all relationships among them thereby capturing the structure of co-methylation network. Then, we train it using the power of Graph Neural Networks, albeit in a resource constraint setting. Thus, we present the first age prediction model, GraphAge, that attempts to leverage all available structural information [34]. And despite our resource constraints, GraphAge was able to “come of age” by outperforming AltumAge in both dimensions of efficacy and interpretatbility. Notably, instead of a pan-tissue training approach, we opted for a tissue specific approach because gene expression follows a tissue specific pattern [35, 36]. For interpretation purposes, we leverage the Graph Neural Network Explainer [37].

For the sake of better comprehension, here we give a brief overview of our methods (Fig. 1); the details can be found in Section 4. We constructed a graph where nodes are CpG sites and edges represent their relationships. We selected CpG sites common across all methylation data platforms, using structural information from NCBI GEO [38]. We transformed tabular methylation data into a graph, with node attributes including methylation beta values, CpG island information, base pair positions, and normalized distances from transcription start sites. Edge features included co-methylation values, same chromosome, and gene indicators. We applied universal, chromosome-specific, and distance-based thresholds for edge filtering, allowing more flexibility for closer CpG sites with lower co-methylation values [30]. Our model used a Principal Neighborhood Aggregation (PNA) [39] convolutional layer in a Graph Neural Network (GNN) followed by a fully connected layer for age prediction.

For interpreting the model, we used GNN explainer [37]. We averaged node and edge importance across age groups for a comprehensive understanding and visualized these using Graphviz [40]. Importance analysis involved averaging node attribute importance across samples and plotting temporal patterns. We performed Enrichr analysis [41, 42, 43] on hypo- and hypermethylating genes, and identified nodes with significant upward or downward trends using linear regression.

We trained GraphAge and, as has been mentioned above, compared it primarily with AltumAge [26] for the reasons already mentioned above. We first split the data into a training and test sets (see Section 4), then conducted a 5-fold cross-validation on the former, and subsequently, we tested the best model on the latter. Fig. 3 reports the model performances. GraphAge overall performance is quite competitive with AltumAge, showing in fact a slight improvement over the latter on the test set in our constrained compute setting (more on this later). In particular, the overall MAE (MSE) of GraphAge is 3̃.207 (2̃5.277) as opposed to 3̃.296 (2̃8.310) of AltumAge registering a modest improvement of 2.7% (10.71%) over the latter. Furthermore, we grouped all samples (in the test set) into age-specific groups and measured the model performance within each group. Fig. 3(b) and Fig. 3(a), illustrate the age error in different age (and sex) groups. We see similar but slightly higher values for MSE in female than male. The trend for MAE is similar across all age groups. A consistent observation across all models is a decline in performance with increasing age. For example, starting from an MSE of 1.53 in male babies aged 0, it grows incrementally to 10.24 for ages 0-20 and 21.95 for ages 20-45. This supports the claim that epigenetic noise increases as chronological age diverges from epigenetic age [3]. Fig. 3(d) demonstrates how both models predict epigenetic age in relation to true chronological age, showing that both are quite similar even in terms of age acceleration (i.e., the direction and value in which the predicted epigenetic age diverts from the chronological age.) on healthy datasets.

Additionally, as an interesting case study, we tested both GraphAge and AltumAge on unhealthy datasets, i.e., in this experiment, the unhealthy samples were considered as unseen test samples for both the models. While this analysis also revealed similar age acceleration across all age groups for both the models (Fig. 3(c)), there are indeed some interesting points that are worth-discussing as follows. For postmenopausal women with ovarian cancer, GraphAge shows, on average, a (slightly) higher age acceleration, i.e., more in the positive direction than AltumAge. To elaborate, for age group 45-55, AltumAge shows an age acceleration value of 0.670 while GraphAge registers a much higher value of 3.452 and for age group 55-65, the values are -2.821574 and -0.783119 respectively. So, for this disease, GraphAge in most cases predicts a higher age acceleration than AltumAge, which is expected in case of a disease (in general).

For the other two diseases, namely, schizophrenia and osteoporosis, a different phenomenon is noticed. In case of schizophrenia most samples are from age group 20-45. Here, the age acceleration values for both models are slightly negative: for AltumAge, it is -0.009141 and for GraphAge, -0.346602. This may be attributed to the fact that schizophrenia is a mental disorder and hence does not show much age acceleration. Finally for osteoporosis, for age group 45-55 (55-65), age acceleration values for AltumAge and GraphAge are -6.786402 (-6.978263) and -2.500751 (-3.728412), respectively. So, both models suggest a lack of age acceleration and GraphAge is more closer to zero. While this may sound surprising and counter-intuitive, this actually is inline with the current knowledgebase as no evidence for accelerated epigenetic age in blood of osteoporotic patients were found as reported in [44].

As has been alluded to above, we had to work in a resource constrained setting. To this end, Fig. 3(e) indicates that GraphAge performance is influenced by the threshold value (see Section 4) used to filter edges from the underlying graph. Lowering the threshold value, which includes more edges in the underlying graph structure, enhances performance, albeit at the cost of increased computational cost. Just to provide an idea, the time for training one epoch increases from 4̃ min to 8̃ min when we decrease the threshold value from 0.8 to 0.7. Informatively, as we work in a constrained computational setting, we could not continue lowering the threshold value beyond the value reported here. However, Fig. 3(e) promises further performance improvement if more computational power could be made available.

We further remark that, in parallel to performance improvement, the lowering of the threshold value also promises to provide more insightful interpretation as more edges are included in the underlying graph. As is evident from the Supplementary Materials (see Section Supplementary Materials), decreasing the threshold value from 0.75 to 0.70 increases the edges by 90% thereby potentially making the graph more informative.

Datasets from NCBI [77] and EBI [78] were used. Datasets were in two different databases - Gene Expression Omnibus (GEO) and Array Express. Because Gene expression patterns are tissue dependent, a tissue type (blood) was selected. They contain both healthy (3707) and unhealthy samples (624). Only healthy data points were used for model training and testing. Among the unhealthy samples, we have postmenopausal women with ovarian cancer (266), schizophrenia (326) and osteoporosis (32).

Missing values were imputed with K-nearest neighbors (with $K=5$). The dataset were created using two different technologies. Therefore, they were not compatible. To fix this issue, a normalization method called Beta-mixture quantile normalization (BMIQ) was used [79].

We trained the model on healthy blood tissue datasets of 3707 samples where we first split the data into a test data of 756, i.e., 20% of the data and the rest where used for 5 fold cross validation. With the final model, we experimented to see the performance of the models in 8 age groups, namely, 0, 0-20, 20-45, 45-55, 55-65, 65-75, 75-80 and 80+ to examine how they performed in various phases of life. Middle ages were more segmented to see how changes occur in those phases of life.

As this is a regression task, we used Mean Absolute Error (MAE) and Mean Squared Error (MSE) as our main evaluation metric. MSE measures the average of the squares of the errors, giving more weight to larger errors and thus being useful for identifying models that produce significant errors. MAE measures the average magnitude of errors in a set of predictions, without considering their direction, providing a clear indication of the average error and being less sensitive to outliers compared to MSE. Based on prior literature, we also used another metric, age acceleration which is defined as predicted age minus the real age, to understand the direction and value in which the predicted epigenetic age diverts from chronological age.

In our graph, nodes are the CpG sites and edges are their relationships. Following the structure of AltumAge [26], we selected CpG sites that were common across all types of methylation data platforms. The structural information for these CpG sites was obtained from the public website NCBI GEO’s supplementary file [38]. Since the information of distance and position is chromosome specific, CpG site information was prepared based on their chromosomal location, ensuring that their positional and distance information were chromosome-specific.

For our model, nodes in the graph included the following attributes:

• •

  Methylation beta values for each CpG site.

• •

  Boolean value for being inside a CpG island (CpG islands are long sequences of Cytosine and Guanine).

• •

  Length of the CpG island.

• •

  One-hot encoding of the next base pair (three such base pair information available).

• •

  Starting base pair position of the island if it is an island, else 0, normalized.

• •

  Ending base pair position of the island if it is an island, else 0, normalized.

• •

  Normalized distance from transcription start site; null values are given 1 to indicate they are the farthest.

• •

  Map_Info, which is the position of the CpG site for both island and non-island sites, is also normalized and added to the node attribute.

The most crucial part of the formulation is encoding the relationships among the edges. Edge features are:

• •

  Co-methylation value among two CpG sites, using methylation values of CpG sites from training data.

• •

  Boolean value if two CpG sites are on the same chromosome.

• •

  Boolean value if two CpG sites are of the same gene.

We applied three thresholds for edge filtering:

• •

  A universal threshold to filter all edges.

• •

  A lower threshold for edges on the same chromosome (Secondary Threshold). We have taken this to be 0.2 lower than the universal threshold.

• •

  A distance-based threshold within the same chromosome to allow even lower co-methylation edges to be included if they are closer (Tertiary Threshold). We have taken this to be 0.4 lower than the universal threshold.

The rationale behind using three thresholds is based on the influence of closer CpG sites, which tend to have a stronger effect on co-methylation [30]. The concept of distance is applicable only within the same chromosomes. This approach allows for more flexibility in selecting edges that are closer but may have slightly lower co-methylation values. These parameters must be chosen carefully, balancing speed and accuracy requirements.

Our model comprises of GNN layers followed by a fully connected layer. So after formulating the graph, we used DataLoaders from PyTorch[80] to create data loaders. The graph was fed into a Graph Neural Network (GNN) in a stochastic manner. We chose the Principal Neighborhood Aggregation (PNA) convolutional [39] layer for our GNN because of multiple aggregators to capture diverse aspects in graphs. After the GNN layer, the output was passed through a fully connected layer to obtain the final age prediction.

We used 5-fold cross-validation to find the best model and then tested the best model on an test dataset. For other models, we followed exactly the same procedure for a fair comparison.

We used the GNN explainer [37] to interpret the model. Two types of explanations were generated:

• •

  Importance of each node attribute.

• •

  Importance of nodes and edges to identify connected subnetworks that work in conjunction.

GNN explainer gives us individual sample explanations which we save. We take these individual sample importance of each site and average them based on the age groups that we used for testing. We then create a graphical representation of the average node and edge importance for each age group. We also color-coded these subnetworks for visualization. Importantly, we generated individual explanations for each sample and averaged the explanations across certain age groups to provide a more comprehensive understanding.

We employed vizgraph[40] for constructing visualizations:

• •

  Red node: hypomethylating CpG sites. With age, the methylation value decreases, and methylation upregulates gene expression.

• •

  Blue node: hypermethylating CpG sites. With age, the methylation value increases, and methylation downregulates gene expression.

• •

  Green circle around node: the importance of that particular CpG site increases with age.

• •

  Yellow circle around node: the importance of that particular CpG site decreases with age.

We annotated the graph for better understanding

After using the GNN explainer, we analyzed the importance of node attributes by averaging them across all samples. We also plotted the importance of these attributes over different age groups to observe their temporal patterns. For node and edge importance, we divided all samples into various age groups and averaged the importance within each group. We then eliminated nodes with zero importance and edges with low importance. The resulting subnetworks were visualized using Graphviz. We also created individual visualizations for subnetworks that contained a certain number of CpG sites. Additionally, we identified genes from these subnetworks and categorized them based on whether they were hypo- or hypermethylating. We performed Enrichr analysis [41, 42, 43] on these two groups of genes separately.

Furthermore, we used linear regression to find the slope of each node’s importance with respect to age. We identified the top 10 nodes with the most positive slopes as upward trending and the top 10 nodes with the most negative slopes as downward trending. Nodes with small slope values were considered to have no significant trend.

We used kaggle free tier version using the GPU P100. We have used most implementation from pytorch geometric library [80]. All the datasets and other information are all present in the following notebooks.

Our code for predicting and compiling explanation: https://www.kaggle.com/code/sakibahmed91/graphage-final-version-for-result-and-explanation

For analysing the explanation : https://www.kaggle.com/code/salehsakibahmed/graph-age-exaplanation-analysis

All codes and result compilation are available here: https://github.com/bojack-horseman91/GraphAge

Our paper includes all results from using threshold 0.7 but in Fig. 9 we show the results for threshold 0.75, secondary_threshold 0.73 and tertiary_threshold of 0.71. One other difference is that the number of epoch to train GNN explainer in this case was 150 while in the explanation of threshold 0.7 epoch was 100 due to longer training time. As can be observed from the figure, particularly in contrast with Fig. 4(b), although node attribute importance remains quite similar the edge importance scores change. We can see from Fig. 10 that by decreasing the threshold from 0.75 to 0.7, on an average for each age group, we see an increase of 812.75 (90.162%) in the number of edges. So by decreasing the threshold we are encoding more interactions/relationships into the graph structure thereby increasing the potential for a more comprehensive understanding of the complexities of aging while decoding (i.e., interpreting) it. Table 1 provides a detailed comparison in terms of MAE among GraphAge, AltumAge, and Horvath’s model including results for multiple threshold values of GraphAge. Also, Table 2 presents a comparison between GraphAge and AltumAge in different age groups.

There are multiple studies in the literature on mice that linked the identified important genes with aging. TSC1, a tumor suppressor gene [81], is shown to be associated with accelerated retinal aging in knockout mice [82]. Moderate lifelong overexpression of TSC1 has been demonstrated to improve health and survival in mice [83]. PDPK1 signaling has great influence on determining the reproductive aging and the length of reproductive life in females [84]. On the other hand, aging affects leptin (LEP) actions differently. There is a decline in hypermetabolic responses but increased sensitivity in lean rats. Older rats are more prone to obesity-induced leptin resistance, while calorie restriction enhances leptin responsiveness, particularly in older rats [85].

These are all the dataset that we used in our entire process. For clarity of user we have added a small description for each of the dataset.

• [label=]

• •

  E-GEOD-51388: Blood samples taken longitudinally from monozygotic twins. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-52588: Blood samples from subjects with or without Down syndrome. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-53128: Whole blood samples. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-53740: Blood samples from progressive supranuclear palsy patients, frontotemporal dementia patients, and healthy controls. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-54399: Cord blood and placenta. All ages were encoded as zero age. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-54690: Blood samples from subjects with or without dietary flavanol supplementation. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-56553: Peripheral mononuclear blood cells of asthmatic patients. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-57484: Whole blood samples from normal and obese children. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  E-GEOD-58045: Whole blood samples. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  E-GEOD-59509: Whole blood, saliva, menstrual blood, vaginal swab, and semen samples. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-59592: Blood from infants exposed to varying degrees of aflatoxin B1. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-62219: Longitudinal peripheral blood leukocyte samples from infants. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-64495: Whole blood from subjects with or without developmental disorder syndrome X. Non-healthy patient samples were separated for further analysis and not included in training, validation, or testing. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-64940: Cord blood samples from newborns. All samples were encoded as zero age. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  E-GEOD-65638: Whole blood samples from twins. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-67444: Neonatal blood samples. All ages were encoded as zero age. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-67705: Blood samples from HIV+ and HIV- subjects. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-71245: Blood samples from different types of blood cells. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-71955: CD4+ and CD8+ T-cell samples from subjects with Graves’ disease or healthy controls. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-72338: Neutrophils and monocytes from patients with tuberculosis and household controls. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-77445: Whole blood samples from subjects with different stress levels. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-79056: Cord blood samples. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-GEOD-83334: Whole blood and cord blood from newborns and infants measured longitudinally. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  E-MTAB-2344: White blood cell samples in patients with stroke and/or obesity and healthy controls. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  E-MTAB-2372: White blood cell samples from obese patients subject to two different diets. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  GSE19711: Whole blood samples from patients with or without ovarian cancer. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE20236: Whole blood samples. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE20242: CD4+ T-cells and CD14+ monocytes samples. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE34257: Blood cord samples and whole blood samples from 9-month-old infants. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE36642: Blood mononuclear cell, human umbilical vascular endothelial cell, and placenta samples from monozygotic and dizygotic twins. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE37008: Peripheral blood mononuclear cell samples. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE41037: Whole blood samples in schizophrenia patients and healthy subjects. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE49904: Blood buffy coat samples. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE56606: CD14+ monocyte samples from monozygotic twins discordant for type 1 diabetes. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE57285: Whole blood samples from women with BRCA1 wild-type or mutants, and with or without breast cancer. The platform used was Illumina’s Infinium 27k Human Methylation Beadchip.

• •

  GSE69176: Umbilical cord blood samples from newborns. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.

• •

  GSE99624: Whole blood samples from osteoporotic and healthy control patients. The platform used was Illumina’s Infinium 450k Human Methylation Beadchip.