Precision prognostication in breast cancer: unveiling a long non-coding RNA-based model linked to disulfidptosis for tailored immunotherapeutic strategies

The procurement and compilation of sample data

A comprehensive dataset, encompassing 1097 instances of breast cancer (BRCA), was acquired from the TCGA database (https://portal.gdc.cancer.gov/). This dataset amalgamated gene expression profiles with relevant clinical particulars for each specimen, encompassing survival status, age, gender, survival duration, and TNM stage. In the initial preprocessing of transcriptomic data, a conversion from FPKM (Fragments Per Kilobase Million) to TPM (Transcripts per million) was executed [23]. Following this, the processed expression values underwent a logarithmic transformation with a base of 2. This transformation facilitated the amalgamation of genes and their corresponding expression values, leading to the creation of comprehensive gene expression profiles for each sample.

Concerning clinical information, a streamlined process of downloading and merging relevant datasets sufficed for the subsequent analysis. This meticulous approach ensures the robust integration of gene expression data and clinical parameters, laying a solid foundation for the in-depth exploration of breast cancer within this cohort.

Development and validation of the model

The dataset sourced from the TCGA database underwent a meticulous division into two sets, maintaining a balanced 1:1 ratio, to fulfill distinct roles in the analytical process. The training set assumed the pivotal role of model construction, while the validation set functioned as a critical evaluator of the model’s performance. To pinpoint noteworthy lncRNAs, a univariate Cox regression analysis was executed using the “glmnet” package in R. Subsequently, a LASSO regression approach, incorporating a tenfold cross-validation strategy, was employed to pinpoint the optimal point of minimal error [24]. This stringent methodology facilitated the discernment of pivotal genes and their corresponding parameters via a multivariate Cox regression analysis, leading to the revelation of eight differentially expressed long non-coding RNAs (DALs). These DALs formed the basis for constructing essential models, wherein each patient’s risk score was meticulously calculated. This score was computed by summing the product of the expression level of each long non-coding RNA (lncRNA) and its respective coefficient, offering a quantitative assessment of risk for each participant in the cohort.

Correlation analysis of DAL with different clinical features

To scrutinize the correlation between DALs and diverse clinical characteristics, a comprehensive exploration was undertaken. Utilizing the “Limma” package within the R software, a rigorous examination of differential gene expression was executed on the training dataset. The primary objective of this analysis was to discern variations in the expression levels of DALs associated with clinical indicators and risks. For a nuanced portrayal of the diversity in DAL expression across distinct clusters, we harnessed the capabilities of the “pheatmap” package in R, thereby generating heatmap representations. These depictions offer a visual elucidation of the disparate expression patterns of DALs concerning various clinicopathological features. Such visual aids serve to illuminate the intricate relationship between DAL expression and diverse clinical characteristics.

Model equations

To scrutinize the correlation between differentially expressed long non-coding RNAs (DALs) and various clinical characteristics, we conducted a thorough analysis of gene expression differences within the training dataset, utilizing the “Limma” package in R. The overarching goal of this analysis was to discern discrepancies in the expression levels of DALs associated with clinical indicators and risks. For a graphical representation of the nuanced variations in DAL expression across distinct clusters, we leveraged the capabilities of the “pheatmap” package in R, facilitating the generation of heatmap depictions. These visual representations vividly illustrate the distinct patterns of DAL expression in correlation with clinicopathological features. These visual aids constitute a valuable method for extracting insights into the intricate connection between the expression of differentially expressed DALs and diverse clinical characteristics.

Development of an independent prognostic model and construction of a nomogram

To establish the independent prognostic relevance of the risk score, we performed both univariate and multivariate Cox regression analyses. These assessments aimed to assess the association between the risk score and patient survival outcomes, concurrently accounting for other clinicopathological variables as covariates. Furthermore, employing the ‘rms’ package in R, we crafted histograms to visually illustrate the distribution of risk scores and clinicopathological variables. These graphical representations provided valuable insights into the interrelation between risk scores and various clinical characteristics. Subsequently, based on the results of the prognostic model, we devised a nomogram. Functioning as a graphical predictive instrument, this nomogram amalgamates the risk score with other pertinent clinical variables to estimate personalized survival probabilities for individuals within the TCGA-BRCA cohort.

Investigation of enriched biological functions

To glean insights into the functional annotation and enrichment pathways linked to the eight differentially expressed genes associated with DALs in BRCA, a comprehensive functional enrichment analysis was undertaken. For Gene Ontology (GO) analysis, the ClusterProfiler package in R was employed. Furthermore, GSVA analysis, utilizing the “c2.cp.kegg.v7.4.symbols.gmt” package from MSigDB and the “GSVA” package in R, was conducted to discern pathway disparities among distinct clusters. The outcomes of these analyses were visually presented using the ‘heatmap’ package in R. In addition, Gene Set Enrichment Analysis (GSEA) was carried out, employing the “c2.cp.kegg.Hs.symbols.gmt” package, to delve deeper into enriched pathways associated with the differentially expressed genes. These analytical endeavors were geared towards unraveling the biological functions and pathways potentially influenced by the identified DALs in the context of BRCA.

Immunological scrutiny of risk attributes

To examine the immunological attributes associated with risk assessment, a suite of robust algorithms, including EPIC, MCPCOUNTER, CIBERSORT, CIBERSORT-ABS, TIMER, XCELL, and QUANTISEQ, was deployed for the analysis of immune infiltration scores. Subsequently, using the “limma” package, we performed Spearman’s correlation analyses aimed at elucidating the relationship between the proportion of immune cells and risk assessment. Additionally, the CIBERSORT algorithm was utilized to categorize BRCA patients based on their immune cell characteristics.

In order to assess the enrichment of genetic features associated with the cancer immune cycle and treatment, we performed GSVA analyses on high-risk cohorts using the “clusterProfiler” package [25]. Finally, the “ggcor” R package facilitated the establishment of connections between risk scores and genetic traits. These comprehensive immunological analyses were performed with the overarching objective of offering profound insights into the immune characteristics intricately linked with the risk levels identified in our study.

Drug sensitivity

The Genomics of Drug Sensitivity in Cancer (GDSC) repository, available at https://www.cancerrxgene.org/ (accessed on 22 March 2023) [26], serves as a valuable database enabling the categorization of BRCA patients into low- and high-risk groups based on their half-maximal inhibitory concentration (IC50). This categorization offers crucial insights into the varying sensitivities of distinct patient populations to diverse treatments. To further scrutinize the response of these patient groups to specific therapeutic interventions, the R package “pRRophetic” proves instrumental. This package empowers the prediction of drug sensitivity grounded in genomic profiles, thereby enabling the personalized selection of treatments and the assessment of treatment responses in individuals affected by BRCA.

Somatic mutations

For the analysis of mutation information within BRCA samples, the TCGA-BRCA mutation database (https://portal.gdc.cancer.gov/, accessed on 22 March 2023) served as the primary data source. Subsequently, the Maftools program [27] was employed for the comprehensive analysis and visualization of mutation profiles. In our investigation, particular emphasis was placed on evaluating the risk of breast cancer through a comparative examination of mutational load, commonly referred to as the tumor mutation burden score. The TMB score was calculated by multiplying the mutation-to-covered bases ratio by 10^6, following the methodology outlined in the investigation by Chalmers and collaborators [28]. This approach enabled the quantification of the overall mutational burden, facilitating an assessment of its potential association with the risk of breast cancer within the analyzed samples. We also used the “limma” package to analyze the relationship between mapping tumor stem cells and risk scores.

Statistical analysis

All bioinformatics analyses were conducted leveraging the capabilities of R software version 4.2.2 and Strawberry Perl version 5.30.0. Spearman‘s correlation analysis method was employed to scrutinize the relationship between immune cell infiltration and risk scores. To compare overall survival (OS) outcomes between the high-risk and low-risk groups, Kaplan-Meier (KM) survival curves, complemented by log-rank tests, were utilized. Statistical significance was established at p-values < 0.05, with the control of false discovery rates (FDR) set at < 0.05.

Data availability statement

The datasets utilized in our investigation are accessible via the TCGA repository (http://cancergenome.nih.gov/). The unprocessed data can be obtained from jianguoyun at the following link: https://www.jianguoyun.com/p/DXU0NlwQo47JCxiEw4MFIAA.

Identification of lncRNAs correlated with disulfidptosis-induced cellular demise

Figure 1 presents a synopsis of the pivotal steps undertaken in this investigation. In the initial phase, an exhaustive analysis was performed on the TCGA-BRCA dataset to sift through protein-coding genes, leading to the discovery of 13,162 lncRNAs. Subsequent to this, co-expression analysis was executed between ten established genes linked to disulfidptosis and the lncRNAs associated with disulfidptosis. The aim was to discern those manifesting noteworthy co-expression relationships with disulfidptosis genes (Figure 2A). Moreover, the expression levels of these lncRNAs were extracted, and their affirmative and adverse regulatory associations were ascertained. Univariate Cox analysis was utilized to identify lncRNAs linked to elevated and diminished risk (Figure 2B). This analysis facilitated the identification of lncRNAs significantly linked to the risk levels. Following this, the LASSO algorithm was employed to choose distinctive lncRNAs, utilizing LASSO regression and observing their regression validation trends to pinpoint the juncture with the lowest cross-validation error (Figure 2C, 2D). Using a multivariate Cox proportional hazard regression model, eight disulfidptosis-related lncRNAs (MIR4435-2HG, AL139035.1, AL451085.2, AP001160.1, NRAV, AC015922.2, YTHDF3-AS1, TP53TG1) were chosen along with their corresponding regression coefficients for model construction. Furthermore, a comprehensive investigation into the correlation between the identified lncRNAs and the genes associated with disulfidptosis was undertaken (Figure 2E). The results unveiled a strong correlation between these long non-coding RNAs associated with disulfidptosis (DALs) and the genes implicated in disulfidptosis (Figure 2F).

Figure 1. The accompanying schematic diagram presents a thorough and all-encompassing summary of the primary cognitive processes conducted throughout the course of this study.

Figure 2. Construction and evaluation of a predictive model for disease onset. (A) A co-expression analysis was executed to explore the intricate interplay between genes associated with disulfidptosis and pivotal lncRNAs, forming the cornerstone of our model. (B) Univariate Cox analysis was employed to pinpoint differentially expressed lncRNAs and assess their correlation with high- and low-risk cohorts. (C) Utilizing the LASSO algorithm, integrated with 10-fold cross-validation, we identified the most significant lncRNAs linked to disulfidptosis. (D) The coefficients obtained from the LASSO algorithm were scrutinized to establish the foundation for our predictive disease model. (E) A correlation heatmap was generated to delve deeper into the intricate relationships between the selected lncRNAs and disulfidptosis-related genes. (F) Subsequent analysis of the interconnections among the chosen lncRNAs provided valuable insights into the underlying mechanisms of disulfidptosis.

Development and validation of the prognostic model

To formulate a resilient prognostic model, the cohort was stratified into two cohorts: the Training group and the Test group, ensuring an equitable distribution of samples in each. The Test group functioned as a validation set to appraise the model’s precision. The expression levels of the eight lncRNAs associated with disulfidptosis were multiplied by their respective regression coefficients for each sample. The resultant values were then aggregated to compute a risk score predicated on the lncRNA expression levels. Following this, samples were arranged in ascending order based on their risk scores, and the median score was employed to categorize them into high-risk and low-risk groups. This methodology was independently applied to all samples, encompassing both the Training and Test groups (Figure 3A–3C).

Figure 3. Model construction and evaluation of disease predictive value. (A–C) The dataset, comprising the overall, training, and testing sets, underwent stratification based on the risk score. Following this, samples were delineated into high-risk (depicted in red) and low-risk (depicted in blue) groups, utilizing the median risk score as the threshold. (D–F) The relationship between the risk score and both survival time and patient status was examined across the entire dataset, training set, and testing set. (G–I) Examination of the expression profiles of each Disulfidptosis-Associated LncRNA (DAL) was conducted, comparing high-risk and low-risk groups within the entire dataset, training set, and testing set. (J–L) Survival curves were juxtaposed between the high-risk and low-risk groups in the entire dataset, training set, and testing set to elucidate differences in survival outcomes. (M–O) Time-related Receiver Operating Characteristic (ROC) curve analysis was executed to appraise the predictive performance of the model across the entire dataset, training set, and testing set.

The correlation analysis between the risk score and survival status revealed a positive association, indicating that with an elevation in the risk score, there was a corresponding increase in the number of fatalities (Figure 3D–3F). Further examination of the expression patterns of the lncRNAs in the high-risk and low-risk groups, encompassing all samples, the Training group, and the Test group, disclosed that MIR4435-2HG, AL139035.1, and YTHDF3-AS1 lncRNAs exhibited an upregulated expression with a rising risk score in the high-risk group. In contrast, AL451085.2, AP001160.1, NRAV, AC015922.2, and TP53TG1 lncRNAs displayed a downregulated expression (Figure 3G–3I).

To evaluate the accuracy of the prognostic model, survival curves were plotted, and receiver operating characteristic (ROC) analysis was conducted at different time points. Consistently, the results demonstrated significantly lower survival rates in the high-risk group compared to the low-risk group at various time points, with statistical significance at P < 0.05. The constructed model effectively differentiated between the high-risk and low-risk groups (Figure 3J–3L). Furthermore, the ROC curves at different time points indicated that a larger area under the curve corresponded to higher predictive accuracy of the constructed model (Figure 3M–3O).

Principal component analysis (PCA) for gene expression

To delve further into gene expression patterns, we performed principal component analysis (PCA) on four distinct categories: all genes, genes associated with disulfidptosis, long non-coding RNAs (lncRNAs) linked to disulfidptosis, and the eight lncRNAs employed in formulating the prognostic model. The PCA analysis unveiled discernible separation among these groups based on their expression profiles. Notably, upon comparing the PCA plots of disulfidptosis-associated genes, all genes, the model lncRNAs, and disulfidptosis-associated lncRNAs, it became apparent that the model lncRNAs exhibited a more pronounced differentiation between the low and high-risk groups. This observation provides compelling evidence that the lncRNAs incorporated into the prognostic model possess the capability to effectively discriminate between low and high-risk sample groups (Figure 4A–4D).

Figure 4. Principal component analysis (PCA) was employed to scrutinize the differentiation between the high-risk and low-risk groups employing diverse gene sets. The four PCA plots depict the sample distribution within the corresponding risk groups based on distinct gene sets: (A) incorporating all genes, (B) comprising disulfidptosis-associated genes, (C) encompassing disulfidptosis-associated lncRNAs, and (D) involving model lncRNAs. The PCA analysis sought to appraise the unique clustering patterns and the potential discriminatory efficacy of these gene sets in distinguishing high-risk from low-risk groups.

Examination of the correlation between DAL and clinical characteristics

To investigate the association between the low and high-risk groups and diverse clinical features, we constructed heatmaps for a visual representation of the relationships. These heatmaps offered a comprehensive overview of the correlation between risk scores and clinical characteristics, encompassing tumor grade, gender, age, T, M, and N (Figure 5A). Moreover, we performed analyses to distinguish the distribution disparities of clinical features, including tumor grade, gender, age, T, M, and N, between the low and high-risk groups. (Figure 5B–5H). These analytical methodologies yielded valuable insights into potential connections between Disulfidptosis-Associated LncRNAs (DALs) and various clinical characteristics, thereby illuminating the potential clinical implications of our findings.

Figure 5. Correlation analysis of disulfidptosis-associated lncRNAs (DALs) with diverse clinical characteristics. (A) Heatmap illustrating the correlation between DALs and risk levels, along with various clinical attributes. (B) Age, (C) gender, (D) survival status, (E) staging, (F) T, (G) N, and (H) M distinctions between the high and low-risk cohorts.

Subgroup analysis of a disulfidptosis-associated lncRNA model in predicting clinical outcomes

To assess the prognostic implications of the DAL model within distinct clinical subgroups, a rigorous survival analysis was undertaken. The low and high-risk groups were stratified based on gender (male and female), age (≤50 and >50 years), clinical stage (I-II and III-IV), and TNM (T1-2 and T3-4, N0-1 and N2-3, M0 and M1). The findings consistently underscored that individuals categorized in the high-risk group experienced markedly diminished survival durations in comparison to their counterparts in the low-risk category across all clinical subgroups (p < 0.05) (Figure 6A–6L). These results affirm the resilience and widespread applicability of the established long non-coding RNA (lncRNA) model associated with disulfidptosis in forecasting clinical outcomes among varied patient cohorts.

Figure 6. The capacity of the DAL model to forecast survival variations in high and low-risk groups among distinct clinical subgroups. (A, B) age (>50 and ≤50 years), (C, D) gender, (E, F) M stage, (G, H) N stage, (I, J) clinical stage, and (K, L) T stage.

Independent prognostic analysis of clinical features and the creation of nomograms

To ascertain the autonomous prognostic efficacy of the formulated model, we executed independent prognostic analyses on clinical attributes. Univariate independent prognostic analysis was employed to scrutinize the association between singular clinical characteristics and the temporal and categorical aspects of survival. Additionally, a multivariate independent prognostic analysis was undertaken to appraise the collective impact of various clinical traits on survival metrics, taking into account plausible interactions between these factors. Our observations unveiled age, staging, and risk score as autonomous prognostic determinants, underscoring their potential to exert influence on patient outcomes irrespective of other variables (Figure 7A, 7B).

Figure 7. Independent prognostic analysis of clinical characteristics and nomogram creation. (A) Assessment of the independent prognostic impact of each clinical factor individually. (B) Evaluation of the independent prognostic influence of multiple factors for each clinical variable. (C) Nomograms were formulated incorporating each clinical indicator and risk scores to present a visual aid for predicting personalized patient outcomes. (D) Calibration curves were generated to appraise the nomograms’ accuracy, comparing predicted survival probabilities with actual survival rates at 1-year, 3-year, and 5-year intervals. (E) Calculation of the C-index was performed to gauge the discriminatory efficacy of gender, age, stage, risk scores, T category, N category, and M category in forecasting patient outcomes. (F) Receiver operating characteristic (ROC) analysis was undertaken for each clinical indicator to ascertain its predictive performance and capacity to differentiate between distinct outcome groups.

Moreover, a nomogram was formulated, amalgamating diverse clinical attributes to forecast patient survival. The calibration curves for 1, 3, and 5 years demonstrated a close alignment with the gray line segment, denoting a commendable concordance between projected and actualized survival probabilities (Figure 7C, 7D). Remarkably, the C-index value and the area under the ROC curve for the risk score surpassed those of alternative clinical traits, indicating the superior predictive accuracy of the devised model in prognosticating patient survival (Figure 7E, 7F). These results highlight the robustness and superior prognostic ability of the developed model, underscoring its potential clinical utility in predicting patient outcomes.

Enrichment analysis of prognostic markers in breast cancer patients

To unravel the functional implications of the identified prognostic markers, we conducted Gene Ontology (GO) enrichment analysis. The differential genes between the two groups were found to participate in diverse biological processes, encompassing ovulatory cycles and endocrine processes. At the cellular level, these genes were associated with blood particle composition and intermediate fibers, while at the molecular level, they exhibited activities such as peptide chain endonuclease activity and hormone functions (Figure 8A).

Figure 8. Enrichment analysis of prognostic markers in breast cancer patients. (A) Illustrates the pathway enrichment of differentially expressed genes concerning cellular components, molecular functions, and biological processes. (B) GSEA enrichment analysis reveals five upregulated functional pathways. (C) GSVA is employed to analyze the differential expression of various KEGG pathways between the two risk groups.

Moreover, GSEA highlighted a notable upregulation of pathways associated with immune rejection and antigen processing and presentation within the high-risk group (Figure 8B). For a holistic comprehension of differentially enriched pathways, GSVA enrichment analysis was executed. The outcomes unveiled heightened activity in pathways such as circadian rhythm in mammals, linoleic acid metabolism, and taurine and hypotaurine metabolism within the low-risk group. Conversely, the remaining pathways exhibited augmented activity in the high-risk group (Figure 8C).

These discoveries offer insightful revelations into the potential biological processes and pathways linked to the identified prognostic markers in individuals with breast cancer. The heightened activation of immune-related pathways in the high-risk cohort implies the engagement of immune mechanisms in disease progression. Simultaneously, the divergent activities across pathways between the low and high-risk cohorts underscore potential molecular signatures associated with distinct risk profiles.

Predictive analysis of the correlation between disulfidptosis-related lncRNA model and immune cell infiltration and tumor microenvironment

To probe the interconnection between the disulfidptosis-related lncRNA model and immune cell infiltration, along with the tumor microenvironment, a battery of algorithms, including QUANTISEQ, XCELL, TIMER, EPIC, CIBERSORT, MCPCOUNTER, and CIBERSORT-ABS, was employed. These analyses unveiled both affirmative and adverse correlations between risk scores and specific immune cell types (Figure 9A). CIBERSORT analysis additionally disclosed noteworthy distinctions in the abundance of plasma cells, Macrophages M2, T cells CD4 memory resting, Macrophages M0, T cells CD8, B cells naive, Tregs, T cells CD4 memory activated, Monocytes, and Mast cells resting between the high and low-risk cohorts (Figure 9B). Furthermore, disparities in immune functions, such as Type II IFN Response, were noted between the two risk groups (Figure 9C, 9D).

Figure 9. Analysis of immune cell infiltration and tumor microenvironment. (A) Bubble plots were employed for the visualization of immune cell infiltration within the tumor microenvironment. (B) Variances in immune response were evident through differential immune cell infiltration between the high and low-risk groups. (C, D) The high and low-risk groups exhibited significant disparities in immune function, implying modifications in immune activity associated with breast cancer risk. (E) Observable distinctions in tumor microenvironment scores between the high and low-risk groups indicated variations in the composition and characteristics of the tumor microenvironment. Statistical significance was denoted by * (p < 0.05), ** (p < 0.01), and *** (p < 0.001).

In addition, studies of the tumor microenvironment found that both stromal and immune cells were more abundant in the high-risk group than in the low-risk group, which may be related to their poor prognosis (Figure 9E).

While further validation is warranted, these findings suggest that the constructed disulfidptosis-related lncRNA model can offer valuable insights into the relationship between immune cell infiltration, risk scores, and immunotherapy pathways in BRCA.

Disparate analysis of drug sensitivity linked to long non-coding RNAs correlated with Disulfidoptosis-induced mortality

The investigation into disparities in drug sensitivity between cohorts classified as high and low risk represents a potential avenue for refining strategies in breast cancer therapeutics. Our scrutiny disclosed noteworthy distinctions in the IC50 values of 32 pharmaceutical agents among these risk-stratified groups, implying that risk scores could function as predictive benchmarks for discerning drug sensitivity and informing judicious clinical dosing decisions. Notably, BI-2536, Erlotinib, KU-55933, and MG-132 exhibited lower IC50 values in the high-risk group, implying enhanced efficacy of these drugs in treating high-risk breast cancer patients. Conversely, Afuresertib, BMS-345541, GNE-317, Ipatasertib, IWP-2, Leflunomide, Mitoxantrone, and MK-2206 demonstrated relatively higher IC50 values in the high-risk group, suggesting potentially reduced effectiveness and the necessity for higher dosages in high-risk patients. These findings propose the potential utility of risk scores in drug selection and dosage optimization for breast cancer treatment. However, it is crucial to note that further validation is necessary to confirm these results and evaluate the clinical applicability of risk scores in predicting drug sensitivity (Figure 10A–10L).

Figure 10. DAL-related drug sensitivity analysis. (A) Afuresertib, (B) BI-2536, (C) BMS-345541, (D) Erlotinib, (E) GNE-317, (F) Ipatasertib, (G) IWP-2, (H) KU-55933, (I) Leflunomide, (J) MG-132, (K) Mitoxantrone, (L) MK-2206.

Comparison of somatic mutations in the low and high-risk groups

We undertook an examination of somatic mutation frequencies in cohorts of breast cancer patients, aiming to discern variances between those deemed low and high risk. Our observations unveiled discernible patterns in the mutation frequencies of pivotal genes. In particular, the incidence of PIK3CA mutations was less frequent in the high-risk category (31%) in contrast to the low-risk subset (38%), while TP53 mutations exhibited a heightened prevalence in the high-risk group (42%) compared to their low-risk counterparts (21%). It is noteworthy that PIK3CA and TP53 emerged as the most frequently mutated genes among individuals with breast cancer (Figure 11A, 11B).

Figure 11. Comparative analysis of somatic mutation rates and influencing factors in high and low-risk groups. (A, B) Comprehensive assessment was undertaken to discern disparities in somatic mutations between the high and low-risk groups. (C) Observable distinctions in tumor mutation burden (TMB) were noted between the high and low-risk groups, signifying variations in the mutational landscape. (D) Investigation into the relationship between TMB and risk scores yielded insights into the connection between mutational burden and breast cancer risk. (E) Examination of the correlation between risk scores and cancer stem cells (CSC_s) provided illumination on the potential involvement of CSCs in breast cancer risk. (F) Comparative analysis of survival curves over time between high and low tumor mutation burdens elucidated the impact of mutational burden on patient survival. (G) Survival curves over time were scrutinized for different combinations of tumor mutation burdens and risk scores, emphasizing the synergistic effect of mutational burden and risk scores on patient prognosis.

Moreover, the cohort identified as high-risk demonstrated a markedly elevated tumor mutation load in comparison to their low-risk counterparts. Additionally, a discernible positive correlation emerged between the risk score and the tumor mutation load (R=0.18, P=2.5e-08) (Figure 11C, 11D). Tumor stem cells (DNAss and RNAss) are crucial regulators in the advancement of tumors. We also investigated the relationship between stem cell scores, derived from transcriptomic data, and patient risk, revealing a positive correlation between the two (P=1.8e-12, R=0.25) (Figure 11E).

Survival analysis based on tumor mutational load demonstrated a gradual decrease in patient survival over time. Individuals with a diminished mutation load exhibited a more favorable prognosis compared to those harboring elevated mutation burdens, whereas those with concurrent high-risk scores and heightened tumor mutation burden demonstrated the most unfavorable prognosis (P<0.01) (Figure 11F, 11G).