Bayesian biclustering of gene expression data

Simulation results

Bayesian biclustering for yeast datasets

We analyzed the same yeast expression data as in [15] using the BBC procedure. This dataset was derived by combining the environmental stress data of [16] and the cell cycle data of [17]. The combined dataset contains 6108 genes and 250 conditions (or samples). We applied the 90% IQRN procedure across all conditions. Since the dataset contains many missing data, we imputed them along with our BBC iterations. The BBC algorithm was asked to search for K biclusters, with K ranging from 30 to 65. We observed that the BIC [18] achieved the optimal value with K = 57. Out of 6108 genes, 6021 were included in one of the clusters, and all conditions were included in at least one cluster.

We analyzed the clustering results from three aspects. First, we identified the significant categories of experimental conditions for each cluster. More precisely, we classified the 250 experimental conditions into 22 categories according to the biological nature of each experiment. Some examples of categories are heat shock stress, amino acid starvation, and α factor synchronization. Then we searched for significant enrichment for each category. Second, we did functional enrichment test of genes in each cluster using functional information from the MIPS database. Third, we searched the promoter sequences (up to 800 bps upstream) of genes in each bicluster for the enrichment of transcription factor binding sites (TFBS). We applied the TFBS scores of 51 known yeast transcription factors used in [15], which measures how likely a promoter sequence contains a TFBS and was computed using ScanACE [19]. A cutoff value of 0.5 was used to make a presence/absence call for a TFBS. The presence frequency of each TFBS in all genes in the dataset was used as null hypothesis. For all the three types of enrichment analysis, we used the criterion of P value< 0.05 after the Bonferroni correction.

Out of the 57 clusters, 36 have significant gene functions enrichments, 26 have significant TFBS enrichments, 51 have significant experimental condition categories enrichments, 22 have all three types of enrichments and 54 have at least one type of enrichment. We named a few of these biclusters and listed them in Table 4. The enriched gene functions and the significant experimental conditions showed strong correlations. For example, in the cell cycle cluster, the enriched gene function terms are cell cycle, DNA processing and cell division, and the significant experimental conditions are all cell-cycle experiments including cln3 and clb2 experiments, α factor, cdc15, cdc28 and elutriation synchronization. Another example is stress response and protein folding cluster, the enriched gene function terms are stress response and protein folding, and the significant conditions are heat shock, diamide and osmolarity stress experiments. The biclustering results are also supported by the TFBS information. In many clusters, the enriched TFBSs correspond to TFs known to be involved in the pathways and biological functions that were found significant from gene functional enrichment analysis. For instance, in the ribosome protein bicluster, the enriched TFBS RAP1 is known to be involved in ribosome protein transcription [20]. The ubiquitin cluster is enriched with TFBS of RPN4, which was shown to be part of the ubiquitin fusion degradation pathway [21]. The nitrogen, sulfur and selenium cluster shows significant over-representation of binding sites for TFs CBF1, GCN4, MET31, and MET4. CBF1 is known to induce sulfate assimilation pathway along with MET4 [22], GCN4 is an activator involved with protein biosynthesis, and MET31 is a known transcriptional regulator of sulfur amino acid metabolism [23]. The G1 phase cluster is enriched with MBP1 and SWI4 binding sites. MBP1 and SWI4 are known to act together to regulate late G1-specific transcription of targets and genes for DNA synthesis [24]. The oxidative stress cluster is enriched by CAD1 and YAP1 binding sites, where YAP1 activates the transcription of anti-oxidant genes in response to oxidative stress [25]. The glycolysis regulation cluster is enriched by TFBS GCR1, which is known to be involved in glycolysis [26].

Bayesian biclustering model

Consider a microarray dataset with N genes and P conditions (or samples), in which the expression value of the i^th gene and j^th condition is denoted as y_ij, i = 1, 2, • • •, N, j = 1, 2, • • •, P. We assume that

where K is the total number of clusters (unknown), μ_k is the main effect of cluster k, and α_ik and β_jk are the effects of gene i and condition j, respectively, in cluster k, ε_ijk is the noise term for cluster k, and e_ij models the data points that do not belong to any cluster. Here δ_ik and κ_jk are binary variables: δ_ik = 1 indicates that row (gene) i belongs to cluster k, and δ_ik = 0 otherwise; similarly, κ_jk = 1 indicates that condition (column) j is in cluster k, and κ_jk = 0 otherwise.

