GOBO: Gene Expression-Based Outcome for Breast Cancer Online

Microarray-based gene expression analysis holds promise of improving prognostication and treatment decisions for breast cancer patients. However, the heterogeneity of breast cancer emphasizes the need for validation of prognostic gene signatures in larger sample sets stratified into relevant subgroups. Here, we describe a multifunctional user-friendly online tool, GOBO (http://co.bmc.lu.se/gobo), allowing a range of different analyses to be performed in an 1881-sample breast tumor data set, and a 51-sample breast cancer cell line set, both generated on Affymetrix U133A microarrays. GOBO supports a wide range of applications including: 1) rapid assessment of gene expression levels in subgroups of breast tumors and cell lines, 2) identification of co-expressed genes for creation of potential metagenes, 3) association with outcome for gene expression levels of single genes, sets of genes, or gene signatures in multiple subgroups of the 1881-sample breast cancer data set. The design and implementation of GOBO facilitate easy incorporation of additional query functions and applications, as well as additional data sets irrespective of tumor type and array platform.

section_content: Breast cancer is one of the most common malignancies in the world. It is evident that breast cancer is a heterogeneous disease characterized by different tumor biological and histological subtypes, age of onset, clinical course and responsiveness to treatment. Although currently used prognostic and treatmentpredictive markers for breast cancer have markedly improved treatment decisions, additional markers are required to further tailor treatment for individual patients. 

Microarray-based gene expression profiling has been used extensively to characterize the transcriptome of breast cancer, resulting in the identification of new molecular subtypes and markers or signatures of potential therapeutic and prognostic importance [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]. Although microarray-based gene expression profiling of breast carcinoma holds promise for future improvement in treatment decisions, current studies aimed at deriving prognostic or treatment predictive gene signatures are often limited by a low number of included patients, typically tens to hundreds, due to cost and availability of suitable tissue specimens. The growing accumulation of data from prognostic and/or treatment predictive gene expression studies on breast cancer in public repositories, e. g. Gene Expression Omnibus [13], provides an opportunity to construct pooled gene expression data sets comprising a larger number of patients with sufficient clinical data. Pooled gene expression data sets of different sizes have been used for prognostic validation of single and multigene predictors [14] [15] [16] [17], identification of gene set modules associated with prognosis [7, [18] [19] [20], distinction of clinical types of breast cancer [21], and dissection of underlying processes in reported prognostic predictors [19, 22]. Although data sets in larger pooled sets are often individually available through public repositories, the combined pooled data sets are commonly not available, thereby representing an obstacle to researchers with limited bioinformatic resources. To address this latter shortcoming, we developed Gene expression-based Outcome for Breast cancer Online (GOBO). GOBO is aimed at providing an online tool for prognostic validation of single genes, sets of genes or simple predictors in a pooled breast cancer data set comprising 1881 cases from eleven public data sets analyzed using Affymetrix U133A arrays. In addition, GOBO allows online investigation of the expression of single genes or sets of genes across a large set of commonly used breast cancer cell lines [23], as well as identification of co-expressed genes in both breast cancer cell lines and breast tumors. Using the CCNB1 gene and a Core Serum Response gene signature as two examples we illustrate the usefulness of GOBO for rapid online analysis in a large breast cancer data set
section_name: Introduction

section_name: Results

section_content: The web interface of GOBO allows precompiled data sets to be queried by the three main applications of GOBO: Gene Set Analysis (GSA), Co-expressed Genes (CG), and Sample Prediction (SP) (Figure 1 ). Currently, the precompiled data sets consist of gene expression data and annotation data for a pooled 1881-sample breast tumor set and 51 previously reported breast cancer cell lines [23]. The 1881-sample breast tumor set comprises 11 public data sets (Table 1 ) analyzed using Affymetrix U133A arrays and processed as described (in [15] and File S1). GSA is further divided into outcome analysis in breast tumors (GSA-Tumor) and expression patterns in breast cancer cell lines (GSA-Cell line). In both GSA applications the input is either a single gene or probe identifier, or a set of gene/probe identifiers (referred to as a gene set hereinafter). CG allows identification of co-expressed genes by provision of a single gene identifier in both the breast tumor data set and the panel of breast cancer cell lines. SP allows users to investigate the association of their classifiers (in certain predefined forms) with outcome in the 1881-sample breast cancer set. 

Pooling of gene expression data sets generated at different institutions and time points may introduce artifacts due to, e. g., different reagent lots, scanner settings, RNA extraction procedures, sample selections and laboratory practices. To investigate whether such data set related artifacts existed in the 1881-sample set we undertook three investigations. First, we performed principal component analysis (PCA) and analyzed associations between principal components and a variety of sample annotations (Figure S1 in File S2). We observed that the different data sets are contributing variation, however, in combination with biological/ clinical features. Importantly, the first three components are significantly associated with molecular subtypes, ER status, histological grade and lymph node status. Since some of the included data sets contain only ER-negative, ER-positive or node negative tumors, respectively, variation associated with data sets cannot be corrected for without the risk of deleteriously affecting biological/clinical features. Second, we examined gene expression levels of Affymetrix probes suggested from literature to provide high sensitivity and specificity in assigning cases to clinical subtypes based on ER and HER2-status. Expression levels were investigated in nine of eleven data sets for samples stratified by ER-status and PAM50 gene expression subtype classes (Figures S2 to S4 in File S2). These analyses show consistency across data sets in the expression of ER and HER2 across samples clinically annotated with ER-status, as well as in what would be expected for the PAM50 gene expression subtypes. Finally, we performed hierarchical clustering of the 1881sample data set using the PAM50 gene set [24], showing that the different data sets are intermixed (Figure S5 in File S2). Taken together, these validations indicate that the 1881-sample set can be used for the types of supervised analyses made available through GOBO. Moreover in the GOBO applications, we provide analysis that reveals the influence of each data set on the result.
section_name: Overview of GOBO applications

section_content: Numerous studies have investigated the prognostic importance of mRNA expression, as measured by microarrays or quantitative PCR, from single genes or sets of genes in different breast cancer subgroups [14, 18, 19, 25]. An online tool was recently reported that allows Kaplan-Meier survival analysis of single genes in a pooled 1809-sample microarray-based breast cancer data set [26]. The GSA-Tumor application of GOBO extends on the application reported by Gyorffy et al. [26] by allowing: 1) expression of gene sets to define patient cohorts based on a larger number of quantiles, 2) gene weights to be used in computation of a combined expression score, 3) time-censoring, 4) multivariate analysis of gene expression groups with other covariates such as estrogen receptor (ER) status, lymph node (LN) status, histological grade, stratified tumor size (#20 mm or. 20 mm), stratified patient age (#50 years or. 50 years), 5) analysis of gene set expression across clinical (ER, LN, and histological grade) and molecular breast cancer subgroups [1, 24], and 6) correlation of expression of individual genes in a gene set with co-expressed gene modules emulating breast cancer specific biological processes (Fredlund et al., Manuscript in preparation). 

In GSA-Tumor a gene set comprises either a single or multiple genes. Weights can be associated with each gene, thus allowing genes with positive and negative relative expression levels to be combined into a gene set. In the case of multiple genes an averaged gene set expression (including weights if supplied) is computed, prior to dividing the data set into patient cohorts based on gene expression quantiles. Association with outcome is investigated for stratified patient cohorts using distant metastasis-free survival (DMFS), relapse-free survival (RFS) or overall survival (OS) in different breast cancer subgroups. GSA-Tumor allows division of the 1881-sample data set into five major subgroups based on clinical annotations available from GEO [13] ; 1) all 1881 tumors, 2) ERpositive tumors (n = 1225), 3) ER-negative tumors (n = 395), 4) systemically untreated patients (n = 927), and 5) patients treated with tamoxifen alone (n = 326, GSE6532 [27] and GSE12093 [28] ). In each of the selected major subgroups, survival analysis including Kaplan-Meier, univariate and multivariate analysis, is performed on subsequently more stratified groups based on gene expression subtypes (basal-like, HER2-enriched, luminal A, luminal B, normal- like) [24, 29], ER-status, LN-status, histological grade, and treatment status (untreated or tamoxifen treated tumors). 

