BACKGROUND
optimizing cancer treatment aims for a cure which kills all cancerous cells in the body with as little detriment to the patient as possible.
cancer is a highly heterogeneous disease with extreme genomic, intra- and inter-tumour heterogeneity; unsurprisingly, patients show a large variety in response to treatment  <cit> .
personalizing treatment is therefore expected to improve treatment response, and thus patient outcome.
for example, in some cases surgical resection of the tumour is curative; additional treatment, which has serious side-effects, is unnecessary.
in contrast, other patients presenting with similar clinical characteristics  could have more aggressive disease, for which adjuvant treatment is required to cure or control disease  <cit> .
without markers to distinguish these patients, all are given the same treatment, resulting in over-treatment in some patients and under-treatment in others.
to address this urgent clinical need, many groups have sought to create transcriptomic biomarkers using microarray-, pcr- or rna-seq-based assessments of mrna abundances.
the resulting multi-gene prognostic biomarkers  can identify patient subgroups that would be particularly likely to derive benefit from more intense therapy  <cit> .
however, there have been numerous challenges in the development of clinically-useful biomarkers; most published biomarkers fail to enter routine clinical practice  <cit> .
in cancer, where heterogeneity plays such an important role, these challenges are magnified; important tumour biomarkers may be missed when using the common practice of a single tumour biopsy to direct treatment.
if faced with uncertainty in biomarkers, these are deemed unsuitable for clinical applications and clinicians prefer to treat without the information and save costs  <cit> .
in order to advance personalized medicine, robust, reproducible biomarkers are required.
we have shown that, at least in lung cancer one of the major sources of biomarker irreproducibility is their sensitivity to relatively subtle changes in pre-processing  <cit> .
we found that analyzing a single biomarker with different pre-processing techniques yielded highly-divergent results, and these could indeed change clinical management for individual treatments  <cit> .
however, we also found tantalizing hints that different ways of analyzing a single biomarker could be integrated: an “ensemble” of pre-processing methodologies out-performed any individual one in a 442-patient cohort of non-small cell lung cancer patients.
it appears that each pre-processing technique removes a different aspect of the underlying noise in a dataset, and thus a large enough collection of them provides a more accurate estimate of the underlying biological signal.
to generalize and extend this finding, we explored the impact of data pre-processing on a micro-environmental biomarker problem: the prediction of tumour hypoxia.
tumor hypoxia  contributes to both inter- and intra-tumour heterogeneity, and can compromise cancer treatment.
it is a result of the uncontrolled growth of tumour cells and the formation of an abnormal tumour vascular network  <cit> , and is related to chemotherapy and radiotherapy resistance, tumour aggressiveness and metastasis  <cit> .
hypoxia is associated with poor prognosis  <cit> , and a marker for hypoxia both identify patients with more aggressive disease and those who might benefit from specific therapeutic options  <cit> .
many different predictors of hypoxia have been generated  <cit> .
to understand pre-processing sensitivity and how ensemble-classification can be best exploited, we evaluate this approach for  <dig> separate biomarkers in  <dig> datasets comprising transcriptomic profiles of  <dig>  primary, treatment-naïve breast cancers.
methods
datasets
the ensemble approach  <cit>  was applied to two separate groups of primary breast cancer datasets.
the first group comprises  <dig> datasets profiled on the affymetrix human genome u133a microarrays , with  <dig>  total patients  <cit> .
the second group is made up of  <dig> datasets profiled on affymetrix human genome u <dig> plus  <dig>  genechip array , comprising a combined  <dig> patients  <cit> .
only datasets reflected similar disease states and profiles were included, for example datasets of metastatic tumours were excluded  <cit> .
all samples included were treatment-naïve.
biomarkers
a series of  <dig> published hypoxia gene biomarkers were evaluated.
the following signatures were included: buffa metagene  <cit> , chi signature  <cit> , elvidge up gene set  <cit> , hu signature  <cit> , the 0% and 2% early seigneuric signatures  <cit> , sorensen gene set  <cit> , winter metagene  <cit>  and starmans clusters  <dig> to  <dig>  <cit> .
descriptions of each biomarker are given in additional file 1: table s <dig> and additional file 2: table s <dig>  the signatures evaluated here only contain up-regulated genes for which high gene expression is associated with poor survival.
pre-processing
all analyses were performed in the r statistical environment .
the first step was to pre-process each dataset in  <dig> different ways: all combinations of  <dig> pre-processing algorithms,  <dig> types of gene annotations and  <dig> approaches for dataset handling.
thus, each pipeline was defined by three factors .
each of these is outlined in detail in the following paragraphs.
the first factor creating pipeline variation for the ensemble classifier was the pre-processing algorithm.
we used robust multi-array average   <cit> , microarray suite  <dig>    <cit> , model-base expression index   <cit> , genechip robust multi-array average   <cit> .
all of which are available in the r statistical environment .
rma and gcrma return data in log2-transformed space whereas mas <dig> and mbei return data in normal space.
it is common practice to log2-transform mas <dig> and mbei pre-processed data, therefore both normal-space and log2-transformed versions of these two methods were included, giving us six pre-processing algorithms.
the second factor was the annotation approach.
a key part of microarray pre-processing involves mapping the  <dig> base-pair oligonucleotide probes to specific parts of the transcriptome .
this is accomplished using a chip description file .
our understanding of the human transcriptome is continually evolving, causing the annotation of individual probesets to change.
these advances are reflected in updated probeset annotation   <cit> .
therefore, we included both the “default" annotation  and updated entrez gene-based “alternative” annotation .
the number of probesets for each annotation is given in table  <dig>
the number of probe sets for each annotation and microarray platform after completion of pre-processing.
the last aspect of pipeline variation considered was dataset handling.
pre-processing was either done on each dataset individually or on all datasets merged into one.
separate dataset handling involves pre-processing of a single dataset as a unit, independent of others.
each separate dataset went through the pipeline and was classified independent of the other datasets.
from all separate datasets, patients classified as having good prognosis were pooled and patients predicted to have poor prognosis were pooled.
alternatively, for merged data handling, the cel files from all datasets were combined during pre-processing and went through the entire pipeline as one dataset.
univariate gene analysis
for each gene represented on both array platforms, patients were median dichotomized into low and high risk groups based on the signal-intensity of that gene across all patients for a single pipeline variant.
cox proportional hazards modeling was used to assess whether survival properties were significantly different between the low risk and high risk patients.
statistical significance was assessed using the wald test  and p-values were false-discovery rate  adjusted to correct for multiple-testing.
linear modeling
a simple linear model of platform, pre-processing algorithm, annotation method and dataset-handling type:
y=v+w+x+∑i=15zi
where y is the number of genes, v is the annotation method, w is the platform, x is the data handling and z is the pre-processing algorithm, was evaluated to determine if the model was a good fit for the data.
second, starting with a complete model of all pairwise interactions and main effects:
y=v+w+x+∑i=15zi+v:w+v:x+w:x+∑i=15v:zi+w:zi+x:zi
where y is the number of genes, v is the annotation method, w is the platform, x is the data handling and z <dig> .z <dig> specify the  <dig> options for the pre-processing algorithm, backwards stepwise refinement was performed using the akaike information criterion .
the linear modelling was constructed with alternative annotation as the baseline for v, hg-u133a as the baseline for w, merged data handling as the baseline for x, and gcrma as the baseline for z.
patient risk group classification
each gene signature was used to classify patients into one of two groups.
the number of genes present on each array for each annotation is shown in additional file 2: table s <dig>  after data pre-processing, a multi-gene signature score was calculated for each patient using all genes on that platform that are in the signature's gene list:
score=∑n=1ngeneexpr,n
where n is the number of genes in a signature and geneexpr,n is the median dichotomized value for the gene expression of the nth gene in the signature compared to the expression levels of that gene from all samples.
if the level of the nth gene is above the median for all samples then geneexpr,n is  <dig>  otherwise − <dig>
after calculating a score for each patient, these scores were used to median dichotomize patients into high and low risk groups for each signature.
ensemble classification
the patient risk group classifications across all pre-processing methods were combined to create an ensemble classification by looking for unanimous agreement between all pipeline variants.
the high risk classification for the ensemble classification is given to the patients who have been classified as high risk in all  <dig> pre-processing pipeline variants; similarly for the low risk grouping.
patients with conflicting classifications between pipeline variants were deemed to have unreliable molecular classifications and were thus excluded from ensemble classification as before  <cit>  as a conservative approach that might be used in the clinic.
individual classification for subset of patients
for better comparison between the ensemble classification and individual classifications, the number of patients classified based on one pre-processing approach was reduced to match the number of patients classified in the ensemble classifier.
instead of median dichotomization, the patients were ordered by their multi-gene signature score.
then the number of patients that the ensemble had classified as high risk was selected from the top of the order as high risk patients and this was equivalently done for the low risk classifications.
classifier evaluation
kaplan-meier survival curves and unadjusted cox proportional hazard ratio modeling  were used to assess survival differences between the low risk and high risk groups.
the wald test was used to determine whether the hazard ratio was statistically different from unity.
in all analyses, the superior classification was defined as the classification with the higher cox proportional hazard ratio.
permutation sampling for variable number of pipelines in the ensemble
in these analyses, the ensemble classification is generally a combination of all  <dig> pipeline variants.
however, we also varied the number of pipeline variants being combined.
to represent a combination of n pipeline variants, we randomly sampled n pipelines  and created an ensemble classifier as outlined above.
this process was repeated with replacement  <dig> times for each value of n ranging from  <dig> to  <dig>
student's t-test methods comparison
the pool of all  <dig> individual methods across the  <dig> signatures was split based on a single aspect of the pipeline .
we compared pipelines only differing on a single aspect using the paired t-test to assess statistical differences between pipelines.
permutation sampling for variable number of pipelines in the ensemble when subgrouping for methods comparison
as part of the method comparison, the pipelines where subgrouped based on a single aspect of the pipeline and then within the subgroups ensembles of a varying number of the pipelines were constructed.
to represent a combination of n pipeline variants, we sampled n pipelines  and created an ensemble classifier.
for each value of n , all possible combinations containing n unique pipeline variants were created.
visualization
all plotting was performed in the r statistical environment  using the lattice , latticeextra , rcolorbrewer  and cluster  packages.
RESULTS
ensemble classification approach
each dataset was pre-processed using  <dig> different pipeline variants.
each biomarker was then applied separately for each pipeline variant, producing an ensemble of  <dig> predictions for each patient and biomarker.
these were analyzed for consistency and combined to form a single ensemble classification.
figure  <dig> outlines the approach used.
we separated our datasets according to the microarray platform used, and tested the two most widely-used platforms at the time of writing according to depositions in the gene expression omnibus: hg-u133a and hg-u <dig> plus  <dig> .
since both platforms are affymetrix arrays and therefore have the same set of potential normalization methods, we can perform inter-platform analysis independent of pre-processing.
univariate gene analysis
we first investigated the univariate performance of individual genes to determine how the prognostic power of these simple biomarkers is influenced by pre-processing differences.
as shown previously for lung cancer  <cit> , the prognostic ability of individual genes varied dramatically across methods.
of the  <dig>  genes represented on both array platforms tested, 74% reached statistical significance after multiple-testing correction in at least 1/ <dig> pipeline variants.
by contrast, only 16% reached significance in at least 12/ <dig> pipelines  and none were significant in all pipelines.
three pipeline variants identified zero genes, while three others found a single gene , which was not identified in the other  <dig> pipelines.
these data clearly indicate that simple union  and intersection  approaches are inappropriate.
interestingly, all six pipelines that resulted in either one or no prognostic genes involved analysis of hg-u133a data , using either the rma or mbei algorithms, along with the “separate” dataset-handling approach.
there is an evident difference between the patterns of significant genes on each platform.
the lowest concordance between pipelines is shown in the inter-platform correlations.
different aspects of the pipeline appear more highly correlated depending on the platform and there is no clear ordering of which aspect is more important without interactions .
we were able to use linear-modeling to show that the choice of pre-processing method is strongly deterministic for the number of statistically-significant genes identified.
we considered a complete model of all pairwise interactions and main effects, then used the akaike information criterion  for backwards stepwise refinement.
a model containing the main effects: platform, pre-processing algorithm, data-handling type and their pairwise interactions resulted , indicating that the relationship is deterministic, not stochastic.
we note that interactions are critical: a simple model of main-effects was not explanatory .
for the linear model, y=w+x+∑i=15zi+w:x+∑i=15w:zi+x:zi where y is the number of genes, w is the platform, x is the data handling and z1…z <dig> are specify the  <dig> options for the pre-processing algorithm, the coefficients that have a p <  <dig>  are shown.
multi-gene signatures
we next focused on multi-gene classifiers, seeking to determine if our single-gene results could be generalized.
we compared the hazard ratios from cox modeling of the ensemble and the  <dig> individual classifications for  <dig> published hypoxia signatures.
for all multi-gene signatures, superior classification was defined as the classification with a higher hazard ratio.
as seen with the single gene classifiers, variation was observed between classifications from the different pipelines and there was not one single variant which consistently resulted in larger risk stratification than the others.
further this analysis identified microarray platform as another possible source for variation.
one pipeline variant  showed the lowest risk stratification of the  <dig> pipelines on one platform  and the largest of the  <dig> pipelines on the other platform  .
as shown in figure  <dig>  ensemble classification performed better than individual pipelines and improved signature performance for both microarray platforms.analyses for all signatures showed that performance was sensitive to pre-processing choices and, in the majority of cases, the ensemble classification improved prognostic ability over individual pipeline variants .
for half of the signatures, ensemble classification resulted in superior risk stratification  compared to classifications from the individual pre-processing pipelines.
moreover the ensemble technique was almost always superior to the “typical” pre-processing techniques, exceeding the median of the  <dig> techniques in 24/ <dig> signature comparisons.
the buffa metagene and the winter metagene showed similar results across pipeline variants, but many of the signatures performed very differently depending on the dataset platform .
overall signatures showed better risk-stratification on hg-u <dig> plus  <dig>  arrays , although this was signature-specific.
some signatures  showed consistently better results on the hg-u <dig> plus  <dig>  dataset compared to the hg-u133a dataset.
conversely, starmans cluster  <dig> and cluster  <dig> performed better in the hg-u133a datasets.
the buffa and the winter metagene were the only signatures which were statistically significant across all pipelines tested.
hu and sorensen, additionally, were other signatures with statistically significant ensemble classifications for both datasets.
in contrast, starmans clusters  <dig>   <dig>   <dig> and seigneuric 0% early signatures did not perform well in either dataset; none of their ensemble classifications were statistically significant.
in general, if a signature performed poorly for single pipeline variants, using the ensemble classification did not improve it.
this was demonstrated by the correlation between the hazard ratios for the ensemble classification and the maximum hazard ratios for classification from the individual pipeline variants .
since previous analyses involved comparing unequal numbers of patients classified, we also compared ensemble classification to classification for the individual pre-processing methods.
in this way, we match patient numbers between the two conditions, removing this potential confounding variable.
in general, this approach yielded fewer statistically significant results , although both the range and the variance of hazard ratios increased for every signature using this classification algorithm.
however the comparison between of ensemble classifications and individual classifications shows that patient-number differences are not the origin of the superior performance of ensemble classification.
for 13/ <dig> signatures, the ensemble classification was superior to all classifications from the individual pre-processing pipelines and in 26/ <dig> signatures the ensemble exceeded the median classification.
signature comparison
to better understand which signatures were more successful, all individual classifications were compared.
unsupervised clustering of the percentage agreement of concordant patient classifications between individual pipeline variants for each signature showed that they mainly clustered by signature, rather than by pipeline composition .
this indicated that, although pre-processing substantially influenced biomarker performance, the genes in the signature characterized the overall partition and determined whether it was a poor or good biomarker.
the buffa metagene had the most consistent patient classifications across pipelines, but hazard ratios still ranged from  <dig>  to  <dig> .
although, we evaluated only hypoxia signatures, patient classifications did not agree across signatures .
signatures of ensemble classifications that were statistically significant generally classified a larger fraction of patients .
what is the optimal ensemble size?
having shown that the ensemble-approach improved classification for most biomarkers and datasets, we explored the limits of its performance.
we wondered if  <dig> distinct pipelines were always necessary, and therefore evaluated the number of pipeline variants required for optimal performance  of the ensemble classifier.
if creating an ensemble of four pipeline variants is equally successful to one from eight variants, then it is not beneficial to introduce the complexity and computational-costs of pre-processing with four extra pipelines.
focusing on signatures with a significant 24-pipeline ensemble, different combinations of pipelines, ranging from combinations of only  <dig> to all  <dig>  were evaluated.
these analyses indicated that in general increasing the number of pipeline variants resulted in an increase in absolute effect size which started to plateau as the number of methods in the ensemble increased .
in parallel, the percentage of patients classified with the ensemble method decreased and plateaued .
most signatures shared the same shape but with different rates of hazard ratio increase.
the sorensen signature on the hg-u133a dataset plateaued at about four pipeline variants.
therefore, in this case, randomly choosing four pipeline variants to combine provided roughly the same risk stratification as using all  <dig> pipelines.
conversely, for the winter metagene signature in either dataset, the mean hazard ratio continued to increase all the way up to  <dig> pipelines, though the curve was steeper at the beginning then in the end.
although the hazard ratio stopped increasing in some cases, stability continued to increase as the number of methods in the ensemble increased.
this is demonstrated in additional file 8: figure s <dig> by the tightening of the hazard ratio range as the number of pipelines is increased.
considering the winter metagene signature in hg-u133a data, the ensembles created from nine or more of the  <dig> pipelines outperformed all single pipeline classifiers .
many ensembles did not require all  <dig> variants to be an improvement over all non-ensemble methods .
even if the ensemble of  <dig> variants was not an improvement over non-ensemble methods, there may still have been an ensemble of a subset of the variants which was superior to the non-ensemble methods .
these data provide a compelling rationale to consider and evaluate ensemble pipelines for all microarray-based biomarkers.
methods comparison
after showing that ensembles are beneficial, we wanted to look at whether we can determine the combination of pipelines that lead to higher hazard ratios in order to add the most benefit for each additional pre-processing pipeline.
there is a clear relationship between the number of patients classified in the ensemble and the gain in hazard ratio, meaning that the ensemble is choosing to exclude the right subset of patients .
methods that produce less-correlated classifications gain more from the ensemble classification.
however, if we look at which methods are diverse by a different metric such as the profiles of prognostic ability of each gene as a single gene classifier, there is only a slight but not obvious increase in hazard ratio from using more diverse pipelines in the ensemble classification .
to help direct pipeline choices, we sought to address whether certain aspects of the pipeline resulted in better or worse performance.
for each aspect of the pipeline , the hazard ratios were grouped per variant of that aspect and compared.
this was done for both platforms separately and combined.
on both platforms there was a significant difference between annotations.
on hg-u133a, alternative annotation had higher hazard ratios .
in direct contrast, hg-u <dig> plus  <dig>  performed better with default annotation .
by contrast, the optimal pre-processing algorithm was similar in both platforms, with rma and mbei performing better than gcrma and mas <dig> .
rma and mbei showed similar results  as did gcrma and mas <dig> .
furthermore, we analyzed the effect of changing the number of variants in the ensemble when creating only ensembles from common pipeline variants .
once again, variant success is not necessarily consistent across signatures.
ensemble of signatures
to further filter out unreliable classifications we investigated combining the classifications from two signatures.
the buffa metagene and the winter metagene performed best across our analyses.
these two signatures share  <dig> genes .
expansion of the ensemble classification to only classify patients that both signatures agreed on  improved risk stratification  compared to ensemble evaluations of both signatures .
to complete the analysis and expand the number of patients classified, we also pooled the unanimous classifications .
this failed to improve risk stratification compared to ensemble evaluations of both signatures; however, prognostic performance was improved over all the signatures’ individual pre-processing methods.
further, many more patients were classified than with the basic ensemble technique , suggesting that ensembles of signatures could be used to further remove noise or to increase the number of patients given confident molecular classifications.
discussion
the purpose of pre-processing is to remove “noise” from the data.
however, since no method is perfect, each pre-processing pipeline removes a somewhat different aspect of the “noise”.
indeed, groups around the world have focused on identifying the “optimal” pre-processing technique for different types of data  <cit> .
the principle of ensemble classification is that by combining pre-processing approaches we can select the parts of the data which are reliable across the multiple approaches.
the central tendency of this pool of methods is thus predicted to lie closer to the “true” value, and thereby to provide a better biomarker.
although different pre-processing methods may cause some variation in the analysis, pre-processing is expected to have a minor effect on the core experimental results and conclusions  <cit> .
our previous work has indicated this is not the case and pre-processing caused major outcome differences in non-small cell lung cancer  <cit> .
here we systematically extend and deepen these analyses to explore the variation caused by algorithmic diversity in pre-processing.
at the single gene level substantial differences in prognostic power were seen in univariate analysis.
therefore pre-processing is part of the reason different studies identify different biomarker genes.
many authors will use public data to show that a given gene is prognostic; however, essentially all genes  can meet that criterion, depending on which platform and pre-processing technique is used.
single genes did not appear to behave the same across pipelines demonstrating variation in classification results are expected and signatures are dependent on the pre-processing platform they were discovered on.
the fluctuation in prognostic ability for each gene partially explains why we get different results for multi-gene signatures and why there is such difficulty validating biomarkers between research groups  <cit> .
in combination with data showing a plurality of gene sets are associated with outcome in breast cancer and other cancers  <cit> , the variation suggests that mining of public data for prognostic significance is very prone to over-fitting and multiple-testing concerns.
therefore robust, permutation-based approaches need to be developed  <cit> .
in non-small cell lung cancer, starmans et al.
<cit>  showed one example of a permutation-based approach where an ensemble biomarker classifier improved survival separation between low risk and high risk patient groups.
here we extended this finding showing that the method replicates on two microarray platforms representing  <dig> separate datasets in breast cancer for a series of  <dig> biomarkers.
across both platforms, there was not a single pipeline that unfailingly outperformed all other pipelines; therefore, the ensemble classification provides a generalized approach to improve biomarkers, both in terms of performance and stability, without determining the actual optimal pre-processing pipeline.
furthermore, in many of the cases, the ensemble classifications outperformed all single pre-processing methods.
the ultimate value of the ensemble classifiers as a concept was demonstrated with the buffa metagene and winter metagene.
for these two signatures, any ensemble classifier comprising at least nine pipelines on hg-u133a or  <dig> pipelines on hg-u <dig> plus  <dig>  arrays generated superior risk stratifications compared to all the classifiers that used only a single pre-processing pipeline.
consequentially, ensemble classifiers are most definitely beneficial and should be used.
the ensemble approach did not improve all biomarkers.
biomarkers with generally bad risk stratifications across pre-processing pipelines still showed poor performance when combined in an ensemble.
the ensemble approach magnified the separation of risk groups rather than corrected for a poor initial biomarker.
therefore the ensemble approach can also be used as a metric to assess the quality of biomarkers, distinguishing between poor and good signatures.
by the statistical significance and consistency in risk stratification improvement across the datasets, the buffa and winter metagenes are shown to be strong, consistent signatures.
by the same metric, seigneuric 0% early and starmans cluster  <dig>   <dig>   <dig> appear to be poor signatures validating previous findings where these signatures did not show prognostic power  <cit> .
a disadvantage of the ensemble classification is that a fraction of patients are not classified.
only patients with robust risk classification across pipelines in the ensemble are assigned to risk groups.
here, using only the statistically significant ensembles 16% to 80% of the patients were not classified.
the signatures showing significance on both platforms tend to have a higher percentage of patients classified  than the signature significant on only on platform  and the signatures not significant on either platform .
nevertheless a patient classified as unreliable with one signature may be robustly classified using a different signature.
this was shown by intersecting the buffa and winter metagenes, which resulted in improved prognostic power compared to the single pre-processing pipelines and classified more patients than with the two ensembles individually.
one might consider taking this into account in biomarker-developement by attempting to construct ensembles that minimize the fraction of unclassified patients in the training dataset, although unclassified patients could resort to standard clinical care.
an important note about our approach to using the ensemble classification is the diversity of the pre-processing methods.
our choice of signature evaluation meant that a log2-transform on the data does not create different classifications.
for example, the multiple methods mas <dig> and log <dig> mas <dig> pre-processing are actually only one pipeline variant and aren't filtering out additional unreliable patients.
here, the array of  <dig> pre-processing methods in reality only gives  <dig> unique classifications so there is less diversity than the numbers indicate.
however, for other signatures, such as the risk-score and clustering methods evaluated previously  <cit> , up to  <dig> unique classifications are possible.
in multiple signatures, the ensemble of all  <dig> pre-processing pipelines is an improvement but not the optimal ensemble classification .
so future analysis to refine the process and finesse which pipelines to use and how they should be combined would be advantageous.
an important future direction is exploring how ensemble methods can be improved by incorporating greater algorithm diversity.
continual addition of diverse methods may increase the optimal ensemble classification or the optimal ensemble may be a certain combination of pipelines and new additions may not lead to an increase.
in some cases, such as the hu signature on hg-u <dig> plus  <dig>  addition of another separate data handling pipeline is not likely to increase the risk stratification but adding a merged data handling pipeline would be advantageous .
CONCLUSIONS
we systematically show that differences in pre-processing create differences when using biomarkers.
this effect of pre-processing is important for the research community to recognize and consider, as accurately accounting for it will advance biomarker discovery, validation and ultimately clinical application.
we found that the buffa metagene is the most consistent biomarker and therefore most clinical useful signature evaluated and we show that application of ensemble classification technique is beneficial for improving risk stratification both in terms of effect size and stability of biomarkers.
competing interests
all authors declare that they have no competing interests.
authors’ contributions
database generation and curation: sh.
performed statistical and bioinformatics analyses: nsf, mhws, pcb.
data interpretation: nsf, mhws, pcb.
wrote the first draft of the manuscript: nsf.
initiated the project: pcb.
supervised research: mhws, pl, pcb.
all authors read and approved the final manuscript.
supplementary material
additional file 1: table s1
prognostic signature descriptions.
click here for file
additional file 2: table s2
gene counts per prognostic signature.
click here for file
additional file 3: figure s1
correlation of gene univariate analysis.
analysis of consistency between methods for the prognostic ability of each gene shown in figure  <dig>  the heatmap shows pairwise comparison of all the pipeline variants where the comparison is spearman's correlation estimate of the fdr-adjusted p-values  for univariate cox proportional hazard ratio modeling analysis of genes analyzed on the set of pipelines.
click here for file
additional file 4: figure s2
platform comparison by signature.
comparison of hazard ratios for the series of prognostic signatures on hg-u133a and hg-u <dig> plus  <dig> .
hazard ratios were derived from cox proportional hazard ratio modeling.
each triangle represents the ensemble classifier's hazard ratio and the circles represent the individual pipeline variants.
the 95% confidence interval is shown for each ensemble.
for the individual pipeline variants, the 95% confidence intervals are shown in additional file 5: table s <dig>
click here for file
additional file 5: table s3
hazard ratio 95% confidence intervals for classifications on the individual pipeline variants.
click here for file
additional file 6: figure s3
risk stratification across classification pipelines and prognostic signatures with equal number of patients classified.
comparison of hazard ratios  and corresponding p-values from cox proportional hazard ratio modeling between ensemble classifications and individual classifications on a subset of patients with the highest and lowest signature scores on  hg-u133a platform,  hg-u <dig> plus  <dig>  platform.
the hazard ratio is represented by the size and colour of the dot and the background shade represents the p-value.
click here for file
additional file 7: figure s4
signature comparison.
analysis of consistency across both significant prognostic signatures and signatures that were not .
heatmaps are shown for the pair-wise comparison  of all the single pipeline classifications for the individual pre-processing methods  and the ensemble scores derived from these individual classifications per patient for each signature .
in b, the signatures are ordered by the number of patients classified unanimously across all the pipeline variants.
from left to right, the number of patients classified in the ensemble for each signature decreases.
click here for file
additional file 8: figure s5
ensemble hazard ratio range.
the range of hazard ratios for ensembles from different number of pipeline variants.
the horizontal pink dashed line shows the highest hazard ratio of the individual methods; all the ensembles above the line are improvements on current pre-processing practice.
the x-axis indicate the number of pipeline variants combined to create ensembles.
the grey background shows the numbers of pipeline variants where all the ensembles created are superior to every single individual method.
the hazard ratio, p-value and number of patients classified for each ensemble shown is provided in additional file 11: table s <dig> and additional file 12: table s <dig>
click here for file
additional file 9: table s <dig>
the hazard ratios, p-values and number of patients classified for all of the classifications on hg-u133a in additional file 8: figure s <dig>  to make the data easier to use, each signature is in a separate a tab delimited table/text file and files for each platform are packaged and compressed separately.
each row in the tables is a patient classification and there are repeated rows since we were sampling the pipelines with replacement.
the first  <dig> columns of each table is whether the pipeline specified in the column name is used in the ensemble classification with  <dig> meaning the pipeline is in the ensemble and  <dig> meaning it is not.
columns  <dig> –  <dig> are the hazard ratio, p-value and the number of patients classified for the classification respectively.
click here for file
additional file 10: table s <dig>
the hazard ratios, p-values and number of patients classified for all of the classifications on hg-u <dig> plus  <dig>  in additional file 8: figure s <dig>  to make the data easier to use, each signature is in a separate a tab delimited table/text file and files for each platform are packaged and compressed separately.
each row in the tables is a patient classification and there are repeated rows since we were sampling the pipelines with replacement.
the first  <dig> columns of each table is whether the pipeline specified in the column name is used in the ensemble classification with  <dig> meaning the pipeline is in the ensemble and  <dig> meaning it is not.
columns  <dig> –  <dig> are the hazard ratio, p-value and the number of patients classified for the classification respectively.
click here for file
additional file 11: figure s <dig>
method correlation effect on hazard ratio.
comparison of the effect of method diversity in ensembles of  <dig> on the increase in hazard ratio from the maximum of the individual classifications for winter metagene classifications on hg-u133a  and hg-u <dig> plus  <dig>  .
part a measures how correlated methods are by their percent agreement between methods  which is also equivalent to the number of patients classified.
part b measures the relatedness of the methods by the spearman's correlation of how prognostic each gene is for a method .
click here for file
additional file 12: figure s <dig>
combining signatures.
prognostic ability of combining the ensemble approach for the winter metagene and the buffa metagene was evaluated with kaplan-meier survival analyses.
hazard ratios and p-values are from cox proportional hazard ratio modeling.
the intersect is using only patients that are in agreement between winter metagene and buffa metagene.
the union is pooling the patients from winter metagene and buffa metagene .
click here for file
acknowledgements
the authors thank nathalie moon for advice on statistical analysis and all members of the boutros lab for helpful suggestions.
funding sources
this study was conducted with the support of the ontario institute for cancer research to pcb through funding provided by the government of ontario.
pcb was supported by a terry fox research institute new investigator award.
this work was supported by the center for translational molecular medicine  , eu imi program , eu 7th framework program  and dutch cancer society  through funding to mhws and pl.

 optimizing cancer treatment aims for a cure which kills all cancerous cells in the body with as little detriment to the patient as possible. cancer is a highly heterogeneous disease with extreme genomic, intra- and inter-tumour heterogeneity; unsurprisingly, patients show a large variety in response to treatment. personalizing treatment is therefore expected to improve treatment response, and thus patient outcome. for example, in some cases surgical resection of the tumour is curative; additional treatment, which has serious side-effects, is unnecessary. in contrast, other patients presenting with similar clinical characteristics could have more aggressive disease, for which adjuvant treatment is required to cure or control disease. without markers to distinguish these patients, all are given the same treatment, resulting in over-treatment in some patients and under-treatment in others. to address this urgent clinical need, many groups have sought to create transcriptomic biomarkers using microarray-, pcr- or rna-seq-based assessments of mrna abundances. the resulting multi-gene prognostic biomarkers can identify patient subgroups that would be particularly likely to derive benefit from more intense therapy. however, there have been numerous challenges in the development of clinically-useful biomarkers; most published biomarkers fail to enter routine clinical practice. in cancer, where heterogeneity plays such an important role, these challenges are magnified; important tumour biomarkers may be missed when using the common practice of a single tumour biopsy to direct treatment. if faced with uncertainty in biomarkers, these are deemed unsuitable for clinical applications and clinicians prefer to treat without the information and save costs. in order to advance personalized medicine, robust, reproducible biomarkers are required. we have shown that, at least in lung cancer one of the major sources of biomarker irreproducibility is their sensitivity to relatively subtle changes in pre-processing. we found that analyzing a single biomarker with different pre-processing techniques yielded highly-divergent results, and these could indeed change clinical management for individual treatments. however, we also found tantalizing hints that different ways of analyzing a single biomarker could be integrated: an “ensemble” of pre-processing methodologies out-performed any individual one in a 442-patient cohort of non-small cell lung cancer patients. it appears that each pre-processing technique removes a different aspect of the underlying noise in a dataset, and thus a large enough collection of them provides a more accurate estimate of the underlying biological signal. to generalize and extend this finding, we explored the impact of data pre-processing on a micro-environmental biomarker problem: the prediction of tumour hypoxia. tumor hypoxia contributes to both inter- and intra-tumour heterogeneity, and can compromise cancer treatment. it is a result of the uncontrolled growth of tumour cells and the formation of an abnormal tumour vascular network, and is related to chemotherapy and radiotherapy resistance, tumour aggressiveness and metastasis. hypoxia is associated with poor prognosis, and a marker for hypoxia both identify patients with more aggressive disease and those who might benefit from specific therapeutic options. many different predictors of hypoxia have been generated. to understand pre-processing sensitivity and how ensemble-classification can be best exploited, we evaluate this approach for separate biomarkers in datasets comprising transcriptomic profiles of primary, treatment-naïve breast cancers. methods datasets the ensemble approach was applied to two separate groups of primary breast cancer datasets. the first group comprises datasets profiled on the affymetrix human genome u133a microarrays, with total patients. the second group is made up of datasets profiled on affymetrix human genome u plus genechip array, comprising a combined patients. only datasets reflected similar disease states and profiles were included, for example datasets of metastatic tumours were excluded. all samples included were treatment-naïve. biomarkers a series of published hypoxia gene biomarkers were evaluated. the following signatures were included: buffa metagene, chi signature, elvidge up gene set, hu signature, the 0% and 2% early seigneuric signatures, sorensen gene set, winter metagene and starmans clusters to. descriptions of each biomarker are given in additional file 1: table s and additional file 2: table s the signatures evaluated here only contain up-regulated genes for which high gene expression is associated with poor survival. pre-processing all analyses were performed in the r statistical environment. the first step was to pre-process each dataset in different ways: all combinations of pre-processing algorithms, types of gene annotations and approaches for dataset handling. thus, each pipeline was defined by three factors. each of these is outlined in detail in the following paragraphs. the first factor creating pipeline variation for the ensemble classifier was the pre-processing algorithm. we used robust multi-array average, microarray suite, model-base expression index, genechip robust multi-array average. all of which are available in the r statistical environment. rma and gcrma return data in log2-transformed space whereas mas and mbei return data in normal space. it is common practice to log2-transform mas and mbei pre-processed data, therefore both normal-space and log2-transformed versions of these two methods were included, giving us six pre-processing algorithms. the second factor was the annotation approach. a key part of microarray pre-processing involves mapping the base-pair oligonucleotide probes to specific parts of the transcriptome. this is accomplished using a chip description file. our understanding of the human transcriptome is continually evolving, causing the annotation of individual probesets to change. these advances are reflected in updated probeset annotation. therefore, we included both the “default" annotation and updated entrez gene-based “alternative” annotation. the number of probesets for each annotation is given in table  the number of probe sets for each annotation and microarray platform after completion of pre-processing. the last aspect of pipeline variation considered was dataset handling. pre-processing was either done on each dataset individually or on all datasets merged into one. separate dataset handling involves pre-processing of a single dataset as a unit, independent of others. each separate dataset went through the pipeline and was classified independent of the other datasets. from all separate datasets, patients classified as having good prognosis were pooled and patients predicted to have poor prognosis were pooled. alternatively, for merged data handling, the cel files from all datasets were combined during pre-processing and went through the entire pipeline as one dataset. univariate gene analysis for each gene represented on both array platforms, patients were median dichotomized into low and high risk groups based on the signal-intensity of that gene across all patients for a single pipeline variant. cox proportional hazards modeling was used to assess whether survival properties were significantly different between the low risk and high risk patients. statistical significance was assessed using the wald test and p-values were false-discovery rate adjusted to correct for multiple-testing. linear modeling a simple linear model of platform, pre-processing algorithm, annotation method and dataset-handling type: y=v+w+x+∑i=15zi where y is the number of genes, v is the annotation method, w is the platform, x is the data handling and z is the pre-processing algorithm, was evaluated to determine if the model was a good fit for the data. second, starting with a complete model of all pairwise interactions and main effects: y=v+w+x+∑i=15zi+v:w+v:x+w:x+∑i=15v:zi+w:zi+x:zi where y is the number of genes, v is the annotation method, w is the platform, x is the data handling and z.z specify the options for the pre-processing algorithm, backwards stepwise refinement was performed using the akaike information criterion. the linear modelling was constructed with alternative annotation as the baseline for v, hg-u133a as the baseline for w, merged data handling as the baseline for x, and gcrma as the baseline for z. patient risk group classification each gene signature was used to classify patients into one of two groups. the number of genes present on each array for each annotation is shown in additional file 2: table s after data pre-processing, a multi-gene signature score was calculated for each patient using all genes on that platform that are in the signature's gene list: score=∑n=1ngeneexpr,n where n is the number of genes in a signature and geneexpr,n is the median dichotomized value for the gene expression of the nth gene in the signature compared to the expression levels of that gene from all samples. if the level of the nth gene is above the median for all samples then geneexpr,n is otherwise − after calculating a score for each patient, these scores were used to median dichotomize patients into high and low risk groups for each signature. ensemble classification the patient risk group classifications across all pre-processing methods were combined to create an ensemble classification by looking for unanimous agreement between all pipeline variants. the high risk classification for the ensemble classification is given to the patients who have been classified as high risk in all pre-processing pipeline variants; similarly for the low risk grouping. patients with conflicting classifications between pipeline variants were deemed to have unreliable molecular classifications and were thus excluded from ensemble classification as before as a conservative approach that might be used in the clinic. individual classification for subset of patients for better comparison between the ensemble classification and individual classifications, the number of patients classified based on one pre-processing approach was reduced to match the number of patients classified in the ensemble classifier. instead of median dichotomization, the patients were ordered by their multi-gene signature score. then the number of patients that the ensemble had classified as high risk was selected from the top of the order as high risk patients and this was equivalently done for the low risk classifications. classifier evaluation kaplan-meier survival curves and unadjusted cox proportional hazard ratio modeling were used to assess survival differences between the low risk and high risk groups. the wald test was used to determine whether the hazard ratio was statistically different from unity. in all analyses, the superior classification was defined as the classification with the higher cox proportional hazard ratio. permutation sampling for variable number of pipelines in the ensemble in these analyses, the ensemble classification is generally a combination of all pipeline variants. however, we also varied the number of pipeline variants being combined. to represent a combination of n pipeline variants, we randomly sampled n pipelines and created an ensemble classifier as outlined above. this process was repeated with replacement times for each value of n ranging from to student's t-test methods comparison the pool of all individual methods across the signatures was split based on a single aspect of the pipeline. we compared pipelines only differing on a single aspect using the paired t-test to assess statistical differences between pipelines. permutation sampling for variable number of pipelines in the ensemble when subgrouping for methods comparison as part of the method comparison, the pipelines where subgrouped based on a single aspect of the pipeline and then within the subgroups ensembles of a varying number of the pipelines were constructed. to represent a combination of n pipeline variants, we sampled n pipelines and created an ensemble classifier. for each value of n, all possible combinations containing n unique pipeline variants were created. visualization all plotting was performed in the r statistical environment using the lattice, latticeextra, rcolorbrewer and cluster packages.  ensemble classification approach each dataset was pre-processed using different pipeline variants. each biomarker was then applied separately for each pipeline variant, producing an ensemble of predictions for each patient and biomarker. these were analyzed for consistency and combined to form a single ensemble classification. figure  outlines the approach used. we separated our datasets according to the microarray platform used, and tested the two most widely-used platforms at the time of writing according to depositions in the gene expression omnibus: hg-u133a and hg-u plus. since both platforms are affymetrix arrays and therefore have the same set of potential normalization methods, we can perform inter-platform analysis independent of pre-processing. univariate gene analysis we first investigated the univariate performance of individual genes to determine how the prognostic power of these simple biomarkers is influenced by pre-processing differences. as shown previously for lung cancer, the prognostic ability of individual genes varied dramatically across methods. of the genes represented on both array platforms tested, 74% reached statistical significance after multiple-testing correction in at least 1/ pipeline variants. by contrast, only 16% reached significance in at least 12/ pipelines and none were significant in all pipelines. three pipeline variants identified zero genes, while three others found a single gene, which was not identified in the other pipelines. these data clearly indicate that simple union and intersection approaches are inappropriate. interestingly, all six pipelines that resulted in either one or no prognostic genes involved analysis of hg-u133a data, using either the rma or mbei algorithms, along with the “separate” dataset-handling approach. there is an evident difference between the patterns of significant genes on each platform. the lowest concordance between pipelines is shown in the inter-platform correlations. different aspects of the pipeline appear more highly correlated depending on the platform and there is no clear ordering of which aspect is more important without interactions. we were able to use linear-modeling to show that the choice of pre-processing method is strongly deterministic for the number of statistically-significant genes identified. we considered a complete model of all pairwise interactions and main effects, then used the akaike information criterion for backwards stepwise refinement. a model containing the main effects: platform, pre-processing algorithm, data-handling type and their pairwise interactions resulted, indicating that the relationship is deterministic, not stochastic. we note that interactions are critical: a simple model of main-effects was not explanatory. for the linear model, y=w+x+∑i=15zi+w:x+∑i=15w:zi+x:zi where y is the number of genes, w is the platform, x is the data handling and z1…z are specify the options for the pre-processing algorithm, the coefficients that have a p <  are shown. multi-gene signatures we next focused on multi-gene classifiers, seeking to determine if our single-gene results could be generalized. we compared the hazard ratios from cox modeling of the ensemble and the individual classifications for published hypoxia signatures. for all multi-gene signatures, superior classification was defined as the classification with a higher hazard ratio. as seen with the single gene classifiers, variation was observed between classifications from the different pipelines and there was not one single variant which consistently resulted in larger risk stratification than the others. further this analysis identified microarray platform as another possible source for variation. one pipeline variant showed the lowest risk stratification of the pipelines on one platform and the largest of the pipelines on the other platform. as shown in figure  ensemble classification performed better than individual pipelines and improved signature performance for both microarray platforms.analyses for all signatures showed that performance was sensitive to pre-processing choices and, in the majority of cases, the ensemble classification improved prognostic ability over individual pipeline variants. for half of the signatures, ensemble classification resulted in superior risk stratification compared to classifications from the individual pre-processing pipelines. moreover the ensemble technique was almost always superior to the “typical” pre-processing techniques, exceeding the median of the techniques in 24/ signature comparisons. the buffa metagene and the winter metagene showed similar results across pipeline variants, but many of the signatures performed very differently depending on the dataset platform. overall signatures showed better risk-stratification on hg-u plus arrays, although this was signature-specific. some signatures showed consistently better results on the hg-u plus dataset compared to the hg-u133a dataset. conversely, starmans cluster and cluster performed better in the hg-u133a datasets. the buffa and the winter metagene were the only signatures which were statistically significant across all pipelines tested. hu and sorensen, additionally, were other signatures with statistically significant ensemble classifications for both datasets. in contrast, starmans clusters and seigneuric 0% early signatures did not perform well in either dataset; none of their ensemble classifications were statistically significant. in general, if a signature performed poorly for single pipeline variants, using the ensemble classification did not improve it. this was demonstrated by the correlation between the hazard ratios for the ensemble classification and the maximum hazard ratios for classification from the individual pipeline variants. since previous analyses involved comparing unequal numbers of patients classified, we also compared ensemble classification to classification for the individual pre-processing methods. in this way, we match patient numbers between the two conditions, removing this potential confounding variable. in general, this approach yielded fewer statistically significant results, although both the range and the variance of hazard ratios increased for every signature using this classification algorithm. however the comparison between of ensemble classifications and individual classifications shows that patient-number differences are not the origin of the superior performance of ensemble classification. for 13/ signatures, the ensemble classification was superior to all classifications from the individual pre-processing pipelines and in 26/ signatures the ensemble exceeded the median classification. signature comparison to better understand which signatures were more successful, all individual classifications were compared. unsupervised clustering of the percentage agreement of concordant patient classifications between individual pipeline variants for each signature showed that they mainly clustered by signature, rather than by pipeline composition. this indicated that, although pre-processing substantially influenced biomarker performance, the genes in the signature characterized the overall partition and determined whether it was a poor or good biomarker. the buffa metagene had the most consistent patient classifications across pipelines, but hazard ratios still ranged from to. although, we evaluated only hypoxia signatures, patient classifications did not agree across signatures. signatures of ensemble classifications that were statistically significant generally classified a larger fraction of patients. what is the optimal ensemble size? having shown that the ensemble-approach improved classification for most biomarkers and datasets, we explored the limits of its performance. we wondered if distinct pipelines were always necessary, and therefore evaluated the number of pipeline variants required for optimal performance of the ensemble classifier. if creating an ensemble of four pipeline variants is equally successful to one from eight variants, then it is not beneficial to introduce the complexity and computational-costs of pre-processing with four extra pipelines. focusing on signatures with a significant 24-pipeline ensemble, different combinations of pipelines, ranging from combinations of only to all were evaluated. these analyses indicated that in general increasing the number of pipeline variants resulted in an increase in absolute effect size which started to plateau as the number of methods in the ensemble increased. in parallel, the percentage of patients classified with the ensemble method decreased and plateaued. most signatures shared the same shape but with different rates of hazard ratio increase. the sorensen signature on the hg-u133a dataset plateaued at about four pipeline variants. therefore, in this case, randomly choosing four pipeline variants to combine provided roughly the same risk stratification as using all pipelines. conversely, for the winter metagene signature in either dataset, the mean hazard ratio continued to increase all the way up to pipelines, though the curve was steeper at the beginning then in the end. although the hazard ratio stopped increasing in some cases, stability continued to increase as the number of methods in the ensemble increased. this is demonstrated in additional file 8: figure s by the tightening of the hazard ratio range as the number of pipelines is increased. considering the winter metagene signature in hg-u133a data, the ensembles created from nine or more of the pipelines outperformed all single pipeline classifiers. many ensembles did not require all variants to be an improvement over all non-ensemble methods. even if the ensemble of variants was not an improvement over non-ensemble methods, there may still have been an ensemble of a subset of the variants which was superior to the non-ensemble methods. these data provide a compelling rationale to consider and evaluate ensemble pipelines for all microarray-based biomarkers. methods comparison after showing that ensembles are beneficial, we wanted to look at whether we can determine the combination of pipelines that lead to higher hazard ratios in order to add the most benefit for each additional pre-processing pipeline. there is a clear relationship between the number of patients classified in the ensemble and the gain in hazard ratio, meaning that the ensemble is choosing to exclude the right subset of patients. methods that produce less-correlated classifications gain more from the ensemble classification. however, if we look at which methods are diverse by a different metric such as the profiles of prognostic ability of each gene as a single gene classifier, there is only a slight but not obvious increase in hazard ratio from using more diverse pipelines in the ensemble classification. to help direct pipeline choices, we sought to address whether certain aspects of the pipeline resulted in better or worse performance. for each aspect of the pipeline, the hazard ratios were grouped per variant of that aspect and compared. this was done for both platforms separately and combined. on both platforms there was a significant difference between annotations. on hg-u133a, alternative annotation had higher hazard ratios. in direct contrast, hg-u plus performed better with default annotation. by contrast, the optimal pre-processing algorithm was similar in both platforms, with rma and mbei performing better than gcrma and mas. rma and mbei showed similar results as did gcrma and mas. furthermore, we analyzed the effect of changing the number of variants in the ensemble when creating only ensembles from common pipeline variants. once again, variant success is not necessarily consistent across signatures. ensemble of signatures to further filter out unreliable classifications we investigated combining the classifications from two signatures. the buffa metagene and the winter metagene performed best across our analyses. these two signatures share genes. expansion of the ensemble classification to only classify patients that both signatures agreed on improved risk stratification compared to ensemble evaluations of both signatures. to complete the analysis and expand the number of patients classified, we also pooled the unanimous classifications. this failed to improve risk stratification compared to ensemble evaluations of both signatures; however, prognostic performance was improved over all the signatures’ individual pre-processing methods. further, many more patients were classified than with the basic ensemble technique, suggesting that ensembles of signatures could be used to further remove noise or to increase the number of patients given confident molecular classifications. discussion the purpose of pre-processing is to remove “noise” from the data. however, since no method is perfect, each pre-processing pipeline removes a somewhat different aspect of the “noise”. indeed, groups around the world have focused on identifying the “optimal” pre-processing technique for different types of data. the principle of ensemble classification is that by combining pre-processing approaches we can select the parts of the data which are reliable across the multiple approaches. the central tendency of this pool of methods is thus predicted to lie closer to the “true” value, and thereby to provide a better biomarker. although different pre-processing methods may cause some variation in the analysis, pre-processing is expected to have a minor effect on the core experimental results and conclusions. our previous work has indicated this is not the case and pre-processing caused major outcome differences in non-small cell lung cancer. here we systematically extend and deepen these analyses to explore the variation caused by algorithmic diversity in pre-processing. at the single gene level substantial differences in prognostic power were seen in univariate analysis. therefore pre-processing is part of the reason different studies identify different biomarker genes. many authors will use public data to show that a given gene is prognostic; however, essentially all genes can meet that criterion, depending on which platform and pre-processing technique is used. single genes did not appear to behave the same across pipelines demonstrating variation in classification results are expected and signatures are dependent on the pre-processing platform they were discovered on. the fluctuation in prognostic ability for each gene partially explains why we get different results for multi-gene signatures and why there is such difficulty validating biomarkers between research groups. in combination with data showing a plurality of gene sets are associated with outcome in breast cancer and other cancers, the variation suggests that mining of public data for prognostic significance is very prone to over-fitting and multiple-testing concerns. therefore robust, permutation-based approaches need to be developed. in non-small cell lung cancer, starmans et al. showed one example of a permutation-based approach where an ensemble biomarker classifier improved survival separation between low risk and high risk patient groups. here we extended this finding showing that the method replicates on two microarray platforms representing separate datasets in breast cancer for a series of biomarkers. across both platforms, there was not a single pipeline that unfailingly outperformed all other pipelines; therefore, the ensemble classification provides a generalized approach to improve biomarkers, both in terms of performance and stability, without determining the actual optimal pre-processing pipeline. furthermore, in many of the cases, the ensemble classifications outperformed all single pre-processing methods. the ultimate value of the ensemble classifiers as a concept was demonstrated with the buffa metagene and winter metagene. for these two signatures, any ensemble classifier comprising at least nine pipelines on hg-u133a or pipelines on hg-u plus arrays generated superior risk stratifications compared to all the classifiers that used only a single pre-processing pipeline. consequentially, ensemble classifiers are most definitely beneficial and should be used. the ensemble approach did not improve all biomarkers. biomarkers with generally bad risk stratifications across pre-processing pipelines still showed poor performance when combined in an ensemble. the ensemble approach magnified the separation of risk groups rather than corrected for a poor initial biomarker. therefore the ensemble approach can also be used as a metric to assess the quality of biomarkers, distinguishing between poor and good signatures. by the statistical significance and consistency in risk stratification improvement across the datasets, the buffa and winter metagenes are shown to be strong, consistent signatures. by the same metric, seigneuric 0% early and starmans cluster appear to be poor signatures validating previous findings where these signatures did not show prognostic power. a disadvantage of the ensemble classification is that a fraction of patients are not classified. only patients with robust risk classification across pipelines in the ensemble are assigned to risk groups. here, using only the statistically significant ensembles 16% to 80% of the patients were not classified. the signatures showing significance on both platforms tend to have a higher percentage of patients classified than the signature significant on only on platform and the signatures not significant on either platform. nevertheless a patient classified as unreliable with one signature may be robustly classified using a different signature. this was shown by intersecting the buffa and winter metagenes, which resulted in improved prognostic power compared to the single pre-processing pipelines and classified more patients than with the two ensembles individually. one might consider taking this into account in biomarker-developement by attempting to construct ensembles that minimize the fraction of unclassified patients in the training dataset, although unclassified patients could resort to standard clinical care. an important note about our approach to using the ensemble classification is the diversity of the pre-processing methods. our choice of signature evaluation meant that a log2-transform on the data does not create different classifications. for example, the multiple methods mas and log mas pre-processing are actually only one pipeline variant and aren't filtering out additional unreliable patients. here, the array of pre-processing methods in reality only gives unique classifications so there is less diversity than the numbers indicate. however, for other signatures, such as the risk-score and clustering methods evaluated previously, up to unique classifications are possible. in multiple signatures, the ensemble of all pre-processing pipelines is an improvement but not the optimal ensemble classification. so future analysis to refine the process and finesse which pipelines to use and how they should be combined would be advantageous. an important future direction is exploring how ensemble methods can be improved by incorporating greater algorithm diversity. continual addition of diverse methods may increase the optimal ensemble classification or the optimal ensemble may be a certain combination of pipelines and new additions may not lead to an increase. in some cases, such as the hu signature on hg-u plus addition of another separate data handling pipeline is not likely to increase the risk stratification but adding a merged data handling pipeline would be advantageous.  we systematically show that differences in pre-processing create differences when using biomarkers. this effect of pre-processing is important for the research community to recognize and consider, as accurately accounting for it will advance biomarker discovery, validation and ultimately clinical application. we found that the buffa metagene is the most consistent biomarker and therefore most clinical useful signature evaluated and we show that application of ensemble classification technique is beneficial for improving risk stratification both in terms of effect size and stability of biomarkers. competing interests all authors declare that they have no competing interests. authors’ contributions database generation and curation: sh. performed statistical and bioinformatics analyses: nsf, mhws, pcb. data interpretation: nsf, mhws, pcb. wrote the first draft of the manuscript: nsf. initiated the project: pcb. supervised research: mhws, pl, pcb. all authors read and approved the final manuscript. supplementary material additional file 1: table s1 prognostic signature descriptions. click here for file additional file 2: table s2 gene counts per prognostic signature. click here for file additional file 3: figure s1 correlation of gene univariate analysis. analysis of consistency between methods for the prognostic ability of each gene shown in figure  the heatmap shows pairwise comparison of all the pipeline variants where the comparison is spearman's correlation estimate of the fdr-adjusted p-values for univariate cox proportional hazard ratio modeling analysis of genes analyzed on the set of pipelines. click here for file additional file 4: figure s2 platform comparison by signature. comparison of hazard ratios for the series of prognostic signatures on hg-u133a and hg-u plus. hazard ratios were derived from cox proportional hazard ratio modeling. each triangle represents the ensemble classifier's hazard ratio and the circles represent the individual pipeline variants. the 95% confidence interval is shown for each ensemble. for the individual pipeline variants, the 95% confidence intervals are shown in additional file 5: table s click here for file additional file 5: table s3 hazard ratio 95% confidence intervals for classifications on the individual pipeline variants. click here for file additional file 6: figure s3 risk stratification across classification pipelines and prognostic signatures with equal number of patients classified. comparison of hazard ratios and corresponding p-values from cox proportional hazard ratio modeling between ensemble classifications and individual classifications on a subset of patients with the highest and lowest signature scores on hg-u133a platform, hg-u plus platform. the hazard ratio is represented by the size and colour of the dot and the background shade represents the p-value. click here for file additional file 7: figure s4 signature comparison. analysis of consistency across both significant prognostic signatures and signatures that were not. heatmaps are shown for the pair-wise comparison of all the single pipeline classifications for the individual pre-processing methods and the ensemble scores derived from these individual classifications per patient for each signature. in b, the signatures are ordered by the number of patients classified unanimously across all the pipeline variants. from left to right, the number of patients classified in the ensemble for each signature decreases. click here for file additional file 8: figure s5 ensemble hazard ratio range. the range of hazard ratios for ensembles from different number of pipeline variants. the horizontal pink dashed line shows the highest hazard ratio of the individual methods; all the ensembles above the line are improvements on current pre-processing practice. the x-axis indicate the number of pipeline variants combined to create ensembles. the grey background shows the numbers of pipeline variants where all the ensembles created are superior to every single individual method. the hazard ratio, p-value and number of patients classified for each ensemble shown is provided in additional file 11: table s and additional file 12: table s click here for file additional file 9: table s the hazard ratios, p-values and number of patients classified for all of the classifications on hg-u133a in additional file 8: figure s to make the data easier to use, each signature is in a separate a tab delimited table/text file and files for each platform are packaged and compressed separately. each row in the tables is a patient classification and there are repeated rows since we were sampling the pipelines with replacement. the first columns of each table is whether the pipeline specified in the column name is used in the ensemble classification with meaning the pipeline is in the ensemble and meaning it is not. columns – are the hazard ratio, p-value and the number of patients classified for the classification respectively. click here for file additional file 10: table s the hazard ratios, p-values and number of patients classified for all of the classifications on hg-u plus in additional file 8: figure s to make the data easier to use, each signature is in a separate a tab delimited table/text file and files for each platform are packaged and compressed separately. each row in the tables is a patient classification and there are repeated rows since we were sampling the pipelines with replacement. the first columns of each table is whether the pipeline specified in the column name is used in the ensemble classification with meaning the pipeline is in the ensemble and meaning it is not. columns – are the hazard ratio, p-value and the number of patients classified for the classification respectively. click here for file additional file 11: figure s method correlation effect on hazard ratio. comparison of the effect of method diversity in ensembles of on the increase in hazard ratio from the maximum of the individual classifications for winter metagene classifications on hg-u133a and hg-u plus. part a measures how correlated methods are by their percent agreement between methods which is also equivalent to the number of patients classified. part b measures the relatedness of the methods by the spearman's correlation of how prognostic each gene is for a method. click here for file additional file 12: figure s combining signatures. prognostic ability of combining the ensemble approach for the winter metagene and the buffa metagene was evaluated with kaplan-meier survival analyses. hazard ratios and p-values are from cox proportional hazard ratio modeling. the intersect is using only patients that are in agreement between winter metagene and buffa metagene. the union is pooling the patients from winter metagene and buffa metagene. click here for file acknowledgements the authors thank nathalie moon for advice on statistical analysis and all members of the boutros lab for helpful suggestions. funding sources this study was conducted with the support of the ontario institute for cancer research to pcb through funding provided by the government of ontario. pcb was supported by a terry fox research institute new investigator award. this work was supported by the center for translational molecular medicine, eu imi program, eu 7th framework program and dutch cancer society through funding to mhws and pl.
