Precious1GPT: multimodal transformer-based transfer learning for aging clock development and feature importance analysis for aging and age-related disease target discovery

Performance of the transformer-based multimodal aging clock

In the current study, we have developed a comprehensive pipeline, as illustrated in Figure 1. The pipeline consists of several key steps, including training a multimodal transformer-based regressor on normal sample data for age prediction and subsequently using the learned weights to fine-tune a transformer-based classifier for distinguishing between case and control samples. Next, we perform gene prioritization by employing the feature importance values obtained from the regressor to rank genes based on their relevance to aging and utilizing the importance values from the classifier to rank genes in terms of their relevance to both aging and disease. Finally, we analyze the resulting gene lists using the PandaOmics TargetID Platform to gain insights into potential targets for age-related diseases. In this study, we implemented a transformer-based architecture to accommodate both numerical and categorical data as input for our predictive model. This strategy, which we call Precious1GPT, enables the construction of multimodal classifiers and regressors that can effectively process diverse data types, such as RNA-seq expression data and epigenomics methylation data taking into account data type and tissue type. Consequently, our model demonstrates proficiency in age prediction and case-control classification, showcasing its versatility in handling multifaceted inputs.

We employed Optuna [23], a hyperparameter optimizer, to optimize the parameters for each model. We optimized L1 and L2 regularization, the activation function, dropout value, batch size, the number of neurons in hidden layers, and the gradient update used. The final regression metrics for the optimized models are shown in Table 1, calculated for the epigenetic and expression sub-datasets and for the whole dataset, respectively. The metrics were calculated by splitting the data to train (80%) and test (20%) set with stratification by sample tissue. For five-fold cross-validation stratified by tissue and data modality, results are shown in Supplementary Table 1. Metrics for individual tissues are shown in Supplementary Table 2. The methylation data subset, expression data subset, and all the test data combined were also calculated for each metric. Learning curves depicting MAE during training are shown in Supplementary Figure 1.

Important genes for age prediction

SHapley Additive exPlanations (SHAP) values [24] are a technique for explaining the output of machine learning models. From the regressor model, we obtained a list of features ranked by their SHAP values representing their importance for age prediction (Supplementary Table 3). Pathway enrichment analysis was subsequently performed for the top-100 genes ranked based on the SHAP values, which showed that these genes are implicated in multiple pathways associated with aging and age-related diseases (Table 2).

Identification of potential targets for age-related diseases through feature importance analysis

Utilizing the feature importance analysis based on SHAP values, we generated a list of genes associated with aging (Supplementary Table 3). We then compared these genes with known drug targets in our in-house database to identify potential therapeutic interventions for 4 selected age-related diseases, namely idiopathic pulmonary fibrosis (Supplementary Table 4), chronic obstructive pulmonary disease (COPD) (Supplementary Table 5), Parkinson’s disease (PD) (Supplementary Table 6) and heart failure (Supplementary Table 7). Evaluation metrics for case-control classifiers are shown in Supplementary Table 8.

We adopted a transfer learning approach to identify genes involved in disease development in the context of aging. We first trained a DL-model as a regressor to predict age using an age dataset. Subsequently, we fine-tuned the model by re-training it as a case-control classifier while keeping the previously learned weights frozen, with the exception of the last layer. The SHAP values generated from this analysis were then used to determine the relative importance of molecular features in driving disease development in the context of aging. This allowed us to identify specific genes that are involved in disease development in the context of aging and determine their relative importance. To establish a baseline, we trained the same classifiers on the complete feature set. For a number of diseases, we observed a slight but significant increase in classification metrics (Supplementary Table 8). These results indicated that the last layer of the neural network, which was trained to predict biological age, contains sufficient information to differentiate between case and control.

To validate the performance of our model and to establish its ability to accurately estimate age based on methylation data, we acquired two methylation datasets from the Gene Expression Omnibus (GEO) repository - one on cells that were reprogrammed to induced pluripotent stem cells (iPSC) (dataset GSE54848) and the other on fibroblasts of the developing fetuses (dataset GSE76641). The selection of the independent dataset on the developing cells, along with reverse aging cells data, allows us to provide another level of validation evidence confirming the potency of the trained model. These datasets were processed using the same methods as our primary dataset on age-related methylation, and we used the model to predict age. In order to increase the robustness of our validation, we used the same processing methods for both the iPSC and fetus datasets as we had for our primary dataset. This ensured consistency and minimized the possibility of any discrepancies or variations in our results due to differences in processing methods. The results were consistent with expectations, as iPSCs become younger during induction (Figure 2, Left) and fetal tissue becomes older during development (Figure 2, Right).

Manual analysis of the resulting targets

