Biclustering reveals breast cancer tumour subgroups with common clinical features and improves prediction of disease recurrence

Biclustering breast cancer gene expression data identifies tumour subgroups with similar clinical features

We analysed a previously published dataset which combines Affymetrix HG-U133A and HG-U133PLUS2 microarray data and matched clinical data from 6 cohorts of patients from independent studies. We selected 437 primary breast tumours from patients who had not received adjuvant or neo-adjuvant chemotherapy, endocrine-based therapy or radiotherapy (see Methods). We restricted our analysis to non-adjuvant treated tumours so that direct associations between molecular information and long-term disease outcomes might be observed without the confounding affects of different treatments. After processing, the dataset consisted of 7756 probe sets (genes) across the 437 tumours. The 7 clinical annotations for each tumour/patient were: lymph node (LN) status; hormone receptor status (estrogen receptor (ER) status and progesterone receptor (PgR) status); tumour grade; tumour subtype; tumour size; patient age. The clinical information associated with these 437 tumours is summarised in Figure 1. More than 90% of these were lymph node negative (LN-) tumours, which is consistent with the fact that a clinical decision had been made not to treat these patients with adjuvant therapy.

Using cMonkey to bicluster the dataset we obtained 44 biclusters with an average number of 20 tumours contributing to each bicluster. To determine whether the tumour subgroups and corresponding gene modules identified by biclustering were biologically and clinically meaningful, we assessed whether tumours allocated to the same bicluster had similar clinico-pathological annotations and clinical outcomes, and whether the genes contributing to these biclusters were associated with the activity of distinct biological pathways. To make sure that the tumour subgroups identified in the biclusters do not simply recapitulate clinical cohort we verified that biclusters contain data derived from a range of cohorts. Then we computed the proportion of tumours in each subgroup in each of the 7 clinical categories. Figure 2 illustrates the proportion of tumour samples in each category of (A) histological grade, (B) ER status, and (C) molecular tumour subtype in each of the biclusters, ranked in order of the proportion of patients who suffered early relapse (Disease Free Survival time ≤ 2 years). These results show that biclusters do not simply recapitulate any of the already available clinical annotations.

Analysis of the clinical annotations of tumours in each bicluster however reveals similar clinico-pathological features. Clinical annotations for each bicluster are summarised in the supplementary material (Additional file 1: Figure S1). For example, a summary of clinical annotations for biclusters 21 and 7 are shown in Figure 3. Both consist of 31 tumour samples collected from all six cohorts. In bicluster 21 most tumours (> 90%) have either Basal-like (42%) or Her2 (52%) intrinsic subtype; most (~87%) are estrogen-receptor negative (ER-), and are grade 2 (32%) or 3 (58%). Therefore, the clinical information for bicluster 21 suggests that tumours in this group are aggressive and have pathological features associated with poor prognosis. By comparison, tumours in bicluster 7 are predominantly Luminal-like tumours, i.e. Luminal A (64%) and Luminal B (35%), and are estrogen-receptor positive (~93%), and histological grade 1 (32%) or 2 (55%). Thus the clinical information of tumours in bicluster 7 suggests that patients in this group have clinical features that are generally associated with less aggressive tumours, associated with good prognosis. Similar patterns can be observed in the remaining biclusters (see Additional file 1).

Overall, these data suggest that biclustering identifies subgroups of tumours selected from multiple cohorts, which do not simply recapitulate existing clinical annotations, but which are similar in clinico-pathological characteristics.

Biclustered subgroups are associated with differential survival outcomes

To test the hypothesis that subgroups of tumours identified in biclustering are associated with differential patient outcomes, we examined the disease free survival (DFS) distribution of patients in each bicluster. Figure 4 shows Kaplan-Meier survival curves for all 44 biclusters (grey) superimposed on the survival curves for histological grade (panel A) and tumour subtype (panel B). Bicluster-associated survival curves show a greater range than the grade-associated and subtype-associated survival curves. To test whether this range of survival distributions might arise by chance, we randomly assigned (sampling with replacement) the 437 tumours into 44 arbitrary groups with the same distribution of sizes as found in the biclusters. The Kaplan-Meier survival curves of these arbitrary groups (see Additional file 1: Figure S2) visually suggest that the bicluster-specific survival curves are more widely separated than expected at random. In order to assess this statistically, we computed the logrank-p values for all possible pairings of biclusters and for all possible pairings of the randomised groups. A two-sample Kolmogorov-Smirnov (KS) test rejected at the 5% significance level the null hypothesis that these distributions were the same (p-value = 1.7193 × 10^-16). These results indicate that membership of the biclusters is associated with differential survival outcomes.

