The gene expression profiles were obtained from the GEO database, including 489 samples and 12,958 genes ( Tables S1 – S3 ). These genes were ranked by the ReliefF algorithm, and then the aging predictor and disease predictor were built using the k-nearest neighbors (kNN; k=3 with the correlation distance) algorithm, optimized by 10-fold cross-validation. The accuracy of the aging predictor in the test set was 0.70455 and 0.7279 in the aging and disease predictors ( Figure 1 ; Table 1 ), respectively. Furthermore, the ROC area under the curve (AUC) for the aging and disease predictor models were 0.7712 and 0.72788 ( Figure 2 ), respectively. As a result, our predictors were sufficiently accurate in both aging and disease models.

Both aging and disease markers have meaningful biological functions. For example, OSBPL1A (oxysterol binding protein-like 1A, ReliefF weight=0.058) was the top aging marker. OSBPL1A is one of a set of intracellular lipid receptors and is closely related to lipid metabolism and cholesterol metabolism (32, 33). TIGD4 (tigger transposable element derived 4, ReliefF weight=0.0253), as the top disease marker, was related to glycogen metabolism. In sum, the abnormal metabolism of lipids, cholesterol and glycogen can lead to T2DM (34). These results indicated the crucial role of energy metabolism in T2DM.

The integrated oxidative aging model was built to explore essential relationships among aging, oxidative and T2DM markers (details are shown in Materials and Methods 5.3, with a total of 11829 “aging-oxidative-disease” triples). The top 10 aging, oxidative and disease markers are shown in Table 2 , including relative experimental details (35–43). For example, ADP-ribosylarginine hydrolase (ADPRH) is the top aging marker, participating in the regulation of various cellular processes, including both immunity and aging (44). ADPRH adversely influences the immune system via CD8+ T cells, hence promoting an imbalance in energy metabolism (45). TPST1 (tyrosyl protein sulfotransferase 1) is the top disease marker, catalyzing the posttranslational sulfation of tyrosine residues within acidic motifs of many polypeptides in all multicellular organisms (46). TPST1 promoted the secretion of some cytokines and then induced the inflammatory response (47). COX5A (cytochrome C oxidase subunit 5A) is the top oxidative marker related to mitochondrial function (48), which induces an imbalance in energy metabolism and insulin resistance (35). In addition, the predictor accuracy calculated by the selected disease markers was 0.7662 ( Table 1 ). In sum, these results indicated that the integrated oxidative aging model could identify essential relationships in T2DM, even with enough prediction ability.

The Markov chain Monte Carlo (MCMC) method was used to evaluate the sensitive relationship between oxidative aging and T2DM. As a result, a series of triples were identified as key components (2501 out of 11829) in the integrated oxidative aging model.

The top 10 sensitive relationships (by calculating the absolute differential frequency) are shown in Table 3 , where the top relationship was “OSBPL7-COX7C-TM6SF1” (difference=-0.03935). Additionally, Table 3 also displayed experimental details of relative oxidative markers (49–55). OSBPL7 (oxysterol binding protein like 7) is an oxysterol-binding protein-like (OSBPL) family member involved in lipid binding and transport and induces cholesterol efflux (56, 57). COX7C (cytochrome C oxidase subunit 7C) is an enzyme in the electron transport chain related to cellular respiration and is also a potential biomarker of diabetes mellitus (58, 59). Transmembrane 6 superfamily member 1 (TM6SF1) participates in regulating transmembrane transport in macrophages (60). Overall, these results indicated that oxidative stress played an important role in the development of T2DM.

The top sensitive aging, disease, oxidative markers (evaluated by the occurrence times, also along with relative experimental details (61–69)) and are also shown in Table 4 . For example, the top aging marker was HPS1 (Hermansky-Pudlak Syndrome 1 gene), inducing the biogenesis of lysosome-associated cellular organelles (70), which regulates the aging process through sphingolipids (71). The top disease marker was PPP1R15A (protein phosphatase 1 regulatory subunit 15A). PPP1R15A plays an important role in insulin resistance via energy metabolism (72, 73). The top oxidative marker was ATOX1 (antioxidant 1 copper chaperone). It has been reported that ATOX1 can regulate the copper level in the cell and maintain the redox balance as a defense antioxidant (74, 75). In short, the sensitivity analysis emphasized the crucial relationship among aging, oxidative stress and T2DM.

