Improved prognostic classification of breast cancer defined by antagonistic activation patterns of immune response pathway modules

Abstract

Background: Elucidating the activation pattern of molecular pathways across a given tumour type is a key challenge necessary for understanding the heterogeneity in clinical response and for developing novel more effective therapies. Gene expression signatures of molecular pathway activation derived from perturbation experiments in model systems as well as structural models of molecular interactions ("model signatures") constitute an important resource for estimating corresponding activation levels in tumours. However, relatively few strategies for estimating pathway activity from such model signatures exist and only few studies have used activation patterns of pathways to refine molecular classifications of cancer.

Methods: Here we propose a novel network-based method for estimating pathway activation in tumours from model signatures. We find that although the pathway networks inferred from cancer expression data are highly consistent with the prior information contained in the model signatures, that they also exhibit a highly modular structure and that estimation of pathway activity is dependent on this modular structure. We apply our methodology to a panel of 438 estrogen receptor negative (ER-) and 785 estrogen receptor positive (ER+) breast cancers to infer activation patterns of important cancer related molecular pathways.

Results: We show that in ER negative basal and HER2+ breast cancer, gene expression modules reflecting T-cell helper-1 (Th1) and T-cell helper-2 (Th2) mediated immune responses play antagonistic roles as major risk factors for distant metastasis. Using Boolean interaction Cox-regression models to identify non-linear pathway combinations associated with clinical outcome, we show that simultaneous high activation of Th1 and low activation of a TGF-beta pathway module defines a subtype of particularly good prognosis and that this classification provides a better prognostic model than those based on the individual pathways. In ER+ breast cancer, we find that simultaneous high MYC and RAS activity confers significantly worse prognosis than either high MYC or high RAS activity alone. We further validate these novel prognostic classifications in independent sets of 173 ER- and 567 ER+ breast cancers.

Conclusion: We have proposed a novel method for pathway activity estimation in tumours and have shown that pathway modules antagonize or synergize to delineate novel prognostic subtypes. Specifically, our results suggest that simultaneous modulation of T-helper differentiation and TGF-beta pathways may improve clinical outcome of hormone insensitive breast cancers over treatments that target only one of these pathways.

Figure 1

Measuring pathway module activation. Flowchart figure showing overall strategy used for inferring pathway module activity in clinical tumor samples from a model (perturbation) signature. A) A gene mRNA signature that represents a perturbed cancer cell phenotype (i.e oncogene overexpression) is combined with mRNA expression data of a large panel of clinical tumor specimens to derive an "expression relevance network" where nodes represent genes from the signature and an edge between two nodes indicates a statistically significant Pearson correlation between the two corresponding genes as measured over the clinical tumor panel. Having constructed the relevance network, the network is first pruned so that network edges that are inconsistent with prior information are removed. Signs on edges between labelled genes indicate the sign of the significant correlation between the two genes, which must be consistent with their directionality as given by the model signature. Modules defined as subnetworks with higher than average edge density are then inferred using a spectral decomposition algorithm (see Methods). B) For a given relatively large module, the module of pathway activation (MPA), pathway activity is then computed using a metric defined over the topology of the module. In the formula, PA_s stands for the estimated pathway module activity in sample s, M is the number of genes in the module, σ_i is a binary weight (1,-1) indicating the directionality of gene expression of gene i (1 = upregulated, -1 = downregulated), z_is is the z-score normalised gene expression value in sample s and A_ij is the adjacency matrix of the module. Effectively, this metric gives more weight to gene interactions that are supported by the data. Color and sign of nodes reflect the directionality of expression in the in-vitro signature (Red = upregulated &σ = 1, Green = downregulated &σ = -1). Pathway activity levels can then be shown as heatmaps (blue = high activity, yellow = low activity).

Figure 2

Clustering analysis over pathway modules. A) Heatmaps of pathway activation (blue = high relative activation, yellow = low relative activation) over the merged ER- and ER+ cohorts [2,24,40-44]. Color bars indicate the intrinsic subtype (Pink = HER2+, green = normal, dark-red = basal, skyblue = luminal A, blue = luminal B) and the cluster inferred using a variational Bayesian method [52]. B) Kaplan Meier plots for distant metastasis free survival (DMFS) for the predicted clusters in ER- and ER+ breast cancer, respectively.

Figure 3

Pathway module activation scores across intrinsic subtypes in ER+ breast cancer. A-H) For the selected module in each pathway, we show boxplots of predicted pathway module activation scores across the major intrinsic subtypes within ER+ breast cancer as estimated in Set1. I-P) Corresponding boxplots as estimated in Set2. Number of samples in each subtype shown above corresponding boxplot. (Color Code: green = normal, skyblue = lumA, blue = lumB, pink = HER2+).

