BACKGROUND
one out of ten women who reaches the age of ninety will have had breast cancer in her lifetime.
most tumors are treated with a combination of surgery, radiation therapy, and adjuvant systemic therapy .
60–80% of tumors express the estrogen receptor er and respond to treatment with hormonal agents such as aromatase inhibitors or tamoxifen  <cit> .
20–40% have amplification of the her <dig> gene  <cit>  which is a marker of increased recurrence rates and poorer prognosis.
the outcome of these her2+ tumors can be improved by the addition of the humanized anti-her <dig> antibody trastuzumab  to their treatment regimen.
10–15% of tumors neither express the estrogen receptor nor harbor her <dig> amplification and have a characteristic gene expression profile  <cit> .
these cancers, called basal-like  <cit> , are high grade aggressive malignancies with poor overall prognosis, and at present there is no targeted therapy for them.
in spite of these classifications and treatment choices, therapy is confounded by the fact that tumors with similar histopathology often have divergent course and varied response to therapy  <cit> .
microarrays have the potential to shed light on this picture because of their ability to provide a snapshot of the genetic state of the cell.
in principal, they should be able to identify the genes and pathways altered in cancer initiation, progression and metastasis.
this promise has resulted in microarray technology being aggressively pursued by researchers, hospitals and pharmaceutical companies to get an improved understanding of the disease process, better diagnostic protocols, new drugs, and new treatment regimens.
however, the success of these efforts has been limited by practical considerations.
the biggest limitation is that the results from microarray studies are sensitive to noise and analysis method  <cit> .
this often leads to ambiguous results and biologically non-intuitive genes and pathways for stratification  <cit> .
efforts to use microarray data to identify the underlying biology of disease progression and help characterize the disease phenotype have met with limited success.
in this paper, we develop and give results from a robust method which addresses the issues outlined above.
we first use principal component analysis   <cit>  to identify the overall structure of clusters in the data and to select the subset of genes that distinguish the clusters.
we then use this set of genes and a new consensus ensemble k-clustering technique, which averages over several clustering methods and many data perturbations, to identify strong, stable clusters.
we also define a simple criterion to find the optimum number of clusters and a method to identify robust markers for disease progression within each cluster.
RESULTS
applied to a breast cancer microarray data set, our method results in stable lists of genes and pathways that distinguish high and low grade tumors.
it also identifies other robust gene sets which mark progression of disease from dcis or ductal carcinoma in-situ to idc or invasive ductal carcinoma.
the clusters paint a portrait of the disease at varying levels of granularity.
when the data is divided into two clusters, the normal samples form one cluster and the disease samples form another.
at the next level of clustering, the low grade and high grade samples separate.
the optimal number of clusters is seven, corresponding to two sub-clusters  of the low grade samples and four  of the high grade samples.
these sub-clusters are well separated by a strong set of markers which are able to distinguish them with sensitivity and specificity in the 80–100% range.
we identify the genes and pathways that mark disease progression in each sub-cluster.
a major result of our analysis is that each sub-cluster contains samples from non-invasive and invasive tumors from the same patient.
this suggests that within each grade of breast cancer, different groups of patients progress to the same final phenotype along different pathways.
this result suggests that the sub-clusters identified here are distinct diseases.
if validated on larger datasets with larger gene-sets, it would have significant implications for disease identification and treatment.
using the genes specific to each cluster and er, pr, her2+ levels in the data, we can place the clusters into the standard categories used to classify breast cancer as defined in  <cit> .
we find that the low grade clusters correspond to one luminal a subtype and one luminal b subtype.
the high grade samples correspond to two additional luminal b subtypes, one her2+ subtype and one basal-like subtype.
description of data
the data was obtained from  <cit>  and consisted of samples from  <dig> breast cancer patients of which  <dig> were diagnosed with at least two out of three pathological stages of disease: atypical ductal hyperplasia or adh, ductal carcinoma in situ or dcis and invasive ductal carcinoma or idc respectively.
the remaining  <dig> patients were diagnosed to have pre-invasive disease  only.
microarray analysis was also done for samples collected from normal breast epithelial tissue extracted from three healthy women during routine mammoplasty.
from the cancer patients, normal as well as disease samples were collected from as many different stages of disease  as possible.
these samples were harvested in triplicate using laser capture micro-dissection , taking care to avoid contamination between cells of different stages from the same patient.
each sample was analyzed in duplicate with a  <dig>  gene cdna microarray.
it was determined that the "normal cells" from cancer patients were highly similar to the normal epithelium of the three disease free patients.
this suggested that the normal samples from cancer patients could be used as a baseline to determine disease state and progression.
the data provided in  <cit>  consisted of the expression levels of  <dig> genes across  <dig> samples.
<dig> of these were from disease free or normal tissue,  <dig> were adh samples,  <dig> were dcis samples and  <dig> were idc samples.
the  <dig> genes  were selected in  <cit>  by their ability to distinguish "normal cells" and each of the disease stages adh, dcis and idc using a linear discriminant function.
the patients were further classified by pathological analysis into  <dig> categories based on the tumor grade: grade i , grade ii  and grade iii .
the mapping of sample labels to stage and grade and a patient identifier is given in additional file  <dig>  this table also includes the classification of the samples into the disease subtypes using the methods described in this paper.
the flow chart of our analysis method is presented in figure  <dig>  first the dataset was normalized and missing entries imputed robustly.
next, pca was used to find the genes which accounted for most of the variation in the data.
the optimal number of clusters kopt in the data was estimated using gap statistics  <cit>  and silhouette scores  <cit> .
next a variety of clustering techniques and data perturbations are averaged to divide the data into  <dig> ,...kopt clusters.
we estimate the number of clusters using the silhouette scores and gap statistics.
this provides a range for the number of clusters.
in our case, the range was 6– <dig> so we chose kmax =  <dig>  this range is compared to the optimum number of clusters provided by the em mixture modeling approach.
.
we then create ensemble consensuses of k clusters for each k =  <dig>   <dig> ..kmax as follows:
for each k, we integrate the clustering results across clustering methods and perturbations into an agreement matrix.
to obtain a global optimum clustering solution, we apply simulated annealing to the agreement matrix at each k to sort the samples into the best k clusters, which appear as blocks along the diagonal in the sorted agreement matrix.
the agreement matrix cost function used in simulated annealing optimizes between the similarity inside the clusters and the dissimilarity between different clusters for each k. once the data is clustered by simulated annealing into k pieces, we compute the silhouette scores for these k clusters  using the fpc package from r.
we stop at k = kopt clusters beyond which the "quality" of the clusters obtained from the agreement matrix begins to deteriorate.
this is assessed by analyzing the accuracy of cluster membership assignments using weighted voting and knn on single and multiple gene markers  which distinguish the clusters.
we stop clustering when the accuracy of such an assignment falls below 75%.
principal component analysis
pca showed that 50% of the variation in the data was represented by the first  <dig> pcs and 85% by the first  <dig> pcs.
we identified  <dig> genes as those with highest absolute value  in the coefficients of the first  <dig> eigenvectors as representative of most of the data variation.
these thresholds were estimated through a calibration step whose aim was to optimize the overall cluster membership assignment for the optimal number of clusters identified in the data restricted to the selected genes.
consensus ensemble k-clustering
gap statistic and the silhouette scores  <cit>  estimated k =  <dig> as the optimal value for the number of clusters in the data.
the data was divided into k =  <dig>   <dig> ..., <dig> clusters by using the  <dig> genes identified by pca and by applying consensus ensemble k-clustering .
the results are shown schematically in figure  <dig>  at k =  <dig>  the samples separated into a "normal"  group, which contained all the normal samples and one adh sample , and a "breast cancer"  group, which contained all the remaining breast cancer samples.
at k =  <dig>  the normal group was unaltered but the bca group split into a low grade  tumor group containing  <dig> samples labeled grade  <dig> and  <dig> samples labeled grade  <dig>  and a high grade  tumor group containing  <dig> samples labeled grade  <dig> and  <dig> samples labeled grade  <dig>  as k increased progressively from  <dig> through  <dig>  the lg group split into  <dig> distinct subgroups  and the hg group split into  <dig> distinct subgroups .
even though the clusters at each level were determined independently, at clustering level k+ <dig>  two clusters always emerged as splits of a parent cluster at level k, while the remaining k- <dig> clusters were inherited unchanged from the previous level k. this shows that the data inherently contains a hierarchy of detail, providing portraits of disease at different levels of clustering resolution.
the separation of samples into "normal" and "disease" at k =  <dig>  the split of the disease samples into "low" and "high" grades at k =  <dig> and so on, strongly suggests that disease progression is a hierarchical process and is readily and robustly identifiable by our clustering procedure.
genes discriminating low and high grade tumors
using a non-stringent signal-to-noise-ratio  test  we found  <dig> gene markers which distinguish the group lg from hg.
a subset of  <dig> markers was selected based on their performance on leave-one-out cross-validation experiments for weighted voting  and k-nearest neighbors  classification models.
the models trained on these  <dig> markers produced only  <dig> false positive error  and  <dig> false negative error  in leave-one-out experiments.
using the snr test and leave-one-out experiments for the wv and knn models, we identified  <dig> markers which distinguish the lg samples from all others  with 90% accuracy.
we find that rbsk, homo sapiens cdna flj <dig> fis, clone nt2rp <dig> and crip <dig> are up-regulated in the lg group, and eya <dig>  anxa <dig>  runx <dig>  dkfzp762a <dig>  gprc5b are down-regulated in the lg group.
for the high grade cluster, the classification accuracy was 97% with  <dig> false positive and zero false negative errors.
the top markers up-regulated in hg are tram, hspc <dig>  tacc <dig>  cdkn <dig>  ube2c, and top markers down-regulated in hg are x <dig>  gng <dig>  sh3bgrl <dig>  loh11cr2a and homo sapiens, clone image: <dig> mrna, partial cds.
low grade substructure
high grade substructure
as seen in table  <dig>  all the samples in the hg <dig> subgroup were er and pr negative while those in the hg <dig> and hg <dig> subgroups were mostly er and pr positive.
the hg <dig> samples had mixed er and pr signatures.
the hg <dig> subgroup, which is the worst prognosis group based on clinical characteristics, had as discriminatory markers the genes bcl <dig> , rad <dig> , gstp <dig>  and rrm <dig> .
hg <dig> markers also include up-regulated bcl <dig>  and down-regulated rrm <dig>
the hg <dig> markers, include a group of down-regulated genes in chromosomal region 17q23- <dig> which harbors the erbb <dig> amplicon 17q  <dig> .
these genes are kpna <dig> , amplified in breast cancer  <dig> , bcl- <dig> inhibitor of transcription , hypothetical protein tanc <dig> , and two proteosome protein psmc <dig>  and psmd <dig> .
this suggests the possibility that patients in the hg <dig> subgroup might have a re-arrangement or deletion of genes around the her <dig> gene leading to loss of regulation or function for these genes which might explain why only 15% of hg <dig> patients are her2+, while 53% are her2- and 15% are undetermined.
since ihc data for the erbb <dig> gene was not made available in  <cit> , the down-regulation of the erbb <dig> amplicon genes in the hg <dig> samples identifies its clinical signature as her2-.
the most notable hg <dig> marker was a down-regulation in the transforming growth factor beta receptor ii .
mutations in this gene have been associated with the development of various types of tumors.
the over-expression of this gene was found to be associated with poor prognosis breast tumors.
overall, gene markers and clinical parameters lead to the conclusion that among the high grade subgroups hg <dig> is probably the best prognosis group composed of grade ii tumors that are all er+ and pr+.
based on these observations, we identify hg <dig> as basal-like  <cit> , hg <dig> as her2+, and hg <dig> and hg <dig> as additional subtypes of luminal b  <cit> .
figures  <dig> and  <dig> present heatmaps using the top  <dig> upregulated markers which classify the tumors by grade and stage respectively.
since the sample sizes are small, the p values for each classification were obtained using permutation experiments and the fdr rates inferred from these.
the fdr values for the genes in figure  <dig> are  <dig>  for lg,  <dig>  for hg and for the genes in figure  <dig> are  <dig>  for lg <dig>   <dig>  for lg <dig>   <dig>  for hg <dig>   <dig>  for hg <dig>   <dig>  for hg <dig> and  <dig>  for hg <dig>  it should be noted that the sample sizes in this study are small.
consequently, these results are to be considered as hypothesis generating.
all these results should be validated on larger data sets.
figure  <dig> maps the genes identified for progression in different grades into pathways for disease progression using the classification of hanahan and weinberg  <cit> .
discussion and
CONCLUSIONS
the use of ensemble consensus clustering is absolutely critical to distinguish the subtypes.
pca by itself could identify a collection of useful markers, but could not identify the rich stratification discovered by consensus ensemble k-clustering.
hierarchical clustering by itself would separate the samples into clusters but the stratification would be very sensitive to bootstrap, indicating that the clusters are unstable to data perturbation.
robustness of clustering is only obtained by averaging over many clustering techniques and data perturbations as is done in the ensemble consensus clustering technique used here.
our results show that progression of the disease from non-invasive to invasive status occurs along different pathways.
progression in the low-grade groups seems to correlate with changes in metabolic and transportation pathways, while in the high grade groups it is related to alterations in cell-cycle and signaling pathways, with distinct subsets of genes involved in each.
the main observation of the original paper of ma et al  <cit>  was that the molecular signature of breast cancer is already present in the early  stage of the disease.
the genes that distinguish adh from normal progressively change their levels away from normal as the disease progresses to dcis and idc.
they also noticed that that breast cancer progression is defined by distinct markers for low and high grade tumors.
our results, particularly the hierarchy we see when the data is grouped into k =  <dig> ,... <dig> clusters  agree with this observation.
our methods identified six different subtypes of breast cancer with distinct patterns of progression.
looking at the histopathology of the samples in the clusters identified by our analysis, four of these subtypes  have a strongly luminal signature ; one subtype  had the triple negative  characteristic of the basal-like subtype, and one subtype  had a predominantly her2+ signature .
the validation of these subtypes on a larger dataset with more genes is currently underway.
at k =  <dig>  each of the six bca clusters always contained samples in both dcis and idc stages from the same patient.
this strong heterogeneity in the genetic signature of subtypes suggests that breast cancer is composed of distinct disease subtypes that develop early and progress along different pathways because progression within a subtype is less distinct than the subtypes themselves.
treatment decisions may benefit by taking account of these subtypes in addition to the current practice of using the markers er, pr and her <dig>

 one out of ten women who reaches the age of ninety will have had breast cancer in her lifetime. most tumors are treated with a combination of surgery, radiation therapy, and adjuvant systemic therapy. 60–80% of tumors express the estrogen receptor er and respond to treatment with hormonal agents such as aromatase inhibitors or tamoxifen. 20–40% have amplification of the her gene which is a marker of increased recurrence rates and poorer prognosis. the outcome of these her2+ tumors can be improved by the addition of the humanized anti-her antibody trastuzumab to their treatment regimen. 10–15% of tumors neither express the estrogen receptor nor harbor her amplification and have a characteristic gene expression profile. these cancers, called basal-like, are high grade aggressive malignancies with poor overall prognosis, and at present there is no targeted therapy for them. in spite of these classifications and treatment choices, therapy is confounded by the fact that tumors with similar histopathology often have divergent course and varied response to therapy. microarrays have the potential to shed light on this picture because of their ability to provide a snapshot of the genetic state of the cell. in principal, they should be able to identify the genes and pathways altered in cancer initiation, progression and metastasis. this promise has resulted in microarray technology being aggressively pursued by researchers, hospitals and pharmaceutical companies to get an improved understanding of the disease process, better diagnostic protocols, new drugs, and new treatment regimens. however, the success of these efforts has been limited by practical considerations. the biggest limitation is that the results from microarray studies are sensitive to noise and analysis method. this often leads to ambiguous results and biologically non-intuitive genes and pathways for stratification. efforts to use microarray data to identify the underlying biology of disease progression and help characterize the disease phenotype have met with limited success. in this paper, we develop and give results from a robust method which addresses the issues outlined above. we first use principal component analysis to identify the overall structure of clusters in the data and to select the subset of genes that distinguish the clusters. we then use this set of genes and a new consensus ensemble k-clustering technique, which averages over several clustering methods and many data perturbations, to identify strong, stable clusters. we also define a simple criterion to find the optimum number of clusters and a method to identify robust markers for disease progression within each cluster.  applied to a breast cancer microarray data set, our method results in stable lists of genes and pathways that distinguish high and low grade tumors. it also identifies other robust gene sets which mark progression of disease from dcis or ductal carcinoma in-situ to idc or invasive ductal carcinoma. the clusters paint a portrait of the disease at varying levels of granularity. when the data is divided into two clusters, the normal samples form one cluster and the disease samples form another. at the next level of clustering, the low grade and high grade samples separate. the optimal number of clusters is seven, corresponding to two sub-clusters of the low grade samples and four of the high grade samples. these sub-clusters are well separated by a strong set of markers which are able to distinguish them with sensitivity and specificity in the 80–100% range. we identify the genes and pathways that mark disease progression in each sub-cluster. a major result of our analysis is that each sub-cluster contains samples from non-invasive and invasive tumors from the same patient. this suggests that within each grade of breast cancer, different groups of patients progress to the same final phenotype along different pathways. this result suggests that the sub-clusters identified here are distinct diseases. if validated on larger datasets with larger gene-sets, it would have significant implications for disease identification and treatment. using the genes specific to each cluster and er, pr, her2+ levels in the data, we can place the clusters into the standard categories used to classify breast cancer as defined in. we find that the low grade clusters correspond to one luminal a subtype and one luminal b subtype. the high grade samples correspond to two additional luminal b subtypes, one her2+ subtype and one basal-like subtype. description of data the data was obtained from and consisted of samples from breast cancer patients of which were diagnosed with at least two out of three pathological stages of disease: atypical ductal hyperplasia or adh, ductal carcinoma in situ or dcis and invasive ductal carcinoma or idc respectively. the remaining patients were diagnosed to have pre-invasive disease only. microarray analysis was also done for samples collected from normal breast epithelial tissue extracted from three healthy women during routine mammoplasty. from the cancer patients, normal as well as disease samples were collected from as many different stages of disease as possible. these samples were harvested in triplicate using laser capture micro-dissection, taking care to avoid contamination between cells of different stages from the same patient. each sample was analyzed in duplicate with a gene cdna microarray. it was determined that the "normal cells" from cancer patients were highly similar to the normal epithelium of the three disease free patients. this suggested that the normal samples from cancer patients could be used as a baseline to determine disease state and progression. the data provided in consisted of the expression levels of genes across samples. of these were from disease free or normal tissue, were adh samples, were dcis samples and were idc samples. the genes were selected in by their ability to distinguish "normal cells" and each of the disease stages adh, dcis and idc using a linear discriminant function. the patients were further classified by pathological analysis into categories based on the tumor grade: grade i, grade ii and grade iii. the mapping of sample labels to stage and grade and a patient identifier is given in additional file this table also includes the classification of the samples into the disease subtypes using the methods described in this paper. the flow chart of our analysis method is presented in figure first the dataset was normalized and missing entries imputed robustly. next, pca was used to find the genes which accounted for most of the variation in the data. the optimal number of clusters kopt in the data was estimated using gap statistics and silhouette scores. next a variety of clustering techniques and data perturbations are averaged to divide the data into,...kopt clusters. we estimate the number of clusters using the silhouette scores and gap statistics. this provides a range for the number of clusters. in our case, the range was 6– so we chose kmax = this range is compared to the optimum number of clusters provided by the em mixture modeling approach.. we then create ensemble consensuses of k clusters for each k =..kmax as follows: for each k, we integrate the clustering results across clustering methods and perturbations into an agreement matrix. to obtain a global optimum clustering solution, we apply simulated annealing to the agreement matrix at each k to sort the samples into the best k clusters, which appear as blocks along the diagonal in the sorted agreement matrix. the agreement matrix cost function used in simulated annealing optimizes between the similarity inside the clusters and the dissimilarity between different clusters for each k. once the data is clustered by simulated annealing into k pieces, we compute the silhouette scores for these k clusters using the fpc package from r. we stop at k = kopt clusters beyond which the "quality" of the clusters obtained from the agreement matrix begins to deteriorate. this is assessed by analyzing the accuracy of cluster membership assignments using weighted voting and knn on single and multiple gene markers which distinguish the clusters. we stop clustering when the accuracy of such an assignment falls below 75%. principal component analysis pca showed that 50% of the variation in the data was represented by the first pcs and 85% by the first pcs. we identified genes as those with highest absolute value in the coefficients of the first eigenvectors as representative of most of the data variation. these thresholds were estimated through a calibration step whose aim was to optimize the overall cluster membership assignment for the optimal number of clusters identified in the data restricted to the selected genes. consensus ensemble k-clustering gap statistic and the silhouette scores estimated k = as the optimal value for the number of clusters in the data. the data was divided into k =..., clusters by using the genes identified by pca and by applying consensus ensemble k-clustering. the results are shown schematically in figure at k = the samples separated into a "normal" group, which contained all the normal samples and one adh sample, and a "breast cancer" group, which contained all the remaining breast cancer samples. at k = the normal group was unaltered but the bca group split into a low grade tumor group containing samples labeled grade and samples labeled grade and a high grade tumor group containing samples labeled grade and samples labeled grade as k increased progressively from through the lg group split into distinct subgroups and the hg group split into distinct subgroups. even though the clusters at each level were determined independently, at clustering level k+ two clusters always emerged as splits of a parent cluster at level k, while the remaining k- clusters were inherited unchanged from the previous level k. this shows that the data inherently contains a hierarchy of detail, providing portraits of disease at different levels of clustering resolution. the separation of samples into "normal" and "disease" at k = the split of the disease samples into "low" and "high" grades at k = and so on, strongly suggests that disease progression is a hierarchical process and is readily and robustly identifiable by our clustering procedure. genes discriminating low and high grade tumors using a non-stringent signal-to-noise-ratio test we found gene markers which distinguish the group lg from hg. a subset of markers was selected based on their performance on leave-one-out cross-validation experiments for weighted voting and k-nearest neighbors classification models. the models trained on these markers produced only false positive error and false negative error in leave-one-out experiments. using the snr test and leave-one-out experiments for the wv and knn models, we identified markers which distinguish the lg samples from all others with 90% accuracy. we find that rbsk, homo sapiens cdna flj fis, clone nt2rp and crip are up-regulated in the lg group, and eya anxa runx dkfzp762a gprc5b are down-regulated in the lg group. for the high grade cluster, the classification accuracy was 97% with false positive and zero false negative errors. the top markers up-regulated in hg are tram, hspc tacc cdkn ube2c, and top markers down-regulated in hg are x gng sh3bgrl loh11cr2a and homo sapiens, clone image: mrna, partial cds. low grade substructure high grade substructure as seen in table all the samples in the hg subgroup were er and pr negative while those in the hg and hg subgroups were mostly er and pr positive. the hg samples had mixed er and pr signatures. the hg subgroup, which is the worst prognosis group based on clinical characteristics, had as discriminatory markers the genes bcl, rad, gstp and rrm. hg markers also include up-regulated bcl and down-regulated rrm the hg markers, include a group of down-regulated genes in chromosomal region 17q23- which harbors the erbb amplicon 17q. these genes are kpna, amplified in breast cancer, bcl- inhibitor of transcription, hypothetical protein tanc, and two proteosome protein psmc and psmd. this suggests the possibility that patients in the hg subgroup might have a re-arrangement or deletion of genes around the her gene leading to loss of regulation or function for these genes which might explain why only 15% of hg patients are her2+, while 53% are her2- and 15% are undetermined. since ihc data for the erbb gene was not made available in, the down-regulation of the erbb amplicon genes in the hg samples identifies its clinical signature as her2-. the most notable hg marker was a down-regulation in the transforming growth factor beta receptor ii. mutations in this gene have been associated with the development of various types of tumors. the over-expression of this gene was found to be associated with poor prognosis breast tumors. overall, gene markers and clinical parameters lead to the conclusion that among the high grade subgroups hg is probably the best prognosis group composed of grade ii tumors that are all er+ and pr+. based on these observations, we identify hg as basal-like, hg as her2+, and hg and hg as additional subtypes of luminal b. figures and present heatmaps using the top upregulated markers which classify the tumors by grade and stage respectively. since the sample sizes are small, the p values for each classification were obtained using permutation experiments and the fdr rates inferred from these. the fdr values for the genes in figure are for lg, for hg and for the genes in figure are for lg for lg for hg for hg for hg and for hg it should be noted that the sample sizes in this study are small. consequently, these results are to be considered as hypothesis generating. all these results should be validated on larger data sets. figure maps the genes identified for progression in different grades into pathways for disease progression using the classification of hanahan and weinberg. discussion and  the use of ensemble consensus clustering is absolutely critical to distinguish the subtypes. pca by itself could identify a collection of useful markers, but could not identify the rich stratification discovered by consensus ensemble k-clustering. hierarchical clustering by itself would separate the samples into clusters but the stratification would be very sensitive to bootstrap, indicating that the clusters are unstable to data perturbation. robustness of clustering is only obtained by averaging over many clustering techniques and data perturbations as is done in the ensemble consensus clustering technique used here. our results show that progression of the disease from non-invasive to invasive status occurs along different pathways. progression in the low-grade groups seems to correlate with changes in metabolic and transportation pathways, while in the high grade groups it is related to alterations in cell-cycle and signaling pathways, with distinct subsets of genes involved in each. the main observation of the original paper of ma et al was that the molecular signature of breast cancer is already present in the early stage of the disease. the genes that distinguish adh from normal progressively change their levels away from normal as the disease progresses to dcis and idc. they also noticed that that breast cancer progression is defined by distinct markers for low and high grade tumors. our results, particularly the hierarchy we see when the data is grouped into k =,... clusters agree with this observation. our methods identified six different subtypes of breast cancer with distinct patterns of progression. looking at the histopathology of the samples in the clusters identified by our analysis, four of these subtypes have a strongly luminal signature; one subtype had the triple negative characteristic of the basal-like subtype, and one subtype had a predominantly her2+ signature. the validation of these subtypes on a larger dataset with more genes is currently underway. at k = each of the six bca clusters always contained samples in both dcis and idc stages from the same patient. this strong heterogeneity in the genetic signature of subtypes suggests that breast cancer is composed of distinct disease subtypes that develop early and progress along different pathways because progression within a subtype is less distinct than the subtypes themselves. treatment decisions may benefit by taking account of these subtypes in addition to the current practice of using the markers er, pr and her 
