BACKGROUND
epigenetics is the study of heritable changes in gene function that cannot be explained by changes in dna sequence  <cit> .
one of the most commonly studied epigenetic alterations is cytosine methylation, which occurs in the context of a cpg dinucleotide.
concentrations of cpgs known as cpg islands, when sufficiently methylated, are associated with transcriptional gene silencing tantamount to "one hit" as part of knudson's two hit hypothesis of tumor suppressor gene inactivation  <cit> .
dna methylation associated gene silencing is a well recognized epigenetic mechanism that often occurs at tumor suppressor gene loci in human cancer.
hundreds of reports of methylation induced silencing at tumor suppressor genes in virtually all types of human cancer have been published  <cit> .
while there has been a tremendous effort to characterize epigenetic alterations in cancer, surprisingly little work has been done in disease-free tissues.
there is a basic need for epigenetic profiling of normal tissues to better understand the contribution of these profiles to tissue specificity, especially in the context of phenotypically important cpgs, where deregulation is associated with human diseases such as cancer.
while efforts to characterize the methylation profiles of normal tissues in humans and mice have begun, and certain themes are slowly becoming apparent, relatively few reports have emerged  <cit> .
most cpgs or cpg regions have been found to have a bimodal distribution of methylation profiles, either hypo- or hypermethylated.
another theme is the disproportionate location of cell or tissue type dependent differentially methylated regions in non-cpg island sites  <cit> .
furthermore, the human epigenome project showed that the human major histocompatability loci have differential methylation across various human tissues, but that differential methylation does not necessarily lead to differential expression  <cit> .
it is therefore critical to first outline the basal-state of phenotypically important epigenetic marks that are known to contribute to cancer in order to have a background for comparison to other normal and diseased tissues.
an approach which characterizes the basal epigenetic state is best suited to foster the discovery of epigenetic profiles that are associated with particular disease states, patient characteristics, exposures or other covariates that may contribute to pathogenesis.
cluster analysis is often used to identify methylation subgroups in data  <cit>  and, in particular, siegmund et al.
argue that model-based clustering techniques are often superior to nonparametric approaches  <cit> .
large-scale methylation arrays are now available for studying methylation genome-wide; the goldengate methylation platform from the manufacturer illumina  simultaneously measures cytosine methylation at  <dig> phenotypically-important loci associated with over  <dig> cancer-related genes.
the result of the array is a sequence of "beta" values, one for each locus, calculated as the average of approximately  <dig> replicates  of the quantity max/, where u is the fluorescent signal from an unmethylated allele on a single bead, m is that from a methylated allele, and q is a constant chosen to ensure that the quantity is well-defined; an absolute value is used in the denominator of the formula to compensate for negative signals due to background subtraction.
note that each beta value is an approximately continuous variable lying between zero and one, where zero represents an unmethylated locus and one represents a methylated locus.
as such, the beta value is appropriately modeled with a beta distribution; we further motivate the choice of beta distribution in the methods section below.
a data set consisting of such sequences produces a high-dimensional data-analysis problem which poses challenges for traditional clustering approaches.
in addition, analysis of heterogeneous tissue data can lead to a large number of clusters, as we demonstrate below, which presents further challenges for clustering techniques.
for example, nonparametric approaches rely on a choice of metric, which may be difficult to justify in the context of high dimensions and numerous clusters.
on the other hand, in model-based clustering, multi-modality of the data likelihood may lead to numerical instability or difficulty in determining the best solution  <cit> .
we propose a novel method for model-based clustering of data of the type produced by illumina goldengate arrays.
our method makes use of a beta mixture model  <cit> .
although one could use bic  to select the number of clusters in the data set, we propose a recursive-partitioning algorithm that provides the number of clusters and a reliable solution in a shorter amount of time than sequential attempts with different numbers of assumed clusters.
this is similar in spirit to the idea of recursive partitioning used in hierarchical ordered partitioning and collapsing hybrid , in which clusters are recursively partitioned using a nonparametric algorithm such as pam  <cit> .
our method is also an unsupervised variant of hierarchical mixtures of experts  <cit> , a fuzzy version of cart  <cit> .
we also propose a method for reducing the number of loci considered in the analysis, and selecting the optimal number using an "augmented" bic statistic.
we also present a simulation study comparing its properties to those of competitor methods.
finally, we demonstrate the methodology on goldengate methylation array data obtained from  <dig> normal tissue samples.
RESULTS
throughout the text, we assume a beta mixture model.
for subject i, let yi =  be a vector of j continuous outcomes falling between  <dig> and  <dig>  and let there be n such vectors.
we posit a mixture model having k classes, so that subject i belongs to class ci ∈ { <dig>  ..., k} and conditional on class membership, each outcome yij is an independent beta-distributed variable with parameters αkj and βkj depending on both class k and dimension j. the objective of our method is to obtain posterior class membership probabilities wik = p. each class has a profile that is obtained from the collection of parameters αkj and βkj, or equivalently the first two moments:
e=μjk=αkj−1andvar⁡=μjk− <dig>
simulation
hc = hierarchical clustering
dyntree = hierarchical clustering with classes determined by dynamic tree cutting
hopach = hopach with 'best' number of classes
hopach = hopach with 'greedy' number of classes
mm = beta mixture model fitting 1– <dig> classes sequentially
rpmm  = recursively partitioned mixture model employing icl-bic
rpmm  = recursively partitioned mixture model employing bic
j = number of loci considered in analysis
dyntree = hierarchical clustering with classes determined by dynamic tree cutting
hopach = hopach with 'best' number of classes
hopach = hopach with 'greedy' number of classes
rpmm = recursively partitioned mixture model employing bic
j = number of loci considered in analysis
in the methods section we propose a variant of bic to compare model fit for different numbers j of loci.
for case i, the mixture models always minimized the "augmented bic" at j =  <dig>  while for case ii, the mixture models always minimized the augmented bic at j =  <dig>  for case i, nearly all  <dig> dimensions were at least somewhat informative; it is interesting to note that j was always minimized at its highest value for this case.
for case ii, the number of informative dimensions was j =  <dig>  so the minimum j was closest to the true number of informative markers among the j considered in this simulation.
in additional simulations that used a finer mesh of j, j was minimized at  <dig>  similar results were obtained when the classes were less balanced .
normal tissue
we applied the recursively-partitioned mixture model algorithm to the normal tissue described below in the methods section.
for this analysis, we attempted to split a node only if the weight assigned to the node was greater than  <dig>  the augbicj statistic, which we propose as a comparative measure of model fit for different numbers j of loci, was minimized at j =  <dig>  and the algorithm found  <dig> classes, whose profiles  are depicted in figure  <dig>  all posterior class membership probabilities were indistinguishable from  <dig> or  <dig> within numerical error.
table  <dig> displays the cross-classification between mixture model latent class and tissue sample type.
blood samples were completely separated from other solid tissue samples.
in addition, adult blood samples were completely separated from newborn blood samples obtained from guthrie cards.
placenta samples were also separated from other tissues aside from a single pleura sample.
for the most part, head and neck tissue and brain were separated from other samples, but were poorly distinguished between each other.
these results were consistent with a random forests analysis  <cit>  , in which we found blood perfectly classified, low classification error for placenta, and some confusion among head and neck tissue and brain tissue.
classes are labeled with the sequence vector representing the terminal node from which the class was derived.
using a permutation test with chi-square statistic, the p value for a hypothesis of no association between class and sample type was less than  <dig> .
thus, our proposed method found clusters relevant to sample type.
in addition, a permutation test using a kruskal-wallis test statistic produced p <  <dig>  for a hypothesis of no difference in mean age among the classes.
interestingly, when the clustering and subsequent hypothesis test was restricted to blood, the p <  <dig>  for a hypothesis of no difference in mean age among latent classes.
between the two classes found among the adult liquid blood samples, age was significantly different , consistent with known associations between age and methylation.
note that the recursion path contains information about relative similarity.
for example, in figure  <dig>  the classes  seem rather distinct from the classes .
there is a band of loci that show high methylation in two classes at the top of the figure, indexed as  <dig> and  <dig>  corresponding to brain tissues, but not in the other classes in the left node , mostly corresponding to blood samples; this band also occurs at the bottom of figure  <dig>  in classes such as  <dig>  which also include brain tissues.
however, as in any hierarchical clustering procedure, the ordering of the children within a node is meaningless, so that classes  and  could have their positions swapped.
there is a wider band of loci to the right, which though the distinction may be subtle, appears to drive the distinction between  and .
however, because the initial levels of recursion may constitute only a crude splitting of the data, the tree representing the recursion path likely has greater meaning only at lower levels of recursion.
we also analyzed the normal tissue methylation data using hopach.
the greedy version of the algorithm produced only  <dig> classes.
the "best" version produced  <dig> clusters, which are cross-classified with tissue type in table  <dig> and with the latent classes obtained from our proposed method in table  <dig>  as table  <dig> shows, the classes found by our proposed method were, for the most part, subsets of the  <dig> classes found using hopach, with a few exceptions that involve minor disagreements in classification.
while the apparent compactness of the hopach classification seems, at first glance, more attractive than the classification produced by our model-based method, we remark that is has a few subtle problems.
it has three singletons, clusters  <dig>   <dig>  and  <dig>  which could be grouped together with cluster  <dig> to comprise a class that entirely represents placental tissue.
while a similar criticism could be made of our proposed method with respect to classification of blood, some of the classes have verifiable meaning; for example, the classes indexed  <dig> and  <dig> distinguish age among blood samples taken from adults.
hopach classification also associates one head and neck tissue sample with blood, and two cervical samples with numerous other tissues in the 5th class.
the 5th class associates numerous tissue types, and comprises  <dig> different classes produced by our proposed method.
table  <dig> shows the cross-classification of tissue type with these  <dig> classes for the subjects comprising the 5th hopach class.
the classification correctly isolates two of three cervix samples, and has a tendency to distinguish pleura from lung samples.
using a permutation test with chi-square statistic, the p value for a hypothesis of no association between class and sample type among the  <dig> classes from our model-based method was less than  <dig> , demonstrating that the classes produced by the mixture model have additional information with respect to tissue type.
in order to compare the predictive ability of the two classification schemes overall, we applied the random forest algorithm to indicator variables representing hopach clusters  and to indicator variables representing our model-based classification .
in the former case we obtained a misclassification rate of  <dig> %, and in the latter case a misclassification rate of  <dig> %, the difference being the misclassification of one less sample using the model-based method.
employing the random forest algorithm in a similar manner to predict age, we obtained a mean-squared-residual of  <dig>  for hopach and  <dig>  for the model-based classification, with variance explained equal to  <dig> % and  <dig> % respectively.
thus, the model-based classification seems to offer modest improvements over hopach in ability to make biological distinctions.
rows represent classes from proposed method, labeled with the sequence vector representing the terminal node from which the class was derived.
columns represent clusters from hopach.
classes are labeled with the sequence vector representing the terminal node from which the class was derived.
discussion
our proposed method is a model-based version of the hopach algorithm  <cit> : we recursively use a beta-mixture model  <cit>  to propose a split of an existing cluster, preserving the split only when it is judged on the basis of bic to better fit the data.
we have found that the integrated-classification-likelihood bic   <cit>  and bic produce practically identical results; this result is not unexpected, since the two differ only in a term involving entropy and a perfect classification results in zero entropy.
a high-dimensional data set will tend to result in approximately perfect classification of most subjects, even though each subject has the opportunity to be portioned out over multiple clusters .
consequently, in the deeper levels of recursion, the algorithm will not depend on the distinction between bic and icl-bic.
for the earlier recursion stages, we assume that there is a strong preference for splitting, which overcomes the entropy penalty in the icl-bic.
siegmund et al.
argue that model-based clustering is preferred in this context over hierarchical clustering  <cit> , a finding that bears out in our simulations.
one reason for the superior performance, at least in a high-dimensional context, is that the metric used to characterize the differences in nonparametric contexts may be relatively insensitive to differences in particular dimensions.
this may play a role in the apparent differences in classification of normal tissue between our proposed method and hopach.
k-means have been used recently to cluster methylation outcomes  <cit> , though the work of van der laan and pollard  seems to suggest that hopach may yield results that are superior to k-means.
in particular, with k-means it is difficult to know how many classes are inherent in the data without resampling-based methods such as the gap statistic  <cit> , with implications for scalability.
also, the "curse of dimensionality" would tend to degrade the performance of procedures such as k-means when there are a large number of clusters and the observed data is of high dimension.
in general, nonparametric methods such as the fanny algorithm  <cit>  rely on tuning parameters that are difficult to optimize without resampling.
an additional problem with non-parametric procedures is that they typically consider only the first moment  of the underlying distributions, ignoring the second-moment  which for dna methylation as measured by the goldengate assay, may play a critical role in distinguishing tissues.
we propose a dimension-reduction strategy which simply ranks candidate dimensions on the basis of some criterion such as variance, fits the top j dimensions in a mixture model, and employs an augmented version of bic to compare model fit across different values of j. this is a departure from the penalized-likelihood methods of the kind described in  <cit> , which would become computationally difficult for truly high-dimensional data.
our approach is similar in spirit to supervised principal components methods such as  <cit> .
interestingly, for the normal tissue data, all  <dig> loci were found to be informative.
the implication is that methylation at even the least variable locus, col6a1_p283_f , contains information about tissue type.
in fact, box-plots showing the distribution of col6a1_p283_f methylation demonstrate great heterogeneity in apparent distribution by tissue type , even though all methylation average beta values were less than  <dig> .
the kruskal-wallis p value for testing the association of col6a1_p283_f methylation with tissue type was p <  <dig> .
while this difference could be a subtle artifact of plate-to-plate variability and heterogeneity in plate sample composition, it also may be that the average beta measured by the goldengate assay is in fact an average of methylation status over different cell types: that is, tissue samples may consist of cells having heterogeneous methylation states, with each tissue type having a unique mixture.
a separate paper describing analysis of normal tissue in detail is under preparation.
we remark that we did not normalize the different plates used in laboratory analysis, for reasons described below.
however, when we used the analytical methods of this paper to classify methylation profiles within plate strata and within tissue type strata, we found much greater variability between tissue types than between plates .
our results are consistent with expectations developed from an understanding of the current literature on dna methylation in normal tissues, though we cannot entirely rule out subtle biases introduced by variation between plates.
it may be possible to incorporate plate effects in our proposed methodology by modeling the mean response , but the computational considerations are beyond the scope of the present paper.
CONCLUSIONS
in summary, our method appears to have good properties both with respect to classification error and computation time.
it achieves these properties by combining the strengths of model-based and hierarchical methods, navigating the underlying clusters quickly through recursive partitioning, but doing so in a way that makes use of a reasonable probability model.
this model is also used to compare different dimensions j of input, thus refining the discriminative ability in a scalable manner.

 epigenetics is the study of heritable changes in gene function that cannot be explained by changes in dna sequence. one of the most commonly studied epigenetic alterations is cytosine methylation, which occurs in the context of a cpg dinucleotide. concentrations of cpgs known as cpg islands, when sufficiently methylated, are associated with transcriptional gene silencing tantamount to "one hit" as part of knudson's two hit hypothesis of tumor suppressor gene inactivation. dna methylation associated gene silencing is a well recognized epigenetic mechanism that often occurs at tumor suppressor gene loci in human cancer. hundreds of reports of methylation induced silencing at tumor suppressor genes in virtually all types of human cancer have been published. while there has been a tremendous effort to characterize epigenetic alterations in cancer, surprisingly little work has been done in disease-free tissues. there is a basic need for epigenetic profiling of normal tissues to better understand the contribution of these profiles to tissue specificity, especially in the context of phenotypically important cpgs, where deregulation is associated with human diseases such as cancer. while efforts to characterize the methylation profiles of normal tissues in humans and mice have begun, and certain themes are slowly becoming apparent, relatively few reports have emerged. most cpgs or cpg regions have been found to have a bimodal distribution of methylation profiles, either hypo- or hypermethylated. another theme is the disproportionate location of cell or tissue type dependent differentially methylated regions in non-cpg island sites. furthermore, the human epigenome project showed that the human major histocompatability loci have differential methylation across various human tissues, but that differential methylation does not necessarily lead to differential expression. it is therefore critical to first outline the basal-state of phenotypically important epigenetic marks that are known to contribute to cancer in order to have a background for comparison to other normal and diseased tissues. an approach which characterizes the basal epigenetic state is best suited to foster the discovery of epigenetic profiles that are associated with particular disease states, patient characteristics, exposures or other covariates that may contribute to pathogenesis. cluster analysis is often used to identify methylation subgroups in data and, in particular, siegmund et al. argue that model-based clustering techniques are often superior to nonparametric approaches. large-scale methylation arrays are now available for studying methylation genome-wide; the goldengate methylation platform from the manufacturer illumina simultaneously measures cytosine methylation at phenotypically-important loci associated with over cancer-related genes. the result of the array is a sequence of "beta" values, one for each locus, calculated as the average of approximately replicates of the quantity max/, where u is the fluorescent signal from an unmethylated allele on a single bead, m is that from a methylated allele, and q is a constant chosen to ensure that the quantity is well-defined; an absolute value is used in the denominator of the formula to compensate for negative signals due to background subtraction. note that each beta value is an approximately continuous variable lying between zero and one, where zero represents an unmethylated locus and one represents a methylated locus. as such, the beta value is appropriately modeled with a beta distribution; we further motivate the choice of beta distribution in the methods section below. a data set consisting of such sequences produces a high-dimensional data-analysis problem which poses challenges for traditional clustering approaches. in addition, analysis of heterogeneous tissue data can lead to a large number of clusters, as we demonstrate below, which presents further challenges for clustering techniques. for example, nonparametric approaches rely on a choice of metric, which may be difficult to justify in the context of high dimensions and numerous clusters. on the other hand, in model-based clustering, multi-modality of the data likelihood may lead to numerical instability or difficulty in determining the best solution. we propose a novel method for model-based clustering of data of the type produced by illumina goldengate arrays. our method makes use of a beta mixture model. although one could use bic to select the number of clusters in the data set, we propose a recursive-partitioning algorithm that provides the number of clusters and a reliable solution in a shorter amount of time than sequential attempts with different numbers of assumed clusters. this is similar in spirit to the idea of recursive partitioning used in hierarchical ordered partitioning and collapsing hybrid, in which clusters are recursively partitioned using a nonparametric algorithm such as pam. our method is also an unsupervised variant of hierarchical mixtures of experts, a fuzzy version of cart. we also propose a method for reducing the number of loci considered in the analysis, and selecting the optimal number using an "augmented" bic statistic. we also present a simulation study comparing its properties to those of competitor methods. finally, we demonstrate the methodology on goldengate methylation array data obtained from normal tissue samples.  throughout the text, we assume a beta mixture model. for subject i, let yi = be a vector of j continuous outcomes falling between and and let there be n such vectors. we posit a mixture model having k classes, so that subject i belongs to class ci ∈ {..., k} and conditional on class membership, each outcome yij is an independent beta-distributed variable with parameters αkj and βkj depending on both class k and dimension j. the objective of our method is to obtain posterior class membership probabilities wik = p. each class has a profile that is obtained from the collection of parameters αkj and βkj, or equivalently the first two moments: e=μjk=αkj−1andvar⁡=μjk− simulation hc = hierarchical clustering dyntree = hierarchical clustering with classes determined by dynamic tree cutting hopach = hopach with'best' number of classes hopach = hopach with'greedy' number of classes mm = beta mixture model fitting 1– classes sequentially rpmm = recursively partitioned mixture model employing icl-bic rpmm = recursively partitioned mixture model employing bic j = number of loci considered in analysis dyntree = hierarchical clustering with classes determined by dynamic tree cutting hopach = hopach with'best' number of classes hopach = hopach with'greedy' number of classes rpmm = recursively partitioned mixture model employing bic j = number of loci considered in analysis in the methods section we propose a variant of bic to compare model fit for different numbers j of loci. for case i, the mixture models always minimized the "augmented bic" at j = while for case ii, the mixture models always minimized the augmented bic at j = for case i, nearly all dimensions were at least somewhat informative; it is interesting to note that j was always minimized at its highest value for this case. for case ii, the number of informative dimensions was j = so the minimum j was closest to the true number of informative markers among the j considered in this simulation. in additional simulations that used a finer mesh of j, j was minimized at similar results were obtained when the classes were less balanced. normal tissue we applied the recursively-partitioned mixture model algorithm to the normal tissue described below in the methods section. for this analysis, we attempted to split a node only if the weight assigned to the node was greater than the augbicj statistic, which we propose as a comparative measure of model fit for different numbers j of loci, was minimized at j = and the algorithm found classes, whose profiles are depicted in figure all posterior class membership probabilities were indistinguishable from or within numerical error. table displays the cross-classification between mixture model latent class and tissue sample type. blood samples were completely separated from other solid tissue samples. in addition, adult blood samples were completely separated from newborn blood samples obtained from guthrie cards. placenta samples were also separated from other tissues aside from a single pleura sample. for the most part, head and neck tissue and brain were separated from other samples, but were poorly distinguished between each other. these results were consistent with a random forests analysis, in which we found blood perfectly classified, low classification error for placenta, and some confusion among head and neck tissue and brain tissue. classes are labeled with the sequence vector representing the terminal node from which the class was derived. using a permutation test with chi-square statistic, the p value for a hypothesis of no association between class and sample type was less than. thus, our proposed method found clusters relevant to sample type. in addition, a permutation test using a kruskal-wallis test statistic produced p < for a hypothesis of no difference in mean age among the classes. interestingly, when the clustering and subsequent hypothesis test was restricted to blood, the p < for a hypothesis of no difference in mean age among latent classes. between the two classes found among the adult liquid blood samples, age was significantly different, consistent with known associations between age and methylation. note that the recursion path contains information about relative similarity. for example, in figure the classes seem rather distinct from the classes. there is a band of loci that show high methylation in two classes at the top of the figure, indexed as and corresponding to brain tissues, but not in the other classes in the left node, mostly corresponding to blood samples; this band also occurs at the bottom of figure in classes such as which also include brain tissues. however, as in any hierarchical clustering procedure, the ordering of the children within a node is meaningless, so that classes and could have their positions swapped. there is a wider band of loci to the right, which though the distinction may be subtle, appears to drive the distinction between and. however, because the initial levels of recursion may constitute only a crude splitting of the data, the tree representing the recursion path likely has greater meaning only at lower levels of recursion. we also analyzed the normal tissue methylation data using hopach. the greedy version of the algorithm produced only classes. the "best" version produced clusters, which are cross-classified with tissue type in table and with the latent classes obtained from our proposed method in table as table shows, the classes found by our proposed method were, for the most part, subsets of the classes found using hopach, with a few exceptions that involve minor disagreements in classification. while the apparent compactness of the hopach classification seems, at first glance, more attractive than the classification produced by our model-based method, we remark that is has a few subtle problems. it has three singletons, clusters and which could be grouped together with cluster to comprise a class that entirely represents placental tissue. while a similar criticism could be made of our proposed method with respect to classification of blood, some of the classes have verifiable meaning; for example, the classes indexed and distinguish age among blood samples taken from adults. hopach classification also associates one head and neck tissue sample with blood, and two cervical samples with numerous other tissues in the 5th class. the 5th class associates numerous tissue types, and comprises different classes produced by our proposed method. table shows the cross-classification of tissue type with these classes for the subjects comprising the 5th hopach class. the classification correctly isolates two of three cervix samples, and has a tendency to distinguish pleura from lung samples. using a permutation test with chi-square statistic, the p value for a hypothesis of no association between class and sample type among the classes from our model-based method was less than, demonstrating that the classes produced by the mixture model have additional information with respect to tissue type. in order to compare the predictive ability of the two classification schemes overall, we applied the random forest algorithm to indicator variables representing hopach clusters and to indicator variables representing our model-based classification. in the former case we obtained a misclassification rate of %, and in the latter case a misclassification rate of %, the difference being the misclassification of one less sample using the model-based method. employing the random forest algorithm in a similar manner to predict age, we obtained a mean-squared-residual of for hopach and for the model-based classification, with variance explained equal to % and % respectively. thus, the model-based classification seems to offer modest improvements over hopach in ability to make biological distinctions. rows represent classes from proposed method, labeled with the sequence vector representing the terminal node from which the class was derived. columns represent clusters from hopach. classes are labeled with the sequence vector representing the terminal node from which the class was derived. discussion our proposed method is a model-based version of the hopach algorithm: we recursively use a beta-mixture model to propose a split of an existing cluster, preserving the split only when it is judged on the basis of bic to better fit the data. we have found that the integrated-classification-likelihood bic and bic produce practically identical results; this result is not unexpected, since the two differ only in a term involving entropy and a perfect classification results in zero entropy. a high-dimensional data set will tend to result in approximately perfect classification of most subjects, even though each subject has the opportunity to be portioned out over multiple clusters. consequently, in the deeper levels of recursion, the algorithm will not depend on the distinction between bic and icl-bic. for the earlier recursion stages, we assume that there is a strong preference for splitting, which overcomes the entropy penalty in the icl-bic. siegmund et al. argue that model-based clustering is preferred in this context over hierarchical clustering, a finding that bears out in our simulations. one reason for the superior performance, at least in a high-dimensional context, is that the metric used to characterize the differences in nonparametric contexts may be relatively insensitive to differences in particular dimensions. this may play a role in the apparent differences in classification of normal tissue between our proposed method and hopach. k-means have been used recently to cluster methylation outcomes, though the work of van der laan and pollard seems to suggest that hopach may yield results that are superior to k-means. in particular, with k-means it is difficult to know how many classes are inherent in the data without resampling-based methods such as the gap statistic, with implications for scalability. also, the "curse of dimensionality" would tend to degrade the performance of procedures such as k-means when there are a large number of clusters and the observed data is of high dimension. in general, nonparametric methods such as the fanny algorithm rely on tuning parameters that are difficult to optimize without resampling. an additional problem with non-parametric procedures is that they typically consider only the first moment of the underlying distributions, ignoring the second-moment which for dna methylation as measured by the goldengate assay, may play a critical role in distinguishing tissues. we propose a dimension-reduction strategy which simply ranks candidate dimensions on the basis of some criterion such as variance, fits the top j dimensions in a mixture model, and employs an augmented version of bic to compare model fit across different values of j. this is a departure from the penalized-likelihood methods of the kind described in, which would become computationally difficult for truly high-dimensional data. our approach is similar in spirit to supervised principal components methods such as. interestingly, for the normal tissue data, all loci were found to be informative. the implication is that methylation at even the least variable locus, col6a1_p283_f, contains information about tissue type. in fact, box-plots showing the distribution of col6a1_p283_f methylation demonstrate great heterogeneity in apparent distribution by tissue type, even though all methylation average beta values were less than. the kruskal-wallis p value for testing the association of col6a1_p283_f methylation with tissue type was p <. while this difference could be a subtle artifact of plate-to-plate variability and heterogeneity in plate sample composition, it also may be that the average beta measured by the goldengate assay is in fact an average of methylation status over different cell types: that is, tissue samples may consist of cells having heterogeneous methylation states, with each tissue type having a unique mixture. a separate paper describing analysis of normal tissue in detail is under preparation. we remark that we did not normalize the different plates used in laboratory analysis, for reasons described below. however, when we used the analytical methods of this paper to classify methylation profiles within plate strata and within tissue type strata, we found much greater variability between tissue types than between plates. our results are consistent with expectations developed from an understanding of the current literature on dna methylation in normal tissues, though we cannot entirely rule out subtle biases introduced by variation between plates. it may be possible to incorporate plate effects in our proposed methodology by modeling the mean response, but the computational considerations are beyond the scope of the present paper.  in summary, our method appears to have good properties both with respect to classification error and computation time. it achieves these properties by combining the strengths of model-based and hierarchical methods, navigating the underlying clusters quickly through recursive partitioning, but doing so in a way that makes use of a reasonable probability model. this model is also used to compare different dimensions j of input, thus refining the discriminative ability in a scalable manner. 