Figure 4

Pathway module activation scores across intrinsic subtypes in ER- breast cancer. A-H) For the selected module in each pathway, we show boxplots of predicted pathway module activation scores across the major intrinsic subtypes within ER- breast cancer as estimated in Set1. I-P) Corresponding boxplots as estimated in Set2. Number of samples in each subtype shown above corresponding boxplot. (Color Code: green = normal, pink = HER2+, red = basal).

Figure 5

Correlation patterns of molecular pathway modules. A) Pearson correlation heatmaps between molecular pathway modules in the ER- and ER+ breast cancer (Set1), respectively (Red = high positive correlation, White = zero or insignificant correlation, Green = high negative correlation). B) Validation of pairwise pathway module Pearson correlations in external set (Set2). Left panel (ER-), right panel (ER+). Overall Pearson correlation (PC) between training (Set1) and validation set (Set2) is given.

Figure 6

Improved prognostic models through non-linear interactions of pathway modules. A) For pathways that correlated with DMFS in ER+ and ER- breast cancer (Table 2), we consider corresponding Boolean interaction Cox regression models describing the pairwise interaction of any two pathways (the best model out of a total of 4, i.e up-up, up-down, down-up, down-down, is shown). y-axis labels the pathway interaction pair and best boolean model, x-axis denotes the log-likelihood of the corresponding model. (Black = log-likelihood of model for first pathway in pair, Grey = log-likelihood of model for second pathway in pair, Red = log-likelihood of the best Boolean interaction model, pink dashed line highlights those Boolean models with improved log-likelihoods). B) Heatmaps of likelihood ratio test (LRT) p-values comparing nested prognostic models. Specifically, LRT p-value for pathway p_y on y-axis and pathway p_x on x-axis is obtained by comparing Cox-regression models with the single pathway p_x plus non-linear Boolean interaction B(p_x, p_y) as predictors against the model with only p_x as predictor. C) As B), but LRT p-value for pathway p_y on y-axis and pathway p_x on x-axis is obtained by comparing Cox-regression models with the single pathway p_x plus non-linear Boolean interaction B(p_x, p_y) as predictors against the model with only B(p_x, p_y) as predictor. Color codes: red (P < 0.01), pink (P < 0.05), white (P > 0.05). All Cox regression were stratified regression using the cohorts as strata to account for variations in the hazard rate between cohorts.

Figure 7

Novel prognostic subtypes in ER- and ER+ breast cancer. A) & B) Kaplan Meier DMFS curves for dichotomised pathway activity levels, for IL12 and TGFB, in ER- breast cancer (Set1). C) Corresponding Kaplan Meier curve for the four subtypes stratified according to up/down activity of the two pathways. Hazard ratio refers to the IL12up-TGFBdn subtype relative to the rest. D) Independent validation in the test cohort (Set2). E) & F) Kaplan Meier DMFS curves for dichotomised pathway activity levels, for MYC and RAS, in ER+ breast cancer (Set1). G) Corresponding Kaplan Meier curve for the four subtypes stratified according to up/down activity of the two pathways. Hazard ratio refers to the MYCup-RASup subtype relative to the rest. H) Independent validation in the test cohort (Set2).

Figure 8

Proposed model of how immune response pathways affect clinical outcome in ER- breast cancer. Figure adapted from [26]. Hypothetical model in which the balance of cytokines in the breast tumour microenvironment determines the relative strength of Th1 and Th2 differentiation. Stronger activation of a Th1 immune response leads to increased production of IL2 and IFNG which mediate formation of M1 macrophages and cytotoxic killer cells, which is tumour inhibitory [26]. Correspondingly we observe that genes that are upregulated in these pathways are associated with good prognosis (DMFS) in ER- breast cancer (significant associations shown in green). Conversely, stronger activation of a Th2 immune response leads to production of IL13 and TGFB cytokines through an M2 macrophage polarization program. The cytokine TGFB is known to suppress the tumour inhibitory role of Th1 [26]. Correspondingly, we observe that genes that are upregulated in these pathways confer poor prognosis (DMFS) (significant associations shown in blue). Genes implicated in the Th1 and Th2 pathways were generally anticorrelated, indicative of an unbalanced differentiation program. It follows from this model that simultaneous high Th1 (IL2, IL12, IFNG) and low TGFB would confer better prognosis than either high Th1 or low TGFB alone, in agreement with our observations.