ICGs expression patterns in TNBC
We first investigated the ICGs expression in TCGA-TNBC cohort. As shown in Fig. 1a, 53 of 53 ICGs were aberrantly expressed in TCGA-TNBC datasets. Among them, 44 ICGs were high expression in TNBC tissues compared to normal tissues. A complex and finely tuned signaling ICGs network was presented by STRING database (Fig. 1b). The prognostic value of ICGs in TNBC was evaluated by Univariate Cox regression analysis and Kaplan–Meier (KM) log-rank test (Table S1). We next combined TCGA-TNBC and METABRIC-TNBC dataset to perform consensus clustering analysis. Unsupervised clustering analysis uncovered two distinct ICGs clusters in TNBC (Fig. 1c and Fig. 1d). Interestingly, TNBC patients with ICGscluster-B had a longer survival time than ICGscluster-A (Fig. 1e). ICGs expression profiles analyses found that ICGscluster-B subgroup exhibited higher expression of checkpoint gene compared to ICGscluster-A (Fig. 1f).
Construction of ICGs-related gene clusters in TNBC
Principal component analysis was conducted to classify TNBC samples into two classes with clear differences after the dimensionality reduction (Fig. 2a). To further explore the biological behaviors between different ICGs cluster, the results of GSVA revealed that ICGscluster-B was highly enriched in various immune-related pathways, such as natural killer (NK) cell mediated cytotoxicity, T cell receptor signaling pathway, interferon gamma response and IL-6 JAK STAT3 signaling (Fig. 2b and Fig. 2c). To further characterize the immune infiltration in different ICGs clusters, we applied the ssGSEA method in TNBC samples. As shown in Fig. 2d, ICGcluster B patients had higher abundance of antitumor immune cells including activated CD8 T cell, activated dendritic cell and so on. We next identified ICGs-related DEGs using “limma” algorithm and 570 DEGs were recognized between two ICGs clusters (Table S2). GO and KEGG enrichment analysis displayed that these ICGs-related DEGs mainly involved in T cell activation, MHC protein complex, immune receptor activity and cytokine−cytokine receptor interaction (Supplement Fig. 1 and 2). We next performed univariate cox regression analysis to screen the prognostic ICGs-related DEGs (Table S3). Consensus clustering analysis further showed two ICGs-related geneClusters based on the expression of prognostic ICGs-related DEGs (Fig. 2e and Fig. 2f). Survival analysis indicated that ICGs geneCluster B was associated with better prognosis (Fig. 2g). The heat plot showed ICGs geneCluster-B patients were largely overlapped with ICGscluster-B subtype (Fig. 2h). As expected, higher expression of ICGs was observed in ICGscluster-B subtype (Fig. 2i).
Generation of an ICGs-related risk model
We next constructed Lasso logistic regression model to assess the immune status for individual TNBC patients. Four ICGs-related genes were selected to generate the best risk score model (Fig. 3a and Fig. 3b). The riskscore was calculated by the formula ICGs-related riskscore = expression level of FGD3 *(-0.22) + expression level of TMEM176A*(-0.19) + expression level of CD1B*(-0.57)+ expression level of MATK*(-0.24). We also noted that TNBC patients with ICGs geneCluster-B subtype or ICGscluster-B subtype had a lower ICGs-related riskscore, indicating ICGs-related risk score was a poor prognostic indicator (Fig. 3c and Fig. 3d). The relationship among ICGs clusters, ICGs-related geneCluster and ICGs-related riskscore was illustrated as Fig. 3e. ICGscluster-B subtype had a high overlap with the ICGs geneCluster-B subtype that showed a low ICGs-related riskscore with a favorable survival prognosis. Finally, we found TNBC patients with low ICGs-related riskscore had a higher expression of ICGs than those with high ICGs-related riskscore (Fig. 3f).
Validation of the efficacy of ICGs-related riskscore model
We then divided TNBC cohort into training cohort (N = 211) and validation cohort (N = 210). Patients were classified into low-risk and high-risk groups based on their median risk scores. In both the training and validation cohort, TNBC patients with high ICGs-related riskscore predicted shorter OS than those with low ICGs-related riskscore (Fig. 4a and Fig. 4b). Time-dependent ROC curves also exhibited excellent stability of our ICGs-related riskscore model in training cohort (5‐year AUC, 0.705; 3‐year AUC, 0.684; 1‐year AUC, 0.746; Fig. 4c ) and validation cohort (5‐year AUC, 0.606; 3‐year AUC, 0.644; 1‐year AUC, 0.605; Fig. 4d). The distribution of ICGs-related riskscore in training cohort and validation cohort was presented in Fig. 4e and Fig. 4f. Moreover, the scatter plot showed that the mortality rate of them increased with the increase of ICGs-related riskscore in both the training cohort and validation cohort (Fig. 4g and Fig. 4h). As shown in Fig. 4i and Fig. 4j, high expression of FGD3, TMEM176A and MATK was observed in low ICGs-related riskscore group.
Construction of Nomogram for predicting OS
PCA and t-SNE analysis were applied to separate between two different risk subgroups. As shown in Fig. 5a-5d, different risk subgroups were divided into two obvious directions. To explore whether the effect of our ICGs-related riskscore was an independent factor, we combined clinicopathological variables including age, stage tumor size and node status with ICGs-related riskscore to perform Cox analyses. Univariable Cox regression analysis showed that ICGs-related riskscore was a high-risk factor for TNBC patients [HR =1.714, 95% CI (1.323, 2.221), p < 0.001, Fig. 5e]. Multivariate Cox analysis further supported that ICGs-related riskscore could serve as the adverse independent prognostic factor [HR=1.544, 95% CI (1.173, 2.003), p =0.002, Fig. 5f]. We next developed a nomogram combining the clinicopathological features and ICGs-related riskscore to help clinicians to predict the survival of TNBC. Calibration curves for OS prediction at 1, 3, and 5 years were as shown in Fig. 5g. In addition, if the total nomogram score of one patient was 138 points, the estimated survival probabilities for this patient were 96.8%, 77.5%, and 66.2% for 1, 3, and 5 years, respectively (Fig. 5h).
Immune Characteristics and genomic features in different risk subgroups
We further investigated the immune cells infiltration in different risk subgroups by CIBERSORT algorithm. As in Fig. 6a, the expression of CD1B was strongly associated with the infiltration of CD8 T cells, T help cells and resting dendritic cells. FBP1 had a strong correlation with infiltration of CD8 T cells, memory CD4 T cells and M1 macrophage. Furthermore, a strong positively correlation was observed between the infiltration of CD8 T cells, memory CD4 T cells and activated NK cells. Overall, ICGs-related riskscore was negatively correlated with infiltration of activated NK cells and CD8 T cells (Fig. 6b and Fig. 6c). M1 macrophage and M2 macrophage cells infiltration were positively correlated with ICGs-related riskscore. However, M0 macrophage cells infiltration was negatively correlated with ICGs-related riskscore (Supplement Fig.3). The results of ESTIMATE algorithm also demonstrated that TNBC patients with low ICGs-related riskscore had a higher level of ImmuneScore and ESTIMATEScore (Fig.6d). In addition, the immune function in low ICGs-related riskscore patients might be attributed to the activation of various immune cells (Fig.6e). High ICGs-related riskscore was proved to be positively associated with the stemness index DNAss and RNAss (Supplement Fig.4). Then, we investigated the genomic features in different ICGs-related riskscore subgroups of TCGA-TNBC cohort. The waterfall plot demonstrated that high ICGs-related riskscore group presented an increased tumor mutation burden compared with the low ICGs-related riskscore group (Fig.6f and Fig.6g). TTN (33% vs.12%) had higher mutation rates in high ICGs-related riskscore group.
Predictive value of ICGs-related riskscore model in ICIs therapy
Associations of clinical parameters and ICGs-related riskscore model were also explored. As the Fig.7a shown, our ICGs-related riskscore was significantly correlated with tumor stage and tumor size. We next adopted IPS score to assess the potential clinical efficacy of ICIs in different riskscore subgroups. IPS score is a valid indicator as a measure of immunogenicity. There was no significant difference for TNBC patients with PD-1 negative expression (Fig.7b and Fig.7d). More remarkably, low ICGs-related riskscore might serve as a valid candidate biomarker for ICIs in patients with PD-1 positive expression (Fig.7c and Fig.7e). Although TNBC patients with high ICGs-related riskscore might not be the most suitable candidates for ICIs therapy, drug sensitivity analysis revealed that TNBC patients with high ICGs-related riskscore were more sensitivity to gemcitabine and methotrexate (Fig.7f and Fig.7g).