BACKGROUND
genetic markers such as snps are the key to dissecting the genetic susceptibility of common complex diseases, such as asthma, diabetes, atherosclerosis and some cancers  <cit> .
the purpose is identifying combinations of genetic determinants which should accumulate among affected subjects.
generally, in such combinations, each genetic variant only exerts a modest impact on the observed phenotype, whereas, in contrast, the interaction between genetic variants and, possibly, environmental factors is determinant.
decreasing genotyping costs now enable the generation of hundreds of thousands of snps, spanning the whole human genome, across cohorts of cases and controls.
this scaling up to genome-wide association studies  makes the analysis of high-dimensional data a hot topic  <cit> .
despite recent technological advances and extensive research effort, the genetic basis of the aforementioned diseases remains to a large extent unknown.
yet, the search for associations between single snps and the variable describing case/control status requires carrying out a large number of statistical tests.
since snp patterns, rather than single snps, are likely to be determinant for complex diseases, a high rate of false positives as well as a perceptible statistical power decrease, not to mention intractability, are severe issues to be overcome.
the simplest type of genetic polymorphism, single nucleotide polymorphism , involves only one nucleotide change, which occurred generations ago within the dna sequence.
to fix ideas, we emphasize that one single individual can be uniquely defined by only  <dig> to  <dig> independent snps and unrelated individuals differ in about  <dig> % of their  <dig>  billion nucleotides  <cit> .
compared with other kinds of dna markers, snps are appealing because they are abundant, genetically stable and amenable to high-throughput automated analysis.
consistently, advances in high-throughput snp genotyping technologies lead the way to various down-stream analyses, including gwass.
exploiting the existence of statistical dependences between neighboring snps, also called linkage disequilibrium , is the key to association study achievement  <cit> .
indeed, a causal variant  may not be a snp.
for instance, insertions, deletions, inversions and copy-number polymorphisms may be causative of disease susceptibility.
nevertheless, a well-designed study will have a good chance of including one or more snps that are in strong ld with a common causal variant.
in the latter case, indirect association with the phenotype, say affected/unaffected status, will be revealed .
interestingly, ld also offers solutions to reduce data dimensionality in gwass.
in the human genome, ld is highly structured into the so-called "haplotype block structure"  <cit> : regions where statistical dependences between contiguous markers  are high alternate with shorter regions characterized by low statistical dependences .
the most likely explanation of this phenomenon is related to the presence of large regions with low recombination rates separated by recombination hotspots   <cit> .
relying on this feature, various approaches were proposed to achieve data dimensionality reduction: testing association with haplotypes   <cit> , partitioning the genome according to spatial correlation  <cit> , selecting snps informative about their context, or snp tags  <cit>  .
recent methods, such as haplobuild  <cit> , have permitted to construct more biologically relevant haplotypes where the "haplotype cluster structure", instead of the "haplotype block structure", is assumed: haplotypes are not constrained by contiguous orientation.
unfortunately, these methods do not take into account all existing dependences since they miss higher-order dependences.
actually, these methods do not consider the fuzzy nature of ld: the ld block boundaries are not accurately defined over the genome .
due to their ability to represent conditional independences between variables, probabilistic graphical models  offer an adapted framework for an accurate modeling of dependences between snps.
a pgm is a probabilistic model relying on a graph representing conditional independences within a set of random variables.
inherently, this model simplifies the description of the joint distribution of the set of variables.
several subclasses of pgms exist such as markov random fields  and bayesian networks .
the main difference between these two subclasses remains in the nature of the graph: in contrast with mrfs, bayesian networks are directed graphs.
although the observed variables  are often sufficient to describe their joint distribution, sometimes, additional unobserved variables, also named latent variables , have a role to play.
only few research works have been dedicated to snp dependence modeling through pgms.
a hard task because of high data dimensionality, tackling this modeling issue through pgms nevertheless offers an attractive lead.
approaches based both on mrfs  <cit>  and bns have been designed.
regarding the latter, some methods only consider observed variables  <cit>  whereas other models include latent variables  <cit> .
in particular, hierarchical bns are the most promising models for ld representation: their hierarchical structure supported by lvs allows flexible information synthesis, thus efficiently reducing the data dimensionality.
to our knowledge, modeling ld through hierarchical bns in order to reduce snp data dimensionality has not yet been designed.
notably, scalability remains a crucial issue for gwass.
in this paper, we emphasize the interest of using a forest of hierarchical latent class models , to reduce the dimension of the data to be further submitted to statistical analyses devoted to the discovery of genetic factors potentially involved in the disease.
such studies encompass single-snp analysis  <cit> , multiple-snp analysis  <cit> , snp-snp interaction analysis  <cit>  and analysis integrating gene expression  <cit> .
an fhlcm is a hierarchical bn with discrete observed and latent variables.
basically, latent variables capture the information borne by underlying markers.
in their turn, latent variables are clustered into groups and, if relevant, such groups are subsequently subsumed by additional latent variables.
iterating this process yields a hierarchical structure.
first, the great advantage to gwass is that further statistical analyses can be chiefly performed on latent variables.
thus, a reduced number of variables will be examined.
second, a model based on a hierarchical structure provides a flexible data mining tool.
for example, different degrees of data dimensionality reduction are available to the statistician.
moreover, the hierarchical structure is meant to efficiently conduct refined association testing: zooming in through narrower and narrower regions in search for a stronger association with the disease ends pointing out the potential markers of interest.
however, most algorithms dedicated to the learning of hierarchical latent class models  fail the scalability criterion when the data describe thousands of variables and a few hundreds of individuals.
in a previous work-on progress paper  <cit> , we designed an algorithm devoted to learning fhlcms.
this algorithm was named cfhlc, which stands for construction of forests of hierarchical latent class models.
the contribution brought in the present extended version is the following:  we advocate the use of fhlcms to model ld;  we provide a detailed description of the main concepts underlying our approach;  using real data, we show that the fhlcm graph is representative of the haplotype cluster structure;  in addition, we compare the haplotype cluster structure obtained through cfhlc with those output by four other algorithms;  relying on both real and simulated data, we demonstrate the ability of fhlcms to concisely model snp dependences, showing that the multiple layers of the model can take into account different ld degrees and haplotype diversity;  finally, we present a thorough study focused on both scalability and impact of adjustment of the input parameters of cfhlc algorithm.
as a prerequisite to further understanding, section preliminaries provides an informal definition of bayesian networks, focusing on latent class models and hierarchical latent class models.
then the section dedicated to the state of the art first points out the few anterior works devoted to hlcm learning in general.
this section ends with a short review of the few attempts to implement probabilistic graphical models for the specific purpose of ld modeling.
section methods motivates the modeling of ld through fhlcms and informally describes such models.
then, the focus is set on the general outline of the method proposed for fhlcm learning.
the next section depicts the sketch of algorithm cfhlc.
the last section is dedicated to experimental results and discussion.
in this section, we first test and discuss the ability of fhlcms to accurately represent the haplotype cluster structure of genetic data.
then, we compare our algorithm to other methods with respect to faithfulness in ld modeling and data dimension reduction.
we end the section with a thorough study centered on scalability and influence of the input parameters of the cfhlc algorithm.
preliminaries
from now on, we will restrain the study to discrete and finite variables .
for readers that are not familiar with pgms, figure  <dig> clarifies the meaning of specific key terms used hereafter.
bayesian networks are probabilistic graphical models.
they are defined by a directed acyclic graph , g, and a set of parameters, θ. the set of nodes x = {x <dig>  ..., xn} represents n random variables and the set of edges e captures the conditional dependences between these variables .
the variables are either observed or latent.
the set of parameters θ is a matrix of conditional probability distributions θi= where paxi denotes node i's parents.
if a node has no parent, then it is described by an a priori probability distribution.
for further understanding, we now briefly introduce the concepts of marginal independence and conditional independence between two variables.
definition  <dig> the marginal independence between two variables xi and xj is defined referring to the joint distribution p: p = p p.
a non-equality implies that xi and xj are marginally dependent.
definition  <dig> more restrictive, the definition of conditional independence between two variables xi and xj given a subset of variables s ⊆ x\{xi , xj } is the following: p = p p.
a non-equality implies that xi and xj are conditionally dependent given s.
a latent class model  is a particular type of bayesian network.
it is defined as containing a unique latent variable connected to each of the observed variables.
the latent variable simultaneously influences all observed variables and hence renders them dependent.
in the lcm framework, an underlying assumption, called local independence , states that the observed variables are pairwise independent, conditional on the latent variable  <cit> .
the intuition behind li is that the latent variable is the only explanation for the dependences between observed variables.
however, this assumption is often violated for observed data.
to tackle this issue, hlcms were proposed as a generalization of lcms.
hlcms are tree-shaped bns where leaf nodes are observed while internal nodes are not.
in a bayesian network, local dependence between variables may be modeled through the use of an additional latent variable .
on a larger scale, multiple latent variables organized in a hierarchical structure allow high modeling flexibility.
additional file  <dig> illustrates the ability of hlcms to depict a large variety of relations encompassing local to higher-order dependences.
state of the art
hlc model learning
various methods have been conceived to tackle hlcm learning.
these approaches differ by the following points:  structure learning;  determination of the latent variables' cardinalities;  learning of parameters, i.e.
a priori and conditional probabilities;  scalability;  main usage.
as for general bns, besides learning of parameters , i.e.
a priori and conditional probabilities, one of the tasks in hlcm learning is structure  inference.
this task generally remains the most challenging due to the complexity of the search space.
to address this issue, two main categories of hlcm learning methods have been developed.
the first category, structural expectation maximization , successively optimizes θ conditional on s and s conditional on θ. amongst a few proposals, greedy search  <cit>  and dynamic programming  <cit>  were designed.
they explore the space of possible graphs guided by a scoring function, such as the bayesian information criterion   <cit> .
when using maximum likelihood estimation, the bic score prevents model overfitting through a penalty term on the number of parameters in the model.
as regards greedy search, the search space of hlcm structures can be visited through two operations: a structure in the neighborhood of the current structure may either result from the addition or the removal of latent nodes or from the addition or the dismissing of states, for existing nodes.
in the other solution, implemented in gerbil algorithm, dynamic programming discovers the best segmentation of a genomic region into blocks of contiguous snps.
then, for each previously learned block, an lcm is learned.
alternative approaches implement ascending hierarchical clustering , which provides clusters within which the snps are not necessarily contiguous.
in the following, we will use the terms "blocks" and "clusters" to distinguish between these two possibilities.
relying on pairwise dependence strength, wang and co-workers first build a binary tree; then they apply regularization and simplification transformations which may result in subsuming more than two nodes through a latent variable  <cit> .
hwang and collaborators' approach confines the hlcm search space to binary trees augmented with possible connections between siblings   <cit> .
to construct the tree, they design an ahc strategy.
first, a partition of the observed variables into clusters of size  <dig> is performed, based on a mutual information criterion.
any such cluster then defines a new lcm  in the upper layer under construction.
second, the parameters of each lcm are learned.
thus missing values of lvs can be imputed.
therefore these lvs can be considered as observed variables for the next step.
a tree is completed through the iteration of these two steps .
parameter learning requires the determination of the lvs' cardinalities, e.g.
the number of possible states  for each lv.
a simple method is to arbitrarily set a small value for the cardinality.
following this idea, hwang and collaborators constrain lvs to binary variables.
this method is very fast but presents several drawbacks: on the one hand, a too small cardinality can lead to a loss of information in the process subsuming child variables into a unique lv; on the other hand, a too large cardinality can entail model overfitting and heavy computational burden.
wang and co-workers propose a regularization step to reduce the cardinality of an lv y, knowing the cardinality of its neighbor variables zi:|y|=Πi=1k|zi|maxi=1k|zi|.
other authors use a greedy search approach, starting with a preset value and incrementing or decrementing it to meet an optimal criterion.
the latter method has the drawback of entailing computational overload because it runs several steps of the expectation maximization  algorithm, implementing an iterative procedure.
usually, the em algorithm is used for parameter learning in the presence of lvs or missing data, but it is computationally expensive and does not guarantee that the global optimum will be reached.
to speed up the em process, hwang and collaborators implemented a heuristic based on partial imputation of binary lvs' missing values.
thus, the em algorithm is actually run on partially imputed data.
for an lcm containing two child variables yj and yk, the heuristic is the following: all individuals showing the most probable configuration of {yj, yk} are assigned an lv value of  <dig>  similarly, the individuals characterized with the second most probable configuration are assigned an lv value of  <dig>  to avoid getting trapped in local optima while running the em algorithm to learn a set of latent models, other authors adapted a simulated annealing approach  <cit> .
hwang and co-workers' approach is the only one we are aware of that succeeds in processing high-dimensional data: in an application dealing with a microarray dataset, more than  <dig> genes have been processed for around  <dig> samples.
to the best of our knowledge, no running time was reported for this study.
nevertheless, the twofold binarity restriction  and the lack of control for information decay as the level increases are severe drawbacks to reach our aims: i.e.
to achieve realistic snp dependence modeling and perform subsequent association study with sufficient power.
graphical models for ld modeling
to address ld modeling through a probabilistic graphical model framework, various models were proposed: hidden markov models , markov random fields and bayesian networks with or without latent variables.
hmms represent simple but efficient models to partition a snp sequence into blocks, because no structure learning step is required  <cit>  and the latent states may represent common haplotypes.
verzilli and co-workers modeled snp dependences using markov random fields  <cit> .
they designed an mcmc  method to sample over the space of possible graphs while exploiting prior biological knowledge.
their approach allows to discover cliques of dependent snps, to further allow the identification of causal relations between markers and the disease status indicator.
to implement a tractable method for genome-wide data, verzilli and co-workers reduce the space of possible graphs by specifying a maximal physical distance between snps belonging to the same clique, as well as a maximal size of  <dig> snps for any clique.
in the family of bayesian networks without lvs, haploblock implements a statistical model of haplotype block variation  <cit> .
this model's advantage lies in integrating population genetics concepts such as recombination hotspots, bottleneck, genetic drift and mutations.
another method, bntagger, was developed for snp tag selection; it exploits conditional independence between variables  <cit> .
to learn the structure, bntagger implements a greedy search with random restarts; then it determines a subset of independent and highly predictive snps.
the two latter methods were only tested on a small number of snps  and the authors reported running times of  <dig> h for  <dig> snps  <cit>  and between  <dig> and  <dig> h for only  <dig> snps  <cit> .
thus, these methods do not seem fitted to gwas data processing.
regarding the family of bayesian networks with lvs, nefian modeled snp dependences through embedded bayesian networks.
her model is indeed a set of lcms augmented with snp-snp dependences and lv-lv dependences  <cit> .
to learn the model, the snp data sequence is split into contiguous windows of fixed common size.
then, for each window, an lcm is created.
the lack of flexibility of the snp partitioning method used remains a severe draw-back.
zhang and ji also proposed to model ld through a set of lcms, using an sem strategy  <cit> .
their method does not require splitting the sequence into fixed-size windows.
nevertheless, the number of lcms has to be specified.
as far as we know, no execution times were reported for the two latter approaches when run on high-dimensional data.
other methods are based on regularization, such as the graphical lasso  <cit> , and have been applied to learning sparse pgms for proteomics or gene expression studies, whose data dimensionality is high  but lower than that of genome-wide data .
the basic idea is to consider that the observations follow a multivariate gaussian distribution with mean µ and covariance matrix Σ .
if the ijth partial correlation coefficient of the precision matrix Σ- <dig> equals zero, variables i and j are conditionally independent, given the other variables.
the use of lasso aims at restraining the learning task to sparse pgms through finding a least-square solution under the following constraint: ∑ν|βν| ≤ t, meaning that the sum, over the whole variable set, of the absolute values of the regression coefficients v has to be inferior or equal to a constant t. this lasso-based approach has been extended to the case where it is reasonable to assume that the variables can be clustered into groups sharing similar correlation patterns  and where sparse block-structured precision matrices are estimated  <cit> .
to our knowledge, verzilli et al's method is the only one whose tractability regarding gwas data is known.
however, in practice, their mrf modeling reveals a drawback.
the ld is modeled through cliques containing a maximum number of  <dig> snps, whereas, generally, several tens or hundreds of snps may be dependent.
furthermore, no dependences between cliques are taken into account.
in contrast, bns with lvs offer a crucial advantage over other models: they provide synthesizing variables useful to reduce data dimensionality.
however, when the number of variables exceeds several hundreds, implementing the sem approach for ld modeling leads to prohibitive computational burden.
when dealing with genome-wide data, the imperious requirement for tractability leads us to choose a hierarchical clustering approach, in the line of hwang and co-workers.
methods
motivation of the fhlc model for gwass
the hlcms offer several advantages for gwass.
first, beside data dimensionality reduction, they allow a simple test of direct dependence between an observed variable and a target variable such as the phenotype, conditional on the latent variable, parent of the observed variable.
note that the phenotype variable is not included in the hlcm.
in the context of gwass, this test helps find the markers which are directly associated with the phenotype, i.e.
causal markers, should there be any.
second, hlcms can deal with the fuzzy nature of ld blocks.
indeed, hlcms can take into account various degrees of ld strength between any two snps, depending on the height of their lowest common lv node ancestor in the tree.
thirdly, the hierarchical structure allows zooming in through narrower and narrower regions in search for a stronger association with the disease, thus offering a data mining tool.
this zooming process ends pointing out the potential markers of interest.
finally, the latent variables may be interpreted in terms of biological meaning.
for instance, in the case of haplotypes, that is, phased genotypes, the latent variables are likely to represent the so-called haplotype block structure of ld.
to a certain extent, an lv might be interpreted as the shared ancestry of the haplotypes defined by the observed variables, namely, the contemporary haplotypes of the tree rooted in the lv.
each state of an lv may represent a group of similar haplotypes.
in the situation of limited ancestral recombination, similar haplotypes tend to share recent common ancestry.
although this situation is not guaranteed along the genome, it is very likely for low-level lvs, since they are expected to cover very small genomic regions showing strong ld.
thus the directed edges, lv → snp, can represent causal effects and provide a biological sense.
besides, it has to be noted that when the latent variables capture dependences between distant snps , they can be viewed as population structure.
however, snp dependences would better be more wisely modeled through a forest of hlcms.
in the case of a forest, higher-order dependences are captured only when relevant, i.e when meeting a strength criterion.
therefore, fhlcms allow to model a larger set of configurations than hlcms do.
typically, an hlcm is limited to represent clusters of close dependent snps.
actually, in this model, variables are constrained to be dependent upon one another, either directly or indirectly.
consequently, hlcms cannot account for potential independence between groups of distant snps or snps located on different chromosomes.
but realistic modeling requires a more flexible framework.
for instance, the ld plot of the  <dig> mb sequence shown in additional file  <dig> reveals that the greatest part of ld is observed between snps in vicinity.
ld rarely exceeds  <dig> kb between snps.
an fhlcm consists of a directed acyclic graph  also called the structure whose nonconnected components are trees, and of θ, the parameters .
figure  <dig> illustrates a possible structure for an fhlcm.
principle of fhlc model construction
our method can process both genotypic  and haplotypic  data.
it takes as an input a matrix dx defined on a finite discrete domain, say { <dig>   <dig>  2} for unphased snps or { <dig>  1} for phased snps, describing n individuals through p variables x = {x <dig>  ..., xp}.
algorithm cfhlc yields an fhlcm, that is a forest structure and θ, the parameters of a set of a priori distributions and local conditional distributions allowing the definition of the joint probability distribution.
two search spaces are explored: the space of directed forests and the probability space.
in addition, the whole set of latent variables h of the fhlcm is output, together with the associated imputed data matrix.
to handle high-dimensional data, our proposal combines two strategies.
the first strategy splits up the genome-scaled data into contiguous regions.
in our case, splitting into  windows is not a mere implementational trick; it satisfies biological grounds: the overwhelming majority of dependences between genetic markers  is observed for close snps.
the user interested in taking into account long-range ld due to the presence of population structure will be faced with the following choices:  adjusting the window size, relying on biological background defining the maximum physical distance between snps in long-range ld ;  slightly diminishing the density of the studied snp sequence.
a combination of these two approaches may be more convenient.
then, an fhlcm is learnt for each window in turn.
within a window, subsumption is performed through an adapted ahc procedure:  at each agglomerative step, a partitioning method is used to identify clusters of variables;  each such cluster is intended to be subsumed into an lv, through an lcm.
for each lcm, parameter learning and missing data imputation  are performed.
a global schema of our method is presented in figure  <dig>
along with a hierarchy-based proposal of hwang and collaborators  <cit>  developed for gene expression studies, our method also implements data subsumption, meeting the two following additional requirements:  a more flexible thus more faithful modeling of underlying reality,  a control of information decay due to subsumption.
node partitioning
following martin and vanlehn  <cit> , ideally, we would propose to associate a latent variable with any clique of variables in the undirected graph of dependence relations .
in the case when introducing an additional lv increases a scoring function such as the bic score  <cit> , the lcm is validated.
however, searching for such cliques is an np-hard task.
moreover, in contrast with these authors' objective, fhlcms do not allow clusters to have more than one parent each: non-overlapping clusters are required for our purpose.
thus, an approximate method solving a clique partitioning problem when provided with pairwise dependence measures is relevant; the clique partitioning problem consists in finding the best partition of a graph into cliques.
an algorithm meeting this purpose has already been described in the literature: bendor and co-authors designed cast, a clique partitioning algorithm devoted to variable clustering  <cit> .
they especially applied cast for gene expression clustering.
as an input, cast requires a binary similarity matrix.
the adaptation of cast to our case is straightforward: the dependence measure between two snps, evaluated through mutual information, is used as a similarity measure.
all mutual information values less than a threshold tmi are assigned a similarity value of  <dig>  whereas the others are assigned a value of  <dig>  as a threshold tmi , the median value  of the mutual information matrix can be used.
then, the cast algorithm constructs the clusters one at a time.
the authors define the affinity a of an element x to be the sum of similarity values between x and the elements present in the current cluster copen.
x is an element of high affinity if it verifies inequality a ≥ tcast|copen|, where tcast is a specified similarity threshold.
otherwise, x is considered an element of low affinity.
to summarize, the algorithm alternates between adding high affinity elements to copen and removing low affinity elements from it.
when the process stabilizes, copen is closed.
a new cluster can be started.
determining cardinalities for latent variables
a steep task is choosing - ideally optimizing - the cardinality of each lcm's latent variable.
this problem cannot be remedied using greedy search because of its intractability regarding high data dimensionality.
although the regularization method of wang and collaborators has the advantage of being very quick , in our context, their method is impracticable.
for instance, let us consider an hlcm learned from genome-wide data.
the first layer of the hlcm contains the majority of the lvs in the model.
in our case, an lv in the first layer can subsume more than  <dig> child ovs .
as the cardinality is the same for all ovs , the resulting cardinality of the lv after regularization remains generally very large.
for example, for an lcm containing  <dig> ovs {x <dig> , ..., x <dig> } of cardinalities equal to  <dig>  the cardinality of the lv h would be : |h| = 310/ <dig> =  <dig> =  <dig>  the simplest solution remains to arbitrarily set a small value for lv cardinalities, but it has several drawbacks .
instead of using an arbitrary constant value common to all latent variables, we propose that the cardinality be estimated for each latent variable through a function of the underlying cluster's size.
the rationale for choosing this function is the following: the more child nodes a latent variable has, the larger the total number of possible combinations is for the values of the child variables and the larger also is the expected number of such combinations observed over all individuals .
therefore, the cardinality of this latent variable should depend on the number of child nodes.
nonetheless, to keep the model complexity within reasonable limits, a maximum cardinality is fixed.
parameter learning and imputation
parameter learning is carried out step by step, each time generating additional latent variables and imputing their values for each individual.
at ith step, this task simply amounts to performing parameter learning for as many lc models as there are clusters of variables identified.
we recall that the nodes in the topology of an lcm are reduced to a unique root and leaves.
therefore, at ith step, each lcm's structure is rooted in a newly created latent variable.
when latent variables are the source nodes in a bn, parameter learning may be performed through a standard em procedure.
this procedure takes as an input the cardinalities of the latent variables and yields the probability distributions, that is, prior distributions for those nodes with no parents and distributions conditional to parents for the remaining nodes.
after imputing the missing data corresponding to latent variables, new data are available to seed the next step of the fhlcm construction: latent variables identified through step i will be considered as observed variables during step i +  <dig>
it has to be noted that designing an imputation method to infer the values of the latent variable for each individual is a matter for investigation.
once the prior and conditional distributions have been estimated for a given lcm, probabilistic inference in bns may be performed.
a straightforward way would consist in imputing the latent variable value for each individual as follows: h*=argmaxh{p}.
however, in the framework of probabilistic models, this deterministic approach is disputable.
in contrast, a more convincing alternative will draw a value h for latent variable h, knowing the probabilities p for each individual.
controlling information decay
conversely to hwang and co-workers' approach, which mainly aims at data compression, information decay control is required: in step i, any candidate latent variable h which does not bear sufficient information about its child nodes must be invalidated.
as a consequence, such child nodes will be seen as isolated nodes in step i +  <dig>
let us consider two variables x and h. basically, the mutual information measures the difference of entropies between the independent model p p and the dependent model p p: ℐ=+ℋ)-+ℋ)=ℋ-.
therefore, the mutual information measures the dependence of the two variables.
the larger the difference between entropies, the higher is the dependence.
now, let us consider a set of child variables x= {x <dig>  x <dig>  ..., xn} and the parent variable h. in our case, we want to compare the two models: p p ... p p and p p ... p p. thus, Δ, the difference of entropies between the two models is: +ℋ)-+ℋ)=∑i=1n-ℋ)=∑i=1nℐ.
Δ corresponds to the sum of mutual information values over all lcm's edges.
normalization through entropy and averaging are performed to provide a more intuitive criterion: c=1sh∑i ∈ clusterℐmin,ℋ), with sh the size of cluster .
c represents the average percentage of information captured by the lv with respect to its child variables.
algorithm
the sketch of cfhlc is presented in algorithms  <dig> and  <dig>  the user may tune seven parameters.
window size s specifies the number of contiguous snps - or variables - spanned per window.
the aforementioned criterion c is meant to estimate information decay, thus allowing information dilution to be constrained to a minimal threshold t. parameters a, b and cardmax participate in the calculus of the cardinality of each latent variable.
finally, parameter partitioningalg enables flexibility in the choice of the method dedicated to clustering highly-correlated variables into non-overlapping groups.
within each window i, the ahc process is initiated from the first layer consisting of univariate models.
each such univariate model is built for any observed variable in the set wi .
the ahc process stops if each cluster identified is reduced to a singleton  or if no cluster of size strictly greater than  <dig> could be validated .
each cluster containing at least two nodes is subject to lcm learning  followed by validation .
in order to simplify the fhlcm learning, the cardinality of the latent variable is estimated as an affine function of the number of variables in the corresponding cluster .
algorithm learn_latent_class_model is plugged into this generic framework .
after validation through threshold t , the lcm is used to enrich the fhlcm associated with the current window : a specific merging process links the additional node corresponding to the latent variable to its child nodes, themselves already present in the fhlcm structure under construction; the prior distributions of the child nodes are replaced with distributions conditional on the latent variable.
the newly created latent variable, ljk, is added to the set of latent variables, whereas its imputed values, d, are stored .
in wi, the variables in cℓjk are now replaced as a whole with the corresponding latent variable; data matrix d is updated accordingly .
in contrast, the nodes in unvalidated clusters are kept isolated for the next step.
finally, the collection of forests, dag, is successively augmented with each forest built within a window .
in parallel, due to assumed independence between windows, the joint distribution of the final fhlcm is merely computed as the product of the distributions associated with the windows .
input:
x, a set of p variables ,
dx, the corresponding observations for n individuals,
s , a window size,
c, a criterion designed to estimate information decay while building the fhlcm,
t, a threshold used to constrain information dilution, based on criterion c,
partitioningalg, an algorithm dedicated to partition a set of variables into non-overlapping clusters of variables,
a; b and cardmax, parameters used to estimate the cardinality of latent variables.
output:
dag and θ, respectively the dag structure and the parameters of the fhlcm constructed,
l, the whole set of latent variables identified through the construction , dl, the corresponding data imputed for n individuals.
1: nbw← p/s /* computation of the number of contiguous windows */
2: dag ←∅; θ ← ∅; l← ∅; dl← ∅
3:
4: for i =  <dig> to nbw
5:    /* processing of layer  <dig> */
6:    wi←{x×s+ <dig> ...,xi×s};d←d
7:    {∪j∈widagunivj,∪j∈wiθunivj}←learn_univariate_models
8:    dagi←∪j∈widagunivj;θi←∪j∈wiθunivj
9:
10:    step ← 1
11:    while true
12:    {cℓ <dig>  ..., cℓ#c}←partition
13:    if all clusters cℓq are singletons then break end if
14:
15:    {cℓj <dig>  ...,cℓj#c2}←identify_clusters_of_size_strictly_greater_than_one
16:    nbvalidclusters ← 0
17:
18:    for k =  <dig> to #c2
19:    card lv ← min + b; cardmax)
20:    {dagjk, θjk, ljk, {dagjk,θjk,ljk,d}←learn_latent_class_model
21:
22:    /* validation of current cluster - see subsection controlling information decay */
23:    if ≥t)
24:       incr
25:       dagi ← merge_structures; θi ← merge_parameters
26:       l←l∪ljk;dl←dl∪d
27:       d←∪d;wi←∪ljk
28:    end if
29: end for
30:
31: if  then break end if
32:
33: incr
34: end while
35:
36: dag ← dag ∪ dagi; θ ← θ × θi
37: end for
algorithm 1: cfhlc
input:
cℓu: a cluster containing at least two nodes,
d: the corresponding observations for n individuals,
cardlv: the cardinality of the latent variable to be created.
output:
a latent class model described by:
dagu and θu, respectively the structure and the parameters of the latent class model,
lu, a latent variable,
d, the data imputed for the latent variable .
1: lu ← create_latent_variable()
2: dagu ← build_structure_of_latent_class_model
3: θu ← run_standard_em
4: d ← impute_data
algorithm 2: learn latent class model
experimental
RESULTS
implementation
algorithm cfhlc has been developed in c++, relying on the probt library dedicated to bns http://bayesian-programming.org.
we have plugged into cfhlc a c++ implementation of cast based on the original implementation provided in java by ben fry http://benfry.com/clustering/.
regarding the visualization of the dags, the software tulip http://tulip.labri.fr/tulipdrupal/ was chosen, meeting both high representation quality and compactness requirements.
cfhlc was run on a standard pc .
experimental protocol
the performance of the fhlcm-based method is evaluated using real phased and unphased genetic data, on the one hand, and simulated phased and unphased genetic data on the other hand.
for real data analysis, the well-known daly et al.
dataset  <cit> , available at http://www-genome.wi.mit.edu/humgen/ibd5/index.html, was used.
this dataset consists of  <dig> trios, each composed of two parents and one child.
for each individual,  <dig> snps are genotyped in the 5q <dig> region and cover  <dig> kb.
we only analyzed the child data.
regarding simulated data, two well known programs were used: hapgen and hap-simu.
using hapgen http://www.stats.ox.ac.uk/~marchini/software/gwas/hapgen.html, we generated  <dig> unrelated individuals  for a several hundreds kb region containing around 20- <dig> snps.
the haplotypes used as references come from the hapmap phase ii and concern u.s. residents of northern and western european ancestry  http://hapmap.ncbi.nlm.nih.gov/.
five sequences showing variable ld degrees ) ranging from  <dig>  to  <dig>  were generated.
hapsimu http://l.web.umkc.edu/liujian/ was used to simulate genotypes with simulation parameters described in additional file  <dig>  three sample sizes were chosen with respect to the number of observed variables:  <dig> k,  <dig> k and  <dig> k snps .
for each sample size, twenty benchmarks were generated.
for these experimentations, imputation of lvs' values was achieved by assigning the most probable values given the observations.
a drawback of this method is the loss of probabilistic relation between a variable and its parent variable.
a definite advantage lies in its running in around half the time required by imputation through simulation .
ld modeling
real data
regarding the daly benchmark, our aim was to evaluate how the forest obtained keeps up with the real structure of the biological data.
moreover, the cfhlc algorithm was compared to four other methods.
we learned fhlcms on both haplotype  and genotype  data.
the corresponding graphs are displayed in figures  <dig> and  <dig>  respectively.
globally, the two graphs are similar: most of snps which are connected through an lv in the haplotype-data graph  are also connected through an lv in the genotype-data graph , e.g.
snp <dig>  snp <dig> and snp <dig>  moreover, a substantial part of these snps share a common parent in both graphs: for instance, in both hdg and gdg, we observe that snp <dig> and snp <dig> are linked by an lv belonging to layer  <dig>  thus, learning fhlcm from genotype data instead of haplotype data leads to similar hierarchical structures.
however, on average, we observe that the snps in the gdg are more connected:  <dig> and  <dig> connected components are identified in the hdg and the gdg, respectively.
for example, the two framed trees  <dig> and  <dig> in the hdg of figure  <dig> are linked by a high-level lv in the gdg of figure  <dig> .
we expect that the fhlcms' graphs will reflect the "haplotype block structure": large blocks of correlated contiguous snps separated by recombination hotspots.
first, we observe that the physical position of snps influences their connection, since close snps tend to be linked by an lv belonging to a low layer, whereas distant snps are generally connected by a high-level lv.
however, strong dependences between distant snps are also observed, e.g.
between snp <dig> and snp <dig> or snp <dig> and snp <dig> .
this characteristic reveals that the ld structure is not only dominated by spatial effects and justifies our haplotype cluster approach .
in addition, the graphs interestingly show trends consistent with biological reality, that is the variation of the recombination rates inferred by software phase v <dig>   <cit>  along the studied sequence .
indeed, most of subtrees rooted in low-level lvs cover regions with low recombination rates .
more than 68% and 94% of lvs from layer  <dig> cover chromosomic segments showing rrs below  <dig> cm/mb and  <dig> cm/mb, respectively.
the same tendency is observed for more than 44% and 66% of lvs from layer  <dig>  respectively.
these results show the relevance of  interpreting low-level lvs as haplotype shared ancestry when cfhlc's input is haplotype data.
we compared the structure obtained by cfhlc with those output by four other approaches: daly et al.
's method  <cit> , gerbil  <cit> , haploblock  <cit>  and zhang et al.
's algorithm  <cit> .
all these methods were detailed in section state of the art.
the three former methods partition the sequence into blocks of contiguous snps.
in contrast, the latter algorithm yields  clusters of non-contiguous snps.
we recall that cfhlc algorithm generates a hierarchical clustering of non-contiguous snps.
in figure  <dig>  we compare the haplotype block- or cluster-structures obtained through all five methods aforecited.
in spite of the fact that these methods differently tackle ld modeling, common trends emerge .
for instance, the last block identified by daly et al.
's method, gerbil and zhang et al.
's algorithm  is also inferred by our algorithm in line  <dig>  slight differences are observed with the two first blocks resulting from daly et al.
's method and gerbil which only form one block for zhang et al.
's algorithm  and cfhlc .
compared to other methods, most divergences with our algorithm remain in its unique ability to take into account the fuzzy delimitations of clusters.
this is illustrated with the central area of the sequence , which actually presents two weak recombination hotspots .
another difference with the other methods is the presence of "unclustered" snps, like snp <dig>  snp <dig> and snp <dig> in our model.
the running times, dimension reduction rates and entropy compression rates of all methods are reported in table  <dig>  results show that gerbil is the fastest algorithm tested, with a running time of  <dig> s. however, cfhlc and zhang et al.
's algorithm, which learn more complex models , achieve their tasks in quite a reasonable time,  <dig> s and  <dig> s, respectively.
compared to others, haploblock is the slowest method, with a running time of  <dig> mn, due to the high complexity of learning models based on population genetics.
for the three methods exhibiting a partition of contiguous snps, we defined the dimension reduction rate  as the ratio of the number of blocks to the number of snps.
as regards zhang et al.
's algorithm, the drr was defined as the ratio of the number of clusters to the number of snps.
in the case of cfhlc, we consider that the information of each fhlcm's tree can be synthesized by its root, providing the best dimension reduction.
therefore, in this case, the drr is defined as the number of roots in the whole forest divided by the number of snps.
haploblock generates the lowest number of blocks with an average of  <dig>  , whereas zhang et al.
's algorithm partitions the sequence in  <dig> clusters , and daly et al.
's method and gerbil both identify  <dig> blocks .
cfhlc presents the lowest dimension reduction with  <dig> trees , due to the presence of  <dig> "unclustered" snps: snp <dig>  snp <dig>  snp <dig>  snp <dig>  snp <dig>  snp <dig> and snp <dig>  as an alternative measure of compression, we defined the entropy compression rate  as the ratio of the sum of block  entropies in a partition to the entropy, assuming no structure .
we observe a different ranking of the methods.
we notice that cfhlc and zhang et al.
's algorithm, which both learn cluster models, provide the best  ecr values , whereas haploblock, gerbil and daly et al.
's method show ecr values of  <dig> ,  <dig>  and  <dig> , respectively.
regarding the ecr criterion, the comparatively better results obtained for the two cluster models are explained by the absence of the constraint for compulsory physical proximity between snps .
moreover, the ecr criterion does not penalize anymore the cfhlc algorithm, since the unclustered snps contribute relatively little to the overall information content.
we ran the last three programs on a standard computer.
as we had no access to daly et al.
's software, we could only compare the dimension reduction rates and entropy compression rates calculated from their results with the dimension reduction rates and entropy compression rates obtained with the other methods.
in subsection motivation of the fhlc model for gwass, we argued that the multiple layers of an fhlcm can describe various degrees of ld strength.
to analyze this property, we plotted the r <dig> squared correlation coefficient of any pair of snps against the level of their most recent common ancestor .
figure  <dig> and  <dig> show such plots drawn for haplotype and genotype data, respectively.
starting from values in the range , the r <dig> correlation coefficient quasi-linearly decreases when the mrca level increases.
we conclude that the layered structure of the fhlcm faithfully reflects ld strength variety.
these encouraging results lead us to visually compare the ld plot and the triangular matrices of the mrca levels for haplotype  and genotype  data, as presented in figure  <dig>  for this purpose, the same color code was used in the ld plot and the triangular matrix of mrca levels.
in the ld plot, the color  of each cell varies with the r <dig> value.
since the median value of r <dig> can be computed for each mrca level, the color of each cell in the triangular matrix of mrca levels is set, relying on the color scale used for the ld plot.
we expected to observe a correspondence between the three plots.
the outstandingly clear correspondence demonstrates the ability of fhlcms to accurately model multiple levels of ld strength: the overall majority of the ld plot dependences are also present in the mrca level matrix.
interestingly, modeling from genotype data leads to quite good results compared to haplotype data.
haplotype diversity is generally very low within haplotype blocks or clusters.
in our hierarchical model, haplotype diversity is expected to be all the larger within a cluster as the level of the lv subsuming this cluster is high.
to check this point, we have relied on the clusters of figure  <dig>  for each lv, haplotype diversity has been calculated as the number of the most common haplotypes observed at level l of the tree rooted in this lv.
figure  <dig> plots the number of the most common haplotypes against the lv level.
the plot shows that the haplotype diversity median remains very low  for the first four layers and dramatically increases to around  <dig> in the fifth layer.
these results confirm our expectation relative to haplotype diversity in fhlcms.
simulated data
finally, the impact of varying ld degrees was studied.
for this purpose, we generated haplotype data with the hapgen software.
five sequences, showing variable ld degrees ) ranging from  <dig>  to  <dig> , were used to learn fhlcms.
figure  <dig> shows the forests obtained.
the forests reveal that increasing ld degrees entails higher graph connectivity as well as a larger number of layers.
indeed, when median equals  <dig> ,  <dig> connected components are identified and the highest lv belongs to the third layer.
conversely, in the case when median is equal to  <dig> , the forest is only composed of  <dig> connected components and the highest lv belongs to layer  <dig>  thus, we conclude that cfhlc can process sequences with various ld degrees and generate fhlcms whose structures reliably reflect linkage disequilibrium.
scalability for gwass
scalability has been studied through the data simulated with hapsimu.
in the hardest case , additional file  <dig> shows that only  <dig> hours are required with a window size s set to  <dig>  for the same dataset processed in the cases "s = 200" and "s = 600", running times are  <dig>  h and  <dig>  h, respectively.
for the same number of ovs , wang et al.
report running times of about two months.
regarding the  <dig> k case, running times are  <dig>  h,  <dig> h and  <dig>  h for "s = 100", "s = 200" and "s = 600", respectively.
further analysis of cfhlc algorithm
finally, many other experimentations are reported with their commentaries in additional files  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and  <dig>  additional file  <dig> focuses on the impact of window size on running time.
in additional file  <dig> data dimension reduction is evaluated from the distribution of lvs over the forest's layers.
additional files  <dig>   <dig>   <dig> and  <dig> study the impact of window size on the number of roots in the forest, the number of lvs, the number of layers and the distribution of lvs over the forest's layers, respectively.
additional files  <dig> and  <dig> analyze how information fades while the layer number increases.
an important result is that cfhlc can achieve a data dimensionality reduction of more than 80% of the number of observed variables .
regarding spatial complexity of cfhlc algorithm, the entire fhlcm does not require to be stored in ram because each window can be saved on the hard disk.
CONCLUSIONS
our contribution in this paper is twofold:  a new framework has been described, which was tested on both real and simulated data and was proven able to concisely model ld and to reduce snp data dimensionality;  cfhlc, an algorithm dedicated to learn fhlcms, has been shown to be efficient when run on genome-scaled benchmarks.
compared to verzilli and co-workers' works, our algorithm provides a more accurate modeling of ld and synthesizes genetic marker information through lvs.
in addition, unlike nefian or zhang and ji, our method does not require to specify the number of lcms and can capture multiple levels of dependences, thus taking into account the fuzzy nature of ld.
to our knowledge, our hierarchical method is the first one shown to achieve fast model learning for genome-scaled data sets, while maintaining satisfying information scores and relaxing the twofold binarity restriction of hwang and collaborators' model .
hwang and collaborators' purpose is only data compression.
we are faced with a more demanding challenge: to make a sufficiently powerful down-stream association analysis possible.
in discussing the biological interpretation of latent variables, we mentioned the potentiality of fhlcms for population substructure description.
in essence, using hierarchical models is highly appealing to take into account the long-range ld expected in substructured populations.
however, this interesting use of such hierarchical models as fhlcms is somewhat precluded by the technical necessity to partition the genome into small regions.
as a first palliative, we indicated two strategies  to cope with this current technical limitation.
however, the strong expectation for faithful substructure modeling through fhlcms advocates further efforts to clear the hurdle on path to realistic long-range ld modeling.
a bottleneck currently lies in the clique partitioning method chosen, which forbids window sizes encompassing more than  <dig> observed variables.
in addition to investigating alternative partitioning methods, a lead to cope with this bottleneck may be to adapt the specific processings at the limits of contiguous windows or use overlapping windows.
in short, fhlcms can be used to resolve several major problems in the gwass' context.
beside flexible data dimension reduction through fhlcms' lvs, fine mapping of causal snps is expected thanks to conditional independence properties encoded in such models.
for instance, fhlcms' lvs can be used to condition tests for independence between a snp and the phenotype.
moreover, due to their hierarchical structure, fhlcms represent an original and appropriate solution to study long-range ld in substructured populations, a recurring problematic in gwass.
finally, genome-wide visualization of ld can be easily achieved with these models using a graph visualization tool and will provide an intuitive representation of snp - snp dependences as well as information synthesis through latent variables.
in this current version of cfhlc, when processing haplotype data, we were not interested in knowing the sequence of each ancestral haplotype.
we just wanted to know from which ancestral haplotype  a contemporary haplotype comes.
nevertheless, it is feasible to infer the sequence of each ancestral haplotype using probabilistic inference.
finally, although our modeling is designed for gwas data, we emphasize that it could be applied to other data presenting spatial dependences between variables, in particular sequential data.
beyond this specific case, it would be interesting to assess the model's generality, in order to determine if it can be applied to generic graphical model learning.
authors' contributions
rm and cs both wrote the manuscript.
rm proposed the modeling of the genetic data, the conception of a tractable algorithm and carried out the implementation and the experiments.
cs designed the study, participated in its coordination and partook in the conception of the algorithm.
pl participated in the design of the study, the conception of the algorithm, helped to improve the manuscript and provided indispensable feedback.
all authors read and approved the final version of the manuscript.
supplementary material
additional file 1
direct and indirect associations between a genetic marker and the phenotype.
the figure included into this additional file illustrates the cases of direct and indirect associations between a genetic marker and the phenotype.
click here for file
additional file 2
linkage disequilibrium plot for a simplified haplotype block structure.
the figure included into this additional file describes a standard representation of pairwise dependences between genetic markers.
click here for file
additional file 3
linkage disequilibrium plot of a real  <dig> kb snp sequence.
the figure presented in this additional file shows the linkage disequilibrium plot of a real  <dig> kb snp sequence.
click here for file
additional file 4
hierarchical latent class model.
the figure presented in this additional file depicts a hierarchical latent class model.
click here for file
additional file 5
linkage disequilibrium plot of a  <dig> mb snp sequence.
the figure included in this additional file describes the linkage disequilibrium plot of a  <dig> mb snp sequence.
click here for file
additional file 6
parameter value adjustment for the generation of simulated genotypic data through software hapsimu.
the table included in this additional file enumerates the values chosen for the parameters of software hapsimu.
click here for file
additional file 7
average running time versus number of variables.
the figure presented in this additional file plots the running time of the cfhlc algorithm versus the number of snps in the dataset.
click here for file
additional file 8
impact of window size on running time.
the figure presented in this additional file plots the running time of the cfhlc algorithm versus the window size.
click here for file
additional file 9
number of variables per layer over the whole fhlc model.
the figure included in this additional file describes the average distribution of the variables over the layers .
click here for file
additional file 10
impact of window size on the number of roots.
the figure included in this additional file depicts the impact of window size on the number of roots.
click here for file
additional file 11
impact of window size on the number of latent variables.
the figure presented in this additional file shows the impact of window size on the number of latent variables.
click here for file
additional file 12
impact of window size on the number of layers.
the figure presented in this additional file describes the impact of window size on the number of layers.
click here for file
additional file 13
impact of window size of the number of latent variables per layer and on the ratio of the number of latent variables per layer to the total number of variables.
the two subfigures included in this additional file depict the impact of window size on the number of latent variables per layer on the one hand and the impact of window size on the number of latent variables per layer to the total number of variables, on the other hand.
click here for file
additional file 14
impact of window size on scaled mutual information, per layer.
the figure presented in this additional file describes the impact of window size on scaled mutual information, per layer, over the whole fhlc model.
click here for file
additional file 15
average scaled mutual information per layer over the whole fhlc model; impact of parameters a and b. the figure presented in this additional file shows the impact of parameters a and b on scaled mutual information, per layer, over the whole fhlc model.
click here for file
acknowledgements
the authors are grateful to two anonymous referees for constructive comments and help in improving their manuscript.
the authors would also like to thank christian dina and jean-jacques schott  for useful discussions regarding genetics and statistical genetics.
yulong zhang  has to be thanked for his precious help in the comparison between his algorithm and cfhlc algorithm.
rm's work was supported by a grant from the bil bioinformatics research project of pays de la loire region, france.

 genetic markers such as snps are the key to dissecting the genetic susceptibility of common complex diseases, such as asthma, diabetes, atherosclerosis and some cancers. the purpose is identifying combinations of genetic determinants which should accumulate among affected subjects. generally, in such combinations, each genetic variant only exerts a modest impact on the observed phenotype, whereas, in contrast, the interaction between genetic variants and, possibly, environmental factors is determinant. decreasing genotyping costs now enable the generation of hundreds of thousands of snps, spanning the whole human genome, across cohorts of cases and controls. this scaling up to genome-wide association studies makes the analysis of high-dimensional data a hot topic. despite recent technological advances and extensive research effort, the genetic basis of the aforementioned diseases remains to a large extent unknown. yet, the search for associations between single snps and the variable describing case/control status requires carrying out a large number of statistical tests. since snp patterns, rather than single snps, are likely to be determinant for complex diseases, a high rate of false positives as well as a perceptible statistical power decrease, not to mention intractability, are severe issues to be overcome. the simplest type of genetic polymorphism, single nucleotide polymorphism, involves only one nucleotide change, which occurred generations ago within the dna sequence. to fix ideas, we emphasize that one single individual can be uniquely defined by only to independent snps and unrelated individuals differ in about % of their billion nucleotides. compared with other kinds of dna markers, snps are appealing because they are abundant, genetically stable and amenable to high-throughput automated analysis. consistently, advances in high-throughput snp genotyping technologies lead the way to various down-stream analyses, including gwass. exploiting the existence of statistical dependences between neighboring snps, also called linkage disequilibrium, is the key to association study achievement. indeed, a causal variant may not be a snp. for instance, insertions, deletions, inversions and copy-number polymorphisms may be causative of disease susceptibility. nevertheless, a well-designed study will have a good chance of including one or more snps that are in strong ld with a common causal variant. in the latter case, indirect association with the phenotype, say affected/unaffected status, will be revealed. interestingly, ld also offers solutions to reduce data dimensionality in gwass. in the human genome, ld is highly structured into the so-called "haplotype block structure": regions where statistical dependences between contiguous markers are high alternate with shorter regions characterized by low statistical dependences. the most likely explanation of this phenomenon is related to the presence of large regions with low recombination rates separated by recombination hotspots. relying on this feature, various approaches were proposed to achieve data dimensionality reduction: testing association with haplotypes, partitioning the genome according to spatial correlation, selecting snps informative about their context, or snp tags. recent methods, such as haplobuild, have permitted to construct more biologically relevant haplotypes where the "haplotype cluster structure", instead of the "haplotype block structure", is assumed: haplotypes are not constrained by contiguous orientation. unfortunately, these methods do not take into account all existing dependences since they miss higher-order dependences. actually, these methods do not consider the fuzzy nature of ld: the ld block boundaries are not accurately defined over the genome. due to their ability to represent conditional independences between variables, probabilistic graphical models offer an adapted framework for an accurate modeling of dependences between snps. a pgm is a probabilistic model relying on a graph representing conditional independences within a set of random variables. inherently, this model simplifies the description of the joint distribution of the set of variables. several subclasses of pgms exist such as markov random fields and bayesian networks. the main difference between these two subclasses remains in the nature of the graph: in contrast with mrfs, bayesian networks are directed graphs. although the observed variables are often sufficient to describe their joint distribution, sometimes, additional unobserved variables, also named latent variables, have a role to play. only few research works have been dedicated to snp dependence modeling through pgms. a hard task because of high data dimensionality, tackling this modeling issue through pgms nevertheless offers an attractive lead. approaches based both on mrfs and bns have been designed. regarding the latter, some methods only consider observed variables whereas other models include latent variables. in particular, hierarchical bns are the most promising models for ld representation: their hierarchical structure supported by lvs allows flexible information synthesis, thus efficiently reducing the data dimensionality. to our knowledge, modeling ld through hierarchical bns in order to reduce snp data dimensionality has not yet been designed. notably, scalability remains a crucial issue for gwass. in this paper, we emphasize the interest of using a forest of hierarchical latent class models, to reduce the dimension of the data to be further submitted to statistical analyses devoted to the discovery of genetic factors potentially involved in the disease. such studies encompass single-snp analysis, multiple-snp analysis, snp-snp interaction analysis and analysis integrating gene expression. an fhlcm is a hierarchical bn with discrete observed and latent variables. basically, latent variables capture the information borne by underlying markers. in their turn, latent variables are clustered into groups and, if relevant, such groups are subsequently subsumed by additional latent variables. iterating this process yields a hierarchical structure. first, the great advantage to gwass is that further statistical analyses can be chiefly performed on latent variables. thus, a reduced number of variables will be examined. second, a model based on a hierarchical structure provides a flexible data mining tool. for example, different degrees of data dimensionality reduction are available to the statistician. moreover, the hierarchical structure is meant to efficiently conduct refined association testing: zooming in through narrower and narrower regions in search for a stronger association with the disease ends pointing out the potential markers of interest. however, most algorithms dedicated to the learning of hierarchical latent class models fail the scalability criterion when the data describe thousands of variables and a few hundreds of individuals. in a previous work-on progress paper, we designed an algorithm devoted to learning fhlcms. this algorithm was named cfhlc, which stands for construction of forests of hierarchical latent class models. the contribution brought in the present extended version is the following: we advocate the use of fhlcms to model ld; we provide a detailed description of the main concepts underlying our approach; using real data, we show that the fhlcm graph is representative of the haplotype cluster structure; in addition, we compare the haplotype cluster structure obtained through cfhlc with those output by four other algorithms; relying on both real and simulated data, we demonstrate the ability of fhlcms to concisely model snp dependences, showing that the multiple layers of the model can take into account different ld degrees and haplotype diversity; finally, we present a thorough study focused on both scalability and impact of adjustment of the input parameters of cfhlc algorithm. as a prerequisite to further understanding, section preliminaries provides an informal definition of bayesian networks, focusing on latent class models and hierarchical latent class models. then the section dedicated to the state of the art first points out the few anterior works devoted to hlcm learning in general. this section ends with a short review of the few attempts to implement probabilistic graphical models for the specific purpose of ld modeling. section methods motivates the modeling of ld through fhlcms and informally describes such models. then, the focus is set on the general outline of the method proposed for fhlcm learning. the next section depicts the sketch of algorithm cfhlc. the last section is dedicated to experimental results and discussion. in this section, we first test and discuss the ability of fhlcms to accurately represent the haplotype cluster structure of genetic data. then, we compare our algorithm to other methods with respect to faithfulness in ld modeling and data dimension reduction. we end the section with a thorough study centered on scalability and influence of the input parameters of the cfhlc algorithm. preliminaries from now on, we will restrain the study to discrete and finite variables. for readers that are not familiar with pgms, figure clarifies the meaning of specific key terms used hereafter. bayesian networks are probabilistic graphical models. they are defined by a directed acyclic graph, g, and a set of parameters, θ. the set of nodes x = {x..., xn} represents n random variables and the set of edges e captures the conditional dependences between these variables. the variables are either observed or latent. the set of parameters θ is a matrix of conditional probability distributions θi= where paxi denotes node i's parents. if a node has no parent, then it is described by an a priori probability distribution. for further understanding, we now briefly introduce the concepts of marginal independence and conditional independence between two variables. definition the marginal independence between two variables xi and xj is defined referring to the joint distribution p: p = p p. a non-equality implies that xi and xj are marginally dependent. definition more restrictive, the definition of conditional independence between two variables xi and xj given a subset of variables s ⊆ x\{xi, xj } is the following: p = p p. a non-equality implies that xi and xj are conditionally dependent given s. a latent class model is a particular type of bayesian network. it is defined as containing a unique latent variable connected to each of the observed variables. the latent variable simultaneously influences all observed variables and hence renders them dependent. in the lcm framework, an underlying assumption, called local independence, states that the observed variables are pairwise independent, conditional on the latent variable. the intuition behind li is that the latent variable is the only explanation for the dependences between observed variables. however, this assumption is often violated for observed data. to tackle this issue, hlcms were proposed as a generalization of lcms. hlcms are tree-shaped bns where leaf nodes are observed while internal nodes are not. in a bayesian network, local dependence between variables may be modeled through the use of an additional latent variable. on a larger scale, multiple latent variables organized in a hierarchical structure allow high modeling flexibility. additional file illustrates the ability of hlcms to depict a large variety of relations encompassing local to higher-order dependences. state of the art hlc model learning various methods have been conceived to tackle hlcm learning. these approaches differ by the following points: structure learning; determination of the latent variables' cardinalities; learning of parameters, i.e. a priori and conditional probabilities; scalability; main usage. as for general bns, besides learning of parameters, i.e. a priori and conditional probabilities, one of the tasks in hlcm learning is structure inference. this task generally remains the most challenging due to the complexity of the search space. to address this issue, two main categories of hlcm learning methods have been developed. the first category, structural expectation maximization, successively optimizes θ conditional on s and s conditional on θ. amongst a few proposals, greedy search and dynamic programming were designed. they explore the space of possible graphs guided by a scoring function, such as the bayesian information criterion. when using maximum likelihood estimation, the bic score prevents model overfitting through a penalty term on the number of parameters in the model. as regards greedy search, the search space of hlcm structures can be visited through two operations: a structure in the neighborhood of the current structure may either result from the addition or the removal of latent nodes or from the addition or the dismissing of states, for existing nodes. in the other solution, implemented in gerbil algorithm, dynamic programming discovers the best segmentation of a genomic region into blocks of contiguous snps. then, for each previously learned block, an lcm is learned. alternative approaches implement ascending hierarchical clustering, which provides clusters within which the snps are not necessarily contiguous. in the following, we will use the terms "blocks" and "clusters" to distinguish between these two possibilities. relying on pairwise dependence strength, wang and co-workers first build a binary tree; then they apply regularization and simplification transformations which may result in subsuming more than two nodes through a latent variable. hwang and collaborators' approach confines the hlcm search space to binary trees augmented with possible connections between siblings. to construct the tree, they design an ahc strategy. first, a partition of the observed variables into clusters of size is performed, based on a mutual information criterion. any such cluster then defines a new lcm in the upper layer under construction. second, the parameters of each lcm are learned. thus missing values of lvs can be imputed. therefore these lvs can be considered as observed variables for the next step. a tree is completed through the iteration of these two steps. parameter learning requires the determination of the lvs' cardinalities, e.g. the number of possible states for each lv. a simple method is to arbitrarily set a small value for the cardinality. following this idea, hwang and collaborators constrain lvs to binary variables. this method is very fast but presents several drawbacks: on the one hand, a too small cardinality can lead to a loss of information in the process subsuming child variables into a unique lv; on the other hand, a too large cardinality can entail model overfitting and heavy computational burden. wang and co-workers propose a regularization step to reduce the cardinality of an lv y, knowing the cardinality of its neighbor variables zi:|y|=Πi=1k|zi|maxi=1k|zi|. other authors use a greedy search approach, starting with a preset value and incrementing or decrementing it to meet an optimal criterion. the latter method has the drawback of entailing computational overload because it runs several steps of the expectation maximization algorithm, implementing an iterative procedure. usually, the em algorithm is used for parameter learning in the presence of lvs or missing data, but it is computationally expensive and does not guarantee that the global optimum will be reached. to speed up the em process, hwang and collaborators implemented a heuristic based on partial imputation of binary lvs' missing values. thus, the em algorithm is actually run on partially imputed data. for an lcm containing two child variables yj and yk, the heuristic is the following: all individuals showing the most probable configuration of {yj, yk} are assigned an lv value of similarly, the individuals characterized with the second most probable configuration are assigned an lv value of to avoid getting trapped in local optima while running the em algorithm to learn a set of latent models, other authors adapted a simulated annealing approach. hwang and co-workers' approach is the only one we are aware of that succeeds in processing high-dimensional data: in an application dealing with a microarray dataset, more than genes have been processed for around samples. to the best of our knowledge, no running time was reported for this study. nevertheless, the twofold binarity restriction and the lack of control for information decay as the level increases are severe drawbacks to reach our aims: i.e. to achieve realistic snp dependence modeling and perform subsequent association study with sufficient power. graphical models for ld modeling to address ld modeling through a probabilistic graphical model framework, various models were proposed: hidden markov models, markov random fields and bayesian networks with or without latent variables. hmms represent simple but efficient models to partition a snp sequence into blocks, because no structure learning step is required and the latent states may represent common haplotypes. verzilli and co-workers modeled snp dependences using markov random fields. they designed an mcmc method to sample over the space of possible graphs while exploiting prior biological knowledge. their approach allows to discover cliques of dependent snps, to further allow the identification of causal relations between markers and the disease status indicator. to implement a tractable method for genome-wide data, verzilli and co-workers reduce the space of possible graphs by specifying a maximal physical distance between snps belonging to the same clique, as well as a maximal size of snps for any clique. in the family of bayesian networks without lvs, haploblock implements a statistical model of haplotype block variation. this model's advantage lies in integrating population genetics concepts such as recombination hotspots, bottleneck, genetic drift and mutations. another method, bntagger, was developed for snp tag selection; it exploits conditional independence between variables. to learn the structure, bntagger implements a greedy search with random restarts; then it determines a subset of independent and highly predictive snps. the two latter methods were only tested on a small number of snps and the authors reported running times of h for snps and between and h for only snps. thus, these methods do not seem fitted to gwas data processing. regarding the family of bayesian networks with lvs, nefian modeled snp dependences through embedded bayesian networks. her model is indeed a set of lcms augmented with snp-snp dependences and lv-lv dependences. to learn the model, the snp data sequence is split into contiguous windows of fixed common size. then, for each window, an lcm is created. the lack of flexibility of the snp partitioning method used remains a severe draw-back. zhang and ji also proposed to model ld through a set of lcms, using an sem strategy. their method does not require splitting the sequence into fixed-size windows. nevertheless, the number of lcms has to be specified. as far as we know, no execution times were reported for the two latter approaches when run on high-dimensional data. other methods are based on regularization, such as the graphical lasso, and have been applied to learning sparse pgms for proteomics or gene expression studies, whose data dimensionality is high but lower than that of genome-wide data. the basic idea is to consider that the observations follow a multivariate gaussian distribution with mean µ and covariance matrix Σ. if the ijth partial correlation coefficient of the precision matrix Σ- equals zero, variables i and j are conditionally independent, given the other variables. the use of lasso aims at restraining the learning task to sparse pgms through finding a least-square solution under the following constraint: ∑ν|βν| ≤ t, meaning that the sum, over the whole variable set, of the absolute values of the regression coefficients v has to be inferior or equal to a constant t. this lasso-based approach has been extended to the case where it is reasonable to assume that the variables can be clustered into groups sharing similar correlation patterns and where sparse block-structured precision matrices are estimated. to our knowledge, verzilli et al's method is the only one whose tractability regarding gwas data is known. however, in practice, their mrf modeling reveals a drawback. the ld is modeled through cliques containing a maximum number of snps, whereas, generally, several tens or hundreds of snps may be dependent. furthermore, no dependences between cliques are taken into account. in contrast, bns with lvs offer a crucial advantage over other models: they provide synthesizing variables useful to reduce data dimensionality. however, when the number of variables exceeds several hundreds, implementing the sem approach for ld modeling leads to prohibitive computational burden. when dealing with genome-wide data, the imperious requirement for tractability leads us to choose a hierarchical clustering approach, in the line of hwang and co-workers. methods motivation of the fhlc model for gwass the hlcms offer several advantages for gwass. first, beside data dimensionality reduction, they allow a simple test of direct dependence between an observed variable and a target variable such as the phenotype, conditional on the latent variable, parent of the observed variable. note that the phenotype variable is not included in the hlcm. in the context of gwass, this test helps find the markers which are directly associated with the phenotype, i.e. causal markers, should there be any. second, hlcms can deal with the fuzzy nature of ld blocks. indeed, hlcms can take into account various degrees of ld strength between any two snps, depending on the height of their lowest common lv node ancestor in the tree. thirdly, the hierarchical structure allows zooming in through narrower and narrower regions in search for a stronger association with the disease, thus offering a data mining tool. this zooming process ends pointing out the potential markers of interest. finally, the latent variables may be interpreted in terms of biological meaning. for instance, in the case of haplotypes, that is, phased genotypes, the latent variables are likely to represent the so-called haplotype block structure of ld. to a certain extent, an lv might be interpreted as the shared ancestry of the haplotypes defined by the observed variables, namely, the contemporary haplotypes of the tree rooted in the lv. each state of an lv may represent a group of similar haplotypes. in the situation of limited ancestral recombination, similar haplotypes tend to share recent common ancestry. although this situation is not guaranteed along the genome, it is very likely for low-level lvs, since they are expected to cover very small genomic regions showing strong ld. thus the directed edges, lv → snp, can represent causal effects and provide a biological sense. besides, it has to be noted that when the latent variables capture dependences between distant snps, they can be viewed as population structure. however, snp dependences would better be more wisely modeled through a forest of hlcms. in the case of a forest, higher-order dependences are captured only when relevant, i.e when meeting a strength criterion. therefore, fhlcms allow to model a larger set of configurations than hlcms do. typically, an hlcm is limited to represent clusters of close dependent snps. actually, in this model, variables are constrained to be dependent upon one another, either directly or indirectly. consequently, hlcms cannot account for potential independence between groups of distant snps or snps located on different chromosomes. but realistic modeling requires a more flexible framework. for instance, the ld plot of the mb sequence shown in additional file reveals that the greatest part of ld is observed between snps in vicinity. ld rarely exceeds kb between snps. an fhlcm consists of a directed acyclic graph also called the structure whose nonconnected components are trees, and of θ, the parameters. figure illustrates a possible structure for an fhlcm. principle of fhlc model construction our method can process both genotypic and haplotypic data. it takes as an input a matrix dx defined on a finite discrete domain, say { 2} for unphased snps or { 1} for phased snps, describing n individuals through p variables x = {x..., xp}. algorithm cfhlc yields an fhlcm, that is a forest structure and θ, the parameters of a set of a priori distributions and local conditional distributions allowing the definition of the joint probability distribution. two search spaces are explored: the space of directed forests and the probability space. in addition, the whole set of latent variables h of the fhlcm is output, together with the associated imputed data matrix. to handle high-dimensional data, our proposal combines two strategies. the first strategy splits up the genome-scaled data into contiguous regions. in our case, splitting into windows is not a mere implementational trick; it satisfies biological grounds: the overwhelming majority of dependences between genetic markers is observed for close snps. the user interested in taking into account long-range ld due to the presence of population structure will be faced with the following choices: adjusting the window size, relying on biological background defining the maximum physical distance between snps in long-range ld; slightly diminishing the density of the studied snp sequence. a combination of these two approaches may be more convenient. then, an fhlcm is learnt for each window in turn. within a window, subsumption is performed through an adapted ahc procedure: at each agglomerative step, a partitioning method is used to identify clusters of variables; each such cluster is intended to be subsumed into an lv, through an lcm. for each lcm, parameter learning and missing data imputation are performed. a global schema of our method is presented in figure along with a hierarchy-based proposal of hwang and collaborators developed for gene expression studies, our method also implements data subsumption, meeting the two following additional requirements: a more flexible thus more faithful modeling of underlying reality, a control of information decay due to subsumption. node partitioning following martin and vanlehn, ideally, we would propose to associate a latent variable with any clique of variables in the undirected graph of dependence relations. in the case when introducing an additional lv increases a scoring function such as the bic score, the lcm is validated. however, searching for such cliques is an np-hard task. moreover, in contrast with these authors' objective, fhlcms do not allow clusters to have more than one parent each: non-overlapping clusters are required for our purpose. thus, an approximate method solving a clique partitioning problem when provided with pairwise dependence measures is relevant; the clique partitioning problem consists in finding the best partition of a graph into cliques. an algorithm meeting this purpose has already been described in the literature: bendor and co-authors designed cast, a clique partitioning algorithm devoted to variable clustering. they especially applied cast for gene expression clustering. as an input, cast requires a binary similarity matrix. the adaptation of cast to our case is straightforward: the dependence measure between two snps, evaluated through mutual information, is used as a similarity measure. all mutual information values less than a threshold tmi are assigned a similarity value of whereas the others are assigned a value of as a threshold tmi, the median value of the mutual information matrix can be used. then, the cast algorithm constructs the clusters one at a time. the authors define the affinity a of an element x to be the sum of similarity values between x and the elements present in the current cluster copen. x is an element of high affinity if it verifies inequality a ≥ tcast|copen|, where tcast is a specified similarity threshold. otherwise, x is considered an element of low affinity. to summarize, the algorithm alternates between adding high affinity elements to copen and removing low affinity elements from it. when the process stabilizes, copen is closed. a new cluster can be started. determining cardinalities for latent variables a steep task is choosing - ideally optimizing - the cardinality of each lcm's latent variable. this problem cannot be remedied using greedy search because of its intractability regarding high data dimensionality. although the regularization method of wang and collaborators has the advantage of being very quick, in our context, their method is impracticable. for instance, let us consider an hlcm learned from genome-wide data. the first layer of the hlcm contains the majority of the lvs in the model. in our case, an lv in the first layer can subsume more than child ovs. as the cardinality is the same for all ovs, the resulting cardinality of the lv after regularization remains generally very large. for example, for an lcm containing ovs {x,..., x } of cardinalities equal to the cardinality of the lv h would be: |h| = 310/ = = the simplest solution remains to arbitrarily set a small value for lv cardinalities, but it has several drawbacks. instead of using an arbitrary constant value common to all latent variables, we propose that the cardinality be estimated for each latent variable through a function of the underlying cluster's size. the rationale for choosing this function is the following: the more child nodes a latent variable has, the larger the total number of possible combinations is for the values of the child variables and the larger also is the expected number of such combinations observed over all individuals. therefore, the cardinality of this latent variable should depend on the number of child nodes. nonetheless, to keep the model complexity within reasonable limits, a maximum cardinality is fixed. parameter learning and imputation parameter learning is carried out step by step, each time generating additional latent variables and imputing their values for each individual. at ith step, this task simply amounts to performing parameter learning for as many lc models as there are clusters of variables identified. we recall that the nodes in the topology of an lcm are reduced to a unique root and leaves. therefore, at ith step, each lcm's structure is rooted in a newly created latent variable. when latent variables are the source nodes in a bn, parameter learning may be performed through a standard em procedure. this procedure takes as an input the cardinalities of the latent variables and yields the probability distributions, that is, prior distributions for those nodes with no parents and distributions conditional to parents for the remaining nodes. after imputing the missing data corresponding to latent variables, new data are available to seed the next step of the fhlcm construction: latent variables identified through step i will be considered as observed variables during step i + it has to be noted that designing an imputation method to infer the values of the latent variable for each individual is a matter for investigation. once the prior and conditional distributions have been estimated for a given lcm, probabilistic inference in bns may be performed. a straightforward way would consist in imputing the latent variable value for each individual as follows: h*=argmaxh{p}. however, in the framework of probabilistic models, this deterministic approach is disputable. in contrast, a more convincing alternative will draw a value h for latent variable h, knowing the probabilities p for each individual. controlling information decay conversely to hwang and co-workers' approach, which mainly aims at data compression, information decay control is required: in step i, any candidate latent variable h which does not bear sufficient information about its child nodes must be invalidated. as a consequence, such child nodes will be seen as isolated nodes in step i + let us consider two variables x and h. basically, the mutual information measures the difference of entropies between the independent model p p and the dependent model p p: ℐ=+ℋ)-+ℋ)=ℋ-. therefore, the mutual information measures the dependence of the two variables. the larger the difference between entropies, the higher is the dependence. now, let us consider a set of child variables x= {x x..., xn} and the parent variable h. in our case, we want to compare the two models: p p... p p and p p... p p. thus, Δ, the difference of entropies between the two models is: +ℋ)-+ℋ)=∑i=1n-ℋ)=∑i=1nℐ. Δ corresponds to the sum of mutual information values over all lcm's edges. normalization through entropy and averaging are performed to provide a more intuitive criterion: c=1sh∑i ∈ clusterℐmin,ℋ), with sh the size of cluster. c represents the average percentage of information captured by the lv with respect to its child variables. algorithm the sketch of cfhlc is presented in algorithms and the user may tune seven parameters. window size s specifies the number of contiguous snps - or variables - spanned per window. the aforementioned criterion c is meant to estimate information decay, thus allowing information dilution to be constrained to a minimal threshold t. parameters a, b and cardmax participate in the calculus of the cardinality of each latent variable. finally, parameter partitioningalg enables flexibility in the choice of the method dedicated to clustering highly-correlated variables into non-overlapping groups. within each window i, the ahc process is initiated from the first layer consisting of univariate models. each such univariate model is built for any observed variable in the set wi. the ahc process stops if each cluster identified is reduced to a singleton or if no cluster of size strictly greater than could be validated. each cluster containing at least two nodes is subject to lcm learning followed by validation. in order to simplify the fhlcm learning, the cardinality of the latent variable is estimated as an affine function of the number of variables in the corresponding cluster. algorithm learn_latent_class_model is plugged into this generic framework. after validation through threshold t, the lcm is used to enrich the fhlcm associated with the current window: a specific merging process links the additional node corresponding to the latent variable to its child nodes, themselves already present in the fhlcm structure under construction; the prior distributions of the child nodes are replaced with distributions conditional on the latent variable. the newly created latent variable, ljk, is added to the set of latent variables, whereas its imputed values, d, are stored. in wi, the variables in cℓjk are now replaced as a whole with the corresponding latent variable; data matrix d is updated accordingly. in contrast, the nodes in unvalidated clusters are kept isolated for the next step. finally, the collection of forests, dag, is successively augmented with each forest built within a window. in parallel, due to assumed independence between windows, the joint distribution of the final fhlcm is merely computed as the product of the distributions associated with the windows. input: x, a set of p variables, dx, the corresponding observations for n individuals, s, a window size, c, a criterion designed to estimate information decay while building the fhlcm, t, a threshold used to constrain information dilution, based on criterion c, partitioningalg, an algorithm dedicated to partition a set of variables into non-overlapping clusters of variables, a; b and cardmax, parameters used to estimate the cardinality of latent variables. output: dag and θ, respectively the dag structure and the parameters of the fhlcm constructed, l, the whole set of latent variables identified through the construction, dl, the corresponding data imputed for n individuals. 1: nbw← p/s /* computation of the number of contiguous windows */ 2: dag ←∅; θ ← ∅; l← ∅; dl← ∅ 3: 4: for i = to nbw 5:    /* processing of layer */ 6:    wi←{x×s+...,xi×s};d←d 7:    {∪j∈widagunivj,∪j∈wiθunivj}←learn_univariate_models 8:    dagi←∪j∈widagunivj;θi←∪j∈wiθunivj 9: 10:    step ← 1 11:    while true 12:    {cℓ ..., cℓ#c}←partition 13:    if all clusters cℓq are singletons then break end if 14: 15:    {cℓj ...,cℓj#c2}←identify_clusters_of_size_strictly_greater_than_one 16:    nbvalidclusters ← 0 17: 18:    for k = to #c2 19:    card lv ← min + b; cardmax) 20:    {dagjk, θjk, ljk, {dagjk,θjk,ljk,d}←learn_latent_class_model 21: 22:    /* validation of current cluster - see subsection controlling information decay */ 23:    if ≥t) 24:       incr 25:       dagi ← merge_structures; θi ← merge_parameters 26:       l←l∪ljk;dl←dl∪d 27:       d←∪d;wi←∪ljk 28:    end if 29: end for 30: 31: if then break end if 32: 33: incr 34: end while 35: 36: dag ← dag ∪ dagi; θ ← θ × θi 37: end for algorithm 1: cfhlc input: cℓu: a cluster containing at least two nodes, d: the corresponding observations for n individuals, cardlv: the cardinality of the latent variable to be created. output: a latent class model described by: dagu and θu, respectively the structure and the parameters of the latent class model, lu, a latent variable, d, the data imputed for the latent variable. 1: lu ← create_latent_variable() 2: dagu ← build_structure_of_latent_class_model 3: θu ← run_standard_em 4: d ← impute_data algorithm 2: learn latent class model experimental  implementation algorithm cfhlc has been developed in c++, relying on the probt library dedicated to bns http://bayesian-programming.org. we have plugged into cfhlc a c++ implementation of cast based on the original implementation provided in java by ben fry http://benfry.com/clustering/. regarding the visualization of the dags, the software tulip http://tulip.labri.fr/tulipdrupal/ was chosen, meeting both high representation quality and compactness requirements. cfhlc was run on a standard pc. experimental protocol the performance of the fhlcm-based method is evaluated using real phased and unphased genetic data, on the one hand, and simulated phased and unphased genetic data on the other hand. for real data analysis, the well-known daly et al. dataset, available at http://www-genome.wi.mit.edu/humgen/ibd5/index.html, was used. this dataset consists of trios, each composed of two parents and one child. for each individual, snps are genotyped in the 5q region and cover kb. we only analyzed the child data. regarding simulated data, two well known programs were used: hapgen and hap-simu. using hapgen http://www.stats.ox.ac.uk/~marchini/software/gwas/hapgen.html, we generated unrelated individuals for a several hundreds kb region containing around 20- snps. the haplotypes used as references come from the hapmap phase ii and concern u.s. residents of northern and western european ancestry http://hapmap.ncbi.nlm.nih.gov/. five sequences showing variable ld degrees ) ranging from to were generated. hapsimu http://l.web.umkc.edu/liujian/ was used to simulate genotypes with simulation parameters described in additional file three sample sizes were chosen with respect to the number of observed variables: k, k and k snps. for each sample size, twenty benchmarks were generated. for these experimentations, imputation of lvs' values was achieved by assigning the most probable values given the observations. a drawback of this method is the loss of probabilistic relation between a variable and its parent variable. a definite advantage lies in its running in around half the time required by imputation through simulation. ld modeling real data regarding the daly benchmark, our aim was to evaluate how the forest obtained keeps up with the real structure of the biological data. moreover, the cfhlc algorithm was compared to four other methods. we learned fhlcms on both haplotype and genotype data. the corresponding graphs are displayed in figures and respectively. globally, the two graphs are similar: most of snps which are connected through an lv in the haplotype-data graph are also connected through an lv in the genotype-data graph, e.g. snp snp and snp moreover, a substantial part of these snps share a common parent in both graphs: for instance, in both hdg and gdg, we observe that snp and snp are linked by an lv belonging to layer thus, learning fhlcm from genotype data instead of haplotype data leads to similar hierarchical structures. however, on average, we observe that the snps in the gdg are more connected: and connected components are identified in the hdg and the gdg, respectively. for example, the two framed trees and in the hdg of figure are linked by a high-level lv in the gdg of figure. we expect that the fhlcms' graphs will reflect the "haplotype block structure": large blocks of correlated contiguous snps separated by recombination hotspots. first, we observe that the physical position of snps influences their connection, since close snps tend to be linked by an lv belonging to a low layer, whereas distant snps are generally connected by a high-level lv. however, strong dependences between distant snps are also observed, e.g. between snp and snp or snp and snp. this characteristic reveals that the ld structure is not only dominated by spatial effects and justifies our haplotype cluster approach. in addition, the graphs interestingly show trends consistent with biological reality, that is the variation of the recombination rates inferred by software phase v along the studied sequence. indeed, most of subtrees rooted in low-level lvs cover regions with low recombination rates. more than 68% and 94% of lvs from layer cover chromosomic segments showing rrs below cm/mb and cm/mb, respectively. the same tendency is observed for more than 44% and 66% of lvs from layer respectively. these results show the relevance of interpreting low-level lvs as haplotype shared ancestry when cfhlc's input is haplotype data. we compared the structure obtained by cfhlc with those output by four other approaches: daly et al.'s method, gerbil, haploblock and zhang et al.'s algorithm. all these methods were detailed in section state of the art. the three former methods partition the sequence into blocks of contiguous snps. in contrast, the latter algorithm yields clusters of non-contiguous snps. we recall that cfhlc algorithm generates a hierarchical clustering of non-contiguous snps. in figure we compare the haplotype block- or cluster-structures obtained through all five methods aforecited. in spite of the fact that these methods differently tackle ld modeling, common trends emerge. for instance, the last block identified by daly et al.'s method, gerbil and zhang et al.'s algorithm is also inferred by our algorithm in line slight differences are observed with the two first blocks resulting from daly et al.'s method and gerbil which only form one block for zhang et al.'s algorithm and cfhlc. compared to other methods, most divergences with our algorithm remain in its unique ability to take into account the fuzzy delimitations of clusters. this is illustrated with the central area of the sequence, which actually presents two weak recombination hotspots. another difference with the other methods is the presence of "unclustered" snps, like snp snp and snp in our model. the running times, dimension reduction rates and entropy compression rates of all methods are reported in table results show that gerbil is the fastest algorithm tested, with a running time of s. however, cfhlc and zhang et al.'s algorithm, which learn more complex models, achieve their tasks in quite a reasonable time, s and s, respectively. compared to others, haploblock is the slowest method, with a running time of mn, due to the high complexity of learning models based on population genetics. for the three methods exhibiting a partition of contiguous snps, we defined the dimension reduction rate as the ratio of the number of blocks to the number of snps. as regards zhang et al.'s algorithm, the drr was defined as the ratio of the number of clusters to the number of snps. in the case of cfhlc, we consider that the information of each fhlcm's tree can be synthesized by its root, providing the best dimension reduction. therefore, in this case, the drr is defined as the number of roots in the whole forest divided by the number of snps. haploblock generates the lowest number of blocks with an average of, whereas zhang et al.'s algorithm partitions the sequence in clusters, and daly et al.'s method and gerbil both identify blocks. cfhlc presents the lowest dimension reduction with trees, due to the presence of "unclustered" snps: snp snp snp snp snp snp and snp as an alternative measure of compression, we defined the entropy compression rate as the ratio of the sum of block entropies in a partition to the entropy, assuming no structure. we observe a different ranking of the methods. we notice that cfhlc and zhang et al.'s algorithm, which both learn cluster models, provide the best ecr values, whereas haploblock, gerbil and daly et al.'s method show ecr values of, and, respectively. regarding the ecr criterion, the comparatively better results obtained for the two cluster models are explained by the absence of the constraint for compulsory physical proximity between snps. moreover, the ecr criterion does not penalize anymore the cfhlc algorithm, since the unclustered snps contribute relatively little to the overall information content. we ran the last three programs on a standard computer. as we had no access to daly et al.'s software, we could only compare the dimension reduction rates and entropy compression rates calculated from their results with the dimension reduction rates and entropy compression rates obtained with the other methods. in subsection motivation of the fhlc model for gwass, we argued that the multiple layers of an fhlcm can describe various degrees of ld strength. to analyze this property, we plotted the r squared correlation coefficient of any pair of snps against the level of their most recent common ancestor. figure and show such plots drawn for haplotype and genotype data, respectively. starting from values in the range, the r correlation coefficient quasi-linearly decreases when the mrca level increases. we conclude that the layered structure of the fhlcm faithfully reflects ld strength variety. these encouraging results lead us to visually compare the ld plot and the triangular matrices of the mrca levels for haplotype and genotype data, as presented in figure for this purpose, the same color code was used in the ld plot and the triangular matrix of mrca levels. in the ld plot, the color of each cell varies with the r value. since the median value of r can be computed for each mrca level, the color of each cell in the triangular matrix of mrca levels is set, relying on the color scale used for the ld plot. we expected to observe a correspondence between the three plots. the outstandingly clear correspondence demonstrates the ability of fhlcms to accurately model multiple levels of ld strength: the overall majority of the ld plot dependences are also present in the mrca level matrix. interestingly, modeling from genotype data leads to quite good results compared to haplotype data. haplotype diversity is generally very low within haplotype blocks or clusters. in our hierarchical model, haplotype diversity is expected to be all the larger within a cluster as the level of the lv subsuming this cluster is high. to check this point, we have relied on the clusters of figure for each lv, haplotype diversity has been calculated as the number of the most common haplotypes observed at level l of the tree rooted in this lv. figure plots the number of the most common haplotypes against the lv level. the plot shows that the haplotype diversity median remains very low for the first four layers and dramatically increases to around in the fifth layer. these results confirm our expectation relative to haplotype diversity in fhlcms. simulated data finally, the impact of varying ld degrees was studied. for this purpose, we generated haplotype data with the hapgen software. five sequences, showing variable ld degrees ) ranging from to, were used to learn fhlcms. figure shows the forests obtained. the forests reveal that increasing ld degrees entails higher graph connectivity as well as a larger number of layers. indeed, when median equals, connected components are identified and the highest lv belongs to the third layer. conversely, in the case when median is equal to, the forest is only composed of connected components and the highest lv belongs to layer thus, we conclude that cfhlc can process sequences with various ld degrees and generate fhlcms whose structures reliably reflect linkage disequilibrium. scalability for gwass scalability has been studied through the data simulated with hapsimu. in the hardest case, additional file shows that only hours are required with a window size s set to for the same dataset processed in the cases "s = 200" and "s = 600", running times are h and h, respectively. for the same number of ovs, wang et al. report running times of about two months. regarding the k case, running times are h, h and h for "s = 100", "s = 200" and "s = 600", respectively. further analysis of cfhlc algorithm finally, many other experimentations are reported with their commentaries in additional files and additional file focuses on the impact of window size on running time. in additional file data dimension reduction is evaluated from the distribution of lvs over the forest's layers. additional files and study the impact of window size on the number of roots in the forest, the number of lvs, the number of layers and the distribution of lvs over the forest's layers, respectively. additional files and analyze how information fades while the layer number increases. an important result is that cfhlc can achieve a data dimensionality reduction of more than 80% of the number of observed variables. regarding spatial complexity of cfhlc algorithm, the entire fhlcm does not require to be stored in ram because each window can be saved on the hard disk.  our contribution in this paper is twofold: a new framework has been described, which was tested on both real and simulated data and was proven able to concisely model ld and to reduce snp data dimensionality; cfhlc, an algorithm dedicated to learn fhlcms, has been shown to be efficient when run on genome-scaled benchmarks. compared to verzilli and co-workers' works, our algorithm provides a more accurate modeling of ld and synthesizes genetic marker information through lvs. in addition, unlike nefian or zhang and ji, our method does not require to specify the number of lcms and can capture multiple levels of dependences, thus taking into account the fuzzy nature of ld. to our knowledge, our hierarchical method is the first one shown to achieve fast model learning for genome-scaled data sets, while maintaining satisfying information scores and relaxing the twofold binarity restriction of hwang and collaborators' model. hwang and collaborators' purpose is only data compression. we are faced with a more demanding challenge: to make a sufficiently powerful down-stream association analysis possible. in discussing the biological interpretation of latent variables, we mentioned the potentiality of fhlcms for population substructure description. in essence, using hierarchical models is highly appealing to take into account the long-range ld expected in substructured populations. however, this interesting use of such hierarchical models as fhlcms is somewhat precluded by the technical necessity to partition the genome into small regions. as a first palliative, we indicated two strategies to cope with this current technical limitation. however, the strong expectation for faithful substructure modeling through fhlcms advocates further efforts to clear the hurdle on path to realistic long-range ld modeling. a bottleneck currently lies in the clique partitioning method chosen, which forbids window sizes encompassing more than observed variables. in addition to investigating alternative partitioning methods, a lead to cope with this bottleneck may be to adapt the specific processings at the limits of contiguous windows or use overlapping windows. in short, fhlcms can be used to resolve several major problems in the gwass' context. beside flexible data dimension reduction through fhlcms' lvs, fine mapping of causal snps is expected thanks to conditional independence properties encoded in such models. for instance, fhlcms' lvs can be used to condition tests for independence between a snp and the phenotype. moreover, due to their hierarchical structure, fhlcms represent an original and appropriate solution to study long-range ld in substructured populations, a recurring problematic in gwass. finally, genome-wide visualization of ld can be easily achieved with these models using a graph visualization tool and will provide an intuitive representation of snp - snp dependences as well as information synthesis through latent variables. in this current version of cfhlc, when processing haplotype data, we were not interested in knowing the sequence of each ancestral haplotype. we just wanted to know from which ancestral haplotype a contemporary haplotype comes. nevertheless, it is feasible to infer the sequence of each ancestral haplotype using probabilistic inference. finally, although our modeling is designed for gwas data, we emphasize that it could be applied to other data presenting spatial dependences between variables, in particular sequential data. beyond this specific case, it would be interesting to assess the model's generality, in order to determine if it can be applied to generic graphical model learning. authors' contributions rm and cs both wrote the manuscript. rm proposed the modeling of the genetic data, the conception of a tractable algorithm and carried out the implementation and the experiments. cs designed the study, participated in its coordination and partook in the conception of the algorithm. pl participated in the design of the study, the conception of the algorithm, helped to improve the manuscript and provided indispensable feedback. all authors read and approved the final version of the manuscript. supplementary material additional file 1 direct and indirect associations between a genetic marker and the phenotype. the figure included into this additional file illustrates the cases of direct and indirect associations between a genetic marker and the phenotype. click here for file additional file 2 linkage disequilibrium plot for a simplified haplotype block structure. the figure included into this additional file describes a standard representation of pairwise dependences between genetic markers. click here for file additional file 3 linkage disequilibrium plot of a real kb snp sequence. the figure presented in this additional file shows the linkage disequilibrium plot of a real kb snp sequence. click here for file additional file 4 hierarchical latent class model. the figure presented in this additional file depicts a hierarchical latent class model. click here for file additional file 5 linkage disequilibrium plot of a mb snp sequence. the figure included in this additional file describes the linkage disequilibrium plot of a mb snp sequence. click here for file additional file 6 parameter value adjustment for the generation of simulated genotypic data through software hapsimu. the table included in this additional file enumerates the values chosen for the parameters of software hapsimu. click here for file additional file 7 average running time versus number of variables. the figure presented in this additional file plots the running time of the cfhlc algorithm versus the number of snps in the dataset. click here for file additional file 8 impact of window size on running time. the figure presented in this additional file plots the running time of the cfhlc algorithm versus the window size. click here for file additional file 9 number of variables per layer over the whole fhlc model. the figure included in this additional file describes the average distribution of the variables over the layers. click here for file additional file 10 impact of window size on the number of roots. the figure included in this additional file depicts the impact of window size on the number of roots. click here for file additional file 11 impact of window size on the number of latent variables. the figure presented in this additional file shows the impact of window size on the number of latent variables. click here for file additional file 12 impact of window size on the number of layers. the figure presented in this additional file describes the impact of window size on the number of layers. click here for file additional file 13 impact of window size of the number of latent variables per layer and on the ratio of the number of latent variables per layer to the total number of variables. the two subfigures included in this additional file depict the impact of window size on the number of latent variables per layer on the one hand and the impact of window size on the number of latent variables per layer to the total number of variables, on the other hand. click here for file additional file 14 impact of window size on scaled mutual information, per layer. the figure presented in this additional file describes the impact of window size on scaled mutual information, per layer, over the whole fhlc model. click here for file additional file 15 average scaled mutual information per layer over the whole fhlc model; impact of parameters a and b. the figure presented in this additional file shows the impact of parameters a and b on scaled mutual information, per layer, over the whole fhlc model. click here for file acknowledgements the authors are grateful to two anonymous referees for constructive comments and help in improving their manuscript. the authors would also like to thank christian dina and jean-jacques schott for useful discussions regarding genetics and statistical genetics. yulong zhang has to be thanked for his precious help in the comparison between his algorithm and cfhlc algorithm. rm's work was supported by a grant from the bil bioinformatics research project of pays de la loire region, france.
the forest of hierarchical latent class models offers several advantages for genome-wide association studies: accurate modeling of linkage disequilibrium, flexible data dimensionality reduction and biological meaning borne by latent variables. in addition, the data dimensionality can be reduced through the latent variables of the model which synthesize the information borne by genetic markers. this model offers an adapted framework to deal with the fuzzy nature of linkage disequilibrium blocks. discovering the genetic basis of common genetic diseases in the human genome represents a public health issue. in order to tackle the learning of both forest structure and probability distributions, a generic algorithm has been proposed. we present an accurate modeling of dependences between genetic markers, based on a forest of hierarchical latent class models which is a particular class of probabilistic graphical models. however, the dimensionality of the genetic data and its complexity make the statistical analysis a challenging task.
 discovering the genetic basis of common genetic diseases in the human genome represents a public health issue. however, the dimensionality of the genetic data and its complexity make the statistical analysis a challenging task.  we present an accurate modeling of dependences between genetic markers, based on a forest of hierarchical latent class models which is a particular class of probabilistic graphical models. this model offers an adapted framework to deal with the fuzzy nature of linkage disequilibrium blocks. in addition, the data dimensionality can be reduced through the latent variables of the model which synthesize the information borne by genetic markers. in order to tackle the learning of both forest structure and probability distributions, a generic algorithm has been proposed. a first implementation of our algorithm has been shown to be tractable on benchmarks describing variables for individuals.  the forest of hierarchical latent class models offers several advantages for genome-wide association studies: accurate modeling of linkage disequilibrium, flexible data dimensionality reduction and biological meaning borne by latent variables. 