We sourced senescence-related genes from the SenMayo dataset and categorized them into three groups based on Dominik Saul et al.’s classification: transmembrane, secreted, and intracellular (23). We validated the expression of these genes in bone biopsies from aged human cohorts. A total of 125 genes constitute the senescence gene set, with 102 genes included in the training set for signature modeling.

In this study, we obtained pan-cancer transcriptomic data from TCGA (n = 8,739). The TCGA data were downloaded from USCS Xena (http://xenabroswer.net/hub). We randomly divided the dataset into a training set and a test set at a ratio of 7:3.

The workflow was as follows: (1) We performed univariate analysis based on the training set. The prognostic genes (p < 0.0001) remained for the following analysis ( Supplementary Table S1 ). (2) We used the LASSO algorithm for further screening of prognostic genes. (3) We used the Cox proportional hazard regression model with stepwise for the final model fitting. We calculated the senescence score for each patient as follows: senescence score = coefficient *expression of genes (SenMayo: 0.026*WNT16 + 0.037*VGF+-0.128*TUBGCP2+-0.068*TNFRSF10C+0.05*SPP1+-0.033*SERPINE2 + 0.135*SERPINE1+-0.107*SCAMP4 + 0.106*PLAUR+-0.041*PLAU+0.221*PIGF+-0.121*NAP1L4+-0.034*MMP9+-0.044*MMP2 + 0.182*MMP14+-0.138*LCP1 + 0.105*ITPKA+-0.063*ITGA2+-0.038*IQGAP2 + 0.033*IL7+-0.045*IL1B+0.03*IL1A+0.034*IGFBP3 + 0.049*IGFBP2+-0.084*ICAM3 + 0.044*ICAM1 + 0.052*HGF+-0.085*ETS2 + 0.067*EREG+0.037*EGFR+-0.035*EDN1 + 0.052*CXCR2 + 0.054*CXCL3+-0.056*CXCL2 + 0.034*CXCL10 + 0.041*CSF1+-0.038*CCL2+-0.037*CCL13 + 0.066*C3 + 0.023*BMP2+-0.038*AREG+-0.027*ANGPTL4. FRIDMAN.SENESCENCE.UP, a gene set of senescence-related genes, originated from Fridman [29]: -0.07*CD44 + 0.17*CDC25B+-0.056*CITED2 + 0.22*CLTB+ -0.054*COL1A2 + 0.13*CCNB1 + 0.198*E2F4 + 0.111*EIF2B1 + 0.088*FN1 + 0.095*HPS5 + 0.081*HTATIP2 + 0.03*F3 + 0.041*IGFBP2 + 0.107*IGFBP3 + 0.076*WISP1+-0.055*IGF1R+0.064*MORF4 + 0.059*MYC+-0.051*NDN+0.092*OPTN+-0.044*CDKN2B+-0.139*CDKN2D+0.082*MAPK14+-0.151*NOLC1 + 0.201*PEA15 + 0.055*RAB31 + 0.111*RAB5B+-0.313*RABGGTA+0.256*RAC1+-0.119*S100A11 + 0.024*SERPINB2+-0.085*STAT1 + 0.039*TERT+0.043*TFAP2A+0.033*THBS1 + 0.049*TNFAIP2+-0.065*TNFAIP3. Tumor samples were classified into high-risk and low-risk groups according to the median senescence score.

The z-score, as introduced by Lee et al. (24), is employed to estimate the activity of specific pathways by amalgamating the expressions of feature genes. We analyzed gene sets associated with glycolysis (sourced from https://www.gsea-msigdb.org/gsea/msigdb/human/geneset/WP_GLYCOLYSIS_IN_SENESCENCE), autophagy (25), and epithelial-to-mesenchymal transition (EMT, sourced from Tian-jian Yu et al.) (26), and we analyzed cell cycle (27) using the GSVA implementation in the R package. The resulting values for each gene set were designated as glycolysis z-score, autophagy z-score, EMT z-score, and cell cycle z-score, respectively.

We defined differentially expressed genes (DEGs) as those with a log2 fold change (|log2FC|) of 1 or greater and a false discovery rate (FDR) p-value of less than 0.05. We used the limma R package for this analysis.

We performed all statistical analysis using R version 4.2.2. We assessed patients’ survival outcomes in the low-risk and high-risk groups using Kaplan–Meier analysis. To assess the clinical utility of each model, we employed decision curve analysis (DCA), a method originally proposed by Vickers and Elkin (28). We evaluated prognostic performance using time-dependent receiver operating characteristic (ROC) curve analysis. We compared the areas under the curves (AUCs) of the senescence score and the nomogram model. For all analyses, we considered a p-value of less than 0.05 to be statistically significant.

The workflow for constructing the senescence-related signature of pan-cancer is shown in Figure 1 . First, we performed a univariate analysis on 102 senescence-related biomarkers, retaining 57 genes with a p-value of less than 0.0001 ( Supplementary Table S1 ). Next, we used these 57 biomarkers for LASSO regression analysis in the TCGA pan-cancer training set ( Figures 2A, B ). Finally, we applied the Cox proportional hazard regression model with stepwise, resulting in the selection of 42 senescence-related biomarkers. We derived the senescence score based on the 42 senescence-related biomarkers’ normalized expression levels.

We classified senescence-related genes into transmembrane, secreted, and intracellular categories based on their functional characteristics (23). To further understand the differences in 42 senescence-related biomarkers between tumor and normal tissues, we investigated their expression levels in pan-cancer versus normal tissues. We discovered that secreted biomarkers (SPP1, SERPINE1, PLAU, and MMP9) are highly expressed in most tumor tissues, suggesting the role of the SASP in driving tumor progression ( Figure 2C ). To further understand the impact of senescence-related genes on tumor prognosis, we performed univariate Cox analysis ( Figure 2D ). Our findings indicate that most senescence-related genes have a detrimental effect on tumor prognosis.

Significant differences in senescence scores are evident across various organs. Tumors originating from the gastrointestinal tract, such as esophageal cancer, stomach cancer, and colorectal cancer, exhibit higher senescence scores. In contrast, certain gland-associated tumors, including thymoma (THYM), prostate adenocarcinoma (PRAD), and thyroid carcinoma (THCA), display lower senescence scores ( Figure 3A ). The glioblastoma multiforme (GBM), malignant mesothelioma (MESO), and pancreatic adenocarcinoma (PAAD) exhibit the highest senescence scores. Moreover, the senescence score was a risky factor in many cancer types ( Figure 3A ). However, the senescence score was a protective factor in DLBC.

Within the TCGA training cohort, patients were classified into high-score or low-score groups based on senescence scores. Patients with high senescence scores across various cancer types experienced multiple adverse survival outcomes, including short disease-specific survival (DSS), overall survival (OS), and progression-free interval (PFI) ( Figure 3B ). We observed similar results in the test cohort, where the high-score group exhibited worse survival outcomes ( Figure 3C ). This suggests that our senescence score may have a universal applicability in pan-cancer.

We then analyzed the relationship between senescence scores and different pathological stages in pan-cancer. We examined 21 tumor types and observed a significant increase in senescence scores in 14 of them as tumor stages advanced. They are adrenocortical carcinoma (ACC), bladder cancer (BLCA), breast cancer (BRCA), colon adenocarcinoma (COAD), head and neck cancer (HNSC), kidney chromophobe (KICH), kidney clear cell carcinoma (KIRC), kidney papillary cell carcinoma (KIRP), liver cancer (LIHC), lung adenocarcinoma (LUAD), PAAD, stomach cancer (STAD), testicular germ cell tumor (TGCT), and THCA ( Figure 4A ). With increasing tumor grades, we also observed a significant increase of senescence score in uterine corpus endometrial carcinoma (UCEC) and lower-grade glioma (LGG) ( Figure 4B ). These findings suggest that our senescence score not only is applicable for assessing advanced pan-cancer but also holds significance for identifying early-stage tumors. To validate our findings, we analyzed another senescence signature, termed FRIDMAN.SENESCENCE.UP (29). Significant differences in senescence scores from different tumor types were also observed. The GBM and MESO also exhibit the highest senescence scores ( Figure 5A ). Compared with gland-associated tumors, like PRAD and pheochromocytoma and paraganglioma (PCPG), gastrointestinal tract tumors exhibit higher senescence scores. Furthermore, the validation dataset also supports the idea that the senescence score was a risky factor in many cancer types, consistent with senescence scores in the SenMayo dataset ( Figure 3A ). In the validation dataset, we classified patients into high- or low-score groups based on their senescence scores. Regardless of whether they were in the training cohort or the test cohort, patients with high scores had worse survival outcomes, including shorter DSS, OS, and PFI ( Figures 5B , C ). In that way, this validation indicates that our senescence score could be applicable across various tumor types.

To enhance the practical application of senescence scores in clinical practice, we have developed a comprehensive nomogram ( Figure 6A ). The nomogram integrates senescence scores with clinical characteristics like age and cancer type. In that way, the prediction of patient survival can be easily conducted. The calibration curves for 3-year and 5-year OS ( Figure 6B ) closely align with the standard curve, indicating that the nomogram accurately predicts actual survival probabilities. This alignment indicates that the nomogram accurately predicts the actual survival probabilities. Furthermore, the AUC values for the nomogram are 0.76 and 0.79 in the training and test sets, respectively, underscoring their efficacy in predicting survival rates ( Figure 6C ). Furthermore, DCA in both training and test cohorts shows that the net benefit of combining individual senescence scores with the nomogram is consistently above zero ( Figure 6D ). Although DCA suggests that the nomogram may not offer a significant advantage over senescence scores alone, its utility in clinical practice lies in its ability to help identify high-risk patients with lower survival probabilities, based on tumor type and patient age. For example, as indicated by the nomogram, patients diagnosed with GBM tend to have a lower survival rate compared to those with testicular germ cell tumors (TCGT). This finding underscores the significance of integrating senescence scores with clinical characteristics to enhance the accuracy of prognostic predictions.

Several biological activities occur during the initiation of malignant tumors. As normal cells transform into a malignant state, they require active EMT, rapid and unlimited proliferation, and increased glycolysis frequency (30). It has been reported that senescent cells exhibit high glycolysis rates, which align with the increased demand for proteins and lipids (31). Conversely, autophagy partly suppresses cellular senescence (32). To explore the relationship between tumor characteristics and the senescence signature, we employed the z-score algorithm to quantify the activity of cancer-related pathways in various tumor types ( Figure 7 ). We found that the senescence score was positively correlated with glycolysis (R = 0.28, p < 0.001), EMT (R = 0.35, p < 0.001), and cell cycle (R = 0.44, p < 0.001). The senescence score was negatively correlated with autophagy (R = 0.28, p < 0.001). In other words, tumors with a strong senescence state exhibit more active glycolysis, EMT, and cell cycle activities, which represent an aggressive phenotype. Meanwhile, these tumors usually show reduced autophagy, which might promote tumor progression. Furthermore, we evaluated the activity of glycolysis, EMT, cell cycle, and autophagy across a wide range of tumor types. In most tumors, the EMT exhibited strong positive correlations with senescence scores ( Figure 7E ). The correlations between glycolysis and senescence scores across different tumor types are shown in Supplementary Figure S1 , while those between the cell cycle and senescence scores are depicted in Supplementary Figure S2 . Most kinds of tumors exhibit positive correlations with senescence scores. Conversely, z-scores of autophagy were mostly negatively correlated with senescence scores ( Supplementary Figure S1B ). In summary, our findings indicate that tumors with high senescence scores exhibit more aggressive cellular characteristics.

Based on our analysis, we know that tumors with high senescence scores exhibit more aggressive characteristics. It is of great significance whether the biological functions involved in the high senescence score group are related to the mechanism of tumor invasion. To investigate this, we conducted functional analyses to discern differences between high and low senescence score groups. This analysis included biological processes (BP), cellular components (CC), molecular functions (MF), and pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG).

In the group with high senescence scores, there was a notable enrichment in molecular functions such as tubulin binding, microtubule binding, cytokine activity, extracellular matrix structural constituent, and growth factor binding. Biological processes like organelle fission, nuclear division, chromosome segregation, mitotic nuclear division, and mitotic sister chromatid segregation were also enriched. Upregulated cellular components included the spindle, chromosome centromeric region, condensed chromosome, and kinetochore. Additionally, KEGG signaling pathways showed significant enrichment in cytokine–cytokine receptor interaction, cell cycle, IL-17 signaling, amoebiasis, and pertussis ( Figure 8A ).

In the high-risk group, we observed a downregulation of genes associated with pathways such as adrenergic signaling, dopaminergic synapse, cholinergic synapse, insulin secretion, and thyroid hormone synthesis. Molecular functions like ATPase-coupled transmembrane transporter activity, ligand-activated transcription factor activity, nuclear receptor activity, transcription coactivator binding, and dinucleotide phosphatase activity were also downregulated. Additionally, biological processes and cellular components related to the metabolic processing of phenol-containing compounds, mammary gland epithelium development, thyroid hormone metabolism, mammary gland duct morphogenesis, thyroid hormone generation, and the apical plasma membrane showed decreased activity in the high-risk group ( Figure 8B ). These findings contribute new insights for future research into related biological mechanisms.