Development of a novel transcriptomic measure of aging: Transcriptomic Mortality-risk Age (TraMA)

Training TraMA

Because we are interested in developing a measure that is accurate and portable to other human datasets, we restricted the set of genes used for training to coding genes with relatively high expression in human venous blood. To accomplish this, of the 50,611 transcripts that were measured and were successfully mapped in HRS, we restricted ourselves to the 19,291 protein coding genes. We further restricted ourselves to genes with a mean count per million greater than 3 in the total HRS sample, leaving 10,964 genes. We also included chronological age in years and sex/gender as features to reduce age and sex/gender bias in the selection of genes (i.e., they are covariates along with the 10,964 genes used to train TraMA).

We then randomly split the HRS sample into training (N = 1801) and testing (N = 1794) subsamples. Descriptive statistics for each subsample are shown in Table 1A and 1B. 222 participants in the training data died during the 4-year follow-up. 197 participants in the testing data died during follow-up. We ran elastic net models predicting 4-year mortality in the training sample. Hyperparameters, including the alpha and lambda penalty terms, were selected using a grid search procedure with 5x cross validation. Log₂ adjusted counts per million (log₂cpm) were used in all analyses, including the Health and Retirement Study (HRS), the Long Life Family Study (LLFS), and publicly available data described below. This procedure selected an alpha of 1 (equivalent to LASSO regression) and a lambda of 0.0198. This model selected (i.e., did not reduce regression coefficients to 0) 35 genes and age. Gene names, Ensembl IDs, and coefficients are shown in Table 2. TraMA scores were transformed to have a mean and variance equivalent to the HRS training set chronological age. In the training data, Harrel’s C-index predicting survival from the TraMA score was 0.835, suggesting very good fit.

Gene ontologies, associated traits, pathway analysis, and functional enrichment analysis

Validation in the health and retirement study testing sample

Validation in the Long Life Family Study

To ensure portability of this measure, and as a further check for robustness, we also validated this measure in an external cohort that includes a large number of older adult humans, the LLFS. Using mixed effect Cox proportional hazards models, we regressed time to death in each sample on TraMA, controlling for batch (Model 1); adding age, race/ethnicity (in HRS; because nearly all LLFS participants were White, we did not control for Race/Ethnicity in LLFS), and sex/gender (Model 2); and adding RNA-based cell type (Model 3). All LLFS models were adjusted for family relatedness as a fixed effect. Results are shown in Figure 1D. TraMA was significantly associated with time to death in both LLFS and HRS with all controls included.

The hazard ratio for LLFS is higher than for HRS. This potential because LLFS participants were selected for their longevity. LLFS also had twice the number of mortality events as HRS. Because the LLFS includes a sample of older adults, their children, and spousal controls, we additionally ran analyses only in the sample of older adults (referred to in the LLFS as the proband generation). These results largely replicated those shown in Model 3 above, with a hazard ratio of 2.48 and p-value less than 0.001 for the proband generation, 4.76 and p-value less than 0.01 for the offspring generation, and 1.63 and p-value less than 0.01 for the spousal controls.

Finally, LLFS has also recently processed RNA data from visit 2 (N = 1,263). We additionally regressed mortality hazard in 2020 on TraMA scores from time 1 and 2, with covariates from Model 3 above. When TraMA scores from both time points were entered in the model, visit 2 TraMA was associated with mortality hazard (HR = 1.66, p < 0.05) and time 1 TraMA was not (HR = 1.13, p > 0.05), suggesting that residual change in TraMA score is associated with mortality hazard, and that longitudinal change in TraMA may track mortality risk.

Validation in small and clinical samples

To be a maximally valuable measure of the aging process, TraMA should be useful not only in large representative samples, but also in small specialty and clinical samples. We therefore additionally validated four publicly available datasets from the Gene Expression Omnibus (GEO) with RNA-seq data from whole blood and information about chronological age. First, to validate expected associations with health behaviors, we estimated TraMA in data from 454 current and 767 former smokers from the COPDGene Study (GEO series GSE171730), including non-Hispanic White and African American people between the ages of 45 and 80 in the US (shown in Figure 2A) [40, 41]. In these data, TraMA was positively associated with smoker status and number of smoking pack years, and negatively associated with lung functioning assessed with forced expiratory volume in one second (FEV1) predicted percentage. Current smokers had an estimated TraMA 2.83 years older than former smokers controlling for age, race, sex/gender, RNA-based cell type, and batch.