Tumours grouped in biclusters are molecularly distinct

To reveal molecular activities of oncogenic pathways associated with the biclusters, we used a univariate Cox PH model to determine the molecular features most strongly associated with clinical outcome for each bicluster, identifying probe sets whose expression was significantly associated with DFS (Cox PH p-value < 0.01). We then used the GeneSetDB web tool [33] to identify the significantly enriched (GeneSetDB p-value < 0.01) functional annotations of these genes (see Additional file 2). For example, RNAs associated with several of the biclusters were enriched for cell cycle annotations (e.g. cell division [GO: 0051301] in bicluster 1, cell cycle checkpoint [GO: 0000075] in bicluster 10, mitotic cell cycle [GO: 0000278] in bicluster 17 and mitotic cell cycle checkpoint [GO: 0007093] in bicluster 26), proliferation (e.g. decreased cell proliferation [MP:0000352] in biclusters 18, 20 and 26, cell proliferation [GO: 0008283] and decreased B cell proliferation [MP: 0005093] in bicluster 23, increased cell proliferation [MP:0000351] in bicluster 33). RNAs associated with several of the biclusters were also enriched for immune response annotations (e.g. immunological synapse [GO: 0001772] and wound healing [GO: 0042060] in bicluster 7). Genes strongly associated with specific biclusters were also found to be significantly enriched for cancer-associated transcription factors, such as MYC (e.g. biclusters 35 and 38), STAT1 (e.g. biclusters 7 and 28), E2F (e.g. V$E2F_Q3_01 and V$E2F1_Q6: E2F-1 in biclusters 25 and 35), SP1 (e.g. biclusters 30 and 33) and NFE2L2 (bicluster 23), and enriched in regulatory pathways, which play an important role in breast cancer progress, such as Regulatory DNA replication (REACTOME Pathway ID: 69304) and p53-Dependent G1/S DNA damage checkpoint (REACTOME pathway ID: 69620 and 69613)) in bicluster 23, G1/S-Specific Transcription (REACTOME pathway ID: 453279 and 69205) in biclusters 18 and 20, Interferon signalling (REACTOME pathway ID: 913531) in biclusters 7 and 10, and SREBP signalling (Wikipathway PW:0000753) in bicluster 23. Additional file 3: Table S1 summarises biological meaningful annotations for biclusters associated with low risk of disease recurrence and for biclusters associated with high risk of relapse.

To compare the gene sets identified for each bicluster we searched for biologically significant genes which were differentially expressed between pairs of biclusters, showing a significant difference (logrank p-value < 0.01; Additional file 1: Figure S3) in survival distribution (see Methods; Additional file 1: Figure S4). As an example, Additional file 3: Table S2 summarises significant differentially expressed genes between bicluster 7 and bicluster 21. We analysed these gene classifiers and found them to be highly enriched in GO categories including response to estrogen stimulus, mammary gland development, defence response and wound healing. (A selected list of biological terms significantly enriched with differentially expressed genes between biclusters 7 and 21 is summarised in Additional file 3: Table S3). The gene classifiers were also highly associated with transcription factor STAT group (TransFac, p-value < 0.001), and were significantly enriched in regulatory pathways such as Downregulated of MTA-3 in ER-negative Breast Tumors (Biocarta, p-value = 1.80 × 10^-4) and Role of ERBB2 in Signal Transduction and Oncology (Biocarta, p-value = 2.50 × 10^-3). Breast abscess (Disease/Phenotype, p-value = 5.90 × 10^-6) and abnormal mammary gland alveolus morphology (Disease/Phenotype, p-value = 7.90 × 10^-5) were also found in genes characterizing the difference in survival outcomes between these two biclusters.

Overall, these results imply that the biclusters are molecularly defined by functionally coherent, but distinct, gene expression modules.

Bicluster gene classifiers are different to commercial Oncotype DX and MammaPrint profiles

To determine whether currently known biomarkers are associated with the biclusters, in Figure 5 (A) we plot gene expression for tumours in each bicluster of genes from the commercialised Oncotype DX and MammaPrint diagnostic tests, the genetic grade signature proposed by Ivshina et al. [34], and other key oncogenic markers. The figure indicates three broad groupings of biclusters, which have significant differences in disease-free survival, shown in Figure 5 (B). The first 10 biclusters (cyan bar) are associated with low risk of disease recurrence (longest DFS). The last 10 biclusters (magenta bar) are associated with high risk of disease recurrence (shortest DFS). The heatmap suggests that these known biomarkers are not able to distinguish patient subgroups identified by biclustering.

