A model selection criterion for model-based clustering of annotated gene expression data

1 Introduction

Genome annotation broadly refers to the set of meta-data associated with the coding regions in the genome, typically including the identification of the location of each gene as well as a determination of the functions related to the gene product (e.g. protein or RNA). In particular, gene annotations correspond to known functions related to the gene product, including molecular functions, biological pathways, or the cellular location of the gene products. A variety of well-known unified databases have been constructed with known functional annotations collected from bibliographic sources across species, including the Gene Ontology (GO) (Ashburner et al., 2000), the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) or the MSigDB (Molecular Signatures) databases (Liberzon et al., 2011). Although such databases contain a rich source of functional information about the genome in a large variety of species (e.g. Arabidopsis thaliana, human, rat, mouse, fly), our knowledge of functional annotations is often far from complete (Tipney and Hunter, 2010).

In recent years, substantial improvements in high-throughput technologies, such as microarrays (Schena et al., 1995) and more recently RNA sequencing (RNA-seq) (Mortazavi et al., 2008), have enabled the simultaneous measurement of the expression levels of tens of thousands of genes. A rich body of work is now available on the appropriate statistical analyses of such gene expression data, including the analysis of differential expression (Smyth, 2004; Anders and Huber, 2010) and co-expression analyses to identify groups of genes with similar profiles across several experimental conditions or over the course of time (Yeung et al., 2001; Rau et al., 2015). The latter is of particular interest in this work, as identifying genes that share the same dynamic patterns of expression may help identify groups of genes that are involved in similar biological processes and generate hypotheses about the functional properties of poorly characterized genes (Eisen et al., 1998; Jiang et al., 2004). Reviews and comparisons of different clustering methods for gene expression data may be found in Datta (2003).

In practice, annotation databases are often used to perform a posteriori validation and interpretation of co-expressed gene clusters through tests of functional enrichment (Steuer et al., 2006). Such functional annotation may instead be directly integrated into the clustering model itself. For example, Tari et al. (2009) incorporate GO annotations as prior knowledge in a fuzzy c-means clustering. Verbanck et al. (2013) proposed a clustering approach based on a distance defined conjointly on the similarity among expression profiles and that among functional profiles. Pan (2006) and Huang et al. (2006) proposed including gene annotation as prior information in a stratified mixture model. However, the inclusion of gene annotation directly in the model itself in this way may be questionable, particularly when they are also used to validate the gene clusters a posteriori. Moreover, as gene annotations tend to be incomplete, biases may be introduced if they are directly incorporated in the model, as unannotated genes (which represent those known to be unassociated with a given function as well as those of unknown function) may be erroneously separated from annotated genes.

One alternative to such approaches is to define a clustering model that accounts for external gene annotations without directly including them in the model itself. To this end model-based clustering provides a convenient framework, as it 1) allows for a large set of clustering models to be fit to the gene expression alone, and 2) facilitates the choice among this set, a parsimonious model that simultaneously provides a good fit to the data and coherence with the external gene annotations. In this work, we address these points by proposing a model selection criterion that accounts for external gene annotations.

The rest of this paper is organized as follows. In Section 2, we present the context of model-based clustering and review classic model selection criteria. Our proposed annotated model selection criterion is presented in Section 3, and numerical illustrations of its behavior are presented on simulated data in Section 4 using Gaussian mixture models. Finally, we illustrate a co-expression analysis of real RNA-seq data in Section 5, and a discussion ends the paper.

2 Model-based clustering and model selection

Let y be the (n×q) matrix of observed gene expression, where n is the number of genes and q the number of biological samples. The vector y_i denotes the expression of gene i (i=1, …, n) across the q samples. In the context of model-based clustering, the data y are assumed to be sampled from a finite mixture density of K random variables, each with parameterized density ϕ(y_i; a_k), k=1, …, K, where the mixture parameters (a₁, …, a_K) are all assumed to be distinct. The density of y may thus be written as

where θ_K=(p₁, …, p_K-1, a₁, …, a_K) are the parameters of the mixture model, and (p₁, …, p_K) are the mixing proportions with p_k∈(0, 1) for all k, ∑k=1Kpk=1.