TraMA was also associated with diabetes status in a sample of 43 healthy participants and 39 participants with type 1 diabetes (GEO series GSE123658; results shown in Figure 2B) [42]. In these data, controlling for age, sex/gender, RNA-based cell type, and batch, participants with diabetes had a predicted TraMA 2.12 years older than healthy controls.

We also validated in data from the Mount Sinai Crohn’s and Colitis Registry (MSCCR; GEO series GSE186507) [43] including 821 participants with inflammatory bowel disease (IBD; 432 with Crohn’s disease and 389 with ulcerative colitis) and 209 healthy controls (results shown in Figure 2C). IBD status was not significantly associated with TraMA after statistically controlling for age, sex/gender, and RNA-based cell type; however, this association was significant without RNA-based cell type (log odds = 0.07, p < 0.001), and Crohn’s disease participants had elevated TraMA compared to healthy controls with all covariates (b = 1.00, p < 0.05), this may be because Crohn’s disease patients were more likely to have active IBD (Harvey-Bradshaw index (HBI) ≥5; p < 0.01). Among all IBD patients, participants with clinically active IBD (according to a physician’s evaluation) had higher TraMA, and TraMA was positively associated with IBD severity using the Simple Endoscopic Score for Crohn’s Disease (SESCD).

TraMA was associated with sepsis, compared to healthy controls, in a sample of 348 sepsis patients and 44 healthy controls (GEO series GSE171730; results shown in Figure 2D) [44] in a regression including age, sex/gender, and RNA-based cell type as covariates. Compared to healthy controls, sepsis patients had an estimated TraMA 10.31 years older. TraMA was also associated with Sequential Organ Failure Assessment (SOFA) scores among sepsis patients. The association between TraMA and mortality among sepsis patients approached significance (p = 0.09), though this association was significant in a model without cell type (log odds = 0.05, p < 0.001).

Validation with other platforms and model species

A truly portable measure of transcriptomic aging would also be valid in human data using other RNA measurement platforms (e.g., array data) and in model species. To that end we also estimated TraMA in GEO data using an Affymetrix array, an Illumina array, and in a mouse (Mus musculus) sample. Because expression values from arrays and from RNAseq have different distributions, we do not expect the means and variances of TraMA calculated from arrays to be meaningful. We therefore use standardized scores in all of these analyses.

We begin by analyzing blood samples from 1,013 human cancer patients and 1,832 control samples with RNA profiled on an Affymetrix array (GEO series GSE203024; results shown in Figure 3A). Participants with cancer had higher TraMA scores compared to those without cancer controlling for age, sex, and RNA-based cell type. TraMA was also associated with lupus diagnosis in a sample of 134 juvenile patients and 36 healthy controls, statistically controlling for age, sex, race, RNA-based cell type, and batch (GEO series GSE65391 [45]; results shown in Figure 3B). Thus, TraMA appears to still perform robustly in a juvenile sample.

Finally, we estimated TraMA in a sample of mice exposed to 1.5, 3, 6 or 10 Gy of gamma-rays or sham irradiated controls and sacrificed at either 1, 2, 3, 5 or 7 days after exposure (GEO series GSE124612) [46]. Greater levels of radiation were associated with more TraMA aging controlling for number of days after exposure, but this pattern largely disappeared after controlling for RNA-based cell type (mice exposed to 1.5 Gy were slightly lower on TraMA compared to controls and no other differences were significant, see Figure 3C). However, there appears to be a joint effect of time and radiation (Figure 3D). Mice exposed to 1.5 Gy were mostly not significantly different on TraMA compared to controls sacrificed on all days, except on day 5 they were slightly lower. Mice exposed to 3 Gy experienced a spike in TraMA age on days 2 and 3 that appears to dissipate by day 5, and they are slightly lower than controls on day 7. Mice exposed to 6 Gy showed a similar pattern of a spike on days 2 and 3 that dissipates. Mice exposed to 10 Gy, showed higher TraMA age on all days after day 1. We also generated volcano plots from regressions of each gene in all of the human GEO datasets on age (see Supplementary Figure 4).