We then looked at the genes and related pathways identified as being differentially expressed between biclusters associated with low and high survival prognosis. Interestingly, many of these gene classifiers are not found in currently available prognostic gene classifiers for breast cancer. For each pair of biclusters showing differentially survival outcomes, we used prediction analysis of microarrays (PAM) [35] to identify those of probe sets which were able to accurately characterise survival outcome differences between the biclusters (Additional file 4). Among the gene classifiers for bicluster 7 and bicluster 21 (see Additional file 3: Table S4), GRB7 and ERBB2 (highly expressed in bicluster 21) and ESR1 and BCL2 (highly expressed in bicluster 7) are members of the Oncotype DX gene list. Overall, however, we found only a small proportion of gene classifiers that coincide with the commercial gene lists.

Biclusters show strong prognostic association with breast cancer patient survival in comparison to conventional clinical markers

We used a univariate Cox PH model to examine the association between conventional prognostic factors (including hormone receptor status, lymph node status, tumour size and age) and disease free survival (DFS) time. Clinical variables were entered as binary data while tumour size and age were discretized into categorical variables (see Methods). The strength of association between clinical variables and the survival outcome, as quantified by Cox PH p-value, was compared to prognostic significance of biclustered tumour groups, using the same univariate Cox PH model. Bicluster membership was found to have a strong association with DFS (univariate Cox PH p = 8.67 × 10^-5), after LN status (p = 1.71 × 10^-5) and tumour size (p = 7.44 × 10^-5), shown in Additional file 3: Table S5 (a). We also found that biclusters give better prognostication of relapse compared to the conventional molecular tumour subtype (univariate Cox PH p = 0.000427).

Next, we used a multivariate Cox PH model to estimate the association of biclusters with DFS in combination with the conventional clinical factors. We found that biclusters were significantly statistically associated with DFS time (Additional file 3: Table S5 (a), multivariate Cox PH p = 0.00103; compared to randomly allocated groups, Additional file 3: Table S5 (f), multivariate Cox PH p = 0.6163).

Finally, we also found that biclusters show stronger prognostic association with disease recurrence within pathology-specific tumour subsets (Additional file 3: Table S5 (b)-(e): LN- patient subgroup, univariate Cox PH p = 0.000183; ER+ patient subgroup, p = 0.00406; Grade 2 tumour subgroup, p = 0.000171; Grade 3 tumour subgroup, p = 0.0153). In each of these pathology-specific tumour subsets, multivariate Cox PH analysis also demonstrated that the tumour biclusters were significantly in associated with DFS compared to molecular tumour subtype when combined with other conventional clinical variables (Additional file 3: Table S5 (b)-(e)).

These univariate and multivariate Cox PH analyses suggest that the classification of tumours using biclustering can, alone or in combination with other clinical factors, be a significant prognostic variable, especially in pathology-specific patient subsets, in comparison to conventional clinical-pathological factors.

Biclusters increase prediction of patient survival

The ultimate goal of tumour classification is to improve prediction of an individual’s risk of relapse, so as to aid the clinician in decision making on treatment strategy (e.g. whether a patient should or should not be given adjuvant therapy). To test whether biclusters improve prediction of patient survival, we used a conditional inference tree model to predict disease recurrence. We chose to examine disease recurrence (relapse) within 2 years of primary surgery. Therefore, we defined a dichotomous clinical outcomes decision tree, in which early relapse equates to disease free survival time ≤ 2 years, and late relapse to DFS time ≥ 8 years (defined in this manner the data set contains 92 early relapse and 211 late relapse patients). Bicluster membership (representing molecular information) and clinical information including hormone receptor status (ER status), lymph node status, histological grade, tumour size and patient age were included as binary variables for each tumour.

