Biological network-driven gene selection identifies a stromal immune module as a key determinant of triple-negative breast carcinoma prognosis

Abstract

Triple-negative breast cancer (TNBC) is a heterogeneous group of aggressive breast cancers for which no targeted treatment is available. Robust tools for TNBC classification are required, to improve the prediction of prognosis and to develop novel therapeutic interventions. We analyzed 3,247 primary human breast cancer samples from 21 publicly available datasets, using a five-step method: (1) selection of TNBC samples by bimodal filtering on ER-HER2 and PR, (2) normalization of the selected TNBC samples, (3) selection of the most variant genes, (4) identification of gene clusters and biological gene selection within gene clusters on the basis of String© database connections and gene-expression correlations, (5) summarization of each gene cluster in a metagene. We then assessed the ability of these metagenes to predict prognosis, on an external public dataset (METABRIC). Our analysis of gene expression (GE) in 557 TNBCs from 21 public datasets identified a six-metagene signature (167 genes) in which the metagenes were enriched in different gene ontologies. The gene clusters were named as follows: Immunity1, Immunity2, Proliferation/DNA damage, AR-like, Matrix/Invasion1 and Matrix2 clusters respectively. This signature was particularly robust for the identification of TNBC subtypes across many datasets (n = 1,125 samples), despite technology differences (Affymetrix© A, Plus2 and Illumina©). Weak Immunity two metagene expression was associated with a poor prognosis (disease-specific survival; HR = 2.68 [1.59–4.52], p = 0.0002). The six-metagene signature (167 genes) was validated over 1,125 TNBC samples. The Immunity two metagene had strong prognostic value. These findings open up interesting possibilities for the development of new therapeutic interventions.

Keywords:

• immune signature
• molecular subtypes
• prognosis
• triple-negative breast cancer

Abbreviations

AR	=	Androgen receptor
BC	=	Breast cancer
BC_CL	=	Breast cancer cell lines
BCSS	=	Breast cancer-specific survival
CCLE	=	Cancer Cell Line Encyclopedia
CDF	=	Consensus distribution function
CGP	=	Cancer Genome Project
EMT	=	Epithelial-mesenchymal transition
ER	=	Estrogen receptor
GE	=	Gene expression
HER2	=	Human epidermal growth factor receptor 2
IHC	=	Immunohistochemistry
IM	=	Immunomodulatory
LAR	=	Luminal androgen receptor
M	=	Mesenchymal
MSL	=	Mesenchymal stem-like
no CT	=	No chemotherapy
NPI	=	Nottingham Prognostic Index
pCR	=	Pathological complete remission
PR	=	Progesterone receptor
RMA	=	Robust multichip average
TILs	=	Tumor-infiltrating lymphocytes
TNBC	=	Triple-negative breast cancer
TNBC_CL	=	Triple-negative breast cancer cell lines.

Introduction

TNBC, defined by the absence of estrogen and progesterone receptor expression and a lack of HER2 overexpression/amplification, is an aggressive disease accounting for 15%–20% of breast cancers. It differs from other molecular subtypes ^1–3 in displaying axillary lymph node involvement, local and regional recurrence, differences in the time lag to metastasis (distant metastatic events occurring within 5 y of diagnosis), high rates of brain, lung and distant nodal metastasis and in its response to neoadjuvant treatment.

TNBC constitutes a major clinical challenge because there has been no substantial improvement in treatment for this subgroup in the recent past. Even if adjuvant chemotherapy has significantly improved outcome, reducing the risk of death by approximately 30%,⁴ but these cancers do not respond to endocrine or targeted therapy. TNBC is, thus, currently the breast cancer subgroup with the worst outcome.⁵ Moreover, the shape of the survival curve for this subgroup differs from that for other BC subtypes: there is a sharp decrease in survival during the first 3–5 y after diagnosis, but distant relapses, occurring after this interval, are much less common.⁵

TNBC is a highly heterogeneous group of tumors differing in terms of their histological features, GE profiles, clinical behavior, overall prognosis ⁶ and sensitivity to systemic treatment.^7-9

