 Methodology article
 Open Access
 Published:
Gene expression trees in lymphoid development
BMC Immunology volume 8, Article number: 25 (2007)
Abstract
Background
The regulatory processes that govern cell proliferation and differentiation are central to developmental biology. Particularly well studied in this respect is the lymphoid system due to its importance for basic biology and for clinical applications. Gene expression measured in lymphoid cells in several distinguishable developmental stages helps in the elucidation of underlying molecular processes, which change gradually over time and lock cells in either the B cell, T cell or Natural Killer cell lineages. Largescale analysis of these gene expression trees requires computational support for tasks ranging from visualization, querying, and finding clusters of similar genes, to answering detailed questions about the functional roles of individual genes.
Results
We present the first statistical framework designed to analyze gene expression data as it is collected in the course of lymphoid development through clusters of coexpressed genes and additional heterogeneous data. We introduce dependence trees for continuous variates, which model the inherent dependencies during the differentiation process naturally as gene expression trees. Several trees are combined in a mixture model to allow inference of potentially overlapping clusters of coexpressed genes. Additionally, we predict microRNA targets.
Conclusion
Computational results for several data sets from the lymphoid system demonstrate the relevance of our framework. We recover wellknown biological facts and identify promising novel regulatory elements of genes and their functional assignments. The implementation of our method (licensed under the GPL) is available at http://algorithmics.molgen.mpg.de/Supplements/ExpLym/.
Background
The study of gene regulatory mechanisms controlling cell proliferation and differentiation is central in developmental biology. Because all hematopoietic cells are easily obtained as individual cells, and due to high clinical interest, the development of lymphocytes is particularly wellstudied [1, 2]. In mammals, all blood cells develop from pluripotent, selfrenewing hematopoietic stem cells (pHSC) of the bone marrow. In the classical model, these pHSC differentiate into common myeloerythroid progenitors and common lymphoid progenitors [3]. The latter give rise to all cells of the adaptive immune system, including T, B and natural killer cells, which are the focus of our work.
Lymphocytes are well characterized; they can be purified by fluorescence activated cell sorting (FACS) exploiting the large variety of cell surface antigens, which appear in specific order during differentiation as the result of a linear sequence of genomic rearrangements at the T and B cell receptor loci [4, 5]. Based on this, lineagespecific expression and roles of transcription factors have been studied extensively [1, 2, 6]. It has been shown, for example, that Gata3 is required for CD4 T cell maturation and that Runx3 silences the CD4 gene in CD8 T cells. Very recently, a new class of regulatory RNAs, microRNAs, have been identified as being involved in lymphocyte cell development [7–9].
Several groups [4, 5, 10–12] have combined FACS mediated cell sorting and mRNA expression profiling to derive a more comprehensive picture of the lymphocytes in distinguishable developmental stages. Our interest focuses on these patterns of gene expression in the distinct stages of the developmental tree, the developmental profiles of genes; see Fig. 1 for a developmental tree. Observing such patterns, the first natural question to ask is whether further genes exhibit the same developmental profile; for example, are there other genes coexpressed with Gata3. It is reasonable to assume that genes with a prescribed pattern of expression, such as "upregulated in proliferating cells", might be relevant for specific functions of cells in a particular stage of differentiation. Clearly, not all relevant developmental profiles are known beforehand, so clustering is the next logical step. Clustering allows us to divide genes into groups of similar developmental profiles, some of which will be irrelevant–genes expressed in all stages–others will differ in distinct branches of the developmental tree and thus indicate relevance for differentiation. Once the gamut of developmental profiles is determined, further questions can be addressed with statistical methods: which regulatory effects might cause differentiation, which subgroups of developmental stages share regulatory patterns or at which developmental stage is the difference in expression between two groups the largest. Prior work in this context relies on classical clustering methods, such as selforganizing maps [4, 5], hierarchical clustering [12], or on performing tests of differential expression between cell types of interest [11]. Further studies concentrated on smallscale data, where selected genes are used to infer regulatory networks. One such study applied a statespace model to infer networks of T cell activation [13]. Troncale and colleagues adopted Petri Nets to model and infer regulatory networks of early pHSC development [14], while Basso and colleagues proposed a novel algorithm for a similar task [15].
Classical clustering relies on distance functions between developmental profiles such as correlation or Euclidean distance, which neglect the dependence structure of the developmental tree (Fig. 1). As a matter of fact, the clustering result does not change if one permutes all the variables. Biology suggests however, that the very sequence of changes does matter as this exact sequence of events is what takes a cell from pluripotent to, say, mature Bcell. Thus we propose dependence tree models–see [16] for the discrete variate version–to model expression during the course of development. Our model assumes that the dependence of gene expression between subsequent stages is the most relevant one for identification of coexpressed genes. We assume that gene expression has been measured for a sufficient number of stages, in particular those relevant for differentiation processes, and that the cell population in a particular stage is sufficiently pure. The disagreement between reality and our assumptions is subsumed as noise, which our method can successfully deal with on simulated data. If we consider all pairwise dependencies between developmental stages our model would be equivalent to a multivariate Gaussian distribution with full covariance matrix. Due to the complexity the estimation of such models is prone to overfitting [17, 18]. The dependence tree model represents a tradeoff between methods assuming independence between variables, such as kmeans and hierarchical clustering, and complex models, such as multivariate Gaussians, which makes estimation more robust.
With one such tree we can find genes with a specified developmental profile, for example similar to the developmental profile of Gata 3, by ranking genes in order of decreasing likelihood under the tree. To cluster developmental profiles we combine several trees with the same topology but with distinct parameters in a classical mixture model [17]; tree topologies are taken from the biological literature. Thus we obtain a robust and flexible statistical model for clustering genomewide mRNA expression data sets, which takes the inherent dependencies between developmental stages explicitly into account. The resulting clusters of genes sharing similar developmental expression profiles are wellsuited for a subsequent search for common regulators such as transcription factors or microRNAs.
Our choice of model class is motivated by the successful application of mixtures of complex statistical models to the analysis of mRNA expression timecourses. There, models that take temporal dependencies into account, such as Splines [19, 20], Autoregressive models [21] or Hidden Markov models [22], outperform simpler models, which assume independence of the variables, for example kmeans, selforganizing maps or hierarchical clustering.
For discrete variates, dependence trees were first proposed by Chow and Liu [16], who showed that efficient computation is possible. Mixtures of trees were first proposed and applied in image recognition problems [23], where more efficient versions of the structure learning algorithm for sparse data sets became necessary. In bioinformatics, mixtures of trees were applied to infer mutation events in HIV strains [24]. We present an extension of the dependence trees to continuous variates, requiring modifications to the densities and provide a framework for robust clustering based on mixtures. To the best of our knowledge, there is no prior work on genomescale mRNA expression analysis in which the developmental tree structure is taken into account. Both the biological application and our approach of combining tree models with mixture estimation for this purpose is novel. However, the main methodological ingredients are wellestablished. Our advanced statistical framework allows us to identify clusters of genes with similar developmental profiles. We detect interesting groups of genes not found using standard techniques, such as selforganizing maps [25], in developing lymphoid cells. Results on simulated data show the conditions under which our method has a technical advantage. From our clustering results we can identify plausible regulatory roles of microRNAs known to be involved in hematopoiesis. We provide a graphical user interface and a web database of clustering results; see [26] for implementations, a tutorial on how to use the tools, and a web database with the results presented below. Our findings suggest that our framework is wellsuited for analysis of genomewide expression data from detailed cell development studies.
Results/Discussion
In the next two sections, we describe the dependence trees and how they are combined in a mixture to find groups of developmental profiles. Subsequently, we present the results of the application of our method to three lymphoid cell datasets. In the last subsection, we analyze the groups of genes, given by our mixture of dependence trees (MixDTrees) results, for common microRNA binding sites patterns, in order to gain insights into regulatory function of microRNAs.
Dependence trees
The main assumption behind the dependence trees (DTree) is that expression levels of a particular developmental stage depend primarily on expression levels of the immediately preceding stage. For example, cf. Fig. 2, we can approximate the joint probability density function (pdf) of four random variables (X_{ A }, X_{ B }, X_{ C }, X_{ D }) by
In other words, we condition the probability of a given variable on its immediate predecessor, in accordance with the tree structure shown in Fig. 2. There, also a cluster of hypothetical genes with similar developmental profiles is depicted (Fig. 2, right). The genes display average expression in stage A, upregulation in stage B, downregulation in stage C and upregulation in stage D. Furthermore, the genes have clearly distinct expression intensities, but similar relative expression changes. Genes strongly overexpressed in B are also strongly underexpressed in C and strongly expressed in D. These dependencies are reflected in the correlation between these stages. For example, A and B (or B and D) are positively correlated, and stages B and C are negatively correlated. A statistical model for such developmental profiles has to include these dependencies between subsequent stages, as it is provided by dependence trees. Let X = (X_{1}, ..., X_{ u }, ..., X_{ L }) be a Ldimensional continuous random vector where the variable X_{ u }denotes the expression values of the developmental stage u and x = (x_{1}, ..., x_{ L }) denotes a realization of X representing a developmental profile of a gene. We represent a tree by its predecessor or parent map, pa {1, ..., L} ↦ {1, ..., L} for which we assume without loss of generality that 1 < pa(u) < u and pa(1) = 1. Then we can write for the probability density function (pdf) of a conditional
We denote the model parameters by θ = (τ_{1}, ..., τ_{ u }, ... τ_{ L }) and the DTree by the tuple (X, pa, θ). Note, that a DTree can also be viewed as an approximation of the joint distribution of a Ldimensional continuous random vector by a product of L  1 second order distributions [16].
We use conditional Gaussian density functions [27] as conditional densities, denoted by p [x_{ u }x_{pa(u)}, τ_{ u }] in Eq. 2. Hence, for a given developmental profile x and a nonroot developmental stage u with pa(u) = v, the pdf takes the form
where τ_{ u }= (μ_{ u }, w_{uv}, ${\sigma}_{uv}^{2}$) are the parameters for one conditional density in the model.
For a given expression data set consisting of measurements for N genes at L developmental stages, let x_{ i }= (x_{i 1}, ..., x_{ iu }, ..., x_{ iL }) be the developmental profile of gene i, and x_{ iu }be the expression value of the gene i in development stage u for 1 ≤ i ≤ N and 1 ≤ u ≤ L. As derived in the Protocol in the Additional data file 1, the maximum likelihood estimates (MLE) for the parameters of the conditional Gaussian are
These terms can be computed from the sufficient statistics
The conditional normal distribution can be seen as estimating a linear fit between X_{ u }and X_{ v }, where w_{uv}> 0 indicates a positive linear correlation and w_{uv}< 0 a negative linear correlation between variables; w_{uv}= 0 if the variables are independent. Furthermore, w_{uv}and ${\sigma}_{uv}^{2}$ are related because the better the linear fit the smaller the variance. For the special case of the root (recall pa(1) = 1), w_{11} is set to zero, and the conditional density is effectively a univariate normal. In total, the model has 3L  1 free parameters.
A very simple, but useful application, is to query the developmental profiles from a data set with a tree model. By defining the model parameters in an interactive manner, we can compute the likelihood (Eq. 2) of all expression profiles x_{ i }. rank them accordingly, and list the m most likely profiles (see [26] for the tool description and tutorial). This interactive tool allows biological experts to find genes following a developmental profile of interest.
Returning to the example in Fig. 2, the model estimates given the tree and developmental profiles are
As expected, w_{BA}and w_{DB}are positive, indicating a linear dependence between these variables. On the other hand w_{CB}is negative.
Mixtures of dependence trees
In order to find clusters of coexpressed genes, we combine several dependence trees (DTree) in a mixture. Each DTree is a representation of a cluster or group of genes with similar developmental profiles; that is, each DTree models distinct patterns of gene expression in the course of development (see Fig. 3 for an example). The differentiation of cells is conveniently represented as a developmental tree and the structure or topology of this tree is wellknown for most data sets under investigation. Consequently, all trees in a mixture share the same topology. A mixture of dependence trees accommodates overlapping clusters while reflecting the inherent dependencies between stages. Throughout this paper we refer to the presented method as well as to the resulting model as MixDTrees.
More formally, we combine a set of K DTrees in a mixture model $f(x\Theta )={\displaystyle {\sum}_{k=1}^{K}{\alpha}_{k}p[x{\theta}_{k}]}$, where Θ = (θ_{1}, ..., θ_{ K }, α_{1}, ..., α_{ K }), θ_{ k }denotes the parameters of the kth DTree and α_{ k }is proportional to the number of developmental profiles assigned to the kth Dtree; as usual α_{ k }≥ 0 and ${\sum}_{k=1}^{K}{\mathcal{\alpha}}_{k}=1$. To avoid overfitting of the tree models, in particular for components with low component priors α_{ k }–that is, a small number of assigned genes–we propose a maximumaposteriori (MAP) approach, which regularizes the estimates from Eq. 5 and Eq. 6. Given this preferable characteristic, MAP estimates are used in all MixDTrees experiments, unless otherwise stated. Note also, the parameters of the mixture are estimated with the ExpectationMaximization (EM) algorithm [28] (see Methods section for EM and MAP details).
As stated in the introduction, the problem approached here is closely related to gene expression timecourse analysis. There is a vast amount of literature on models and clustering methods suitable for timecourses [18–22, 29, 30]. Lately, attention has been given to the fact that these timecourses have usually few time points [31, 32], a characteristic previously ignored. This aspect is also essential to our application, since the number of distinguishable developmental stages is usually small, for example at most seven in our data sets. Note that a single chain of subsequent development stages, such as the stages of Bcell differentiation in Fig. 1, is by definition a tree. While dependence trees are indeed also suitable for timecourses, the complex dependency structures necessary due to branching of the developmental tree into distinct lineages prevents the use of timecourse models, as there is no effective way of incorporating the necessary extensions into these models [19, 22]. In the context of mixtures, our method represents an alternative to the parameterization of the covariance matrix of a mixture of multivariate Gaussians [17]. With MLE, the dependence tree model essentially imputes zeros in the covariance matrix reducing the number of parameters to the order of L. If we would consider all the covariances between observations for L developmental stages; it would be straightforward to represent the data distribution by a Lvariate Gaussian model with full covariance matrix. However, the estimates for the L^{2} parameters are often unreliable even for small values of L and the parameter estimation is prone to overfit to outliers often found in noisy and scarce data. In fact, mixtures of Gaussians with full covariance matrix were outperformed by simpler parameterizations of the covariance matrices in the context of gene expression time courses [18].
Application in lymphocyte cell development
We apply our method to obtain MixDTrees for the data sets TCell, BCell, LymphoidTree, and SIM (see Methods section for details) and compare our clustering results to previous work. Our data is complemented with information from OMIM [33], the Gene Ontology database [34] and from literature. For TCell and BCell, we use the same number of clusters as Hoffmann and colleagues (20) [4, 5, 35] and for LymphoidTree we apply the BIC criterion [36] (see Fig. S4 in Additional data file 2), which also resulted in an optimal choice of 20 clusters. As discussed in Dependence trees section, a simple way to check for similarities in the expression between developmental stages is to compute the correlation matrix of the data set at hand (see Mixtures of dependence trees estimation section).
T cell development (TCell)
TCell is a gene expression data set from seven differentiation stages of the T cell development (see Methods section and Fig. 1 for details). The only branch in this tree is the final differentiation of DPS precursors into CD4 single positive SP4 cells and CD8 single positive SP8 cells. Most clusters show a distinctive pattern of differential expression along the developmental path but do not differ between SP4 and SP8 cells (clusters 4, 7, 11, 13, 14, 15, 16, 19, and 20). The most drastic changes occur at the DPL stage in which the cells are proliferating and subsequently start to rearrange the TCRαlocus. This is also reflected in the overall correlation matrix (Table S1 in Additional data file 3). Although the expression values of all neighboring stages are positively correlated, the correlation between the DPL stage and the DPS stage is much smaller in comparison to the double negative stages, all of which are relatively highly correlated. The correlation matrix suggests that SP4 and SP8 cells are more similar to each other than to their precursor DPS cells, which is expected since the two types of mature T cells share many cellular functions [4]. The largest differences with respect to SP4 and SP8 are found in clusters 5 and 18 (Fig. 4). In cluster 5, cellcycle genes are clearly enriched. In contrast, cluster 18 mainly contains regulatory proteins involved in transcription and signaling (see Fig. 4).
Hoffmann and colleagues used selforganizing maps (SOM) to cluster the expression profiles [4, 5, 35]. From now on, we refer to Hoffman and colleagues' results simply as SOM. In our analysis we observe clusters with similar developmental profiles, which we define as the average over the gene expression profiles of a cluster. As expected, there is not a onetoone relationship between the two clusterings. While the single gene profiles are similar since we used analogous normalization and filtering procedures (see Methods section), the actual gene clusterings differ (see Table S12 in Additional data file 3). An objective assessment of clustering quality on developmental data is impossible due to lack of benchmarking data. Furthermore, there is no agreement in the literature on a methodology to validate clustering results [37]. In order to demonstrate that our method is able to extract additional biological information, we concentrate our discussion on clusters of distinct developmental profiles that could not be detected by SOM [4]. For such a cluster we assign functions to genes using the GO term annotation and complementary literature. Ideally, the functions of all genes of the cluster would match the cellular processes of the particular developmental stage at which these genes are overexpressed. Additionally, if some of these genes are of unknown function then the developmental profile can help to generate hypotheses about their functional role. In our analysis we find that genes of cluster 8 are overexpressed in DN3 and DN4 cells, a developmental profile that has not been previously discovered (Fig. 4). With SOM, the genes of this cluster are dispersed over the two clusters (see Table S12 in Additional data file 3). Out of the 30 genes of cluster 8 seven are related to vesicle transport, or to the Golgi/ER system. Additionally, we find five cellcycle related genes, three involved in mitochondrial function, and seven genes of other functions, which are mainly involved in signaling. These findings agree with the functions of DN3 and DN4 cells, which is the transport of precursor receptor molecules to the cell surface membrane and the initiation of proliferation. This demonstrates that our method is able to identify functionally relevant gene sets even if the expression changes are not as large as for the DPL stage, for example. The complete results, including gene expression plots, analysis of GOterm and microRNA enrichment, can be found in our web database [26].
B cell development (BCell)
In a similar approach to the TCell study, we investigated gene expression for five consecutive stages during B cell development (see Methods section and [4, 5] for details). The correlation matrix of BCell suggests dependencies between gene expression values of successive stages, with the largest correlation between preBI and large preBII cells and between immature and mature B cells (see Table S2 in Additional data file 2). When we compare, as in the TCell set, our clustering results to those of Hoffmann and colleagues [4], we observe similar average developmental profiles although the contingency table indicates differences in the cluster compositions (Table S13 in Additional data file 3). Clusters 3, 5 and 6, for example, contain genes that are upregulated in preBI and large preBII cells and downregulated in later developmental stages (Fig. 5). Consistent with the phenotype of these cells, the function assigned to the genes of this cluster are mainly related to proliferation. GO categories that are associated with mitosis, cellcycle and chromatin remodeling are clearly overrepresented in these clusters (see our web database [26]).
Cluster 20 shows an average developmental profile that was not detected with SOM [4, 5]. The genes of this cluster are downregulated in preBI cells, in which the first rearrangement of the D^{H}and J^{H}segments on the H chain loci has taken place, and upregulated in all the following developmental stages (Fig. 5). With SOM [4], these 23 genes are found distributed over the four clusters 11, 13, 14 and 17 (Table S12 in Additional data file 3). The most palpable common function of many cluster 20 genes is the regulation of survival and apoptosis during B cell development. The gene products Nfkbia, Traf5 and the Srcfamily protein tyrosine kinases Lyn and Syk are known regulators of NFkappa B activity, which in turn has been found to be involved in B cell fate decision and survival [38–40]. Similarly, Krupellike factor 2 (Klf2) protects cells against TNFalpha induced apoptosis [41]. Furthermore, Icam2 and Rhoh, whose encoding genes are two other members of cluster 20, regulate the adhesiveness of primary B cells depending on their activation state and protect them from apoptosis [33, 42].
Lymphoid tree (LymphoidTree)
LymphoidTree combines data sets of several studies [10–12], and the resulting tree contains expression measurements from lymphoid cells of six developmental stages, namely hematopoietic stem cells, proB, preB, and immature B cells, mature SP4 T cells, and natural killer (NK) cells. This integration of data is possible because the studies were carried out on the same array platform. Although the developmental tree is far less detailed compared to TCell and BCell, we still gain insights on differences between the cell lineages. As expected, the correlation matrix shows that the expression patterns of the three B cell stages are more highly correlated among each other then expression patterns of different lineages. Moreover, the overall expression of SP4 cells and NK cells is positively correlated. The resulting clusters provide a basis to hypothesize about early developmental decisions and suggest target genes for further investigations. For example cluster 11 contains genes that are strongly upregulated in NK cells, weakly induced in the SP4 cells and not expressed in the precursor B cells (Fig. 6). Many of the cluster 11 genes are well known to be expressed in NK cells, as for example the cell surface receptor genes Cd244, Klra1, and Crtam [33, 43]. Among the lesser known genes is the one that codes for the Pu.1 related transcription factor SpiC, which has already been found to be temporarily expressed during B cell development [44]. In contrast, cluster 19 contains genes that are upregulated in SP4 cells and in all B cell precursors but not in NK cells (Fig. 6). Important functions during B and T cell maturation are reflected by genes in this cluster, like the bruton tyrosine kinase Btk, the transcription factor Pou2af1, which is involved in immunoglobulin gene regulation, and the DNA repair genes Trp53bp1 and Pnkp [33].
Simulated data (SIM)
We demonstrate with simulated data that our novel method outperforms established methods, such as SOM, kmeans and mixture of Gaussians, when inferring tree components in complex mixtures for varying levels of dependence between the individual variates. The dependence is reflected in the magnitude of w_{uv, k}(Eq. 5) of a tree. By sampling these parameters from different intervals, [ε, ε ], [0.5, 0.5], [1, 1], [1.0, 0.5] ∪ [0.5, 1] and [1, 1 + ε] ∪ [1 – ε, 1], we can create mixtures with components ranging from independent models to highly dependent ones. We generate a data set for each sampled mixture. We used MixDTrees, mixture of Gaussians, kmeans and SOM to compute clusters, which we can compare to the classes used in data generation to compute specificity and sensitivity of the clustering solutions. Method performance is evaluated with a paired ttest. Details are given in Methods section.
We observe that the MixDTrees with MAP estimates (MixDTreesMAP) have a higher specificity and sensitivity than kmeans and SOM in all experimental settings (Fig. 7 top) (pvalue below 0.005). In the (almost) independent case (w_{uv, k}∈ [ε, ε]), this is not expected, since the data agrees well with the assumptions of kmeans and SOM. This also explains the large standard deviations of MixDTreesMAP in that case. As expected, the MixDTreesMAP clearly improves the cluster recovery in settings with pronounced dependence structure, while the performance of kmeans and SOM deteriorates slightly. In comparison to others mixture model methods (Fig. 7 bottom), MixDTreesMAP also obtains a significantly higher specificity and sensitivity in almost all experimental settings. The mixture of Gaussians with diagonal covariance matrices performs well in the independent case (1), which meets the model assumptions, but it has poor results in experiments with higher dependence (pvalues below 0.05 for settings 3, 4 and 5). The mixture of Gaussians with full covariance matrix (MGFull) has a reasonable sensitivity in all settings, but very poor specificity (pvalue below 0.05 in settings 3, 4 and 5 for specificity and in all settings for specificity). The reason for these results is that MGFull tends to populate some clusters with few data points, a problem known as spurious local maxima [17]. Note that we use a MAP estimate of MGFull to mitigate this problem. Even though there are other methods for detection of spurious local maxima in MGFull, which could lead to better specificity, this would require extensions of the EM method, and consequently slower convergence [17]. On the other hand, MixDTrees, which has a lower computational running time than MGFull, achieves good results without the need of any extension. MixDTrees with MLE estimates (MixDTreesMLE) has good overall performance, but is outperformed by MixDTreesMAP in all cases, except experimental settings 1 and 5 (pvalue below 0.05 for settings 2, 3 and 4). In experimental setting 5, where data is highly dependent, by definition, both methods work similarly well. Nevertheless, such high dependency would never be found in real data sets, since noise in the data obfuscates dependencies between variables. Additionally, we performed further experiments with simulated data to evaluate the robustness of the method with respect to noise (see Additional data file 1). There, MTMAP maintains good sensitivity and specificity of cluster recovery even for high noise levels.
This demonstrates that the MixDTrees is a superior alternative to SOM and kmeans in all cases. In relation to other mixture models, MixDTrees represents a good tradeoff between a complex model class such as multivariate Gaussian with full covariance matrices and the simple Gaussian with diagonal covariance matrices. Furthermore, MAP estimates of the MixDTrees represent a more robust alternative to the MLE counterpart.
MicroRNA target discovery
LympMIR contains a set of 17 microRNAs that are potentially involved in lymphocyte cell development (for details see Methods section). It has been proposed that microRNAs bind target mRNAs specifically via base pairing, which subsequently leads to interference with the translational machinery or mRNA degradation, and thus can control whole groups of genes simultaneously [45]. Recent microarray studies have demonstrated that the microRNA expression negatively correlates with mRNA target expression in a tissue specific manner [46–48].
Having identified a cluster of coexpressed genes during lymphoid development we ask whether a certain microRNA could be a potential regulator of this cluster (see Fig. 8). For this task we first obtain lists of potential target genes for each microRNA from the miRBase Targets database [49], which contains predictions made by sequence based methods. Given our clustering results, we use the statistic of the ChiSquare Test [50] to obtain a list of microRNAs, whose potential targets are overrepresented in a cluster. This is an analogous approach to finding Gene Ontology [51] terms overrepresented in a cluster of genes. Given a set of n genes, we count the number c of genes in a given cluster, the number t of genes identified as targets for a given microRNA and the number h of genes that are both in the cluster and are targets of the microRNA. The resulting pvalue reflects the statistical significance of observing a count h, given n, c and t. A lower pvalue indicates a higher "microRNA enrichment", and, consequently, a better result. By choosing a pvalue cutoff, we can construct a list of enriched microRNAs for each cluster as well as a list of target genes related to the enriched microRNAs. Note, that the statistics for microRNAbinding are not well developed; intricate dependencies introduced by sequence similarities among the microRNAs and the target genes exist and complicate the analysis. As we also consider a manually selected set of microRNAs, we choose a somewhat relaxed pvalue cutoff, foregoing multiple testing corrections [52], followed by a careful biological evaluation. For the following discussions we restrict our result set to clusters that contain at least four target genes in total.
In summary, in TCell our target prediction scheme detects significant enrichment for eleven out of the 17 initial microRNAs in four out of the 20 clusters (Table 1). In these four clusters we detect in total 35 candidate target genes, which is a considerable reduction of the set of 229 targets that have been predicted by sequence based methods alone [49]. For BCell these numbers are respectively, eleven out of the 17 microRNAs, four out of the 20 clusters, and 29 out of the 273 predicted targets (Table 1). In particular, we find the five microRNA families miR15, miR181, miR221, miR26, and miR1423p to be enriched in both TCell and BCell by our criterion. See Table S6 in Additional data file 3 for microRNA enrichment in LymphoidTree and Table S7, Table S8, Table S9, for pvalues of microRNA enrichment of all data sets. As mentioned earlier, the BCell clusters 3, 5, and 6 show a similar expression profile. We find that cluster 5 of the results of the TCell set overlaps substantially with clusters 3 and 5 of BCell (Table 1). In TCell cluster 5 we find miR15a, miR181a, miR26a, miR24, and miR221 as potential regulators and 20 potential target genes, seven of which are also present among the 18 BCell candidate genes of clusters 3 and 5. The developmental profiles of the clusters of both lineages show strikingly analogous phenotypical features, namely upregulation in the proliferating large cell populations (DN4, DPL, large preBII) and from then on strict downregulation. In TCell cluster 5 there are eight genes and in the BCell clusters 3 and 5 there are nine target genes that are known to be involved in DNA metabolism, cellcycle and mitosis (Table 1). This suggests a regulatory role for the identified microRNAs in reducing the transcript levels of genes that are important for cell proliferation. This is supported by the fact that a similar role for microRNA was found in Drosophila germline stem cells [53].
At the individual gene level we identify some candidate microRNA targets for further detailed analysis: the three known genes (H2Eb1, Ltb, Tap2) of BCell cluster 19 are all involved in the antigen presentation by MHC class II molecules [33, 54]. In the context of the cell cycle, Chek1 (clusters TCell 5 and BCell 5) and Cdc25a (cluster TCell 5) are important for the transition between G1/S and G2/M phases [55].
Furthermore, both genes are candidate targets of the same microRNA, miR15a, which is related to apoptosis in chronic lymphoid leukemia cells [56]. Another interesting gene codes for the nuclear factor Y (Nfyb; cluster BCell 5), which regulates Hoxb4 [57], Cdc34 [58] and the major histocompatibility complex in mice [59]. These are all important genes for lymphoid development. The mRNA of the growth factor independence1 transcription factor (Gfi1; cluster TCell 10) is a potential target of miR1423p with a function in the restriction of cell proliferation and maintenance of the functional integrity of lymphocyte cells [60]. Moreover, Gfi1 is implicated in the transition from CD4/CD8 double negative to double positive T cells [61].
In order to relate our approach with [4, 5], we also perform a microRNA enrichment analysis with the results of SOM (see Table S4 and S5 in the Additional data file 3). In TCell there is little overlap between the microRNA targets, with the exception of SOM cluster 6, which is a subset of targets genes from cluster 5 from MixDTrees. We also compare the pvalues obtained by both methods in a procedure similar to the one performed in [31]. For TCell, MixDTrees results in lower pvalues in nine out of 14 microRNAs (see Fig. S5 in Additional data file 2). In BCell, gene targets found with SOM are partially a subset of the ones encountered with MixDTrees; 14 out of 24 targets genes in BCell SOM are also detected by MixDTrees (Table S5 in the Additional data file 3). For BCell, (Fig. S6 in Additional data file 2), MixDTrees obtains lower pvalues in 8 out of 14 microRNAs. Even though SOM obtains lower pvalues for microRNAs found to be enriched with both methods, MixDTrees detects seven enriched microRNA not significantly enriched in SOM. An inspection of the cumulative distribution function of these pvalues also reinforces the view that MixDTrees is more sensitive in detecting enriched microRNAs than SOM in BCell (Fig. S8 in Additional data file 2). Overall, the results suggests a higher sensitivity of MixDTreesMAP in finding groups of microRNA targets sharing similar expression patterns compared to SOM. Additionally, we performed microRNA enrichment pvalue comparison between MixDTreesMAP and MixDTreesMLE for both data sets (see Additional data file 2 Fig. S9 and S10). For TCell, MixDTreesMAP achieves a higher enrichment for nine out of 14 microRNAs; while for BCell, six out of 13 microRNAs. In summary, clusters computed according to MAP have an increased enrichment for TCell and a slightly lowered enrichment for BCell. A manual inspection of the contingency table comparing the clusters from MAP and MLE (Additional data file 3 Table S15) and in the cluster size distributions (Additional data file 2 Fig. S11) shows that MixDTreesMLE has a tendency to produce spurious, small clusters as a result of overfitting, a known disadvantage of MLE estimates [17]. Note that the resulting pvalues decrease drastically as a function of the cluster size, making a clustering which joins clusters appear preferable. Enrichment analysis should be used cautiously to compare clusterings, if the cluster size distributions are not similar, as it is the case for the MLE results. This and the results on simulated data supports our preference of MixDTreesMAP over MixDTreesMLE.
Conclusion
The regulatory processes behind cell proliferation and differentiation are of central interest to developmental biologists and clinicians alike and are frequently the focus of largescale studies to investigate gene expression along paths of differentiation. To make full use of this data in a principled manner we present a novel statistical framework which models gene expression in the course of development. By combining the dependence trees in a classical mixture model, we facilitate interactive querying and visualization of data and, more importantly, the detection of possibly overlapping clusters of coexpressed genes, which provide a basis for the identification of key players in the regulatory mechanism and their mode of action.
In particular, we detect interesting groups of genes not found by classical clustering methods such as SOM. By incorporating microRNA binding data, we show how to identify complex regulatory relationships. Compared to an analysis based only on sequence, we predict a manageable number of plausible microRNA targets. Moreover, our method offers some insights into the biological role of predicted microRNAs, by the inspection of the developmental profiles of gene targets associated with one microRNA. A comparison with SOM indicates that our approach is more sensitive for finding coexpressed genes on which the same microRNA can have a regulatory effect.
Extensions to accommodate further types of data are straightforward. Binding sites of transcription factors can be analyzed completely analogous to the microRNA analysis. If expression levels of microRNAs in developmental stages investigated in TCell or BCell were available, we could incorporate a target prediction framework [62]. Furthermore, we can simply apply established techniques [63–66] to extend our mixture model to integrate heterogeneous data–sequence information, protein interaction, genotype, phenotype data–and semisupervised extensions to mixture estimation can be applied to make use of biological knowledge about functional similarities and regulatory relationships [22, 67, 68]. This is of highest relevance, because the identification of regulatory modules is actually feasible compared to the automated inference of regulatory networks [69]. Once a statistical model is obtained, further detailed questions about the significance of differences, or the most likely stage, at which differentiation occurs can be easily answered.
Fascinating extensions are possible, even when one only considers gene expression data and the basic method. None of the currently publicly available data sets offers both a tree with a large number of branches and a detailed view of all, in particular early, development stages ([70] concentrates on mature and immature cells in final development stages); combining data from several microarray platforms suffers from the usual problems. Hence, we concentrate on two smaller but detailed studies covering several stages of T cell and B cell development [4, 5], and a tree containing three lineages of lymphoid cells. Note that in the latter several cell types of intermediary development stages are not measured. Nevertheless, our analysis indicates that our method takes advantage of the tree structure information in detecting relevant differences of gene expression in these lineages. This also reinforces the importance of the creation of expression compendia, such as the one in [70], where many intermediary stages of differentiation of the developmental tree are also present. Such data will be of great value as computational methods can exploit characteristics intrinsic to cell development.
Lastly, developmental biologists are still redrawing developmental trees with the discovery of new intermediary stages and "alternative" paths of development [1–3]; a particular developmental stage might also be formed by a mixture of distinct cell types not well characterized yet. As an example of an alternative path, there has been evidence that DN1 T cells can be originated not only from the lymphoid progenitor as depicted in Fig. 1, but also from the earlier multipotent progenitor cells [3]. It is an exciting prospect to infer branches and stages of a developmental tree from gene expression data, ideally per functional module. This structure learning (see [16] for discrete data) can be incorporated in the EMbased parameter estimation. In conclusion, our results suggest that the mixture of dependence trees provides a natural and powerful representation of developmental gene expression data. Furthermore, our results reinforce the importance of the creation of detailed and heterogeneous data sets for helping elucidate the regulatory mechanisms of development.
Methods
Data
Our work concentrates on two detailed studies covering several stages of the B and T cell development [4, 5] and a tree containing three lineages of lymphoid cells [10–12]. All gene expression data sets analyzed are deposited at the Gene Expression Omnibus [71]. Their accession entries are: GDS44 and GDS52 for BCell, GDS237 and GDS257 for TCell, and GDS1077 (HSC), GSE2227 (Bcells) and GDS828 (NK and SP4) for the LymphoidTree data. Final normalized and filtered data sets are found in [26]. Furthermore, we also use simulated data sets in order to evaluate the method. Finally, we describe a set of microRNAs that are used in our study.
T cell development (TCell)
This data set contains measurements of gene expression during the development of T cells in mouse [4]. Based on cell surface markers seven stages have been distinguished: CD4 and CD8 double negatives (DN2, DN3, DN4), large double positives (DPL), small double positives (DPS), single positive CD4 (SP4) and single positive CD8 (SP8) (see Fig. 1 for the corresponding tree, and the original publication for details [4]). Affymetrix MU11k chips with four or five replicates were used to measure the expression levels of 13,104 mouse genes. We performed variance stabilization [72] on all chips, and computed the median values of replicates. To facilitate comparisons, we restrict the set to the same list of 1318 differentially expressed genes that was used by Hoffmann and colleagues [4]. Furthermore, we normalize the expression levels separately for each gene to mean zero and standard deviation one, as is routine in gene expression analysis. Finally, we map each probe set to a gene symbol if it exists in the respective chip platform annotation provided by the GEO database [73]. The final dataset is found at Additional data file 4.
B cell development (BCell)
This data set contains expression levels of five consecutive stages of the B cell lineage, PreBI, large PreBII, small PreBII, immature, and mature B cells [5]. This study was conducted on Affymetrix MU11k chips also, and we preprocess the data exactly as it is described for TCell. The final dataset is found at Additional data file 5.
Lymphoid tree (LymphoidTree)
We combine the data of the wildtype control measurements of three studies [10–12] based on the Affymetrix U74 platform to obtain a development tree with distinct lymphoid lineages. This results in expression values of a hematopoietic stem cell (pHSC) from [10], of Natural Killer cells (NK) and of SP4 cells from [11], and of three B cell stages from [12], which are proB, preB and immature B cells. We preprocess the data exactly as it is described for TCell. Additionally, we remove genes which do not display a 2fold change in expression at least once. The final dataset is found at Additional data file 6.
Simulated data (SIM)
We use MixDTrees with random parameterizations to generate simulated data. For the tree structure given in Fig. 2, we randomly chose the μ_{uv, k}from the range [1.5, 1.5] and ${\sigma}_{uv,k}^{2}$ from [0, 1]. We create five experimental settings to inspect the performance of the method in the presence of distinct levels of dependence. For these five settings, we sample w_{uv, k}uniformly from [ε, ε ] (independent data), [0.5, 0.5], [1, 1], [1.0, 0.5] ∪[0.5, 1] and [1, 1 + ε ] ∪ [1 – ε, 1] (tree dependent data) respectively, where ε = 0.01. We chose K = 5 and mixture coefficients equal to α = (0.1, 0.15, 0.2, 0.2, 0.35). For each experimental setting, we generate ten such mixtures, and sample 500 development profiles from each (see Additional data file 1 for more results on simulated data and Additional data file 7 for datasets). To evaluate the results we compare the class information from the data generation to compute sensitivity, $\frac{\#TP}{\#TP+\#FN}$, and specificity, $\frac{\#TP}{\#TP+\#FP}$, where, for a given clustering result and the class information, TP denotes the number of pairs of objects in the same cluster and same class. The remaining three types of pairs are counted as FP (same cluster, distinct class), TN (distinct cluster and class) and FN (distinct cluster, same class). For each method, we compute the sensitivity and specificity on all 10 data sets of an experimental setting and take the mean (see Fig. 7). To compare MixDTreesMAP with other methods, we apply a one tailed paired ttest to evaluate the hypothesis that two methods have the same mean specificity (or sensitivity) in a given experimental setting. Low pvalues indicate that the equal means hypothesis was rejected and that mean specificity (or sensitivity) was significantly higher in MixDTreesMAP. For brevity, in the Simulated data section, we simply state–MixDTreesMAP had a higher sensitivity than method X (pvalue below 0.05)–when the null hypothesis is rejected.
Lymphoid development related microRNAs (LympMIR)
We collect 17 microRNAs that have been found to be involved in Lymphoid development or at least differentially expressed between distinguishable lymphocyte cell types [7–9, 56, 74]: mmumiR24, mmumiR26a, mmumiR1423p, mmumiR146, mmumiR150, mmumiR155, mmumiR181a, mmumiR181b, mmumiR181c, mmumiR191, mmumiR221, mmumiR222, mmumiR223 and mmumiR342. Additionally, we include mmumiR15a, mmumiR15b, and mmumiR16 because, according to recent papers, they participate in the regulation of cell proliferation and apoptosis [75, 76]. Since we refer exclusively to microRNAs of the mouse in this work, the species prefix mmu is omitted throughout the text. The lists of candidate targets of these were obtained in the miRBase Targets database [49] (Version 2.0), which uses the Miranda algorithm [77] to search for possible microRNA binding sites in the gene sequences.
Mixtures of dependence trees estimation
We combine K DTrees in a mixture $f(x\Theta )={\displaystyle {\sum}_{k=1}^{K}{\mathcal{\alpha}}_{k}p[x{\theta}_{k}]}$, where Θ = (θ_{1}, ..., θ_{ K }, α_{1}, ..., α_{ K }), θ_{ k }denotes the parameter set of the k th Dtree and αk ≥ 0, ${\sum}_{k=1}^{K}{\alpha}_{k}=1$, are the mixture weights or component priors. By introducing a discrete hidden variable Y = {y_{ i }} for 1 ≤ i ≤ N, which indicates which DTree generated which developmental profile x_{ i }, we can formulate a complete loglikelihood function and estimate the parameters with the ExpectationMaximization (EM) algorithm [28]. Given an initial parameterization Θ^{0}, EM iterates two steps: first estimating the posterior probabilities $P[{y}_{i}=k{x}_{i},{\theta}_{k}^{m}]$ (E Step), and second the computation of the maximumlikelihood parameters Θ^{m + 1}(Mstep), as defined in Eq. 4, Eq. 5 and Eq. 6. We refer the reader to [36] for details of the EMalgorithm.
To avoid overfitting the models, in particular for components with low component priors α_{ k }–that is, a small number of assigned genes–we propose maximumaposteriori (MAP) approach. We assume that w_{uv, k}~ N(0, α_{ k }β_{uv, k}, ${\sigma}_{uk}^{2}$) [78]. Consequently, the estimates take the form.
For the sake of simplicity we omit the coefficients k which indicates a tree in a given mixture from formulas in the Dependence tree section. See Protocol for exact MLE and MAP formulas in the mixture context. When β → ∞, we obtain a noninformative prior, for which the MAP and MLE estimates are equal. As β → 0, w → 0 and we have a univariate Gaussian. As in [78], we use a empirical Bayes approach to estimate the value of the hyperparameter β_{uv, k}as
where r_{ ik }is equal to the posterior probability P [y_{ i }= kx_{ i }, θ_{ k }] calculated in the E step. This term can be interpreted as the inverse of the linearity evidence. It penalizes components with low responsibilities and larger variances, enforcing lower w_{uv, k}values (see Protocol in Additional data file 1 for derivations of all formulas).
The last step after the mixture estimation is the assignment of genes to groups. This is done by assigning genes to the component that maximizes the posterior of the ith gene, which is y_{ i }= argmax_{1 ≤ k ≤ K}(r_{ ik }). Note, that more refined assignment schemes [22] (i.e., decoding a mixture) which increase the robustness of the clustering method can also be used.
Application in lymphoid development
We perform the following steps on each of the sets TCell, BCell, LymphoidTree, and SIM. The mixture estimation method is initialized with K random DTrees, which are obtained by choosing random values uniformly and in [0, 1] independently for each r_{ ik }and estimating the individual models. Subsequently, we train the mixture model using the EMalgorithm and MAP estimates. To avoid the effect of the initialization, all estimations are repeated 15 times, and the one with highest likelihoods is selected (a similar procedure is applied for kmeans and SOM). The implementation of our method (licensed under the GPL) and MS Windows binaries are available at [26]. There you can also find a web database–generated with our MixDTrees Report tool–with results of all analyses described in this article.
On TCell and BCell, we used the SOM results as given by [4, 5]. For SOM experiments on SIM data, we used the default parameters of the implementation [25], which uses a set of heuristics to select the values. Furthermore, we performed a clustering of SOM nodes with kmeans as it is a common practice [79]. In order to facilitate the comparison between our clustering results and the clusters of the original work we reorder our clusters accordingly. Dependence between developmental stages is measured as the correlation between variables. Given two stages, X_{ u }and X_{ v }the correlation is defined as
where 1 ≤ ρ_{u, v}≤ 1 and ρ_{u, v}= 0 indicates independence of variables.
Abbreviations
 BCell:

B cell development data
 DTree:

dependence tree
 DN:

CD4/CD8 double negative cells
 DPL:

CD4+/CD8+ double positive large cells
 DPS:

CD4+/CD8+ double positive small cells
 FACS:

fluorescence activated cell sorting
 LympMIR:

hematopoiesis related microRNAs data
 LymphoidTree:

lymphoid tree data
 MAP:

maximumaposteriori
 MLE:

maximum likelihood estimates (MLE)
 MixDTrees:

mixtures of dependence trees
 MixDTreesMAP:

mixtures of dependence trees with MAP estimates
 MixDTreesMLE:

mixtures of dependence trees with MLE estimates
 NK:

natural killer cells
 pHSC – pluripotent:

selfrenewing hematopoietic stem cells
 SIM:

simulated data
 SOM:

selforganizing maps
 SP4:

single positive CD4
 SP8:

single positive CD8
 TCell:

T cell development data
References
 1.
Matthias P, Rolink AG: Transcriptional networks in developing and mature B cells. Nat Rev Immunol. 2005, 5 (6): 497508.
 2.
Rothenberg EV, Taghon T: Molecular genetics of T cell development. Annu Rev Immunol. 2005, 23: 601649.
 3.
Bhandoola A, Sambandam : From stem cell to T cell: one route or many?. Nature Reviews Immunology. 2006, 6: 117126.
 4.
Hoffmann R, Bruno L, Seidl T, Rolink A, Melchers F: Rules for gene usage inferred from a comparison of largescale gene expression profiles of T and B lymphocyte development. J Immunol. 2003, 170 (3): 13391353.
 5.
Hoffmann R, Seidl T, Neeb M, Rolink A, Melchers F: Changes in gene expression profiles in developing B cells of murine bone marrow. Genome Res. 2002, 12: 98111.
 6.
Warren LA, Rothenberg EV: Regulatory coding of lymphoid lineage choice by hematopoietic transcription factors. Curr Opin Immunol. 2003, 15 (2): 166175.
 7.
Chen CZ, Li L, Lodish HF, Bartel DP: MicroRNAs modulate hematopoietic lineage differentiation. Science. 2004, 303 (5654): 8386.
 8.
Monticelli S, Ansel KM, Xiao C, Socci ND, Krichevsky AM, Thai TH, Rajewsky N, Marks DS, Sander C, Rajewsky K, Rao A, Kosik KS: MicroRNA profiling of the murine hematopoietic system. Genome Biol. 2005, 6 (8): R71
 9.
Ramkissoon SH, Mainwaring LA, Ogasawara Y, Keyvanfar K, McCoy JP, Sloand EM, Kajigaya S, Young NS: Hematopoieticspecific microRNA expression in human cells. Leuk Res. 2005
 10.
Bystrykh L, Weersing E, Dontje B, Sutton S, Pletcher MT, Wiltshire T, Su AI, Vellenga E, Wang J, Manly KF, Lu L, Chesler EJ, Alberts R, Jansen RC, Williams RW, Cooke MP, de Haan G: Uncovering regulatory pathways that affect hematopoietic stem cell function using 'genetical genomics'. Nat Genet. 2005, 37 (3): 225232.
 11.
Poirot L, Benoist C, Mathis D: Natural killer cells distinguish innocuous and destructive forms of pancreatic islet autoimmunity. Proc Natl Acad Sci USA. 2004, 101 (21): 81028107.
 12.
Tze LE, Schram BR, Lam KP, Hogquist KA, Hippen KL, Liu J, Shinton SA, Otipoby KL, Rodine PR, Vegoe AL, Kraus M, Hardy RR, Schlissel MS, Rajewsky K, Behrens TW: Basal immunoglobulin signaling actively maintains developmental stage in immature B cells. PLoS Biology. 2005, 3 (3): e82
 13.
Rangel C, Angus J, Ghahramani Z, Lioumi M, Sotheran E, Gaiba A, Wild DL, Falciani F: Modeling Tcell activation using gene expression profiling and statespace models. Bioinformatics. 2004, 20 (9): 13611372. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/20/9/1361]
 14.