For parameter estimation, the mixture model in Equation (1) may be thought of as an incomplete data structure model where z is the (n×K) matrix of unknown mixture labels, where z_ik=1 if gene i is from group k and 0 otherwise. Note that this matrix defines a partition of the genes.

Using the mixture labels z, the completed density of y may be written as follows:

The maximum likelihood estimate θ^K of the mixture parameters is computed through the expectation-maximization algorithm (Dempster et al., 1977) by replacing the unknown labels z in Equation (2) with z^, defined as:

z^ik={1if arg maxℓτiℓ(θ^K)=k0otherwise,

where τiℓ(θ^K) denotes the conditional probability given y_i of the ℓ th mixture component under θ^K:

τiℓ(θ^K)=p^ℓϕ(yi; a^ℓ)∑t=1Kp^tϕ(yi; a^t).

In the context of model-based clustering, one important task is the choice of an appropriate model, most notably the revelant number of clusters K. To this end, a standard model selection criterion is the Bayesian information criterion (BIC) (Schwarz, 1978):

BIC(K)=logf(y; K^, θ^K)−νK2log(n),

where θ^K is the maximum likelihood estimator of the mixture parameters and v_K the number of free parameters in the model with K components. This criterion is an asymptotic approximation of the logarithm of the integrated likelihood:

f(y; K)=∫θKf(y; K, θK)π(θK)dθK,

where π(θ_K) is a weakly informative prior distribution on θ_K.

An alternative to the BIC is the integrated completed likelihood (ICL) criterion (Biernacki et al., 2000):

where Ent(K) is the estimated mean clustering entropy

Note that the ICL is a BIC-like approximation of the logarithm of the completed integrated likelihood:

f(y, z; K)=∫θKf(y, z; K, θK)π(θK)dθK.

Because of the additional entropy term defined in Equation (4), the ICL favors models that lead to data partitions with the greatest evidence in terms of classification.

More recently, Baudry et al. (2014) proposed an ICL-like criterion that takes advantage of the potential explicative ability of external categorical variables u=(u¹, …, u^R) where uiℓr=1 indicates that the gene i is in category ℓ for the r^th external categorical variable and 0 otherwise. The idea is to choose a classification z based on y that is coherent with u. Assuming that y and u are conditionally independent given z, the Supervised Integrated Completed Likelihood (SICL) criterion is an asymptotic approximation of the logarithm of the integrated completed likelihood:

f(y, u, z; K)=∫f(y, u, z; K, θK)π(θK)dθK.

The SICL criterion is defined as follows:

where U_r is the number of levels of the variable u^r,

nkℓr=card{i:zik=1 and uiℓr=1},

and nk.=∑ℓ=1Urnkℓr. The last additional term in Equation (5) quantifies the strength of the link between the categorical variables u and the classification z.

3 Taking genome annotations into account

As previously stated, the objective of this work is to make use of external gene annotations to choose a model for which clusters may be meaningfully interpreted both with respect to their expression profiles and their functional properties. To do so, we propose a novel model selection criterion that highlights the association between the clusters of expression profiles and the functional annotations associated with a subset of genes. Since gene annotations are binary variables (i.e. a gene is either annotated or unannotated), it may seem natural to directly use the SICL defined in Equation (5). However, in contrast to the situation considered by Baudry et al. (2014), gene annotation information is often incomplete. More precisely, for each of the G annotation terms, indexed by g, the available information u^g is as follows:

uig={1if gene i is known to be implicated in function g,0if gene i is not known to be implicated in function g.

Note that uig=0 can indicate that information is missing (i.e. gene i has not yet been identified for annotation g) or that gene i is known to be unrelated to annotation g. As such, uig=0 does not represent the null level of variable and thus represents an incomplete binary variable. For this reason, the SICL criterion is not an appropriate measure of the link between an external annotation u^g and a classification z, and a specific criterion must be defined to incorporate the gene annotation information into the model selection step. To this end, we propose the integrated completed annotated likelihood (ICAL) criterion as follows.

For each gene annotation u^g, we first define the random matrix b^g of latent variables indicating the allocation of the annotations among the K clusters:

Each row of the matrix b^g is a random vector following a multinomial distribution with parameters uig and (p1g, …, pKg) if uig>0, and is the null vector 0 if uig=0.

For the sake of simplicity, we first derive ICAL when a single external annotation b¹ is available. ICAL aims to select the clustering model that maximizes the logarithm of the integrated annotated likelihood:

As for the definition of the SICL, the variables y and b¹ are assumed to be conditionally independent given z. Using Bayes formula, we have

f(y, z, b1; K, θK)=f(y, z; K, θK)f(b1|y, z; K, θK).

Note that since y and b¹ are assumed to be independent given z, the conditional distribution of b¹ given z does not depend on y or the mixture parameters. Thus, as f(b¹|y, z; K, θ_K)=f(b¹|z; K), it follows that:

The last term in Equation (8) can be approximated with ICL(K) from Equation (3), and the first term may be approximated with

logf(b1|z^; K)=∑k=1Knk1lognk1n1,

where n1=card{i:ui1=1} and nk1=card{i:z^ik=1 and ui1=1}. Finally, an asymptotic approximation of the expression in (7) leads to the ICAL criterion:

ICAL(K)=ICL(K)+∑k=1Knk1lognk1n1.

The generalization of this criterion to the case where G>1 gene annotations are available is straightforward. The aim is now to maximize the logarithm of the integrated annotated likelihood:

logf(y, z, b1, …, bG; K)=log∫θKf(y, z, b1, …, bG; K, θK)π(θK)dθK.

Assuming that b¹, …, b^G and y are conditionally independent given z, we have

logf(y, z, b1, …, bG; K)=logf(b1, …, bG; z, K)+log∫θKf(y, z; K, θK)π(θK)dθK.

Assuming in addition that b¹, …, b^G are independent and that gene annotations are missing at random, we can write

leading to the generalized ICAL criterion:

Comparing ICAL and SICL If we ignore the uncertainty associated with uig=0 (i.e. that gene i could either be unassociated with function g or that this information is missing), the SICL criterion could be considered to choose the model dimension K. In this case, using the notation from Section 2 and defining n_k the size of the cluster k, the SICL may be written as follows:

SICL(K)=ICL(K)+penSICL,

where

penSICL=∑g=1G∑k=1Knk1glognk1gnkg+∑g=1G∑k=1Knk0glognk0gnkg,=∑g=1G∑k=1Knk1glognk1g+∑g=1G∑k=1Knk0glognk0g−G∑k=1Knklognk.

On the other hand, using the notation from Section 2 and defining n.1g=∑k=1Knk1g, the ICAL may be written as follows:

ICAL(K)=ICL(K)+penICAL,

where

penICAL=∑g=1G∑k=1Knk1glognk1gn.1g,=∑g=1G∑k=1Knk1glognk1g−∑g=1Gn.1glogn.1g.

We note that the last term in the equation above is a constant independent of K. Finally, we can rewrite ICAL as a function of SICL:

From Equation (11), we note that the SICL takes into account both modalities (0 and 1) of the external variables u, while the ICAL discards the null modality (the −∑g=1G∑k=1Knk0glognk0g term). Moreover, it can be seen that the ICAL penalizes a large number of clusters, while the SICL does not (the G∑k=1Knklognk term). As such, the ICAL tends to select parsimonious models with a relatively small number of clusters, as compared to SICL.

It is also helpful to consider the behavior of the ICAL and SICL criteria in extreme conditions. If the number of clusters K equals 1, the ICAL penalty pen_ICAL equals zero whereas SICL penalty pen_SICL is not null (∑g=1Gn1glogn1gn+∑g=1Gn0glogn0gn). In contrast, if the number of clusters K is equal to the number of observations, with one gene per cluster, the SICL penalty pen_SICL equals zero whereas the ICAL penalty is not null (∑g=1Gn1glogn1g). In general, ICAL tends to merge clusters to group genes annotated for the same function, reducing the number of optimal clusters K with respect to the optimal number of clusters selected by ICL. SICL tends to split clusters in order to obtain clusters made up only of annotated genes, increasing the number of optimal clusters with respect to the optimal number of clusters selected by ICL. In other words, SICL tends to select more complex models than ICL while ICAL tends to favor more parsimonious models than ICL. Note that this behavior of ICAL and SICL is a general trend, not a rule: ICAL does not always merge clusters and SICL does not always split them since clusters for different solutions are not necessarily nested in each other.

Code to implement our method is available in the R package ICAL, which may be found at the following website: https://github.com/Gallopin/ICAL.

4 Numerical illustrations

4.1 Simulation settings

To illustrate the behavior of the proposed ICAL criterion, we consider a numerical example. We simulate 200 observations from a mixture of four bivariate Gaussian distributions, 100 independent times (see parameters in Table 1). The first two components are close to one another while the third and fourth are clearly distinct from the first two and also distinct from each other. For a given model indexed by K, the estimation of parameters is performed with the R package Rmixmod (Biernacki et al., 2006; Lebret et al., 2015) for a Gaussian mixture model with diagonal variance matrix (that is, the p_kL_kB_k model in the notation of the Rmixmod package, corresponding to clusters with variable proportions, variable volumes, variable shapes and vertical or horizontal orientation). We estimate the parameters for models with the number of clusters K varying from 1 to 10 and perform model selection to select the most appropriate number of clusters. Over the 100 replicated datasets, the BIC most frequently selects four clusters (81 times). Indeed, we note that these clusters correspond to the simulated Gaussian components. The ICL criterion selects either three (54 times) or four clusters (46 times), as it tends to merge the two similar components (1 and 2 from Table 1).

Table 1

Parameters of simulated datasets: the first two components are close to one another while the third and fourth are clearly distinct from the first two and also distinct from each other.

Component	Mixing proportions	Component distribution
1	0.25	N((00) , (1001.7))
2	0.25	N((05) , (1001.7))
3	0.25	N((98) , (0.1001))
4	0.25	N((90) , (1000.8))

We illustrate the potential utility of accounting for external gene annotations in model selection by simulating such annotations and performing model selection with the corresponding SICL and the ICAL criteria. We simulate three types of functional annotations: u_A, u_B and u_C (see Figure 1). The genes annotated for the first function u_A are shared by the two closest mixture components (components 1 and 2 from Table 1). This annotation is designed to be associated to the components in the sense that it suggests the interest of merging the two clusters, as they share similar joint distributions and external annotations. The genes annotated for the second function u_B are shared only by the two clearly distinct components (components 3 and 4 from Table 1). This annotation is designed to be unassociated with the components: although the components share a similar function, their joint distributions are too distinct to be merged from a modelling point of view. Finally, the genes annotated for the third function u_C are randomly spread over the four components: meaning the annotation is mixed (half associated/half unassociated). For each function, we simulate the annotation using binomial random variables, with parameters fixed to yield on average n_annot annotated genes over 200 possible genes. In the following, we tested two levels of annotation densities dannot=nannotn: dannot=0.05 and d_annot=0.25. Since the ICAL and SICL criteria can be used with more than one external annotation, we also illustrate the potential utility of including more than one annotation. For each set of annotations used in our simulation, the number of annotations varies from 1 to 12. Each annotation in the set is simulated independently.

Figure 1:

Illustration of a simulated dataset and three annotation patterns. For each figure, the 200 observations are drawn from a mixture of Gaussian bivariate components whose parameters are defined in Table 1: circles, triangles, inverted triangles and diamonds correspond to components 1–4. The three figures correspond to three annotation patterns: associated annotation u_A (A), unassociated annotation u_B (B) and mixed annotations u_C (C). For each annotation, the 20 annotated genes are represented by colored bold crosses.

4.2 Simulation results

All penalized criteria (BIC, ICL, SICL and ICAL) versus the number of clusters for one simulated dataset with d_annot=0.05 are displayed in Figure 2A. Over the 100 simulated datasets, ICAL selects three clusters 72 times, merging the two closest components 1 and 2 (Table 2). This three cluster solution is meaningful with respect to the information provided by u_A, as all annotated genes are attributed to the same cluster. In this case, the external information provided by the associated annotation u_A reinforces the model selection. Using the same pattern as u_A (annotations shared by components 1 and 2 only), we simulate 12 independent associated annotations with density d_annot=0.05. The evidence to merge clusters increases with the number of annotations (one versus 12 annotations). The number of annotations in the set (12 annotations) is chosen so that the evidence to merge clusters was sufficiently strong: over the 100 simulated datasets, ICAL systematically selects a three cluster solution (Table 2). For this set of 12 external annotations, the peak of the ICAL displayed in Figure 2B is much sharper than the peak of ICL. In contrast, SICL more frequently selects a four- or even five-cluster solution, as it leads to a preference of smaller clusters containing only annotated genes (i.e. a high specificity of annotation within each cluster). This demonstrates the utility of the ICAL criterion over the SICL, as it does not correctly take into account the specificity of gene annotation.

Figure 2:

BIC, ICL, SICL and ICAL information criteria versus the number of clusters on one simulated dataset for the informative annotations: u_A (A) and uA1, …, uA12 (B). Triangles indicate the maximum value attained by each criterion.

Table 2

Number of simulated datasets for which each model (K=1, …, 10) was selected by BIC, ICL, SICL and ICAL for several external annotations with density d_annot=0.05 over 100 independent datasets simulated with parameters detailed in Table 1.

		K	1	2	3	4	5	6	7	8	9	10
		BIC			19	81
		ICL			53	47
Associated annotations	uA	SICL			51	49
		ICAL			72	28
	uA1, …, uA12	SICL			27	70	3
		ICAL			100
Unassociated annotations	uB	SICL			53	46	1
		ICAL			53	47
	uB1, …, uB12	SICL			53	44	3
		ICAL			53	47
Mixed annotations	uC	SICL			50	50
		ICAL			67	33
Multiple annotations	uA, uB, uC	SICL			48	51	1
		ICAL			81	19

The model most commonly selected for each criterion is highlighted in bold.

For unassociated annotation u_B with d_annot=0.05, the behavior of the information criteria versus the number of clusters for one simulated dataset is displayed in Figure 3A. We note that the ICAL criterion behaves similarly to the ICL. Over the 100 simulated datasets, ICAL, as ICL, leads to some uncertainty as to whether a three cluster solution (53 times) or a four cluster solution (47 times) is best (Table 2). In this case, the annotation u_B is not related to the components and has no impact on the resulting clustering, even if the number of annotations in the set is increased to 12, each simulated with the same pattern as u_B as displayed in Figure 3B.

Figure 3:

BIC, ICL, SICL and ICAL information criteria versus the number of clusters on one simulated dataset for the non informative annotations: u_B (A) and uB1, …, uB12 (B). Triangles indicate the maximum value attained by each criterion.

Finally, for the mixed annotation u_C with d_annot=0.05, ICAL most frequently selects three clusters (67 times) or four clusters (33 times). Because the annotation u_C is mixed, there is less evidence to merge the two clusters on the left than in the case of the informative annotation u_A. Using the three types of annotations at the same time (u_A, u_B, u_C), the ICAL criterion almost systematically selects three clusters (Table 2).

The potential utility of accounting for external gene annotations in model selection is highlighted in the numerical results summarized in Table 2. First, these results illustrate that the SICL is not well-adapted to account for gene annotations in model selection; at best SICL behaves like ICL and at worst, erroneously splits clusters that should be merged. However, if the external information is associated to the components, even partially so, the use of the ICAL criterion improves model selection in terms of functional interpretability. If the external information is unassociated to the components, the ICAL criterion simply behaves like the ICL.

To assess the influence of annotation density on our results, we repeated the experiment for a higher density (d_annot=0.25), as summarized in Table 3. Generally speaking, we note that when the number of annotations in the set increases, the evidence to merge clusters is stronger; for a higher density of annotations, a smaller number of annotations in the set is needed to merge clusters. In particular, for a lower density of annotations (d_annot=0.05), increasing the number of associated low density annotations in the set (from one to twelve) also increases the level of confidence to merge clusters; clusters 1 and 2 are merged 72% of the time for a single annotation, and systematically merged in all simulated datasets for a set of twelve annotations (Table 2). On the other hand, for a higher density of annotations (d_annot=0.25), clusters 1 and 2 are systematically merged even for a set of one single annotation (Table 3). If the annotation density is high (i.e. where more than half of the genes in clusters 3 and 4 are annotated), multiple unassociated annotations can lead to overly parsimonious solutions, such as the two cluster solutions in Table 3.

Table 3

Number of simulated datasets for which each model (K=1, …, 10) was selected by BIC, ICL, SICL and ICAL for several external annotations with density d_annot=0.25 over 100 independent datasets simulated with parameters detailed in Table 1.

		K	1	2	3	4	5	6	7	8	9	10
		BIC			19	81
		ICL			53	47
Associated annotations	uA	SICL			52	48
		ICAL			100
	uA1, uA2, uA3	SICL			53	47
		ICAL			100
Unassociated annotations	uB	SICL			53	47
		ICAL			53	47
	uB1, uB2, uB3	SICL			53	47
		ICAL		6	47	47
Mixed annotations	uC	SICL			52	48
		ICAL			94	6
Multiple annotation	uA, uB, uC	SICL			51	49
		ICAL			100

The model most commonly selected for each criterion is highlighted in bold.

5 RNA-seq data analysis

5.3 Model selection

We compare the results of model selection performed by the four different criteria presented in Sections 2 and 3: BIC selects 28 clusters, ICL and SICL select 23 clusters while ICAL selects 20 (see Figure 4). Figures 5 and 6 are heatmaps of the resulting clusters for the ICL and ICAL solutions respectively. The approximate correspondences between clusters in the ICAL and ICL solutions are displayed in Table 5. Although the result of the former is not perfectly nested in the latter, in many cases the attribution of genes to clusters in the ICAL solution is a result of collapsing or partially collapsing several clusters from the ICL solution. For example, the ICAL merges most of cluster 2 and parts of clusters 5 and 18 from the ICL solution because they share similar expression profiles and functional annotations, as illustrated in Figure 7. This suggests that the ICL favors a slightly more complex solution, as expected; we next investigate whether the more parsimonious solution of the ICAL appears to be coherent given the set of CP used.

Figure 4:

BIC, ICL, SICL and ICAL information criteria (respectively A–D) versus the number of clusters for the pig RNA-seq data for 10 independent initializations, represented by 10 gray circles for each number of clusters K. Solid lines link the maximum of the criteria over the 10 initializations over the collection of models. The red crosses correspond to the maximum of the criteria, which correspond to the selected models.

Figure 5:

Heatmap of the ICL clusters. The RNA-seq data, expressed in log-cpm is centered and scaled. Colors on the heatmap reflect the level of normalized expression: low in red, high in yellow. The 23 clusters are indicated by colors on the left side of the heat map.

Figure 6:

Heatmap of the ICAL clusters. The RNA-seq data, expressed in log-cpm is centered and scaled. Colors on the heatmap reflect the level of normalized expression: low in red, high in yellow. The 20 clusters are indicated by colors on the left side of the heat map.

Table 5

Approximate composition of the 20 clusters of the ICAL solution with respect to the 23 clusters of the ICL solution.

ICAL clusters	ICL clusters
Cluster 1	129
Cluster 2	15
Cluster 3	10
Cluster 4	11
Cluster 5	125	+	12	+	1220	+	1222
Cluster 6	2	+	125	+	1318
Cluster 7	8
Cluster 8	4	+	9	+	16
Cluster 9	3	+	126
Cluster 10	7	+	126
Cluster 11	13	+	1420	+	1222
Cluster 12	23
Cluster 13	7	+	17	+	1318
Cluster 14	21
Cluster 15	19
Cluster 16	1
Cluster 17	14	+	1318
Cluster 18	1
Cluster 19	1420	+	125
Cluster 20	16

Lines in bold correspond to clusters of the ICAL solution that are formed by several clusters or parts from clusters of the ICL solution. For example, cluster 5 of the ICAL solution is approximately made of clusters 12 and parts of clusters 5, 20 and 22 of the ICL solution.

Figure 7:

Heatmap of cluster 6 of the ICAL solution. The RNA-seq data, expressed in log-cpm is centered and scaled. Colors on the heatmap reflect the level of normalized expression: low in red, high in yellow. The colors and numbers on the left side of the heatmap indicate the corresponding clusters from the ICL solution that have been merged in the ICAL solution.

For the ICL and ICAL solutions, we examine associations between clusters and CP using Fisher’s exact test. Significant p-values are summarized in Table 6. The ICAL criterion yields a clustering that maximizes the number of genes annotated in each cluster for each CP while still only grouping genes that share sufficiently similar expression profiles. For example, we note that CP8 is associated with two different clusters in the ICL solution, while it is associated with a single cluster in the ICAL solution; similarly, CP10 is associated with three clusters in the ICL solution and only two clusters in the ICAL solution. On the other hand, although clusters 10 and 17 in the ICAL solution both share annotations for CP10, these clusters are not collapsed into one using the proposed criterion, as their expression dynamics are too different. As such, the ICAL solution appears to enable the identification of more biologically interpretable clusters than the ICL, while still ensuring that the clustered genes share sufficiently similar expression dynamics.

Table 6

Table of associations between clusters and CP for the ICL solution (A) and the ICAL solution (B).

	Size	CP1	CP2	CP3	CP4	CP5	CP6	CP7	CP8	CP9	CP10
(A) [ICL solution]
Cluster 2	58		★	★	★
Cluster 5	203								★
Cluster 6	47										★★
Cluster 7	258	★					★				★
Cluster 8	96					★★
Cluster 10	287									★
Cluster 14	225										★★
Cluster 22	144			★★					★★★
(B) [ICAL solution]
Cluster 3	297									★
Cluster 5	379			★★					★★★
Cluster 6	156		★★	★
Cluster 7	92					★
Cluster 10	267	★					★★				★★
Cluster 17	235										★★

Associations are detected using Fisher’s exact tests: the number of stars indicates the value of the p-value (★ below 0.01, ★★ below 0.001, ★★★ below 0.0001).

Finally, we note that the ICAL solution exhibits two clusters of particular interest with respect to the biological processes studied: Cluster 5 (379 genes) is associated with CP3 (reactome homeostasis, p=0.0002) and CP8 (KEGG ECM receptor interaction, p=0.00001). Cluster 10 (297 genes) is associated with CP1 (reactome metabolism of lipids and lipoproteins, p=0.002), CP6 (reactome fatty acid triacylglycerol and ketone body metabolism, p=0.005) and CP10 (KEGG peroxisome, p=0.0001), all of which correspond to fatty acid metabolism. Both clusters 5 and 10 contain unknown genes that may be good candidates for follow-up studies to determine whether they may be implicated in the corresponding canonical pathways.

6 Discussion

In this paper, we present a novel way to incorporate functional annotations into model-based clustering of gene expression data. To this end, we have developed a model selection criterion, the Integrated Completed Annotated Likelihood (ICAL) which is designed to select the model that jointly maximizes the goodness-of-fit to the data and the association of clusters and annotations. From a biological point of view, the ICAL criterion aims to select models with more interpretable clusters than those selected by BIC or ICL. It is important to note that the functional annotations are not directly included in the clustering model and are only used to select the best model. This approach is a good compromise between two opposite strategies: including functional annotations directly in the clustering model (Morlini, 2011) or excluding them altogether and using them only to validate clusters a posteriori. Since we do not include annotations in the clustering model, we detect associations between annotations and clusters with a stronger evidence than if we had included the external annotations in the clustering model. In particular, the ICAL criterion is a good way to include prior biological expertise without according it too much importance, which is a good balance between what can be observed in the data and what experts expect to see in the data.

As illustrated in numerical simulations, the model selected by ICAL depends on the quality of the annotation information provided. Selecting the appropriate annotations to include in ICAL is an important step and should be done based on expert knowledge. We also suggest the use of gene annotation databases that are curated manually by experts, such as the gene sets collection from the MSigDB database (Liberzon et al., 2011). However, the choice of annotations should reflect the biological functions of interest to a particular study, making it difficult to provide general guidelines about how such annotations should be chosen in practice. We note that if the annotations chosen are not relevant to the cluster patterns present in the data, they do not contribute any information and ICAL tends to behave like the ICL criterion.

In this work, we applied the ICAL using the framework of Gaussian mixture models, but the extension to other mixture models is straightforward; including Poisson (Rau et al., 2015) or Dirichlet multinomial mixture models (Holmes et al., 2012). In addition, this model selection strategy may be useful for other types of data which may also be associated with incomplete external annotations (e.g. sociology, marketing).