clustering analysis of microarray data is often criticized for giving ambiguous results because of sensitivity to data perturbation or clustering techniques used. in this paper, we describe a new method based on principal component analysis and ensemble consensus clustering that avoids these problems. further, the low grade cluster splits into two subclusters and the high grade cluster into four. we illustrate the method on a public microarray dataset from breast cancer patients of whom were diagnosed with at least two of three pathological stages of disease, ductal carcinoma in situ and invasive ductal carcinoma. the first split separates the normal samples from the cancer samples. an intriguing observation is that the disease phenotype is distinguishable in adh and progresses along distinct pathways for each subtype. we confirm that the cancer phenotype can be identified in early stage because the genes altered in this stage progressively alter further as the disease progresses through dcis into idc. our method identifies an optimum set of genes and divides the samples into stable clusters which correlate with clinical classification into luminal, basal-like and her2+ subtypes. the genetic signature for disease heterogeneity across subtypes is greater than the heterogeneity of progression from dcis to idc within a subtype, suggesting that the disease subtypes have distinct progression pathways. we identify six subtypes of disease which have distinct genetic signatures and remain separated in the clustering hierarchy. next, the cancer cluster splits into low grade and high grade while the normal cluster is unchanged. our findings suggest that the heterogeneity of disease across subtypes is higher than the heterogeneity of the disease progression within a subtype, indicating that the subtypes are in fact distinct diseases. our analysis reveals a hierarchical portrait of breast cancer progression and identifies genes and pathways for each stage, grade and subtype.
 clustering analysis of microarray data is often criticized for giving ambiguous results because of sensitivity to data perturbation or clustering techniques used. in this paper, we describe a new method based on principal component analysis and ensemble consensus clustering that avoids these problems.  we illustrate the method on a public microarray dataset from breast cancer patients of whom were diagnosed with at least two of three pathological stages of disease, ductal carcinoma in situ and invasive ductal carcinoma. our method identifies an optimum set of genes and divides the samples into stable clusters which correlate with clinical classification into luminal, basal-like and her2+ subtypes. our analysis reveals a hierarchical portrait of breast cancer progression and identifies genes and pathways for each stage, grade and subtype. an intriguing observation is that the disease phenotype is distinguishable in adh and progresses along distinct pathways for each subtype. the genetic signature for disease heterogeneity across subtypes is greater than the heterogeneity of progression from dcis to idc within a subtype, suggesting that the disease subtypes have distinct progression pathways. our method identifies six disease subtype and one normal clusters. the first split separates the normal samples from the cancer samples. next, the cancer cluster splits into low grade and high grade while the normal cluster is unchanged. further, the low grade cluster splits into two subclusters and the high grade cluster into four. the final six disease clusters are mapped into one luminal a, three luminal b, one basal-like and one her2+.  we confirm that the cancer phenotype can be identified in early stage because the genes altered in this stage progressively alter further as the disease progresses through dcis into idc. we identify six subtypes of disease which have distinct genetic signatures and remain separated in the clustering hierarchy. our findings suggest that the heterogeneity of disease across subtypes is higher than the heterogeneity of the disease progression within a subtype, indicating that the subtypes are in fact distinct diseases. 