When multiple biclusters are allowed, the original plaid model usually finds biclusters greatly overlapping with each other. This effect is quite artificial and is likely due to the nonidentifiability problem caused by the additive assumption made for overlapping clusters. We solve this problem by allowing biclusters to overlap only in one direction, either the gene or the condition direction, but not both. This results in two versions of the BBC model: non-overlapping gene biclustering and non-overlapping condition biclustering. In non-overlapping condition biclustering, a condition can be in one or none of the clusters, but a gene can be assigned to multiple clusters. Mathematically, this constraint can be written as ${\displaystyle {\sum }_{k=1}^K{\kappa }_{jk}}\le 1$. In non-overlapping gene biclustering, a condition can be assigned into multiple clusters, while a gene can only be in no more than one cluster. This corresponds to ${\displaystyle {\sum }_{k=1}^K{\delta }_{ik}\le 1}$. Note that in either of these two versions, different biclusters do not overlap. Without loss of generality, we focus our discussions on the non-overlapping gene biclustering in this paper. Thus the priors of the indicators κ and δ are set so that a condition can be in multiple clusters and a gene be in no more than one cluster, i.e.,

                                                                           k_ij ~ Bernoulli(q_k)

where p_k and q_k are set to be constant. We tested different values for p_k and q_k, and found out that different values do not affect the results much. We used q_k =0.1 and

for the yeast dataset.

We assume a priori that

The hyperpriors for the ${\sigma }_{\mu k}^2,{\sigma }_{\alpha k}^2,{\sigma }_{\beta k}^2,{\sigma }_{\in k}^2,{\sigma }_e^2$ are set to be inverse Gamma distributed.

In our model, an observation y_ij can belong to either one or none of the biclusters. Thus, we can rewrite the probability distribution of y_ij conditional on the cluster indicators. If y_ij belongs to cluster k,(k = 1,2,… ,K), then

If y_ij belongs to none of the clusters, then

With Gaussian zero-mean priors on the effects parameters, we get the marginal distribution of the y_ij conditional on the indicators as:

where Σ is the covariance matrix of Y, and Y = (Y₀, Y₁, Y₂, …, Y_K)^T with Y_k = {y_ij : δ_ikκ_jk = 1}, k ≥ 1, being the vector composed of the data points belonging to cluster k, and Y₀ being the vector of data points belonging to no cluster. More specifically, Σ is a sparse matrix of the form

where Σ_k = cov(Y_k, Y_k) is the covariance matrix of all data points belonging to cluster k whose entries are:

Gibbs sampling for Bayesian biclustering

In order to make inference from the BBC model, we implement a Gibbs sampling method [27] to draw samples from the posterior distribution of the indicator variables, which can be derived by combining equation (8) with a prior distribution on the δ and κ. Initializing from a set of randomly assigned values of δ's and κ's, we sample the column (condition) indicators κ by calculating the following log-probability ratio:

Since each data point belongs to no more than one cluster, we can therefore divide data points into two sets given the current parameters except κ_jk. The first set contains data points not in cluster k, i.e., V₁ = {y_il : δ_ik = 0 or κ_lk =0,l ≠ j}. The second set contains data points that are or can in cluster k, i.e., V₂ = {y_il : δ_ik = 1,κ_lk = 1,l ≠ j}U{y_ij : δ_ik = 1}.

Two data points are independent if they belong to different clusters, therefore we can write the joint likelihood of Y as a product of the joint likelihood for data in V₁ and V₂, respectively. As a consequence, the log-posterior probability ratio can be simplified as

In order to calculate the likelihood term in the above ratio, we need to take the inverse and determinant of the covariance matrices for the vector V₂ in both cases. The dimensions of these covariance matrices are huge in practice (in the order of thousands), so a brute force calculation would be expensive. Since the covariance matrices have the special structures as shown in equation (9), we can simplify the likelihood ratio term. The final simplified form only involves multiplications and additions of matrices with dimension I_k × I_k, where I_k is the number of genes in cluster k given current parameters.

Similarly we can obtain the log- posterior probability ratio for gene indicators δ_ik,

Since our model requires that ${\displaystyle {\sum }_{k=0}^K{\delta }_{ik}\le 1}$, then the gene indicators δ_ik are correlated. We thus need to calculate the log-posterior probability ratio for every δ_ik, k = 1, 2, • • • K, and sample them jointly.

We also sample the effect parameters based on the indicators δ and κ. The gene effects α_ik and condition effects β_jk can serve as scores for genes and conditions in a cluster. Moreover the BBC model is very convenient and coherent in handling missing data in microarray datasets: just treat them as additional unknown variables and iteratively impute them in the Gibbs sampling iterations. Suppose at step t, we have sampled both the indicator variables and effects parameters, we can impute the missing data, say y_ij, by sampling from distribution

In the above procedure, we preset the value K for the total number of clusters. However, the information of K is not available in general. In practice, we search biclusters for a number of K's and select the best K based on the Beyesian information criterion (BIC) [18].

An executable program for the BBC algorithm is available at http://www.people.fas.harvard.edu/~junliu/BBC