Troncale S, Tahi F, Campard D, Vannier JP, Guespin J: Modeling and simulation with Hybrid Functional Petri Nets of the role of interleukin6 in human early haematopoiesis. Pac Symp Biocomput. 2006 : 427438.
 15.
Basso K, Margolin AA, Stolovitzky G, Klein U, DallaFavera R, Califano A: Reverse engineering of regulatory networks in human B cells. Nat Genet. 2005, 37 (4): 382390.
 16.
Chow C, Liu C: Approximating discrete probability distributions with dependence trees. IEEE Trans Info Theory. 1968, 14 (3): 462467.
 17.
McLachlan GJ, Bean RW, Peel D: A mixture modelbased approach to the clustering of microarray expression data. Bioinformatics. 2002, 18 (3): 413422.
 18.
Yeung KY, Fraley C, Murua A, Raftery AE, Ruzzo WL: Modelbased clustering and data transformations for gene expression data. Bioinformatics. 2001, 17 (10): 977987.
 19.
BarJoseph Z, Gerber GK, Gifford DK, Jaakkola TS, Simon I: Continuous representations of timeseries gene expression data. J Comput Biol. 2003, 10 (3–4): 341356.
 20.
Luan Y, Li H: Clustering of timecourse gene expression data using a mixedeffects model with Bsplines. Bioinformatics. 2003, 19 (4): 474482. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/19/4/474]
 21.
