Development of a prognostic model based on four genes related to exhausted CD8+ T cell in triple-negative breast cancer patients: a comprehensive analysis integrating scRNA-seq and bulk RNA-seq

Abstract

Low immune infiltration is closely associated with poor clinical results and an unfavorable response to therapy in triple-negative breast cancer (TNBC). T-cell exhaustion (TEX) is a significant risk factor for tumor immunosuppression and invasion. Although improving TEX and enhancing effector function are promising strategies for strengthening immunotherapy, their role in the pathogenesis of TNBC remains unclear. This study's objective was to develop a prognostic model for TNBC based on exhausted CD8+ T-cell (CD8+ Tex)-related differentially expressed genes (DEGs) and to investigate its clinical and immune relevance. Initially, 398 CD8+ Tex-related genes were screened utilizing single-cell RNA sequencing (scRNA-seq) data from TNBC patients. Pseudotime analysis confirmed that CD8+ Tex mainly clustered at the end of the differentiation pathways, making them a critical subset in TNBC progression. By analyzing the TCGA cohort, ten CD8+ Tex-related DEGs were identified as significantly correlated with overall survival (OS) in TNBC patients, and a prognostic model containing four biomarkers (GBP1, CTSD, ABHD14B, and HLA-A) was constructed. The model demonstrated robust predictive capability in both the TCGA cohort and an external cohort, with the low-risk group exhibiting elevated expression of immunological checkpoint molecules and immune cell infiltration, as well as better responses to immunotherapy and chemotherapy. Furthermore, these four biomarkers were found to be highly expressed on CD8+ Tex and were associated with cellular communication efficiency. Therefore, this model is expected to be a new method for forecasting TNBC patients' prognosis and effectiveness of treatment, providing new insights for clinical decision-making.

Keywords: Cellular communication; Exhausted CD8+ T cell; Single-cell RNA sequencing; Triple-negative breast cancer; Tumor immune microenvironment.

Fig. 1

Single-cell data analysis. A Post-batch correction clustering colored plot of different samples. B UMAP visualization colored to distinguish between various clusters. C UMAP plot colored to represent diverse cell types. D Bubble chart illustrating the expression patterns of marker genes associated with distinct cell types in TNBC. E Bar graph depicting the percentage of various cell types present within each individual sample. F UMAP visualization colored based on clusters of CD8+ T cells. G UMAP plot colored according to the CD8+ T cell types. H Bubble plot showing marker genes for different CD8+ T cell types. I Plot depicting the trajectory of CD8+ T cells across various states. J Visualization illustrating the differentiation pathways of CD8+ T cells over time. K Overall trajectory analysis of CD8+ T cell subgroups. L Hierarchical trajectory plot outlining the subgroups of CD8+ T cells

Fig. 2

Prognostic model building. A Volcano plot representing the results of differential analysis comparing TNBC tumor tissue with normal tissue. B Intersection of CD8+ Tex marker genes from scRNA-seq data and DEGs from the TCGA cohort. C Analysis of GO enrichment for the overlapping genes. D Enrichment analysis of KEGG pathways for the overlapping genes. E Univariate Cox regression analysis of intersection genes. F and G LASSO regression analysis. H Multivariate Cox stepwise regression analysis

Fig. 3

Prognostic model validation. A Distribution of risk score, survival status, and the expression profile of model genes for high-risk and low-risk patient groups in the training set. B KM curves for survival analysis of patients in the training set. C ROC curves of patients in the training set at 1, 3, and 5 years. D Distribution of risk score, survival status, and the expression profile of model genes for high-risk and low-risk patient groups in the validation set. E KM curves for survival analysis of patients in the validation set. F ROC curves of patients in the validation set at 1, 3, and 5 years

Fig. 4

Prognostic model accuracy verification and pathway analysis. A Univariate Cox regression analysis of risk score and clinical characteristics. B Multivariate Cox regression analysis of risk score and clinical characteristics. C Nomogram constructed to predict 1, 3, and 5-year OS. D Calibration curves for the nomogram. E DCA for the nomogram and clinical characteristics. F Donut plot showing the distribution of patients at different clinical stages. G and H GSEA pathway analysis results

Fig. 5

Immunoinfiltration analysis. A Box plot showing the correlation between high- and low-risk groups with ESTIMATEScore, StromaScore, TumorPurity, and ImmuneScore. B Box plot of immune cell enrichment differences. C Heatmap of multiple algorithms to assess the abundance of immune cell infiltration. (ns p ≥ 0.05, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001)

Fig. 6

Immunotherapy and chemotherapy response analysis. A Violin plot of IPS scores. B Violin plot of TIDE scores. C Violin plot of immune checkpoint expression. D Violin plot showing differences in sensitivity to 12 clinically common chemotherapy drugs (ns p ≥ 0.05, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001)

Fig. 7

Single-cell level validation. A Feature plot and violin plot showing the expression levels of four prognostic genes at the single-cell level. B Bar plot and feature plot based on cell scoring using the prognostic model gene set. C Bar plot showing differences in the number (left) and strength (right) of cellular communications between high- and low-exhaustion groups. D Network plot showing differences in the number (left) and strength (right) of cellular communications between CD8+ T cell subsets (Blue represents the low-exhaustion group with stronger communication abilities, and red represents the high-exhaustion group with stronger communication abilities). E Heatmap showing differences in the number (left) and strength (right) of cellular communications between CD8+ T cell subsets. F Proportional plot (left) and numerical comparison bar plot (right) showing differences in enriched signaling pathways between high- and low-exhaustion groups. G Heatmap showing the overall signaling pathway levels in the low-exhaustion group (left) and high-exhaustion group (right). H Bubble plot showing differences in total ligand-receptor pair communication probabilities between high- and low-exhaustion groups