Meta-analysis of prostate cancer gene expression data identifies a novel discriminatory signature enriched for glycosylating enzymes

Background Tumorigenesis is characterised by changes in transcriptional control. Extensive transcript expression data have been acquired over the last decade and used to classify prostate cancers. Prostate cancer is, however, a heterogeneous multifocal cancer and this poses challenges in identifying robust transcript biomarkers. Methods In this study, we have undertaken a meta-analysis of publicly available transcriptomic data spanning datasets and technologies from the last decade and encompassing laser capture microdissected and macrodissected sample sets. Results We identified a 33 gene signature that can discriminate between benign tissue controls and localised prostate cancers irrespective of detection platform or dissection status. These genes were significantly overexpressed in localised prostate cancer versus benign tissue in at least three datasets within the Oncomine Compendium of Expression Array Data. In addition, they were also overexpressed in a recent exon-array dataset as well a prostate cancer RNA-seq dataset generated as part of the The Cancer Genomics Atlas (TCGA) initiative. Biologically, glycosylation was the single enriched process associated with this 33 gene signature, encompassing four glycosylating enzymes. We went on to evaluate the performance of this signature against three individual markers of prostate cancer, v-ets avian erythroblastosis virus E26 oncogene homolog (ERG) expression, prostate specific antigen (PSA) expression and androgen receptor (AR) expression in an additional independent dataset. Our signature had greater discriminatory power than these markers both for localised cancer and metastatic disease relative to benign tissue, or in the case of metastasis, also localised prostate cancer. Conclusion In conclusion, robust transcript biomarkers are present within datasets assembled over many years and cohorts and our study provides both examples and a strategy for refining and comparing datasets to obtain additional markers as more data are generated. Electronic supplementary material The online version of this article (doi:10.1186/s12920-014-0074-9) contains supplementary material, which is available to authorized users.


Background
Alterations in transcriptional programmes are often involved in neoplastic transformation and progression and defining these changes will help to understand the underlying biology of the malignancies. Gene Expression Microarray Analysis and more recently high-throughput RNA sequencing (RNA-seq) are commonly used techniques when trying to acquire an unbiased view of the expression levels of large numbers of genes. In order to define more compact and manageable expression modules that might predict risk or prognosis, various approaches have been used across several studies. These include the identification of consensus profiles across multiple datasets [1] and identifying biologically categorised gene sets with enriched representation of deregulated genes [2,3]. Furthermore, smaller expression modules have also been identified using hierarchical clustering methods to generate clusters containing genes with similar expression profiles across glioblastoma samples [3]. The high degree of prostate tissue heterogeneity, however, represents a challenge for transcriptomics since the relative prevalence of each cell type within a given sample determines the overall expression profile. This makes it difficult to compare prostate samples that have very different epithelial and stromal contents. Many studies have compared tumor tissue with benign hyperplastic tissue, or with nontumoral prostate tissues that were not precisely characterised in terms of location or epithelial representation. Therefore, the outcomes of these analyses were possibly biased because the comparisons included tissues of diverse histological or embryological origins. Various approaches have been used to overcome this issue including in silico corrections to compensate for variable epithelial representations in different samples [4], and laser microdissection combined with in vitro linear amplification [5]. The laser capture microdissection study of Tomlins et al. yielded several informative molecular concepts (multi-gene modules), which provide a rich source of data for further refinement and follow-up as well as distinguishing between stromal and epithelial cancer signatures [5]. It is, however, not clear how detectable those concepts might be in material extracted from heterogeneous whole tissue sections, an important point given the time and expense associated with laser capture microdissection.
In this study, we have therefore set out with a number of goals. First and foremost amongst these was to determine whether we could identify gene signatures that were statistically significant in datasets generated from both whole tissue sections and laser capture microdissected material. If so, this might indicate that with the right filtering approach, sample heterogeneity might not be a completely confounding challenge to transcriptomic analysis. Secondly, if we were able to identify such signatures, we then wanted to be able to refine them to a point that the signature and any pathway or process enriched within it could be easily validated by other experimental and clinical research groups. Here, we report a concise 33-gene signature with biological enrichment for glycosylation, which discriminates between benign tissue and prostate cancer (PCa) across multiple transcript detection platforms and sample types.

Description of datasets
Five datasets were downloaded and used in this study.