Initially we generated a decision tree based on conventional clinical covariates only. Figure 6 (A) demonstrates that lymph node (LN) status was found to be the most important clinical variable for estimation of patient clinical outcomes (p < 0.001) in comparison to the remaining clinical variables, followed by Tumour Grade (p < 0.001) and ER status (p = 0.088). Since ~90% of the tumours in our dataset are LN negative tumours, we constructed a decision tree for the LN− branch, shown in Figure 6 (B). In LN− tumours, histological grade was the most significant variable (p < 0.001), followed by tumour size (p = 0.027) and ER status (p = 0.183). Next we constructed a decision tree for LN− tumours using both biclusters and clinical variables. Figure 7 shows that the bicluster 36 (p = 0.001) and bicluster 21 (p < 0.001) were significant molecular covariates in LN− tumours compared to the other clinical variables. We determined similar decision trees for continuous DFS time without dichotomisation (data not shown).

To assess the performance of the conditional inference tree predictors, we performed a cross-validation by generating 50 random splits of the data into training (80%) and test (20%) sets, ensuring that both training and test sets contained similar proportion of early and late relapse patients as the full data set. Biclusters derived by applying cMonkey to the training data only were combined with clinical variables and used to generate the decision tree. Performance was assessed by comparing relapse time predicted by the decision tree to the actual survival outcomes of patients in the test set. Tumours in the test set were assigned to biclusters by calculating correlation with the first principal component of the gene expression profile of each bicluster (see Methods). To assess the contribution of biclusters to prediction of relapse times, we compared these results with relapse times predicted using decision trees constructed from clinical variables only. The average rate of correctly predicted early relapse increased from 41% to 58% by incorporating biclusters into the predictors.

To assess the performance of the predictors for a clinically-determined subset of tumours we estimated of early relapse for patients with moderately differentiated grade 2 and poorly differentiated grade 3 tumours. We found 17% and 20% improvements in prediction of early relapse in these patients by including bicluster information (increase from 39% to 56% for patients with grade 2 tumours and from 46% to 66% for patients with grade 3 tumours).

We note that the conditional inference tree prediction model performed consistently with different choice of survival time boundaries (data not shown): no significant deviation in the performance was observed when different cut-off times for relapse were chosen. We also considered an alternative prediction model. Bøvelstad et al. [36, 37] performed a comparative study of six Cox regression-based survival prediction methods, including univariate selection, supervised principal component analysis, partial least squares regression, ridge regression and lasso, finding that ridge regression presented the best overall performance. We repeated the above analysis using a Cox regression-based ridge regression model. A higher median area under the ROC curve (AUC) [38] was observed, in particular for DFS time less than 2 years, across 50 random data splits when the model incorporated bicluster information (Figure 8), indicating better prediction of early relapse compared to the model with clinical covariates only, qualitatively consistent with the results from the conditional inference tree.

Taken together, these results suggest that combining conventional clinical information with bicluster membership is able to provide a better prediction of early relapse in breast cancer patients.

Non-adjuvant treated breast cancer dataset

This dataset contains the expression profiles of more than 20,000 probe sets across 960 patient samples, measured using Affymetrix HG-U133A and HG-U133PLUS2 GeneChip arrays, quantile normalized across all patient samples, as described in Lasham et al. [41]. This study described gene expression microarray data and matched clinical data for primary breast tumour patients assembled from multiple published studies (GSE3494, GSE1456, GSE6532, GSE7390, GSE4922, GSE36771, and GSE36772).

In order that direct associations between molecular information and long-term disease outcomes might be observed without the confounding affects of different treatments, we restricted our analysis to patients who had not received adjuvant or neo-adjuvant chemotherapy, endocrine-based therapy or radiotherapy, resulting in a dataset for analysis containing 437 tumours from the following cohorts (Figure 1): GSE1456 (33 samples), GSE3494 (122), GSE6532_U133A (44), GSE7390 (198), NZ_U133A (30 samples; part of GSE36771) and Singapore (10 samples; 6 from GSE4922 and 4 from GSE36772).

We removed all genes without significant variation across samples by estimating a noise threshold based on the mean signal of the BioB 3’ probe set across all chips, and retaining all probes that were expressed above this noise level in at least 10% of chips. Since Affymetrix HG-U133A microarrays include multiple probe sets annotated with same official gene symbol (OGS), to avoid inclusion of outlying probe sets and to reduce computational cost, probe sets without OGS annotation in version 2.6.3 of the Bioconductor hgu133a.db annotation package were also removed. We retained 7756 probe sets for analysis – one probe set each OGS (whose expression had the highest correlation to the median expression of all probe sets with the same annotation).

Finally, to remove non-biological variation associated with cohort, ComBat [42] processing was then applied to adjust for cohort-correlated batch effect across the non-adjuvant treated tumour data set.

Biclustering using cMonkey to identify molecular-based tumour subgroups

