We downloaded the gene expression files of 370 HCC and 50 adjacent normal tissues from the TCGA database (https://portal.gdc.cancer.gov/) using the “TCGAbiolinks” package (Colaprico et al., 2016). The corresponding clinicopathologic characteristics and their survival information including overall survival (OS), disease-free survival (DFS), progression-free survival (PFS), and disease-specific survival (DSS) were also gathered. The one-class logistic regression (OCLR)-based mRNAsi and mDNAsi for evaluating the differentiation degree from bulk sequencing were obtained as previously described (Malta et al., 2018). For survival-related analysis, 336 of 370 patients with the survival time of >1 month were included (Cai et al., 2020; Guo et al., 2022), which were divided into the training and validation datasets with the ratio of 2:1 (222 vs. 114). The transcriptome profiles and survival data of 238 HCC patients from the International Cancer Genome Consortium (ICGC) database (ICGC-LIRI-JP) were acquired to serve as the external validation set. All data were normalized and processed as previously reported (Guo et al., 2022; Tang et al., 2022). 1,259 cellular senescence-related genes (CSGs) including 525 positive CSGs and 734 negative CSGs were retrieved from a recently published paper (Wang et al., 2022). All the patient data used for model construction in the current study were summarized in Supplementary Table S1.

Differentially expressed genes (DEGs) between HCC and adjacent normal tissues were screened by using the “limma” package (Ritchie et al., 2015), with the criteria cutoff of adj.P.Val<0.01 and | logFC|>1. To balance the samples between the tumor and normal groups, we utilized the subset-based method as we previously described (Guo et al., 2021). Generally, five subsets were randomly generated from the tumor group with the sample size of 74, yielding a ratio of 1.48:1 for tumor and normal. The common DEGs in the five subset-based analyses were obtained by plotting a Venn diagram, which were further intersected with 1,259 CSGs (Wang et al., 2022) to acquire the DECSGs in HCC. Subsequently, Kaplan-Meier survival analysis and log-rank test were used to identify the prognostic DECSGs for OS. To enhance the stability, we adopted the “multi-split” approach as we reported with some modifications (Tang et al., 2020; Guo et al., 2022; Zhang et al., 2022). Next, the prognostic DECSGs were utilized to perform unsupervised clustering for the patients with complete survival information, and the R package “ConsensusClusterPlus” was used to explore novel CSG patterns of HCC (Wilkerson and Hayes, 2010; Wang et al., 2021; Guo et al., 2022).

The package “limma” was utilized to screen the DEGs between two CSG subtypes of HCC, followed by the enrichment analysis of biological functions and pathways. The “clusterProfiler” R package was applied to conduct the Gene Ontology (GO) enrichment and the most significant terms were visualized by using the “ggplot2” package. Meanwhile, the online DAVID tool (https://david.ncifcrf.gov/tools.jsp) was used to perform the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, and a network showing the most significant enriched pathways was drawn by the package “ggraph.”

The Gene set variation analysis process was completed to calculate the enrichment scores of 50 hallmark gene sets of the molecular signature database (Liberzon et al., 2015) for the TCGA-LIHC cohort by the “GSVA” package (Hanzelmann et al., 2013). Differential analysis was then performed by the “limma” package to reveal the relative activities of these hallmark pathways in the two CSG subtypes, which was determined by the criterion of p < 0.01 combined with absolute t > 2.5 (Lambrechts et al., 2018; Zhang et al., 2022). GSEA was further employed to verify the GSVA results using the GSEA software with default settings. The cutoff of FDR <0.25 was considered significant as recommended (Zhang et al., 2022).

In this section, all 336 samples from the TCGA-LIHC cohort were randomly classified into the training set (n = 222) and the validation set (n = 114) at a ratio of 2:1, and the ICGC-LIRI-JP dataset (n = 238) was selected as an external validation set. In the training set, two feature selection algorithms suitable for high dimensionality were applied to screen the most subtype-relevant features: LASSO regression and the Random Forest. The expression of the top 208 DEGs between two CSG subtypes with the cutoff of |logFC| > 1.5 and adj.P.Val<0.01 was used as the input variable and the CSG subtype status was used as the binary outcome (0 or 1). For the LASSO selection operator, the “multi-split” strategy was adopted to minimize the effect of arbitrary choice in the random sample split process (Xu et al., 2017; Luo et al., 2020). We subsampled 75% of the training dataset 1,000 times and only those repeatedly occurred more than 700 times were selected. Moreover, random forest was carried out for variable elimination with the OOB error as the minimization criterion. The intersecting genes were further narrowed down by comparing them with the hub genes in HCC that were reported before (Tang et al., 2020; Zhang et al., 2021). Afterward, two overlapping critical genes were used to fit a binary prediction model with multivariate logistic regression analysis (Xu et al., 2017; Luo et al., 2020; Wang et al., 2021). The discriminative performance of the CSG subtype predictor was investigated by the receiver operating characteristic (ROC) curves and the areas under the ROC curve (AUCs) in the training and validation datasets. The optimal cutoff value derived from the training set according to the maximum Youden’s index was applied to the training and validation datasets to generate the predicted outcomes and confusion matrixes. Finally, the CSG subtypes for the external dataset ICGC-LIHC-JP were predicted in a similar way and the prognostic value of the CSG subtype was further validated.

We examined the utility of the CSGs for the construction of a prognostic CSGs scoring system with the aforementioned prognostic DECSGs. Based on the TCGA training dataset, all the 238 robust prognostic DECSGs were subjected to LASSO-COX and stepwise regression algorithms to generate the best subset of prognostic genes, which were subsequently used to build the prognostic CSGs signature: CSGscore = Σ(coef (β)*EXP_β), where β indicates each selected prognostic DECSG. Next, the CSGscores for the HCC patients in the training, validation, and external dataset were computed, and all the patients were categorized into the high- and low-risk groups depending on the median CSGscores in the training set. Kaplan-Meier survival curves and log-rank tests were used to compare the divergent survival outcomes in both groups. Time-dependent receiver operating characteristic (tROC) curves were drawn to evaluate the predictive power of CSGscore. Besides, stratified survival analysis was used to explore its additional prognostic value in subgroups divided by clinicopathologic variables.

Univariate and multivariate analyses were conducted to verify the independent prognostic capacity of CSGscore. Based on the results of Univariate analysis, a CSGscore-integrated nomogram was established with the “rms” package. Calibration curves were plotted to estimate the consistency between predicted probabilities of 3- and 5-year survival and actual ones. The decision curve analysis (DCA) curves at 1-, 3-, and 5- years were depicted to compare the net benefits of the nomogram and that of the pathologic stage or tumor burden status. Kaplan-Meier analysis was further applied to test the clinical relevance of the nomogram for OS, DFS, PFS, and DSS.

Based on the whole TCGA-LIHC cohort, the potential role of CSGscore in predicting the immunotherapeutic efficacy of HCC was preliminarily assessed by comparing the gene expression of 50 immune checkpoint genes (ICGs) between distinct CSGscore groups. Immunophenoscore (IPS) representing the tumor immunogenicity of each sample was acquired from The Cancer Immunome Atlas (TCIA) (https://tcia.at/home), and a lower IPS stands for worse sensitivities to immunotherapy (Charoentong et al., 2017). The tumor immune dysfunction and exclusion (TIDE) score, which was designed to predict the influences on survival and immunotherapeutic responses based on two mechanisms of immune evasion: T-cell exclusion and T-cell dysfunction (Jiang et al., 2018), was computed with the expression profiles of the whole TCGA-LIHC cohort and the ICGC-LIRI-JP dataset. Wilcoxon test was used to compare the TIDE scores between the high- and low- CSGscore groups, and Chi-square test was used to analyze the ratio of responders or non-responders in both groups. Meanwhile, the “pRRophetic” package was used to calculate the semi-inhibitory concentration (IC₅₀) values of 138 chemo/targeted drugs (Geeleher et al., 2014).

For the immune infiltration landscape, we utilized the CIBERSORT algorithm (Newman et al., 2015) to estimate the relative percentage of 22 tumor immune cell types based on the whole TCGA-LIHC cohort. The tumor immune cell infiltration for the high- and low-risk groups were computed and presented with a heatmap. Spearman correlation analysis was used to investigate the correlation of immune cell infiltration and CSGscore as well as the expression of CSGscore genes, followed by the visualization of a correlation heatmap.

The innovative OCLR machine-learning algorithm-derived stemness indices including mRNAsi and mDNAsi have correlated with multiple clinical observations in malignant tumors such as poor survival, tumor metastasis, and therapeutic resistance (Malta et al., 2018; Wang et al., 2021; Huang et al., 2022). Thus, we used the whole TCGA-LIHC cohort to compare the mRNAsi and mDNAsi in the two CSGscore risk groups by Wilcoxon test. Spearman correlation analysis was then performed to validate the relationship between each CSGscore gene and mRNAsi or mDNAsi. The differential expression patterns of these CSGscore genes were further verified by paired HCC and non-cancerous tissue samples (n = 50). Finally, the correlation of CSGscore genes’ expression and clinicopathologic parameters including tumor grade, pathologic stage, tumor burden status, and CSG subtype was also investigated via boxplots.

Kaplan-Meier and log-rank tests were used to compare the survival outcomes for OS, DFS, PFS, and DSS in different groups. Wilcoxon test or Kruskal-Wallis test was used for the comparisons of continuous variables in two or three groups. Correlation coefficients and statistical significance were calculated by Spearman correlation analysis. Pearson Chi-square test was utilized to examine the distribution differences of categorical variables. The R package “glmnet” and “randomForest” were used to perform the LASSO and random forest algorithms for dimensionality reduction. The “pROC” package was used to assess the predictive accuracies for the CSG subtype predictor and the “timeROC” package was used to determine the prognostic reliability of the CSGs scoring system. All data visualization was completed in RStudio (Version 1.1.383, https://www.rstudio.com/).

Applying the “limma” package with the screening criteria of |logFC| > 1 and adj.P.Val <0.01, we extracted 1997 overlapping DEGs among five randomly subsampled subsets from TCGA-LIHC dataset using the subset-based approach (Figures 2A,B). Next, we intersected the previously reported 1,259 cellular senescence genes (CSGs) with these 1997 DEGs to obtain 331 DEGs associated with CS (DECSGs) (Figure 2C). Currently, a growing body of evidence has linked CS with cancer progression and metastasis, which in turn affects the prognosis (Wang et al., 2022). Therefore, we screened the robust DECSGs significantly correlated to the prognosis of HCC patients by the “multi-split” method, and 238 DECSGs were observed to be associated with prognosis in more than 990 of 1,000 subsamples (Supplementary Table S2), and by those genes, two subtypes of HCC were recognized by unsupervised clustering analysis (Figures 2D–F, Supplementary Figure S1). Furthermore, by utilizing Kaplan-Meier analysis, we found that patients in the two subtypes had distinct clinical outcomes, and patients in subtype 2 had a significantly shorter survival time than the patients in subtype 1 (Figure 2G). We also exhibited the gene profiles of the 238 DECSGs in the two HCC subtypes along with clinical traits (Figure 2H).

Differential analysis of clinicopathologic characteristics between the two HCC subtypes from the TCGA-LIHC cohort was performed. Notably, Chi-square test results showed significantly differential sample distribution for pathologic stage, tumor grade, OS status, and progression status of HCC concerning the defined CSG subtype. (Figure 3A). More patients with advanced stage (35.5% vs. 18.8%), high grade (50.7% vs. 25.9%), and poor progression status (61.1% vs. 44.9%) were clustered into subtype 2 than subtype 1, corresponding to the poorer prognosis of subtype 2.

GO enrichment analysis revealed that the most significant terms enriched by the DEGs between the two subtypes were the biological process (BP) of DNA replication, cellular component (CC) of chromosomal region, and molecular function (MF) of catalytic activity acting on DNA (Figure 3B). For the KEGG enrichment analysis, DEGs between the two CSG subtypes mostly participated in the pathways of metabolic pathways and cell cycle (Figure 3C). GSVA regarding hallmark gene sets was performed in the two HCC subtypes. As Figure 3D shows, the CSG subtype 2 was significantly associated with several cell cycle-related pathways such as G2M checkpoint, E2F targets, and mitotic spindle pathway while those from subtype 1 were associated with immune- and metabolism-related pathways such as coagulation, xenobiotic metabolism, and fatty acid metabolism pathway. Meanwhile, GSEA confirmed the differences in hallmark pathways enriched in the two HCC subtypes. Similar to the results of GSVA analysis, we observed that multiple cell cycle- and cancer-related pathways were more frequently enriched in subtype 2 (Figure 3E).

Preliminary work inspired us to investigate potential biomarkers with indicative roles in the classification of CSGs-related HCC subtypes to achieve a precise and concise model for CSG subtype categorization of HCC. We first picked out the top 208 of 784 DEGs between the two HCC subtypes with the cutoff of |logFC| > 1.5 (Supplementary Table S3). Next, we conducted LASSO regression and random forest with the 208 DEGs and screened out nine and 14 variables, respectively. Five genes (GINS1, TOP2A, KIF11, KIF2C, and MELK) were found in common by both algorithms. Interestingly, TOP2A and KIF11 were widely recognized as hub genes of HCC, which were also presented to be strongly associated with the CSGs-related subtype classification. Figure 4A outlined the screening process of these two critical hub genes (TOP2A and KIF11), which were further used to construct the CSGs-related subtype predictor. Confusion matrixes showed that the expression of TOP2A and KIF11 in HCC patients could precisely predict the CSGs-related subtype in the TCGA training cohort (totally 217 of 222 samples were correctly predicted) and the TCGA validation cohort (totally 109 of 114 samples were correctly predicted) and ROC curves also demonstrated the reliability of the two-gene subtype predictor (Figures 4B,C). Moreover, the confusion matrix and ROC curve for the whole TCGA-LIHC cohort were shown in Supplementary Figure S2. Additionally, we performed CSG subtype prediction with the subtype predictor for the ICGC-LIRI-JP dataset and 238 patients were successfully classified into two groups with distinct OS status (Figures 4D,E).

To evaluate the capability of CSGs in the prognosis prediction of HCC, the abovementioned 238 robust prognostic DECSGs were applied as the input variables for the LASSO-COX regression analysis to construct a prognostic model, which resulted in 10 potential genes for the next incorporation into the stepwise regression (Figure 5A). The final model consisting of PSRC1, SOCS2, TMEM45A, CCT5, and STC2 generated a minimal AIC value (AIC = 709.63). Therefore, the final prognostic CSGscore risk model was constructed with the five prognostic CSGs, whose coefficients and significance levels in the model were shown in Figure 5B. To evaluate the prognostic value of the CSGscore model for HCC, all patients in the TCGA training set were classified into high- and low-risk groups according to the median value of CSGscore (0.9824545), and the CSGscore distribution, survival status, and expression profile of HCC patients were exhibited Figure 5C. As shown in Figure 5D, Kaplan-Meier analysis showed that high-risk patients had exceedingly lower OS rates relative to low-risk patients in the TCGA training cohort. The time-dependent receiver operating characteristic curve showed the area under the ROC curve (AUC) of 0.871, 0.801, and 0.784 in 1-year, 3-year, and 5-year OS, respectively, suggesting the high accuracy of CSGscore for OS prediction (Figure 5E). Moreover, we also found the high-risk group had a worse prognosis for DFS, PFS, and DSS (Figures 5F–H). Notably, stratified analysis indicates the CSGscore signature could recognize different risk status significantly in all subgroups divided by clinicopathologic parameters, which denotes its good potential and clinical application to exert additional prognostic value to existing risk factors (Supplementary Figure S3).

To validate the stability and robustness of CSGscore for OS prediction, we adopted the TCGA-validation dataset and the whole TCGA cohort as internal validation sets and chose the ICGC-LIRI-JP dataset as an external set. Using the same cutoff value that was derived from the training set (0.9824545), all patients in each validation set were categorized into high- and low-risk groups. As Figure 6 shown, the CSGscore was consistent in predicting the prognosis of HCC in all of them, and CSGscore-defined high-risk patients had a significantly worse prognosis than the low-risk group.

The correlation between CSGscore and clinicopathologic features was subsequently explored. As presented in Figure 7A and Supplementary Figure S4, the CSG subtype, tumor burden, stage, and tumor grade were significantly related to CSGscore, while other clinical characteristics showed no significant association with CSGscore.

The accurate preoperative prediction of OS might help surgeons to make better clinical decisions and comprehensive nomograms can be used for the combined diagnosis or disease status prediction with multiple indicators. In this situation, accuracy was of the essence, so we were prompted to construct a novel nomogram combing the CSGscore model and multiple clinicopathologic traits to constitute a quantitative tool for predicting the clinical outcomes of HCC patients. While univariate analysis identified the CSGscore classifier, stage, CSG subtype, and tumor burden were prognostic variables, in multivariate analysis, only CSGscore and tumor stage were proved to be independent prognostic indicators for HCC (Figure 7B). Afterward, by integrating the prognostic factors including CSGscore, stage, tumor burden, and CSG subtypes, we established an integrative nomogram to predict the 1-, 3-, and 5-year OS of HCC patients (Figure 7C). Well-fitted calibration curves demonstrated that the nomogram had excellent consistency between the predicted and observed 3- and 5-year OS (Figure 7D). Moreover, The 1-, 3-, and 5-year DCA curves showed that the CSGscore-integrated nomogram could yield better net benefits compared with other indices for the vast majority of threshold probabilities (Figures 7E–G). Importantly, Kaplan-Meier analyses of OS, DFS, PFS, and DSS for low- and high-risk HCC patients classified by the median value of the nomogram score suggested that the CSGscore-integrated nomogram was a utilizable tool in evaluating the prognosis of HCC patients (Figure 7H).

Ongoing efforts are being made to decipher the interplay between cellular senescence and immune response in the tumor microenvironment to determine whether senescence changes the role of the immune system from anti-tumor to pro-tumor response. In this context, we eagerly mined the potential correlation between CSGscore and the tumor immune microenvironment of HCC. We found most of the immune checkpoint genes (ICGs) were highly expressed in the CSGscore-defined high-risk group compared to the low-risk group (Figure 8A). Moreover, Immunophenoscore (IPS), an index to measure the overall immunogenicity of tumors with machine learning, was designed to predict patients’ response to immunocheckpoint inhibitors (ICIs) therapy, and a higher IPS was representative of a more immunogenic tumor. As Figure 8B shows, in the whole TCGA-LIHC cohort, the low-risk patients had higher scores of IPS, IPS-CTLA4, and PD1/PD-L1/PD-L2 blocker than high-risk patients, which demonstrated that low-risk patients had higher sensitivity to ICIs therapy. Also, we performed TIDE scoring on patients in the TCGA-LIHC cohort and the ICGC-LIRI-JP dataset, both of which showed the low-risk patients had significantly lower TIDE scores and a higher proportion of predicted responders (Figures 8C,D). By the CIBERSORT algorithm, we calculated the correlation matrix of the expression levels of the five CSGscore genes and 22 immune cells, respectively. Interestingly, we found that the CSGscore genes were all significantly associated with the M0 macrophage infiltration, positively or negatively in HCC, corresponding to their hazard ratio of them inside the model (Figure 8E). The relative abundance of the 22 immune cell types for the whole TCGA cohort was also shown in Supplementary Figure S5.

In addition, we evaluated the potential application of CSGscore for predicting the HCC patients’ responses to chemo-/targeted therapy. The IC₅₀ values of 138 chemo-/targeted drugs were calculated by “pRRophetic” algorithm for HCC patients and 10 common chemo/targeted drugs were selected to show the different IC₅₀ values of them in low- and high-CSGscore groups (Supplementary Figure S6).

Significant correlations of cellular senescence levels with tumor stemness have been reported previously (Wang et al., 2022), thus, we suppose the statistical linkage between CSGscore and stemness indices. As shown in Figure 9A and Figure 9B, mRNAsi and mDNAsi both showed significantly different between the low- and high-risk groups, and we noticed that patients in the high-risk group had higher tumor stemness indices, suggesting a worse prognosis in patients with a higher stemness level. Specifically, four CSGscore genes (STC2, SOCS2, CCT5, and PSRC1) were found significantly correlated with mRNAsi while three CSGscore genes (SOCS2, STC2, and TMEM45A) had significant correlations with mDNAsi (Figures 9C,D). Further, we validated the differential expression of the five signature genes with paired boxplots using the TCGA-LIHC dataset (Figure 9E). To seek the clinical relevance of the CSGscore, we analyzed their expression levels according to tumor grade, pathologic stage, tumor status, and CSG subtype. Consequently, most of them showed statistical differences in different subgroups of tumor grade and pathologic stage (Figures 9F,G), implying that they may good indicators for tumor progression. While only TMEM45A was associated with tumor burden status, all of the five signature genes were strongly related to CSG subtypes (Supplementary Figure S7), suggesting their close relationships with the CSG patterns of HCC.