Cohorts

The Health and Retirement Study (HRS) is an ongoing panel study of older adults since 1992 that is designed to be representative of older US adults when weighted. As part of 2016 data collection, venous blood was collected from a subsample of the HRS. 2.5 ml of blood was collected in PAXgene tubes from about 4000 participants. Total RNA extraction was performed on the QIACube semi-automated method using the PAXgene Blood miRNA Kit. Assays used 200-500 ng of RNA for each sample. All RNA species were extracted and stored for future use. RNA was extracted from only half a PAXgene tube to ensure RNA storage in a variety of formats. Ribosomal RNA and globin reduction was performed using the TruSeq stranded Total Library Prep Gold kit - Ribozero Gold kit. RNAseq was performed on a NovaSeq (Illumina Inc.) using 50 bp paired end reads. All samples were sequenced to a minimum depth of 20 M reads. RNA-Seq was successfully performed on 3685 participants. The HRS pipeline closely mirrors the TOPMed/GTEX RNA-Seq analysis pipeline with minor modifications. More information about RNAseq pipelines is available elsewhere [49] including the HRS website (https://hrs.isr.umich.edu/about).

The Long Life Family Study (LLFS) is a longitudinal sample of nearly 5000 participants from 539 families that were selected because of their exceptional longevity. There have been three waves of data collected 6-8 years apart. The first and second waves of data included blood collection. We use data from the first wave to align with HRS. More information is available at the LLFS website https://longlifefamilystudy.com/.

RNA sequencing for Visit 1 was performed using RNA extracted from PAXgene™ Blood RNA tubes, processed with the Qiagen PreAnalytiX PAXgene Blood miRNA Kit. Library preparation, quality control, and sequencing were carried out by the Division of Computation and Data Sciences at Washington University, using the nf-core/rnaseq 3.14.0 pipeline for read alignment, duplicate marking, and transcript quantification. Genes with low expression (fewer than 4 counts per million in at least 98.5% of samples) and those with significant intergenic overlap were filtered out. This resulted in a final dataset of 1,810 samples and 16,418 genes. For this study, we utilized RNAseq data from the LLFS dataset, with the filtered raw counts converted to a Log₂CPM (counts per million) scale for further analysis.

The COPDGene study (GSE171730) that is publicly available includes 454 current and 767 former smokers, including non-Hispanic White and African American men and women between the ages of 47 and 86 in the US. RNAseq was performed on whole blood using the Illumina HiSeq 2000 platform. More information is available on the COPDGene website (https://copdgene.org/). Information about current smoker status, pack years, and forced expiratory volume in one second (FEV1) predicted percentage, race, sex/gender, and batch are available.

GSE123658 is a sample of 43 healthy donors and 39 type 1 diabetes patients between ages 19 and 73. RNAseq was assessed in whole blood using Illumina NextSeq 500 or HiSeq 4000 platforms. Information about diabetes status, age, and sex/gender are available.

The Mount Sinai Crohn’s and Colitis Registry (GSE186507) includes 821 irritable bowel disease (IBD) patients and 209 healthy controls aged 19 to 82 recruited during an endoscopy appointment from December 2013 to September 2016. RNAseq was assessed in whole blood using the Illumina HiSeq 2500 platform. Information about IBD status, active IBD status (Harvey-Bradshaw index (HBI) ≥ 5), disease severity Simple Endoscopic Score for Crohn’s Disease (SESCD), age, and sex/gender are available.

GSE185263 is a sample of 348 sepsis patients and 44 healthy controls aged 18 to 96 from countries, including Australia, Colombia, the Netherlands, and Canada (sites in Toronto and Vancouver). RNAseq was assessed in whole blood using the Illumina HiSeq 2500 platform. Information about sepsis status, severity using Sequential Organ Failure Assessment (SOFA) scores, mortality, age, and sex/gender are available.

GSE203024 includes blood samples from 1,013 human cancer patients with 11 different types of cancer or colorectal polyps and 1,832 control samples without a cancer diagnosis with RNA profiled on Affymetrix U133 Plus 2.0 GeneChips. Expression values were log₂ transformed and missing values were set to the median. Affymetrix IDs were matched to ensembl IDs using biomaRt [37]. If an ensembl ID matched more than one probe, the mean of the values was taken. 34 of the 35 TraMA genes were available.

GSE65391 is a sample of pediatric lupus patients. Expression values were log₂ transformed and missing values were set to the median. Illumina IDs were matched to ensembl IDs using biomaRt [37]. If an ensembl ID matched more than one probe, the mean of the values was taken. 32 of the 35 TraMA genes were available.

GSE124612 is a sample of male C57BL/6 mice exposed to 1.5, 3, 6 or 10 Gy of gamma-rays or sham irradiated controls and sacrificed at either 1, 2, 3, 5 or 7 days after exposure. 10 mice were in each experimental group except for 10 Gy on day 7, which only had 8 mice. RNA was profiled using the Agilent-026655 Whole Mouse Genome Microarray. Expression values were log₂ transformed and missing values were set to the median. Mus musculus genes were matched to homologous Homo sapiens genes using biomaRt [37]. If an ensembl ID matched more than one probe, the mean of the values was taken. 15 of the 35 TraMA genes were available. Because the ages of all of the mice were equal, age was arbitrarily set 8 to calculate TraMA. Two of the genes used to indicate cell type were not available, so we only use CD3D, CD19, CD4, CD8A, and CD14 for these analyses.

Measures

Time to death in the HRS was assessed using information about date of interview and date of death from the HRS tracker file. We use 4-year mortality with participants known to be deceased to the HRS. Time to death was calculated as the difference between the 2016 interview month and the known month of death. For participants who survived (i.e., were not known to have died), time at risk was calculated as the time between the 2016 interview month and the most recent interview month available. This measure was used in elastic net models.

Mortality in the HRS for logistic regression. We created a binary indicator of 4-year mortality used in logistic regression models.

Time to death in the LLFS. Because relatively few of the LLFS participants died in 4 years compared to HRS, we used 8-year mortality for the replication analysis using participants known to be deceased by the LLFS. Time to death was calculated as the difference between 2006 blood sample collection date and the known date of death censored on 31 December 2014 (about 8 years after Wave 1 data collection).

Multimorbidity was calculated in HRS as the sum of diseases that a participant had ever been told by a doctor that they had, including high blood pressure, diabetes, cancer, lung disease, heart disease, stroke, and arthritis.

Cognitive dysfunction was assessed in HRS using the Telephone Interview for Cognitive Status (TICS). To make this a measure of dysfunction, we used errors (27 minus the sum of these scores) in immediate recall (10 words), delayed recall (the same 10 words after about 5 minutes of other survey questions), serial 7s (participants were asked to subtract 7 from 100 and continue subtracting for 5 trials), and backwards counting from 20 (participants were asked to count backward from 20 to 10 and given 2 points for a correct first try and 1 for a correct second try).

ADLs and IADLs in HRS are the sum of self-reported difficulty with walking across a room, dressing, bathing, eating, getting in and out of bed, using the toilet, using a map, using the phone, taking medications, managing money, shopping for groceries, and preparing a meal.

Other biological aging measures used in the current study include two epigenetic aging measures produced by HRS [50], PhenoAge [3] and GrimAge [9], and Expanded Biological Age (ExpandedAge) [6]. PhenoAge and GrimAge are both so-called second-generation epigenetic clocks. They are indices of cytosine-phosphate-guanine (CpG) sites where differential methylation is associated with age-related biomarkers and mortality. These have been widely used in past research and have been of extraordinary value in advancing geroscience and in clarifying the biological processes underlying social, psychological, and demographic differences in health and aging [48, 51–53]. We use so-called principal component versions of these clocks which have been shown to be more reliable [54]. ExpandedAge is an index of 22 biomarker of phenotypic aging that has been linked to mortality and other health outcomes [6].

Demographic factors used in regression analyses include chronological age in years, sex/gender (female as the reference group), and race/ethnicity (non-Hispanic Black, Hispanic, non-Hispanic other race, and non-Hispanic White as the reference group).

Socioeconomic factors include years of education as reported to HRS (0-11 years, 12 years (the typical number of years for a high school degree in the US), 13-15 years, and 16 or more years as the reference group), as well as total wealth as calculated by RAND for HRS [55].

Health behaviors include smoker status as reported by respondents (never smoked, past smoker, and current smoker as the reference group), BMI split into five categories (underweight (less than 18.4), overweight (25 to 29.9), obese (30 to 34.9), morbidity obese (greater than or equal to 35), and normal weight (18.5 to 24.9) as the reference group), alcohol use based on number of drinks per day drinking (non-drinker, five or more drinks, and one to four drink as the reference group), and an index of physical activity (the sum of respondent reported light, moderate, and vigorous activity, each ranging from 1 (never) to 5 (every day)).

Analytic plan

Machine learning analyses were conducted in R 4.4.0 “Puppy Cup” [56] using the tidyverse [57], glmnet [58, 59], and lubridate [60] packages. We restricted the set of genes used for training to coding genes with relatively high expression in human venous blood. Of the 50,611 genes that were measured and were successfully mapped in HRS, we restricted ourselves to the 19,291 protein coding genes and to genes with a mean count per million greater than 3 in the total HRS sample, leaving 10,964 genes.

As sensitivity analysis, we also estimated models using the full 50,611 genes in HRS and the 19,291 protein coding genes in HRS and evaluated them in the HRS testing set. These elastic net models had similar C-index in testing data (0.805 and 0.804, respectively), suggesting that the subset of sites captures similar variance. We use the reduced gene set for greatest portability.

We ran elastic net models using Cox regression to predict 4-year mortality hazard using these 10,964 genes, chronological age, and sex in a N = 1802 training set randomly selected from HRS [58, 59]. Mortality hazard was assessed using month of death of participants known by HRS to have died. Hyperparameters, including the alpha and lambda penalty term, were selected using 5x cross validation with a grid search procedure. 11 alpha values were tested with more values close to 0 (viz., 0.000, 0.001, 0.008, 0.027, 0.064, 0.125, 0.216, 0.343, 0.512, 0.729, 1.000). The alpha and lambda values that produced the lowest mean square error were selected. We adjusted the mean and variance to be the same as the mean and variance of HRS age in 2016 to make this an age-like variable.

$$\text{TraMA}=-18.07571+\text{pred}\text{.prob}\times 10.02709(1)$$

Where pred.prob is the log hazard of mortality from the Cox elastic net. We then tested this surrogate score in the N = 1794 testing set. If a gene was missing in a validation cohort, the mean value of that gene from the HRS sample was imputed for all samples (n.b., LLFS had all necessary genes). The R code used to train the TraMA measure, as well as code needed to reproduce this measure in other data is available on Githib at https://github.com/etklopack/TraMA.

Validation in the HRS was conducted in R 4.4.0 “Puppy Cup” [56] using the tidyverse [57] and survey [61] packages. Validation regressions and descriptive statistics in Table 1A, 1B and Figures 1 and 2 used survey weights provided by HRS for use in the RNAseq subsample (vbsi16wgtra). The total testing sample was 1794 participants. Because some participants were missing data on individual outcome variables and covariates, for mortality analyses N = 1791, for multimorbidity analyses N = 1791, for cognitive dysfunction N = 1791, and for ADLs and IADLs analyses N = 1588.

Validation in the LLFS was conducted in R 4.4.0 “Puppy Cup” [56] using coxme and dplyr packages. The total sample size was 1920 participants belonging to visit 1 in LLFS. TraMA was calculated using the algorithm developed in HRS. Mortality was analyzed using mixed-effect Cox proportional hazards regression models with adjustment for family effects.

Validation in public datasets was conducted in R 4.4.0 “Puppy Cup” [56] using the tidyverse [57] and MASS [62] packages. For GSE171730, smoker status was analyzed using logistic regression, and pack years and FEV1 predicted percentage were assessed with linear regression. For GSE123658, diabetes status was assessed using logistic regression. For GSE186507, IBD status and active IBD status were assessed using logistic regression, and severity level was assessed using ordinal logistic regression. For GSE185263, sepsis status and mortality were assessed using logistic regression, and severity was assessed using linear regression.