We used the cMonkey biclustering algorithm [32] to identify subgroups of non-adjuvant treated tumours with coherent patterns of gene expression over a subset of genes. cMonkey performs 2D clustering on both genes and conditions (microarrays) by using a Markov chain process that iteratively adds or removes genes and conditions to or from a cluster, initially with a single gene, until the optimal cluster is achieved [43]. Using this approach, genes and tumour samples can be allocated into more than one bicluster, consistent with the principle that genes may participate in more than one biological process or gene expression module [44].

We note that we used cMonkey to identify clusters of tumours conditional on common expression profiles across a set of genes, rather than the standard usage (to identify co-expressed genes across a subset of conditions), and therefore the data matrix was presented in a transposed form in comparison to normal usage. Since the genome information for Homo sapiens was not available for the new R package version of cMonkey (version 4.8.4) [45], the motif and network information was switched off in the input of the cMonkey command line.

We chose to use cMonkey because of its computational efficiency in dealing with high-dimensional microarray data, and its ability to detect co-expression patterns by considering the membership of tumours and genes simultaneously without large overlap between biclusters. However, a limitation is the relatively narrow distribution in the number of genes used to determine each bicluster – a known current limitation of cMonkey. While this biases the analysis towards biclusters defined by gene sets of similar size, the gene membership was shown to be different between biclusters.

Differences in survival distribution – logrank test

We used the Kaplan-Meier plot to estimate disease free survival (DFS) distributions. The logrank test was used to assess the differences in survival distribution between any two sets of survival curves (smaller p-value indicates stronger statistical evidence against the null hypothesis that two curves were drawn from the same survival distribution). The Kaplan-Meier estimation and logrank test were performed using MatLab scripts kmplot and logrank[46].

PAM to determine molecular classifiers for each bicluster

Prediction analysis of microarrays (PAM) [35] was used to determine the minimal set of genes classifying each pair of biclusters associated with low and high risk of recurrence, with the lowest misclassification rate. Initially, we ran the PAM algorithm with all 7756 probe sets as input in the low risk group to high risk group comparison and obtained a minimal set of probe sets, which gave the lowest misclassification error. All probes in the minimal set are differentially expressed between the low and high risk-associated biclusters. This was executed using the pamr package in R.

Cox proportional hazards regression

The prognostic effect of molecular and clinical-pathological covariates on patient survival was assessed within the framework of Cox proportional hazard (PH) regression. A univariate Cox PH model was used to examine the individual association between each covariate and survival time. A multivariate Cox PH model was used to test the prognostic effect of combinations of covariates. The significance of association was quantified by the Cox PH p-value: smaller p-values indicate the stronger contributions of molecular and clinical variables to survival outcomes. The computation of the association between gene signatures and biclusters was implemented using the coxphfit function in MatLab and prognostic association was assessed using the survival package in R.

Cross validation: Association of test data with biclusters

In order to determine the association between a tumour from the ‘test’ data set and the biclusters calculated from the ‘training’ set, we used the first principal component (PC1) as a representative gene expression profile for each tumour sample in each bicluster. We computed the Pearson’s correlation between these profiles and the gene expression profiles of tumours in the ‘test’ set for each bicluster. ‘Test’ tumours were associated with biclusters for which the correlation coefficient was larger than a pre-defined threshold value (> 0.9), such that each ‘test’ tumour sample could be allocated into more than one bicluster.

Survival prediction model

To detect whether classification of tumours using biclustering can be used to improve the prediction of survival with censoring, we used two different approaches to estimate a regression relationship between survival time and the status of multiple covariates including both clinical and molecular variables (biclusters): (1) a conditional inference tree [47] constructed by recursive binary partitioning process and (2) ridge regression method based on Cox regression model [48]. For these models the covariates were entered as categorical variables, i.e. ER status (ER+, ER-); PgR status (PgR+, PgR-); LN status (LN+, LN-); Histological Tumour Grade (G1, G2, G3); Tumour Size (size < 10 mm, 10 mm ≤ size 20 mm, 20 mm ≤ size <40 mm, and size ≥40 mm); Patient Age (age < 40 years, 40 years ≤ age < 50 years, 50 years ≤ age < 60 years, and age ≥ 60 years); Bicluster Membership (0,1). The conditional inference tree was implemented through ctree function in the R package party. MatLab code for the Cox regression based ridge regression method is provided by Bøvelstad et al. [36].

The computational methods used in this study will be made available on request to the authors.