Tumor microenvironment characterization in triple-negative breast cancer identifies prognostic gene signature

Different TME clusters of TNBC

Immune cell populations modulate different immune responses through the infiltration of TNBC TME and lead to antitumor effects. To investigate whether there was immune-related TME heterogeneity in TNBC, we collected 29 immune-related terms from published reports. Then, single-sample Gene Set Enrichment Analysis (ssGSEA) method was applied to assess the enrichment status of these 29 immune-related terms for each TNBC patient from the GSE58812 cohort (training set consisting of 107 TNBC samples). Unsupervised clustering of the enrichment scores of the 29 immune-related terms identified two clusters in the TNBC patients (Figure 1A): Cluster 1, which included 65 cases of TNBC, was defined as a “low immune infiltration” cluster because of its low enrichment in the immune-related terms; whereas Cluster 2, which included 42 cases of TNBC, was defined as a “high immune infiltration” cluster because of its high enrichment in the immune-related terms. This result reflected different patterns of infiltration of the immune cells based on the adaptive and innate immune systems of each patient. To validate the results of the TME clusters, we used two independent algorithms, TIMER2.0 and CIBERSORTx, to estimate the TME cell abundance. We observed similar results as most of the tumor infiltrating immune cells were significantly increased in Cluster 2 than in Cluster 1 (TIMER2.0: Supplementary Figure 1; CIBERSORTx: Supplementary Figure 2). These results indicated that Cluster 2 was associated with high immune infiltration. Further analysis showed that the patients in the two clusters manifested different outcomes. For instance, the patients in Cluster 2 displayed better metastasis-free survival (MFS; Figure 1C; P = 0.012) and overall survival (OS; Figure 1D; P = 0.0054) compared with the patients in Cluster 1. To validate the identified TME clusters, we independently performed ssGSEA and unsupervised clustering for the validation METABRIC cohort (validation set consisting of 299 TNBC patients). The results showed that all the patients were also categorized into two heterogeneous clusters with 131 and 168 patients each (Figure 1B). Significant differences in OS and pattern of cellular infiltration were observed between these two TME clusters as well (P = 0.0034; Figure 1E). Further, in Cluster 1, microenvironmental cell permeability was relatively low. In contrast, Cluster 2 exhibited both innate and adaptive immune responses with high abundance of plasmacytoid dendritic cells (pDCs), immature dendritic cells (iDCs), macrophages, B cells, and a high penetration of T cells, cytotoxic cells. We systematically compared the abundance of immune-related terms in the two clusters. The results likewise confirmed that the abundance of immune-related terms in Cluster 2 was significantly greater than in Cluster 1 (Figure 1F).

Figure 1. TME Clusters of Triple-negative breast cancer. (A) TNBC TME clusters in training set. Unsupervised clustering of tumor microenvironment immune cells for 107 TNBC patients from training set. Clinicopathological information including age, MFS, OS, as well as TME clusters, is shown in annotations above. (B) TNBC TME clustering in the validation set. Unsupervised clustering of tumor microenvironmental immune cells from 299 cases of TNBC from validation set. The OS and the number of months of survival are shown in the comments above. (C) Kaplan–Meier curves for Metastasis-free survival (MFS) stratified by TME clusters in the GEO cohort. (D) Kaplan–Meier curves for overall survival (OS) stratified by TME clusters in the training set. Hierarchical clustering was performed using Euclidean distance and ward linkage. (E) Kaplan–Meier curves for overall survival (OS) stratified by TME clusters in the validation set. (F) The abundance of immune-related terms in the Cluster 1 and 2.

Clinical features of different TME clusters of TNBC

We compared the clinical characteristics, including age, Nottingham prognostic index, cellularity, chemotherapy, laterality, tumor grade, tumor size, tumor stage, and positive lymph nodes between the two distinct TME clusters (Table 1). Cluster 1 (131 patients) was predominantly associated with low-immune infiltration and Cluster 2 (168 patients) was characterized by high immune infiltration. Cluster 2 had significantly smaller tumor sizes, while Cluster 1 had larger tumor sizes (P<0.001). However, the differences were not significant with respect to age, tumor grade, and tumor stage (age, P = 0.359; tumor grade, P = 0.725; tumor stage, P = 0.773), indicating that the TME clusters were independent of age, and tumor grade and stage.

Differential gene and pathway analysis of different TME clusters of TNBC