the reproducibility of transcriptomic biomarkers across datasets remains poor, limiting clinical application. while further research is required, this approach has the potential to become a standard for transcriptomic biomarkers. importantly, ensemble classification improves biomarkers with initially good results but does not result in spuriously improved performance for poor biomarkers. to test this hypothesis, we systematically assessed the effects of pre-processing on biomarker classification using different pre-processing methods and distinct signatures of tumour hypoxia in datasets. importantly, exploiting different pre-processing techniques in an ensemble technique improved classification for a majority of signatures.
 the reproducibility of transcriptomic biomarkers across datasets remains poor, limiting clinical application. we and others have suggested that this is in-part caused by differential error-structure between datasets, and their incomplete removal by pre-processing algorithms. methods to test this hypothesis, we systematically assessed the effects of pre-processing on biomarker classification using different pre-processing methods and distinct signatures of tumour hypoxia in datasets.  we confirm strong pre-processing effects for all datasets and signatures, and find that these differ between microarray versions. importantly, exploiting different pre-processing techniques in an ensemble technique improved classification for a majority of signatures.  assessing biomarkers using an ensemble of pre-processing techniques shows clear value across multiple diseases, datasets and biomarkers. importantly, ensemble classification improves biomarkers with initially good results but does not result in spuriously improved performance for poor biomarkers. while further research is required, this approach has the potential to become a standard for transcriptomic biomarkers. 