Prostate cancer dependent expression changes
To generate an initial broad progression dependent gene set, we used the prostate progression expression dataset GSE3325 (NCBI GEO database). We quantile normalised probe level intensity values and generated probe set signal estimates using RMA (1,2). We first characterised reporters with a coefficient of variance of less than or equal to 0.05 as uninformative and removed them from further analysis. Reporters having intensities below the 10 th quantile (3.91) in more than 75% of the samples were also removed. We identified progression associated expression changes by linear model. Primary tumour versus benign and metastatic versus primary contrasts were run and differential reporters identified using a 0.01 FDR threshold. Reporters were further filtered selecting those with a differential effect size of greater than or equal to 2-fold. This resulted in a progression signature set of 4662 reporters, 3021 genes (121 primary, 2900 metastatic, primary ∩ metastatic = 102). Signatures derived from this primary dataset were subsequently applied to two additional datasets, one prostate dataset generated by Tomlins et al. [5] in a laser capture microdissection study (GSE6099) and another generated by Ramaswamy et al. [6] and representing multi-tissue primary tumours and metastases (accessible through the Broad Institute data repository: http:// www.broadinstitute.org/cgi-bin/cancer/datasets.cgi).

Identifying correlated gene modules
We clustered our progression gene set using hierarchical clustering with a Ward agglomerative method designed to minimize intra-cluster variance (hclust, Bioconductor) and a 1 -Pearson correlation coefficient dissimilarity measure. We found this method produced a more highly correlated clustering structure when compared to other methods leading to more compact sub-clusters (Additional file 1). We characterised correlated gene modules by cutting the cluster dendrogram at branch lengths ranging from log10(0.05) to log10(3000) giving 39 equal intervals across the log scale. We removed clusters containing less than 3 members from further analysis. We selected modules defined at branch lengths 0f 0.6, 0.8, 1.1, 1.9, 2.5, 4.5, 10.6, 24.7 and 101.6 for further analysis since these gave a broad range of cluster numbers. Since a smaller branch length threshold does not always sub divide a parent module modules can be duplicated at different thresholds. These were removed from further analysis assigning them to the largest branch threshold at which they appears. We assigned Gene Ontology classifications to modules by testing for enrichment at GO nodes using a hyper-geometric distribution and a 0.01 p-value threshold. We carried out this analysis at the gene level by translating chip reporter probeset ids to Entrez gene ids. All reporters from the progression signature with assigned Entrez gene ids were used as background. Analysis was carried out using the GoStats package, Bioconductor.

Phenotype dependent transcript module expression changes
To determine differential regulation of modules within other expression datasets we first identified phenotype dependent expression changes for each sample using an absolute fold change filter of greater than 2. To generate fold changes against which we could filter each gene was scaled to a baseline intensity value. In the case of dataset GSE3325 each signal intensity from the primary tumour samples was scaled to the corresponding median gene signal intensity across the benign tumour samples. Likewise all metastatic samples were scaled to the median across all primary tumour samples. Prior to mapping modules to dataset GSE6099 we background corrected each sample using a normexp method and print tip loess normalised (normalizeWithinArrays(), Bioconductor). We then scaled PIN samples to EPI_ADJ_PCA control samples, PCA samples to PIN, MET_HNF to PCA and MET_HR to PCA. To identify hormone refractory dependent expression changes MET_HR samples were scaled to the median across the non-refractory samples MET_HN. To determine module induction or repression within the scaled samples we tested for enrichment of module genes within the sample associated expression changes using a hypergeometric distribution, <= 0.05 fdr. Mapping was achieved across array platforms using NCBI Entrez gene ids.
Modules with an intersection of less than 3 were discarded from the analysis.

Phenotype segregation
To determine if any of the enriched modules were capable of segregating samples on phenotype we built contingency tables across clinical conditions from each of the data sets (Tomlins [5] and Ramaswamy [6]) for induced and repressed modules and tested for sample enrichment using a Fisher's Exact test. Here we tested each phenotypic group against all others from each data set.

Cluster analysis
To determine the best clustering method and branch length thresholds to apply to our analysis we clustered our 4662 reporter prostate tumour progression signature (Additional file 2: Table S1) using single, average, complete, Ward's 1 minimum variance method and mcquitty 2 agglomerative hierarchical clustering methods along with a divisive method. The hclust function from Bioconductor 3 was used for the agglomerative techniques and the diana function from the cluster package from Bioconductor was used to run the divisive method. We used the cutree function at branch length thresholds ranging from 0.05 to 2 in increments of 0.05 to derive groups of correlated genes. In the case of the Ward agglomerative method where the branch scales [ it is unclear how Ward calculates its branch lengths, need to find out ] branch length thresholds ranged from log 10 (0.05) to log 10 (3000) in increments of log 10 (3000/0.05)/ 39. Branch length threshold intervals were chosen to produce a broad range of cluster numbers.

Cluster correlation
To assess the extent to which genes assigned to clusters are correlated we calculated a within-cluster dissimilarity value for each cluster 4 . This is given by where d x i , x j is the dissimilarity between genes i and j across all samples, i,j = 1,2,…N where N is the total number of cluster members and k = 1,2,…,K where K is the total number of clusters. In our case the dissimilarity measure is 1 -Pearson correlation coefficient between 2 genes across all samples. W(C) dissimilarity values across the array of branch length thresholds can be seen plotted against cluster number in Additional file 1: Figure S1. As observed the Ward agglomerative method out performs all other methods producing clusters that are less dissimilar and therefore more highly correlated than those generated from other methods relative to the number of clusters produced. These results provide a justification for using hierarchical clustering with a ward agglomerative method to generate sets of co-regulated genes.

Cluster gene ontology entropy
To quantify the information content of our clusters from a biological perceptive we assigned GO terms to cluster members. This was achieved by mapping GO terms via reporter entrez gene id assignments using the GO.db annotation package from Bioconductor. To quantify the GO information content of a cluster we calculated Shannon Entropy bit values given by: where x i is a cluster associated GO term, p(x i ) is the probability of choosing x i from all cluster GO terms , i = 1,2,…N where N is the total number of unique cluster GO terms, k = 1,2,…K where K is the total number of clusters and b = 2. H(X) bit values for different branch length thresholds can be seen plotted against cluster number for the different clustering techniques in Additional file 3: Figure S2. As observed the Ward clustering method produces clusters with higher GO bit values when compared to other methods. This implies greater uncertainly in the GO term mappings for clusters generated by the ward method thus indicating the production of clusters more GO information rich when compared to other methods. This provides further justification for using hierarchical clustering with a ward agglomerative method to generate sets of co-regulated genes.

Visualization of gene signatures through heatmaps
For visualization, sample groups were averaged using the mean prior to high level mean and variance normalization using the freely available software J-Express 2012 (http://jexpress.bioinfo.no/site/). Subsequently, both sample groups and genes were hierarchically clustered using complete linkage and Euclidian distance using the freely available software Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/software.htm). Heatmaps were produced using Java Tree View (http://jtreeview.sourceforge.net/). To evaluate the prediction performance 33 gene signature we analysed data from a microarray experiment of Grasso et al as available from GEO (GSE35988) http:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE35988.

Evaluating gene signature specificity and sensitivity
The dataset includes measurements on two different microarray platforms GPL6480 and GPL6848 and includes three different tissue types (benign, localised and metastatic castrate resistant prostate cancer).
We used the samples typed on platform GPL6848 as trainings data to derive the weights of the genes in the multi-gene signature. First, we replaced missing values using the k nearest neighbor algorithm as implemented in the R package impute. The gene PCA3 was not measured on the microarray and the gene CRISP3 was not observed in more than half of these samples. Thus, these two genes were excluded from the signature. We estimated the weights of each gene using an L2 regularized logit regression model [7] with the R package glmnet.
Then we used the samples typed on platform GPL6480 as test data to evaluate the prediction performance of the gene-signature. Per sample we computed one score as the weighted average over the 31 proposed genes where the weights were defined by the independent training data. Finally, we computed ROC statistics and report the area under the curve (AUC) of the ROC curve (R package ROCR). The AUC indicates the ability of a marker to distinguish between two groups, where a value 0.5 is random and a value of 1 represents a perfect distinction between the groups. Additionally, we were looking at the AUC of specific genes, in particular KLK3, ERG and AR.

Results and discussion
In order to define our starting signatures, we selected a dataset published by Varambally et al. about a decade ago and consisting of a small number of whole tissue sections [8]. This constitutes the smallest and oldest dataset used in our meta-analysis. It was, however, extensively validated at both the transcript and protein level in the original study and therefore provides a high-degree of confidence in data quality. We chose to define our starting signatures using this dataset in order to assess how much information could be derived despite the limitations in size and age. Within these data, we firstly identified transcripts that were differentially expressed in localised prostate versus benign tissue or in metastatic disease versus localised cancer using a conventional linear model approach. This approach identified 121 genes differentially expressed in localised primary cancers (primary versus benign, 0.01 FDR) and 2900 genes associated with metastatic status (metastatic versus primary, 0.01 FDR), which were covered by 4662 probes in total (Additional file 2: Table S1). To further refine these gene lists into discrete signatures, we constructed a gene coexpression network using Pearson correlation coefficients and hierarchical clustering using the Ward agglomerative method (See Methods section).
A number of different correlation or dissimilarity metrics have been employed when constructing co-expression networks. To determine the correlation between the genes we used a Pearson correlation coefficient to construct a dissimilarity matrix across all affected samples in the prostate tumour progression dataset and all genes identified in the preliminary analysis. We then used hierarchical clustering to group the genes. There are a number of available agglomeration methods available each producing their own clustering structure. To determine the best agglomeration method to apply in constructing our expression modules, we clustered our prostate tumour progression signature using single, average, complete, the Ward [9] minimum variance method and the Mcquitty [10] agglomerative hierarchical clustering method along with a divisive method. The performance of these clustering methods by using an algorithm to determine the extent to which genes assigned to clusters are correlated generating a within-cluster dissimilarity value for each cluster (Methods section -'Cluster Correlation'). In addition, we assessed the information content in gene ontology terms associated with clusters generated using each method by calculating Shannon Entropy bit values (Methods section -'Cluster Gene Ontology Entropy'). Shannon entropy and coefficient of variation are well known in a great many application domains, from theoretical physics to computational chemistry to materials science. They have been applied in bioinformatics as well, most notably in statistical genetics and molecular biology. Shannon entropy is derived from information theory [11]. Most relevant for this study the approach has previously been used as a measure of the robustness of gene regulatory networks [12], to accelerate feature elimination when classifying microarray expression data [13]. By these measures the Ward clustering method provided both more tightly associated coexpressed gene clusters as well as clusters with higher GO bit values when compared to other methods, indicative of greater information content in the ontologies derived for coexpression clusters generated using the Ward approach than using the other approaches. Additional file 4: Table S2 provides a complete list of coexpressed genes signatures generated used the Ward approach at all branching thresholds.  -06 CAV2, CAV1, MYL5, MYL2, TNC, PTEN, MYL9, VCL, IGF1R, LAMB3,  LAMB2, ITGB8, ILK, ITGB6, PDGFC, PAK1, THBS1, THBS4, COL4A4,  PRKCA, ACTB, MET, ITGA1, ACTN1, IGF1, HGF, COL4A6, FLNA,  LAMA4, ITGA6, CCND2, ITGA5, JUN, ITGA8, COL1A2,        Four large gene signatures were generated at the least stringent cut-point consisting of 1334 genes referred to as signature 1 (annotated as 101.6.1: Additional file 5: Table S3), 652 genes referred to as signature 2 (annotated as 101.6.2: Additional file 6: Table S4), 836 genes referred to as signature 3 (annotated as 101.6.3: Additional file 7: Table S5) and 357 genes referred to as signature 4 (annotated as 101.6.4: Additional file 8: Table S6). Signatures 1 and 2 contained genes that predominantly discriminated between localised PCa and benign tissue. Using DAVID ontology enrichment search (http://david.abcc.ncifcrf.gov/) to determine whether KEGG pathways were enriched within these signatures, we identified focal adhesions (hsa04510: Focal adhesion, p-value 7.21x10 −6 ) ( Table 1) as the most significant pathway for signature 1and complement and coagulation cascades for signature 2 (hsa04610: Complement and coagulation cascades, p-value 0.002) ( Table 2). The 38 genes associated with the focal adhesion annotation (p-value 7.21x10 −6 ) in signature 1 are listed in Table 1 and were all significantly downregulated in metastatic samples relative to benign and localised PCa. Half of these genes were laminins (eg. laminin alpha subunit-4 (LAMA4)), integrins (eg. integrin, alpha 1 (ITGA1) and five others), thrombospondins (thrombospondins 1 and 4 (THBS1/4), collagens, actins and myosins which may reflect the remodelling of the extracellular matrix and loss of stroma in particular during the transition to metastasis. The enrichment for complement and coagulation cascades in signature 2 (p value 0.002) included complement (eg. C1R, C1QA, C3) and plasma factors as well as serpin peptidase inhibitor as listed in Table 2 and were also predominantly downregulated in metastatic cases versus benign tissue and localised PCa. Collectively, these pathway enrichments might reflect a combination of extracellular matrix changes and the contribution of infiltrating immune cells and the inflammatory response. However, given that the Varambally dataset consists of whole-tissue sections it is not possible in this meta-analysis to precisely attribute these signatures to a particular biological process.
By contrast, signatures 3 and 4 contained genes that predominantly discriminated between metastatic cases and benign tissue samples. The dominant pathway for signature 3 was cell cycle regulation (hsa04110: Cell cycle, p-value 9x10 −20 ) and the enrichment arose from the overexpression of a total of 36 genes linked to this process in the metastatic cases versus benign tissue. The genes are listed in Table 3 and included E2F transcription factors, DNA replication licensing factors, cyclin-dependent kinase inhibitors, cell division cycle genes and components of the mitotic spindle checkpoint control apparatus. Many of these overexpressed genes also constitute a prognostic cell cycle progression gene signature, which has been validated at the transcript level in biopsy samples [14]). For signature 4, steroid biosynthesis was the most enriched pathway (hsa00100: Steroid biosynthesis, p-value 0.03squalene epoxidase (SQLE), farnesyl-diphosphate farnesyltransferase 1 (FDFT1), sterol-C4-methyl oxidase-like gene (SC4MOL). In this case the enrichment was due to the differential expression of three genes that are functionally tightly linked in some cases on consecutive steps in the cholesterol biosynthesis pathway. FDFT1 was overexpressed, SC4MOL was downregulated and SQLE showed a switch in expression in which one probe on the array was repressed and another was overexpressed (Table 4). Downregulation occured predominantly in localized PCa relative to benign tissue and expression seemed higher in metastatic cases than localized prostate cancers. FDFT1 overexpression, and increases in the expression of one probe for SQLE, were most significant in the metastatic cases compared to benign tissue and localised disease. These are enzymes associated with cholesterol biosynthesis in particular and collectively catalyse 3 out of 4 consecutive reactions in the conversion of farnesyl pyrophosphate to lathosterol via squalene. FDFT1 catalyses the production of squalene from farnesyl pyrophosphate, SQLE catalyses the conversion of squalene (See figure on previous page.) Figure 1 Gene signatures capable of discriminating between prostate cancer subgroups and classify metastatic disease. Gene signatures generated using the Varambally dataset and found to be significant discriminators of metastatic disease and primary/localised cancers (Additional file 10: Table S8) when applied to the Tomlins and Rawaswamy datasets were used to cluster samples in these datasets in a heatmap. The gene signatures represented are those capable of characterising samples from at least one progression stage (Fischer's exact < = 0.05). Gene signatures are rows and samples are columns. The colour coded bar at the base of the heatmap indicates the clinical grouping for each sample as also defined in the key. Metastatic hormone refractory, metastatic hormone naïve and hormone refractory vs. naïve represent prostate cancer cases from the Tomlins dataset, as do PIN (prostatic intraepithelial neoplasia) and primary carcinoma. The other categories (metastatic and primary) are samples from the Rawaswamy dataset and are metastatic and primary cancers from multiple organ sites, not simply the prostate gland. The blue bar graph on the right-hand side of the heatmap depicts the number of genes in each signature which are differentially expressed and contribute to the sample clustering in this analysis. For signature 1 (dist 101.6.1 and Additional file 5: Table S3) this is 1748 genes in total as highlighted and other bars are numbers of genes relative to this. The colour scale represents the mean log2 fold change for differential gene signatures (> = abs log2 (2)). Red indicates module induction, green repression. Gene signatures significant in both directions are indicated in yellow. Using the mean module log2 fold change we clustered the samples and modules using hierarchical clustering with euclidean distance as a measure of dissimilarity. Data points that contained both induced and repressed values have been excluded from the clustering.  to 2,3-epoxysqualene and SC4MOL catalyses the conversion of lanosterin to lathosterol. The two metabolites consecutively further downstream in the pathway are dehydrocholesterol and cholesterol. FDFT1 overexpression has previously been associated with aggressive PCa [15]). This is particularly intriguing since metastatic PCa is characterised by increases in the proliferative index of tumours [16] and the ability to produce autocrine steroid hormones from cholesterol in order to maintain androgen receptor activity [17]. Consequently the observation of increased levels of these enzymes in metastatic cases may hypothetically imply enhanced cholesterol biosynthesis to sustain its use for steroid hormone biogenesis by the tumours. Discrimination between cancer and benign control tissue and also between metastatic disease and other clinical cases represents an important goal of biomarker research. Thus, we used these gene signatures to classify clinical samples in prostate cancer samples and metastatic tissue samples in two additional datasets. One consisted of prostate cancer samples isolated by laser A B capture microdissection generated by Tomlins et al. [5] and the other contained expression array data from primary and metastatic tumours from multiple tissue sites generated by Ramaswamy et al. [6]. The Tomlins dataset consisted of various refined subgroups based on isolation of cell sub-populations including stromal fractions, epithelial fractions, localised prostate cancer and hormone-naïve and refractory metastatic disease. The Ramaswamy et al. dataset consisted of cancers from 14 organ sites with paired normal samples as well as normal tissues. In each dataset, we asked whether our signatures and sub-signatures could discriminate between the sample groups. To determine this, we first assessed the mean foldchange in the expression of each gene signature in each sample group in both datasets (Additional file 9: Table S7). We then performed a Fischer's Exact test to identify signatures that were capable of discriminating between localised prostate cancer, metastatic prostate cancer and the other sample groups defined in all three published studies -  Table S8). Gene ontologies were assigned to these statistically significant gene clusters and the clustering is represented in a heatmap for the classifying modules combining both the Tomlins and Ramaswamy sample sets (Figure 1 and Additional files 10 and 11: Tables S8 and S9 for gene ontology annotations). The smallest gene signature (dist.0.6.34) capable of subclustering localised prostate cancer from other samples in all three datasets consisted of 71 genes (Additional file 12). This small signature was a sub-component of the original signature 1 (101.6.1). The most significantly enriched biological process associated with these genes was vascular smooth muscle contraction (hsa04270: Vascular smooth muscle contraction, p-value 2x10 −3 ) ( Table 5). The four genes within this signature that were individually most significantly overexpressed in localised prostate cancers compared to benign tissues and metastatic cases were an oncogenic transcription factor, v-myc avian myelocytomatosis viral oncogene homolog (MYC), a proteoglycan capable of sequestering transforming growth factor beta called fibromodulin (FMOD), a mitochondrial enzyme associated with fatty acid metabolism called glycine N-acyltransferaselike protein 1 (GLYATL1) and an extraneuronal monoamine transporter called solute carrier family 22 member 3 (SLC22A3). MYC has been shown to be overexpressed in prostate cancer [18] and to drive tumourigenesis in a transgenic model of the disease [19]. Fibromodulin has not been widely studied in cancer and has not been implicated in prostate cancer. It is, however, known to be significantly overexpressed in chronic lymphocytic leukemia (CLL) versus normal B lymphocytes [20] and associated with a resistance signature to DNA damage-induced apoptosis in CLL [21]. Furthermore, the expression of fibromodulin is known to be induced in leiomyoma in response to TGF-beta through Smad and MAP kinase signalling [22]. GLYATL1 has not been associated with cancers. SLC22A3 has been reported to be overexpressed in localised prostate cancer at the transcript level when compared to benign tissue [5].
The other genes within this coexpression signature were downregulated in prostate cancers versus benign tissue and the majority were myosins, such as myosin, heavy polypeptide 11, smooth muscle (MYH11), myocardin (MYOCD), and myosin, light chain 9, regulatory (MYL9)thus accounting for the pathway enrichment for vascular smooth muscle contraction. As prostate cancer progresses to more advanced stages there is a depletion of stromal cells from the tissue and this perhaps explains the dominant contribution from downregulated muscleassociated genes to the signature and also other features of pathway enrichments particularly of the focal adhesion classification [23]. In order to determine whether our signature was consistent across more recent datasets, we downloaded an exon-array dataset generated by Taylor et al., and also The Cancer Genome Atlas (TCGA) data recently generated using high-throughput transcript sequencing of prostate cancers [24] (data generated by the data generated by the TCGA Research Network: http:// cancergenome.nih.gov/). MYC and GLYATL1 remain significantly overexpressed features (>1.3 fold) within these signatures in both datasets (Figure 2) with the vast majority of other gene transcripts downregulated including those enriched in the KEGG pathway analysis for vascular smooth muscle contraction. Whilst our 71 gene signature mainly contains differentially expressed genes that are downregulated in cancers versus benign tissues, most prostate cancer biomarkers that are currently under evaluation are overexpressed transcripts and proteins in the disease state. Consequently, we next sought to evaluate genes that were overexpressed in localised prostate cancers in signatures 1-4 more thoroughly in other datasets. There were 97 annotated gene transcripts in total overexpressed (Additional file 13). We had previously performed pathway analyses on signatures 1-4 which included both up-and downregulated genes (Tables 1, 2, 3 and 4). We now repeated this solely for the 97 overexpressed genes and this yielded pathway enrichment for O-glycan biosynthesis (hsa00512: O-glycan biosynthesis, p-value 0.009) as the most significant KEGG enrichment (Table 6).
To further refine this gene set, we then interrogated the Oncomine compendium of expression array data to determine which of these 97 genes are significantly overexpressed in at least three additional independent prostate cancer datasets when a Top 1% overexpression threshold was applied together with a p-value threshold of 1x10 −4 [25]. Thirty three annotated genes, around one-third of the 97-gene set fulfilled these criteria. This included 3 of the 4 overexpressed genes (MYC, GLYATL1 and SLC22A3) in the 71 gene signature subgrouping prostate cancers in the Varambally, Tomlins and Ramaswamy studies (highlighted in bold in Additional file 13). This 33 gene set also included four of the five glycosylating enzymes (UDP N-acetylglucosamine pyrophosphorylase 1 (UAP1), glucosaminyl (N-acetyl) transferase 1, core 2 (GCNT1), beta-1,3-glucuronyltransferase 1 (B3GAT1) and RAP1 GTPase activating protein 2 (RAP1GAP2/ GARNL4)) contributing to the ontology enrichment for glycan biosynthesis in the larger 97 gene set. Notably, others and we have recently reported that UAP1 and GCNT1 are overexpressed in prostate cancer tissue using immunohistochemistry. In addition, an aminosugar conjugate, O-linked N-acetylglucosamine (O-GlcNAc), is also significantly elevated in prostate cancer [26,27]. Furthermore, the UAP1 transcript has also been reported to be detectable in urine and plasma samples as a component of a multi-gene signature [28]. Additionally, UDP sugar conjugates have been identified as elevated in prostate cancers through metabolomics and O-linked N-acetylglucosamine is an overexpressed prostate cancer tissue biomarker, which can be conjugated to a variety of proteins to affect their stability and activity including c-Myc [29,30]. Consequently, the presence of these genes encoding glycosylating enzymes within this signature has been partly validated in tissue at the proteins level and suggests that more systematic profiling of glycoproteins may reveal new biomarkers.
Biologically, it is interesting to consider what might contribute to the increased expression of these genes in prostate cancers. Prostate cancer is driven by the dysregulated expression and activity of a number of transcription factors. The most notable example is the androgen receptor but others are overexpressed through chromosomal rearrangements and gene fusions as well as copy number variation as prostate cancer develops and progresses. This in turn has a significant impact on the expression of gene targets for these transcription factors and makes it plausible that a proportion of overexpressed genes reflect changes in transcription factor expression and activity. In this context, it is noteworthy that a total of five transcription factors were present in this group of 33 annotated genes ((single-minded family bHLH transcription factor 2 (SIM2), MYC, distal-less homeobox 1 (DLX1), homeobox C6 (HOXC6) and v-ets avian erythroblastosis virus E26 oncogene homolog (ERG)).
c-Myc is a well-established oncogenic transcription factor, which is overexpressed through chromosomal amplification on 8q24 but also through post-translational events, which may include glycosylation of the N-terminal transactivation and concomitant antagonism of proteasomal degradation [18,31,32]. ERG is part of a highly prevalent gene fusion affecting chromosome 21 and driven by the activity of the AR [33]. It is overexpressed in around 50% of prostate cancers through a chromosomal rearrangement, which fuses it to the upstream androgen receptordependent regulatory element controlling TMPRSS2 expression. SIM2 overexpression is associated with changes in transcriptional control affecting other loci on chromosome 21 [34,35]. Target genes for MYC and ERG have been extensively explored in clinical and cell-line datasets using expression array profiling with targeted knockdown and overexpression in prostate cancer cells. These approaches have linked MYC to processes including ribosome biogenesis and splicing and ERG to cell motility and migration, respectively [36][37][38]. Whilst the 33 genes did not include significant number of established MYC target genes, several reported ERG target were present including B3GAT1, phospholipase A1 member A (PLA1A) and collagen, type IX, alpha 2 (COL9A2) [39]. In addition, there were a number of direct AR targets including UAP1 and GCNT [30,40]. Importantly, whilst the AR is the principal transcription factor driving all stages of prostate cancer development, its target genes cannot be easily inferred by coexpression with the AR in contrast to ERG relative to ERG target genes. Target genes for HOXC6, SIM2 and DLX1 are less well defined in prostate cancers but given the presence of ERG and AR target genes within this geneset it is highly likely that they also contribute, being transcription factors, to the expression of some of these genes. A more systematic understanding of the interplay between these transcription factors and dependent gene networks will emerge in future studies. This will require targeting the expression of the transcription factors in experimental model systems and profiling concomitant changes in transcription factor recruitment, chromatin architecture and gene expression.
In the interim, however, it was possible to infer codependency based on co-clustering of genes in clinical samples. We did so in two additional datasets, an exon-array dataset generated by Taylor et al. and a transcriptomic dataset generated for prostate cancer through highthroughput sequencing by the The Cancer Genome Atlas (TCGA) (Figure 3). In both datasets, we were able to firstly reconfirm the ability of these 33 genes to discriminate between localised prostate cancer and benign tissue samples ( Figure 3). Secondly, ERG co-clustered within these 33 genes with bona fide target genes such as B3GAT1 and PLA1A corroborating a contribution at least from ERG to this prostate cancer-specific overexpression signature [39]. Intriguingly, another transcription factor, DLX1, also coclustered with ERG raising the possibility of a transcription factor hierarchy in which early emergence of an ERG gene fusion may trigger aberrant expression of other developmental transcription factors.
Currently prostate-specific antigen (PSA)/kallikrein 3 (KLK3) is the most widely used protein biomarker for prostate cancer. The androgen receptor (AR) is the most significant transcription factor driving prostate cancer, but is also expressed at high levels in untransformed luminal epithelial cells and therefore is predominantly used as a transcript biomarker associated with metastatic disease and concomitant copy-number amplification [41]. Gene fusions have been detected which significantly elevate the transcript levels of ETS transcription factors and the most prevalent example in prostate cancer is the TMPRSS2-ERG fusion [33]. Detection of the fusion has been reported in biological fluids including urine samples [42].  To assess the performance of the 33-gene signature in comparison to KLK3/PSA, AR and ERG we interrogated an additional independent expression array dataset generated by Grasso et al., and consisting of benign tissue, localised prostate cancer and metastatic cases [43]. This dataset was generated using two different array platforms on distinct sets of samples (Methods section). Cysteine-rich secretory protein-3 (CRISP3) was excluded from the signature due to missing values in the datasets for this gene and prostate cancer antigen 3 (PCA3) was not represented on the arrays leaving a 31-gene signature for evaluation. In the first phase of the signature evaluation we assessed the weighted contribution of each gene in the signature using a logistical regression model on a training dataset consisting of the samples profiled on an Agilent oligo microarray platform. We then used the samples profiled on a second platform, the 4x44K Agilent microarray to evaluate the performance of the signature and compared this to KLK3, ERG and AR. Three pairwise sample comparisons were undertakenbenign versus localised prostate cancer, benign versus metastatic cases and localised prostate cancers versus metastatic cases. Whilst all three transcripts and the signature discriminated between metastatic samples and benign tissue with good specificity and sensitivity as reflected in an area-under-the-curve (AUC) ranging from 0.83 for the AR to 0.99 for the signature, only the signature provided an AUC of ≥ 0.95 for all three pairwise comparisons (Table 7). Since both the AR and KLK3 are expressed in both untransformed prostate cells and prostate cancer it is perhaps not surprising that neither yielded an AUC of >0.65 in discriminating between localised prostate cancer and benign tissue samples. By contrast ERG expression is driven by a cancer-associated gene fusion and the AUC was 0.81 (Table 7). AR is amplified and overexpressed in metastatic prostate cancers and this likely explains the higher AUC for this marker (0.84) in discriminating metastatic cases from localised prostate cancers [41]. KLK3/PSA was also higher, 0.87, in this context. ERG by contrast whilst consistently overexpressed in the majority of localised prostate cancers is of variable utility as a prognostic marker according to the study cohort examined associating variously positively or negatively with progression and metastasis [44][45][46]. In our evaluation the AUC for ERG in discriminating localised prostate cancer from metastatic cases was 0.61, performing more poorly than as a discriminator of localised prostate cancer from benign tissue samples. The AUC differences between the markers and the signature in each pairwise comparison of the sample sets was also visualised in receiver operating characteristic (ROC) curves ( Figure 4). These comparisons highlight the importance of using a multi-gene signature since no single gene provides robust discrimination at all stages of the disease, no doubt reflecting changes in the underlying biological drivers during disease progression. We provide in addition AUC values for each individual gene and array probe for each of the pairwise sample comparisons in the test set (Additional file 14: Table S12) and the validation set (Additional file 15: Table S13). Although beyond the scope of this paper we hope that this will assist in further evaluation of the signature by researchers in the field.

Conclusions
In conclusion, in this study we have used a multi-step approach to refine gene signatures derived from diverse transcript detection platforms and sample types in order to arrive at a robust gene signature able to discriminate between PCa and benign tissue ( Figure 5). This is the first time that this has been attempted and demonstrates that value exists in transcript signatures generated from amongst the earliest microarray studies right through to high-throughput sequencing. In brief, beginning with a small expression array dataset consisting of 13 macrodissected samples, we have been able to derive gene signatures capable of subclustering localised PCa and metastases in a larger microdissected and a multi-cancer dataset ( Figure 5). This highlights that there are valuable gene transcript signatures that can be robust despite cellular heterogeneity in (See figure on previous page.) Figure 5 Workflow for the identification of robust gene signatures and gene sets for clustering prostate cancer cases. In step 1, we identified all statistically significant differentially expressed Affymetrix array probes in a small dataset consisting of 13 macrodissected clinical samples encompassing localised benign prostatic hyperplasia, localised prostate cancer and metastatic disease (GSE3325). We then generated gene signatures from these based on gene coexpression at varying stringency thresholds. These gene signatures were then applied to two additional datasets, a microdissected dataset (Tomlins et al.) and a multi-tissue site cancer and metastatic dataset (Ramaswamy et al.). A large number of the coexpression gene signatures clustered localised prostate cancers from metastatic disease and prostate metastases from other sample sets. The most compact gene signature able to do so consisted of 71 genes (A) and we assessed its expression pattern in two additional datasets, an exon-array dataset (Taylor et al.) and in a RNA-sequenced dataset (TCGA-PRAD). Few of the genes in the significant coexpression gene signatures were overexpressed genes in localised prostate cancers. In the second phase of the study, we abstracted all of the overexpressed genes and refined this down to a set of 33 genes based on significant overexpression in additional publicly available prostate cancer microarray datasets housed within the Oncomine database (B). These genes also effectively clustered benign versus cancer cases in an exon-array dataset (Taylor et al.) an expression microarray dataset (Grasso et al.) and a RNA-sequenced dataset (TCGA-PRAD) (C and D). In conclusion, it is possible to generate gene classifiers of clinical prostate cancer from a small dataset of macrodissected samples with the capacity to classify larger sequenced and microdissected datasets based on clinical characteristics.
PCa and the evolution of transcript detection technologies. In addition, we have discovered that gene transcripts that are significantly overexpressed within these signatures are also overexpressed in much more recently acquired exonarray and sequence-based TCGA data, transcript detection platforms that were unavailable when the Varambally, Tomlins and Ramaswamy studies were undertaken ( Figure 5). Finally, we have evaluated the performance of these transcripts as a signature in discriminating between benign tissue samples, localised PCa and metastatic disease in an additional dataset generated by Grasso et al. ROC curves reveal that the signature exceeds the performance of ERG, KLK3 or the AR as a classifier. Intriguingly, one third of these genes are glycosylating enzymes and transcription factors. PCa is significantly driven by a transcription factor, the AR, but there is increasing evidence of contributions by others and of interplays between them and indeed our signature does not include the AR itself. However, it includes both established examples (MYC and ERG) but also others that have so far been less studied (SIM2, DLX1 and HOXC6). Mechanistically, future work will investigate this transcriptional co-dependency in more detail and clinically these signatures will be further evaluated in clinical cohorts.