To identify the potential biological characteristics of the two TME clusters, differentially expressed genes (DEGs) were identified between TME Cluster 1 and Cluster 2. In TME Cluster 1, 37 genes were significantly upregulated and 778 genes were significantly downregulated when compared with the expression in Cluster 2 (Figure 2A; P < 0.05). DEG expression between the two clusters was visualized using a heatmap (Figure 2B). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the DEGs revealed enrichment of immune-related pathways, such as “leukocyte migration” and “leukocyte proliferation” in TME Cluster 2, which was supportive of the high immune cell infiltration patterns in this cluster (Figure 2C). In TME Cluster 1, significant upregulated genes were associated with development (Figure 2D), which were not immune-related. We further examined the expression of several immune checkpoint molecule genes between the two clusters, and found that four genes, PD-1, PD-L1, CTLA-4, and TIM-3, were significantly upregulated in TME Cluster 2 than in TME Cluster 1 (Figure 2E).

Figure 2. Differential gene and pathway Analysis of Different TME Cluster of TNBC. (A) The volcano plot showing the differentially expressed gene (DEGs) between TME Cluster 1 and Cluster 2. Blue is highly expressed in Cluster 2 and red is highly expressed in Cluster1. (B) Heatmap shows the expression of DEGs between the two clusters. (C) KEGG functional enrichment analyses of the up-regulated genes in Cluster 2. (D) KEGG functional enrichment analyses of the up-regulated genes in Cluster 1. (E) Expression of four immune checkpoint molecules genes in the two clusters.

Establishment and validation of TME signature

To establish a prognostic prediction model by TME-related characteristics, we applied univariate Cox regression to the survival data of the TNBC patients and expression profiles of 815 DEGs. Genes with P values < 0.1 in the univariate Cox regression were further subjected to multivariate Cox regression and stepwise regression [16, 17]. The final TME signature consisted of eight genes, PPFIA4, TNFRSF1B, ARHGAP9, ZNF831, CTLA4, BLK, ANKRD22, and CLEC4E. The hazard ratio (HR) and P values of the eight genes in the TME signature are shown in Figure 3A. Next, TME scores were obtained according to the TME signature, and the patients were categorized into high- and low-risk groups according to the median TME score. Survival analysis demonstrated that the TME scores exhibited strong power to distinguish good outcomes from poor in the TNBC patients in the training set (P < 0.001; Figure 3B). Patients at high-risk had significantly shorter survival compared to those at low-risk. To confirm our findings, we validated the TME signature in the validation set. Using the same model in the training set, we calculated the TME scores for each patient in the validation set, and evaluated the relationship between the TME score and their survival status. As shown in Figure 3C, patients at high-risk had significantly shorter survival compared to those at low-risk in the validation set. The distribution of the risk scores, survival statuses, and prognostic gene expression in the training and validation sets is shown in Figures 3D, 3E. To assess the prediction accuracy of the prognostic prediction model, receiver operating characteristic (ROC) curve was used. The ROC curve showed that the TME score demonstrated good discrimination in the training set with an area under the curve (AUC) of 0.84, and an acceptable discrimination in the validation set with an AUC of 0.64 (Figure 3F). Together, these results indicated that the prognostic prediction model based on the TME scores could be used to predict the survival of the TNBC patients. Since proteins seldom affect in isolation, it is important to know the interactions between the genes included in TME signature to identify the hub genes. We constructed a protein-protein interaction (PPI) network using STRING database. The PPI network comprised 8 nodes and 10 edges (Supplementary Figure 3). TNFRSF1B and BLK had the maximum neighboring genes, and were identified as the hub genes. We used immunohistochemistry (IHC) to analyze the correlation of TNFRSF1B and BLK expression with the prognosis of TNBC patients using an in-house paraffin embedded sample. Representative IHC figures are shown in Figure 4A. TNBC patients with high expression of either TNFRSF1B or BLK showed poorer prognosis than those with low expression of these hub genes (Figure 4B, 4C). These results were also consistent with the TME signature (HRs of TNFRSF1B and BLK were >1; Figure 3A), suggesting that TNFRSF1B and BLK were risk factors for poor prognosis in TNBC patients.

Figure 3. TME signature is a promising marker of survival in TNBC patients. (A) the HR and P-value from the eight genes in the prognostic model. (B) The survival curves for patients with high-risk and low-risk in the training set. (C) The survival curves for patients with high-risk and low-risk in the validation set. (D) Distribution of risk scores, survival profiles, and heat maps showing characteristic expressions of the low and high risky groups in the training set. Top panels indicate risk scores of patients. Middle panels depict survival status and survival time of patients distributed by risk score. Bottom panels display heatmap of expression for eight predictive factors distributed by risk score. (E) Distribution of risk scores, survival profiles, and heat maps showing characteristic expressions of low- and high-risk groups in the validation set. (F) Comparison of the predictive accuracy of the training set (GSE58812) and the validation set (METABRIC).

