A total of three publicly available datasets derived from The Cancer Genome Atlas (TCGA), Molecular Taxonomy of Breast Cancer International Consortium (METABRIC) and Gene Expression Omnibus (GEO) databases served as sources of data for study cohorts designated A, B and C, respectively. A survival prognosis model was developed based on the 303 TNBC patients of cohort A, of whom some patients were assigned to the early-stage group (I~II, 83 patients) for training and the remaining patients were assigned to the late-stage group (IIA~IV, 220 patients) for internal validation (31). To validate the established model, cohort B gene expression data (derived from METABRIC) derived from 154 paracancerous samples and 209 BC samples were employed (32). Thereafter, the survival prognosis model was validated for different TNBC subtypes by performing virtual microdissection analysis of cohort C data (accession number GSE16446) (33, 34) followed by survival analysis performed using Kaplan-Meier estimation analysis.

First, bulk RNA-seq-based gene expression profiles of the 691-patient training cohort were subjected to virtual microanatomical analysis using NMF in R using one of the most critical NMF parameters, the decomposition level r, to define the number of clusters. When r was set to 4, the highest co-occurrence correlation coefficient was attained and the TCGA training cohort dataset was efficiently deconstructed ( Figure 1A ). Based on this result, r was set to 4. Following the method used in a recently reported study, immune and stromal enrichment scores were determined using single-sample gene set enrichment analysis (ssGSEA) (35, 36), which was bundled in GSVA to enable the identification of immune-related and stromal-related expression profiles (15). Once the immune and stromal cell enrichment scores were combined with the four NMF-identified clusters (numbered 1-4), we noticed that cluster 1 yielded stronger enrichment scores than the other clusters. As a result, cluster 1 was considered to be the “immune stromal cluster” within this context (37).