To further explore the underlying mechanisms between oxidative aging and T2DM, the shortest path between each pair of oxidative aging and disease markers was identified, and then enrichment analysis was performed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and biological process (BP) terms in Gene Ontology (GO). As a result, relative enrichment results were summarized in Figures 3 , S1 , as well as Tables 5 (75–88), 6 (89–102) and S4 (76–81, 83–85, 103, 104, 111, 113), S5 (89–91, 94, 96, 103, 105–110).

The top 10 KEGG pathways are shown in Tables 5 , S3 . The most enriched KEGG pathway was “Parkinson’s Disease” (enriched in 1213 shortest paths). It has been reported that Parkinson’s disease (PD) and T2DM have common pathological mechanisms (76–78, 111). For example, oxidative stress and mitochondrial dysfunction are involved in both T2DM and PD pathogenesis (77). Strikingly, there are also a series of common biological pathways in T2DM, PD and cancer, such as mitochondrial dysfunction and protein homeostasis (112). Furthermore, the most significant KEGG pathway with the minimum FDR was “Leishmania Infection” (FDR=0.0000629) ( Figure 3B ), indicating the inflammatory response in the immune system (82, 113, 114). Notably, the inflammatory response is also often closely related to cancer (113). The classical aging pathway, the “mTOR signaling pathway” was also enriched in shortest pathway ( Figure S2 ), indicating the interrelationship between oxidative aging and T2DM.

The top 10 BP terms are shown in Tables 6 , S4 . For example, the top enriched BP term was “Regulation of aerobic respiration” (enriched in 29 shortest paths), which was related to energy and mitochondrial function (96). In addition, reactive oxygen species (ROS) are byproducts of aerobic respiration that control various cellular functions (97). The BP term with the minimum FDR was “Regulation of DNA binding” (FDR=0.0000282) ( Figure 3A ), which is vital to T2DM by dysregulating mitochondria and energy metabolism (100). Obviously, the accumulation of DNA damage is also a hallmark of cancer (115). Overall, these results identified various aspects of risk factors for T2DM, such as oxidative stress, aging, energy metabolism and immune systems.

Network markers were identified by calculating the betweenness in the shortest path of each ‘‘oxidative-disease’’ pair, where the top markers are shown in Table 7 . For example, the top network marker was SCD (stearyl-coenzyme A desaturase), which is mainly expressed in adipose tissue and can catalyze the synthesis of monounsaturated fatty acids (116). In addition, SCD can affect lipid metabolism and mediate steroidogenesis, playing an important role in insulin resistance (117, 118). Furthermore, SCD participates in mediating the inflammatory reaction, which promotes the progression of cancer (119). Moreover, there were also a series of shortest paths through SITR1 ( Figure S3 , where permutation p-value=0.002 and 0, before and after sensitive analysis), which was as a clssical aging marker. Thus, network markers indicate the crucial role of oxidative stress dysfunction, along with energy metabolism, in T2DM.

Pan-cancer analysis was used to further verify the relative functions of T2DM oxidative aging markers in cancer. For example, oxidative aging markers in the integrated model were used to evaluate the survival index across different cancer types. There were 9 out of 15 cancer types with significant results (including COAD, ESCA, KIRC, LIHC, LUAD, LUSC, PRAD, THCA and UCEC, shown in Figure 4 ). These results suggest that oxidative aging markers can also be used as relative risk factors in cancer.

