The task of clustering a set of objects based on multiple sources of data arises in several modern applications. We propose an integrative statistical model that permits a separate clustering of the objects for each data source. These separate clusterings adhere loosely to an overall consensus clustering, and hence they are not independent. We describe a computationally scalable Bayesian framework for simultaneous estimation of both the consensus clustering and the source-specific clusterings. We demonstrate that this flexible approach is more robust than joint clustering of all data sources, and is more powerful than clustering each data source separately. This work is motivated by the integrated analysis of heterogeneous biomedical data, and we present an application to subtype identification of breast cancer tumor samples using publicly available data from The Cancer Genome Atlas. Software is available at http://people.duke.edu/~el113/ software.html. Several fields of research now analyze multi-source data (also called multimodal data), in which multiple heterogeneous datasets describe a common set of objects. Each dataset represents a distinct mode of measurement or domain. Table 1 gives examples of multi-source data from very diverse research areas. While the methodology described in this article is broadly applicable, our primary motivation is the integrated analysis of heterogeneous biomedical data. The diversity of platforms and technologies that are used to collect genomic data, in particular, is expanding rapidly. Often multiple types of genomic data, measuring various biological components, are collected for a commen set of samples. For example, The Cancer Genome Atlas (TCGA) is a large-scale collaborative effort to collect and catalog data from several genomic technologies. The integrative analysis of data from these disparate sources provides a more comprehensive understanding of cancer genetics and molecular biology. In Section 6 we present an analysis of RNA expression, DNA methylation, microRNA expression, and proteomic data from TCGA for a common set of breast cancer tumor samples. Separate analyses of each data source may lack power and will not capture inter-source associations. At the other extreme, a joint analysis that ignores the heterogeneity of the data may not capture important features that are specific to each data source. Exploratory methods that simultaneously model shared features and features that are specific to each data source have recently been developed as flexible alternatives [10, 9, 20, 16] . The demand for such integrative methods motivates a very dynamic area of statistics and machine learning. This article concerns integrative clustering. Clustering is a very widely used exploratory tool to identify similar groups of objects (for example, clinically relevant disease subtypes). Hundreds of general algorithms to perform clustering have been proposed in the literature. However, our work is motivated by the need for an integrative clustering method that is computationally scalable and robust to the unique features of each data source. 2 2 Related Work Most applications of clustering multi-source data follow one of two general approaches: 1. Clustering of each data source separately, potentially followed by a post hoc integration of these separate clusterings. 2. Combining all data sources to determine a single "joint" clustering. Under approach (1) the level of agreement between the separate clusterings may be measured by the adjusted rand index [6] or a similar statistic. Furthermore, consensus clustering can be used to determine an overall partition of the objects that agrees the most with the data-specific clusterings. Several objective functions and algorithms to perform consensus clustering have been proposed (for a survey see Nguyen & Caruana [13]). These methods are most commonly used to establish consensus among multiple clustering algorithms, or multiple realizations of the same clustering algorithm, on a single dataset. Consensus clustering has also been used to integrate multi-source data [4, 1] . Such an approach is attractive in that it models source-specific features yet still determines an overall clustering, which is often of practical interest. However, the two stage process of performing entirely separate clusterings followed by post hoc integration limits the power to identify and exploit shared structure (see Section 5.2 for an illustration of this phenomenon). Approach (2) effectively exploits shared structure, at the expense of failing to recognize features that are specific to each data source. Within a modelbased statistical framework, one can find the clustering that maximizes a joint likelihood. Assuming that each source is conditionally independent given the clustering, the joint likelihood is the product of the likelihood functions for each data source. This approach is used by Kormaksson et al. [8] in the context of integrating gene expression and DNA methylation data. The iCluster method [18, 11] performs clustering by first fitting a Gaussian latent factor model to the joint likelihood; clusters are then determined by K-means clustering of the factor scores. Rey & Roth [17] propose a dependency-seeking model in which the goal is to find a clustering that accounts for associations across the data sources. Perhaps most similar to our approach in spirit and motivation is the Multiple Dataset Integration (MDI) method [7], which uses a statistical framework to cluster each data source separately while simultaneously modeling dependence between the clusterings. By explicitly modeling dependence, 3 MDI permits borrowing strength across data sources. The fundamental difference between MDI and the approach we describe in the following concerns how dependence is modeled. Specifically, MDI models the dependence between each pair of data sources rather than adherence to an overall clustering. We elaborate on this distinction in Sections 3.1 and 5.2. Here we briefly describe the finite Dirichlet mixture model for clustering a single dataset, with the purpose of laying the groundwork for the integrative model given in Section 3. Given data X n for N objects (n = 1, ..., N ), the goal is to partition these objects into at most K clusters. Typically X n is a multi-dimensional vector, but we present the model in sufficient generality to allow for more complex data structures. Let f (X n |θ) define a probability model for X n given parameter(s) θ. For example, f may be a Gaussian density defined by the mean and variance θ = (µ, σ 2 ). Each X n is drawn independently from a mixture distribution with K components, specified by the parameters θ 1 , . . . , θ K . Let C n ∈ {1, . . . , K} represent the component corresponding to X n , and π k be the probability that an arbitrary object belongs to cluster k: Then, the generative model is We further assume that Π = (π 1 , . . . , π K ) has a Dirichlet distribution parameterized by a K-dimensional vector β of positive reals. This allows some π k to be small, and therefore N objects may not represent all K clusters. Letting K → ∞ gives a Dirichlet process. Under a Bayesian framework one can put a prior distribution on Π and the parameter set Θ = (θ 1 , . . . , θ K ). Standard computational methods such as Gibbs sampling can then be used to approximate the posterior distribution for Π, Θ, and C = (C 1 , . . . , C N ) (for an overview see Neal [12]). We extend the Dirichlet mixture model to accommodate data from M sources X 1 , . . . , X M . Each data source is available for a common set of N objects, where X mn represents data m for object n. Each data source requires a probability model f m (X n |θ m ) parametrized by θ m . Under the general framework 4 presented here each X m may have disparate structure. For example X 1n may give an image where f 1 defines the spectral density for a Gaussian random field, while X 2n may give a categorical vector where f 2 defines a multivariate probability mass function. We assume there is a separate clustering of the objects for each data source, but that these adhere loosely to an overall clustering. Formally, each X mn n = 1, . . . , N is drawn independently from a K-component mixture distribution specified by the parameters θ m1 , . . . , θ mK . Let L mn ∈ {1, . . . , K} represent the component corresponding to X mn . Furthermore, let C n ∈ {1, . . . , K} represent the overall mixture component for object n. The sourcespecific clusterings L m = (L m1 , . . . , L mN ) are dependent on the overall clustering C = (C 1 , . . . , C N ): where α m adjusts the dependence function ν. The data X m are independent of C conditional on the source-specific clustering L m . Hence, C serves only to unify Throughout this article, we assume ν has the simple form where , 1] controls the adherence of data source m to the overall clustering. More simply α m is the probability that L mn = C n . So, if α m = 1 then L m = C. The α m are estimated from the data together with C and L 1 , . . . , L m . In practice we estimate each α m separately, or assume that α 1 = . . . = α M and hence each data source adheres equally to the overall clustering. The latter is favored when m = 2 for identifiability reasons. More complex models that permit dependence of the α m 's are also potentially useful. Let π k be the probability that an object belongs to the overall cluster k: We assume Π = (π 1 , . . . , π K ) follows a Dirichlet(β) distribution. The probability that an object belongs to a given source-specific cluster follows directly: 5 Moreover, a simple application of Bayes rule gives the conditional distribution of C: where ν is defined as in (1). Although the number of possible clusters K is the same for L 1 , . . . , L M and C, the number of clusters that are actually represented may vary. Generally the source-specific clusterings L m will represent more clusters than C, rather than vice-versa. This follows from Equation (2) and is illustrated in AppendixrefSec2. Intuitively if object n is not allocated to any overall cluster in data source m (i.e., L mn / ∈ C) then X mn does not conform well to any overall pattern in the data. Data source m X mn Data for object n, source m f m Probability model for source m θ mk Parameters for f m , cluster k p m Prior distribution for θ mk C n Overall cluster for object n π k Probability that C n = k L mn Cluster specific to X mn ν Dependence function for C n and L mn α m Probability that L mn = C n  Integrating over the overall clustering C gives the joint marginal distribution of L 1 , . . . , L M : 6 Under the assumption that α 1 = . . . = α M the model simplifies: where t k is the number of clusters equal to k and U = (K−1)α 1 1−α 1 ≥ 1. This marginal form facilitates comparison with the MDI method for dependent clustering. In the MDI model φ ij > 0 control the strength of association between the clusterings L i and L j : whereπ mk = P (L mn = k). For M = 2 andπ 1· =π 2· it is straightforward to show that (4) and (5) are functionally equivalent under a parameter substitution (see Appendix C). There is no such equivalence for M > 2, regardless of restrictions onΠ and Φ. This is not surprising, as MDI gives a general model of pairwise dependence between clusterings rather than a model of adherence to an overall clustering. Here we present a general Bayesian framework for estimation of the integrative clustering model. We employ a Gibbs sampling procedure to approximate the posterior distribution for the parameters introduced in Section 3. The algorithm is general in that we do not assume any specific form for the f m and the parameters θ mk . We use conjugate prior distributions for α m , Π, and (if possible) θ mk . • . By default we choose a m = b m = 1, so that the prior for α m is uniformly distributed between 1 K and 1. • Π ∼ Dirichlet(β 0 ). By default we choose β 0 = (1, 1, ..., 1), so that the prior for Π is uniformly distributed on the standard (M − 1)-simplex. • The θ mk have prior distribution p m . In practice, one should choose p m so that sampling from the conditional posterior Markov chain Monte Carlo (MCMC) proceeds by iteratively sampling from the following conditional posterior distributions: • L m |X m , Θ m , α m , C ∼ P (k|X mn , C n , θ mk , α m ) for n = 1, . . . , N , where ), where τ m is the number of samples n satisfying L mn = C n . • C|L m , Π, α ∼ P (k|Π, {L mn , α m } M m=1 ) for n = 1, . . . , N , where • Π|C ∼ Dirichlet(β 0 + ρ), where ρ k is the number of samples allocated to cluster k in C. This algorithm can be suitably modified under the assumption that α 1 = . . . = α M (see Appendix A.2). Each sampling iteration produces a different realization of the clusterings C, L 1 , · · · , L m , and together these samples approximate the posterior distribution for the overall and source-specific clusterings. However, a point estimate may be desired for each of C, L 1 , · · · , L m to facilitate interpretation of the clusters. In this respect methods that aggregate over the MCMC iterations to produce a single clustering, such as that described in [3], can be used. It is possible to derive a similar sampling procedure using only the marginal form for the source-specific clusterings given in Equation (3). However, the overall clustering C is also of interest in most applications. Furthermore, incorporating C into the algorithm can actually improve computational efficiency dramatically, especially if M is large. As presented, each MCMC iteration can be completed in O(M N K) operations. If the full joint marginal distribution of L 1 , . . . , L M is used the computational burden increases exponentially with M (this also presents a bottleneck for the MDI method). For each iteration, C n is determined randomly from a distribution that gives higher probability to clusters that are prevalent in {L 1n , . . . , L mn }. In this sense C is determined by a random consensus clustering of the sourcespecific clusterings. Hence, we refer to this approach as Bayesian consensus clustering (BCC). BCC differs from traditional consensus clustering in three key aspects. 8 1. Both the source-specific clusterings and the consensus clustering are modeled in a statistical way that allows for uncertainty in all parameters. 2. The source-specific clusterings and the consensus clustering are estimated simultaneously, rather than in two stages. This permits borrowing of information across sources for more accurate cluster assignments. 3. The strength of association to the consensus clustering for each data source is learned from the data and accounted for in the model. We have developed software for the R environment for statistical computing [15] to perform BCC on multivariate continuous data sources (available at http://people.duke.edu/~el113/software.html). Specifically, we use a Normal-Gamma conjugate prior distribution to model cluster-specific means and variances. See Appendix A for more details. In theory the specified maximum number of clusters K can be large, for example K = N . The number of clusters realized in C and the L m may still be small. However, we find that this is not the case for high-dimensional structured data such as that used for the genomics application in Section 6. The model tends to select a large number of clusters even if the Dirichlet prior concentration parameters β 0 are very small. From an exploratory perspective we would like to identify a small number of interpretable clusters. For this reason Dirichlet process models that allow for an infinite number of clusters are also not appropriate. We therefore consider an alternative heuristic measure that selects the value of K that gives maximum adherence to an overall clustering. The adjusted adherence parameters α * m ∈ [0, 1] given by are computed for each candidate value of K. We then select the K giving the highest mean adjusted adherencē In practice, we find that this method selects a small number of clusters that reveal shared structure across the data sources. 9 Here we present applications of BCC to simple simulated datasets. To the best of our knowledge there is no existing method that is directly comparable to BCC in terms of the motivating model. Nevertheless, we do illustrate the flexibility and potential advantages of BCC over other model-based clustering methods for multi-source data. We find that with reasonable signal the α m can generally be estimated with accuracy and without substantial bias. To illustrate, we generate simulated datasets X 1 : 1 × 200 and X 2 : 1 × 200 as follows: 1. Let C define two clusters, where C n = 1 for n ∈ {1, . . . , 100} and C n = 2 for n ∈ {101, . . . , 200} 2. Draw α from a Uniform(0.5, 1) distribution. 3. For m = 1, 2 and n = 1, . . . , 200, generate L mn ∈ {1, 2} with probabilities 4. For m = 1, 2 draw values X mn from a Normal(1.5, 1) distribution if L mn = 1 and from a Normal(−1.5, 1) distribution if L mn = 2 We generate 100 realizations of the above simulation, and estimate the model via BCC for each realization. We assume α 1 = α 2 in our estimation and use a uniform prior; further computational details are given in Section 4 of the supplemental article. Figure 1 displaysα, the best estimate for both α 1 and α 2 , versus the true α for each realization. The point estimate displayed is the mean over MCMC draws, and we also display a 95% credible interval based on the 2.5 of 97.5 percentiles of the MCMC draws. The estimatedα are generally close to α, and the credible interval contains the true value in 91 of 100 simulations. See Section 4 of the supplemental article for a more detailed study, including a simulation illustrating the effect of the prior distribution onα. To illustrate the flexibility and advantages of BCC in terms of clustering accuracy, we generate simulated data sources X 1 and X 2 as in Section 5.1 but with Normal(1,1) and Normal(-1,1) as our mixture distributions. Hence, the signal distinguishing the two clusters is weak enough so that there is substantial overlap within each data source. We generate 100 simulations and compare the results for four model-based clustering approaches: 1. Separate clustering, in which a finite Dirichlet mixture model is used to determine a clustering separately for X 1 and X 2 . 2. Joint clustering, in which a finite Dirichlet mixture model is used to determine a single clustering for the concatenated data X 1 X 2 . 3. Dependent clustering, in which we model the pairwise dependence between each data source, in the spirit of MDI. The full implementation details for each method are given in Section 5 of the supplemental article. We consider the relative error for each model in terms of the average number of incorrect cluster assignments: where 1 is the indicator function. Note that for joint clusteringL 1 =L 2 . The relative error for each clustering method on M = 3 datasets is shown in Figure 2 (top panel). Smooth curves are fit to the data using LOESS [2] and display the relative clustering error for each method as a function of α. Not surprisingly, joint clustering performs well for α ≈ 1 (perfect agreement) and separate clustering is superior when α ≈ 0.5 (no relationship). BCC and dependent clustering learn the level of cluster agreement, and hence serve as a flexible bridge between these two extremes. The bottom panel of Figure 2 displays the results for a similar simulation in which M = 3 data sources are generated. Again, BCC is competitive with separate clustering when α ≈ 0.5 and joint clustering when α ≈ 1. Dependent clustering does not perform as well in this example. The pairwisedependence model does not assume an overall clustering and hence it has less power to learn the underlying structure for M > 2. We apply BCC to multi-source genomic data on breast cancer tumor samples from TCGA. For a common set of 348 tumor samples, our full dataset includes • RNA expression (GE) data for 645 genes. • DNA methylation (ME) data for 574 probes. • miRNA expression (miRNA) data for 423 miRNAs. • Reverse Phase Protein Array (RPPA) data for 171 proteins. These four data sources are measured on different platforms and represent different biological components. However, they all represent genomic data for the same samples and it is reasonable to expect some shared structure. These data are publicly available from the TCGA Data Portal. See http: //people.duke.edu/~el113/software.html for R code to completely reproduce the analysis, including instructions on how to download and process these data from the TCGA Data Portal. Previously, four comprehensive sample subtypes were identified based on a multi-source consensus clustering of these data [1] . These subtypes were shown to be clinically relevant, as they may be used for more targeted ther-13 apies and prognoses. We select K = 3 clusters for BCC, based on the heuristic described in Section 4.1. Full computational details, as well as charts that illustrate mixing over the MCMC draws, are available in Appendix F. Table 3 shows a matching matrix comparing the overall clustering C with the comprehensive subtypes defined by TCGA. The two partitions have a significant but weak association, and this suggests they may not be driven by the same biological signal. Similar tables that compare the source-specific clusterings L mn with source-specific subtypes are available in Appendix F.3. TCGA subtype 1 13 6 20 2 66 2 4 3 3 80 78 4 0 3 73 Figure 3 provides a point-cloud view of each dataset given by a scatter plot of the the first two principal components. The overall and source-specific cluster index is shown for each sample, as well as a point estimate for the adherence parameter α (alternatively, clustering heatmaps for each data source are given in Appendix F.4). The four data sources show different structure, and the source-specific clusters are more well distringuished than the overall clusters in each plot. However, the overall clusters are clearly represented to some degree in all four plots. Hence, the flexible yet integrative approach of BCC seems justified for these data. This work was motivated by the perceived need for a general, flexible, and computationally scalable approach to clustering multi-source data. We propose BCC, which models both an overall clustering and a clustering specific to each data source. We view BCC as a form of consensus clustering, with advantages over traditional methods in terms of modeling uncertainty and the ability to borrow information across sources. The BCC model assumes a very simple and general dependence between data sources. When an overall clustering is not sought, or when such a clustering does not make sense as an assumption, a more general model of Figure 3 : PCA plots for each data source. Sample points are colored by overall cluster; cluster 1 is black, cluster 2 is red, and cluster 3 is blue. Symbols indicate source-specific cluster; cluster 1 is '•', cluster 2 is '+', and cluster 3 is ' * ' . cluster dependence (such as MDI) may be more appropriate. Furthermore, a context-specific approach may be necessary when more is known about the underlying dependence of the data. For example, [14] exploit functional covariance models for time-course data to determine overall and time-specific clusters. Our implementation of BCC assumes the data are normally distributed and models cluster-specific mean and variance parameters. It is straightforward to extend this approach to more complex clustering models. In particular, models that assume clusters exist on a sparse feature set [19] or allow for more general covariance structure [5] are growing in popularity. While we focus on multi-source biomedical data, the applications of BCC are potentially widespread. In addition to multi-source data, BCC may be used to compare clusterings from different statistical models for a single homogeneous dataset. 15
REF proposes a Bayesian consensus clustering method.
Bayesian Consensus Clustering
Bioinformatics 29 (2013) 2610-2616
Computer Science Medicine Mathematics