BACKGROUND
detecting genes-by-samples patterns in expression data
after removing genes that exhibit little variance, the standard script for exploring microarray data applies two-way hierarchical clustering , followed by a visual search for patterns displayed in a red-green heatmap  <cit> .
for breast cancer in particular, this procedure has proven immensely productive.
it can be credited with the discovery <cit>    <cit>  of the basal subtype, and, more broadly the identification of subtypes of breast cancer that hold out the potential to inform clinical practice  <cit>
despite its utility, the standard script suffers from several limitations, in particular the instability of the binary tree of the clusters found  <cit> .
perturbation and re-sampling techniques are available to gauge the robustness of the clusters defined by subtrees  <cit> .
but, small changes in the selection of genes or choice of samples can result in disconcertingly large changes in the overall configuration of the tree, which calls into question any typology defined on such tree-based partitions  <cit> .
a different problem stems from the disproportionate impact of large, tightly coordinated clusters on the overall arrangement of the tree <cit> .
because the largest gene sets, for example estrogen or immune response, will dominate the branching of the tree, smaller sets may be broken up and redistributed.
the problem of large dense sets of genes occluding smaller sets arises in the simple or one-way application of hc; it is compounded when two trees are visually crossed in the two-way hc used in the standard script.
the question then becomes one of what can be faithfully represented in the two-dimensional arrangement.
the brief answer, as spelled out by hartigan in the context of "direct" clustering, is that clusters jointly defined by two trees can be rendered as contiguous regions only if they are either disjoint or nested  <cit> .
an alternative to two-way hierarchical clustering, biclustering, seeks to find submatrices in the array of expression values that satisfy some defining criteria.
essential to biclustering is the notion that, given a matrix of expression values, the pattern of coordinated expression for a group of genes may be confined to only a subset of samples.
that is, the pattern is "local", and may not be detectable using a "global" measure of the similarity between pairs of genes, e.g., pearson correlation coefficient computed on vectors of expression values  <cit> .
motivated by this concern with detecting "local" structure, several algorithms for detecting biclusters have been advanced.
for surveys, see  <cit> .
because of the computational complexity inherent in the biclustering problem, which is provably np-hard under several formal descriptions <cit> , these programs abandon exhaustive search and either resort to heuristics or impose bounds on the size of the genes-by-samples submatrices that can be recovered.
the method applied in this paper reframes the biclustering problem such that an exhaustive enumeration of all of the submatrices that meet a certain definition can be achieved.
included in the resulting inventory are not only such familiar sets as "estrogen" and "erbb2/17q12", but additional sets highly prognostic of survival.
some of these newly defined sets may serve alongside "estrogen" etc., as building blocks for future models and classifiers of breast cancer.
methods
defining gene patterns by partitions on samples
hc is an agglomerative method that recursively joins sets of genes or, separately, sets of samples  <cit> .
the starting point is the notion that similar genes are assigned to the same set  using a measure of similarity such as correlation <cit> .
a wholly different approach can be predicated on the notion that two genes belong to the same set provided that they induce the same partition on the set of samples.
as an example, the erbb2/17q <dig> gene set in figure  <dig> appears to satisfy such a requirement.
here the samples have been ordered by column sum.
it is apparent that the largest  values of each of these six genes align to a considerable degree.
that is, they indicate, or select for, much the same subset of columns .
while the notion of genes inducing partitions like this may be intuitive, rendering the notion such that it can be implemented in a search algorithm requires further specifying what constitutes an admissible partition.
crucially important is the partition's size and the degree to which the samples induced by any two genes must align for the partition to be considered the "same".
these two parameters  control a matching rule that can serve to explicitly define what is meant by "gene set".
setting up the enumeration
the two parameters that define what counts as the "same partition" suggest a partition-based strategy for finding all sets of genes  in a matrix of expression values, namely, try all combinations of size and tolerance.
because this collects all the genes-by-samples patterns, it constitutes an enumeration.
since the size and quality of the partitions to be discovered are not known in advance , to enumerate all of the partitions in a dataset, and the genes that support them, requires trying all combinations of size and tolerance.
this sets up a two dimensional grid search with partition size  and tolerance  as dimensions.
given the number of samples , and the number of draws, that is genes, bounds can be set for the parameters  and  by computing p-values with a counting rule .
though the rule for bounding t works well enough with smaller datasets, we find that with data sets of the size of the nki <dig> <cit> , these p-values are so impossibly small that they provide little guidance in controlling t. an alternative strategy relies on a heuristic suggested by results of an initial spiral search in which s and t are stepped in large increments until gene sets are found, at which point s and t, are incremented, separately, and in small steps.
this preliminary scan of the four datasets finds:
1) for most partition sizes s, small  values of t fail to detect any gene sets at all, while large t settings yield a single set composed of most, if not all of the genes.
between these two extremes, for any given s, there is a relatively narrow range of t that returns multiple distinct gene sets.
2) in size and composition, these sets vary smoothly as s and t are stepped.
3) the number of unique components  is surprisingly small, for example, between thirty and forty in the nki <dig> data.
in light of these observations from the spiral searches, rather than bounding t by a combinatorial counting rule, we choose to set t such that specific gene sets remain distinct.
in the enumerations reported in this paper, the set we wish to preserve consists of interferon genes, a clinically significant set that the initial spiral search encountered repeatedly.
the resulting "stop rule" says, in effect, that for each partition size, which ranges from a small constant , to n minus that same constant, increment t, starting at  <dig>  until the interferon set disappears as a distinct object as it merges with other immune sets.
the motivation for devising this stop rule will be developed further in the context of the gene sets that decompose "immune response".
as the algorithm steps through its search it typically finds multiple versions of the same set, which may differ only marginally in the number and composition of the genes.
this means that rather than a single list of genes, a "gene set" in fact refers to a set of sets, like those listed for the erbb <dig> set in table  <dig>  each column in that table catalogues the genes found for the erbb2/17q <dig> gene set in the uppsala data as detected at ten combinations of the "size" and "tolerance" parameters.
the list of genes in each column is like a snapshot of the same object in gene space, viewed at a slightly different resolution.
the partition size  and tolerance  are listed at the head of each column.
the first instance of this gene set is detected at s =  <dig>  t =  <dig>  as s and t are stepped, the number of erbb2/17q <dig> genes increases from size to a maximum of  <dig>  at s =  <dig>  t =  <dig>  from that point the number of genes decreases.
it is important to note that the end product of the enumeration is both the list of genes in each set and the partition induced on the samples by that set.
while the number of genes in the sets can vary from six  to as many as  <dig> , regardless of the number of genes in the set there can be only one partitioning pattern per gene set.
it is in this context that the issue of genes with multiple probe sets can be addressed.
as would be expected, multiple probe sets spotted for the same gene will, to the extent that they induce the same partition on the samples, be assigned to the same gene set by the algorithm.
a case in point is the small stromal gene set in figure  <dig> which consists principally of decorin  and fibulin  <dig> .
decorin is spotted four times on the affymetrix hg-u133a and the algorithm finds that all four decorin probe sets induce the same partition pattern on the  <dig> samples in the miller dataset  <cit> .
that same pattern is also induced by two copies of fbln <dig>  as well as by single instances of glt8d <dig>  ctsk, prrx <dig>  and spon <dig>  and several ests.
all in all, in the miller data, this partitioning pattern is found at five combinations of s and t settings.
for the purpose of making the survival analysis manageable, we summarize  these five versions of this gene set into one composed of those probe sets that appear in at least half of the original versions.
the result is the set pictured in the heatmap in the figure, to which metastasis events have been attached.
as described more fully in the results section, the χ <dig> for a partition at the median value of a vector of column sums is  <dig>  .
the χ <dig> for first versus last quartile is  <dig>  .
since each of the genes in a gene set realize the same partitioning pattern, the end result, in terms of partitioning and survival analysis would be virtually the same if the four decorin probe sets were replaced by any one .
this issue of multiple probe sets for the same gene does become a problem for the detection of the smallest gene sets in the case in which the partitioning pattern is realized exclusively in probe sets all of which spot the same gene.
this occurs occasionally in the affymetrix data, and we consider these instances to be artefacts of the chip design.
detecting patterns that were heretofor undetectable
the virtues of hc and the standard script include the fact that it is intuitive, simple to implement, and effective .
moreover, the software that implements this method of discovery is free and well maintained  <cit> , and the output is not only immanently useful, but a thing of beauty.
replacing this standard script requires both a conceptual shift and considerable computational apparatus.
for these reasons, such a switch can only be justified if the alternative not only reproduces the standard results, but significantly extends those results.
that is, it must find additional patterns that are biologically meaningful and clinically significant.
to show that this is in fact the case we enumerate the gene sets in four well-studied breast cancer microarray data sets listed in table  <dig>  we refer to these as: uppsala, stockholm, transbig, and nki data sets.
three of these use affymetrix u133a chips, a platform that has attained the status of an industry standard, and has been used perhaps in more studies than any other microarray platform.
to identify gene sets conserved across platforms, we add the fourth data set, generated by the nederlands kanker instituut , which used a custom agilent chip.
the uppsala and stockholm cohorts  <cit>  are population studies .
the transbig dataset is a validation series for a gene-based classifier  <cit>  and, accordingly, consists of younger  node-negative patients  <cit> .
the nki data was collected in conjunction with a second classifier  <cit> , and consists of younger patients, both node-positive and node-negative.
for purpose of survival analysis in this paper, the end point for the uppsala data is breast cancer specific survival ; the end points in the other three data sets are distant metastasis free survival .
the datasets were selected, in the first instance, because of the quality of the data, both expression and accompanying clinical data.
but they were also chosen to set up explicit comparisons between the gene sets detected by the enumerations in this paper and gene clusters and signatures induced on the same data by other methods and programs.
data sets were downloaded from the geo website  <cit> , or, in the case of the nki data, from the stanford website  <cit>  that accompanies  <cit> .
missing values were estimated using k nearest neighbours,   <cit> .
because the algorithm takes the full expression matrix as input, e.g., in the case of the uppsala data, a matrix of size  <dig>  ×  <dig>  there is no prior filtering/selection of genes.
RESULTS
in each of four breast cancer datasets in table  <dig>  the rank-based bi-clustering algorithm finds between thirty and forty sets of genes.
many of these sets are conserved across the four studies despite differences in platform: affymetrix u133a , and custom agilent .
several of the largest and most prominent of the sets closely resemble gene clusters reported many times in microarray studies of breast cancer.
these include "estrogen", "proliferation", erbb <dig>  and "adipose"  <cit> .
other sets, a number of which are relatively small , appear to be novel.
these include "histone", "hemoglobin" "amplicon 8p11-12".
overall, nearly a third of the gene sets detected are clinically significant as measured by association with survival, and some are very significant with a log-rank χ <dig> values exceeding  <dig> .
for purpose of presentation, in the course of enumerating the gene sets in each of the four data sets, we proceed in three steps.
the gene sets detected in each of the three affymetrix data sets are first tabulated separately.
this is followed by a comparison of these sets to identify gene sets conserved on the affymetrix u133a platform.
this, then, is followed by an enumeration of the gene sets in the nki <dig> data set which uses a custom agilent chip.
these sets are then compared to the affymetrix sets for the purpose of identifying gene sets conserved across the two platforms.
regarding the clinical significance of the gene sets, since a gene set, in general, is detected multiple times , the survival analysis of a gene set would yield multiple log-rank values.
to simplify the analysis we have replaced these multiple versions of each set with a "core" set comprised of those genes that appear in at least half of the instances detected.
five log-rank χ <dig> values are reported for each core gene set.
the first reflects all the samples in a data set partitioned at the median value of the column sums of the genes in the set.
the second partitions all samples by lowest and highest quartiles.
the third restricts the samples to those which are er positive, and partitions at the median.
the fourth also restricts the samples to er positive, partitioning by quartiles.
the fifth restricts the samples to er negative, partitioned at the median.
because of small class size, the er negative samples are not partitioned by first and last quartile.
dataset 1: uppsala cohort
the proliferation set is strongly prognostic of recurrence in the uppsala data, as it is in the other three data sets, which accords with a number of microarray-based studies of breast cancer  <cit> .
other sets strongly prognostic of recurrence are the histones, the metallothioneins/16q <dig>  as well as gapdh, cd <dig>  affx-m27830_ <dig>  and gnas.
the histone set includes hist1h2bf hist1h2be hist1h2bh h2bfs hist1h2bk hist1h2bd, while the metallothioneins gene set consists of seven, or perhaps eight, isoforms and two ests: mt1e, mt1f, mt1g, mt <dig> h, mt1m, mt1x, mt2a, loc <dig>  up-regulation of mt2a, the most abundant of these, is associated with more aggressive breast cancer and poor prognosis  <cit> , and it is reported that the over-expression of metallothioneins predicts resistance to doxorubicin  <cit> .
gnas is located on the long arm of chromosome  <dig> at 20q <dig> , lying just outside the interval of the 20q <dig> amplicon investigated by ginestier et al.
<cit> , but falling within one of the two 20q <dig> amplicons analyzed by yao et al  <cit> .
it is reported that increased copy numbers for 20q <dig> amplicon genes occur in 12% of primary breast tumors, and are associated with more aggressive disease  <cit> .
gapdh is often used as a housekeeping gene  <cit> .
for example it is one of the five genes used to normalize recurrence score in oncogenedx  <cit> .
but, gapdh is reported to be up-regulated in some cancers, e.g., by a factor of  <dig> to  <dig> in non small cell lung cancer compared to normal lung tissue  <cit> .
in a study of gapdh expression in breast cancer, revillion et al., <cit>  conclude with the warning that it should not be used as a control rna.
valenti et al reach a similar conclusion  <cit> .
the fact that the gapdh partitions the uppsala samples into good and poor prognostic groups is added evidence that it is not only unsuitable for scaling expression values, but may in fact be a candidate oncogene.
data set 2: stockholm cohort
similar to the results for the uppsala data, for the stockholm cohort three stromal gene sets are strongly associated with increased survival as reported in table  <dig>  also, again, the adipose set proves significant.
in addition, among the sets associated with increased survival is a hemoglobin set consisting of: hbh <dig>  hbb, hba <dig>  hbb, hba <dig>  hba <dig>  hba <dig>  hbb, hba <dig>  hbg <dig> .
a complete list of the genes in all  <dig> gene sets for the stockholm cohort is available in additional file  <dig>  though er and pr status are reported in the aggregate in the original article <cit> , er status was not available for this data set.
consequently the log-rank values tabulated in table  <dig> represent all  <dig> samples partitioned at the median, and by first and last quartiles.
among the sets significantly associated with decreased survival, as in the uppsala data, proliferation is most prominent.
also, again, the histone set, cd <dig>  and gapdh figure among the sets associated with recurrence.
in addition, among the gene sets that have a negative impact on survival in this data set include actg <dig>  a ribosomal set, and a set we have labelled ezrin, because it includes the gene vil <dig>  ezrin expression is associated with metastasis in a number of cancers  <cit> .
in a murine cell line elliot et al show that ezrin is required for metastasis  <cit> , and li et al demonstrated that ezrin silencing reverses cell migration and invasion in a metastatic breast cancer cell line  <cit> .
with regard to survival in breast cancer, bruce et al report that ezrin expression is associated with poor outcome  <cit> .
the ezrin gene set also contains arf <dig>  which is reported to modulate migration and proliferation in breast cancer cell lines via the regulation of the pi3k pathway  <cit> .
data set 3: transbig
while the transbig data employs the same affymetrix platform as the uppsala and stockholm datasets, it differs in terms of sample selection.
the first two data sets are population-based, while this data set, as a validation series, reflects the sampling criteria used in the development of wang et al.
's 76-gene classifier  <cit> .
as a consequence, patients were younger , node-negative, with smaller tumour grade .
despite these differences, as is apparent in table  <dig>  the gene sets detected are largely the same as those found in the uppsula and stockholm cohorts.
a complete list of the genes in each of the  <dig> sets discovered in the transbig data are available in additional file  <dig>  as in the previous two data sets, stromal is associated with increased survival.
distinctive of this third data set are the three ribosomal sets associated with survival, a relationship which for two of the sets  and ribosomal) is particularly strong for er positive samples.
with regard to gene sets prognostic of decreased survival, again, as with the uppsala and stockholm cohorts, histone and proliferation figure prominently, as do the metallothioneins.
also among the sets associated with decreased survival are the basal set and nktr.
comparing gene sets across the three affymetrix u133a data sets
factors that impact the ability of the algorithm to detect sets that induce a common partition include the design of the microarray platform.
some genes, for example cd <dig> or lst <dig>  are rendered by multiple probe sets.
to the extent that they take similar expression values for the same samples, the algorithm will identify them as a  gene set in their own right.
that is, all gene sets that satisfy the matching criteria are reported in the enumeration, including several as small as six genes in size.
while some of these will prove to be artefacts of chip design, others may prove to be significant, for example the eight-gene 8p11- <dig> gene set detected in the nki data.
a virtue of the method over hc is this capacity to find even the smallest sets that induce important patterns on the samples.
data set 4: nki 2002
of the  <dig> gene sets discovered in the enumeration of the nki <dig>  available in additional file  <dig>  eight are significantly associated with increased survival.
of these, the two largest  are the estrogen set, well-documented in microarray-based studies of breast cancer  <cit> , and a second set that we have labelled "foxa1".
these sets are disjoint except for the pivotal foxa <dig>  gene which appears in both, and which is represented by two probes on the custom agilent chip used in the nki data.
while the estrogen and foxa <dig> gene sets are both significantly associated with increased survival, the two sets induce substantially different partitions on the  <dig> tumor samples, possibly suggesting different positive mechanisms at work.
the gene foxa <dig> is known to correlate strongly with estrogen receptor alpha  <cit> , and the mechanism that accounts for this association has been established, namely, er binding requires foxa <dig> binding in close proximity  <cit> .
but, it has also been shown that foxa <dig> expression is largely independent of estrogen receptor status.
a substantial proportion of er-negative tumors express foxa <dig>  <cit> , while, in one study, more than 40% of er-positive tumors were down-regulated for foxa <dig>  <cit> .
additional support for the notion that foxa <dig> is independent of er, is the fact that most of the known foxa1-response genes are not erα-response genes  <cit> .
so, paradoxically, it would appear that foxa <dig> is both correlated and uncorrelated with er.
the estrogen and foxa <dig> sets detected in the enumeration, support both propositions.
it should be noted that at least two other genes in the foxa <dig> set, tloc <dig> and sdccag <dig>  are reported to act as tumor suppressors in their own right  <cit> .
two stromal sets, described in the next section, are shown in table  <dig> to be associated with increased survival.
the genes in the 8p11- <dig> gene set essentially tile a region of 8p11- <dig> amplicon between  <dig> m and  <dig>  m. the set largely agrees with the amplicon as refined by haverty et al.
<cit> .
seven of the eight genes map to one of the four 8p11- <dig> sub-amplicons that gelsi-boyer et al identified using array-cgh and reyal et al.
's method for correlating genes within 20-gene windows  <cit> .
studies have found that the expression of each of these genes is significantly correlated with copy number  <cit> .
it is interesting to observe that immune sets are associated with survival both positively t-cell and immuneigg) and negatively mhc-i and immuneinterferon).
as with the three affymetrix data sets, the histone and proliferation gene sets are strongly prognostic of poor outcome in the nki data.
two histone sets, h2b_histone and h3_histone are detected in this data set, but it is only the h2b, which corresponds to the histone set in the affymetrix data, that is significant.
gene sets conserved across the agilent and affymetrix platforms
discussion
while the tables in the results section report the gene sets found and their clinical significance, here the purpose is to compare gene sets, in particular those related to immune response and stromal signalling, to gene clusters and signatures identified across a number of microarray studies.
for many of the sets detected the correspondence is immediate.
this is the case in particular for "estrogen", "erbb2", "basal", and "proliferation", all of which have become virtual fixtures in genomic studies of breast cancer  <cit> .
other sets, including several small sets found only in the affymetrix data, may be trivial in the sense that they consist of housekeeping genes, or reflect the design decision to spot single genes multiple times.
in these cases, the resulting patterns, as detected by the algorithm, are real, but may be of only technical interest in so far as they concern issues of normalization and quality control.
of the remaining sets detected by the enumerations, several are easily described and readily interpreted in biological terms.
examples include:
• hemoglobin: hba <dig>  hbb, hba <dig>  hbg <dig> ..
• histone: hist1h2bf, hist1h2be, hist1h2bfs, h2bfs, hist1h2bk, ...
• metallotheioneins: mt1g, mt <dig> h, mt1x, mt2a, mt1f, mt1e, mt1m
of special interest are the ten sets labelled "immune", and the six or more identified as "stromal".
a natural factoring of immune response
eight immune sets are detected in the uppsala data, five in the stockholm, seven in transbig, and six in nki.
all told, across the four data sets, ten unique immune sets are induced.
of these, two are questionable: one consisting essentially of five copies of casp <dig>  and the other, a small set in the nki containing tnfrsf <dig>  which is most likely a subset of the gene set labelled "immuneigg".
setting these aside, the broad category of "immune-related genes" appears to factor naturally into eight distinct subcategories.
to marshal evidence for, and against, this division, and to aid in their biological interpretation, these eight sets can be matched against clusters and lists of immune genes identified by others using standard methods, generally hierarchical clustering.
the correspondence is one-one between the immune gene sets and eight of loi's et al's "pclusts"  <cit>  and seven of rody et al.
's immune-related metagenes  <cit> .
examining this correspondence in greater detail, seven of the eight "immune" sets correspond almost exactly to the gene clusters used by rody et al.
to compute immune "metagenes" in their study of lymphocytic infiltrates in breast cancer  <cit> .
in that study of twelve affymetrix breast cancer data sets, hierarchical clustering repeatedly yielded a cluster of approximately  <dig> immune-related genes.
applying hierarchical clustering a second time to this cluster revealed seven distinct sub-clusters.
the concordance between the genes in these seven clusters and seven of the immune sets as detected in the enumeration of the uppsala data is available in additional file  <dig>  the one minor discrepancy between the enumerated sets and these seven clusters involves the five lst <dig> probe sets in the monocyte cluster .
in the enumerations of the uppsala and stockholm cohorts, lst <dig> forms its own small gene set.
the survival analysis on the two largest datasets, uppsala and nki, confirm rody et al.
's principal finding  <cit> , namely t-cell genes are prognostic of increased survival time for estrogen-negative patients, with χ <dig> =  <dig> , p =  <dig>  for the miller data, and χ <dig> =  <dig> , p =  <dig>  for the van de vijver data.
but on these same datasets, contrary to rody et al.
who find no association between b-cell/immunoglobulin genes and survival, immuneigg is a significant predictor of increased survival in the nki data, χ <dig> =  <dig> , p =  <dig> , and marginally so in the uppsala data, χ <dig> =  <dig> , p =  <dig> .
the genes in immuneigg closely match the b-cell gene cluster that schmidt et al.
find to be prognostic of increased survival among highly proliferating tumor samples <cit> .
contrary to our results, and to that of rody et al., for tumor samples stratified by proliferation, schmidt et al.
fail to find any association between survival and the expression of their t-cell gene cluster.
rody et al.
attribute this discrepancy to possible differences in cohorts and/or treatments, but the difference in results may also be due to the composition of the respective t-cell gene clusters or sets.
schmidt et al.
's t-cell cluster contains our immunet-cell genes as a proper subset, but it also contains genes that belong to immunestat <dig> and immunecomplement, neither of which we find to be significantly associated with survival.
hence, the significance of the t-cell genes as a predictor of survival may be attenuated by genes from other sets.
immediately relevant to the question of t-cell genes versus igg genes as predicators of survival, calabro et al., assembled eighteen genes from the literature to measure the presence of lymphocytic infiltrate  <cit> .
they found that this list is somewhat associated with diminished survival time for estrogen-positive samples, but is strongly prognostic of increased survival for estrogen-negative samples.
the genes in this list range from ccl <dig>  cd <dig>  cde <dig> ..to ighg <dig> and igj, which effectively merges our immunet-cell and immuneigg sets.
therefore, the separate positive associations with survival between immunet-cell and immuneigg for estrogen-negative samples in the uppsala and nki datasets appear to confirm calibro et al.
's results.
the positive effect of the genes in immunet-cell on survival may in part explain the performance of finak et al.
's stromal-derived classifier  <cit> .
in that study, hierarchical clustering applied to individually matched tumor and normal stroma yielded three clusters of samples that are starkly distinguished by outcome.
from the genes that best discriminate between pairs of these clusters, finak et al.
construct a 23-gene classifier.
ten of these  <dig> are associated with the good outcome cluster of samples.
of these ten, eight belong to the immunet-cell set: gzma, cd8a, cd <dig>  cd <dig>  cd <dig>  plek, runx <dig>  and gimap <dig>  this suggests that the stromal-derived classifier, in assigning samples to the good, poor, and mixed groups, may be powered substantially by the association between t-cell genes and increased survival.
if this is the case, then the stromal-derived classifier provides yet more evidence for the association between t-cell genes and survival for estrogen-negative breast cancer observed in the uppsala and nki data.
the capacity of teschendorff et al.
's seven gene immune response module to identify good outcome samples from among estrogen-negative tumors may constitute yet more support for this association.
those seven genes span at least four of our immune sets with immunedt-cell represented by ly <dig>  immuneigg by iglc <dig> and tnfrsf <dig>  immunemhc-i by hla-f, and immunecomplement by c1qa.
it is interesting to note that in the attempt to validate this seven-gene module on independent data, only four of the seven genes prove significant, and of these, three belong to either immunet-cell , or to immuneigg .
if these are the genes that are driving the performance of the teschendorff et al.
immune response module, then the effectiveness of that signature for predicting increased survival among estrogen-negative samples is consistent with the immunet-cell and immuneigg survival results for estrogen-negative samples in the uppsala and nki data.
overall, the strong association between the t-cell gene set and survival, and the milder association between the immunogloblulin/b-cell gene set and survival, as identified in the enumerations, appears to converge with the important results of each of these several studies despite large differences in approach and research design.
immune/interferon
the immune/interferon set closely resembles the cluster of "interferon response" genes identified by buess et al.
<cit> .
in that experiment fibroblasts were co-cultured with several breast cancer cell lines to investigate cell-cell signalling between stroma and malignant epithelial cells.
in comparing the gene expression of co-cultures to matched monocultures, the starkest difference involved interferon-related genes which were preferentially expressed in co-cultures of fibroblasts with estrogen negative cell lines.
of the twelve genes in the immune gene set, as realized for example in the nki data, five belong to this "interferon response": oas <dig>  oas <dig>  mx <dig>  mx <dig>  and ifit <dig>  the remaining seven include: ifit <dig>  isg <dig>  os <dig>  mtap4a, usp <dig>  g1p <dig>  and gs <dig>  buess et al show that these interferon genes are significantly associated with shorter survival time in the nki data.
the enumeration of the nki <dig> confirms this finding.
the stromal gene sets factor stromal signaling
efforts to delineate the interaction between stromal cells and epithelial tumor are challenged by the complexity of the microenvironment, which has been defined as all components of the mammary gland other than luminal and/or tumor epithelial cells  <cit> .
the effective division of the stroma-related genes into six gene sets by the enumerations may be relevant to this problem.
unlike the immune sets which appear to be distinct, monolithic entitites, the stromal sets resemble a constellation with a core entity, which we designate "stromal", accompanied by several satellite sets, "stromal, stromal, stromal and stromal.
each of these is detected as a self-contained and independent set under at least one combination of partition "size" and "tolerance" parameters, but, as the algorithm increases the partition size, and relaxes the stringency of what qualifies as a match, these sets tend to quickly merge into an omnibus stromal set.
despite this, the importance of distinguishing the smaller sets becomes apparent in the survival analysis: five of the six stromal gene sets are associated with increased survival, several exceptionally so with log-rank χ <dig> values in excess of  <dig>
stromal = core dtf stromal signature
the stromal set, induced in each of the four enumerations, is composed of many, if not all of the collagen and ecm remodelling enzymes featured in west et al's desmoids-type fibromatosis stromal signature   <cit> .
prominent genes include: sparc, cspg <dig>  fbln <dig>  fbn <dig>  and type-i, type-iii, and type-vi collagen genes.
that signature was devised as a proof of concept for a larger, on-going program that exploits the mono-cellular property of soft tissue tumors to inductively define subtypes  of fibroblastic stroma cells  <cit> .
the original dtf signature, comprised of genes differentially expressed between two types of soft tissue tumors, was refined for the purpose of identifying a distinctive stromal response in breast cancer.
in five datasets, including three of the four used in this paper, the dtf stromal response identifies a subset of breast cancer patients who experience increased survival.
the two versions of the dtf signature contain  <dig> and  <dig> genes, respectively.
in an analysis of the functional relations among the proteins that correspond to these genes, this list is reduced further to a protein-protein network of  <dig> genes  <cit> .
the close relationship between the dtf signature and the stromal gene set as realized in each of the four data sets can be conveyed with respect to these twenty essential genes .
stromal/col11a1
among the six stromal sets, stromal/col11a <dig> is the only one that is negatively related to survival, though this is apparent only when inspecting the entire family of sets detected.
in the enumerations, this gene set, consisting exclusively of col11a <dig> and fn <dig> probe sets, is almost always subsumed in a larger stromal set composed essentially of stromal/dtf genes.
nevertheless it appears to be a distinct entity under some of the parameter combinations that control the stringency of the match.
it is interesting to observe that the survival value of the stromal set, as measured by log-rank, and which is generally positive, abruptly goes to zero as coll11a <dig> merges with that set.
this scenario is played out in all four datasets.
the implication is that important relationships between stromal expression and outcome can be lost unless the stromal genes are decomposed into their constituent sets.
stromal/lama2
the stromal gene set is detected in the three affymetrix data sets and is characterized by lama <dig> and col14a <dig>  the genes in this gene set tend to merge with the genes in stromal/darc as the size of the partition is increased and the stringency of the match is relaxed.
stromal/lama <dig> is significantly associated with increased survival: χ <dig>  <dig>   for the uppsala <dig>  χ <dig>  <dig> , p =  <dig>  for stockholm <dig>  and χ <dig>  <dig> , p =  <dig>  for transbig <dig>
stromal/darc
the stromal/darc gene set is detected in the uppsala and stockholm cohorts and is positively associated with survival in both: χ <dig>  <dig>   for the uppsala <dig>  χ <dig>  <dig>   for the stockholm <dig>  genes in this set include: darc, tnxb, ccl <dig>  ldb <dig>  lhfp, and c <dig>  duffy antigen receptor for chemokines  is believed to suppress tumor metastasis through two mechanisms: by sequestering angiogenic chemokines  <cit>  and by inducing kai1/cd <dig> tetraspandin in tumor cells, thereby causing senescence  <cit> .
the merozoite of p. vivax malaria enters red blood cells via darc, consequently individuals who lack darc are resistant to that strain of malaria <cit> .
the association between the lack of expression of darc in a large percentage of black men and increased rates of aggressive prostate cancer, compared to whites, has been established  <cit> .
a similar relationship appears to hold between black women and the increased incidence of more aggressive breast cancer.
it is estimated that 70% of blacks of west african descent lack the expression of darc, which is a population that suffers high rates of both prostate and breast cancer  <cit> .
a second gene in the stromal gene set, tenascin b , is down-regulated during the tumor progression of neurofibromatosis  <cit> .
stromal/adamts5
among the smallest of the stromal sets, stromal is comprised of adamts <dig>  znf <dig> , and four agilent probes that lack gene names.
it is significantly associated with increased survival in the nki <dig> data: χ <dig>  <dig> , p =  <dig> , when partitioned at the median, χ <dig>  <dig> , p =  <dig> , when partitioned by first and last quartiles.
stromal/decorin
stromal consists almost exclusively of dcn and fbln <dig> probe sets.
like the other small stromal sets, with the exception of stromal/adamts <dig>  this set tends to merge with stromal, but conceptually and empirically it should be treated as a distinct entity.
while stromal is virtually defined by the dtf signature <cit>  neither dcn nor fbln <dig> are found among the  <dig> genes in the original list of genes that discriminate between dtf and sft  <cit> .
in contrast, dcn and fbln <dig> are both prominent on the list of myoepihelial genes identified by allinen et al.
in a sage-based sequential purification of cell types  <cit> .
therefore, stromal might better be labelled "myoepithelial".
in keeping with reports of the positive effect of myoepithelial cells in co-cluture experiments  <cit> , stromal/decorin is highly significant in all four of the datasets, though in the nki <dig> it is apparent that it merges with genes from stromal.
this may be partially explained by the fact that dcn is spotted only once on the agilent chip, while it is represented by four probe sets on the affymetrix u133a.
as an experiment, if we used the single agilent decorin probe as a surrogate for stromal.
defined in this way the stromal/decorin gene set proves exceptionally significant as a prognosticator of increased survival in the nki <dig> .
across the four datasets the substantive finding is that stromal, principally decorin, appears to be strongly associated with increased survival for patients at risk of early onset metastasis.
in sum, of the six stromal gene sets, five are significantly associated with increased survival in one or more of the data sets.
this appears to be further evidence of the normalizing effect of stroma.
coculture/coinjection experiments have convincingly shown that aberrant stromal cells are required to promote tumor formation in epithelial cells, and the reverse has also been shown, namely that tumorigenic epithelial cells can revert to normal in the presence of normal stroma  <cit> ).
decomposing stromal signalling into its constituent gene sets  may contribute to an understanding of the complex of cell types and signals that comprise the mircroenvironment which is an active participant in the initiation and progression of cancer  <cit> .
CONCLUSIONS
a research program of hierarchical clustering and data visualization has proven immensely productive for more than a decade  <cit> .
but continued reliance on this research script dependent on hierarchical clustering may inhibit the further exploration of large genomic data sets.
here we offer an alternative program for the unsupervised exploration of microarray data, one that delivers all that the standard script delivers plus considerably more.
the paper is merely a proof of the concept that an enumeration of all the genes-by-samples sets in breast cancer is computationally feasible and substantively useful.
by enumerating the gene sets in three data sets that use the industry standard affymetrix u133a, we identify gene sets that are conserved across experiments on a single platform.
by enumerating the gene sets in a forth dataset, the nki data spotted on a custom agilent chip, we identify gene sets that are conserved across platforms.
in terms of substantive results, nearly 40% of the sets detected prove to be significantly associated with survival.
these include subsets of immune response and stromal signalling genes involved in the complex interactions between the epithelial tumor and its microenvironment.
from the first microarray-based studies of cancer, a fundamental challenge has been that of data reduction.
the task is to filter and factor these large matrices such that statistical modelling is possible.
over the course of a decade, expression-based analysis has progressed from single-factor designs, e.g., regressing groups  on survival  <cit> , to two-factor designs, as exemplified by studies that first stratify on estrogen status , then proceed to identify sets or modules of immune genes that correlate with survival  <cit> .
using the gene sets induced by the enumerations, at this point the stage may be set for exploring the effects and interactions of multiple sets of genes as determinants of the progress and outcome of disease.
these sets, or more specifically the partitions they induce, can be conveniently incorporated into  decision trees.
in this way, the best results to date regarding, for example, estrogen status and immune response, can be further qualified and extended via-a-vis important additional factors, such as stromal-status.
the rank-based matching algorithm , can be applied to any matrix of expression values.
as such, it should be useful for the exploration of a broad range of biological and medical data.
limitations of the algorithm are addressed in additional file  <dig>  describing the method as a "direct" approach is a direct reference to hartigan's original work with two-way clustering  <cit> .
input is a matrix of real numbers; output is a set of submatrices that can be read off the original data matrix .
because it works with the entire data set, it dispenses with the standard data reduction steps, e.g., restricting analysis to "differential" genes.
multiple testing is effectively controlled by classical counting rules.
the submatrices  discovered can be immediately interpreted using the original variables, namely, in the present context of breast cancer, in terms of genes and tissue samples.
because the clusters are "two-way", they reveal an association between possible biological mechanisms embodied in the gene sets and subsets or subclasses of breast cancer.
authors' contributions
dm is responsible for concept and original implementation.
dm and jm wrote the paper.
both authors read and approved the final manuscript.
supplementary material
additional file 1
gene sets significant for survival in the uppsala data set.
click here for file
additional file 2
<dig> gene sets detected in the uppsala  data set.
click here for file
additional file 3
<dig> gene sets detected in the stockholm  data set.
click here for file
additional file 4
<dig> gene sets detected in the transbig  data set.
click here for file
additional file 5
<dig> gene sets detected in the nki  data set.
click here for file
additional file 6
seven immune gene sets mapped to seven immune metagenes .
click here for file
additional file 7
stromal genes mapped to dtf genes .
click here for file
additional file 8
limitations of the enumerations as implemented.
click here for file
acknowledgements
we are grateful to two anonymous reviewers who have become virtual co-authors.
we appreciate financial support from wendell mefford who helped defray the cost of publication.
dm thanks martin crane for timely encouragement.

 detecting genes-by-samples patterns in expression data after removing genes that exhibit little variance, the standard script for exploring microarray data applies two-way hierarchical clustering, followed by a visual search for patterns displayed in a red-green heatmap. for breast cancer in particular, this procedure has proven immensely productive. it can be credited with the discovery of the basal subtype, and, more broadly the identification of subtypes of breast cancer that hold out the potential to inform clinical practice despite its utility, the standard script suffers from several limitations, in particular the instability of the binary tree of the clusters found. perturbation and re-sampling techniques are available to gauge the robustness of the clusters defined by subtrees. but, small changes in the selection of genes or choice of samples can result in disconcertingly large changes in the overall configuration of the tree, which calls into question any typology defined on such tree-based partitions. a different problem stems from the disproportionate impact of large, tightly coordinated clusters on the overall arrangement of the tree. because the largest gene sets, for example estrogen or immune response, will dominate the branching of the tree, smaller sets may be broken up and redistributed. the problem of large dense sets of genes occluding smaller sets arises in the simple or one-way application of hc; it is compounded when two trees are visually crossed in the two-way hc used in the standard script. the question then becomes one of what can be faithfully represented in the two-dimensional arrangement. the brief answer, as spelled out by hartigan in the context of "direct" clustering, is that clusters jointly defined by two trees can be rendered as contiguous regions only if they are either disjoint or nested. an alternative to two-way hierarchical clustering, biclustering, seeks to find submatrices in the array of expression values that satisfy some defining criteria. essential to biclustering is the notion that, given a matrix of expression values, the pattern of coordinated expression for a group of genes may be confined to only a subset of samples. that is, the pattern is "local", and may not be detectable using a "global" measure of the similarity between pairs of genes, e.g., pearson correlation coefficient computed on vectors of expression values. motivated by this concern with detecting "local" structure, several algorithms for detecting biclusters have been advanced. for surveys, see. because of the computational complexity inherent in the biclustering problem, which is provably np-hard under several formal descriptions, these programs abandon exhaustive search and either resort to heuristics or impose bounds on the size of the genes-by-samples submatrices that can be recovered. the method applied in this paper reframes the biclustering problem such that an exhaustive enumeration of all of the submatrices that meet a certain definition can be achieved. included in the resulting inventory are not only such familiar sets as "estrogen" and "erbb2/17q12", but additional sets highly prognostic of survival. some of these newly defined sets may serve alongside "estrogen" etc., as building blocks for future models and classifiers of breast cancer. methods defining gene patterns by partitions on samples hc is an agglomerative method that recursively joins sets of genes or, separately, sets of samples. the starting point is the notion that similar genes are assigned to the same set using a measure of similarity such as correlation. a wholly different approach can be predicated on the notion that two genes belong to the same set provided that they induce the same partition on the set of samples. as an example, the erbb2/17q gene set in figure appears to satisfy such a requirement. here the samples have been ordered by column sum. it is apparent that the largest values of each of these six genes align to a considerable degree. that is, they indicate, or select for, much the same subset of columns. while the notion of genes inducing partitions like this may be intuitive, rendering the notion such that it can be implemented in a search algorithm requires further specifying what constitutes an admissible partition. crucially important is the partition's size and the degree to which the samples induced by any two genes must align for the partition to be considered the "same". these two parameters control a matching rule that can serve to explicitly define what is meant by "gene set". setting up the enumeration the two parameters that define what counts as the "same partition" suggest a partition-based strategy for finding all sets of genes in a matrix of expression values, namely, try all combinations of size and tolerance. because this collects all the genes-by-samples patterns, it constitutes an enumeration. since the size and quality of the partitions to be discovered are not known in advance, to enumerate all of the partitions in a dataset, and the genes that support them, requires trying all combinations of size and tolerance. this sets up a two dimensional grid search with partition size and tolerance as dimensions. given the number of samples, and the number of draws, that is genes, bounds can be set for the parameters and by computing p-values with a counting rule. though the rule for bounding t works well enough with smaller datasets, we find that with data sets of the size of the nki, these p-values are so impossibly small that they provide little guidance in controlling t. an alternative strategy relies on a heuristic suggested by results of an initial spiral search in which s and t are stepped in large increments until gene sets are found, at which point s and t, are incremented, separately, and in small steps. this preliminary scan of the four datasets finds: 1) for most partition sizes s, small values of t fail to detect any gene sets at all, while large t settings yield a single set composed of most, if not all of the genes. between these two extremes, for any given s, there is a relatively narrow range of t that returns multiple distinct gene sets. 2) in size and composition, these sets vary smoothly as s and t are stepped. 3) the number of unique components is surprisingly small, for example, between thirty and forty in the nki data. in light of these observations from the spiral searches, rather than bounding t by a combinatorial counting rule, we choose to set t such that specific gene sets remain distinct. in the enumerations reported in this paper, the set we wish to preserve consists of interferon genes, a clinically significant set that the initial spiral search encountered repeatedly. the resulting "stop rule" says, in effect, that for each partition size, which ranges from a small constant, to n minus that same constant, increment t, starting at until the interferon set disappears as a distinct object as it merges with other immune sets. the motivation for devising this stop rule will be developed further in the context of the gene sets that decompose "immune response". as the algorithm steps through its search it typically finds multiple versions of the same set, which may differ only marginally in the number and composition of the genes. this means that rather than a single list of genes, a "gene set" in fact refers to a set of sets, like those listed for the erbb set in table each column in that table catalogues the genes found for the erbb2/17q gene set in the uppsala data as detected at ten combinations of the "size" and "tolerance" parameters. the list of genes in each column is like a snapshot of the same object in gene space, viewed at a slightly different resolution. the partition size and tolerance are listed at the head of each column. the first instance of this gene set is detected at s = t = as s and t are stepped, the number of erbb2/17q genes increases from size to a maximum of at s = t = from that point the number of genes decreases. it is important to note that the end product of the enumeration is both the list of genes in each set and the partition induced on the samples by that set. while the number of genes in the sets can vary from six to as many as, regardless of the number of genes in the set there can be only one partitioning pattern per gene set. it is in this context that the issue of genes with multiple probe sets can be addressed. as would be expected, multiple probe sets spotted for the same gene will, to the extent that they induce the same partition on the samples, be assigned to the same gene set by the algorithm. a case in point is the small stromal gene set in figure which consists principally of decorin and fibulin. decorin is spotted four times on the affymetrix hg-u133a and the algorithm finds that all four decorin probe sets induce the same partition pattern on the samples in the miller dataset. that same pattern is also induced by two copies of fbln as well as by single instances of glt8d ctsk, prrx and spon and several ests. all in all, in the miller data, this partitioning pattern is found at five combinations of s and t settings. for the purpose of making the survival analysis manageable, we summarize these five versions of this gene set into one composed of those probe sets that appear in at least half of the original versions. the result is the set pictured in the heatmap in the figure, to which metastasis events have been attached. as described more fully in the results section, the χ for a partition at the median value of a vector of column sums is. the χ for first versus last quartile is. since each of the genes in a gene set realize the same partitioning pattern, the end result, in terms of partitioning and survival analysis would be virtually the same if the four decorin probe sets were replaced by any one. this issue of multiple probe sets for the same gene does become a problem for the detection of the smallest gene sets in the case in which the partitioning pattern is realized exclusively in probe sets all of which spot the same gene. this occurs occasionally in the affymetrix data, and we consider these instances to be artefacts of the chip design. detecting patterns that were heretofor undetectable the virtues of hc and the standard script include the fact that it is intuitive, simple to implement, and effective. moreover, the software that implements this method of discovery is free and well maintained, and the output is not only immanently useful, but a thing of beauty. replacing this standard script requires both a conceptual shift and considerable computational apparatus. for these reasons, such a switch can only be justified if the alternative not only reproduces the standard results, but significantly extends those results. that is, it must find additional patterns that are biologically meaningful and clinically significant. to show that this is in fact the case we enumerate the gene sets in four well-studied breast cancer microarray data sets listed in table we refer to these as: uppsala, stockholm, transbig, and nki data sets. three of these use affymetrix u133a chips, a platform that has attained the status of an industry standard, and has been used perhaps in more studies than any other microarray platform. to identify gene sets conserved across platforms, we add the fourth data set, generated by the nederlands kanker instituut, which used a custom agilent chip. the uppsala and stockholm cohorts are population studies. the transbig dataset is a validation series for a gene-based classifier and, accordingly, consists of younger node-negative patients. the nki data was collected in conjunction with a second classifier, and consists of younger patients, both node-positive and node-negative. for purpose of survival analysis in this paper, the end point for the uppsala data is breast cancer specific survival; the end points in the other three data sets are distant metastasis free survival. the datasets were selected, in the first instance, because of the quality of the data, both expression and accompanying clinical data. but they were also chosen to set up explicit comparisons between the gene sets detected by the enumerations in this paper and gene clusters and signatures induced on the same data by other methods and programs. data sets were downloaded from the geo website, or, in the case of the nki data, from the stanford website that accompanies. missing values were estimated using k nearest neighbours,. because the algorithm takes the full expression matrix as input, e.g., in the case of the uppsala data, a matrix of size × there is no prior filtering/selection of genes.  in each of four breast cancer datasets in table the rank-based bi-clustering algorithm finds between thirty and forty sets of genes. many of these sets are conserved across the four studies despite differences in platform: affymetrix u133a, and custom agilent. several of the largest and most prominent of the sets closely resemble gene clusters reported many times in microarray studies of breast cancer. these include "estrogen", "proliferation", erbb and "adipose". other sets, a number of which are relatively small, appear to be novel. these include "histone", "hemoglobin" "amplicon 8p11-12". overall, nearly a third of the gene sets detected are clinically significant as measured by association with survival, and some are very significant with a log-rank χ values exceeding. for purpose of presentation, in the course of enumerating the gene sets in each of the four data sets, we proceed in three steps. the gene sets detected in each of the three affymetrix data sets are first tabulated separately. this is followed by a comparison of these sets to identify gene sets conserved on the affymetrix u133a platform. this, then, is followed by an enumeration of the gene sets in the nki data set which uses a custom agilent chip. these sets are then compared to the affymetrix sets for the purpose of identifying gene sets conserved across the two platforms. regarding the clinical significance of the gene sets, since a gene set, in general, is detected multiple times, the survival analysis of a gene set would yield multiple log-rank values. to simplify the analysis we have replaced these multiple versions of each set with a "core" set comprised of those genes that appear in at least half of the instances detected. five log-rank χ values are reported for each core gene set. the first reflects all the samples in a data set partitioned at the median value of the column sums of the genes in the set. the second partitions all samples by lowest and highest quartiles. the third restricts the samples to those which are er positive, and partitions at the median. the fourth also restricts the samples to er positive, partitioning by quartiles. the fifth restricts the samples to er negative, partitioned at the median. because of small class size, the er negative samples are not partitioned by first and last quartile. dataset 1: uppsala cohort the proliferation set is strongly prognostic of recurrence in the uppsala data, as it is in the other three data sets, which accords with a number of microarray-based studies of breast cancer. other sets strongly prognostic of recurrence are the histones, the metallothioneins/16q as well as gapdh, cd affx-m27830_ and gnas. the histone set includes hist1h2bf hist1h2be hist1h2bh h2bfs hist1h2bk hist1h2bd, while the metallothioneins gene set consists of seven, or perhaps eight, isoforms and two ests: mt1e, mt1f, mt1g, mt h, mt1m, mt1x, mt2a, loc up-regulation of mt2a, the most abundant of these, is associated with more aggressive breast cancer and poor prognosis, and it is reported that the over-expression of metallothioneins predicts resistance to doxorubicin. gnas is located on the long arm of chromosome at 20q, lying just outside the interval of the 20q amplicon investigated by ginestier et al., but falling within one of the two 20q amplicons analyzed by yao et al. it is reported that increased copy numbers for 20q amplicon genes occur in 12% of primary breast tumors, and are associated with more aggressive disease. gapdh is often used as a housekeeping gene. for example it is one of the five genes used to normalize recurrence score in oncogenedx. but, gapdh is reported to be up-regulated in some cancers, e.g., by a factor of to in non small cell lung cancer compared to normal lung tissue. in a study of gapdh expression in breast cancer, revillion et al., conclude with the warning that it should not be used as a control rna. valenti et al reach a similar conclusion. the fact that the gapdh partitions the uppsala samples into good and poor prognostic groups is added evidence that it is not only unsuitable for scaling expression values, but may in fact be a candidate oncogene. data set 2: stockholm cohort similar to the results for the uppsala data, for the stockholm cohort three stromal gene sets are strongly associated with increased survival as reported in table also, again, the adipose set proves significant. in addition, among the sets associated with increased survival is a hemoglobin set consisting of: hbh hbb, hba hbb, hba hba hba hbb, hba hbg. a complete list of the genes in all gene sets for the stockholm cohort is available in additional file though er and pr status are reported in the aggregate in the original article, er status was not available for this data set. consequently the log-rank values tabulated in table represent all samples partitioned at the median, and by first and last quartiles. among the sets significantly associated with decreased survival, as in the uppsala data, proliferation is most prominent. also, again, the histone set, cd and gapdh figure among the sets associated with recurrence. in addition, among the gene sets that have a negative impact on survival in this data set include actg a ribosomal set, and a set we have labelled ezrin, because it includes the gene vil ezrin expression is associated with metastasis in a number of cancers. in a murine cell line elliot et al show that ezrin is required for metastasis, and li et al demonstrated that ezrin silencing reverses cell migration and invasion in a metastatic breast cancer cell line. with regard to survival in breast cancer, bruce et al report that ezrin expression is associated with poor outcome. the ezrin gene set also contains arf which is reported to modulate migration and proliferation in breast cancer cell lines via the regulation of the pi3k pathway. data set 3: transbig while the transbig data employs the same affymetrix platform as the uppsala and stockholm datasets, it differs in terms of sample selection. the first two data sets are population-based, while this data set, as a validation series, reflects the sampling criteria used in the development of wang et al.'s 76-gene classifier. as a consequence, patients were younger, node-negative, with smaller tumour grade. despite these differences, as is apparent in table the gene sets detected are largely the same as those found in the uppsula and stockholm cohorts. a complete list of the genes in each of the sets discovered in the transbig data are available in additional file as in the previous two data sets, stromal is associated with increased survival. distinctive of this third data set are the three ribosomal sets associated with survival, a relationship which for two of the sets and ribosomal) is particularly strong for er positive samples. with regard to gene sets prognostic of decreased survival, again, as with the uppsala and stockholm cohorts, histone and proliferation figure prominently, as do the metallothioneins. also among the sets associated with decreased survival are the basal set and nktr. comparing gene sets across the three affymetrix u133a data sets factors that impact the ability of the algorithm to detect sets that induce a common partition include the design of the microarray platform. some genes, for example cd or lst are rendered by multiple probe sets. to the extent that they take similar expression values for the same samples, the algorithm will identify them as a gene set in their own right. that is, all gene sets that satisfy the matching criteria are reported in the enumeration, including several as small as six genes in size. while some of these will prove to be artefacts of chip design, others may prove to be significant, for example the eight-gene 8p11- gene set detected in the nki data. a virtue of the method over hc is this capacity to find even the smallest sets that induce important patterns on the samples. data set 4: nki 2002 of the gene sets discovered in the enumeration of the nki available in additional file eight are significantly associated with increased survival. of these, the two largest are the estrogen set, well-documented in microarray-based studies of breast cancer, and a second set that we have labelled "foxa1". these sets are disjoint except for the pivotal foxa gene which appears in both, and which is represented by two probes on the custom agilent chip used in the nki data. while the estrogen and foxa gene sets are both significantly associated with increased survival, the two sets induce substantially different partitions on the tumor samples, possibly suggesting different positive mechanisms at work. the gene foxa is known to correlate strongly with estrogen receptor alpha, and the mechanism that accounts for this association has been established, namely, er binding requires foxa binding in close proximity. but, it has also been shown that foxa expression is largely independent of estrogen receptor status. a substantial proportion of er-negative tumors express foxa, while, in one study, more than 40% of er-positive tumors were down-regulated for foxa. additional support for the notion that foxa is independent of er, is the fact that most of the known foxa1-response genes are not erα-response genes. so, paradoxically, it would appear that foxa is both correlated and uncorrelated with er. the estrogen and foxa sets detected in the enumeration, support both propositions. it should be noted that at least two other genes in the foxa set, tloc and sdccag are reported to act as tumor suppressors in their own right. two stromal sets, described in the next section, are shown in table to be associated with increased survival. the genes in the 8p11- gene set essentially tile a region of 8p11- amplicon between m and m. the set largely agrees with the amplicon as refined by haverty et al.. seven of the eight genes map to one of the four 8p11- sub-amplicons that gelsi-boyer et al identified using array-cgh and reyal et al.'s method for correlating genes within 20-gene windows. studies have found that the expression of each of these genes is significantly correlated with copy number. it is interesting to observe that immune sets are associated with survival both positively t-cell and immuneigg) and negatively mhc-i and immuneinterferon). as with the three affymetrix data sets, the histone and proliferation gene sets are strongly prognostic of poor outcome in the nki data. two histone sets, h2b_histone and h3_histone are detected in this data set, but it is only the h2b, which corresponds to the histone set in the affymetrix data, that is significant. gene sets conserved across the agilent and affymetrix platforms discussion while the tables in the results section report the gene sets found and their clinical significance, here the purpose is to compare gene sets, in particular those related to immune response and stromal signalling, to gene clusters and signatures identified across a number of microarray studies. for many of the sets detected the correspondence is immediate. this is the case in particular for "estrogen", "erbb2", "basal", and "proliferation", all of which have become virtual fixtures in genomic studies of breast cancer. other sets, including several small sets found only in the affymetrix data, may be trivial in the sense that they consist of housekeeping genes, or reflect the design decision to spot single genes multiple times. in these cases, the resulting patterns, as detected by the algorithm, are real, but may be of only technical interest in so far as they concern issues of normalization and quality control. of the remaining sets detected by the enumerations, several are easily described and readily interpreted in biological terms. examples include: • hemoglobin: hba hbb, hba hbg.. • histone: hist1h2bf, hist1h2be, hist1h2bfs, h2bfs, hist1h2bk,... • metallotheioneins: mt1g, mt h, mt1x, mt2a, mt1f, mt1e, mt1m of special interest are the ten sets labelled "immune", and the six or more identified as "stromal". a natural factoring of immune response eight immune sets are detected in the uppsala data, five in the stockholm, seven in transbig, and six in nki. all told, across the four data sets, ten unique immune sets are induced. of these, two are questionable: one consisting essentially of five copies of casp and the other, a small set in the nki containing tnfrsf which is most likely a subset of the gene set labelled "immuneigg". setting these aside, the broad category of "immune-related genes" appears to factor naturally into eight distinct subcategories. to marshal evidence for, and against, this division, and to aid in their biological interpretation, these eight sets can be matched against clusters and lists of immune genes identified by others using standard methods, generally hierarchical clustering. the correspondence is one-one between the immune gene sets and eight of loi's et al's "pclusts" and seven of rody et al.'s immune-related metagenes. examining this correspondence in greater detail, seven of the eight "immune" sets correspond almost exactly to the gene clusters used by rody et al. to compute immune "metagenes" in their study of lymphocytic infiltrates in breast cancer. in that study of twelve affymetrix breast cancer data sets, hierarchical clustering repeatedly yielded a cluster of approximately immune-related genes. applying hierarchical clustering a second time to this cluster revealed seven distinct sub-clusters. the concordance between the genes in these seven clusters and seven of the immune sets as detected in the enumeration of the uppsala data is available in additional file the one minor discrepancy between the enumerated sets and these seven clusters involves the five lst probe sets in the monocyte cluster. in the enumerations of the uppsala and stockholm cohorts, lst forms its own small gene set. the survival analysis on the two largest datasets, uppsala and nki, confirm rody et al.'s principal finding, namely t-cell genes are prognostic of increased survival time for estrogen-negative patients, with χ =, p = for the miller data, and χ =, p = for the van de vijver data. but on these same datasets, contrary to rody et al. who find no association between b-cell/immunoglobulin genes and survival, immuneigg is a significant predictor of increased survival in the nki data, χ =, p =, and marginally so in the uppsala data, χ =, p =. the genes in immuneigg closely match the b-cell gene cluster that schmidt et al. find to be prognostic of increased survival among highly proliferating tumor samples. contrary to our results, and to that of rody et al., for tumor samples stratified by proliferation, schmidt et al. fail to find any association between survival and the expression of their t-cell gene cluster. rody et al. attribute this discrepancy to possible differences in cohorts and/or treatments, but the difference in results may also be due to the composition of the respective t-cell gene clusters or sets. schmidt et al.'s t-cell cluster contains our immunet-cell genes as a proper subset, but it also contains genes that belong to immunestat and immunecomplement, neither of which we find to be significantly associated with survival. hence, the significance of the t-cell genes as a predictor of survival may be attenuated by genes from other sets. immediately relevant to the question of t-cell genes versus igg genes as predicators of survival, calabro et al., assembled eighteen genes from the literature to measure the presence of lymphocytic infiltrate. they found that this list is somewhat associated with diminished survival time for estrogen-positive samples, but is strongly prognostic of increased survival for estrogen-negative samples. the genes in this list range from ccl cd cde..to ighg and igj, which effectively merges our immunet-cell and immuneigg sets. therefore, the separate positive associations with survival between immunet-cell and immuneigg for estrogen-negative samples in the uppsala and nki datasets appear to confirm calibro et al.'s results. the positive effect of the genes in immunet-cell on survival may in part explain the performance of finak et al.'s stromal-derived classifier. in that study, hierarchical clustering applied to individually matched tumor and normal stroma yielded three clusters of samples that are starkly distinguished by outcome. from the genes that best discriminate between pairs of these clusters, finak et al. construct a 23-gene classifier. ten of these are associated with the good outcome cluster of samples. of these ten, eight belong to the immunet-cell set: gzma, cd8a, cd cd cd plek, runx and gimap this suggests that the stromal-derived classifier, in assigning samples to the good, poor, and mixed groups, may be powered substantially by the association between t-cell genes and increased survival. if this is the case, then the stromal-derived classifier provides yet more evidence for the association between t-cell genes and survival for estrogen-negative breast cancer observed in the uppsala and nki data. the capacity of teschendorff et al.'s seven gene immune response module to identify good outcome samples from among estrogen-negative tumors may constitute yet more support for this association. those seven genes span at least four of our immune sets with immunedt-cell represented by ly immuneigg by iglc and tnfrsf immunemhc-i by hla-f, and immunecomplement by c1qa. it is interesting to note that in the attempt to validate this seven-gene module on independent data, only four of the seven genes prove significant, and of these, three belong to either immunet-cell, or to immuneigg. if these are the genes that are driving the performance of the teschendorff et al. immune response module, then the effectiveness of that signature for predicting increased survival among estrogen-negative samples is consistent with the immunet-cell and immuneigg survival results for estrogen-negative samples in the uppsala and nki data. overall, the strong association between the t-cell gene set and survival, and the milder association between the immunogloblulin/b-cell gene set and survival, as identified in the enumerations, appears to converge with the important results of each of these several studies despite large differences in approach and research design. immune/interferon the immune/interferon set closely resembles the cluster of "interferon response" genes identified by buess et al.. in that experiment fibroblasts were co-cultured with several breast cancer cell lines to investigate cell-cell signalling between stroma and malignant epithelial cells. in comparing the gene expression of co-cultures to matched monocultures, the starkest difference involved interferon-related genes which were preferentially expressed in co-cultures of fibroblasts with estrogen negative cell lines. of the twelve genes in the immune gene set, as realized for example in the nki data, five belong to this "interferon response": oas oas mx mx and ifit the remaining seven include: ifit isg os mtap4a, usp g1p and gs buess et al show that these interferon genes are significantly associated with shorter survival time in the nki data. the enumeration of the nki confirms this finding. the stromal gene sets factor stromal signaling efforts to delineate the interaction between stromal cells and epithelial tumor are challenged by the complexity of the microenvironment, which has been defined as all components of the mammary gland other than luminal and/or tumor epithelial cells. the effective division of the stroma-related genes into six gene sets by the enumerations may be relevant to this problem. unlike the immune sets which appear to be distinct, monolithic entitites, the stromal sets resemble a constellation with a core entity, which we designate "stromal", accompanied by several satellite sets, "stromal, stromal, stromal and stromal. each of these is detected as a self-contained and independent set under at least one combination of partition "size" and "tolerance" parameters, but, as the algorithm increases the partition size, and relaxes the stringency of what qualifies as a match, these sets tend to quickly merge into an omnibus stromal set. despite this, the importance of distinguishing the smaller sets becomes apparent in the survival analysis: five of the six stromal gene sets are associated with increased survival, several exceptionally so with log-rank χ values in excess of stromal = core dtf stromal signature the stromal set, induced in each of the four enumerations, is composed of many, if not all of the collagen and ecm remodelling enzymes featured in west et al's desmoids-type fibromatosis stromal signature. prominent genes include: sparc, cspg fbln fbn and type-i, type-iii, and type-vi collagen genes. that signature was devised as a proof of concept for a larger, on-going program that exploits the mono-cellular property of soft tissue tumors to inductively define subtypes of fibroblastic stroma cells. the original dtf signature, comprised of genes differentially expressed between two types of soft tissue tumors, was refined for the purpose of identifying a distinctive stromal response in breast cancer. in five datasets, including three of the four used in this paper, the dtf stromal response identifies a subset of breast cancer patients who experience increased survival. the two versions of the dtf signature contain and genes, respectively. in an analysis of the functional relations among the proteins that correspond to these genes, this list is reduced further to a protein-protein network of genes. the close relationship between the dtf signature and the stromal gene set as realized in each of the four data sets can be conveyed with respect to these twenty essential genes. stromal/col11a1 among the six stromal sets, stromal/col11a is the only one that is negatively related to survival, though this is apparent only when inspecting the entire family of sets detected. in the enumerations, this gene set, consisting exclusively of col11a and fn probe sets, is almost always subsumed in a larger stromal set composed essentially of stromal/dtf genes. nevertheless it appears to be a distinct entity under some of the parameter combinations that control the stringency of the match. it is interesting to observe that the survival value of the stromal set, as measured by log-rank, and which is generally positive, abruptly goes to zero as coll11a merges with that set. this scenario is played out in all four datasets. the implication is that important relationships between stromal expression and outcome can be lost unless the stromal genes are decomposed into their constituent sets. stromal/lama2 the stromal gene set is detected in the three affymetrix data sets and is characterized by lama and col14a the genes in this gene set tend to merge with the genes in stromal/darc as the size of the partition is increased and the stringency of the match is relaxed. stromal/lama is significantly associated with increased survival: χ for the uppsala χ, p = for stockholm and χ, p = for transbig stromal/darc the stromal/darc gene set is detected in the uppsala and stockholm cohorts and is positively associated with survival in both: χ for the uppsala χ for the stockholm genes in this set include: darc, tnxb, ccl ldb lhfp, and c duffy antigen receptor for chemokines is believed to suppress tumor metastasis through two mechanisms: by sequestering angiogenic chemokines and by inducing kai1/cd tetraspandin in tumor cells, thereby causing senescence. the merozoite of p. vivax malaria enters red blood cells via darc, consequently individuals who lack darc are resistant to that strain of malaria. the association between the lack of expression of darc in a large percentage of black men and increased rates of aggressive prostate cancer, compared to whites, has been established. a similar relationship appears to hold between black women and the increased incidence of more aggressive breast cancer. it is estimated that 70% of blacks of west african descent lack the expression of darc, which is a population that suffers high rates of both prostate and breast cancer. a second gene in the stromal gene set, tenascin b, is down-regulated during the tumor progression of neurofibromatosis. stromal/adamts5 among the smallest of the stromal sets, stromal is comprised of adamts znf, and four agilent probes that lack gene names. it is significantly associated with increased survival in the nki data: χ, p =, when partitioned at the median, χ, p =, when partitioned by first and last quartiles. stromal/decorin stromal consists almost exclusively of dcn and fbln probe sets. like the other small stromal sets, with the exception of stromal/adamts this set tends to merge with stromal, but conceptually and empirically it should be treated as a distinct entity. while stromal is virtually defined by the dtf signature neither dcn nor fbln are found among the genes in the original list of genes that discriminate between dtf and sft. in contrast, dcn and fbln are both prominent on the list of myoepihelial genes identified by allinen et al. in a sage-based sequential purification of cell types. therefore, stromal might better be labelled "myoepithelial". in keeping with reports of the positive effect of myoepithelial cells in co-cluture experiments, stromal/decorin is highly significant in all four of the datasets, though in the nki it is apparent that it merges with genes from stromal. this may be partially explained by the fact that dcn is spotted only once on the agilent chip, while it is represented by four probe sets on the affymetrix u133a. as an experiment, if we used the single agilent decorin probe as a surrogate for stromal. defined in this way the stromal/decorin gene set proves exceptionally significant as a prognosticator of increased survival in the nki. across the four datasets the substantive finding is that stromal, principally decorin, appears to be strongly associated with increased survival for patients at risk of early onset metastasis. in sum, of the six stromal gene sets, five are significantly associated with increased survival in one or more of the data sets. this appears to be further evidence of the normalizing effect of stroma. coculture/coinjection experiments have convincingly shown that aberrant stromal cells are required to promote tumor formation in epithelial cells, and the reverse has also been shown, namely that tumorigenic epithelial cells can revert to normal in the presence of normal stroma ). decomposing stromal signalling into its constituent gene sets may contribute to an understanding of the complex of cell types and signals that comprise the mircroenvironment which is an active participant in the initiation and progression of cancer.  a research program of hierarchical clustering and data visualization has proven immensely productive for more than a decade. but continued reliance on this research script dependent on hierarchical clustering may inhibit the further exploration of large genomic data sets. here we offer an alternative program for the unsupervised exploration of microarray data, one that delivers all that the standard script delivers plus considerably more. the paper is merely a proof of the concept that an enumeration of all the genes-by-samples sets in breast cancer is computationally feasible and substantively useful. by enumerating the gene sets in three data sets that use the industry standard affymetrix u133a, we identify gene sets that are conserved across experiments on a single platform. by enumerating the gene sets in a forth dataset, the nki data spotted on a custom agilent chip, we identify gene sets that are conserved across platforms. in terms of substantive results, nearly 40% of the sets detected prove to be significantly associated with survival. these include subsets of immune response and stromal signalling genes involved in the complex interactions between the epithelial tumor and its microenvironment. from the first microarray-based studies of cancer, a fundamental challenge has been that of data reduction. the task is to filter and factor these large matrices such that statistical modelling is possible. over the course of a decade, expression-based analysis has progressed from single-factor designs, e.g., regressing groups on survival, to two-factor designs, as exemplified by studies that first stratify on estrogen status, then proceed to identify sets or modules of immune genes that correlate with survival. using the gene sets induced by the enumerations, at this point the stage may be set for exploring the effects and interactions of multiple sets of genes as determinants of the progress and outcome of disease. these sets, or more specifically the partitions they induce, can be conveniently incorporated into decision trees. in this way, the best results to date regarding, for example, estrogen status and immune response, can be further qualified and extended via-a-vis important additional factors, such as stromal-status. the rank-based matching algorithm, can be applied to any matrix of expression values. as such, it should be useful for the exploration of a broad range of biological and medical data. limitations of the algorithm are addressed in additional file describing the method as a "direct" approach is a direct reference to hartigan's original work with two-way clustering. input is a matrix of real numbers; output is a set of submatrices that can be read off the original data matrix. because it works with the entire data set, it dispenses with the standard data reduction steps, e.g., restricting analysis to "differential" genes. multiple testing is effectively controlled by classical counting rules. the submatrices discovered can be immediately interpreted using the original variables, namely, in the present context of breast cancer, in terms of genes and tissue samples. because the clusters are "two-way", they reveal an association between possible biological mechanisms embodied in the gene sets and subsets or subclasses of breast cancer. authors' contributions dm is responsible for concept and original implementation. dm and jm wrote the paper. both authors read and approved the final manuscript. supplementary material additional file 1 gene sets significant for survival in the uppsala data set. click here for file additional file 2 gene sets detected in the uppsala data set. click here for file additional file 3 gene sets detected in the stockholm data set. click here for file additional file 4 gene sets detected in the transbig data set. click here for file additional file 5 gene sets detected in the nki data set. click here for file additional file 6 seven immune gene sets mapped to seven immune metagenes. click here for file additional file 7 stromal genes mapped to dtf genes. click here for file additional file 8 limitations of the enumerations as implemented. click here for file acknowledgements we are grateful to two anonymous reviewers who have become virtual co-authors. we appreciate financial support from wendell mefford who helped defray the cost of publication. dm thanks martin crane for timely encouragement.