Figure 4. The correlation between TNFRSF1B and BLK expression levels and TNBC prognosis. (A) Representative immunohistochemical figure showed the high and low expression of TNFRSF1B and BLK. (B) The correlation between TNFRSF1B expression levels and TNBC prognosis. (C) The correlation between BLK expression levels and TNBC prognosis.

Correlation of TME signature with different TME clusters

Relationship between the TME signature and TME clusters were evaluated by a correlation analysis. A TME cell network, depicting a comprehensive landscape of tumor-immune cell interactions with TME scores, was constructed. There was an almost all-round positive correlation between the abundance of the 29 immune-related terms. This phenomenon was most likely due to co-infiltration effects (Figure 5A). Interestingly, the TME scores were negatively correlated with all immune-related terms, indicating that the TME scores reflected the “cold” intratumoral microenvironment of TNBC. The association of TME scores with TNBC in different clusters was also analyzed, and the results showed that the high-infiltration group had a lower score than the low-infiltration group in the training set (Figure 5B). The same trend was also observed in the validation set (Figure 5C). Since the expression of immune checkpoint molecules, such as PD-1, were promising predictive factors for immune treatment response, we explored the association between immune checkpoint molecule expression and TME scores. Interestingly, TME scores were negatively correlated with all four immune checkpoint molecules (PD-1, PD-L1, CTLA-4, and TIM-3; P < 0.05; Figure 5D). To determine whether the patients with high- or low-risk scores were enriched for genes in previously defined biological pathways, GSEA analysis was performed. The results showed that in the low-risk score group, “intestinal immune network for IgA production,” “Th17 cell differentiation,” “B cell receptor signaling pathway,” “T cell receptor signaling pathway,” and “inflammatory bowel disease” were significantly upregulated (Figure 5E). In contrast, in the high-risk score group, “steroid hormone biosynthesis,” “synaptic vesicle cycle,” “biosynthesis of amino acids,” “ECM-receptor interaction,” and “glucagon signaling pathway” were the predominantly upregulated pathways (Figure 5F). Since the 29 immune-related terms contained 16 types of immune cells, we analyzed the association between the TME score and tumor-infiltrating immune cells. Th1 cells, CD8+ T cells, Treg cells, T cells, and 12 other immune-related cell types were shown in Figure 6. A negative correlation was found between the TME score and the tumor-infiltrating immune cells, indicating that the TME scores reflected the intrinsic characteristics of the TME.

Figure 5. Correlation of TME signature with different TME clusters. (A) The correlation between TME signature and infiltrating immune cells. (B) Distribution of TME scores in different TME clusters of training set. (C) Distribution of TME scores in different TME clusters of the validation set. (D) The correlation between TME scores and immune checkpoint molecules. (E) Pathways enriched in low-risk group. (F) Pathways enriched in high-risk group.

Figure 6. Correlation between TME score and Tumor-infiltrating immune cells. The correlation between TME scores with 16 kinds of tumor-infiltrating immune cells.

Correlation of TME signature with clinical parameters

We further investigated the relationship between the TME signature and clinical parameters. The patients were categorized into two groups according to the median of the TME scores. High TME scores were significantly associated with larger tumor sizes (P < 0.001; Supplementary Table 1). Univariate Cox analysis indicated that age, tumor stage, and TME signature were significantly associated with a poor prognosis of TNBC (all P < 0.05; Table 2). Multivariate Cox analysis indicated that age, tumor stage, and TME signature were independent prognostic factors for TNBC (Table 2). We constructed a nomogram based on three independent prognostic factors: age, tumor stage, and TME signature (Figure 7A) to obtain the estimated 3- and 5-year survival probabilities. To assess the calibration of the nomogram, we compared the predicted 3- and 5-year survival probabilities to the actual 3- and 5-year survival probabilities, and found that the calibration curve revealed good concordance between the predicted and observed probabilities (Figure 7B, 7C). These results indicated that the nomogram had proper calibration, and the ROC curve demonstrated that the nomogram had good discrimination with an AUC of 0.88, which was superior to that of the tumor stage (AUC = 0.82) and age (AUC = 0.55) (Figure 7D).