Robust classifiers are urgently required, to improve our understanding of the molecular basis of TNBC and to define novel therapeutic interventions. Lehmann et al. recently published a classification of six molecular subtypes of TNBC¹⁰ and developed a website (http://cbc.mc.vanderbilt.edu/tnbc/) ¹¹ for the classification of TNBC samples on the basis of their GE profiles. This classification has been shown to be relevant, as it identifies the main biological component and pathways of TNBC. However, the large number of genes defining this TNBC molecular classification (2,188 genes) constituted a potential source of instability.^12,13

We developed a two-step biological network-driven gene selection process: (1) identification of the most variant genes displaying highly-correlated patterns of expression, (2) direct connection of these genes within known biological networks. This method has been reported to be efficient for the construction of molecular signatures.^14,15 We defined a robust TNBC molecular subtype classification, providing considerable biological insight, with great potential for use in the development of therapeutic interventions. We also identified a stromal immune module GE profile strongly correlated with TNBC prognosis.

Results

TNBC gene expression profiles identify six main gene clusters

GE profiles were obtained from 21 publicly available datasets, containing data for 3,247 primary human breast cancer samples. These profiles were processed according to the flow chart in Fig. 1. The training set included samples hybridized on HGU-133A Affymetrix© arrays (12 datasets, n = 1,995), to eliminate cross-platform discrepancies and to ensure robust normalization. The validation set included samples hybridized on HGU-133Plus2 Affymetrix© arrays (9 datasets, n = 1,014). We filtered out 42 outlier samples from the training set and 17 from the validation set.

Figure 1. Methodology flow chart.

We also collected two large datasets, for the validation of our classification: the Ignatiadis set (n = 996) and the METABRIC set (n = 1,992). The processing of these two datasets has been described elsewhere.^16,17

Bimodal filtering on ER-PR and HER2 GE identified 262, 295, 314 and 254 TNBC samples in the training set, the validation set, the Ignatiadis set and the METABRIC set, respectively.

We developed a gene selection process based on biological networks, to decrease the intrinsic instability of molecular classification methods.

We identified the 830 most variant genes (SD > 0.8) in the training set (n = 262). A consensus clustering method and hierarchical clustering identified four main gene clusters. Further increases in cluster number yielded no significant increase in the consensus distribution function (CDF) area (Fig. S1 and Materials and Methods).

The various gene clusters were associated with different gene ontologies (Fig. S2). The clusters were thus named as follows (Fig. S3A): Immunity cluster (145 genes), Proliferation/DNA damage cluster (397 genes), AndrogenReceptor(AR)-like cluster (139 genes) and Matrix/Invasion cluster (149 genes).

The Immunity cluster was the most homogeneous, with strong correlations between the GE profiles of most of the genes within this cluster (Fig. S3B).

We used String© database software to analyze our gene selection, with the aim of decreasing the heterogeneity of each main gene cluster. We retained the genes from our initial selection that (1) had high String© database gene connection indexes (greater than 0.7, Fig. S4), (2) had similar patterns of expression to other genes within the same biological network (correlation coefficient of at least 0.5). We selected a final set of 167 genes [Immunity cluster (80), Proliferation/ DNA damage (15), AR-like(15), Matrix/Invasion (57)] (Fig. S5).

Following biological network-driven gene selection, it became clear that the original Immunity and Matrix/Invasion clusters were more accurately described by splitting them into two subclusters displaying minor differences [Immunity1 (33), Immunity2 (47), Matrix/Invasion1 (43), Matrix2 (14)] (Fig. S6A). This approach yielded an increase in the area under the CDF curve (Fig. S7).

For each of the six gene clusters identified in this way, we defined a metagene. The Immunity1 and Immunity2 metagenes displayed similar patterns of expression, with a Pearson correlation coefficient of 0.58; the Pearson correlation coefficient for the expression patterns of Matrix/Invasion1 and Matrix2 was 0.48. The Proliferation/DNA damage and Matrix metagenes displayed the strongest inverse correlation (coefficients of −0.43 and −0.60 for Matrix/Invasion1 and Matrix2, respectively) (Fig. S6B).

We validated this six-gene cluster classification, by applying hierarchical clustering based on the 167 genes selected to the validation set (n = 295). Clustering was highly consistent between the training and validation gene sets (concordance: 93–100%).

The six gene clusters identify six stable TNBC subgroups

Hierarchical clustering was performed on the four TNBC datasets [training set (262), validation set (295), Ignatiadis (314) and METABRIC (254)]. For Affymetrix© arrays, we used the 167 selected genes. For the Illumina© platform, we used 153 common genes. We identified six reproducible subgroups of TNBC, for which GE patterns were similar in the training set and in the three validation sets (total of 1,125 samples). The corresponding heatmaps are shown in Fig. 2. The Pearson correlation coefficients for the relationships between each sample subgroup centroid in the three validation sets and the corresponding subgroup centroid in the training set are shown in Fig. 2.

Figure 2. Heatmaps of the selected genes in the TNBC training set (upper left) and the TNBC validation sets (upper right: validation, lower left: Ignatiadis, lower right: METABRIC).

We illustrated the dynamic links between genes within a biological network, as defined by the String© database, by showing GE levels for a “prototype sample” (Fig. S8).

We compared our sample classification with those reported by Lehmann et al. and Curtis et al. (Fig. S9). Our classification appears very different from that of Lehmann at first glance (χ² test p value = 0.05), but the samples assigned to Centroids one and six (with high-levels of Matrix/Invasion 1 and Matrix 2 gene expression, respectively) tended to be classified as Mesenchymal (M) or Mesenchymal stem-like (MSL), the samples in Centroid 5 (strong expression of Immunity2 genes) tended to be classified as Immunomodulatory (IM), and the samples in Centroid 4 (strong expression of AR-like genes) tended to be classified as of the Luminal androgen receptor (LAR) subtype (Fig. S10A and Fig. S10B). Curtis et al. aimed at defining a new classification across all cancer subtypes, not specific to TNBC subtypes. In this classification, the TNBC samples were mostly classified as IntClust10 or IntClust4, with an even distribution.

Prognostic value of the Immunity2 metagene in TNBC

The prognostic value of the 167-gene TNBC signature was assessed with the METABRIC dataset. The 254 TNBC samples were split into two subgroups: a subgroup treated by chemotherapy (n = 139) and a subgroup not treated by chemotherapy (n = 115). The chemotherapy-naive (noCT) population and the chemotherapy-treated population were significantly different (Table S1). The patients in the noCT population were older (mean age of 61.5 y vs. 50.1 y, p < 1.2^10-11), more likely to be postmenopausal (77% vs. 47%, p = 5.38^10-5), and their tumors were of lower grade (p = 0.01), with less lymph node involvement (81% vs. 17%, p < 2.2 ^10-16), a lower Nottingham Prognostic Index (NPI < 3.4, 17% vs. 2%, p = 2.57^10-5), and less cellularity (p = 0.03).

Univariate analysis identified three factors significantly correlated with a poor outcome (distant disease-free survival) in the chemotherapy-treated population: NPI > 5.4 (HR = 2.15 [1.28–3.60], p = 0.003); p53 mutation (HR = 2.42 [1.15–5.09], p = 0.02); and weak Immunity2 metagene expression (HR = 2.59 [1.54–4.34], p = 0.0002) (Table 1A, Fig. 3A). We did not include p53 mutation status in the multivariate model, due to missing data (n = 79). A NPI > 5.4 and low-levels of Immunity2 metagene expression were retained in the multivariate model and were significantly associated with a poor outcome (HR = 2.30 [1.36–3.89], p = 0.002; HR = 2.68 [1.59–4.52], p = 0.0002, respectively) (Table 1A). The combined variable, NPI score/Immunity2 metagene expression was found to be of particular interest. In a first model, a NPI score greater than 5.4 was associated with a worse prognosis: HR = 3.98 [2.00–7.92], p = 8.72^10–5. For patients with NPI scores of 5.4 or below, Immunity2 metagene expression discriminated between two groups of patients with different outcomes (HR = 2.90 [1.51–5.56], p = 0.001). In a second model, NPI3 patients can also be split into two groups on the basis of Immunity2 metagene expression. The NPI3 group with high-levels of Immunity2 metagene expression had a prognosis similar to that of the NPI1/2 group with low-levels of Immunity2 metagene expression (Table 1B, Fig. 3A).

Table 1 A. Survival analysis (disease-specific survival). Chemotherapy-treated population. Univariate and multivariate analysis.

		Univariate analysis		Multivariate analysis
139 triple-negative breast cancer patients		DS-survival HR [95% CI]	p value	DS-survival HR [95% CI]	p value
Menopausal status	Pre	1
	Post	1.56 [0.95–2.55]	0.08
Tumor size (mm)	<20 mm	1
	>20 mm	1.03 [0.58–1.82]	0.92
Tumor grade	II	1
	III	1.23 [0.45–3.39]	0.69
Lymph node status	0	1
	1	0.84 [0.42–1.65]	0.61
NPI score	<5.4	1		1
	>5.4	2.15 [1.28–3.60]	0.003	2.30 [1.36–3.89]	0.002
Cellularity	Low	1
	Moderate	0.57 [0.22–1.46]	0.24
	High	0.59 [0.25–1.39]	0.23
P53 status	Wild-type	1
	Mutant	2.42 [1.15–5.09]	0.02
Immunity1 metagene expression	High	1
	Low	0.97 [0.60–1.58]	0.91
Immunity2 metagene expression	High	1		1
	Low	2.59 [1.54–4.34]	0.0002	2.68 [1.59–4.52]	0.0002
Proliferation/DNA damage metagene expression	High	1
	Low	1.13 [0.69–1.84]	0.63
AR-like metagene expression	High	1
	Moderate	1.07 [0.59–1.94]	0.82
	Low	0.98 [0.50–1.94]	0.96
Matrix/Invasion1 metagene expression	High	1
	Low	1.23 [0.76–2.01]	0.40
Matrix2 metagene expression	High	1
	Low	0.99 [0.61–1.61]	0.96

Abbreviations: NPI, Nottingham Prognostic Index; AR, androgen receptor; HR, hazard ratio; CI, confidence interval.

Table 1 B. Survival analysis (disease-specific survival). Chemotherapy-treated population. Two univariate models. Combination of NPI score and Immunity2 metagene expression

139 triple-negative breast cancer patients		DS-survival HR [95% CI]	p value
NPI score/Immunity2 metagene expression	NPI1-2/HighI2	1
	NPI1-2/LowI2	2.90 [1.51–5.56]	0.001
	NPI3	3.98 [2.00–7.92]	8.72 ^10-5
NPI score/Immunity2 metagene expression	NPI1-2/HighI2	1
	NPI1-2/LowI2	2.91 [1.51–5.59]	0.001
	NPI3/HighI2	2.31 [0.96–5.57]	0.06
	NPI3/LowI2	6.30 [2.89–13.78]	3.87^10-6

Abbreviations: NPI, Nottingham Prognostic Index; I2, Immunity2; HR, hazard ratio; CI, confidence interval.

NPI = [0.2 × S] + N + G.

S: tumor size (cm).

N: number of lymph nodes involved (0 = 1, 1−3 = 2, >3 = 3).

G: tumor grade according to Elston and Ellis (Grade I=1, Grade II = 2, Grade III = 3).

Figure 4. String Software connections between our Immunity1 and Immunity2 genes and the genes of eight previously published prognostic immune signatures. Stronger associations between genes are represented by thicker lines. Associations between genes with a coefficient < 0.9 are shown in green. Associations between genes with a coefficient ≥0.9 are shown in red. Associations between genes with a coefficient between 0.4 to 0.7 are not shown.

Univariate analysis identified four factors significantly correlated with poor outcome in the noCT population: tumor size >20 mm (HR = 2.36 [1.01–5.48], p = 0.04), lymph node-positive status (HR = 3.66 [1.65–8.11], p = 0.001), NPI score >5.4 (HR = 10.69 [2.74–41.76], p = 0.001) and low-levels of Immunity2 metagene expression (HR = 2.33 [1.09–4.95], p = 0.03) (Table 2A, Fig. 3B). Two of these factors were retained in the multivariate model: NPI score >5.4 (HR = 12.03 [3.05–47.50], p = 0.0004) and low-levels of Immunity2 metagene expression (HR = 2.42 [1.13–5.16], p = 0.02) (Table 2A). As in the chemotherapy-treated subpopulation, the combined variable, NPI score/Immunity2 metagene expression discriminated between two groups of patients with different outcomes in this noCT population (Table 2B, Fig. 3B). The chemotherapy-naive group contained only seven patients classified as NPI3. Stratification of this subgroup defined on the basis of treatment was therefore not considered methodologically relevant.

Table 2 A. Survival analysis (disease-specific survival). Chemotherapy-naive population. Univariate and multivariate analysis.

		Univariate analysis		Multivariate analysis
115 triple-negative breast cancer patients		DS-survival HR [95% CI]	p value	DS-survival HR [95% CI]	p value
Menopausal status	Pre	1
	Post	1.31 [0.56-3.06]	0.53
Tumor size (mm)	<20 mm	1
	>20 mm	2.36 [1.01-5.48]	0.04
Tumor grade	I-II	1
	III	1.33 [0.51-3.49]	0.56
Lymph node status	0	1
	1	3.66 [1.65-8.11]	0.001
NPI score	<3.4	1		1
	3.4-5.4	1.36 [0.47-3.96]	0.57	1.55 [0.53-4.51]	0.43
	>5.4	10.69 [2.74-41.76]	0.001	12.03 [3.05-47.50]	0.0004
Cellularity	Low	1
	Moderate	1.91 [0.54-6.71]	0.31
	High	1.42 [0.41-4.90]	0.58
P53 status	Wild-type	1
	Mutant	0.90 [0.17-4.63]	0.90
Immunity1 metagene expression	High	1
	Low	1.56 [0.76-3.19]	0.22
Immunity2 metagene expression	High	1		1
	Low	2.33 [1.09-4.95]	0.03	2.42 [1.13-5.16]	0.02
Proliferation/DNA damage metagene expression	High	1
	Low	1.14 [0.56-2.32]	0.72
AR-like metagene expression	High	1
	Moderate	0.96 [0.42-2.20]	0.92
	Low	0.74 [0.28-2.00]	0.56
Matrix/Invasion1 metagene expression	High	1
	Low	0.48 [0.23-1.01]	0.06
Matrix2 metagene expression	High	1
	Low	1.31 [0.64-2.66]	0.46

Abbreviations: NPI, Nottingham Prognostic Index; AR, androgen receptor; HR, hazard ratio; CI, confidence interval.

Table 2 B. Survival analysis (disease-specific survival). Chemotherapy-naive population. Univariate analysis. Combination of NPI score and Immunity2 metagene expression.

115 triple-negative breast cancer patients		DS-survival hazard Ratio [95% CI]	p value
NPI score/Immunity2 metagene expression	NPI1-2/HighI2	1
	NPI1-2/LowI2	2.13 [0.95-4.78]	0.07
	NPI3	12.89 [4.07-40.82]	1.37^10−5

Abbreviations: NPI, Nottingham Prognostic Index; I2, Immunity2; HR, hazard ratio; CI, confidence interval.

We compared the prognostic value of the Immunity2 metagene with that of eight previously published immune signatures,^18-25 using the METABRIC dataset.

We generated a heatmap (Fig. S11) of the GE profiles of each of the above prognostic signatures applied to the METABRIC dataset. The samples were ordered according to our classification of low/high Immunity2 metagene expression. Expression patterns were very similar between the Immunity2 GE signature and all the other GE signatures, with the exception of the Bianchini, Karn and Burstein (BLIS) gene-expression signatures.

We first performed a univariate analysis of the prognostic value of the eight-GE signatures, as described in the corresponding original manuscripts. The Rody, Sabatier, Teschendorff, Desmedt, Gu-Trantein Tfh, Gu-Trantien Th1 and Burstein signatures were significantly correlated with the prognosis of TNBC. The Bianchini and Karn GE signatures were not correlated with the prognosis of TNBC (Fig. S12, Table S2). We then performed a multivariate analysis. We included NPI score, the Immunity2 metagene and each of the Rody, Sabatier, Teschendorff, Desmedt, Gu-Trantein Tfh, Gu-Trantien Th1, and Burstein signatures, one-by-one, in the model. In all comparisons the only significant variables remaining in the multivariate model were NPI score and the Immunity2 metagene (Table S2).

The Immunity2 metagene corresponds to B-cell and T-cell pathways

String database connections between the Immunity1 or Immunity2 genes and the genes of the eight published prognostic immune signatures ^18-25 are provided in Fig. 4. The gene intersection was poor, but our immune signature nevertheless appears to be strongly correlated with other published signatures (Supplementary data), suggesting the use of similar immune pathways. The Immunity2 metagene was strongly correlated with the expression metagenes of the above signatures (coefficient greater than 0.8), except for the Bianchini, Karn and BLIS metagenes (Fig. S13).

Figure 3. (A) Kaplan–Meier plots. Disease-specific survival of the chemotherapy-treated population (n = 139). NPI score. Immunity2 metagene. NPI score/Immunity2 metagene. (B) Kaplan–Meier plots. Disease-specific survival of the noCT population (n = 115). NPI score. Immunity2 metagene. NPI score/Immunity2 metagene.

We explored the pathways relating to the immune metagenes in detail, by analyzing the correlation between the expression of the Immunity1 and Immunity2 metagenes and the metagenes defined by Gatza et al.²⁶ (IFN-α, IFNγ, STAT3, TGF-β, TNF-α) and Palmer et al.²⁷ (LB, LT, CD8⁺, GRANS, LYMPHS). This analysis was performed on the METABRIC dataset published by Curtis et al.¹⁷

We showed that the Immunity2 metagene was highly correlated with the B-cell, T-cell and CD8⁺ cell metagenes (Pearson correlation scores: 0.93, 0.91, 0.87, respectively) (Fig. S14). The Immunity1 metagene was highly correlated with the interferon alpha and gamma pathways (Pearson correlation scores: 0.97, 0.94, respectively).

Furthermore, in cancer cell lines (CCLE and CGP datasets), the Immunity2 metagene displayed very low-levels of expression, similar to those of the CD8⁺ metagene (Fig. S15). This was true for all cell lines and BC_CLs tested.

Moreover, the IFN-y, IFN-gamma, STAT3, TGF-β, TNF-α, LB, LT, GRANS metagenes were more strongly expressed in TN BC_CLs than in HER2-positive and luminal BC_CLs (Fig. S16).

We investigated Immunity2 GE in white blood cell populations (Palmer et al.²⁷), by performing a consensus clustering of the Immunity2 genes on Palmer's dataset. This analysis identified four stable clusters of the genes of the Immunity2 signature. Some genes were more strongly expressed in B cells (GZMA, GZMB, CCR7, LY96, MS4A1, CD74 for example), others in T cells (CD3D, CCL2, CD14, CD2, LCK, IL7R), and still others in granulocytes (named Pax cells) (CD163, MNDA, NCF2, CSF2RB, FGL2) (Fig. S17). These findings suggest that, even if the “Immune2” signal is highly homogeneous within tumor samples (the entire set of genes being coordinately either over- or under-expressed), different subpopulations of cells express different subsets of these genes in the periphery.

The Immunity2 metagene is probably expressed by stromal cells

In TNBC cell lines (TNBC_CL), genes from the Immunity2 module displayed very low medians and narrow ranges of expression, suggesting that they were expressed only in the tumor stromal compartment. A similar trend was observed for all BC_CLs. The Immunity1 module genes had higher median expression levels and a broader range of expression in TNBC_CL and in all BC_CL, suggesting that Immunity1 genes were expressed by the tumor cells (Figs. 5A and B).

Figure 5. (A) Boxplots of gene expression for the Immunity1 and Immunity2 metagenes, in each breast cancer cell line subtype from the CCLE. (B) Boxplots of gene expression for the Immunity1 and Immunity2 metagenes in each breast cancer cell line subtype from the CGP.

Furthermore, we explored the contributions of stromal and cancer cells to Immunity1 and Immunity2 expression in detail, by comparing our gene lists to the “stromal contribution to global GE evaluated in PDX RNAseq data”, as defined by Isella et al.²⁸ The Immunity2 metagene had a very high stromal fraction, as for the Matrix/Invasion1 and Matrix2 metagenes. The Immunity1 metagene had a very low stromal fraction, like the AR-like and Proliferation/DNA damage metagenes (Fig. S18).

The Immunity2 metagene opens up interesting new possibilities for therapeutic interventions

To highlight the new opportunities for therapeutic intervention provided by this study, we represented the existing drugs (with or without US Food and Drug Administration approval) for each metagene (Fig. S19 and Supplementary data). Some are undergoing clinical investigation in patients with TNBC.

We explored the links between PD1, PDL1, CTLA4 (and their respective metagenes) and the Immunity2 metagene. We compared the Immunity2 metagene with the TILs signature defined by Schalper et al.,²⁹ who showed that PD-L1 mRNA synthesis was associated with increases in the expression of TILs and recurrence-free survival. This analysis was performed on the METABRIC dataset. The PD1 and CTLA-4 metagenes were constructed from the genes most strongly correlated with the PD1 and CTLA-4 genes, respectively (Pearson correlation score >0.8). The PDL1 metagene was defined by Sabatier et al.³⁰

The Immunity2 metagene was highly correlated with the PD1, PDL1 and CTLA-4 metagenes (Pearson correlation: 0.90, 0.96, 0.91, respectively). The coefficient of correlation between the Immunity2 metagene and the TILs signature was up to 0.90 (Fig. S20).

In cell lines, the PD1, PDL1, CTLA-4 and TILs metagenes were very weakly expressed, like the Immunity2 metagene (Fig. S21).

Using the METABRIC dataset, we compared the prognostic value of these metagenes (PD1, PDL1, CTLA4 and TILs) with that of the Immunity2 metagene. In univariate analysis, high-levels of PD1, PDL1, CTLA-4 and TILs metagene expression were associated with a good prognosis (Fig. S22, Table S3). In multivariate analysis, we included NPI score, the Immunity2 metagene and each of the PD1, PDL1, CTLA4 and TILs metagenes, one-by-one, in the model. In all comparisons, the only significant variables remaining in the multivariate model were NPI score and the Immunity2 metagene.

Discussion

New tools for classifying TNBCs are urgently required, to improve our understanding of the molecular basis of TNBC and to identify potentially useful novel therapeutic interventions. By analyzing the GE profiles of 1,125 TNBCs, we identified a six-metagene signature (167 genes) in which the various metagenes were enriched in different gene ontologies: two clusters were enriched in immunity genes, one in proliferation/DNA damage genes, one in AR pathway genes, and two in matrix/invasion genes. This signature appeared to be particularly robust for identifying TNBC subtypes across different datasets, independently of the gene chip technology used to generate the data. Furthermore, one metagene (Immunity2) was found to be of strong prognostic value for TNBC samples.

Lehmann et al.¹⁰ recently developed a classification of TNBCs in which a 2,188-gene signature was used to classify tumors. They suggested that this classification could also be used to classify xenografts or cell lines. They also developed a website (http://cbc.mc.vanderbilt.edu/tnbc/) for the classification of TNBC samples.¹¹ This study provided important biological insight into the molecular drivers of TNBC, but it also raised several key concerns. First, the normalization process involved data from different platforms. Several studies have shown that large discrepancies in signature composition and absences of concordance concerning outcome may be due to differences in the array platform and preprocessing method used.¹² Second, Lehmann et al. used a very large number of genes (2,188 genes) to establish their molecular signature, and this may have constituted a source of instability, due to the noise introduced.^12,13 As shown by Weigelt et al.,³¹ microarray-based single-sample predictors do not allocate individual samples to a given molecular subtype reproducibly, probably because the use of large numbers of genes leads to instability of the classification when new samples are added. Third, it would be unwise to transpose this classification to various in vitro and in vivo breast cancer models (primary tumor xenografts, cell lines, cell line-derived xenografts), because the stromal environment and the original tumor are very different.^32,33 We found that genes from the Immunity compartment (Immunity2 module) were highly relevant for the classification of TNBC samples and that these genes were not expressed in BC_CLs. The observed lack of reproducibility between classifiers may reflect major differences in the methodology and aims of the studies concerned. Further validation will be required before these models can be used in routine clinical practice.

We developed a strategy for the definition of a GE signature based on the analysis of biological networks for the most variant genes. Within these networks, we then analyzed GE parameters, to select the genes with the most strongly correlated patterns of expression. The validation process showed that our gene matrix identified similar GE patterns across 1,125 TNBC samples. This first step in biological network analysis, which is probably less sensitive to sample fluctuations than other methods, made it possible to capture strong biological signals that might be concealed by the noise present in microarray data. Several studies have reported that the incorporation of network information improves the stability of gene selection and the biological interpretability of biomarker signatures for a given prediction accuracy.^14,15,34

The Immunity2 module was identified as a strong prognostic factor for disease-specific survival (strong expression of this metagene is correlated with a good outcome), regardless of the characteristics of the tumor (NPI score, tumor size, tumor grade and lymph node status). It clearly suggest the presence of an hemopoietic infiltrate, composed of activated cytotoxic T cells, B cells, myeloid cells, natural killer cells and neutrophils. This module includes adhesion molecule-associated genes (SELL, ITGB2), and genes encoding proteins involved antigen processing and presentation (CD74 or ligand, HLA-DRA), B-lymphocyte cell surface molecules (PTPRC, ITGB2, HLA-DRA), the caspase cascade (CASP1), complement pathway (C1QA, C1QB), CTL-mediated immune responses to target cells (ITGB2, CD3D, GZMB), dendritic cell regulation of Th1 and Th2 development (CD2, IL7R), granzyme-mediated apoptosis (GZMA, GZMB), IL12-mediated signaling events (CD3D, HLA-DRA, GZMA, LCK), the IL2 signaling pathway (LCK), interleukin-3, 5 and GM-CSF signaling (HCK, BLNK, CSF2RB), T-cell surface molecules (PTPRC, CD3D, CD2, ITGB2), and the T-cell receptor signaling pathway (PTPRC, CD3D, HLA-DRA, LCK).

Burstein et al.²⁵ identified four different TNBC subtypes (LAR, MES, BLIS, BLIA) with the identification of similar pathways and a prognostic value for the BLIA subgroup similar to that for the signature identified in our study. This subgroup displays an upregulation of B-cell, T-cell, and natural killer cell immune-regulating pathways and an activation of STAT transcription factor-mediated pathways. The authors showed that the prognosis was worse for basal-like immune-suppressed tumors than for basal-like immune-activated tumors, for both disease-free survival (p = 0.04) and disease-specific survival (p = 0.039).

Several recent studies have demonstrated the importance of tumor-infiltrating lymphocytes (TILs) in controlling the clinical progression of various epithelial cancers.³⁵ In breast cancer, recent advances in GE profiling have revealed an association between immune signatures and favorable outcomes.^29,36 A gene signature enriched in cytotoxic CD8⁺ T-cell genes and genes associated with natural killer cell activity has been reported.³⁷ However, the ability of CD8⁺ T cells to control human breast cancer is probably counteracted by the presence of immunosuppressive cells, CD4⁺ T-regulatory cells or macrophages: immunohistochemistry (IHC) analysis of tissue microarray data for 179 treatment-naive breast tumors revealed that high-levels of macrophages and CD4⁺ T cells were correlated with poor overall survival, whereas a combination of high-levels of CD8⁺ T cells and low-levels of macrophages and CD4⁺ T cells was correlated with higher overall survival.³⁸ Intratumoral B cells have also been associated with a favorable prognosis in breast cancer.³⁹ In ER-negative breast cancers, a STAT1 signaling metagene,¹⁶ and a B-cell metagene¹⁹ were found to be associated with better outcomes. Another group identified an immune response-based prognostic gene module (C1QA, XCL2SPP1, TNFRSF17, LY9, IGLC2, HLA-F) associated with a better prognosis than for other ER-negative breast cancers, regardless of lymph node status and lymphocytic infiltration.⁴⁰ According to Bertucci et al.,⁴¹ the IM subtype (overlapping with medullary breast cancers, a rare form of TNBC with a prominent lymphocytic reaction) is associated with a favorable prognosis. The two immune modules identified in this study had many biological connections with other eight immune prognosis signatures published for TNBC.^18-25

Neoadjuvant chemotherapy is increasingly being used for TNBC, because these tumors have a poor prognosis, are assumed to be chemosensitive and no alternative specific systemic treatment is available. Patients with a complete pathologic response (pCR) after neoadjuvant chemotherapy have a better outcome than those with residual disease, and pCR is a good surrogate for long-term survival and cure in this specific subgroup.^9,42

The Immunity2 metagene was not found to be predictive of response to neoadjuvant chemotherapy in TNBC (272 fine needle aspirations of TNBC samples for which information about pCR or its absence was available from the eight datasets previously published by Ignatiadis et al. ¹⁶) (data not shown). This lack of relationship may have resulted from the use of fine needle aspiration biopsy samples. The Immunity2 genes, which are largely expressed in the stromal environment, were less strongly expressed in fine needle aspiration samples than in tumor samples (Fig. S23).

However, intratumoral immune responses are known to be correlated with clinical outcomes in TNBC. This may reflect the role of immune cells in the activity of cytotoxic chemotherapeutic agents. Some chemotherapeutic drugs, such as anthracyclines, act not only through direct cytotoxic effects, but also by activating CD8⁺ T-cell responses. Conflicting results have been published on the ability of other immune-based classifiers to predict outcome in TNBC. High-intratumoral levels of CD8⁺ T cells ⁴³ or TILs ^36,44 are associated with better clinical responses to anthracycline-based chemotherapy. West et al. ⁴⁵ reported that high-levels of lymphocyte GE were associated with a high rate (74%) of complete pathological responses to neoadjuvant anthracycline-based chemotherapy. In 2011, Sabatier et al. ²⁰ showed, by gene-expression profiling, that “Immune High” patients (59%) were more likely to present pCR than “Immune Low” patients (43%), but this difference was not significant (p = 0.29). In 2014,⁴⁶ they showed that “PDL1 mRNA expression high” (57%) patients presented higher rates of pCR than “PDL1 mRNA expression low” (43%) patients (p < 0.001). Wimberley et al. ⁴⁷ showed that PDL1 protein levels in the epithelium and stroma were correlated with pCR only in hormone receptor-positive and HER2-amplified breast cancers. Denkert et al. ⁴⁴ demonstrated the importance of TIL and immune GE signatures for predicting pCR in breast carcinoma. However no significant difference in pCR rate was detected between lymphocyte-predominant breast cancer (LPBC) and no-LPBC in the anthracycline-taxane subgroup.

However, the results of these studies suggest that clinical outcomes in ER-negative breast cancers, including TNBC in particular, are strongly influenced by tumor immune responses and are, thus, highly responsive to immunotherapies. The possible use of immunotherapy approaches to treat TNBC (tumor vaccine approaches, immune-checkpoint inhibitors, antagonists of immunosuppressive molecules and adoptive cell therapies) should be investigated in detail.⁴⁸

The other metagenes studied had no significant prognostic or predictive value. However, they identified sound biological networks providing opportunities for therapeutic intervention. The Immunity1 metagene included genes involved in the interferon α/β signaling pathway or cytokine signaling (STAT1, IRF7, IRF27, OAS1, OAS2, PMSB8, XAF1, IFIT1, IFITM1, ISG15, IGS20, IF6, MX1), the Toll-like receptor signaling pathway (STAT1, CXCL9, CXCL10, CXCL11, CCL5, IRF7), cell-cycle checkpoints and DNA synthesis (PSMB8, PSMB9). Patients displaying strong expression of this metagene often also had high-levels of Immunity2 and Proliferation/DNA damage metagene expression, suggesting the possible existence of common pathways. The IDO1 (indoleamine 2, 3-dioxygenase 1) gene is a particularly interesting potential target. It encodes a tryptophan-degrading enzyme known to suppress antitumor CD8⁺ T cells and it contributes to the inhibition of anticancer immune responses.⁴⁸ This immunosuppressive enzyme is actually investigated as a promising candidate target in cancer immunotherapy.

A subset of TNBC tumors strongly expresses AR-regulated genes.⁴⁹ AR expression has been reported to be lower in triple-negative breast tumor cells than in other types of breast cancer. The overall frequency of AR expression in carcinoma cells varies considerably between studies (0–53%).^50,51 We identified strong expression of AR pathway genes in 25% of our population. The biological role of androgens in TNBC remains a matter of debate. Immunohistochemical studies investigating the presence of AR in tumor cells have reported conflicting results for clinical outcome; some studies have suggested that AR expression is advantageous for survival,^52-54 whereas others found no significant effect.⁵⁵ Lehmann et al. found that the LAR subtype of TNBC displayed the lowest frequency of pCR (10%). The presence of AR in a subset of TNBC patients suggests that androgenic pathways in tumor cells could be targeted in at least some TNBC patients. The widespread availability of agents targeting the AR also makes this approach potentially appealing, as it would be straightforward to incorporate such treatment into clinical practice.

The Matrix/Invasion1 metagene included genes associated with β1 integrin cell surface interactions, ECM-receptor interaction or integrin family cell surface interactions (NID1, TGFBI, COL5A1, COL5A2, COL6A3, COL3A1, COL1A1, COL1A2, COL11A1, FN1, FBN1, THBS1, THBS2), the TGF β signaling pathway (DCN, COMP, THBS1, THBS2), the inhibition of matrix metalloproteases (MMP2, TIMP3), and the AP-1 transcription factor network (DCN, COL1A2, MMP2). Metalloproteinases (MMPs) and their tissue inhibitors are involved in several key pathways of tumor growth, invasion and metastasis.^56,57 The expression and activity of MMPs has been linked to advanced stages of breast cancer, greater tumor invasion and the construction of metastatic formations.^58,59,60 Some studies have highlighted the importance of matrix MMP expression by stromal cells as a prognostic factor in the TNBC subtype.⁶¹ These molecules are thus attractive targets for drug development.⁶²

The Matrix2 metagene included genes associated with the AP-1 transcription factor network (FOS, EGR1, FABP4, DUSP1), the EGR receptor signaling pathway (FOS, DUSP1, EGR1), the Wnt or ALK signaling pathway (CAV1), the MAPK signaling pathway (FOS, DUSP1) or Trk receptor signaling mediated by the MAPK pathway (FOS, EGR1), the mTOR signaling pathway (IGF1), the PPAR signaling pathway (ADIPOQ, CD36, FABP4), and androgen-mediated signaling (FOS, EGR1). These pathways may contribute to cell motility and tumor cell invasion⁶³ and play a prominent role in epithelial-mesenchymal transition (EMT) and in stem cells. These metagenes are strongly expressed in mesenchymal cells and metaplastic breast cancers.⁴ Metaplastic breast cancers have lineage plasticity, including spindle cell foci, and display osseous or cartilaginous differentiation.⁶⁴ Some drugs targeting the pathways relating to the metagenes identified here may be of particular interest for the treatment of TNBC (PI3K/mTOR inhibitor, Wnt/β catenin inhibitor).

Conclusion

In conclusion, our 167-gene TNBC molecular signature, consisting of six metagenes, appears to be particularly robust for the identification of TNBC subtypes. Furthermore, expression of the Immunity2 metagene was strongly correlated with prognosis, and many biological targets have been identified within the corresponding biological network. These findings open up interesting new possibilities for the development of new therapeutic interventions.

Patients and Methods

Data normalization and quality control

We collected 21 publicly available datasets (described in the supplementary data) containing raw GE data from microarray analyses (Affymetrix© Gene Chip Human Genome HG-U133A and HG-U133Plus2) of 3,247 primary human breast cancer samples. The data were normalized by the robust multichip average (RMA) procedure from the EMA R package.⁶⁵ The datasets were split into training (HGU-133A Affymetrix© arrays, 12 datasets, n = 1,995) and validation (HGU-133Plus2 Affymetrix© arrays, (9 datasets, n = 1,014) sets. We also collected two large datasets, to validate our classification: The Ignatiadis dataset (n = 996) and the METABRIC dataset (n = 1,992). Data processing for these two datasets has been described elsewhere.^16,17

Determination and preprocessing of triple-negative breast cancer samples

We identified the TNBC samples in each dataset, using a bimodal mixture of two Gaussian distributions for ER and HER2 gene expression, and the median value for PR expression.

The training, validation and Ignatiadis datasets

Batch effects were eliminated by the median centering of each probe-set across arrays and by a, independent quantile normalization of all arrays for each dataset. We controlled for outliers with the Array Quality Metrics R package.

The METABRIC set

We fitted a linear model (limma R package) to remove the batch effect and probes were filtered according to three criteria: probe quality,⁶⁶ GC content and presence in more than 5% of the samples. We centered expression values, using the R function scale().

Gene selection process

Consensus clustering was applied to the training set, to determine the optimal number of robust gene clusters for the most variant genes (standard deviation > 0.8). We investigated the enrichment of each gene cluster in particular types of genes. We then identified known biological networks, for each gene cluster separately, using String© database software version 9.1 (http://string-db.org/).⁶⁷

We then applied a two-step selection process: (1) we selected strong biological networks by retaining only genes for which connection scores of at least 0.7 were obtained with String© database software, (2) within each biological network, we selected groups of genes with for which expression levels were correlated, with a correlation coefficient of at least 0.5.

For each dataset (the training, validation, Ignatiadis and METABRIC sets), we applied a hierarchical clustering procedure to the TNBC GE profiles, using the selected genes to visualize the optimal number of stable TNBC subtypes.

Prognostic analysis

Prognostic analysis was performed on the METABRIC set published by Curtis et al.¹⁷

Expression data were summarized by a metagene for each gene cluster (details in the supplementary material). The clinical and pathological variables available for each dataset are described in the supplementary data. Qualitative variables were compared in χ² tests or Fisher's exact tests, as appropriate. Quantitative variables were analyzed in Student's t-tests. Survival analyses were performed separately for patients with and without chemotherapy. Survival analyses were performed, with the Kaplan–Meier estimate of the survival function. The endpoint of these analyses was breast cancer-specific survival (BCSS). Survival curves were compared in log rank tests. Hazard ratios were estimated with Cox's proportional hazard model.

Expression of the gene signature in human triple-negative breast cancer cell lines

We downloaded the GE profiles of the human cancer cell lines from the Cancer Cell Line Encyclopedia (CCLE)⁶⁸ of Novartis/the Broad Institute and the Cancer Genome Project (CGP)⁶⁹ of the Sanger Institute. We normalized all the cell lines from different tissues together.

All statistical analyses were performed with R software (www.cran.r-project.org). P -values < 0.05 were considered statistically significant.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.

Acknowledgments

The authors thank Sergio Roman-Roman for the reviewing of the study and the manuscript.

Funding

This work was supported by Institut Curie, INCa (the french National Cancer Institute) Grant INCa-DGOS-4654, ANR-10-IDEX-0001-02 PSL, ANR-11-LABX-0043, CIC IGR Curie 1428 and Fundation ARC (association for Research Against Cancer).