Research and experimental verification on the mechanisms of cellular senescence in triple-negative breast cancer

Data downloading and pre-processing

Single-cell sequencing dataset, GSE176078, comprising nine samples was extracted from the NCBI data Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). The statistical power (Therneau, Hart & Kocher, 2023) of this experimental design, calculated in RNASeqPower is 1.0. Single-cell data were subjected to filtering, quality control, clustering, and dimensionality reduction using the “Seurat” (Gribov et al., 2010) package. Microarray data of GSE58812 was extracted from GEO, and probes were transformed to “symbols” using the annotation file. In total, 107 qualified tumor samples and 16,416 genes were obtained after excluding tumor samples and normal tissues without information on clinical follow-up and overall survival (OS).

Clinical phenotypic data of TNBC were extracted from The Cancer Genome Atlas (TCGA) database. Samples without information on survival time and statuses were removed, and all those with patient survival time >0 days were retained. Subsequently, expression profile data were downloaded from TCGA. Finally, 113 para-cancerous and 113 tumor samples were obtained. The statistical power of this experimental design, calculated in RNASeqPower is 1.0. Furthermore, copy number variations (CNVs) on the Masked Copy Number Segment type of TNBC were obtained from TCGA database and integrated using “gistic2” (Mermel et al., 2011) software. Finally, the single nucleotide variants (SNVs) in TCGA-TNBC dataset processed by “mutect2” software, were obtained from TCGA database. Cellular senescence-related pathways were extracted from the Molecular Signature Database (MSigDB: https://www.gsea-msigdb.org/gsea/index.jsp). Sangerbox platform (http://vip.sangerbox.com/) was introduced to assist bioinformatics analysis in this study (Shen et al., 2022).

Single-cell clustering and dimensionality reduction

First, 38,582 cells were obtained after filtering the single-cell RNA-seq data matrix according to the cell criterion (<250 transcripts/cell) and gene criterion (<3 cells/gene). The “PercentageFeatureSet” function was then used to determine the percentages of mitochondria and rRNAs. Ultimately, 38,007 cells that met the following inclusion criteria were obtained: One cell could express between 100 and 6,000 genes, the mitochondrial content was less than 25%, and each cell had a unique molecular identifier (Yasuda et al., 2021) of at least 100.

Variance stabilization transformation (VST) was used to identify variable characteristics after the aforementioned data were log-normalized. High-variance genes were then discovered by employing the “FindVariableFeatures” tool. Further, the “FindIntegrationAnchors” function was utilized to remove batch effects using the canonical correspondence analysis (CCA) method for 10 samples. The data were integrated using the “IntegrateData” function. The anchor (dim = 30) was determined by principal component analysis (PCA) and dimensionality reduction after normalization using the “ScaleData” function for all genes. The cells were clustered using “FindNeighbors” and “FindClusters” functions with resolution = 0.1; this yielded 11 subpopulations. Finally, the TSNE dimensionality reduction analysis was performed on 38,582 cells using the “RunTSNE” function.

With threshold parameters of logfc = 0.5 (fold-change), Minpct = 0.35 (minimum percentage of differential gene expression), and adjusted P0.05, marker genes of subpopulations were screened using the “FindAllMarkers” tool. Subsequently, changes in cellular CNVs in single-cell data were predicted using the “copycat” (Gao et al., 2021) package to distinguish tumor cells. Although cancerous and normal tissues were differentiated during sampling, cancerous tissues might contain normal cells. Therefore, a distinction between the two was necessary.

Gene set enrichment analysis (GSEA) and annotation

Cellular senescence-related pathways were obtained from the GSEA website (https://www.gsea-msigdb.org/gsea/index.jsp). We separately calculated the corresponding scores of non-malignant and malignant cells in cellular senescence-related pathways by ssGSEA using the “GSVA” package and normalized the enrichment scores for each pathway by z-score. The scores of each sample of normal and tumor tissues in the bulk RNA-seq dataset were analyzed by ssGSEA for these senescence-related pathways (Subramanian et al., 2005). Finally, the significance of each senescence-related pathway in cancer and para-cancerous tissues was assessed using “wilcox.test”.

The genes related to the G1/S phase (p15-Cell cycle G1/S, MDM2-p21-Cell cycle G1/S, and p27-Cell cycle G1/S; 27 genes in total) were extracted from the Kyoto Encyclopedia of Genes and Genomes (KEGG) official website. Subsequently, the scores for G1/S genes were calculated for each sample in TCGA dataset by “ssGSEA”. Genes in the HALLMARK_G2M_CHECKPOINT pathway were downloaded from MSigDB via GSEA and the scores of G2 checkpoints were calculated for each sample in TCGA dataset by “ssGSEA”. We downloaded REACTOME_TELOMERE_EXTENSION_BY_TELOMERASE using GSEA and the primary role of this pathway was “Telomere Extension By Telomerase”. Finally, GSEA-derived enrichment scores in TCGA dataset for epithelial-mesenchymal transition (EMT) were computed by ssGSEA for 200 genes in the HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION pathway for each sample (Liberzon et al., 2015).

Univariate Cox and LASSO regression analyses

Univariate Cox analysis was conducted using the “coxph” function in the survival package of R to screen the significant (p < 0.05) prognosis-related genes. LASSO-Cox regression analysis was conducted for prognosis-related genes using the “glmnet” package in R (Friedman, Hastie & Tibshirani, 2010). The changing trajectory of each independent variable was analyzed. Increased lambda was positively related to increased number of independent variable coefficients tending to zero. Finally, by 10-fold cross-validation, a model was constructed; for each lambda, the confidence intervals were calculated.

The risk-related prognostic score (Riskscore) was calculated for each sample according to the equation for the sample risk score: Riskscore = Σβi×Expi, where Expi denotes the level of gene expression of the corresponding gene signature and β is the corresponding Cox regression coefficient. ROC analysis and z-score for RiskScore were computed using the “timeROC” package in R. Samples with RiskScore lower than zero comprised the low-risk group, while those with RiskScore greater than zero after z-score normalization were classified into the high-risk group. Finally, the Kaplan–Meier curves were plotted.

Consensus clustering

Consensus clustering was conducted using the “ConsensusClusterPlus” package (Wilkerson & Hayes, 2010). In total, 500 bootstraps were conducted using the “hc” algorithm with “canberra” as the metric distance; in the training set, 80% of the patients were in each bootstrap process. After the number of clusters was selected to between 2 and 10, “ConsensusClusterPlus” in the TCGA clustered the 113 TNBC samples dataset. The optimal classification was confirmed based on the consensus matrix and cumulative distribution function (CDF).

Mutation analysis

Using hg38 as the reference genome, the CNV-related TCGA-TNBC results were integrated by “gistic2″software at a confidence level of 0.9. The downloaded SNV data from TCGA were analyzed using the “maftools” (Mayakonda et al., 2018) package.

Immune cell scoring and immunotherapy

The scores of the relevant immune cells were calculated using the “MCPcounter.estimate” function of the “MCPcounter” package (Becht et al., 2016), Tumor Immune Estimation Resource (TIMER) (Li et al., 2020) and Estimating the Proportion of Immune and Cancer cells (EPIC) (Zhou et al., 2020) methods. Significant differences were determined using “kruskal.test”. Immune cell infiltration was assessed using the “Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression” (ESTIMATE) (Yoshihara et al., 2013) which provides information on tumor purity, scores of immune cell infiltration, and stromal cell levels in tumor tissues. Subsequently, the TIDE (Jiang et al., 2018) software was used to assess the putative immunotherapeutic efficacy in the defined molecular subtypes (http://tide.dfci.harvard.edu/). Higher TIDE score indicates less immunotherapy benefit as it suggests a greater possibility of immune escape. The correlation and significance of the risk score with the immune cell score were separately calculated using the “Hmisc” package’s “rcorr” function based on Pearson’s method.

Differential expression analysis

Differential expression analysis of cluster 1, 2 and 3 was performed in both GSE and Target datasets in “limma” package (Ritchie et al., 2015). Finally, differentially expressed genes (DEGs) were filtered under “—log2 (Fold Change)— > 1 and p < 0.05”.

Quantitative reverse transcription-polymerase chain reaction

Using TRIzol reagent (Thermo Fisher, Waltham, MA, USA), total RNA was collected from MDA-MB-468, MDA-MB-231 and MCF10A cell lines with 260/280 the value was between 1.8–2.0. High-Capacity cDNA Reverse Transcription Kit (4368814, ThermoFisher, Waltham, MA, USA) was conducted for reverse transcription. Using a LightCycler 480 PCR System (Roche, Indianapolis, IN, USA), the RNA from each sample (2 µg) was subjected to qRT-PCR with FastStart Universal SYBR^®Green Master (Roche, Indianapolis, IN, USA). The cDNA was a template with a reaction volume of 20 µl (0.5 µl of forward and reverse primers, 2 µl of cDNA template, 10 µl of PCR mixture, and appropriate volume of water). For the PCR reactions, cycling began with an initial DNA denaturation phase for 30 s (s) at 95 °C, then a total of 45 cycles were run for 15 s at 94 °C, for 30 s at 56 °C, and for 20 s at 72 °C. Each sample was performed for three separate analyses. Data from the threshold cycle (CT) were normalized to the level of GAPDH by 2^−ΔΔCT. The mRNA expression was compared in the normal and control tissues. Sequences of primer pairs targeting (Shanghai Gemma Gene Co., LTD, Shanghai, China) the genes are listed as follows:

Cell culture and transient transfection

Breast cancer cell lines MDA-MB-468, MDA-MB-231, and non-tumorigenic epithelial cell lines MCF10A were suspended in serum-free cell freeze and immediately frozen in liquid nitrogen tanks. MDA-MB-468, MDA-MB-231 and MCF10A cell lines were commercially obtained from Beijing Bena Biotechnology Co. (Beijing, China). Cells were cultured in DEME F-12 medium (Gibco, Waltham, MA, USA). Lipofectamine 3000 (Invitrogen, Waltham, MA, USA) was used for transfecting the ACP5 siRNA (Sagon, China) and negative control (NC) into the cells. The target sequences for ACP5 siRNA were TCCTAAATCAAGCATCTTTCTGT (si ACP5).

Transwell assay

Migration and invasion of MDA-MB-468 and MDA-MB-231cell lines were detected. Cell (5 × 10⁴) inoculation onto chambers coated (for invasion) or uncoated with Matrigel (BD Biosciences, USA) (for migration) was performed. The upper layer was added with serum-free medium, while the lower layer was added with complete DMEM medium. After incubation for 24 h, migrated or invaded cells were fixed using 4% paraformaldehyde and dyed by 0.1% crystalline violet. Cell number was counted using a light microscope. The experiment was conducted in our own lab.

Statistical analysis

All statistical analysis was done with R program (v4.2.0). The Wilcoxon test evaluated two-group difference. Differences among three groups were analyzed by the Kruskal–Wallis test. Cox analysis and survival analysis both used the log-rank test. For Wet experiments, the number of each group was 3 independent experiments. Assay carried out by investigator’s lab. Consumables in experiments we use are RNAase and DNAase free.

Nine annotated cell types displayed different senescence characteristics

A cellular senescence-related classifier was constructed for TNBC. The prognoses of patients with TNBC were verified using bulk RNA-seq and scRNA-seq data. The flow chart of the study design is presented in Fig. S1.

First, 38,582 cells were obtained by screening the data from nine samples (Table S1, Fig. 1A); these were clustered and 11 subpopulations were obtained (Fig. 1B). Subsequently, these cells were subjected to TSNE dimensionality reduction analysis, and the 11 subpopulations were annotated based on some classical markers of immune cells (Fig. 1C and Table S2). Among them, subpopulation C1 comprised CD8 T cells expressing CD8A, CD3D, GZMA, and CD8B; C6 comprised CD4 T cells expressing CD4 and CD3D; C2 and C9 comprised macrophages expressing CD163 and CD68; C3 and C10 comprised monocytes expressing S100A8; C5 comprised B cells expressing CD19, CD79A, and MS4A1; C4 comprised plasma cells expressing CD79A and JSRP1; C0 comprised epithelial cells expressing EPCAM; C3 and C8 contained fibroblasts expressing ACTA2, FAP, PDGFRB, and NOTCH3, and C7 contained endothelial cells expressing PECAM1. We also screened marker genes among these subpopulations and analyzed their corresponding expression (Fig. 1D). For instance, CD24 and KRT19 were highly expressed in epithelial cells in C0 and were the marker genes of C0.

The proportion of these nine subpopulations in each sample was further analyzed (Fig. 1E). Subsequently, changes in CNVs in single cells were predicted to distinguish normal from tumor cells in each sample (Fig. 1F). In total, 7,709 cancer cells and 30,298 normal cells were identified. Subsequently, the proportion of malignant and non-malignant cells in each sample was calculated (Fig. 1G). The cells in samples CID3946, CID44041, and CID4465 were all non-malignant, while the other samples comprised both malignant and non-malignant cells.

The senescence characteristics of single cells were further assessed. Cellular senescence-related pathway scores were higher in malignant relative to the non-malignant cells (Fig. 2).

Validation of abnormal cellular senescence based on bulk RNA-seq data

Cellular senescence-related pathways were highly expressed in malignant cells at the single cell level relative to the non-malignant cells, with aberrant expression. Therefore, their expression in tumor and normal tissue samples was further analyzed using a bulk RNA-seq dataset. GSEA results suggested that REACTOME_CELLULAR_SENESCENCE, and P53_SIGNALING_PATHWAY were significantly enriched in tumor tissues in TCGA dataset (Fig. 3A). ssGSEA of the bulk RNA-seq dataset (Fig. 3B) revealed that the enrichment scores of REACTOME_DNA_DAMAGE_TELOMERE_STRESS_INDUCED_SENESCENCE, P53_SIGNALING_PATHWAY, REACTOME_CELLULAR_SENESCENCE, REACTOME_DNA_DAMAGE_TELOMERE_STRESS_INDUCED_SENESCENCE were higher in some cancer tissues relative to the para-cancerous tissues. Therefore, these three cellular senescence-related pathways may exhibit regulatory complexity in cancerous tissues.

Construction of cellular senescence-related subtypes

Three cellular senescence-related pathways, including REACTOME_CELLULAR_ SENESCENCE, and P53_SIGNALING_PATHWAY, REACTOME_DNA_DAMAGE_ TELOMERE_STRESS_INDUCED_SENESCENCE, were significantly enriched in TNBC tissues. Therefore, we further analyzed the genes enriched in these three pathways.

The above pathways comprised 253 genes, of which 186 were present in TCGA dataset. Subsequently, a univariate analysis of these 186 genes yielded seven prognosis-related genes, including ETS2, SERPINE1, FOS, SHISA5, IL1A, TP53AIP1, and IGFBP7 (Fig. 4A). Mutations in these seven genes in TNBC were further examined. Among SNVs, ETS2 showed the highest mutation frequency (primarily mutation type: missense) (Fig. 4B). Among CNVs, IGFBP7 had the highest frequency of acquired CNVs (CNV_gain), whereas SHISA5 had the highest frequency of deletion mutations (CNV_loss) (Fig. 4C). The levels of these seven genes in normal versus tumor tissues were analyzed by the “wilcox.test” (Fig. 4D). The expressions of SHISA5, SERPINE1, and IL1A were markedly high in tumor tissues, whereas those of TP53AIP1, ETS2, and FOS were high in the normal tissues.

Consensus clustering of 113 TNBC samples in TCGA dataset was performed based on these seven key prognosis-related genes. Stable clustering results were achieved at k = 3 (Figs. 4E–4F). Therefore, three subtypes (clusters) were obtained (Fig. 4G). Significant differences in the prognostic characteristics among the subtypes were found (Fig. 4H). Overall, cluster 1 showed the best prognosis, followed by clusters 2 and 3. These seven genes showed the highest expression in cluster 3 (Fig. 4I).

Differential analysis of cellular senescence-related subtypes

The expression of genes in these three senescence-related pathways among the three clusters was compared by “kruskal.test” (Fig. S2). Among the three pathways, genes in cluster 2 showed the lowest expression, while those in cluster 3 showed the highest expression.

The distribution of clinical characteristics across molecular subtypes was compared by chi-square test to assess the differences among the three subtypes in TCGA. The survival statuses of these patients differed significantly among the three subtypes. In particular, the proportion of dead patients was high in cluster 3.

Differences in mutations among cellular senescence-related subtypes were analyzed. TCGA data suggested some differences in CNVs among the three subtypes (Figs. 5A and 5C). For instance, cluster 1 had the lowest degree of gene gain and lose. While cluster 3 exhibited the highest degree of CNV gain. Additionally, clusters 2 and 3 had similar CNV loss on chromosomes 2, 3, 4 and 5. Furthermore, TP53, TTN, MUC16, SYNE1, and FAT3 had the highest SNV frequency (Fig. 5B).

Biological characteristics of cellular senescence-related subtypes

Cancer cells can induce cellular senescence through the inhibition of the cell cycle, and cell cycle protein-dependent kinases (CDK inhibitory proteins), such as InK4a and p21, are upregulated in senescent cells, which causes cell cycle arrest (Wang, Lankhorst & Bernards, 2022). A previous study Cuzick et al. (2011) identified 31 genes associated with cell cycle progression (CCP). We computed the CCP scores of each sample in TCGA dataset by ssGSEA. clusters 1 and 2 had slightly higher CCP scores relative to cluster 3 (wilcox.test) (Fig. 6A). No significant differences were observed among the three subtypes in the G1/S phase (wilcox.test) (Fig. 6B), whereas the scores for G2 checkpoints were slightly lower in cluster 3 relative to clusters 1 and 2 (wilcox.test) (Fig. 6C). Prognostic analysis suggested that patients in cluster 3 showed a worse prognosis. Therefore, the cell cycle may be one of the factors influencing cellular senescence. Other mechanisms may also regulate cellular senescence along with the cell cycle.

Similarly, inhibition of telomerase induces cellular senescence (Wang, Lankhorst & Bernards, 2022). Cancer cells usually circumvent telomere attrition by activating telomerase activity. Telomere elongation scores for telomerase in cluster 3 were lower than those in clusters 1 and 2, with the latter showing the highest score. Therefore, other mechanisms may co-regulate cellular senescence with telomerase activity (wilcox.test) (Fig. 6D).

However, in addition to conveying the message “please kill me”, the factors secreted by senescent cells also affect surrounding cells (Yu et al., 2021). This effect can promote tumor migration and metastasis by promoting EMT. Senescent tumor cells can promote the production of blood and lymphatic vessels by recruiting specialized macrophages that provide oxygen and nutrients necessary for growth to other tumor cells, thus promoting tumor growth and metastasis. We calculated the EMT scores. Cluster 2 had a lower EMT score than clusters 1 and 3, indicating that patients in cluster 3 were more likely to experience metastasis (wilcox.test) (Fig. 6E).

Based on genes in the HALLMARK_HYPOXIA pathway, the hypoxia scores of the samples were computed by ssGSEA. Angiogenesis scores of samples were also analyzed based on 24 genes from the literature (Masiero et al., 2013) (wilcox.test) (Figs. 6F–6G). Enrichment scores for 10 pathways associated with tumors in each sample in TCGA dataset were computed by the ssGSEA method described previously (Sanchez-Vega et al., 2018). Only two of the ten pathways (NOTCH and TGF-beta signaling) exhibited significant differences (wilcox.test).

Finally, genes in the HALLMARK_INFLAMMATORY_RESPONSE pathway from GSEA-based MSigDB were analyzed for inflammation-related scores by ssGSEA. Cluster 2 had a significantly lower inflammation-related score relative to clusters 1 and 3 (wilcox.test) (Fig. 6H).

Immune characteristics of cellular senescence-related subtypes

The relationship between cellular senescence-related subtypes and immunity was analyzed to investigate the immune characteristics of these subtypes. The immune scores of patients in TCGA dataset were calculated by “ESTIMATE”. Significant differences were observed among the three subtypes, with a high degree of immune infiltration in clusters 1 and 3 (kruskal.test) (Fig. 7A).Subsequently, the scores for the relevant cells were calculated. Significant differences were noted for all cells except CD8 T cells and cytotoxic lymphocytes (kruskal.test) (Fig. 7B). Similarly, T cells CD4, also displayed distinct difference among three groups in TIMER and EPIC analysis (Figs. 7C–7D). Although the results were obtained from different algorithms, these results all together suggested higher immune infiltration status in clusters 3. The expression of immune checkpoint genes, including NRP1, CD200, BTLA, LAIR1, and TNFRSF14 differed significantly among the three subtypes (kruskal.test) (Fig. 7E).

Finally, the putative clinical efficacy of immunotherapy was assessed among these molecular subtypes. High TIDE scores in cluster 3 in TCGA cohort suggested that patients of this subtype were less likely to benefit from immunotherapy owing to a higher likelihood of immune escape (wilcox.test) (Fig. 7F).

Identification of prognosis-related genes and construction of the risk model

Differential analysis was separately performed between cluster 1 and non-cluster 1; cluster 2 and non-cluster 2, and cluster 3 and non-cluster 3. Finally, 2 down-regulated and 34 up-regulated DEGs were screened between cluster 1 and non-cluster 1; 293 down-regulated and 9 up-regulated DEGs between cluster 2 and non-cluster 2, and 26 down-regulated and 218 up-regulated DEGs between cluster 3 and non-cluster 3. Thus, a total of 391 DEGs were screened.

Univariate Cox analysis of these 391 DEGs yielded 69 genes with high prognostic significance (p < 0.05), including 63 “Risk” and six “Protective” genes (Fig. S3A). These 69 hub genes were subsequently screened by LASSO regression (Fig. S3B). Further, using 10-fold cross-validation, a model was constructed and confidence intervals for each lambda were calculated (Fig. S3C). The model was optimum for lambda = 0.0782, whereby four genes were obtained, including MMP28, CT83, ACP5, and KRT6A; these were selected as the target genes for further analyses.

The four genes’ expression levels were utilized to determine the risk score for each sample using the TCGA dataset as the training set. Then, using ROC analysis, the RiskScore prognostic classification was carried out (Fig. 8A). The prognostic prediction classification efficiency was analyzed separately for 1–5 years with obtained AUC values 0.91, 0.93, 0.74, 0.83 and 0.84, respectively , and the AUC values reached 0.7, indicating good predictive performance. Samples with RiskScore greater than zero after z-score normalization were classified into the high-risk group, whereas those with Riskscore lesser than zero comprised the low-risk group, with a highly significant (p < 0.05) difference between them (Fig. 8B).

To validate the model’s robustness, the GSE58812 dataset was used for verification. Specifically, the risk model was constructed using these four genes and ROC analysis for prognostic classification was performed based on the RiskScore (Fig. 8C). The prognostic predictive classification efficiency was analyzed for 1–5 years with obtained AUC values 0.98, 0.71, 0.62, 0.66 and 0.65, respectively. All AUC values were over 0.6, displaying good predictive ability. Furthermore, a significant difference in prognosis was found between the risk groups (p < 0.05) (Fig. 8D).

Validation of association between prognostic model signature and TNBC malignant phenotype

In this study, molecular assays revealed that the expression of MMP28, ACP5 and KRT6A was increased in MDA-MB-468 and MDA-MB-231 cell lines. In contrast, the expression of CT83 was decreased in MDA-MB-468 and MDA-MB-231 cell lines (Figs. 9A–9D). This result is consistent with previous analyses that high expression of prognostic risk factors and low expression of protective factors associated with higher RiskScore. Given that ACP5 is a key gene influencing high RiskScore, we inhibited the expression of ACP5 in MDA-MB-468 and MDA-MB-231 cell lines, and by transwell assay, we could observe the ability of ACP5 to promote the proliferation of triple-negative breast cancer cell lines. After we inhibited the expression of ACP5 in MDA-MB-468 and MDA-MB-231 cell lines, the invasion and migration ability of the cell lines were decreased (Figs. 9E–9G). This study not only confirms the existence of ACP5 as a prognostic risk factor, but also demonstrates that this gene regulates the TNBC malignant phenotype at the cellular level, thereby influencing TNBC malignant progression.

Immunological characteristics of RiskScore subgroups and sensitivity analysis based on drug IC50

To elucidate differences in the immune microenvironment of patients between the risk groups, differences in the relative cell abundances predicted by “MCPcounter” were compared (wilcox.test) (Fig. 10A). T cells, cytotoxic lymphocytes, and monocytic lineage differed significantly between the risk groups. Immune cell infiltration was also assessed (wilcox.test) (Fig. 10B). “ImmuneScore” was higher in the high-risk group relative to the low-risk group, suggesting that patients in the former had high immune cell infiltration. Additionally, TIMER and EPIC analysis also revealed higher immune infiltration status in high risk group (Figs. 10C–10D).

Subsequently, differences in response to immunotherapy between the risk groups were analyzed in TCGA cohort. First, differences in immune checkpoint expression were compared. Most immune checkpoint genes, including TNFRSF14, LAIR1, and CD244 showed differential expression between the risk groups, with high expression in the high-risk group (wilcox.test) (Fig. 10E).

A significant positive correlation between the Riskscore and T cells, cytotoxic lymphocytes, and monocytic lineage (Fig. 10F) was observed. Risk scores also correlated positively with “ImmuneScore”, “StromalScore”, and “ESTIMATEScore” (Fig. 10G). Drug sensitivity prediction revealed that low-risk group showed higher sensitivity to six traditional drugs (MG-132, Sorafenib, A77004, Bortezomib, Shikonin and AZ628) by wilcox.test (p < 0.05, Fig. 10H).

Combination of RiskScore and clinicopathological characteristics improves the performance of the prognostic model and its survival predictive accuracy

Univariate regression analysis based on RiskScore and other clinicopathological characteristics suggested that the latter was the most significant prognostic factor (p-value = 0.007, Fig. 11A) and remained so (p-value = 0.003, Fig. 11B). A nomogram was constructed to quantify the survival probability of patients and risk assessment by combining RiskScore and other clinicopathological characteristics (Fig. 11C). The results confirmed the greatest impact of RiskScore on survival prediction. Subsequently, the model’s prediction accuracy was assessed using a calibration curve (Fig. 11D). The findings indicated that the calibration curves for 1-, 3-, and 5-year predictions nearly overlapped with the standard curves, pointing to the nomogram’s strong predictive ability. Furthermore, decision curve analysis (DCA) was used to evaluate the model’s dependability. DCA results highlighted that both RiskScore and nomogram had significantly more benefits than the extreme curves with strong survival predictive power relative to other clinicopathological features, especially for the 5-year survival (Figs. 11E–11F).