Development of a CD8+ T cell-based molecular classification for predicting prognosis and heterogeneity in triple-negative breast cancer by integrated analysis of single-cell and bulk RNA-sequencing

Abstract

Background: Triple-negative breast cancer (TNBC), although the most intractable subtype, is characterized by abundant immunogenicity, which enhances responsiveness to immunotherapeutic measures.

Methods: First, we identified CD8⁺ T cell core genes (TRCG) based on single-cell sequence and traditional transcriptome sequencing and then used this data to develop a first-of-its-kind classification system based on CD8⁺ T cells in patients with TNBC. Next, TRCG-related patterns were systematically analyzed, and their correlation with genomic features, immune activity (microenvironment associated with immune infiltration), and clinicopathological characteristics were assessed in the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC), the Cancer Genome Atlas (TCGA), GSE103091, GSE96058 databases. Additionally, a CD8⁺ T cell-related prognostic signature (TRPS) was developed to quantify a patient-specific TRCG pattern. What's more, the genes-related TRPS was validated by polymerase chain reaction (PCR) experiment.

Results: This study, for the first time, distinguished two subsets in patients with TNBC based on the TRCG. The immune microenvironment and prognostic stratification between these have distinct heterogeneity. Furthermore, this study constructed a novel scoring system named TRPS, which we show to be a robust prognostic marker for TNBC that is related to the intensity of immune infiltration and immunotherapy. Moreover, the levels of genes related the TRPS were validated by quantitative Real-Time PCR.

Conclusions: Consequently, this study unraveled an association between the TRCG and the tumor microenvironment in TNBC. TRPS model represents an effective tool for survival prediction and treatment guidance in TNBC that can also help identify individual variations in TME and stratify patients who are sensitive to anticancer immunotherapy.

Keywords: CD8+ T cell; Heterogeneity; Triple-negative breast cancer; Tumor immune microenvironment; single cell sequencing.

Fig. 1

Schematic summary of the workflow.

Fig. 2

Single cell RNA data analysis; Umap diagram of 16 samples(A); umap distribution diagram of annotation cell subgroups(B); umap of four T cell subpopulations after cluster analysis (C); Analysis of the weighted co-expression network in metabric-TNBC cohort. Sample clustering of dataset metabric-TNBC(D). Volcano Plot illustrated the DEGs between CD8⁺ T cells and no CD8⁺ T cells based on single cell datasets, the genes of CD8⁺ T cell related classifier are exhibited(E). Analysis of correlations between modules and CD8⁺ T cell characteristic(F). Different modules are produced and shown in different colors by aggregating genes with strong correlations into a same module (G). Identification of optimal thresholds, which is 3(H). Scatter plot of module eigengenes in the turquoise module (I). Heatmap describing the topological overlap matrix among genes based on co-expression modules (J).

Fig. 3

Patients stratified into two clusters based on the CD8⁺ T cells hub genes. Consensus matrices of patients in the metabric-TNBC(A) and merge-TNBC(C) cohort via the unsupervised consensus clustering method (K-means). Survival analysis of the heterogeneous clusters in the metabric-TNBC (B) and merge-TNBC (D) cohort. Analysis of the hub genes–immune response relationships of TNBC in METABRIC data(E). Heatmap of the clinicopathological manifestations among the TRCG-related patterns (F).

Fig. 4

Comparison of the two subtypes in TNBC patients. Expression of chemokines, receptors and MHC molecules(A) and Immune cell infiltration levels Immune cell infiltration levels(B). Evaluation of the TME in the two TTK related patterns(C). Boxplots depicting the difference of the TME related signatures(D), the cancer–immunity cycle(E) and TIS scores(F) amongs the two TRCG related patterns via ssGSEA. Response prediction to immunotherapy (anti-PD-1 and anti-CTLA4) amongs the two TRCG patterns based on SubMap algorithms(G).

Fig. 5

TRSP in metabric-TNBC cohort analysis. Kaplan–Meier curves (A), time-dependent ROC analysis (B) and risk score(C). The association between TRSP and special clinicopathological traits(D). Kaplan–Meier curves of TRSP related genes in METABRIC-TNBC and merge - TNBC cohort (E, F, G).

Fig. 6

The comparison of the TRSP -high and -low subgroups. Evaluation of the TME in the TRSP -high and -low subgroups(A). The heterogeneity of Immune cell infiltration levels amongs TRSP -high and -low subgroups in METABRIC(B) cohort by MCPcounter. Boxplots depicting the difference of the TME related signatures(C) amongs the two TRCG related patterns via ssGSEA. Correlation between TRSP and immune checkpoint-related genes(D). A comparison of the relative sensitivity of responding to anti-PD-1/PD-L1 as well as anti-CTLA-4 treatment in the TRSP high and low subgroups(E). The distinction of TIS scores amongs TRSP -high and -low subgroups(F). Response prediction to immunotherapy (anti-PD-1 and anti-CTLA4) amongs the TRSP – high and – low subgroups based on TIDE(G) and SubMap algorithms(H). The TRSP related genes are all low expressed in TNBC cells compared to normal breast cancer cell by Student's t-test. 2−ΔΔCt is used to present the fold change in qRT-PCR experiment(I). (*P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001).

Fig. 7

The exploration of TRSP relaed genes. UMPA (A) plots showing different TRSP related genes expression distribution based on scRNA-Seq and Violin plots revealing the difference of the expression of N4BP2L1, IL18RAP, KIR3DL3 in all cell subgroup (B) and T cell subgroup (C), respectively. Trajectory reconstruction of all T cells in TNBC, with a color code for pseudo-time (D), clusters (E), cell subtypes (F), respectively. The branched heatmap indicates the dynamics of the expression of N4BP2L1, IL18RAP, KIR3DL3 during T cells transdifferentiation(G), the redder the color, the higher the expression. The correlation between N4BP2L1, IL18RAP, KIR3DL3 and the level of immune cell infiltration (H).

figs1

figs2

figs3

figs4

figs5

figs6

figs7

figs8