Second, numbers of certain types of immune cells present within immune stromal cluster tumors were studied by collecting data signals representing various immune cell subtypes while enrichment scores were determined according to ssGSEA-derived expression profiles. To evaluate the number of individual immune cells within the immune stromal clusters, enrichment fractions of 28 immune cells were merged with the clusters. Next, absolute fraction data of 22 infiltrating immune cells that were predicted by the CIBERSORT algorithm (according to gene expression patterns) were collected from the TIMER database (http://timer.cistrome.org/infltration estimation for tcga.csv.gz). Thereafter, leukocyte fraction data (TCGA all leuk estimate.masked.20170107.tsv) obtained from Thorsson et al. (https://gdc.cancer.gov/about-data/publications/panimmune) were used to determine DNA methylation-based signatures (38). Leukocyte fraction data were then obtained here based on analysis of images of TCGA tumors (including TNBC tumors) that can be found in the Supplementary Table ( Table S1 ) of the Saltz study (39). To confirm lymphocyte enrichment within the immune stromal clusters, the immune-stromal class data were compared to that of the remaining clusters.

Finally, multiple inhibitory receptor expression profiles were analyzed that revealed immune stromal clusters that overexpressed multiple inhibitory receptors, of which a large proportion of signatures associated with T cell depletion were found to be enriched. As a result, the patients within the immune stromal cluster was designated as the IDC, while the remaining population was designated as the resting class.

To confirm the presence of an immune-depleted state in advanced TNBC patients, the abovementioned methodologies for evaluating RNA-seq-based bulk gene expression profiles of the 619 early-stage TNBC samples were used to analyze corresponding profiles for the 210-sample late-stage cohort. Ultimately, four clusters were identified for the late-stage TNBC cohort (as for the early-stage TNBC cohort). Of these, cluster 2 yielded higher immune cell, stroma and TEX-related signature enrichment scores after signature enrichment scores from late-stage samples of the early TNBC cohort were incorporated while calculating the late-stage cohort signature enrichment score. As a result, cluster 2 was selected to be the IDC of the advanced TNBC group. After proportions of immunological cells and white blood cells and expression levels of various inhibitory receptors were compared between IDC and resting groups, gene set enrichment analysis (GSEA) was conducted to assess enrichment levels of markers and pathways then IDC scores were calculated based on ssGSEA of enrichment scores of 157 deficient immune function-related genes obtained during the training phase. The predictive potential of IDC scores was assessed using receiver operating characteristic (ROC) analysis.

Patient ICB treatment responses were predicted using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm. To investigate IDC patient responses to ICB treatment, programmed death-ligand 1 (PD-L1) expression was compared between IDC and resting classes. Higher TIDE prediction scores were generally associated with worse ICB responses.

According to Mariathasan et al., the cytokine transforming growth factor-β (TGF-β encoded by TGFB1) inhibits antitumor immunotherapies (40). By contrast, therapeutic co-administration of anti-TGF-β blockade-inducing antibody and anti-PD-L1 antibody reduced TGF-β signaling in stromal cells, promoted T cell infiltration into the tumor center and stimulated strong antitumor immunity and tumor regression. Taken together these observations suggest that TGF-β acts by limiting T cell infiltration of the TME to thereby suppress antitumor immunity. Therefore, detection of TGF-β in TNBC patient tumors correlates with resistance to antitumor immunotherapies.

The R language LIMMA package was used to identify differentially expressed genes between immune cell-deficient TNBC and other TNBC samples and univariate COX regression analysis was performed (41). Thereafter, least absolute shrinkage and selection operator (LASSO) regression analysis was conducted to identify significant genes associated with survival for use in building a risk prediction model. The formula for calculating the risk score was as follows:

where coefi is the coefficient and X is the normalized count for each gene. Based on the median risk score, we assigned patients to high-risk and low-risk survival-based groups. To test the model, follow-up TNBC data obtained from TCGA and metabric databases were used as training and test datasets, respectively. The reliability of the risk scoring model was assessed by within-group validation based on survival curves plotted based on training and test datasets then 1-, 3- and 5-year survival rates of patients were predicted based on the model. Thereafter, risk scores and clinical characteristics were evaluated together to obtain clinical information related to TNBC patient survival and prognosis followed by the creation of a forest plot of clinical prognostic factors.

Drug response and drug-targeting pathway information were obtained by searching the Cancer Drug Sensitivity Genomics (GDSC) website (https://www.cancerrxgene.org/). Next, sensitivity data of two TNBC cell lines to different drugs were obtained then the pRRophetic package in R language was used to predict drug sensitivities of different TNBC cell phenotypes as based on gene expression data (42). These results were then used to generate risk scores for different drugs as a basis for clinical drug selection.

Human breast cancer paracancerous cell line (MCF-10A) and human triple-negative breast cancer cell line (MDA-MB-453, MDA-MB-231) were purchased from Shangcheng North Na Chuanglian Biotechnology Co., LTD. MCF-10A was combined with MDA-MB-453 and MDA-MB-231 respectively constructed control groups to analyze the expression differences of PDCD1 and TGFB1 in different triple-negative breast cancer cell lines and normal breast cell lines. Total RNA was extracted using the Redzol kit from Beijing SBS Gene Technology Co., Ltd., and was extracted according to the instructions.

The forward primer sequence is TGFB1: F-5’-TTGACTTCCGCAAGGACCTC-3’, the reverse primer sequence is TGFB1: R-5’-ATCCGCAGTCCTCTCTCCAT-3’, and the product length is 421bp; PDCD1: F-5’-TGACTTCCACATGAGCGTGG-3’, the reverse primer sequence is PDCD1: R-5’-GCTCCTATTGTCCCTCGTGC-3’, and the product length is 294bp. qRT-PCR was performed using the SYBR^® Premix Ex Taq™ II (Takara, Shiga, Japan) Kit and a StepOnePlus Real-Time PCR instrument. Briefly, the mixture contained SureScript RTase Mix (20×) 1 μL, SureScript RT Reaction Buffer (5×) 4.0 μL, Total RNA 1 μg and dd HO (RNase/DNase free) supplemented to 20 μL, then sealed the cDNA with a transparent sealing film Array, mix well, spin off for 5s, 5 min at 25°C, 15 min at 42°C, 5 min at 85°C, keep at 4°C, and store the product at -20°C. qRT-PCR reactions were performed on a LightCycler 96 Real-Time PCR System (Roche Diagnostics, Indianapolis, IN). The reaction mixture was activated at 50°C for 2 min, pre-denatured at 95°C for 10 min, and then subjected to 40 cycles of amplification reactions at 95°C for 15 s and 60°C for 30 s. Finally, LightCycler 96 software (version 1.1.0.1320, Roche) was used for the collection and analysis of qRT-PCR data. With β-actin as an internal reference gene, the relative expression of mRNA was calculated by 2^-ΔΔCt method.

R (version 4.2.1, http://www.rproject.org) was used for statistical discrete analysis. The Wilcoxon rank sum test for continuous data was used to correlate IDC and the remaining categories with tumor-infiltrating lymphocyte (TIL) percentage, mutation number and neoantigen number. Overall survival (OS) data were analyzed using Kaplan-Meier estimates and log-rank testing. To discover variable combinations, we included all clinicopathological factors in the Cox model. P-values of ≤0.05 were considered statistically significant. Pearson’s correlation was used to assess the strengths of two-variable linear relationships. Maftools, which enables the visualization and analysis of somatic mutations and calculation of total somatic mutation numbers (43), was used to determine differences in numbers of mutations between the IDC and other classes. Genomic mutation data of all 303 TNBC tumor samples were obtained from the TCGA database.

NMF analysis was performed on the large number of RNA-seq-based gene expression profiles that were obtained for 691 TNBC samples in the training cohort then TEX-associated TME transcriptome signals were extracted. The training cohort dataset was clearly divisible into four patient clusters ( Figure 1A ). TNBC patients in cluster 1 had high immune and stromal enrichment fractions as assessed by ssGSEA and batch RNA-seq-based gene expression profiling, indicating considerable enrichment of immune cell and stromal feature gene expression signatures. This cluster, hereafter referred to as the immunological stromal cluster ( Figures 1B, C ), was subsequently found to possess a high density of immune cell signatures ( Figure 1D ) associated with immune cell subsets such as M2 macrophages, B cells, CD8⁺ T cells and eosinophils. Furthermore, the immune stromal cluster possessed a considerably greater proportion of leukocytes, as estimated from DNA methylation data, than did the other clusters ( Figure 1E ). Numbers of TILs detected during pathological examinations of images of tumors were also considerably higher in the immune stromal class than in the other clusters ( Figure 1F ). After the enrichment of immune cells in the immune stromal class was confirmed, absolute proportions of immune cells in the immunological stromal class and remaining clusters were compared using CIBERSORT analysis of batch RNA-seq data ( Figure 1G ).

To investigate TEX signaling function in TNBC, expression levels of several inhibitory receptors of immune system receptors, including CTLA-4, PDCD1 (also known as PD-1), BTLA, LAG3, TIGIT, HAVCR2 (also known as TIM-3), SIGLEC7, IDO1 and VISTA, were measured. Based on results of inhibitory receptor expression analyzes and enrichment scores for the TEX signaling gene set, a novel IDC subpopulation of the immune stromal cluster was identified that accounted for 21.4% of the population of the training cohort ( Figure 1H ). The resting class was defined as the resting subpopulation of the training cohort.

To define IDC-associated molecular markers, GSEA based on gene expression profiles of the training cohort was conducted that led to the identification of six marker gene sets that were enriched in the IDC molecular signature, with especially high enrichment noted for immunocytokine-related pathways and markers. Previously reported research results indicated that exhausted T cells can undergo apoptosis, as consistent with results obtained here showing that apoptotic markers were considerably elevated in the IDC of the early-stage TNBC cohort as evidence of severe T cell depletion. Results obtained using NMF analysis of IDC gene expression patterns led to the identification of 273 metagenes (sets of genes with coordinated expression) that after ClueGO network analysis were grouped into Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway networks that revealed their biological functions ( Figure 1I ). The majority of genes within these metagenes were functionally linked to cytokines and associated molecules that mediate immunological and inflammatory responses. These 76 immune function-related genes, which included members of the immunoglobulin superfamily (CDH9 and CDHR3), interleukin receptor molecules (IL9, IL19, IL13RA1 and IL12RB2) and cytokines (CD70 and CD276) were then utilized to validate TNBC patient immunological exhaustion ( Supplementary Tables S2 , S3 ). Taken together, these results demonstrated that we effectively identified an IDC gene expression pattern indicating severe T cell depletion within the TME of some TBNC patients.

We conducted NMF analysis of bulk RNA-seq-based gene expression profiles of TCGA late-stage samples to confirm the presence of reduced immune classes in the TNBC TME, with results of NMA leading to identification of four clusters ( Figure 2A ). When these clusters were combined with ssGSEA-calculated feature enrichment scores and bulk RNA-seq-based gene expression profiles, cluster 1 was found to possess higher stromal and immune enrichment scores than the other clusters ( Figure 2D ). Along with other TEX-related cytokine indicators, these features were considerably enriched ( Figures 2B, C ). As a result, cluster 1 was designated as the IDC of the TGCA samples, as consistent with observed co-upregulation within this cluster of numerous inhibitory receptors ( Figure 2E ) that were enriched for Gene Ontology (GO) functional terms leukocyte migration, positive regulation of the MAPK cascade, cell chemotaxis, positive regulation of cell adhesion, cellular divalent inorganic cation homeostasis, positive regulation of cell adhesion, calcium ion homeostasis, cytokine-mediated signaling pathway, cellular calcium ion homeostasis and positive regulation of cytokine production ( Figures 2D, E ). Notably, ROC curve-based results ( Figure 2F ) confirmed that the 167 depleted immune function-related genes may be useful for predicting the presence of IDC immunosuppressive signatures in TMEs of individual TNBC patients.

In order to explore the utility of reduced immune function status for predicting TNBC patient prognosis, we correlated categories with clinicopathological factors. Previous studies had linked high-density TILs with better outcomes, such as improved overall survival (OS) (44, 45). Here, IDCs of both early and advanced TNBC patients possessed higher proportions of TILs than did IDCs of other patient categories. However, Kaplan-Meier estimates of early and advanced TNBC OS rates revealed significantly lower OS rates in patients with TME IDC immunosuppressive signatures than in TMEs of the resting group (P< 0.0001; Figures 3A, B ), while in both early and advanced TNBC, multivariate survival analysis conducted using Cox regression models showed that the IDC immunosuppressive signature was an independent predictor of OS (P< 0.001). ( Figures 3E, F ). Finally, the prognostic value of IDC immunosuppressive signatures in TCGA and GEO cohorts of TNBC patients in all disease stages were studied. As expected, IDC patients had worse OS than the resting groups ( Figures 3C, D ). Taken together, these findings suggest that although T cells are abundant in IDCs, most of these T cells have lost effector functions for limiting tumor growth, resulting in a state of immunosuppression that leads to continued tumor growth and poor prognosis. Survival data also confirmed that IDC signatures of patients with late-stage TNBC exhibited more pronounced T cell depletion than those of early-stage TNBC patients.

Treatment of cancer patients with therapeutic antibodies targeting the PD-L1 pathway can elicit long-lasting, robust responses. However, efficacies of such treatments are frequently reduced due to the emergence of drug resistance. To study the response of the IDC group to ICB treatment, we compared PD-L1 expression levels between IDC and resting patient groups and found that PD-L1 expression levels were higher in the early IDC group than in the resting group ( Figures 4A, B ). Next, use of the Tumor Immune Dysfunction and Rejection (TIDE) algorithm to predict ICB treatment responses revealed that both early and advanced TNBC IDCs had higher TIDE prediction scores than did the other categories ( Figures 4C, D ), such that a higher TIDE prediction score indicated a worse ICB response. Importantly, TGFβ1 was expressed at a lower level in resting classes than in IDCs ( Figures 4E, F ), which is consistent with results obtained by Mariathasan et al. that indicated the cytokine TGF-β (encoded by TGFβ1) suppresses effects of antitumor immunotherapy (46). Taken together, these results suggest that the IDC immunosuppressive signature has a ICB therapy resistance.

Considering that tumors that escape immune editing in cancer patients often express molecules that suppress the antitumor immune response, such as PD-L1 (also known as B7-H1), indoleamine 2,3-dioxygenase (IDO) and others, here expression levels of these molecules were correlated with TGF-β expression using Pearson correlation analysis. As expected, expression levels of these molecules were positively correlated with TGFβ1 mRNA expression levels ( Figures 5A–I ), thus suggesting that IDC patients remain in an ICB-unresponsive state, despite higher PD-L1 expression levels and that anti-PD-L1 therapies. In addition, these results suggest that in vivo antitumor response molecules play key roles in this process.

Due to the fact that genomic mutations occurring in tumors are strongly associated with immunotherapy outcomes, we explored the landscape of commonly mutated genes between TNBC IDC and resting class molecular signatures. In the TCGA cohort, censorship mutations were the most common TNBC mutations detected ( Figure 6A ), of which most mutations were SNPs ( Figure 6B ) that most often contained C>T base substitutions ( Figure 6C ). The ten genes with the highest mutation frequencies were TP53, USH2A, AHNAK, PIK3CA, NOTCH1, MUC16, DNAH5, DNAH11, AHNAK2 and AKAP9 ( Figure 6D ), with no significant differences in individual gene mutations observed between IDC and resting group signatures ( Figure 6E ). These results thus suggest that somatic mutations are not significantly associated with immunosuppressive TME status.

Ultimately, a total of 211 differentially expressed genes associated with survival were identified ( Figure 7A , Supplementary Table S4 ). The univariate Cox regression algorithm was used to initially obtain 150 genes related to TNBC patient prognosis then hazard ratios (HRs) and P values of these genes were calculated ( Supplementary Table S5 ). Next, we constructed risk models using the LASSO algorithm that ultimately led to the identification of 45 prognostically relevant genes ( Figure 7B ). These genes were then used to construct risk score-based models based on training (n = 303) and test (n = 211) datasets obtained from TCGA-TNBC and metabric TNBC patients, respectively. Results of survival analysis revealed that higher risk scores in training and test sets corresponded to poorer survival (P< 0.0001) ( Figures 7C–F ). Time-dependent ROC curves were generated and used to assess the sensitivity of the prognostic model. Results obtained for the areas under ROC curves (AUCs) revealed 3-, 5- and 10-year AUCs for the training set of 0.768, 0.794 and 0.746, respectively ( Figure 7G ), and corresponding AUCs for the test set of 0.599, 0.656 and 0.689, respectively ( Figure 7H ). Additionally, multivariate Cox regression analysis was conducted to assess whether clinical characteristics, such as age, cancer stage, TNM stage and risk scores were independent factors related to TNBC patient prognosis ( Supplementary Figure S2 ). We found that age, tumor stage and risk score were independent prognostic factors for both training and test sets of TNBC patients. Notably, more high-risk patient molecular signatures contained IDC signatures ( Figure 7I ), suggesting that immune system failure may be responsible for poor survival of high-risk patients.

The R package “pRRophetic” was used evaluate high-risk and low risk patient IDC, and resting class groups for chemotherapeutic responses and resistance. Figure 8 shows drug sensitivity results for three types of triadic breast cancer cell lines to eight anti-cancer therapy drugs (sorafenib, gefitinib, bleomycin, bosutinib, etoposide, lenalidomide, camptothecin, methotrexate). Statistical analysis showed that, except bleomycin, the IC50 level of resting grade patients was higher than that of IDC patients. Sensitive ( Figures 8A–H ). In addition, for the risk group Sorafenib, the low-risk group had a higher score, while for Bosutinib and Camptothecin, the low-risk group had a lower score ( Figures 8I–P ).

Expression levels of selected genes were determined using qRT-PCR to confirm their reliability. As shown in Figures 9A, B , the results for all mRNA transcripts showed that they were expressed at significantly higher levels in the two TNBC tumor cell lines (MDA-MB-231, MDA-MB-453) than in the adjacent normal cell lines.