Proliferation has been recognized as one of the key prognostic factors in breast cancer, and has been found to be the major constituent of several prognostic gene expression signatures [19, 22]. Aberrant cell proliferation requires deregulation of several cellular functions or pathways, including cell cycle progression. Cyclins and cyclin-dependent kinases (CDKs) are key regulators of the cell cycle, controlling the progression and initiation of different phases of the cell cycle [30]. One of the key components of the transition between the G2 and M phase in the cell cycle is the protein complex formed by cyclin B1 (CCNB1) and CDK1 [31]. CCNB1 protein levels have been shown to differ between breast cancer subgroups [32], tumors of different histological grade [32, 33], and to be associated with breast cancer outcome [32] [33] [34]. In addition, CCNB1 is also included in several prognostic gene signatures such as the 21-gene recurrence score [10] and two genomic grade signatures [35, 36]. Using GSA-Tumor we stratified the 1881-tumor set into three quantiles based on CCNB1 gene expression (File S1). GSA-Tumor analysis of CCNB1 gene expression shows, in agreement with previous protein-based studies, that elevated expression is associated with the basal-like and luminal B gene expression subgroups (Figure 2A ), higher histological grade (Figure 2B ), and worse clinical outcome predominantly in ER-positive breast cancers (Figures 2C and D ) supported by multivariate analysis (Figure 2E ) [32] [33] [34]. Gene set expression in breast cancer cell lines (GSA-Cell line)