Additionally, both the commonality and specificity across 15 cancer types were investigated based on enrichment analysis. The top 10 common KEGG pathways are shown in Figures 5 , S4 , where “Alzheimer’s Disease” was the top KEGG pathway. Alzheimer’s disease (AD) and cancer share common risk factors. For example, aging is one of the greatest risk factors for the development of Alzheimer’s disease, and the risk of cancer also increases with increasing age (120). In addition, some cancer patients may have a higher risk of Alzheimer’s disease (121). Figures 6 , S2 showed the top 10 common BP terms in 15 cancers. “Regulation of cellular respiration” was the top BP term, indicating the key role of energy metabolism in cancer (122). Cellular respiration participates in energy metabolism and is also a hallmark of many cancers (123). The specific enrichment results within each cancer are also summarized in Tables 8 , 9 , S6 , S7 (112, 120–163), indicating a series of oxidative aging-related risk factors in cancer, such as the inflammatory response, energy metabolism and mitochondrial function. Overall, our results highlighted a series of crucial functions related to oxidative aging, which can also be used to study potential mechanisms in cancer.

All gene expression data were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), including GSE362, GSE15790, GSE18732, GSE29221, GSE29226, GSE29231, GSE37171, GSE38642, GSE76894, and GSE182120. These datasets were from eight different platforms: GPL96, GPL97, GPL8450, GPL9486, GPL6947, GPL570, GPL6244, and GPL17586.

The gene expression profiles were processed as follows:

As a result, a total of 489 samples were obtained, including 208 samples of healthy aged people (age > 50 years old, 145 training datasets + 63 test datasets), 131 samples of healthy young people (age ≤ 50, 90 + 41), 110 samples of T2DM aged people (age > 50, 75 + 35) and 40 samples of T2DM young people (age ≤ 50, 25 + 15), containing 12958 gene symbols ( Tables S1 – S3 ).