while it has delivered important insights, including a typology of breast cancer subtypes, it suffers from instability in the face of gene or sample selection, and an inability to detect small sets that may be dominated by larger sets such as the estrogen-related genes in breast cancer. many of these sets are readily identified in biological terms, e.g., "estrogen", "erbb2", and 8p11- and several are clinically significant as prognostic of either increased survival or diminished survival. two-way hierarchical clustering, with results visualized as heatmaps, has served as the method of choice for exploring structure in large matrices of expression data since the advent of microarrays. it delivers results comparable to two-way hierarchical clustering, and much more. applied systematically across a range of parameter settings, it enumerates all the partition-inducing gene sets in a matrix of expression values. of special interest are the sets that effectively factor "immune response" and "stromal signalling". applied to four large breast cancer datasets, this alternative exploratory method detects more than thirty sets of co-regulated genes, many of which are conserved across experiments and across platforms. in this regard these inventories confirm and consolidate findings from microarray-based work on breast cancer over the last decade. the sets induced are robust, biologically meaningful, and serve to reveal a finer structure in existing breast cancer microarrays. but, the enumerations also identify gene sets that have not been studied as of yet, some of which are prognostic of survival.
 two-way hierarchical clustering, with results visualized as heatmaps, has served as the method of choice for exploring structure in large matrices of expression data since the advent of microarrays. while it has delivered important insights, including a typology of breast cancer subtypes, it suffers from instability in the face of gene or sample selection, and an inability to detect small sets that may be dominated by larger sets such as the estrogen-related genes in breast cancer. the rank-based partitioning algorithm introduced in this paper addresses several of these limitations. it delivers results comparable to two-way hierarchical clustering, and much more. applied systematically across a range of parameter settings, it enumerates all the partition-inducing gene sets in a matrix of expression values.  applied to four large breast cancer datasets, this alternative exploratory method detects more than thirty sets of co-regulated genes, many of which are conserved across experiments and across platforms. many of these sets are readily identified in biological terms, e.g., "estrogen", "erbb2", and 8p11- and several are clinically significant as prognostic of either increased survival or diminished survival. of special interest are the sets that effectively factor "immune response" and "stromal signalling".  the gene sets induced by the enumeration include many of the sets reported in the literature. in this regard these inventories confirm and consolidate findings from microarray-based work on breast cancer over the last decade. but, the enumerations also identify gene sets that have not been studied as of yet, some of which are prognostic of survival. the sets induced are robust, biologically meaningful, and serve to reveal a finer structure in existing breast cancer microarrays. 