Cancer cell lines, derived from immortalized cancer cells, are commonly employed in cancer research as model systems for functional assessment of aberrant gene function and prediction or indication of response to therapeutic substances. Furthermore, it has been shown that breast cancer cell lines in many aspects mirror the heterogeneity in transcriptional and genomic aberrations found in breast tumors [23, 37], underlining the importance of selecting the appropriate cell line for a specific experiment. The GSA-Cell line application provides an interface for exploring mRNA expression levels of gene sets across 51 previously reported breast cancer cell lines [23], reporting both relative gene expression levels across the cell line panel, as well as summarized average intensity values. Similarly as for GSA-Tumor, if a gene set consists of multiple genes an average expression is computed for the total gene set, taking consideration to gene weights if supplied. 

Using the GSA-Cell line application, we investigated CCNB1 mRNA expression levels across the breast cancer cell line panel. Highest gene expression was observed in the basal B subgroup associated with a more stem-cell like phenotype and recently also the claudin-low subtype [23, 38], followed by the basal A and luminal subgroups as defined by Neve et al. [23] (Figures 3A and B ). In line with previous reports analyzing CCNB1 protein expression we found that triple-negative cell lines displayed highest CCNB1 mRNA expression levels (Figure 3C ) [32].
section_name: Association of gene sets with breast cancer outcome (GSA-Tumor)

section_content: Utilization of gene expression microarrays has made it clear that for a given gene investigated in a large enough data set other genes exist that display high similarity in expression pattern across samples. This co-expression of multiple genes is anticipated, given that functional processes and molecular pathways in an organism involve a highly organized interplay between multiple genes. Moreover, the co-expression of genes suggests that genes in a gene signature are often interchangeable as illustrated in a reanalysis of the van't Veer et al. [12] data set [39], indicating that it is the functional processes captured by a gene signature and not the individual genes that is of importance [22]. Recently, the approach of constructing metagenes or gene expression modules aimed at mimicking key biological processes in breast cancer has been utilized to identify gene sets with prognostic or treatment predictive power in breast cancer subgroups [4, 18, 19, 40]. These metagenes or gene expression modules are often based on identification of co-expressed genes to a specific prototype gene believed to represent the biological process of interest. 

The CG application of GOBO allows single genes to be used as prototype genes for identification of co-expressed genes by Pearson or Spearman correlation analysis in both subsets of the breast tumor data set (similar as GSA-Tumor) and the breast cancer cell line data set. The list of co-expressed genes can be further submitted to iterative correlation analysis, identifying which of the co-expressed genes that show tight co-expression with each other in an iterative process across the selected data set. The result of such an analysis could be a smaller cluster of genes that all are highly correlated with each other, representing a more distinct metagene or gene module (Fredlund et al., Manuscript in preparation). Using the CG application, we searched for genes highly co-expressed with CCNB1 in a 1751-sample subset of the 1881-sample breast cancer data set (excluding Chin et al. [41] cases) (File S1). Using a Pearson correlation cut-off. 0. 6 combined with a gene expression log2 standard deviation cut-off. 0. 9, we identified 34 genes to be positively co-expressed with CCNB1 (Table S1 ). When subjected to iterative correlation analysis one by one, all 35 genes showed at least 5 connections (Pearson correlation. 0. 6) with other highly co-expressed genes (Figure 4A, visualized using Cytoscape [42] ). Functional analysis of the 35 genes using DAVID [43, 44] confirmed the strong association with cell cycle processes, specifically M-phase and mitosis (Table S2 ). CCNB1 mRNA expression has been reported to be significantly and positively correlated with mRNA expression levels of CENPE, AURKB, PLK1, and PLK4 in both breast tumors and breast cancer cell lines [32]. None of these genes were present in the list of 34 genes from the CG analysis, explained by lower correlations (CENPE: 0. 47, AURKB: 0. 42, PLK1: 0. 39, PLK4 0. 38) to CCNB1 in our analysis. However, when performing a CG analysis in the Neve et al. cell line panel we found PLK1, AURKB and CENPE to be correlated (Pearson correlation. 0. 5) with CCNB1 given a gene expression log2 standard deviation cut-off. 0. 4. Association with outcome for a derived co-expressed gene set may be further investigated by subsequent GSA-Tumor or SP analysis (Figure 4B, File S1), while association with biological processes could be made through usage of public functional annotation tools.
section_name: Co-expressed Genes