Utilizing a transfer learning approach, we have built aging-aware case-control classifiers and extracted feature importance values from them. The lists of top-200 genes ranked by expression classifiers were retrieved (Supplementary Tables 4–7) and considered as a starting point for further target identification and prioritization techniques offered by the AI-powered PandaOmics platform to propose a list of promising novel targets for age-related diseases. According to PandaOmics TargetID platform, APLNR was ranked top-20 for all 4 diseases, while IL23R was ranked top-20 for COPD, PD, and heart failure (Figure 3, Supplementary Figures 2–4). APLNR and IL23R were therefore selected as the most promising targets for treating multiple age-related diseases. In general, APLNR, a receptor for Apelin and Elabela peptide ligands, is involved in regulating several important physiological processes, including cardiovascular function, fluid balance, and metabolism. IL23R is a receptor for the pro-inflammatory cytokine IL-23 and is associated with chronic inflammation, which is considered to be one of the hallmarks of aging [25]. Accumulating evidence demonstrated that the effectiveness of our approach in addressing various age-related diseases, as documented in existing literature, provides further validation for our approach. Therefore, our unique approach with the amplification of PandaOmics allows us to identify various potential targets associated with essential aging-driven tissue dysfunction, which may be useful for the delay and treatment of multiple age-related diseases.

Data sources for multimodal aging clock development

To train our models, we used several datasets, including publicly available and in-house-built ones. For training age prediction models based on the epigenetic status of tissue, we employed 450 k Illumina Methylation array data from EWAS Data Hub [47] (number of samples = 8,374).

For building age prediction models based on the transcriptomics status of tissue, we employed RNA-Sequencing data from the GTEx project [48] (number of samples = 12,453). We trained our models in a tissue-agnostic fashion. For methylation, data distribution of samples across ages is shown in Supplementary Figure 5A and across tissues in Supplementary Figure 5B. For expression data, distribution of samples across ages is shown in Supplementary Figure 6A, across tissues in Supplementary Figure 6B.

For assessing prediction results and predicting disease targets using age-pretrained models, we used custom datasets from the PandaOmics software [49] for 4 selected age-related diseases: idiopathic pulmonary fibrosis, COPD, PD, and heart failure from where we have obtained samples annotated as carrying disease (case samples) and healthy ones (control samples).

To evaluate possible interventions to prevent senescence development, we have constructed datasets containing only features corresponding to genes for which approved drugs exist. For this, we have used an in-house constructed database of approved drugs and their targets based on the information from [50].

As input features for our age prediction models were either beta values averaged across CpG probes annotated as TS200 region from Illumina 450k Methylation Array for epigenomic data or TPM values for protein coding (genes) for expression datasets accordingly.

For the DNA methylation data, we obtained the ß-values from the CNCB data hub, where the raw data were obtained using the GMQN package developed by the CNCB [47].

In our study, we have opted for the TSS200 region as the most interpretable for age prediction. Following the aggregation of corresponding beta-values, we are left with approximately 14,000 features representing average methylation of proximal promoter regions.

This choice of region is based on its potential to offer a more accurate and comprehensive assessment of age-related changes in methylation patterns. The TSS200 region, situated within 200 base pairs upstream of the transcription start site, is known to play a crucial role in gene regulation [51]. As such, it is expected to exhibit significant age-related changes in methylation patterns, providing a robust basis for predicting biological age.

To further enhance the interpretability of our age prediction model, we utilized machine learning techniques to identify and select the most informative features from the 14,000-feature dataset. This allowed us to refine our model, ensuring that it captures the most relevant age-related methylation changes in the TSS200 region. In addition, the selected features can potentially shed light on the molecular mechanisms underlying aging, as well as the development of age-related diseases.

By focusing on the TSS200 region, we aim to not only improve the accuracy of our aging clock but also gain a deeper understanding of the complex relationship between DNA methylation, aging, and age-dependent diseases. This knowledge can then be used to develop targeted therapeutic interventions aimed at mitigating the impact of age-related diseases and improving overall health and quality of life in aging populations.

For gene expression models, we have used TPM values which were back-corrected using ComBat [52] and quantile-normalized using qnorm [53]. For performance evaluations, we employed a shuffle split stratified by sample tissue, leaving 20% of all data for the test set.

Transformer-based model for multimodal aging clock

Given the abundance of omics data available for various experimental conditions and the distinct challenges of predicting one omic data type from another, we propose the development of a transformer-based model to estimate sample age across different sample types. This multi-tissue, multi-omics transformer-based age prediction model aims to harness the power of deep learning to effectively integrate diverse data sources and improve age prediction accuracy.

Transformers have demonstrated remarkable success in various applications, particularly in natural language processing tasks. Their capacity to model long-range dependencies and capture complex relationships among features makes them well-suited for multi-omics data integration. By leveraging the transformer architecture, our proposed model is able to identify and exploit relevant information from different omics data types, such as genomics, transcriptomics, proteomics, and metabolomics, as well as different tissue types.

By developing a multi-tissue, multi-omics transformer-based age prediction model, we hope to enhance the accuracy and generalizability of aging clocks, ultimately contributing to a better understanding of the aging process and facilitating the development of targeted therapies for age-related diseases.

