Bréhélin et al. BMC Genomics 2010, 11:35 http://www.biomedcentral.com/1471-2164/11/35 Even if better results can be expected from data pro- duced by identical assays, in practice this is often impos- sible. Rather, the aim here is to find gene clusters likely involved in analogous pathways. The approach Our approach involves the detection of clusters of coex- pressed genes in the query species whose coexpression is conserved in the reference species. We implemented this approach with a probabilistic model specially designed for modeling the cluster of coexpressed genes in each species as well as the conservation of coexpres- sion between the two species. The model In each species, coexpression is modeled with prototypes of gene expression profiles. This is an approach similar to that used in classical clustering algorithms such as SOM [24] or K-MEANS [25,26]. The probabilistic model has the following structure (see Figure 1). We use K prototypes of gene expression profiles for P. falci- parum, and K’ prototypes for the reference species. In the experiments below, we use 100 prototypes for each species. Each prototype (of P. falciparum and the refer- ence species) is modeled with a probabilistic, multivari- ate Gaussian distribution defining a probability function over the gene expression profiles of the species. As all prototypes are different, the profile of a given gene has different probabilities depending on the prototype con- cerned, and hence some prototypes are more likely than others for the gene. Genes are associated with their most likely prototype, and all genes associated with the same prototype are considered to be coexpressed. There is a prior-probability distribution on the prototypes of each species, which reflects the proportion of genes that are associated with each prototype in the species. When the coexpression of a set of genes associated with prototype k of P. falciparum is conserved, the homologs of these genes in the reference species also tend to share a common prototype k’ in this species. This is modeled by a transition between k and k’, labeled with a probability P(k’ |k) representing the proportion of genes associated with k whose homologs are associated with k’. Finally, an additional mute state M is used for modeling all non-conserved coexpressions: it is possible to reach the M state from every prototype of P. falci- parum, and reach every prototype of the reference spe- cies from the M state. Prototypes of P. falciparum are labeled with their prior probability P(k), k = 1 ... K. The prior probabilities of the prototypes of the reference P(k’ ), k’ = 1 ... K’ are on the outgoing transitions from the M state. Usually, one prototype is linked with zero or with one prototype of the other species. However, as illustrated by the last P. falciparum prototype in Figure 1, one pro- totype can sometimes be linked with two or more Page 3 of 14 prototypes. This may occur for several reasons. For bio- logical reasons first: since some genes are part of several pathways, the same pathway is not necessarily active in the two species—either due to the species’ specificities, or because the microarray data monitor different condi- tions in the two species. For experimental reasons next: the experimental noise in microarray data can artificially break down a set of coregulated genes into two (or more) different clusters. Conversely, two clusters with a specific but close signature of expression profiles may appear to be similar and associated with the same prototype. This model defines a probabilistic model of pairs of gene expression profiles. From a computational stand- point, it is a Hidden Markov Model (HMM) [27,28], therefore we benefit from all classical algorithms designed for using and training these models (see refer- ence [27] for a description of these algorithms). Given the expression profiles of a pair of query-reference genes, this HMM can be used to compute the probabil- ity density of the profiles under two different hypoth- eses. Under the hypothesis of dependence of the expression profiles—that is, homologous genes belong- ing to a group of conserved coexpression— this is the probability of generating the two profiles by a path in the HMM that uses only direct transitions between pro- totypes. Under the independence assumption of the expression profiles, this is the probability of generating the two profiles by a path that uses the M state: this way, the density only depends on the prototype’s prior probabilities and not on conditional probabilities between prototypes. Learning the model Learning the model involves 1) building the structure (i. e. deciding on the direct transitions between prototypes); 2) training the model to estimate the numerical para- meters of the HMM (i.e. Gaussian distributions, transi- tion probabilities, and prototype prior probabilities). Building the structure involves identifying the conserved coexpressed clusters between the two species. For this purpose, we designed an algorithm which takes as input a set of gene pairs with high sequence similarity selected with the Reciprocal Best Hit procedure (RBH) between the two species (see Methods). An initial gene clustering of each species in K and K’ clusters is computed with the K-MEANS algorithm [25]. Next, a statistical proce- dure computes, for each cluster pair (k, k’), the number of gene pairs with high sequence similarity between k and k’, and tests if this can be expected by chance. When this is not the case, a transition between k - k’ is added to the model. By the end of the process, the structure of the HMM has been built, and the classical Baum-Welch algorithm [29] is then run to train the model.