section_content: A multitude of microarray-based gene signatures have been derived stratifying breast cancer into subtypes [1, 2, 24], identifying patients at risk [5, 6, 8, 9, 11, 12], and predicting response to therapeutic agents [3, 4]. However, it has become evident that gene signatures often represent the same biological processes [22], are associated with outcome only in subgroups of breast cancer [18, 19], and that validation of gene signatures needs to be performed in large independent data sets [45]. The SP application of GOBO is aimed at providing an interface for prognostic validation of simple gene signatures in similar subsets of the 1881sample breast cancer data set as for GSA-Tumor. Currently, SP supports three forms of predictor designs; 1) PAM-clustering [46], 2) stratification by expression quantiles, and 3) correlative centroid prediction. SP analysis includes, similar to GSA-Tumor, Kaplan-Meier analysis, univariate and multivariate analysis of predicted groups, as well as correlation of predictor genes to functional modules (Figure 1B ). 

Cancer invasion and metastasis represent the often-fatal steps in tumor progression, but remain poorly understood. Clearly, modifications in the tumor microenvironment involving not only tumor cells but also surrounding tissue, e. g., tumor-associated fibroblasts, are required for a tumor to become invasive and able to metastasize. Observations of histological similarities between the tumor microenvironment and normal wound healing have led to the hypothesis that tumors represent ''wounds that do not heal'' [47]. Fibroblasts represent a key component in wound healing, and in an attempt to identify a canonical gene expression signature of the fibroblast serum response to wound healing Chang et al. [48] reported a Core Serum Response (CSR) signature that was subsequently found to be associated with outcome in breast cancer. Identification of CCNB1 co-expressed genes and their association with outcome as a gene set. (A) GOBO analysis using the Coexpressed genes application identified 34 genes to be highly correlated with CCNB1 in a 1751-sample subset of the full combined breast cancer data set (excluding Chin et al. cases). Iterative correlation analysis of the 35 genes showed that all genes were highly correlated to each other with at least 5 connections, as visualized by a Cytoscape V2. 6. 3 spring embedded network. Each connection is visualized as a line between two genes. CCNB1 is highlighted in yellow. (B) PAM clustering into two groups of samples using the Sample Predictor (SP) application with the 35 genes in (a), followed by Kaplan-Meier analysis for 21 subgroups of the 1881-tumor set using DMFS as endpoint and 10-year censoring for 1379 cases with DMFS follow-up. Logrank P-values are shown as 2log10(P-value). doi:10. 1371/journal. pone. 0017911. g004 patients [8]. In order to investigate the prognostic association of the CSR signature in subgroups of breast cancer we submitted the CSR activated fibroblast centroid [8] to SP analysis (File S1). Analysis of the association for the CSR signature with different clinical and molecular variables showed that the activated fibroblast class of tumors was predominantly associated with typical high-proliferative poor outcome classes in breast cancer, such as ER-negative tumors, histological grade 3 tumors, basal-like and luminal B classified tumors (Figure 5A and B ). In contrast, the non-activated fibroblast class of tumors was associated with ERpositive tumors, histological grade 1 and 2 tumors, luminal A and normal-like classified tumors (Figure 5A and B ). As expected from these results the CSR signature was associated with outcome for all tumors (Figure 5C ). SP analysis also reveals if a signature is associated with outcome within subgroups of breast cancer stratified using clinical and molecular variables. Association with outcome for the CSR signature was mainly observed for ER-positive, LN-positive, LNnegative, untreated tumors, tamoxifen (TAM) alone treated tumors, or ER-positive/LN-negative tumors (Figure 5C ). In further support, associations with outcome for CSR in ER-positive (Figure 5D ), LN-positive and TAM tumors were supported by multivariate analysis using ER-status, LN-status and stratified histological grade as covariates when applicable. Taken together, the CSR SP and CCNB1 GSA-Tumor analyses are consistent with results from meta-analyses of different prognostic gene signatures, including CSR, identifying proliferation as the major constituent of several signatures [19, 22]. These meta-analyses have shown the signatures to be predominantly associated with outcome in ERpositive breast cancers, whereas the signatures identify the majority of ER-negative breast cancers as poor outcome and cannot be used to further stratify ER-negative tumors into groups with different outcome. These findings reinforce that prognostic gene signatures in breast cancer need to be evaluated in a subtypespecific manner in sufficiently large sample sets [18, 19, 22, 49]. Interestingly, a majority of genes in the CSR centroid did not show clear co-expression with any of eight different co-expressed gene expression modules, including two proliferation modules (Checkpoint and M-phase) (Fredlund et al. Manuscript in preparation) (Figure 6A ). Moreover, in exploratory analysis we found that genes in the CSR centroid did not show extensive co-expression when tested pair-wise in the 1751-sample subset used by the CG application (Figure 6B ). Similarly, when dissecting several reported gene signatures Wirapati et al. found that when using only a subset of signature genes associated with proliferation, performance was identical or even improved for certain signatures [22]. This observation suggests that many prognostic gene signatures contain only a smaller fraction of genes actually associated with outcome, and a significant proportion of passenger genes that may perturb or even have a detrimental effect on the prognostic performance.
section_name: Sample Prediction