We also obtained paired gene expression (RNAseq) profiles (“Batch effects normalized mRNA data”) and clinical data from the TCGA database through the xena platform (https://xenabrowser.net/hub/). Cancer types with ≥10 adjacent normal samples were retained. As a result, there were 15 cancer types used in this work: BLCA (408 cancer samples and 19 adjacent normal samples), BRCA(1102 + 113), COAD(451 + 41), ESCA(185 + 11), KICH(66 + 25), KIRC(534 + 72), KIRP(291 + 32), LIHC(376 + 50), LUAD(517 + 59), LUSC(504 + 51), PRAD(498 + 52), READ(160 + 10), STAD(414 + 35), THCA(513 + 59) and UCEC(533 + 22). The tumor expression profiles from the same patient were averaged. Genes with missing values ≥30% were deleted.

After randomization as well as a random disorder, the healthy population samples were divided into a training dataset and a test dataset. The ratio of training dataset samples to test dataset samples was close to 2:1. The ReliefF algorithm was used to select key features, and then the first 500 models were studied to train predictors. The optimal model was selected by 10-fold cross-validation. To verify the accuracy of the aging predictor, the selected model was verified in the test dataset.

The integrated oxidative aging model was built to identify the essential relationship among aging, oxidative stress and T2DM. The computational pipeline was referred to by Mendelian randomization (MR), although it was not as strict as MR ( Figure 1C ).

In this model, the aging-related oxidative stress markers were considered oxidative aging markers, where the relative aging/disease markers were identified in “Methods 5.2”. As a result, the essential relationships among aging, oxidative stress and disease (T2DM) markers were identified as key “aging-oxidative-disease” triples in T2DM.

MR is a statistical method for assessing the causal relationship between risk factors and outcomes based on observational data (181, 182). The causal relationships between the instrumental variables, risk factors, and outcome variables were assessed as follows.

Here, the aging marker was used as the auxiliary variable (similar to the instrumental variable in MR), and the oxidative stress markers were used as the candidate risk factor. Then, aging-related oxidative (“oxidative aging”) markers were identified as the risk factor, and disease markers were used as the outcome variable. That is, the integrated oxidative aging model aimed to explore essential relationships among aging, oxidative stress and disease markers in T2DM. This model was performed as follows:

(1) Oxidative markers were obtained as candidate risk factors based on Biological Processes (BP) of Gene Ontology (GO) through the Gene Set Enrichment Analysis (GSEA) platform (http://www.gsea-msigdb.org/gsea/downloads.jsp, “OXIDATIVE” was taken as the keyword). As a result, 310 candidate oxidative markers were selected.

(2) The correlation (differential coexpression) pattern was used to select aging markers that strongly correlated with candidate oxidative stress markers with the help of the Kruskal−Wallis test. Here, the differential coexpression was calculated as follows:

where the phenotype could be defined as 1 (T2DM) and 0 (control).

Furthermore, both a p-value<0.05 and Benjamini−Hochberg false discovery rate (FDR)<0.1 were used to select strongly correlated aging markers.

(3) To reduce the correlation between the auxiliary variable (aging marker) and confounding factors, as well as further select a strong correlation between the aging marker and the candidate oxidative marker, a permutation test was performed by generating the simulated aging markers from the same number of randomly selected markers to each candidate oxidative marker; this process was repeated 1000 times, and then the p-value was calculated as the proportion of occurrence times (larger than the real mean difference) of the absolute difference between T2DM and control in 1000 permutations. The relationship between each aging marker and the candidate oxidative marker was retained if the permutation P<0.05.

(4) Correlation (differential coexpression) was used to select oxidative markers that strongly correlated with disease markers with the help of the Kruskal−Wallis test. Here, the differential coexpression was calculated as follows:

where the phenotype could be defined as 1 (T2DM) and 0 (control).

Furthermore, both a p-value<0.05 and Benjamini−Hochberg false discovery rate (FDR)<0.1 were used to select strongly correlated oxidative markers.

(5) To reduce the correlation between the risk factor (oxidative marker) and confounding factors, as well as further select a strong correlation between the oxidative marker and the disease marker, a permutation test was performed by generating the simulated oxidative markers from the same number of randomly selected markers to each disease marker; this process was repeated 1000 times, and then the p-value was calculated as the proportion of occurrence times (larger than the real mean difference) of the absolute difference between T2DM and control in 1000 permutations. The relationship between each aging marker and the candidate oxidative marker was retained if the permutation P<0.05.

(6) The direct relationships for any other factors (genes) were found to reduce the correlation between the auxiliary variable (aging marker) and confounding factors. If there was another factor (gene) that was directly correlated (differentially coexpressed) to both the aging marker and the disease marker, then the relationship from aging to disease was deleted.

where the phenotype could be defined as 1 (T2DM) and 0 (control).

Furthermore, both a p-value<0.05 and Benjamini−Hochberg false discovery rate (FDR)<0.1 were used to filter out any direct relationships.

(7) To filter out the effect of horizontal pleiotropy, the aging–disease relationship was further examined by comparing the correlation between each aging and disease marker, through the oxidative marker or otherwise. Herein, steps ①–③ were used to calculate the correlations between auxiliary variables and outcome variables without the background of the risk factor, and step ④ was used to calculate the correlations between auxiliary variables and outcome variables with the context of the risk factor.

① The residual of each disease marker (“residual A”) was calculated based on the oxidative marker:

where b₁ is the regression coefficient.

② The residual of each aging marker (“residual B”) was calculated based on the oxidative marker:

where b₂ is the regression coefficient.

③ The abovementioned two residuals were further compared, and the residual of the disease marker was calculated (as “residual C”):

where b₃ is the regression coefficient.

④ The residual of the disease marker (“residual D”) was calculated based on the aging marker.

⑤ The difference (between “residual C” and “residual D”) was tested between the T2DM and control subgroups using the Kruskal–Wallis test (P<0.05 and FDR<0.1).

Finally, the essential relationship among the aging marker, oxidative marker and disease marker was retained. Thus, 11829 “aging-oxidative-disease” triples were identified, including 105 aging markers, 83 oxidative markers and 282 disease markers. Thus, these 83 oxidative markers were used as oxidative aging markers (risk factors), and 282 disease markers were also used to discriminate the T2DM phenotype.

To further explore the relationship among aging, oxidative stress and T2DM, sensitivity analysis was performed based on the Markov chain Monte Carlo (MCMC) method, where “aging-oxidative-disease” triples identified by MR were further evaluated as a candidate relationship. The MCMC method is used to sample certain posterior distributions in a high-dimensional space based on a given probabilistic background. The key step of MCMC is to construct a Markov chain whose equilibrium distribution is equal to the target probability distribution. The steps were as follows:

(1) Constructing the transfer cores of the ergodic Markov chain. The prior distribution of each parameter was normally distributed based on all identified markers in each group (i.e., T2DM and control), respectively.

(2) Simulate the chains until equilibrium is reached. The Metropolis−Hastings sampling method was used to determine whether the new sample (θ *) was acceptable based on the α value.

where P (θ ⁿ | X) and P (θ * | X) are the posterior probability of the nth accepted sample, the new sample q (θ ⁿ → θ *) is the transition probability from the nth accepted sample to the new sample, and q (θ * → θ ⁿ ) is the transition probability from the new sample to the n-th accepted sample.

In this work, the disease score was used to evaluate the simulated samples, with 1000 random samples used as candidate samples for each group (i.e., T2DM or control). The disease score was calculated by comparing the distance between normal and T2DM training samples based on the 282 disease markers identified by the integrated oxidative aging model:

(3) Performing the global sensitivity analysis

The correlation index was used to evaluate each “aging-oxidative-disease” triple in the accepted samples (including both T2DM and control):

As a result, the correlation index was calculated in each “aging-oxidative-disease” triple for all accepted samples. Then, the Kruskal–Wallis test was used to evaluate each correlation index in each “aging-oxidative-disease” triple, where p-value<0.05 and FDR<0.1 were set as the threshold. Finally, 2501 “aging-oxidative-disease” triples were identified as sensitive relationships, including 41 aging markers, 37 oxidative markers and 61 disease markers.

To further reveal the relationship between “oxidative aging” and T2DM, a differential coexpression network was constructed by the following steps:

The gene functions were further explored by enrichment analysis of the shortest pathway. Gene Ontology (GO) terms and KEGG pathways for the GSEA platform were obtained from gene set enrichment analysis (http://software.broadinstitute.org/gsea/downloads.jsp, version 7.5). The hypergeometric distribution was used to test the degree of enrichment of the GO BP and KEGG pathways. Hypergeometric test formula:

where N is the total number of genes in the gene set, M is the number of known genes (such as KEGG pathway or BP terms), which is the number of genes identified in each shortest pathway, and k is the number of common genes between known genes and candidate genes identified in each “oxidative-disease” shortest pathway. The p-value of each path was controlled using the Benjamin-Hochberg method. Finally, pathways with p<0.05 and FDR<0.1 were retained.

The subnetwork with the shortest pathways among the selected “oxidative-disease” pairs was constructed, and genes in the subnetwork were sorted by their betweennesses in descending order. To test whether the top betweenness genes were hubs in the background network, we ran a permutation to count the occurrence time of the top genes in the shortest paths between randomly selected genes (containing the same numbers of “oxidative-disease” pairs, based on the identified “aging-oxidative-disease” triples) when they had greater betweennesses than those in our study. We repeated this process 1000 times, and the p-value was calculated as the proportion of occurrence times of the top betweenness genes in 1000 permutations.

The survival analysis was performed based on the oxidative aging markers (identified by the integrated oxidative aging model in 5.3) for each cancer using the Kaplan−Meier method. The tumor samples of each cancer were divided into two groups based on the mean value of the oxidative aging markers. Then, the Kaplan−Meier method was used to evaluate the survival difference between these two groups, and the significance was estimated by the log-rank test. A p-value<0.05 was considered statistically significant.

Genes were considered differentially expressed if they satisfied the following criteria:

Then, the differential expression networks were constructed for each cancer, where the details were also the same as 5.5. As a result, each shoreat pathway was selected from each pair of oxidative aging markers and differentially expressed genes (as disease markers in cancer) using the Dijkstra algorithm. Furthermore, enrichment analysis was performed by the “oxidative-disease” shortest pathway for each cancer type, where both p<0.05 and FDR<0.1 were used.