More formally about the transformer model. Let X be the input matrix of size (N, D) containing epigenetic (methylation) and expression data, and Z be the input matrix of size (N, M) containing categorical data. Let Y be the output matrix of size (N, 1) containing the predicted age values. Since we used TabTransformer (LINK), we represented the model as a function F that takes X and Z as input and returns Y as output: Y = F (X, Z). The model consists of multiple linear and attention layers, each of which applies a set of learnable parameters to the input and produces an output. We represent each layer as a function G that takes an input matrix A and a set of learnable parameters W and b, and produces an output matrix:

$$\begin{array}{l}B:B=G\text{(}A,W,b\text{)}=\\ \mathrm{Re}LU\text{(}AW+b\text{)}B=G\text{(}A,W,b\text{)}=\\ \mathrm{Re}LU\text{(}AW+b\text{)}\end{array}$$

The TabTransformer [54] model consists of multiple layers, including self-attention layers and feedforward layers. Let A be the output of the previous layer, and let W_q, W_k, and W_v be learnable weight matrices of size (D, d_k). The self-attention layer computes the attention matrix A’ as follows:

$$A^{\prime }=soft\max \left(\frac{A{W}_q{(W_k)}^T}{\sqrt{d_k}}\right)W_{v}O=A^{\prime }{W}_v$$

where W_q, W_k, and W_v are weight matrices of size (D, d_k). The feedforward layer computes the output matrix H as follows:

$$H=G\left(O,W_1,b_1\right)=ReLU\left(O{W}_1+b_1\right)$$

where W₁ is a weight matrix of size (h, d_h), and b₁ is a bias vector of size (1, d_h). We repeat the self-attention and feedforward layers multiple times to create a deep TabTransformer model. Next, we apply a feedforward layer with weight matrix W₁ and bias b₁ to the output of the self-attention layer: H = ReLU (OW₁ + b₁). Finally, we apply a linear layer with weight matrix W₂ and bias b₂ to the output of the last feedforward layer to produce the final output matrix Y:

$$Y=H{W}_2+b_2$$

To train the TabTransformer model, we used mean squared error (MSE) loss that measures the difference between the predicted age values and the true age values. Let $\widehat{Y}$ be the predicted age values and Y_true be the true age values:

$$L=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\text{(}{\widehat{Y}}_l-Y_{true,i}\text{)}}^2}$$

Overall, the experiment involves training a TabTransformer model on a dataset containing epigenetic (methylation) and expression data, as well as categorical data (tissue and dataset type), to predict age values. The model consists of multiple layers, including self-attention and feedforward layers, and is trained using a loss function such as MSE.

In the present study, we utilized PyTorch Tabular [55], which is built on top of PyTorch, for all the work with the transformer. PyTorch Tabular provides a highly optimized and efficient way of handling tabular data with PyTorch. We used PyTorch Tabular’s various functions and modules to preprocess the input data, construct the transformer architecture, and train it on the data.

In the model, all the hyperparameters, such as learning rate, dropout rate, number of hidden layers, and activation functions, were determined through the use of Optuna [23], a hyperparameter optimization framework. The values of these parameters were chosen based on their performance during multiple rounds of training and validation. Based on our experiments, the optimal hyperparameters for TabTransformer are a model architecture with hidden layers of size 128, 2048, and 128, dropout probability of 0, “ELU” activation function, a learning rate of 0.00023, “AdamW” optimizer, weight decay of 0, and a batch size of 96.

Transfer learning for case-control classification

To build models which take into account both senescent and clinical status, we employed a transfer learning strategy. Initially, we trained a model as a regressor to predict age in an age dataset. Subsequently, we froze the model weights, excluding the final layer, and re-trained the resulting model as a case-control classifier. The derived SHAP values indicate the significance of each molecular feature in disease development concerning aging. This approach enabled us to determine the relative importance of various molecular features in driving disease development within the aging context. The pipeline is depicted in Figure 1.

Feature importance analysis

SHAP values, a mathematical method of feature importance analysis that constitutes a robust and interpretable technique for elucidating the contributions of individual features in complex predictive models, is based on the Shapley value from game theory and involves computing the contribution of each feature to explain the final predictions of machine learning models [24]. SHAP values explain individual predictions and identify the most important features in the model. Additionally, the employment of SHAP values for feature analysis offers a considerable advantage in multimodal settings, where discerning the interplay between various factors is indispensable for the development of accurate and efficacious aging clocks.

Pathway enrichment analysis

Pathway enrichment analysis was performed on the list of top-100 genes ranked by their SHAP values obtained from the regressor model with the pathways available on the Reactome database. R package Enrichr was used to calculate the enrichment levels and p-values. Pathways with p-value < 0.05 were considered as significantly enriched.

Manual analysis of resulting targets

Combination of aging clocks and target ID represents an interesting approach to identifying targets for aging-associated diseases. To illustrate the applicability of this approach for target identification, we have investigated 4 diseases associated with aging: idiopathic pulmonary fibrosis, chronic obstructive pulmonary disease (COPD), Parkinson’s disease (PD), and heart failure in PandaOmics. Gene lists of top-200 genes ranked by expression classifiers were used in Target ID projects for the mentioned diseases, along with a filter for small molecules to identify potential druggable proteins across these lists.