The raw data were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) and Cancer Genome Atlas (TCGA) databases. Three TNBC datasets [GSE96058 (Brueffer et al., 2018), GSE86166 (Prabhakaran et al., 2017), and GSE103091 (Jezequel et al., 2015)] and two datasets related to immunotherapy [GSE35640 (Ulloa-Montoya et al., 2013) and GSE78220 (Hugo et al., 2016)] in the GEO database were used for analysis. The pan-cancer data involving 17 cancer types in TCGA were downloaded from the UCSC XENA database (https://xenabrowser.net/datapages/) (Goldman et al., 2020). We extracted TNBC data from all TCGA datasets for the principal analysis, and other tumors were used for validation. Moreover, the profiles of the IMvigor210 cohort were obtained according to official guidelines (http://research-pub.Gene.com/imvigor210corebiologies) (Mariathasan et al., 2018). All of the information of the public datasets is summarized in Supplementary Table S1.

In addition, we used a cohort constructed by the West China Hospital breast cancer specialist research team as an external validation cohort, including 80TNBC biopsies, and this experiment was approved by the Ethics Committee of West China Hospital. Total RNA was extracted and purified following the manufacturer’s protocol. After synthesizing first- and second-strand cDNA using random hexamer primers, DNA polymerase I and RNase H, the library fragments were purified with an AMPure XP system (Beckman Coulter, Beverly, MA, United States) as described in the NEBNext UltraTM Directional RNA Library following the manufacturer’s recommendations. The libraries were then sequenced on the Illumina HiSeq X ten platform (Novogene Bioinformatic Technology Co., Ltd., China) following a 150 bp paired-end read protocol. Eventually, the raw sequencing data from this study have been deposited in the Genome Sequence Archive (GSA) in BIG Data Center (https://bigd.big.ac.cn/) (Zhang et al., 2021a), Beijing Institute of Genomics (BIG), Chinese Academy of Sciences, under the accession number HRA002256.

Single-sample gene set enrichment analysis (ssGSEA) is a well-known method to derive the absolute enrichment scores of previously experimentally validated gene signatures conducted by the R package “GSVA,” a nonparametric and unsupervised method commonly employed to estimate the variations in the pathway and biological process activity of a single sample (Li et al., 2017). Here, we preferred to use ssGSEA to assess the relative abundance of immune cell infiltration levels in a single sample. Two validated immune cell signatures published, labeled immune cell signatures 1 and 2 in this study, were used in this research, containing 24 (Bindea et al., 2013) and 23 (Charoentong et al., 2017) types of immune cells, respectively. The markers of these two signatures are listed in Supplementary Tables S2, S3. To further validate the results from ssGSEA, the CIBERSORT algorithm (Newman et al., 2015), which is a deconvolution algorithm, was employed to infer cell-type proportions with bulk tumor sequence data. Moreover, the third method, called Estimation of Stromal and Immune cells in malignant tumors using Expression data (ESTIMATE) (Yoshihara et al., 2013), was also used to infer the fraction of stromal and immune cells in tumor samples.

Using ssGSEA described previously, GSVA (Hanzelmann et al., 2013) was used to assess pathway activation levels in a single sample with the gene set “c5.all.v6.2. symbols” downloaded from the MSigDB database in GSEA website (Mootha et al., 2003) and another published pathway gene set summarized in Supplementary Table S4 (Mariathasan et al., 2018). GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted using the R package and the online website Database for Annotation, Visualization, and Integrated Discovery (DAVID) (david.ncifcrf.gov) (Dennis et al., 2003).

Unsupervised clustering analysis was used to classify patients based on the immune cell infiltration levels with the ConsensuClusterPlus package (Wilkerson and Hayes, 2010). Differentially expressed gene (DEGs) analysis was conducted by the “limma” R package, with the criterion of adjusted p value < 0.05. The differentially expressed mRNAs were visualized by the “pheatmap” package.

A total of 113 ferroptosis regulators were extracted from the online website FerrDb (http://www.zhounan.org/ferrdb/), and the specific information of these genes is shown in Supplementary Table S5. To describe the ferroptosis level, the ferroptosis index (FPI) was established based on the expression data of genes in ferroptosis, including positive and negative components. The enrichment score (ES) was calculated using ssGSEA, and the FPI to roughly assess ferroptosis trends was calculated as follows:

FPI = ES (positive) − ES (negative) (Liu et al., 2020)

To assess the stemness of cancer cells, a one-class logistic regression algorithm known as mRNA-based stemness index (mRNAsi) was used to calculate the stemness index for each sample under the direction of the workflow available on a previously established database (https://bioinformaticsfmrp.github.io/) (Malta et al., 2018).

The chemotherapeutic response for TNBC was predicted according to the data involved in the Genomics of Drug Sensitivity in Cancer (GDSC) with the “pRRophetic” package (Geeleher et al., 2014). The Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard.edu/) was employed to predict the immunotherapy response of TNBC (Jiang et al., 2018), and the default cutoff value was 0.

The overlapping DEGs among the four TNBC datasets were regarded as TME gene signatures. Principal component analysis (PCA) was used to calculate the TRG score to quantize the TME patterns in TNBC. We summed PC1 and PC2 of genes i by PCA as described before by us (Liu et al., 2021). The TRG score was calculated as follows:

TRG score = Σ (PC1i + PC2i)

Logistic least absolute shrinkage and selection operator (LASSO) regression analysis can construct a prognostic signature to minimize the risk of overfitting (Simon et al., 2011). However, LASSO relies heavily on seeds when it allows. Iteration LASSO was independent of the seed once the roots, the optimal lambda, and the resulting feature were changed (Sveen et al., 2012). The features retained at high frequency can be considered the most influential factors. Genes included under consensus were generated by iteration of LASSO, and AUC further selected the minimum combination of genes associated with survival. The formula of patients’ risk scores was established:

Risk score = Σ (each gene’s expression × corresponding coefficient).

Receiver operating characteristic (ROC) curves and survival curves with the Kaplan–Meier method were used to judge the prediction efficiency of the signature. The best cutoff value of genes in survival analysis was searched by the “survminer” R package. The signature genes obtained from iterated LASSO analysis were used for nomogram construction using logistic and Cox regression analyses. Calibration curves were used to assess the predictive accuracy of the nomogram.

Correlation coefficients and p values were calculated by Spearman correlation analysis among several defined groups. Wilcoxon tests were used to compare differences between the two groups. The asterisks represent the statistical p values (*p < 0.05, **p < 0.01, and ***p < 0.001) in the panels.

The flowchart of this study is depicted in Figure 1. To explore the tumor microenvironment patterns in four independent TNBC cohorts, consensus cluster analysis was used to classify patients with TME conditions (Supplementary Figures S1A–D). By integrating the clustering results of each dataset, two distinct TME subtypes were eventually identified using unsupervised clustering in each cohort, labeled as subtypes 1 and 2 (Figure 2A). Here, we used immune cell signature 1 to perform cluster analysis. At the same time, we found that the infiltration of the levels of immune cells was significantly different using immune cell signature 2 in all four cohorts (Figure 2B). Among them, subtype 1 was enriched with immune cells compared with subtype 2, meaning that subtype 1 was an immune-activating subtype with higher immune cell infiltration levels, the same as the conception of a “hot” tumor. By the CIBERSORT algorithm, we found that some antitumor immune cells, such as CD8+ T cells, activated CD4+ T cells, and M1-like macrophages, were elevated in subtype 1. In contrast, tumor-associated immune cells, such as M2-like macrophages, were more elevated in subtype 2 (Figures 2C,D). Given these differences in the TME for these two subtypes, survival analysis showed that the overall survival of subtype 1 in the four cohorts was better than that of subtype 2 (Figure 2E).

To further investigate the differences between the two TME subtypes, we considered analyzing the biological function variation in the conception of signaling pathways. GSVA showed that all immune-related pathways, such as the IL-2/STAT5, IL-6/STAT3, and interferon response pathways, were enriched in subtype 1, while the TGF-β-, NOTCH-, PI3K/AKT-, and EMT-related pathways were enriched in subtype 2 (Figures 3A,B). ssGSEA with curated signaling pathway signatures showed that the CD8 T effector- and immune checkpoint-related pathways were activated in subtype 1. In contrast, tumor progression-related pathways such as WNT and EMT were activated in subtype 2 (Figures 3C,D). Based on ATAC-seq data from TCGA, differentially expressed peaks were identified between subtypes 1 and 2 (Figure 3E). GO analysis was processed on these differentially expressed peaks annotated by ChIPseeker, and the results showed that genes correlated with T cell activation had higher chromatin activities in subtype 1. In comparison, genes correlated with the regulation of GTPase and cell morphogenesis possessed higher chromatin activities in subtype 2 (Figure 3F).

Moreover, traditional GSEA was also conducted between subtypes in the four cohorts, which was consistent with the abovementioned results (Figure 4A). Using the ESTIMATE method, scores of stromal and immune cells were also higher in subtype 1 (Figure 4B). The expression levels of MHC molecules and immune checkpoint inhibitors (ICIs) are correlated with the activation of the antitumor immune response and the efficacy of immunotherapy. Most MHC molecules and ICIs were significantly different between the two subtypes and were especially higher in subtype 1 (Figure 4C). The abovementioned analysis showed that TME subtype 1 was highly correlated with immune-related phenotypes, while TME subtype 2 was positively associated with tumor progression and metastasis phenotypes. Thus, more comprehensive analyses containing FPI and mRNAsi were employed to analyze the ferroptosis level and the stemness index of single tumor tissue. As we have noticed that the initiation or metastasis of a malignant tumor might be highly correlated with cancer stem cells, we aimed to use mRNAsi to evaluate the differences between two TME subtypes. Moreover, ferroptosis was also a novel and vital phenotype that aroused our interest in further investigating the relationship with immune subtypes; FPI was employed here to satisfy our intention. Eventually, we found that the ferroptosis index (FPI) was higher in subtype 1 than in subtype 2, while the mRNA-based stemness index (mRNAsi) was higher in subtype 2 than in subtype 1 (Figures 4D,E). However, no significant difference was found in tumor mutation burden (TMB) (Figure 4F).

To further investigate the underlying mechanisms between the two TME subtypes, differentially expressed gene (DEG) analysis was conducted in four TNBC cohorts. Taking the intersection of DEGs in four cohorts (Figures 5A,B), 236 TME-related genes (TRG) were identified between TME subtypes, and all of them were upregulated in subtype 1 (Supplementary Table S6; Figure 5C). GO analysis showed that DEGs were highly enriched in T-cell activation and cell adhesion pathways (Figure 5D). For further analysis, a continuous variable called the TRG score by PCA was generated to quantify the different levels of TME in individual patients. The TRG score could well reflect the differences in TME subtypes in TNBC cohorts, and the TRG score was lower in subtype 1 (Figure 5E). Patients with low TRG score demonstrated a greater survival benefit than patients with high TRG score (Figures 5F,G). ssGSEA calculated with immune cell signature 1 showed that the infiltration levels of most immune cells were highly negatively associated with the TRG score (Figure 5H), and ssGSEA calculated with immune cell signature 2 also verified that most of the immune cells were higher in the low TRG score groups (Figure 5I). CIBERSORT analysis showed that as the TRG score was reduced, the percentage of cytotoxic T cells increased (Figure 5J). GSVA showed that immune-related pathways, such as the IL-2/STAT5, IL-6/STAT3, and interferon response pathways, were negatively correlated with the TRG score.

In contrast, glycolysis, the NOTCH signaling pathway, and protein secretion were positively correlated with the TRG score (Figure 6A). ssGSEA with curated pathway signatures verified that the TRG score was negatively linked with antigen processing machinery, CD8 T effector, and immune checkpoint and positively associated with WNT target pathways (Figure 6B). We show the genes involved in the above-curated pathway signatures with statistical significance in Figure 6C. Most of the genes involved in immune-related pathways were highly negatively correlated with the TRG score. The stromal and immune scores calculated by ESTIMATE were undoubtedly negatively correlated with TRG score in all TNBC cohorts (Figure 6D). The low TRG score group still had a higher FPI than the high TRG score group, but the mRNAsi and TMB showed no significant differences (Figures 6E–G). Due to differences in FPI between TME subtypes and TRG score groups, correlation analysis was conducted between TRG score and the expression of ferroptosis-related genes. We found that the expression of TNFAIP3, SOCS1, IFNG, ATM, ALOX5, PML, ISCU, and GCH1 was significantly negatively correlated with TRG score in four TNBC cohorts (Figure 6H). The TRG score showed no significant differences between the AJCC_T, AJCC_N, and stage groups, meaning that the TRG score was a novel factor regardless of clinical traits (Figure 6I).

To explore the association between the TRG score and drug response, we evaluated the estimated IC50 value of 138 drugs included in the GDSC database in four TNBC cohorts. Correlation analyses were conducted between the TRG score and predicted IC50 values (Supplementary Figure S2A). Drugs with significant differences in more than three cohorts were regarded as potential therapeutic drugs; we found that eight drugs were sensitive to the high TRG score group, and 49 drugs were sensitive to the low TRG score group (Figure 7A). The TRG score might logically be related to the efficacy of immunotherapy due to its apparent association with immune cell infiltration and activation. TIDE was utilized to predict the immunotherapy response of TNBC patients, and the TRG score was lower in the immunotherapy response group (Figure 7B). Moreover, TIDE analysis showed that the TRG score was apparently negatively correlated with markers of immunotherapy response and positively correlated with CAFs, myeloid-derived suppressor cells (MDSCs), and TAM M2 (Figure 7C). Lacking TNBC datasets that received immunotherapy, we selected three cohorts that received anti-PDL1, anti-PD1, and anti-MAGE-A3 therapy in bladder cancer (BLCA) and skin melanoma (SKCM) to verify the immunotherapy response prediction value of the TRG score. First, TRG scores were calculated across cancers in TCGA. TRG scores were prognostic risk factors (Supplementary Figure S2B) and were negatively correlated with immune cell infiltration levels in most cancer types, especially in BLCA and SKCM (Figure 7D). Then, we calculated the TRG score in the three immunotherapy cohorts. Interestingly, we found that the TRG score was also a risk factor in IMvigor210 (Supplementary Figure S2C), and patients with a high TRG score and low TMB presented the worst survival advantage (Figure 7E). Correlation analysis further validated that the TRG score was negatively correlated with the expression of MHC, costimulatory, adhesion molecules (Supplementary Figure S2D), and immune cell infiltration levels (Figure 7F). Moreover, a higher TRG score was associated with disease progression (PD), indicating that a higher TRG score might indicate poor response after immunotherapy (Figures 7G,H). As most solid tumors exhibited one of three distinct immunological phenotypes, immune inflamed, immune excluded, or immune desert, studies in the IMvigor210 cohort classified each sample into one of these immune phenotypes (Mariathasan et al., 2018). The immune inflamed phenotype was thought to be rich in immune cell infiltration and sensitive to immunotherapy, while the immune desert was on the contrary. Immune excluded phenotype was surrounded by many immune cells, but the immune cells were confined to the periphery of the tumor cell matrix. We found that a higher TRG score was associated with desert-resistant phenotypes, while inflamed phenotypes possessed a lower TRG score than desert and excluded phenotypes (Figure 7I). In the anti-MAGE-A3 cohort, the TRG score was also negatively correlated with immune cell infiltration levels (Figure 7J) and was lower in the response group (Figure 7K). Similar results could be seen in the anti-PD1 cohort (Figure 7L), although the differences among response groups showed no significance (Supplementary Figure S2E).

The TRG score was calculated as described earlier in 80 triple-negative breast cancer samples from the West China Hospital (TNBC_WC) cohort. It was found that a higher TRG score was related to the disease progression rate (Figure 8A) and poor survival probability (Supplementary Figure S2F) of TNBC patients. ssGSEA also showed a strong correlation between the TRG score and immune cell infiltration levels in TNBC_WC (Figure 8B), as well as the results of the ESTIMATE score (Figure 8C). Some MHC molecules and ICI targets were also negatively correlated with the TRG score, especially PDL1 and PDCD1LG2 (Figure 8D). For pathway analysis, pathways that were associated with the TRG score in TNBC_WC were almost the same as the results in the four training cohorts (Figure 8E). Correlation analysis showed that the TRG score was positively correlated with the FPI and mRNAsi (Figure 8F). The response group predicted by TIDE analysis showed a lower TRG score than the no response group (Figure 8G). Additionally, we predicted the drug IC50 by GDSC analysis and performed correlation analysis with the TRG score (Figure 8H) and intersection drugs with the results in the training cohorts, as shown in Figure 8I. These results illustrated that the TRG score was a novel and robust method to measure immune cell infiltration levels and therapy efficacy.

Considering the accessibility of the TRG score, we aimed to shrink the members of the TRG score and simplify the formula modes to predict the prognosis of TNBC patients. First, survival analysis was processed for 236 TME-related DEGs in TCGA cohorts; 84 genes with a p value < 0.05 were selected for further research (Supplementary Figure S3A). Here, iteration LASSO was then used to simplify the members of the TRG score; after multiple attempts to reach the highest 5-year AUC, we finally constructed a prognostic signature with 20 members from TRG score members (Supplementary Table S7; Figure 9A). We could see that there were 20 genes with the most frequencies of occurrence in 1,000 operation iterations in LASSO algorithms, and prognostic signatures with these 20 genes could reach a high area under the curve (AUC) of ROC for 5 years of survival in the TCGA cohort. To provide a convenient approach for predicting the survival probability of a patient with TNBC, we constructed predictive nomograms with the 20 genes generated previously. We developed a nomogram based on the Cox regression model to predict the 5- and 8-year survival probability for TNBC patients (Supplementary Figure S3B). The calibration plots for the 5- and 8-year survival showed an optimal agreement between the nomogram-predicted and observed OS, which was used to evaluate the accuracy of the prediction signature (Figure 9B). For validation of the prognostic value of the 20-gene signature, the patients in the high-risk group showed a worse prognosis than those in the low-risk group, and the same condition could be seen in other TNBC cohorts (Supplementary Figure S3C). Before the prognostic signature was constructed, we conducted a correlation analysis between the TRG score and the expression of 236 TME-related DEGs in four TNBC cohorts. In view of most genes highly correlated with the TRG score, we are supposed to simplify the TRG score by these 20 genes. Surprisingly, the simplified TRG score (sTRG score) calculated based on the expression of these 20 genes was highly positively correlated with the TRG score in the TNBC_WC cohort (Figure 9C), and patients in the high sTRG score group also showed worse DFS than those in the low sTRG score group (Figure 9D). Not unexpectedly, the correlation coefficients between the TRG score and sTRG score were almost close to 1 in other cohorts, which means that they were virtually interchangeable (Figure 9E). However, the risk score showed no significance with the TRG score, indicating that the risk score was a novel factor generated by the iterative LASSO regression model. Eventually, we set up a coexpression network for 20 genes, and we found strong correlations among them (Figure 9F). The visualization of attribute changes in individual patients using an alluvial diagram indicated that the TRG score might be a powerful method to direct therapeutic efficacy or prognostic risk for TNBC patients (Figure 9G).