epigenetics is the study of heritable changes in gene function that cannot be explained by changes in dna sequence. we propose a novel model-based recursive-partitioning algorithm to navigate clusters in a beta mixture model. we present simulations that show that the method is more reliable than competing nonparametric clustering approaches, and is at least as reliable as conventional mixture model methods. we also show that our proposed method is more computationally efficient than conventional mixture model approaches. our proposed recursively-partitioned mixture model is an effective and computationally efficient method for clustering dna methylation data. one of the most commonly studied epigenetic alterations is cytosine methylation, which is a well recognized mechanism of epigenetic gene silencing and often occurs at tumor suppressor gene loci in human cancer. we demonstrate our method on the normal tissue samples and show that the clusters are associated with tissue type as well as age. arrays are now being used to study dna methylation at a large number of loci; for example, the illumina goldengate platform assesses dna methylation at loci associated with over cancer-related genes. model-based cluster analysis is often used to identify dna methylation subgroups in data, but it is unclear how to cluster dna methylation data from arrays in a scalable and reliable manner.
 epigenetics is the study of heritable changes in gene function that cannot be explained by changes in dna sequence. one of the most commonly studied epigenetic alterations is cytosine methylation, which is a well recognized mechanism of epigenetic gene silencing and often occurs at tumor suppressor gene loci in human cancer. arrays are now being used to study dna methylation at a large number of loci; for example, the illumina goldengate platform assesses dna methylation at loci associated with over cancer-related genes. model-based cluster analysis is often used to identify dna methylation subgroups in data, but it is unclear how to cluster dna methylation data from arrays in a scalable and reliable manner.  we propose a novel model-based recursive-partitioning algorithm to navigate clusters in a beta mixture model. we present simulations that show that the method is more reliable than competing nonparametric clustering approaches, and is at least as reliable as conventional mixture model methods. we also show that our proposed method is more computationally efficient than conventional mixture model approaches. we demonstrate our method on the normal tissue samples and show that the clusters are associated with tissue type as well as age.  our proposed recursively-partitioned mixture model is an effective and computationally efficient method for clustering dna methylation data. 