Figure 7. Nomogram construction and evaluation. (A) Nomogram for predicting 3- and 5-year survival probabilities of TNBC. Three points are allocated for age, tumor stage and TME signature. Draw a vertical straight line from the variable value to the axis labeled “Points”. Then calculate all variables’ points. The total points on the bottom scales that correspond to the 3- and 5-y survival were showed apparently. (B, C) Calibration curves for predicting 3-y (B) and 5-y (C) OS. Predicted survival produced by nomogram is plotted on the x-axis, and actual survival is plotted on the y-axis. Dashed lines represent an identical calibration model in which predicted OS approximate to actual OS. (D) Comparison of the predictive accuracy of nomogram, tumor stage and age.

Clinical specimens

A total of 41 specimens of TNBC tissues were obtained from patients at the Guangxi Medical University Cancer Hospital, People’s Republic of China between 2011 and 2019. TNBC was pathologically confirmed in all patients who were not on chemotherapy or radiotherapy before the collection of the tissues. Written informed consent was obtained from all patients, and the study was approved by the Ethics and Human Subject Committee of Guangxi Medical University Cancer Hospital. All experiments and methods were performed according to relevant guidelines and regulations formulated by the Guangxi Medical University Cancer Hospital.

Data acquisition and processing

We retrieved two transcriptomic expression datasets of BC cohorts from public databases: a microarray dataset from Gene Expression Omnibus (GEO; https://www. ncbi.nlm.nih.gov/geo/) database with accession number GSE58812 [30] and a dataset of clinical and mRNA expression from the Molecular Taxonomy of Breast Cancer International Consortium (METABRIC; cBioportal of Cancer Genomics (https://www.cbioportal.org/) [31]. For GSE58812, the probes in the microarray data were annotated according to Affymetrix Human Genome U133 Plus 2.0, batch effects caused by different studies were removed by “ComBat” (empirical Bayes method) from the “sva” package [32, 33], and background noise was cancelled via data normalization using the “limma” package [34]. For the METABRIC cohort, Illumina HT-12v3 platform was used to measure mRNA expression. METABRIC data was formerly median-centered and log-transformed. GEO and METABRIC datasets were utilized for model development and validation, respectively.

IHC analysis

All tissue specimens obtained from the patients were fixed in 10 % neutral formalin and then embedded using paraffin. The embedded specimens were cut into 5 μM serial sections just before the staining. IHC was performed using SPlink Detection Kits SP-9000 (ZSGB-BIO, China) according to the manufacturer’s instructions. The sections were processed by deparaffinization, rehydration, and blocking of endogenous peroxidase activity, and then incubated with primary antibodies (BLK, 1:200; ABCEPTA and TNFRSF1B, 1:200, ALS17735, ABCEPTA) overnight at 4° C. Next, the sections were incubated with a secondary antibody, i.e., a biotinylated donkey anti-goat IgG and horseradish peroxidase-labeled streptomycin working solution, and then treated with diaminobenzidine hydrochloride to visualize the immunoreactivity. The immunohistochemical scoring was performed independently by two authors who were blinded to the clinicopathological parameters. The score was composed of two parts. One was the extent of staining (percentage of positive tumor cells: 0 % = 0; 1–25 % = 1; 26–50 % = 2; 51–75 % = 3; and 75-100 % = 4), and the other was the intensity of staining (no staining = 0; weak staining = 1; moderate staining = 2; and strong staining = 3). The total scores that ranged from 0 to 4 were considered “low expression,” while those that ranged from 5 to 7 were considered “high expression.”

TME cell abundance calculation

To investigate the TME immune cell abundance in TNBC, we methodically examined published studies, adopted the gene signatures proposed by Bindea et al [35], and obtained the marker gene sets for different immune cell types. Gene signatures signifying microenvironment cell subsets of both adaptive and innate immunity were also obtained. TNBC with qualitatively different immune cell infiltration patterns were grouped by hierarchical agglomerative clustering, depending upon Ward's linkage and Euclidean distance. ssGSEA, from the R package “GSVA,” was applied to identify the enrichment scores of each immune-related term to quantify the infiltration levels of different immune cell types in TME [32]. In addition, TIMER2.0 and CIBERSORTx, were employed to validate the TME cell abundance calculated by ssGSEA using default parameters. TIMER2.0 integrates multiple state-of-the-art algorithms for estimating six types of tumor infiltrating immune cells [36]. Meanwhile, CIBERSORTx is a deconvolution algorithm based on support vector regression, which uses a set of reference gene expression values corresponding to a minimal representation for each cell type to infer the cell type proportions in data from bulk tumor samples with mixed cell types [37]. CIBERSORTx can sensitively and specifically discriminate between 22 human immune cell phenotypes.

Unsupervised cluster analysis

We identified characteristic expression patterns based on the infiltration of immune cells by using unsupervised clustering analysis based on Euclidean distance and Ward's linkage. The patients were categorized for further analysis. We established the best number of “TMEclusters” based on the percentage of variance of the data by applying the “ConsensusClusterPlus” package with 1000 repeats [38]. We categorized the patients into two groups based on the infiltrating degrees of the immune cells. Finally, we calculated the correlation matrix using Spearman's test and generated images using R package “ggplot2.”

Identification of DEGs and construction of PPI networks

In order to identify the DEGs between the two distinct TME clusters, we used the R package “limma” [39]. P values were adjusted by Benjamini-Hochberg correction, and DEGs with adjusted P value < 0.05 and |log2 (fold change)| > 1 were retained. PPI networks were constructed using STRING (https://string-db.org/) with confidence of interaction set at 0.15 [40]. The networks were visualized by Cytoscape version 3.6.0 [41].

Enrichment analysis

DEGs were subjected to KEGG and Gene Ontology (GO) pathway enrichment analyses, and only terms with adjusted P < 0.05 were considered significant. The “GOSemSim” R package was used to assess the similarity between the significant GO terms of the two groups by referring to the annotation data GO.db. Similar GO terms between two groups were graphically represented by a heatmap and tree diagram [42]. The median TME signature was used as a cutoff value to classify the patients into high- and low-risk groups. Then, gene set enrichment analysis (GSEA) was performed to test whether the genes in the high- or low-risk groups were enriched in the predefined KEGG gene sets. After 1000 permutations, the gene sets with P < 0.05 and false discovery rate < 0.05 were considered significant.

TME signature construction

Unsupervised clustering was used to categorize the patients into two TME clusters (Cluster 1 and Cluster 2). Based on the DEGs, a univariate Cox regression analysis was performed to identify the prognostic DEGs by applying the “survival” R package. To build the TME signature, DEGs with P value < 0.1 in the univariate Cox regression were considered the prognostic genes. Next, we used a stepwise multivariate Cox regression analysis to screen for key prognostic genes. Then, we computed the TME signature for each patient using the following formula:

TME signature = β1 × exprG1 + β2 × exprG2 + ... βn × exprG

where exprG is the expression level of the vital prognostic genes and β is the regression coefficient computed from the multivariate Cox regression model. The GSE58812 cohort was used as the training set to develop the TME signature, and the METABRIC cohort was used as the validation set to confirm the stability of the TME signature. A forecast value was calculated for the patients in the validation set according to the risk model that was made in the training set. The ROC curve and AUC were used to calculate the predictive risk model’s discrimination ability.

Construction of the nomogram

Nomograms have recently attracted increased attention as a user-friendly tool for predicting prognosis with strong clinical utility [43, 44]. TME signature and clinical parameters were subjected to univariate Cox proportional hazards analyses. Features with P values < 0.05 were subjected to multivariate Cox proportional hazards analysis. Features with P values < 0.05 after multivariate analysis were incorporated into nomograms that were constructed to predict the 3- and 5-year survival rates [45]. The nomogram was based on three independent prognostic factors: age, tumor stage, and TME signature. Each factor corresponded to a specific point by drawing a line straight upwards to the points axis. The sum of the three factor points was defined as total points. By drawing a perpendicular line from the total point axis to the two-outcome axes, the estimated 3- and 5-year survival probabilities were obtained. The observed 3- and 5-year survival rates were compared with the predicted 3- and 5-year survival rates to further verify the predictive performance of the nomogram. We assessed the goodness-of-fit of the nomogram using calibration plots [46].

Statistical analysis

For comparison between two groups, non-parametric Mann-Whitney U test and parametric unpaired Student’s t-test were used. For comparison between more than two groups, non-parametric Kruskal-Wallis test and parametric one-way analysis of variance (ANOVA) test were used. Spearman and distance analysis were used to perform correlation analysis. The survival curves concerning the prognostic analysis were prepared by Kaplan-Meier method, and log-rank tests were employed to identify the significance of the survival differences. We used the "LR forward" stepwise approach and Cox proportional risk model for the univariate and multivariate analyses. To assess the accuracy of survival prediction by the prognostic model, we used a time-related ROC analysis. All statistical analyses were performed in R version 3.5.0. Two-sided P value < 0.05 was regarded as significant.

Ethics approval and consent to participate

This study was approved by the Ethics and Human Subject Committee of Guangxi Medical University Cancer Hospital. All patients provided written informed consent.

Availability of data and materials

The dataset supporting the conclusions of this article is included within the article.