Ramoni MF, Sebastiani P, Kohane IS: Cluster analysis of gene expression dynamics. Proc Natl Acad Sci USA. 2002, 99 (14): 91219126.
 22.
Schliep A, Costa IG, Steinhoff C, Schonhuth A: Analyzing Gene Expression TimeCourses. IEEE/ACM Trans Comput Biol Bioinform. 2005, 2 (3): 179193.
 23.
Meila M, Jordan MI: Learning with mixtures of trees. J Mach Learn Res. 2001, 1: 148.
 24.
Beerenwinkel N, Rahnenfuhrer J, Daumer M, Hoffmann D, Kaiser R, Selbig J, Lengauer T: Learning multiple evolutionary pathways from crosssectional data. RECOMB '04: Proceedings of the eighth annual international conference on Research in computational molecular biology. 2004, New York, NY, USA: ACM Press, 3644.
 25.
Vesanto J, Himberg J, Alhoniemi E, Parhankangas J: SOM Toolbox for Matlab. Tech rep. 2000, [http://citeseer.ist.psu.edu/vesanto00som.html]
 26.
Supplementary Material. [http://algorithmics.molgen.mpg.de/Supplements/ExpLym/]
 27.
Lauritzen SL, Spiegelhalter DJ: Local computations with probabilities on graphical structures and their application to expert systems. J Royal Statis Soc B. 1988, 50: 157224.
 28.
Dempster A, Laird N, Rubin D: Maximum likelihood from incomplete data via the EM algorithm. JRSSB. 1977, 39: 138.
 29.
BarJoseph Z: Analyzing time series gene expression data. Bioinformatics. 2004, 20 (16): 2493503.
 30.
Yeung KY, Medvedovic M, Bumgarner RE: Clustering geneexpression data with repeated measurements. Genome Biol. 2003, 4 (5): R34
 31.
Ernst J, Nau GJ, BarJoseph Z: Clustering short time series gene expression data. Bioinformatics. 2005, 21 (suppl 1): i159168. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/21/suppl_1/i159]
 32.
MollerLevet C, Klawonn F, Cho K, Wolkenhauer O: Fuzzy clustering of short timeseries and unevenly distributed sampling points. Advances in Intelligent Data Analysis V, Lecture Notes in Computer Science. 2003, Springer Verlag, 2810: 330340.
 33.
Johns Hopkins University M Baltimore: Online Mendelian Inheritance in Man, OMIM (TM). World Wide Web. 2006, [http://www.ncbi.nlm.nih.gov/omim/]
 34.
Ashburner M: Gene Ontology: tool for the unification of biology. Nat Genet. 2000, 25: 2529.
 35.
Hoffmann R, Melchers F: A genomic view of lymphocyte development. Curr Opin Immunol. 2003, 15 (3): 239245.
 36.
McLachlan G, Peel D: Finite Mixture Models. 2000, Wiley Series in Probability and Statistics, Wiley, New York
 37.
D'haeseleer P: How does gene expression clustering work?. Nat Biotechnol. 2005, 23 (12): 14991501.
 38.
Aizawa S, Nakano H, Ishida T, Horie R, Nagai M, Ito K, Yagita H, Okumura K, Inoue J, Watanabe T: Tumor necrosis factor receptorassociated factor (TRAF) 5 and TRAF2 are involved in CD30mediated NFkappaB activation. J Biol Chem. 1997, 272 (4): 20422045.
 39.
Nakano H, Sakon S, Koseki H, Takemori T, Tada K, Matsumoto M, Munechika E, Sakai T, Shirasawa T, Akiba H, Kobata T, Santee SM, Ware CF, Rennert PD, Taniguchi M, Yagita H, Okumura K: Targeted disruption of Traf5 gene causes defects in CD40 and CD27mediated lymphocyte activation. PNAS. 1999, 96 (17): 98039808. [http://www.pnas.org/cgi/content/abstract/96/17/9803]
 40.
Hoffmann A, Levchenko A, Scott ML, Baltimore D: The IkappaBNFkappaB signaling module: temporal control and selective gene activation. Science. 2002, 298 (5596): 12411245.
 41.
Glynne R, Ghandour G, Rayner J, Mack DH, Goodnow CC: Blymphocyte quiescence, tolerance and activation as viewed by global gene expression profiling on microarrays. Immunol Rev. 2000, 176: 216246.
 42.
Perez OD, Kinoshita S, Hitoshi Y, Payan DG, Kitamura T, Nolan GP, Lorens JB: Activation of the PKB/AKT pathway by ICAM2. Immunity. 2002, 16: 5165.
 43.
Boles KS, Barchet W, Diacovo T, Cella M, Colonna M: The tumor suppressor TSLC1/NECL2 triggers NKcell and CD8+ Tcell responses through the cellsurface receptor CRTAM. Blood. 2005, 106 (3): 779786.
 44.
Carlsson R, Hjalmarsson A, Liberg D, Persson C, Leanderson T: Genomic structure of mouse SPIC and genomic structure and expression pattern of human SPIC. Gene. 2002, 299 (1–2): 271278.
 45.
Bartel DP: MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004, 116 (2): 281297.
 46.
Legendre M, Ritchie W, Lopez F, Gautheret D: Differential Repression of Alternative Transcripts: A Screen for miRNA Targets. PLoS Computational Biology. 2006, 2 (5): e43
 47.
Lim LP, Lau NC, GarrettEngele P, Grimson A, Schelter JM, Castle J, Bartel DP, Linsley PS, Johnson JM: Microarray analysis shows that some microRNAs downregulate large numbers of target mRNAs. Nature. 2005, 433 (7027): 769773.
 48.
Sood P, Krek A, Zavolan M, Macino G, Rajewsky N: Celltypespecific signatures of microRNAs on target mRNA expression. PNAS. 2006, 103 (8): 27462751. [http://www.pnas.org/cgi/content/abstract/103/8/2746]
 49.
GriffithsJones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ: miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006, D140D144. 34 Database
 50.
Sokal FR, Rohlf : Biometry. 1995, New York: W. H. Freeman and Company
 51.
Beissbarth T, Speed TP: GOstat: find statistically overrepresented Gene Ontologies within a group of genes. Bioinformatics. 2004, 20 (9): 14641465.
 52.
Westfall P, Young S: ResamplingBased Multiple Testing: Examples and Methods for pValue Adjustment. 1993, WileyInterscience
 53.
Hatfield SD, Shcherbata HR, Fischer KA, Nakahara K, Carthew RW, RuoholaBaker H: Stem cell division is regulated by the microRNA pathway. Nature. 2005, 435 (7044): 974978.
 54.
Powis SH, Mockridge I, Kelly A, Kerr LA, Glynne R, Gileadi U, Beck S, Trowsdale J: Polymorphism in a second ABC transporter gene located within the class II region of the human major histocompatibility complex. Proc Natl Acad Sci USA. 1992, 89 (4): 14631467.
 55.
Busino L, Donzelli M, Chiesa M, Guardavaccaro D, Ganoth D, Dorrello NV, Hershko A, Pagano M, Draetta GF: Degradation of Cdc25A by betaTrCP during S phase and in response to DNA damage. Nature. 2003, 426 (6962): 8791.
 56.
Cimmino A, Calin GA, Fabbri M, Iorio MV, Ferracin M, Shimizu M, Wojcik SE, Aqeilan RI, Zupo S, Dono M, Rassenti L, Alder H, Volinia S, Liu CG, Kipps TJ, Negrini M, Croce CM: miR15 and miR16 induce apoptosis by targeting BCL2. Proc Natl Acad Sci USA. 2005, 102 (39): 1394413949.
 57.
Gilthorpe J, Vandromme M, Brend T, Gutman A, Summerbell D, Totty N, Rigby PWJ: Spatially specific expression of Hoxb4 is dependent on the ubiquitous transcription factor NFY. Development. 2002, 129 (16): 38873899.
 58.
Radomska HS, Satterthwaite AB, Taranenko N, Narravula S, Krause DS, Tenen DG: A nuclear factor Y (NFY) site positively regulates the human CD34 stem cell gene. Blood. 1999, 94 (11): 37723780.
 59.
Zhu XS, Linhoff MW, Li G, Chin KC, Maity SN, Ting JP: Transcriptional scaffold: CIITA interacts with NFY, RFX, and CREB to cause stereospecific regulation of the class II major histocompatibility complex promoter. Mol Cell Biol. 2000, 20 (16): 60516061.
 60.
Kazanjian A, Gross EA, Grimes HL: The growth factor independence1 transcription factor: New functions and new insights. Crit Rev Oncol Hematol. 2006
 61.
Schmidt T, Karsunky H, Rodel B, Zevnik B, Elsasser HP, Moroy T: Evidence implicating Gfi1 and Pim1 in preTcell differentiation steps associated with betaselection. EMBO J. 1998, 17 (18): 53495359.
 62.
Huang J, Morris Q, Frey B: Detecting MicroRNA Targets by Linking Sequence, MicroRNA and Gene Expression Data. Lect. Notes . Comput. Sci. 2006, 3909: 11429. [http://www.springerlink.com/openurl.asp?genre=article&id=doi:10.1007/11732990_11]
 63.
Segal E, Yelensky R, Koller D: Genomewide discovery of transcriptional modules from DNA sequence and gene expression. Bioinformatics. 2003, 19 (Suppl 1): i273i282. [http://bioinformatics.oxfordjournals.org/cgi/content/abstract/19/suppl_1/i273]
 64.
Segal E, Wang H, Koller D: Discovering molecular pathways from protein interaction and gene expression data. Bioinformatics. 2003, 19 (Suppl 1): i264i271.
 65.
Segal E, Shapira M, Regev A, Pe'er D, Botstein D, Koller D, Friedman N: Module networks: identifying regulatory modules and their conditionspecific regulators from gene expression data. Nat Genet. 2003, 34 (2): 166176.
 66.
Yeang CH, Jaakkola T: Time Series Analysis of Gene Expression and Location Data. Third IEEE Symposium on BioInformatics and BioEngineering (BIBE'03). 2003, 305
 67.
Bair E, Tibshirani R: Semisupervised methods to predict patient survival from gene expression data. PLoS Biol. 2004, 2 (4): E108
 68.
Pan W, Shen X, Jiang A, Hebbel RP: Semisupervised learning via penalized mixture model with application to microarray sample classification. Bioinformatics. 2006, 22 (19): 23882395.
 69.
Friedman N: Inferring cellular networks using probabilistic graphical models. Science. 2004, 303 (5659): 799805.
 70.
Hyatt G, Melamed R, Park R, Seguritan R, Laplace C, Poirot L, Zucchelli S, Obst R, Matos M, Venanzi E, Goldrath A, Nguyen L, Luckey J, Yamagata T, Herman A, Jacobs J, Mathis D, Benoist C: Gene expression microarrays: glimpses of the immunological genome. Nat Immunol. 2006, 7 (7): 686691.
 71.
Gene Expression Omnibus. [http://www.ncbi.nlm.nih.gov/projects/geo/]
 72.
Huber W, Heydebreck AV, Sultmann H, Poustka A, Vingron M: Variance stabilization applied to microarray data calibration and to the quantification of differential expression. Bioinformatics. 2002, 18 (Suppl 1): S96S104.
 73.
Edgar R, Domrachev M, Lash AE: Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucl Acids Res. 2002, 30: 207210. [http://nar.oxfordjournals.org/cgi/content/abstract/30/1/207]
 74.
Felli N, Fontana L, Pelosi E, Botta R, Bonci D, Facchiano F, Liuzzi F, Lulli V, Morsilli O, Santoro S, Valtieri M, Calin GA, Liu CG, Sorrentino A, Croce CM, Peschle C: MicroRNAs 221 and 222 inhibit normal erythropoiesis and erythroleukemic cell growth via kit receptor downmodulation. Proc Natl Acad Sci USA. 2005, 102 (50): 1808118086.
 75.
Calin GA, Liu CG, Sevignani C, Ferracin M, Felli N, Dumitru CD, Shimizu M, Cimmino A, Zupo S, Dono M, Dell'Aquila ML, Alder H, Rassenti L, Kipps TJ, Bullrich F, Negrini M, Croce CM: MicroRNA profiling reveals distinct signatures in B cell chronic lymphocytic leukemias. PNAS. 2004, 101 (32): 1175511760. [http://www.pnas.org/cgi/content/abstract/101/32/11755]
 76.
Croce CM, Calin GA: miRNAs, cancer, and stem cell division. Cell. 2005, 122: 67.
 77.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS: MicroRNA targets in Drosophila. Genome Biol. 2003, 5: R1
 78.
Minka TP: Bayesian Linear Regression. Tech rep. 2001, MIT
 79.
Vesanto J, Alhoniemi E: Clustering of the SelfOrganizing Map. IEEE Transactions on Neural Networks. 2000, 11 (3): 586[http://citeseer.ist.psu.edu/article/vesanto00clustering.html]
Acknowledgements
We would like to express our gratitude to Fritz Melchers and Roland Krause (MPI for Infection Biology, Berlin) for helpful discussions, encouragement, and valuable comments about the manuscript. We also thank Christoph Hafemeister for his work on the software, and Benjamin Georgi and Ruben Schilling for revising the manuscript. The first author would like to acknowledge funding from the CNPq(Brazil)/DAAD.
Author information
Affiliations
Corresponding author
Additional information
Competing interests
The author(s) declares that there are no competing interests.
Authors' contributions
IC implemented the approach and performed the experiments. IC and SR evaluated the results. IC, SR and AS designed this study and wrote the manuscript. All authors read and approved the final manuscript.
Ivan G Costa, Stefan Roepcke contributed equally to this work.
Electronic supplementary material
12865_2007_139_MOESM1_ESM.pdf
Additional data file 1: Protocol. This file contains information on software implementations, derivations of estimation formulas and additional experiments with simulated data. (PDF 121 KB)
Supplementary Figures
Additional data file 2: . Figures 1, 2 and 3 contains all clusters results from MixDTrees on BCell, TCell and LymphoidTree, and Figure 4 contains BIC results from LymphoidTree. Figures 5 and 6 contain comparisons between microRNA enrichment with MixDTreesMAP and SOM in TCell and BCell, Figures 7 and 8 depict the empirical cumulative distribution function (cdf) of microRNA enrichment pvalues from TCell and BCell, and Figures 9 and 10 contain comparisons between microRNA enrichment with MixDTreesMAP and MixDTreesMLE in TCell and BCell. Figure 11 describes the cluster size distribution of clustering results in TCell and BCell. (PDF 430 KB)
Supplementary Tables
Additional data file 3: . Tables 1, 2 and 3 contains correlation matrices from BCell, TCell and LymphoidTree datasets; Tables 4, 5 and 6 contains enriched microRNA and gene targets from SOM results on TCell and BCell and from MixDTreesMAP results on LymphoidTree; Tables 7, 8, 9 contains microRNA enrichment pvalues for BCell, TCell and LymphoidTree on MixDTreesMAP results; Tables 10 and 11 contains microRNA enrichment pvalues for BCell and TCell on SOM results; Tables 12 and 13 contain the contingency tables comparing clusters from MixDTreesMAP and SOM with BCell and TCell datasets; and Tables 14 and 15 contain the contingency tables comparing clusters from MixDTreesMAP and MixDTreesMLE with BCell and TCell datasets. (PDF 67 KB)
12865_2007_139_MOESM4_ESM.txt
Additional data file 4: TCell Dataset. Data set after filtering and normalization procedures. The second column indicates the cluster assignment found by the MixDTrees. (TXT 89 KB)
12865_2007_139_MOESM5_ESM.txt
Additional data file 5: BCell Dataset. Data set after filtering and normalization procedures. The second column indicates the cluster assignment found by the MixDTrees. (TXT 51 KB)
12865_2007_139_MOESM6_ESM.txt
Additional data file 6: LymphoidTree Dataset. Data set after filtering and normalization procedures. The second column indicates the cluster assignment found by the MixDTrees. (TXT 79 KB)
12865_2007_139_MOESM7_ESM.zip
Additional data file 7: SIM Datasets. Data sets from simulated MixDTrees. See readme.txt for file descriptions. The first column indicates the true label of the sample. (ZIP 358 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Costa, I.G., Roepcke, S. & Schliep, A. Gene expression trees in lymphoid development. BMC Immunol 8, 25 (2007). https://doi.org/10.1186/14712172825
Received:
Accepted:
Published:
Keywords
 Additional Data File
 Dependence Tree
 microRNA Target
 Developmental Profile
 Developmental Tree