section_content: GOBO is a convenient and user-friendly online tool for preliminary analysis of association with outcome for gene expression levels of single genes, sets of genes or gene signatures in a large public breast cancer microarray data set. Moreover, GOBO offers the possibility of investigation of gene expression levels in breast cancer subgroups and breast cancer cell lines for gene sets, as well as creation of potential metagenes based on iterative correlation analysis to a prototype gene. We have here demonstrated the usefulness of GOBO using CCNB1 and the CSR gene signature as two examples for rapid online extraction of gene expression patterns, co-expressed genes, and survival analysis in a large breast cancer data set. 

followed by Kaplan-Meier survival analysis in 21 subgroups using 1379 cases with DMFS follow-up. Logrank P-values are shown as 2log10(P-value). (D) Corresponding multivariate analysis for ER-positive tumors (n = 554) using lymph node status and stratified histological grade (histological grade 1 and 2 vs. 3) as covariates, DMFS as endpoint and 10-year censoring. ER status is omitted from the multivariate analysis since all investigated cases are ER-positive. doi:10. 1371/journal. pone. 0017911. g005
section_name: Conclusions

section_content: GOBO is designed in three modules; a data set module, a web interface module, and a data processing module (Figure 1 ). The data set module contains normalized gene expression data sets, and corresponding annotation data compiled into R-data objects. The GOBO web-interface module is written using Perl CGI. In the data processing module, application specific Perl scripts are used to parse parameters passed from web-interface, and to launch application specific R-scripts for analysis. R-scripts perform all analyses using standard R packages, such as survival for outcome analysis. Upon completion, R-scripts return analysis results as files available from the web-server for download. In order to simplify extensive analysis comparisons, storage of analysis results and to provide plots as scalable vector graphics supporting incorporation in publications output from GOBO analyses are reported as text and pdf files.
section_name: Design and Implementation

section_content: The present study was supported by the Swedish Cancer Society, the Foundation for Strategic Research through the Lund Strategic Centre for Translational Cancer Research (CREATE Health), and the Swedish Research Council. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

section_name: Availability and Future Directions

section_content: GOBO is freely accessible from the GOBO website http://co. bmc. lu. se/gobo. The design and implementation of GOBO facilitate easy incorporation of additional query functions and applications, as well as additional data sets irrespective of tumor type and array platform in the form of precompiled R-data sets. Future integration of data sets with additional information such as mutation status, copy number alterations or methylation data may further enhance the usefulness of GOBO.
section_name: Availability and Future Directions

section_name: Supporting Information