Identification of immunosuppressive signature subtypes and prognostic risk signatures in triple-negative breast cancer

Abstract

Purpose: Immune checkpoint blockade (ICB) therapy has transformed the treatment of triple-negative breast cancer (TNBC) in recent years. However, some TNBC patients with high programmed death-ligand 1 (PD-L1) expression levels develop immune checkpoint resistance. Hence, there is an urgent need to characterize the immunosuppressive tumor microenvironment and identify biomarkers to construct prognostic models of patient survival outcomes in order to understand biological mechanisms operating within the tumor microenvironment.

Patients and methods: RNA sequence (RNA-seq) data from 303 TNBC samples were analyzed using an unsupervised cluster analysis approach to reveal distinctive cellular gene expression patterns within the TNBC tumor microenvironment (TME). A panel of T cell exhaustion signatures, immunosuppressive cell subtypes and clinical features were correlated with the immunotherapeutic response, as assessed according to gene expression patterns. The test dataset was then used to confirm the occurrence of immune depletion status and prognostic features and to formulate clinical treatment recommendations. Concurrently, a reliable risk prediction model and clinical treatment strategy were proposed based on TME immunosuppressive signature differences between TNBC patients with good versus poor survival status and other clinical prognostic factors.

Results: Significantly enriched TNBC microenvironment T cell depletion signatures were detected in the analyzed RNA-seq data. A high proportion of certain immunosuppressive cell subtypes, 9 inhibitory checkpoints and enhanced anti-inflammatory cytokine expression profiles were noted in 21.4% of TNBC patients that led to the designation of this group of immunosuppressed patients as the immune depletion class (IDC). Although IDC group TNBC samples contained tumor-infiltrating lymphocytes present at high densities, IDC patient prognosis was poor. Notably, PD-L1 expression was relatively elevated in IDC patients that indicated their cancers were resistant to ICB treatment. Based on these findings, a set of gene expression signatures predicting IDC group PD-L1 resistance was identified then used to develop risk models for use in predicting clinical therapeutic outcomes.

Conclusion: A novel TNBC immunosuppressive tumor microenvironment subtype associated with strong PD-L1 expression and possible resistance to ICB treatment was identified. This comprehensive gene expression pattern may provide fresh insights into drug resistance mechanisms for use in optimizing immunotherapeutic approaches for TNBC patients.

Keywords: TNBC; immunogenomics; immunosuppressive cytokines; t cell exhaustion; tumor microenvironment.

Figure 1

Identification and molecular characterization of IDC. (A) Heatmap of gene expression clusters and distinct expression patterns of 691 early-stage (III) TNBC samples from unsupervised NMF. (B) Matrix and immune enrichment scores for clusters of four expression patterns. High and low gene enrichment scores are depicted in purple and orange, respectively. (C) Enrichment fractions of gene signatures identify immune stroma and other clusters of immune cells. (D) CIBERSORT-inferred absolute fraction of TME cells compared between the two classes. (E) Box plot showing the difference in Leukocyte fraction between the two classes. (F) Box plots showing differences in the expression levels of various inhibitory receptors in rest clusters and in the immune and stromal cluster. (G) Histogram of the percentage of immune cells in each sample. (H) GSEA analysis reveals that IDC shows significant enrichment of marker gene sets related to immune cell metabolic processes. (I) KEGG pathway functional grouping network by ClueGO/CluePedia. Colorless and colored nodes represent metagene-specific genes and KEGG pathway terms, respectively. Node colors represent different functional groups. Node size represents the importance of KEGG pathways. The more important the KEGG pathway, the larger the highlighted node. All statistical differences between the two categories were compared using the Wilcoxon rank sum test. ns, >0.05; *P< 0.05; 0.01; ***P< 0.001.

Figure 2

Internal validation of IDC in late stage TCGA-TNBC samples. (A, B) Consensus clustering for the late stage TCGA-TNBC. (C) The comparison of the absolute fractions of TME cells inferred by CIBERSORT between two classes. (D) Kegg pathway enrichment visualization in immune depletion Classifer Genes. (E) Boxplots shows the different expression levels of multiple inhibitory receptors between two classes. (F) ROC curve evaluated the predictive capacity of 157 immune depletion classifer genes in late stage TCGA-TNBC samples. (G) Box plots show the differences of leukocyte fraction between two classes. All statistical differences of two groups were computed by Wilcoxon rank-sum test; *P< 0.05; **P< 0.01; ***P< 0.001.

Figure 3

Prognostic analysis of different stages of IDC and rest types of TNBC. (A, C) Kaplan-Meier estimates of overall survival for IDC and rest classes of advanced, early, and persistent TNBC. (D) Kaplan-Meier estimates of overall survival for IDC and rest classes of TNBC in the gse16446 dataset. P-values were calculated using the log-rank test. (E, F) Multivariate and univariate Cox regression analysis (group, tumor stage, and age) for full-stage TNBC.

Figure 4

Prediction of resistance to ICB therapy. (A, B) Different expressions of PD-L1 expression levels in different stages of IDC patients. (C, D) Different expression of TIDE predictive score for each stage of ICB treatment. (E, F) Different expression of TGFb1 in IDC and rest classes.*: *P< 0.05; ***P< 0.001.

Figure 5

Expression correlation of TGFβ1 and immune checkpoint-related genes in TNBC patient samples. (A) CTLA4, (B) PDCD1, (C) LAG3, (D) BTLA, (E) TIGIT, (F) HAVCR2, (G) IDO1, (H) SIGLEC7, (I) VISTA.

Figure 6

Association of IDC with somatic mutations. (A, B) Summary of mutation information for TNBC samples in the TCGA database. (C) Summary of variant Classification within TNBC. (D) Top 10 triple negative breast cancers mutated genes. (E) The landscape of most frequently mutated genes between the IDC and the rest class in TNBC.

Figure 7

Construction and validation of the risk score. (A) Volcano plot of differentially expressed genes between IDC and rest class samples, blue indicates down-regulated expression and red indicates up-regulated expression. (B) Lasso regression analysis and multivariate stepwise Cox regression analysis for identification of the immune risk signature. (C, E) Kaplan-Meier curves of training set (p<.001, log-rank test) and test set (p<.001, survival rate comparison). (D, F) Association between patient survival and increased risk score. (G, H) Time-dependent receiver operating characteristic (ROC) of training and test sets. (I) Proportion of IDC and rest class in high and low risk patients.

Figure 8

(A–H) The violin plot shows that drug sensitivity prediction score(Sorafenib, Gefitinib, Bleomycin, Bosutinib, Etoposide, Lenalidomide, Camptothecin, Methotrexate) are distributed differently among groups. (I–P) The violin plot shows that drug sensitivity prediction score(Sorafenib, Gefitinib, Bleomycin, Bosutinib, Etoposide, Lenalidomide, Camptothecin, Methotrexate) are distributed differently among risk groups; ns: no significance; *: P < 0.05; ***: P < 0.001.

Figure 9

The mRNA expression of TGFB1 (A) and PDCD1 (B) in a TNBC cell lines and the adjacent cell lines. ***P< 0.001.