This research was conducted on publicly available database through a series of bioinformatics methods. Transcript profiles and corresponding clinical information of six cohorts containing a total of 694 TNBC patients were acquired to construct and validate the risk signature after removing the samples with unknown survival time and outcome. Details were as follows: microarray dataset GSE103091 (107 samples) was downloaded from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) and was selected as the training cohort on account of the optimal sample size. Another three microarray datasets named GSE16446 (107 samples), GSE20685 (225 samples) and GSE20711 (78 samples) were also obtained from GEO database and used as the validation cohorts, in addition, TNBC RNA sequencing datasets which downloaded separately from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov) and Molecular Taxonomy of Breast Cancer International Consortium (METABRIC, https://www.cbioportal.org/) were used as another two validation cohorts, respectively. Summary information of above cohorts was listed in Table 1. Besides, three real-world immunotherapy cohorts (GSE91061, GSE135222, IMvigor210) were chosen to verify the ability of the risk signature on prediction of immunotherapy response. All the raw data were normalized and log₂ transformed.

The Tumor Immunophenotype (TIP) database (http://biocc.hrbmu.edu.cn/TIP) is a webtool that can assess the immune microenvironment on the base of the cancer-immunity cycle (Xu et al., 2018). The marker genes were retrieved from the TIP database and employed to classify TNBC patients into diverse clusters in the training cohort through “ConsensusClusterPlus” R package (Wilkerson and Hayes, 2010). Pam algorithm and “spearman” were used as the metric distance. Each bootstrap process including 80% of the training cohort of patients and was repeated by 500 times. The number of clusters was set to be 2 to 10, and the optimal classification was determined by calculating the consistency matrix and consistency cumulative distribution function (CDF).

The stromal score, immune score and estimate score of training cohort were calculated by ESTIMATE algorithm and were used to compare the immune infiltration between different subtypes and different risk models (Danilova et al., 2019). Then, the c2. cp.kegg.v7.5.1 gene set was downloaded from Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/) and employed to quantify the pathways through “ssGSEA” method. The infiltration level of 10 immune cells were also quantified through “MCPcounter” algorithm. Next, we also calculated the relative infiltration level of 22 kinds of immune cells by CIBERSORT method. Then, the characteristic genes of 28 immune cells which obtained from previous study (Charoentong et al., 2017) were used to calculate the degree of infiltrating immune cells between different risk groups in TNBC.

Differential expression analysis between diverse molecular subtypes and risk models were conducted by “limma” package and visualized through volcano plot. The selection criterion was |log₂FC| > 2 and FDR < 0.05 for molecular subtypes, and |log₂FC| > 1.5 and FDR < 0.05 for diverse risk models, respectively. The “WebGestaltR” package was used to further investigate the functions involved in differential expressed genes (DEGs), and the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed. The Gene set enrichment analysis (GSEA) was then conducted to further analysis the difference of biological functions between different groups based on the Hallmark gene set through “clusterProfiler” package.

To develop a consensus model with high accuracy and stable performance, we integrated 10 machine learning algorithms and 101 algorithm combinations. The integrative algorithms included random survival forest (RSF), elastic network (Enet), Lasso, Ridge, stepwise Cox, CoxBoost, partial least squares regression for Cox (plsRcox), supervised principal components (SuperPC), generalised boosted regression modelling (GBM), and survival support vector machine (survival-SVM). The signature generation procedure was as follows: (a) Firstly, univariate Cox regression was employed to identify the prognostic genes in training cohort; (b) Then, 101 algorithm combinations were performed on the prognostic genes to fit prediction models based on the leave-one-out cross-validation (LOOCV) framework in the training cohort; (c) All models were detected in five validation datasets (GSE20685, METABRIC, TCGA-TNBC, GSE16446, and GSE20711); (d) For each model, the Harrell’s concordance index (C-index) was calculated across all validation datasets, and the model with the highest average C-index was considered optimal. The risk score was calculated as following formula: R i s k   s c o r e = ∑ 1 n Exp * c o e f f i c i e n t .

The “maftool” package was used to explore the somatic mutations in TCGA-TNBC patients and the top 10 mutated genes were presented in waterfall plot. Besides, the copy number variation (CNV) data of TCGA-TNBC was also downloaded and used to display the proportion of deletion and amplification of genes according to the risk model.

In order to explore the superiority of the risk signature, the time-dependent area under curves (tAUC) of the signature and other clinicopathological features were analyzed and compared in METABRIC cohort and TCGA-TNBC cohort, respectively. Then, the TIDE score of three immunotherapy cohorts (GSE91061, GSE135222 and IMvigor210) were calculated through the online tool (http://tide.dfci.harvard.edu/) for immune treatment effect evaluation. The tAUC of risk signature and TIDE score were analyzed in the three cohorts and the comparison between these two indicators were also performed to distinguish the sensitivity and specificity to immunotherapy response.

To further investigate potential therapeutic target drugs in the high-/low-risk group, we used the drug-sensitive cell lines in the CCLE database (https://portals.broadinstitute.org/ccle) as a training set. The drug sensitivity of each patient in the GSE103091 cohort was predicted by CTRP and PRISM methods. Then screening for potential regulation of drugs through the setting as |cor| > 0.3.

R software (https://www.r-project.org, version 4.1.3) and GraphPad Prism 8.0 (GraphPad Software Inc., San Diego, CA, United States) were used for all statistical analysis and visualization. Univariate Cox regression analysis was performed to evaluate the significant prognostic genes. Quantitative data were compared between different groups through Wilcoxon rank sum test. Relationships between risk scores and expression levels of different genes were examined by Spearman’s correlation analysis. Unless otherwise specified, p < 0.05 was considered as statistically significant.

The workflow of our research was presented in Figure 1. The research contents mainly included three parts: 1) Identifying different prognostic immune types and their related DEGs and pathways; 2) Building risk models based on machine learning methods according to these DEGs and exploring the regulatory pathways of different risk models, as well as the relationship with immune cells and chemokines; 3) Analysis of potential targeting drugs for different risk groups.

A total of 166 marker genes were collected in the seven stages of the cancer-immunity cycle, including checkpoints, cytotoxic factors, chemokines, and major histocompatibility complex (MHC) molecules (Figure 2A). The CDF delta area curve indicated that k = 2 could gather relatively stable clustering results which named Cluster 1 (C1) and Cluster 2 (C2) (Figures 2B–D). Further analysis of the prognostic characteristics of these two molecular subtypes revealed significant overall survival (OS) difference between them in the training cohort (Figure 2E). In general, C1 subtype had a poor prognosis compared with C2 subtype (p < 0.05). Similar results were observed in GSE20685 cohort (Figure 2F, p < 0.05). Then, PCA analysis was conducted based on Neutrophils marker genes, and PCA dimension reduction distributions of the two subtypes were shown in Figure 2G. The results demonstrated an obvious batch effect between the two cluster samples.

Above results indicated that the patients in C2 showed a better prognosis than the patients in C1. Next, the study explored the differences in immunity between these two clusters. Obviously, C2 displayed a higher immune score, stomal score, and ESTIMATE score compared with the C1 (Figure 3A). In addition, the two clusters showed significant differences in most immune-related pathways, including JAK-STAT signaling pathway, NF-kappa B signaling pathway, Toll-like receptor pathway, B cell receptor signaling pathway, T cell receptor signaling pathway and inflammatory response (Figure 3B). Besides, the quantitative infiltration levels of most immune cells were much higher in C2 than C1, suggesting the patients in C2 may act more immune activity, detailed information were presented in Figures 3C, D. The difference of KEGG pathways were visualized by heatmap and a coincident result was obtained, that is the C2 exhibit higher activity in tumor immunity related pathways, such as apoptosis and JAK-STAT pathway, et al. (Figure 4A). GSEA analysis further suggested the C1 showed positively correlation with cancer-related pathways, including G2M checkpoint and E2F targets, et al., while the C2 presented positively relationship with immune-related pathways, including INF-alpha response and INF-gamma response, et al. (Figure 4B). Finally, the results of “ssGSEA” score showed that five tumor-related pathways were significantly different between two clusters, including WNT, TP53, PI3K, NRF1, and HIPPO, which have been linked to the development and progression of cancer and have great potential in predicting the prognosis of TNBC patients (Figure 4C).

In order to further investigate and verify the differences of biological functions between the two clusters, the differential analysis was performed to complete this task. The heatmap was showed in Figure 5A. A total of 590 DEGs were collected, among which contains 11 upregulated genes and 579 downregulated genes (Figure 5B). KEGG enrichment indicated the DEGs were mainly participate cytokine-cytokine receptor pathway, chemokine signaling pathway, cell adhesion molecules, hematopoietic cell lineage and viral protein interaction with cytokine and cytokine receptor (Figure 5C). GO biological process results showed these DEGs were enriched in T cell activation, leukocyte mediated immunity, leukocyte cell-cell adhesion, regulation of T cell activation and lymphocyte mediated immunity (Figure 5D). The top five cellular component (CC) and molecular functions (MF) were showed in Figures 5E, F, respectively.

A total of 454 genes were collected from the intersection of DEGs, TCGA, METABRIC, GSE20685, GSE20711 and GSE16446 data sets (Figure 6A). Then, univariate Cox analysis was performed to calculate the relationship between 454 DEGs and TNBC prognosis in the training cohort, and 30 prognostic genes were finally screened after filtering p < 0.05, among which were all protective factors. Next, these 30 genes were used to develop a consistent prognostic model through our integrated program based on machine learning approach. In brief, 101 prediction models were filtered through the LOOCV framework, meanwhile, the C-index for each model was also calculated in the training cohort and validation cohorts to select the most outstanding candidate. Interestingly, the optimal model was a combination of CoxBoost and RSF, with the highest average C-index equal to 0.622 (Figure 6B). Finally, seven key genes were screened to establish the prognostic signature (Figure 6C). The risk score was calculated as above mentioned and subsequently normalized by the “scale” method. The Z-score equal to zero was selected as the cut-off value to separate the cohorts into high-risk group and low-risk group. Survival analysis indicated that the patients in low-risk group presented significant longer OS compared with the high-risk patients both in the training cohort and the validation cohorts (all p < 0.05, Figure 6D). So, this is considered a robust model and worthy to further study.

A total of 1,145 DEGs were collected between high-risk group and low-risk group in the training cohort (Figure 7A). Functional enrichment analysis indicated these DEGs may play a vital role in immune-related pathways and biological functions (Figure 7B). GSEA results showed that immune-related pathways were significantly enriched in the low-risk group, including innate immune system, adaptive immune system, signaling by GPCR, hemostasis, and cytokine signaling in immune system (Figure 7C). On the other hand, only three pathways were enriched in the high-risk group and most of them were related to cell proliferation process (Figure 7D). Due to the strong correlation between risk characteristics and immune-related biological pathways, we further investigated the association between risk scores and tumor-infiltrating immune cells. Firstly, we use an estimation algorithm to quantify the overall somatic immune cells based on the TCGA sequence. From Figure 8A, we found that the risk score and the immune score presented a strong negative correlation (p < 0.001), indicating that the low-risk group which evaluated based on our model had a higher immune infiltration. We further analyzed the differences in the distribution of somatic immune cells between the two subpopulations and found significant differences in T cells and three types of macrophages in the low-risk group (Figure 8B). Then, using the characteristic genes of 28 immune cells obtained from previous study (Charoentong et al., 2017), the infiltration scores of 28 immune cells were calculated by “ssGSEA” method, and 9 out of 12 T cells showed significant differences in the two risk groups (Figure 8C). Furthermore, we analyzed the correlation between risk score and these 12 types of T cells (Figure 8D), and the results showed that there was a strong negative correlation between risk score and T cell scores. It was also found that M1 macrophage score was negatively correlated with risk score, while M0 and M2 showed an opposite trend (Figure 8E). The scores of three macrophage-related pathways were also significantly negatively correlated with risk scores (Figure 8F). It can be seen that 14 out of 40 chemokines were significantly different between the two risk groups, suggesting that different risk groups may have different degrees of immune cells infiltration, and these differences may directly affect the progress of tumor and the effect of immunotherapy (Figure 9A). In addition, we calculated and compared the expression of chemokine receptor genes in the different risk groups and found that there were significant differences in the expression of chemokine receptor genes, including CCR1, CCR2, CCR5, CCR6, CCR7, CCR8, CXCR2, CXCR5 and CXCR6 (Figure 9B). Finally, further analysis indicated the risk score was significant negatively correlated with these genes (Figure 9C). Thus, our study identified and validated two robust immune subtypes based on comprehensively bioinformatics methods.

To further explore the mutational landscape between diverse risk groups, the somatic mutational data was downloaded from TCGA database and used to complete the procession. As a result, top 10 mutated genes were shown in waterfall plot, including TP53, TTN, MUC16, SYNE1, FAT3, SPTA1, CSMD3, DMD, DYNC2H1 and PIK3CA (Supplementary Figure S1A). CNV analysis presented the proportion of deletion and amplification of these seven genes were remarkable changed, especially the CD160 (Supplementary Figure S1B). These findings suggested that these genes with significant mutational differences may play an important role in different immune scores.

In order to verify the prognostic performance of risk signature, we conducted tAUC analysis to compare the specificity and sensitivity with other clinicopathological features. Results showed the risk score played a significantly strong survival prediction ability in METABRIC cohort (Figure 10A). Similar results were viewed in the TCGA cohort (Figure 10B). In addition, we calculated the AUC values of the risk model and TIDE score in IMvigor210 cohort, GSE135222 and GSE91061, respectively. Besides, the prognostic value of risk signature and TIDE score were also compared in these three immunotherapy cohorts. All the results indicated the risk signature displayed better ability in prognosis prediction and immunotherapy response (Figures 10C–K).

To assess the usefulness of risk models in clinical treatment, we analyzed chemotherapy drug sensitivity in low- and high-risk patients. We used the CCLE database of drug-sensitive cell lines as the training set and the GSE103091 data set as the validation set. In the end, a total of 18 CTRP (Figure 11B) and 26 PRISM (Figure 11D) compounds were obtained. The results showed that the high-risk group had higher IC50, indicating that the high-risk group was not sensitive to chemotherapy. Results showed that high-risk groups were more sensitive to MR-1220, GSK2110183 and temsirolimus. Therefore, high-risk samples should be sensitive to these compounds, which may be new options for future TNBC treatment.