Bréhélin et al. BMC Genomics 2010, 11:35 http://www.biomedcentral.com/1471-2164/11/35 achieved when using sequence similarity alone. It also outperforms the accuracy of classical GBA approaches. Results achieved on P. falciparum allowed us to pool and review series of genes involved in particular pro- cesses such as ribosomal complex elaboration and bio- genesis. Notably, it allowed us to highlight potential drug targets. We observed high representation of certain biological functions and an apparent lack of other func- tions. The potential of such an analysis appears to depend on the choice of both the species compared and of the biological conditions screened by the microarray data. Nonetheless, results on P. falciparum show that the approach has high accuracy and potentially allows a function to be inferred for several unknown genes. Methods The model HMMs are probabilistic models that are widely used in biological sequence modeling to recognize protein families and predict gene models [28]. Here they are used for modeling pairs of gene expression profiles. Our HMM has a predefined structure with a mute state M, as depicted in Figure 1, and a set of numerical para- meters Θ. Each prototype is parameterized by a multi- variate normal distribution of dimension equal to the profile sizes—i.e., the number of time points (or experi- ments) of the gene expression profiles. In our experi- ments, we restricted the covariance matrices to be diagonal and equal for all prototypes of one species (homoscedastic assumption). Moreover, there is a prior probability distribution on the prototypes of the query species, and a probability distribution on the outgoing transitions from each of these prototypes and from the M state. Assessing functional conservation We assessed the transfer of function of numerous gene pairs between S. cerevisiae/D. melanogaster and P. falci- parum. We used the BLASTP software of NCBI with default parameters to measure sequence similarity between genes. All gene pairs between the query and the reference species with an e-value below 0.1 were considered. For query genes possessing several homologs in the reference species, the BLAST bit scores (used to compute the e-value) were compared, and gene pairs with scores under 95% of the largest one were removed. This allowed us to associate at least one gene in S. cere- visiae and D. melanogaster with 3 664 and 3 115 genes of P. falciparum, respectively. As some of these P. falci- parum genes are actually associated with more than one gene in each reference species, the total number gene pairs considered was 5 345 between P. falciparum and S. cerevisiae, and 3 904 between P. falciparum and D. melanogaster. For each gene pair, the classical Viterbi algorithm [27] was used to search the prototypes k and Page 11 of 14 k’ that have the highest probability of generating the associated gene expression profiles. Then, each pair for which the ratio P(k’ |k)/P(k’ ) was higher than the threshold ξ = 10 was considered eligible for a functional transfer. Learning the model The learning algorithm takes the number of prototypes K and K’ for the two species as parameters—in the experiments we used K = K’ = 100. It is applied to a learning set L of gene pairs selected by Reciprocal Best Hit (RBH). The procedure involves associating with each gene of the query (respectively reference) species the gene in the reference (resp. query) species that most resembles it, and identifying the gene pairs g - g’ such that g is the best hit for g’, and g’ is the best hit for g. We obtained 1 427 and 1 855 aligned pairs of sequences between P. falciparum and S. cerevisiae, and P. falci- parum and D. melanogaster, respectively. Learning the HMM involves to learn both its structure (i.e. the presence/absence of direct transitions between prototypes), and parameters Θ (i.e. parameters of the Gaussian models, prior probabilities of the prototypes, and transition probabilities). Given a predefined struc- ture, the second problem can be handled with the classi- cal Baum-Welch algorithm [29] used for training the HMMs. Solving the former problem involves deciding, for each prototype pair (k, k’ ), if it involves conserved coexpressed genes. To this end, an initial clustering in K and K’ clusters is computed by running the K-MEANS algorithm [25] on the expression data of each species. Then, for each cluster pair (k, k’ ), the number of gene pairs of L involving genes assigned to clusters k and k’ by the K-MEANS algorithm is computed, and a statisti- cal test is applied to test if this number can be expected by chance (see below). When this is not the case (the hypothesis is rejected), the two clusters are consid- ered to reveal conserved coexpression between the two species, and a transition is added in the HMM from prototype k to prototype k’. At the end of the process, the HMM structure has been built. The Baum-Welch algorithm is then run to train this structure. It is applied to the pairs of expression profiles in L and aims at maximizing the model likelihood: P( | )   P(g, g | ), (g,g) (1) with P(g, g’ |Θ) being the probability of generating the expression profiles associated with genes g and g’ with the HMM with parameters Θ. Thus, the Baum-Welch algorithm searches for the parameters Θ that maximizes Expression (1). It is an iterative algorithm, which starts from an initial set of parameter Θ0, and iteratively rees- timates the parameters at each step of the process. We