BACKGROUND
the chromatin immunoprecipitation  coupled with microarray  analysis, provides the researchers an efficient way of mapping protein-dna interactions across a whole genome.
the chip-chip technology has been used in a wide range of biomedical studies, such as identification of human transcription factor binding sites, investigation of dna methylation, and investigation of histone modifications in animals and plants  <cit> .
data from chip-chip experiments encompass dna-protein interaction measurements on millions of short oligonucleotides  which often tile one or several chromosomes or even the whole genome.
the data analysis consists of two steps:  identifying the bound regions where dna and the protein are cross-linked in the experiments; and  identifying the binding sites through sequence analysis of the bound regions.
the goal of this paper is to develop an effective method for the first step analysis.
analysis of the chip-chip data is very challenging, due to the large amount of probes and the small number of replicates.
the existing methods in the literature can be roughly grouped into three categories, the sliding window methods  <cit> , the hidden markov model  methods  <cit> , and the bayesian methods  <cit> .
other methods have been suggested, e.g., by zheng  <cit> , huber  <cit>  and reiss  <cit> , but are less common.
the sliding window methods are to test a hypothesis for each probe using the information from the probes within a certain genomic distance sliding window, and then try to correct for the multiple hypothesis tests.
the test statistics used are varied.
cawley  <cit>  used wilcoxon's rank sum test, keles  <cit>  used a scan statistic which is the average of t-statistics within the sliding window, and ji and wong  <cit>  used a scan statistic which is the average of empirical bayesian t-statistics within the sliding window.
since each test uses information from neighboring probes, the tests are not independent, rendering a difficult adjustment in the multiple hypothesis testing step.
we note that the power of the sliding window tests is usually low, especially for the tests for the regions where the probe density is low.
this is because there will be only very limited neighboring information available for those tests.
since, in the chip-chip experiments, the dna samples hybridized to the microarrays are prepared by pcr, which is known to perform independently of the form of dna, the far probes should have similar intensity patterns as long as they are of similar positions to their nearest bound regions.
this provides a basis for us to devise powerful methods that make use of information from all probes.
the hmm methods have the potential to make use of all probe information, where the model parameters are estimated using all available data.
however, in most of the existing implementations of hmms, the model parameters are estimated in an ad hoc way.
for example, li  <cit>  estimated the model parameters using previous results on affymetrix snps arrays.
an exception is tilehmm  <cit> , where the model parameters are estimated using the baum-welch and viterbi training algorithms  <cit> .
however, as pointed out by humburg et al.
<cit> , the baum-welch algorithm tends to converge to a local maximum of the likelihood function, and the viterbi training algorithm even fails to converge to a local maximum of the likelihood in some cases.
this renders that the parameter estimates and thus followed inference often suboptimal to the problem.
bayesian methods have also the potential to make use of all probe information.
like the hmm methods, the bayesian methods estimate the model parameters using all available data.
however, these methods usually require multiple replicates or some extra experimental information to parameterize the model.
for example, the joint binding deconvolution model  <cit>  requires one to know the dna fragment lengths, measured separately for each sample via extrophoretic analysis; and the hierarchical gamma mixture model  <cit>  requires one to first divide the data into genomic regions containing at most one bound region, but such information is, in general, unavailable.
the bayesian hierarchical model  <cit>  models the probe intensities using essentially a mixture of normal distributions, and models the spatial structure of the probes using a gaussian intrinsic auto-regression model  <cit> .
gottardo  <cit>  developed a software for the model, bayesian analysis of chip-chip .
using bac  <cit>  does not need extra experimental information, but it is extremely slow, roughly  <dig> hours for a dataset with  <dig>  probes on a personal computer.
one reason for the slow speed is the use of mcmc simulations.
in this paper, we propose a bayesian latent variable model for tiling array data.
our method differs from the existing bayesian methods, such as the joint binding deconvolution model  <cit> , the hgmm  <cit> , and the bayesian hierarchical model  <cit> , in several respects.
firstly, it works on the difference between the averaged treatment and control samples.
this enables the use of a simple model for the data, which avoids the probe-specific effect and the sample  effect.
as a consequence, this enables an efficient mcmc simulation of the posterior distribution of the model, and also makes the model rather robust to the outliers.
secondly, it models the neighboring dependence of probes by introducing a latent indicator vector.
thirdly, it does not require multiple replicates or extra experimental information.
as described below, it can work on a single intensity measurement for the probes.
the bayesian latent model has been successfully applied to several real and ten simulated datasets, with comparisons with some of the existing bayesian methods, hidden markov model methods, and sliding window methods.
the numerical results indicate that the bayesian latent model can outperform the others, especially when the data contain outliers.
our method is also computationally efficient; it takes about  <dig> minutes for a dataset with  <dig>  probes on a personal computer.
the remainder of this paper is organized as follows.
in section  <dig>  we describe the bayesian latent model and its mcmc implementation.
in section  <dig>  we test the bayesian latent model on real and simulated datasets with comparisons with tilehmm, bac and some sliding window methods.
in section  <dig>  we discuss possible extensions of our methods and provide an explanation why our method outperforms tilehmm and bac via a detailed comparison of the models used by them.
in section  <dig>  we conclude the paper.
methods
consider a chip-chip experiment with two conditions, treatment and control.
let x <dig> and x <dig> denote, respectively, the samples measured under the treatment and control conditions.
each sample has ml, l =  <dig>   <dig>  replicates providing measurements for n genomic locations along a chromosome or the genome.
suppose that these samples have been normalized and log-transformed.
in this paper, we summarize the measurements for each probe by
where  is the intensity measurement of probe i averaged over ml replicates.
the underlying assumption for the summary statistic in  is that the intensity measurements for each probes has a variance independent of its genomic position.
the rationale is that the dna samples used in the experiments are prepared by pcr, which is known to perform independently of the form of dna, and that the amount of the dna samples provides the main sources for the variation of probe intensities.
we note that a similar assumption has also been made in other bayesian software, e.g., tilehmm  <cit> .
otherwise, yi can be adjusted by its standard error to a shrinkage t-statistic  <cit>  or an empirical bayes t-statistic  <cit> , depending on the estimate of the standard error.
note that both the adjustments are toward the constant variance of probes.
even with the adjustments, the bayesian latent model developed in this paper can still work reasonably well, as the normality assumption approximately holds for the modified t-statistics.
the bayesian latent model
suppose that the data consists of a total of k bound regions, and that region k consists of nk  consecutive probes.
for convenience, we call all the non-bound regions by region  <dig> and denote by n <dig>  the total number of probes contained in all the non-bound regions, although the probes in which may be non-consecutive.
thus, we have .
let z =  be a latent binary vector associated with the probes, where zi =  <dig> indicates that probe i belongs to a bound region and  <dig> otherwise.
given z, we can re-index , a realization of , by ykj, k =  <dig> ...,k, j =  <dig> ...,nk.
then ykj can be modeled as follows,
where μ <dig> is the overall mean, which models the difference of sample effects ; ν <dig> =  <dig> and νk >  <dig>  k =  <dig> ...,k accounts for the difference of probe intensities in different bound regions; ϵkjs are random errors independently and identically distributed as n. in words, model  assumes that, conditioning on the latent vector z, ykjs are mutually independent and also identically distributed within the same bound region.
we are aware that for the tiling array data, the probe intensities tend to form a peak around the true binding site.
since, given z, the order of probes is meaningless to us, the model  is appropriate if ignoring the order of the probes.
we note that a similar assumption has also been used in the hgmm and hmm methods.
conditioning on z, the likelihood of the model can be written as
to conduct a bayesian analysis for the model, we specify the following prior distributions for the model parameters:
where ig denotes an inverse gamma distribution, u denotes a uniform distribution, and α, β, νmin, νmax are hyperparameters.
in this paper, we set α = β =  <dig> , which form a vague prior for σ2; and set νmin = 2sy and νmax = maxi yi, where sy is the sample standard error of yi.
different values of νmin, e.g., sy and  <dig> sy, have also been tried in our simulations, and the results are similar.
the sensitivity issue of the bayesian latent model to the hyperparameters will be further discussed in section  <dig>  in addition, we assume that the latent vector z follows a truncated poisson distribution,
where k, denoting the total number of bound regions specified by z, and is thus a function of z; λ is a hyperparameter; kmax is the largest number of bounded regions allowed by the model; and
which makes the prior  a proper distribution.
the rationale behind this prior can be explained as follows.
since the length of each bound region is very short comparing to the chromosome or the whole genome, it is reasonable to view each bound region as a single point, and thus, following the standard theory of poisson process, the total number of bound regions can be modeled as a poisson random variable.
conditioning on the total number of bound regions, as implied by , we put an equal prior probability on all possible configurations of z, i.e., assuming a non-informative prior for z. the prior  penalizes a large value of k, where the parameter λ represents the strength of penalty.
we do not recommend to use a large value of λ, as the number of true bound regions is usually small and a large value of λ will lead to discovery of too many false bound regions.
our experience shows that a value of λ around  <dig>  usually works well for the chip-chip data.
in this paper, we set λ =  <dig>  in all simulations.
the parameter kmax is usually set to a large number.
we set kmax =  <dig> in all simulations of this paper.
as long as the value of kmax has been reasonably large, increasing it further would have a negligible effect on simulations.
finally, we would like to point out that the bound region identification problem can also be viewed as a change-point identification problem that has been widely studied in statistics.
for the change-point identification problem, the same truncated poisson prior has been used for modeling the total number of change-points by many authors, see, e.g., phillips and smith  <cit> , dension et al.
<cit> , liang and wong  <cit> , and liang  <cit> .
if ν <dig> ...,νk ∈ , combining the likelihood and prior distributions, integrating out σ <dig>  and taking the logarithm, we get the following log-posterior density function
otherwise, the posterior is equal to  <dig>
due to the design of chip-chip experiments, it is obvious that the intensity measurements of the neighboring probes are positively dependent.
to model this dependence, we use a latent indicator vector z. this makes our model different from the existing models, such as the joint binding deconvolution model  <cit> , the hgmm  <cit> , and the bayesian hierarchical model  <cit> .
both the joint binding deconvolution model and the bayesian hierarchical model model the mean of probe intensities through the gaussian random field , although their formulations may not be in the standard form of the gmf.
like the bayesian latent model, the hgmm models the mean of probe intensities by a piece-wise constant function.
the difference is that the hgmm requires one to first divide the data into genomic regions containing at most one bound regions, and thus it allows different non-bound regions to have different means.
considering the physical property of pcr, which performs independently of the form of dna, allowing different non-bound regions to have different mean values may not be necessary.
mcmc simulations
to simulate from the posterior distribution , we used the metropolis-within-gibbs sampler  <cit> ; see appendix for the details.
note that when a component of z is updated, the sum of square terms in the posterior density can be calculated in a recursive manner, and this simplifies the computation of the posterior density greatly.
inference of bound regions
let pi = p be the marginal posterior probability that probe i belongs to a bound region.
since the bound regions are expected to consist of several consecutive probes with positive ip-enrichment effects, the regions which consists of several consecutive probes with high marginal posterior probabilities are likely to be bound regions.
to identify such regions, we follow gottardo  <cit>  to consider the joint posterior probability
where i is the index of the probes, w is a pre-specified half-window size, and m is the minimum number of probes belonging to the bound region.
as explained in gottardo  <cit> , the purpose of introducing the joint posterior probability is to remove the false bound regions, which usually consists of only few isolated probes with large enrichment effects.
we found that the choice w =  <dig> and m =  <dig> works well in practice.
this choice of w is consistent with the moving window size used in other work, such as ji and wong  <cit> , keles  <cit> , and gottardo  <cit> .
the choice of m is chosen for robustness to false bound regions.
it also reflects our belief that a bound region should consist of at least five consecutive probes with large enrichment effects.
note that estimation of ρi is trivial based on the samples simulated from the posterior distribution.
the value of ρi depends on a lot of parameters, such as w, m and the hyperparameters of the model.
however, we found that the orders of ρi are rather robust to these parameters.
this suggests us to treat ρi as a conventional testing p-value, and to control the false discovery rate  of the bound regions using a fdr control method, e.g., the empirical bayes method  <cit>   or the stochastic approximation-based empirical bayes method  <cit> .
both the methods allow for the dependence between testing statistics and an empirical determination of the density of the testing statistics.
although a strict control of fdr is important to the detection of bound regions, it is not the focus of this paper.
in this paper, we will follow other bayesian methods, such as bac, to simply set a cut- off value of ρi.
we classify probe i as a probe in bound regions if ρi ≥  <dig> , and classify probe i as a probe in nonbound region otherwise.
as we will see in the numerical examples, the joint posterior probability can lead to a good detection of true bound regions.
RESULTS
the estrogen receptor data
the er data were generated by carroll  <cit> , which mapped the association of the estrogen receptor on chromosomes  <dig> and  <dig>  here we just used a subset of the data to illustrate how the bayesian latent model works.
the subset we used is available from the bac software  <cit> .
it consists of intensity measurements for  <dig> probes under the treatment and control conditions with three replicates each.
the same subset has been used by bac for a demonstration purpose.
the bayesian latent model was first applied to the dataset.
the algorithm was run  <dig> times.
each run consisted of  <dig> iterations, and cost about  <dig>  minutes cpu time on a personal computer .
all computations of this paper were done on this computer.
based on the gelman-rubin diagnostic plot  <cit>  , we discarded the first  <dig> iterations for the burn-in process, and used the remaining  <dig>  iterations for further inference.
figure  <dig> shows the estimates of the joint posterior probabilities resulted from one run.
for comparison, bac and tilehmm  were also applied to this dataset.
both bac and tilehmm produced a probability measure for each probe, similar to ρi, on how likely it belongs to a bound region.
the results were shown in figures  <dig> and  <dig>  respectively.
the comparison shows that all the three methods produced very similar results for this dataset.
however, the results produced by the bayesian latent model are neater; the joint posterior probabilities produced by it tend to be dichotomized, either close to  <dig> or close to  <dig>  this gives the user a clear classification for the bound and non-bound regions.
to provide some numerical evidence for this statement, we calculated the ratio #{i : pi >  <dig> }/#{i : pi >  <dig> }, where #{i : pi > a} denotes the number of probes with pi greater than a. here pi refers to the joint posterior probability for the bayesian latent model and bac, and the conditional probability for tilehmm.
the ratios resultant from the bayesian latent model, bac and tilehmm are  <dig> ,  <dig>  and  <dig> , respectively.
later, we assessed the sensitivity of the bayesian latent method to the values of the hyperparameters νmin and λ with other parameters fixed, α = β =  <dig>  and νmax = maxi yi.
the cross settings { <dig> sy,  <dig> sy,  <dig> sy, 2sy,  <dig> sy, 3sy} ×  for  were tried for this dataset.
for each setting, the algorithm was run  <dig> times, and each run consisted of  <dig>  iterations.
to measure the similarity of the bound regions resultant from different settings of the hyperparameters, we propose to use the adjusted rand index  <cit> .
the adjusted rand index is usually used in the literature of clustering, which measures the degree of agreement between two partitions of the same set of observations even when the comparing partitions having different numbers of clusters.
it is obvious that the problem of bound region identification can also be viewed as a clustering problem; where the genome was partitioned into a series of segments, non-bound or bound regions, and each of the segments forms a cluster.
the adjusted rand index is defined as follows.
let Ω denote a set of n observations, let c = {c <dig> ...,cs} and  represent two partitions of Ω, let nij be the number of observations that are in both cluster ci and cluster , let ni.
be the number of observations in cluster ci, and let n.j be the number of observations in cluster .
the adjusted rand index is
a higher value of r means a higher correspondence between the two partitions.
when the two partitions are identical, r is  <dig>  when a partition is random, the expectation of r is  <dig>  under the generalized hyper-geometric model, it can be shown  <cit>  that
refer to hubert and arabie  <cit>  for the theoretical justification of r.
in calculations of the adjusted rand indices for the sensitivity experiments, we used the result shown in figure  <dig> as the standard; that is, if a partition is identical to that partition, r will be  <dig>  the results are summarized in figure  <dig>  where the adjusted rand index is plotted as a function of log for different setting of νmin.
figure  <dig> shows that, for each value of νmin, the adjusted rand index varies between  <dig>  and  <dig>  as λ runs from  <dig>  to  <dig> .
this indicates that the performance of the bayesian latent model is rather robust to the choices of νmin and λ.
finally, we examined the robustness of the bayesian latent model to different choice of w and m with other parameters fixed at α = β =  <dig> , λ =  <dig> , and νmin = 2sy.
the cross settings { <dig>   <dig>   <dig>  10} × { <dig>   <dig>  7} for  were tried for this dataset.
again, the adjusted rand index is used as the similarity criterion and the result shown in figure  <dig> as the standard.
the results were summarized in table  <dig>  which indicates, for this dataset, the bayesian latent model is quite robust to the choices of w and m. in practice, to achieve robustness to outlying probes, we suggest to avoid choosing a small m. in all the following simulations, we set m =  <dig>
the average of adjusted rand indices  is calculated based on  <dig> independent runs.
the entries for the cells with 2w <m are not available .
the robustness of the results with respect to changes of α, β and νmax are not studied in the paper.
the reason is that νmax is completely determined by the data, and the values of α and β we used form a vague prior for the variance σ <dig>
p <dig> data
in a chip-chip experiment, cawley  <cit>  mapped the binding sites of four human transcription factors sp <dig>  cmyc, p53-fl, and p53-do <dig> on chromosomes  <dig> and  <dig>  the experiment consisted of  <dig> treatment and  <dig> input control arrays, and the chromosomes spanned over three chips a, b and c. refer to cawley  <cit>  for the details of the experiment.
for the testing purpose, p53-fl data on chips a, b and c were used in this paper, which contains  <dig> quantitative pcr verified regions.
as in cawley  <cit> , the data were pre-processed by filtering out the local repeats, quantile-normalized  <cit> , rescaled to have a median feature intensity of  <dig> for the purpose of adjusting batch effect, and then log-transformed.
since the normalization is not the focus of this paper, we skipped the details.
the bayesian latent method was first applied to the p <dig> data.
the data on chip a, chip b, and chip c were analyzed separately.
each run consisted of  <dig>  iterations.
diagnostic plot for the convergence of these runs indicates that they can converge within several hundreds of iterations, even the data on each chip consists of more than  <dig>  probes.
accordingly, the first  <dig> iterations were discarded for the burn-in process, and the samples from other iterations are used for further analysis.
for comparison, bac and tilehmm were also applied to this example.
given the posterior probabilities, a cutoff of  <dig>  was used for all methods to detect bound regions.
all resultant bound regions having less than  <dig> probes or  <dig> bps were considered to be spurious and removed, and those regions separated by  <dig> bps or less were merged together to form a predicted bound regions following the approach taken by cawley  <cit> .
the results were summarized in table  <dig>  although tilehmm detected all the  <dig> validated regions, it essentially fails for this example.
it identified a total of  <dig> bound regions, which should contain too many false bound regions.
we suspect that the failure of tilehmm for this example is due to its training algorithm; it is very likely that tilehmm converged to a local maximum of the likelihood function.
this have been noted by humburg et al.
<cit> , tilehmm may converge to a local maximum of the likelihood function with either the baum-welch algorithm or the viterbi training algorithm, rendering an ineffective inference for the model.
both the total numbers of regions and quantitative pcr verified ones detected by each method on each chip are reported.
the columns under "p53" summarize the results on chips a, b and c. the number in the parentheses is the number of clusters needed to cover all  <dig> experimentally validated bound regions.
both the bayesian latent method and bac work well for this example.
at a cutoff of  <dig> , bac identified  <dig> bound regions, which cover  <dig> out of  <dig> experimentally validated bound regions.
the bayesian latent method works even better.
at the same cutoff, it only identified  <dig> bound regions, but which also cover  <dig> out of  <dig> experimentally validated bound regions.
for further comparison of the bayesian latent method and bac, we relaxed the cutoff value and counted the total number of regions needed to cover all experimentally validated regions.
we found that the bayesian latent method only needs to increase the total number of regions to  <dig>  while bac needs to increase to  <dig> regions.
note that the bac and tilehmm's results reported here may be a little different from those reported by other authors, due to the difference of the normalization methods.
simulated data
to have a careful assessment of the performance of the bayesian latent model, we simulated  <dig> datasets based on the sp <dig> data of cawley's  <cit>  experiment.
each dataset consists of  <dig>  probes, two conditions , and six replicates under each condition.
the probe genomic coordinates we used in simulations were the first  <dig>  genomic positions used in the sp <dig> data.
each dataset consisted of  <dig> bound probes, forming  <dig> bound regions.
as in gottardo  <cit> , the bound regions were assumed to describe a peak with the intensity function given by a exp{-42/b2}, where a is the amplitude of the peak, b controls the width of the peak, c represents the center of the peak, and gi is the genomic position of probe i. we also followed gottard  <cit>  to generate the centers of the bound regions randomly across the set of possible coordinates while imposing a separation of at least  <dig> bps between peaks; and to generate the values of parameter b uniformly between  <dig> and  <dig> bps.
the values of parameter a were generated uniformly between  <dig> and  <dig>  the variance of the probe intensity was estimated from the sp <dig> data.
firstly, the performance of different models is assessed using the area under the receiving operating characteristic  curve  <cit>  and the error rate.
the former is a standard measure for the performance of a multiple hypothesis testing method, which shows the true positive discovery rate  against the false positive discovery rate  at probe level.
the later is a standard measure for the performance of a classification method, which shows the proportion of totally incorrect probe calls, including both false positives and false negatives, against different cutoff values.
all the three methods, bayesian latent method, bac and tilehmm, were applied to the  <dig> full datasets.
the averaged roc curve and error rate across a range of cutoffs are obtained and plotted in figure  <dig>  as indicated by figure  <dig>  the bayesian latent method and tilehmm have very similar performances on these datasets, and both are much better than bac.
by further examining the plot on the right, which provides a closer view of the area enclosed by the dotted line and axis on the left, it is easy to see that the bayesian latent method is better than tilehmm for this example.
next, we checked the error rate for each model.
the results were shown in figure  <dig>  again, the bayesian latent method and tilehmm perform very well and both are much better than bac.
from the right plot of figure  <dig>  we can see that the optimal cutoff for tilehmm is close to  <dig> , while it is close to  <dig>  for the bayesian latent method.
figure  <dig> also suggest that both the bayesian latent method and tilehmm are robust to the choice of cutoff values, ranging from  <dig>  to  <dig> , while bac is not.
later, based on the true bound regions which are known for these  <dig> simulate datasets, we use the adjusted rand index r to assess the quality of the results produced by the above three algorithms.
in addition, we calculated p-values of the two-sample t-tests, h0: rbl = ro vs h1: rbl > ro, where rbl denotes the r-value produced by the bayesian latent model, and ro denotes the r-value produced by the other method.
the results were summarized in table  <dig>  the tests indicate that the bayesian latent model can lead to more accurate identifications of true bound regions than bac and tilehmm.
"total" denotes the average number of bound regions identified for each of the  <dig> datasets, nd denotes the number of true bound regions that are not discovered by the algorithm, fd denotes the number of false bound regions discovered by the algorithm, r is the adjusted rand index, the number in the parentheses is the standard error, and "eb t-scan" refers to the empirical bayesian t-scan method proposed by ji and wong  <cit> .
for a thorough comparison, we also applied the sliding window methods, including the wilcoxon rank sum test method  <cit> , t-scan statistic  <cit>  and empirical bayesian t-scan statistic  <cit> , to the  <dig> datasets.
for the testing purpose, we identified the most significant  <dig> probes, which is the same as the true number of bound probes, as the bound probes for each of the datasets and each of the sliding window methods.
we note that this cutoff number should be determined by a multiple hypothesis test in practice, and this choice makes the comparison a little favorly biased toward the sliding window methods.
the results were summarized in the lower panel of table  <dig>  which indicate that the bayesian latent model also outperforms the sliding window methods.
discussion
the bayesian latent model can be generalized in a few ways.
firstly, it can be generalized to allow different bound regions to have different variances.
this generalization has been implemented by us.
the numerical results are very similar to those reported in the paper.
secondly, it can be generalized to work on the multiple replicates directly.
this can be simply done by modifying  to multivariate normals.
this generalization will certainly slow down the simulations, but the results may not be improved significantly.
the reason is that under the assumption of constant variances for probe intensities, the statistic  is sufficient for the mean intensity of probes, while the latter has been designed in the experiment as the main measure for differentiating bound and non-bound regions.
the reason why the bayesian latent method outperforms tilehmm and bac can be explained as follows, through a detailed comparison of the models used by them.
tilehmm implemented a standard two-state hidden markov model, with the emission distribution of state si, i =  <dig>   <dig>  being modeled as a t-distribution.
tilehmm and the bayesian latent model are mainly different in two respects.
• tilehmm is a non-bayesian method, where maximum likelihood estimates are used for all model parameters and inference for the bound regions are based on the conditional probability of the hidden states.
tilehmm is trained using the baum-welch algorithm and the viterbi algorithm.
it is known that the baum-welch algorithm is an em algorithm implemented in the context of hmm, and that it tends to converge to a local maximum of the likelihood function.
the viterbi algorithm provides a fast alternative to the baum-welch algorithm, but may not converge to a local maximum.
the bayesian latent method is a bayesian method, where inference for bound regions is based on the posterior distribution of the latent variable.
the posterior distribution is simulated using the metropolis-within-gibbs sampler, which is known to converge to its target distribution when the number of iterations becomes large.
• tilehmm models all bound regions to have the same mean value, while the bayesian latent model allows different bound regions to have different mean values.
our model fits the real data better.
the mixed performance of tilehmm on the simulated and real datasets indicates that the inferiority of tilehmm is mainly due to its training algorithm.
in addition, as indicated by our simulated examples, tilehmm tends to misidentify the bound regions with relatively low probe intensities, because it models all bound regions to have the same mean value.
bac models the probe intensity using a mixed-effect model:
where c =  <dig> denotes the control sample, c =  <dig> denotes the treatment sample, r is the index of replicates; μp is a random probe effect distributed as ; γcp is the probe enrichment effect with γ1p = 0; and ϵcpr is the random error distributed as .
the authors further modeled the probe enrichment effect by a mixture of a point mass at zero and a truncated gaussian distribution, i.e.,
where tn+() denotes a truncated gaussian distribution truncated at zero, and wp is the a priori proportion of probes belonging to nonbound regions.
the a priori proportion depends on a latent markov random field prior θ = {θp,  <dig> ≤ p ≤ p }, through a logistic transformation
and a gaussian intrinsic autoregressive model  <cit>  for θ,
where ∂p corresponds to the probes p' immediately adjacent to p, np is the cardinality of ∂p, κ is a smoothing parameter, and n is the number of neighboring probes used.
the model is trained using a mcmc algorithm.
the main difference between the bac and the bayesian latent methods is that bac models the control and treatment samples jointly, while the bayesian latent method models the difference between the averaged treatment and control samples.
since bac models the treatment and control samples jointly, it has to include the probe-specific effect in the model and assume a complicated structure for the random error, assuming the variance depends on both the probe and the type of samples .
by working on the difference between the averaged treatment and control samples, the bayesian latent method eliminates the probe effect in the model and the dependence of the random error on the probe and the type of samples.
this simplifies the model greatly and enables an efficient mcmc simulation from the the posterior distribution.
in addition, due to the complicated structure of the model, bac includes too many parameters, and this makes the model potentially overfitted, especially when the number of replicates is small.
this explains why bac always tends to identify too more bound regions than does the bayesian latent model.
on the other hand, the simplicity of the bayesian latent model makes it rather robust to outlying probes.
as indicated by our examples, it work well for all examples studied in this paper.
CONCLUSIONS
we have proposed a bayesian latent model for the chip-chip experiments.
the new model mainly differs from the existing bayesian models, such as the joint deconvolution model, the hierarchical gamma mixture model, and the bayesian hierarchical model, in two respects.
firstly, it works on the difference between the averaged treatment and control samples.
this enables the use of a simple model for the data, which avoids the probe-specific effect and the sample  effect.
as a consequence, this enables an efficient mcmc simulation of the posterior distribution, and also makes the model fairly robust to the outliers.
secondly, it models the neighboring dependence of probes by introducing a latent indicator vector.
a truncated poisson prior distribution is assumed for the latent indicator variable, with the rationale being justified at length.
the bayesian latent model has been successfully applied to the er, p <dig>  and some simulated datasets, with comparisons with bac, tilehmm, and some sliding window methods.
the numerical results indicate that the bayesian latent model can outperform others, especially when the dataset contains outlying probes.
availability and requirements
an r software package called latentchip, which implements the bayesian latent model under linux operating system, is available at the author's web-page  <cit> .
authors' contributions
mw developed the r software package, analyzed the data, and drafted the paper.
fl conceived this work, developed the model, and drafted the paper.
yt participated in the discussion of chip-chip technology.
all authors read and approved the final manuscript.
appendix
the scheme for simulating samples from the posterior distribution:
conditioned on z, updating    using the metropolis-hastings  algorithm, where t indexes the number of iteration cycles.
conditioned on   , updating each component of z according to the following rule: given   : change    to    using the mh algorithm.
when a component of z is updated in step , the sum of square terms in the posterior density function can be calculated in a recursive manner, i.e., only the terms related to zi need to be re-calculated.
supplementary material
additional file 1
convergence study of bayesian latent model.
this file provides a convergence study of bayesian latent model for the er dataset.
click here for file
additional file 2
comparison results for simulated data.
this file provides a comparison study for bac, tilehmm and the latent model on simulated data.
click here for file
acknowledgements
liang's research was supported in part by the grant  made by the national science foundation and the award  made by king abdullah university of science and technology .
tian's research is supported in part by eso <dig>  the authors thank the editor, the associate editor, and the referee for their comments which have led to significant improvement of this paper.

 the chromatin immunoprecipitation coupled with microarray analysis, provides the researchers an efficient way of mapping protein-dna interactions across a whole genome. the chip-chip technology has been used in a wide range of biomedical studies, such as identification of human transcription factor binding sites, investigation of dna methylation, and investigation of histone modifications in animals and plants. data from chip-chip experiments encompass dna-protein interaction measurements on millions of short oligonucleotides which often tile one or several chromosomes or even the whole genome. the data analysis consists of two steps: identifying the bound regions where dna and the protein are cross-linked in the experiments; and identifying the binding sites through sequence analysis of the bound regions. the goal of this paper is to develop an effective method for the first step analysis. analysis of the chip-chip data is very challenging, due to the large amount of probes and the small number of replicates. the existing methods in the literature can be roughly grouped into three categories, the sliding window methods, the hidden markov model methods, and the bayesian methods. other methods have been suggested, e.g., by zheng, huber and reiss, but are less common. the sliding window methods are to test a hypothesis for each probe using the information from the probes within a certain genomic distance sliding window, and then try to correct for the multiple hypothesis tests. the test statistics used are varied. cawley used wilcoxon's rank sum test, keles used a scan statistic which is the average of t-statistics within the sliding window, and ji and wong used a scan statistic which is the average of empirical bayesian t-statistics within the sliding window. since each test uses information from neighboring probes, the tests are not independent, rendering a difficult adjustment in the multiple hypothesis testing step. we note that the power of the sliding window tests is usually low, especially for the tests for the regions where the probe density is low. this is because there will be only very limited neighboring information available for those tests. since, in the chip-chip experiments, the dna samples hybridized to the microarrays are prepared by pcr, which is known to perform independently of the form of dna, the far probes should have similar intensity patterns as long as they are of similar positions to their nearest bound regions. this provides a basis for us to devise powerful methods that make use of information from all probes. the hmm methods have the potential to make use of all probe information, where the model parameters are estimated using all available data. however, in most of the existing implementations of hmms, the model parameters are estimated in an ad hoc way. for example, li estimated the model parameters using previous results on affymetrix snps arrays. an exception is tilehmm, where the model parameters are estimated using the baum-welch and viterbi training algorithms. however, as pointed out by humburg et al., the baum-welch algorithm tends to converge to a local maximum of the likelihood function, and the viterbi training algorithm even fails to converge to a local maximum of the likelihood in some cases. this renders that the parameter estimates and thus followed inference often suboptimal to the problem. bayesian methods have also the potential to make use of all probe information. like the hmm methods, the bayesian methods estimate the model parameters using all available data. however, these methods usually require multiple replicates or some extra experimental information to parameterize the model. for example, the joint binding deconvolution model requires one to know the dna fragment lengths, measured separately for each sample via extrophoretic analysis; and the hierarchical gamma mixture model requires one to first divide the data into genomic regions containing at most one bound region, but such information is, in general, unavailable. the bayesian hierarchical model models the probe intensities using essentially a mixture of normal distributions, and models the spatial structure of the probes using a gaussian intrinsic auto-regression model. gottardo developed a software for the model, bayesian analysis of chip-chip. using bac does not need extra experimental information, but it is extremely slow, roughly hours for a dataset with probes on a personal computer. one reason for the slow speed is the use of mcmc simulations. in this paper, we propose a bayesian latent variable model for tiling array data. our method differs from the existing bayesian methods, such as the joint binding deconvolution model, the hgmm, and the bayesian hierarchical model, in several respects. firstly, it works on the difference between the averaged treatment and control samples. this enables the use of a simple model for the data, which avoids the probe-specific effect and the sample effect. as a consequence, this enables an efficient mcmc simulation of the posterior distribution of the model, and also makes the model rather robust to the outliers. secondly, it models the neighboring dependence of probes by introducing a latent indicator vector. thirdly, it does not require multiple replicates or extra experimental information. as described below, it can work on a single intensity measurement for the probes. the bayesian latent model has been successfully applied to several real and ten simulated datasets, with comparisons with some of the existing bayesian methods, hidden markov model methods, and sliding window methods. the numerical results indicate that the bayesian latent model can outperform the others, especially when the data contain outliers. our method is also computationally efficient; it takes about minutes for a dataset with probes on a personal computer. the remainder of this paper is organized as follows. in section we describe the bayesian latent model and its mcmc implementation. in section we test the bayesian latent model on real and simulated datasets with comparisons with tilehmm, bac and some sliding window methods. in section we discuss possible extensions of our methods and provide an explanation why our method outperforms tilehmm and bac via a detailed comparison of the models used by them. in section we conclude the paper. methods consider a chip-chip experiment with two conditions, treatment and control. let x and x denote, respectively, the samples measured under the treatment and control conditions. each sample has ml, l = replicates providing measurements for n genomic locations along a chromosome or the genome. suppose that these samples have been normalized and log-transformed. in this paper, we summarize the measurements for each probe by where is the intensity measurement of probe i averaged over ml replicates. the underlying assumption for the summary statistic in is that the intensity measurements for each probes has a variance independent of its genomic position. the rationale is that the dna samples used in the experiments are prepared by pcr, which is known to perform independently of the form of dna, and that the amount of the dna samples provides the main sources for the variation of probe intensities. we note that a similar assumption has also been made in other bayesian software, e.g., tilehmm. otherwise, yi can be adjusted by its standard error to a shrinkage t-statistic or an empirical bayes t-statistic, depending on the estimate of the standard error. note that both the adjustments are toward the constant variance of probes. even with the adjustments, the bayesian latent model developed in this paper can still work reasonably well, as the normality assumption approximately holds for the modified t-statistics. the bayesian latent model suppose that the data consists of a total of k bound regions, and that region k consists of nk consecutive probes. for convenience, we call all the non-bound regions by region and denote by n the total number of probes contained in all the non-bound regions, although the probes in which may be non-consecutive. thus, we have. let z = be a latent binary vector associated with the probes, where zi = indicates that probe i belongs to a bound region and otherwise. given z, we can re-index, a realization of, by ykj, k =...,k, j =...,nk. then ykj can be modeled as follows, where μ is the overall mean, which models the difference of sample effects; ν = and νk > k =...,k accounts for the difference of probe intensities in different bound regions; ϵkjs are random errors independently and identically distributed as n. in words, model assumes that, conditioning on the latent vector z, ykjs are mutually independent and also identically distributed within the same bound region. we are aware that for the tiling array data, the probe intensities tend to form a peak around the true binding site. since, given z, the order of probes is meaningless to us, the model is appropriate if ignoring the order of the probes. we note that a similar assumption has also been used in the hgmm and hmm methods. conditioning on z, the likelihood of the model can be written as to conduct a bayesian analysis for the model, we specify the following prior distributions for the model parameters: where ig denotes an inverse gamma distribution, u denotes a uniform distribution, and α, β, νmin, νmax are hyperparameters. in this paper, we set α = β =, which form a vague prior for σ2; and set νmin = 2sy and νmax = maxi yi, where sy is the sample standard error of yi. different values of νmin, e.g., sy and sy, have also been tried in our simulations, and the results are similar. the sensitivity issue of the bayesian latent model to the hyperparameters will be further discussed in section in addition, we assume that the latent vector z follows a truncated poisson distribution, where k, denoting the total number of bound regions specified by z, and is thus a function of z; λ is a hyperparameter; kmax is the largest number of bounded regions allowed by the model; and which makes the prior a proper distribution. the rationale behind this prior can be explained as follows. since the length of each bound region is very short comparing to the chromosome or the whole genome, it is reasonable to view each bound region as a single point, and thus, following the standard theory of poisson process, the total number of bound regions can be modeled as a poisson random variable. conditioning on the total number of bound regions, as implied by, we put an equal prior probability on all possible configurations of z, i.e., assuming a non-informative prior for z. the prior penalizes a large value of k, where the parameter λ represents the strength of penalty. we do not recommend to use a large value of λ, as the number of true bound regions is usually small and a large value of λ will lead to discovery of too many false bound regions. our experience shows that a value of λ around usually works well for the chip-chip data. in this paper, we set λ = in all simulations. the parameter kmax is usually set to a large number. we set kmax = in all simulations of this paper. as long as the value of kmax has been reasonably large, increasing it further would have a negligible effect on simulations. finally, we would like to point out that the bound region identification problem can also be viewed as a change-point identification problem that has been widely studied in statistics. for the change-point identification problem, the same truncated poisson prior has been used for modeling the total number of change-points by many authors, see, e.g., phillips and smith, dension et al., liang and wong, and liang. if ν...,νk ∈, combining the likelihood and prior distributions, integrating out σ and taking the logarithm, we get the following log-posterior density function otherwise, the posterior is equal to due to the design of chip-chip experiments, it is obvious that the intensity measurements of the neighboring probes are positively dependent. to model this dependence, we use a latent indicator vector z. this makes our model different from the existing models, such as the joint binding deconvolution model, the hgmm, and the bayesian hierarchical model. both the joint binding deconvolution model and the bayesian hierarchical model model the mean of probe intensities through the gaussian random field, although their formulations may not be in the standard form of the gmf. like the bayesian latent model, the hgmm models the mean of probe intensities by a piece-wise constant function. the difference is that the hgmm requires one to first divide the data into genomic regions containing at most one bound regions, and thus it allows different non-bound regions to have different means. considering the physical property of pcr, which performs independently of the form of dna, allowing different non-bound regions to have different mean values may not be necessary. mcmc simulations to simulate from the posterior distribution, we used the metropolis-within-gibbs sampler; see appendix for the details. note that when a component of z is updated, the sum of square terms in the posterior density can be calculated in a recursive manner, and this simplifies the computation of the posterior density greatly. inference of bound regions let pi = p be the marginal posterior probability that probe i belongs to a bound region. since the bound regions are expected to consist of several consecutive probes with positive ip-enrichment effects, the regions which consists of several consecutive probes with high marginal posterior probabilities are likely to be bound regions. to identify such regions, we follow gottardo to consider the joint posterior probability where i is the index of the probes, w is a pre-specified half-window size, and m is the minimum number of probes belonging to the bound region. as explained in gottardo, the purpose of introducing the joint posterior probability is to remove the false bound regions, which usually consists of only few isolated probes with large enrichment effects. we found that the choice w = and m = works well in practice. this choice of w is consistent with the moving window size used in other work, such as ji and wong, keles, and gottardo. the choice of m is chosen for robustness to false bound regions. it also reflects our belief that a bound region should consist of at least five consecutive probes with large enrichment effects. note that estimation of ρi is trivial based on the samples simulated from the posterior distribution. the value of ρi depends on a lot of parameters, such as w, m and the hyperparameters of the model. however, we found that the orders of ρi are rather robust to these parameters. this suggests us to treat ρi as a conventional testing p-value, and to control the false discovery rate of the bound regions using a fdr control method, e.g., the empirical bayes method or the stochastic approximation-based empirical bayes method. both the methods allow for the dependence between testing statistics and an empirical determination of the density of the testing statistics. although a strict control of fdr is important to the detection of bound regions, it is not the focus of this paper. in this paper, we will follow other bayesian methods, such as bac, to simply set a cut- off value of ρi. we classify probe i as a probe in bound regions if ρi ≥, and classify probe i as a probe in nonbound region otherwise. as we will see in the numerical examples, the joint posterior probability can lead to a good detection of true bound regions.  the estrogen receptor data the er data were generated by carroll, which mapped the association of the estrogen receptor on chromosomes and here we just used a subset of the data to illustrate how the bayesian latent model works. the subset we used is available from the bac software. it consists of intensity measurements for probes under the treatment and control conditions with three replicates each. the same subset has been used by bac for a demonstration purpose. the bayesian latent model was first applied to the dataset. the algorithm was run times. each run consisted of iterations, and cost about minutes cpu time on a personal computer. all computations of this paper were done on this computer. based on the gelman-rubin diagnostic plot, we discarded the first iterations for the burn-in process, and used the remaining iterations for further inference. figure shows the estimates of the joint posterior probabilities resulted from one run. for comparison, bac and tilehmm were also applied to this dataset. both bac and tilehmm produced a probability measure for each probe, similar to ρi, on how likely it belongs to a bound region. the results were shown in figures and respectively. the comparison shows that all the three methods produced very similar results for this dataset. however, the results produced by the bayesian latent model are neater; the joint posterior probabilities produced by it tend to be dichotomized, either close to or close to this gives the user a clear classification for the bound and non-bound regions. to provide some numerical evidence for this statement, we calculated the ratio #{i: pi > }/#{i: pi > }, where #{i: pi > a} denotes the number of probes with pi greater than a. here pi refers to the joint posterior probability for the bayesian latent model and bac, and the conditional probability for tilehmm. the ratios resultant from the bayesian latent model, bac and tilehmm are, and, respectively. later, we assessed the sensitivity of the bayesian latent method to the values of the hyperparameters νmin and λ with other parameters fixed, α = β = and νmax = maxi yi. the cross settings { sy, sy, sy, 2sy, sy, 3sy} × for were tried for this dataset. for each setting, the algorithm was run times, and each run consisted of iterations. to measure the similarity of the bound regions resultant from different settings of the hyperparameters, we propose to use the adjusted rand index. the adjusted rand index is usually used in the literature of clustering, which measures the degree of agreement between two partitions of the same set of observations even when the comparing partitions having different numbers of clusters. it is obvious that the problem of bound region identification can also be viewed as a clustering problem; where the genome was partitioned into a series of segments, non-bound or bound regions, and each of the segments forms a cluster. the adjusted rand index is defined as follows. let Ω denote a set of n observations, let c = {c...,cs} and represent two partitions of Ω, let nij be the number of observations that are in both cluster ci and cluster, let ni. be the number of observations in cluster ci, and let n.j be the number of observations in cluster. the adjusted rand index is a higher value of r means a higher correspondence between the two partitions. when the two partitions are identical, r is when a partition is random, the expectation of r is under the generalized hyper-geometric model, it can be shown that refer to hubert and arabie for the theoretical justification of r. in calculations of the adjusted rand indices for the sensitivity experiments, we used the result shown in figure as the standard; that is, if a partition is identical to that partition, r will be the results are summarized in figure where the adjusted rand index is plotted as a function of log for different setting of νmin. figure shows that, for each value of νmin, the adjusted rand index varies between and as λ runs from to. this indicates that the performance of the bayesian latent model is rather robust to the choices of νmin and λ. finally, we examined the robustness of the bayesian latent model to different choice of w and m with other parameters fixed at α = β =, λ =, and νmin = 2sy. the cross settings { 10} × { 7} for were tried for this dataset. again, the adjusted rand index is used as the similarity criterion and the result shown in figure as the standard. the results were summarized in table which indicates, for this dataset, the bayesian latent model is quite robust to the choices of w and m. in practice, to achieve robustness to outlying probes, we suggest to avoid choosing a small m. in all the following simulations, we set m = the average of adjusted rand indices is calculated based on independent runs. the entries for the cells with 2w <m are not available. the robustness of the results with respect to changes of α, β and νmax are not studied in the paper. the reason is that νmax is completely determined by the data, and the values of α and β we used form a vague prior for the variance σ p data in a chip-chip experiment, cawley mapped the binding sites of four human transcription factors sp cmyc, p53-fl, and p53-do on chromosomes and the experiment consisted of treatment and input control arrays, and the chromosomes spanned over three chips a, b and c. refer to cawley for the details of the experiment. for the testing purpose, p53-fl data on chips a, b and c were used in this paper, which contains quantitative pcr verified regions. as in cawley, the data were pre-processed by filtering out the local repeats, quantile-normalized, rescaled to have a median feature intensity of for the purpose of adjusting batch effect, and then log-transformed. since the normalization is not the focus of this paper, we skipped the details. the bayesian latent method was first applied to the p data. the data on chip a, chip b, and chip c were analyzed separately. each run consisted of iterations. diagnostic plot for the convergence of these runs indicates that they can converge within several hundreds of iterations, even the data on each chip consists of more than probes. accordingly, the first iterations were discarded for the burn-in process, and the samples from other iterations are used for further analysis. for comparison, bac and tilehmm were also applied to this example. given the posterior probabilities, a cutoff of was used for all methods to detect bound regions. all resultant bound regions having less than probes or bps were considered to be spurious and removed, and those regions separated by bps or less were merged together to form a predicted bound regions following the approach taken by cawley. the results were summarized in table although tilehmm detected all the validated regions, it essentially fails for this example. it identified a total of bound regions, which should contain too many false bound regions. we suspect that the failure of tilehmm for this example is due to its training algorithm; it is very likely that tilehmm converged to a local maximum of the likelihood function. this have been noted by humburg et al., tilehmm may converge to a local maximum of the likelihood function with either the baum-welch algorithm or the viterbi training algorithm, rendering an ineffective inference for the model. both the total numbers of regions and quantitative pcr verified ones detected by each method on each chip are reported. the columns under "p53" summarize the results on chips a, b and c. the number in the parentheses is the number of clusters needed to cover all experimentally validated bound regions. both the bayesian latent method and bac work well for this example. at a cutoff of, bac identified bound regions, which cover out of experimentally validated bound regions. the bayesian latent method works even better. at the same cutoff, it only identified bound regions, but which also cover out of experimentally validated bound regions. for further comparison of the bayesian latent method and bac, we relaxed the cutoff value and counted the total number of regions needed to cover all experimentally validated regions. we found that the bayesian latent method only needs to increase the total number of regions to while bac needs to increase to regions. note that the bac and tilehmm's results reported here may be a little different from those reported by other authors, due to the difference of the normalization methods. simulated data to have a careful assessment of the performance of the bayesian latent model, we simulated datasets based on the sp data of cawley's experiment. each dataset consists of probes, two conditions, and six replicates under each condition. the probe genomic coordinates we used in simulations were the first genomic positions used in the sp data. each dataset consisted of bound probes, forming bound regions. as in gottardo, the bound regions were assumed to describe a peak with the intensity function given by a exp{-42/b2}, where a is the amplitude of the peak, b controls the width of the peak, c represents the center of the peak, and gi is the genomic position of probe i. we also followed gottard to generate the centers of the bound regions randomly across the set of possible coordinates while imposing a separation of at least bps between peaks; and to generate the values of parameter b uniformly between and bps. the values of parameter a were generated uniformly between and the variance of the probe intensity was estimated from the sp data. firstly, the performance of different models is assessed using the area under the receiving operating characteristic curve and the error rate. the former is a standard measure for the performance of a multiple hypothesis testing method, which shows the true positive discovery rate against the false positive discovery rate at probe level. the later is a standard measure for the performance of a classification method, which shows the proportion of totally incorrect probe calls, including both false positives and false negatives, against different cutoff values. all the three methods, bayesian latent method, bac and tilehmm, were applied to the full datasets. the averaged roc curve and error rate across a range of cutoffs are obtained and plotted in figure as indicated by figure the bayesian latent method and tilehmm have very similar performances on these datasets, and both are much better than bac. by further examining the plot on the right, which provides a closer view of the area enclosed by the dotted line and axis on the left, it is easy to see that the bayesian latent method is better than tilehmm for this example. next, we checked the error rate for each model. the results were shown in figure again, the bayesian latent method and tilehmm perform very well and both are much better than bac. from the right plot of figure we can see that the optimal cutoff for tilehmm is close to, while it is close to for the bayesian latent method. figure also suggest that both the bayesian latent method and tilehmm are robust to the choice of cutoff values, ranging from to, while bac is not. later, based on the true bound regions which are known for these simulate datasets, we use the adjusted rand index r to assess the quality of the results produced by the above three algorithms. in addition, we calculated p-values of the two-sample t-tests, h0: rbl = ro vs h1: rbl > ro, where rbl denotes the r-value produced by the bayesian latent model, and ro denotes the r-value produced by the other method. the results were summarized in table the tests indicate that the bayesian latent model can lead to more accurate identifications of true bound regions than bac and tilehmm. "total" denotes the average number of bound regions identified for each of the datasets, nd denotes the number of true bound regions that are not discovered by the algorithm, fd denotes the number of false bound regions discovered by the algorithm, r is the adjusted rand index, the number in the parentheses is the standard error, and "eb t-scan" refers to the empirical bayesian t-scan method proposed by ji and wong. for a thorough comparison, we also applied the sliding window methods, including the wilcoxon rank sum test method, t-scan statistic and empirical bayesian t-scan statistic, to the datasets. for the testing purpose, we identified the most significant probes, which is the same as the true number of bound probes, as the bound probes for each of the datasets and each of the sliding window methods. we note that this cutoff number should be determined by a multiple hypothesis test in practice, and this choice makes the comparison a little favorly biased toward the sliding window methods. the results were summarized in the lower panel of table which indicate that the bayesian latent model also outperforms the sliding window methods. discussion the bayesian latent model can be generalized in a few ways. firstly, it can be generalized to allow different bound regions to have different variances. this generalization has been implemented by us. the numerical results are very similar to those reported in the paper. secondly, it can be generalized to work on the multiple replicates directly. this can be simply done by modifying to multivariate normals. this generalization will certainly slow down the simulations, but the results may not be improved significantly. the reason is that under the assumption of constant variances for probe intensities, the statistic is sufficient for the mean intensity of probes, while the latter has been designed in the experiment as the main measure for differentiating bound and non-bound regions. the reason why the bayesian latent method outperforms tilehmm and bac can be explained as follows, through a detailed comparison of the models used by them. tilehmm implemented a standard two-state hidden markov model, with the emission distribution of state si, i = being modeled as a t-distribution. tilehmm and the bayesian latent model are mainly different in two respects. • tilehmm is a non-bayesian method, where maximum likelihood estimates are used for all model parameters and inference for the bound regions are based on the conditional probability of the hidden states. tilehmm is trained using the baum-welch algorithm and the viterbi algorithm. it is known that the baum-welch algorithm is an em algorithm implemented in the context of hmm, and that it tends to converge to a local maximum of the likelihood function. the viterbi algorithm provides a fast alternative to the baum-welch algorithm, but may not converge to a local maximum. the bayesian latent method is a bayesian method, where inference for bound regions is based on the posterior distribution of the latent variable. the posterior distribution is simulated using the metropolis-within-gibbs sampler, which is known to converge to its target distribution when the number of iterations becomes large. • tilehmm models all bound regions to have the same mean value, while the bayesian latent model allows different bound regions to have different mean values. our model fits the real data better. the mixed performance of tilehmm on the simulated and real datasets indicates that the inferiority of tilehmm is mainly due to its training algorithm. in addition, as indicated by our simulated examples, tilehmm tends to misidentify the bound regions with relatively low probe intensities, because it models all bound regions to have the same mean value. bac models the probe intensity using a mixed-effect model: where c = denotes the control sample, c = denotes the treatment sample, r is the index of replicates; μp is a random probe effect distributed as; γcp is the probe enrichment effect with γ1p = 0; and ϵcpr is the random error distributed as. the authors further modeled the probe enrichment effect by a mixture of a point mass at zero and a truncated gaussian distribution, i.e., where tn+() denotes a truncated gaussian distribution truncated at zero, and wp is the a priori proportion of probes belonging to nonbound regions. the a priori proportion depends on a latent markov random field prior θ = {θp, ≤ p ≤ p }, through a logistic transformation and a gaussian intrinsic autoregressive model for θ, where ∂p corresponds to the probes p' immediately adjacent to p, np is the cardinality of ∂p, κ is a smoothing parameter, and n is the number of neighboring probes used. the model is trained using a mcmc algorithm. the main difference between the bac and the bayesian latent methods is that bac models the control and treatment samples jointly, while the bayesian latent method models the difference between the averaged treatment and control samples. since bac models the treatment and control samples jointly, it has to include the probe-specific effect in the model and assume a complicated structure for the random error, assuming the variance depends on both the probe and the type of samples. by working on the difference between the averaged treatment and control samples, the bayesian latent method eliminates the probe effect in the model and the dependence of the random error on the probe and the type of samples. this simplifies the model greatly and enables an efficient mcmc simulation from the the posterior distribution. in addition, due to the complicated structure of the model, bac includes too many parameters, and this makes the model potentially overfitted, especially when the number of replicates is small. this explains why bac always tends to identify too more bound regions than does the bayesian latent model. on the other hand, the simplicity of the bayesian latent model makes it rather robust to outlying probes. as indicated by our examples, it work well for all examples studied in this paper.  we have proposed a bayesian latent model for the chip-chip experiments. the new model mainly differs from the existing bayesian models, such as the joint deconvolution model, the hierarchical gamma mixture model, and the bayesian hierarchical model, in two respects. firstly, it works on the difference between the averaged treatment and control samples. this enables the use of a simple model for the data, which avoids the probe-specific effect and the sample effect. as a consequence, this enables an efficient mcmc simulation of the posterior distribution, and also makes the model fairly robust to the outliers. secondly, it models the neighboring dependence of probes by introducing a latent indicator vector. a truncated poisson prior distribution is assumed for the latent indicator variable, with the rationale being justified at length. the bayesian latent model has been successfully applied to the er, p and some simulated datasets, with comparisons with bac, tilehmm, and some sliding window methods. the numerical results indicate that the bayesian latent model can outperform others, especially when the dataset contains outlying probes. availability and requirements an r software package called latentchip, which implements the bayesian latent model under linux operating system, is available at the author's web-page. authors' contributions mw developed the r software package, analyzed the data, and drafted the paper. fl conceived this work, developed the model, and drafted the paper. yt participated in the discussion of chip-chip technology. all authors read and approved the final manuscript. appendix the scheme for simulating samples from the posterior distribution: conditioned on z, updating using the metropolis-hastings algorithm, where t indexes the number of iteration cycles. conditioned on, updating each component of z according to the following rule: given: change to using the mh algorithm. when a component of z is updated in step, the sum of square terms in the posterior density function can be calculated in a recursive manner, i.e., only the terms related to zi need to be re-calculated. supplementary material additional file 1 convergence study of bayesian latent model. this file provides a convergence study of bayesian latent model for the er dataset. click here for file additional file 2 comparison results for simulated data. this file provides a comparison study for bac, tilehmm and the latent model on simulated data. click here for file acknowledgements liang's research was supported in part by the grant made by the national science foundation and the award made by king abdullah university of science and technology. tian's research is supported in part by eso the authors thank the editor, the associate editor, and the referee for their comments which have led to significant improvement of this paper.
the new model mainly differs from the existing bayesian models, such as the joint deconvolution model, the hierarchical gamma mixture model, and the bayesian hierarchical model, in two respects. they usually require multiple replicates or some extra experimental information to parametrize the model, and long cpu time due to involving of mcmc simulations. the numerical results indicate that the bayesian latent method can outperform other methods, especially when the data contain outliers. firstly, it works on the difference between the averaged treatment and control samples. the chip-chip technology has been used in a wide range of biomedical studies, such as identification of human transcription factor binding sites, investigation of dna methylation, and investigation of histone modifications in animals and plants. in this paper, we propose a bayesian latent model for the chip-chip data. various methods have been proposed in the literature for analyzing the chip-chip data, such as the sliding window methods, the hidden markov model-based methods, and bayesian methods. the bayesian latent method is successfully applied to real and ten simulated datasets, with comparisons with some of the existing bayesian methods, hidden markov model methods, and sliding window methods. secondly, it models the neighboring dependence of probes by introducing a latent indicator vector. a truncated poisson prior distribution is assumed for the latent indicator variable, with the rationale being justified at length. as a consequence, this enables an efficient mcmc simulation of the posterior distribution of the model, and also makes the model more robust to the outliers. this enables the use of a simple model for the data, which avoids the probe-specific effect and the sample effect. although, due to the integrated consideration of uncertainty of the models and model parameters, bayesian methods can potentially work better than the other two classes of methods, the existing bayesian methods do not perform satisfactorily.
 the chip-chip technology has been used in a wide range of biomedical studies, such as identification of human transcription factor binding sites, investigation of dna methylation, and investigation of histone modifications in animals and plants. various methods have been proposed in the literature for analyzing the chip-chip data, such as the sliding window methods, the hidden markov model-based methods, and bayesian methods. although, due to the integrated consideration of uncertainty of the models and model parameters, bayesian methods can potentially work better than the other two classes of methods, the existing bayesian methods do not perform satisfactorily. they usually require multiple replicates or some extra experimental information to parametrize the model, and long cpu time due to involving of mcmc simulations.  in this paper, we propose a bayesian latent model for the chip-chip data. the new model mainly differs from the existing bayesian models, such as the joint deconvolution model, the hierarchical gamma mixture model, and the bayesian hierarchical model, in two respects. firstly, it works on the difference between the averaged treatment and control samples. this enables the use of a simple model for the data, which avoids the probe-specific effect and the sample effect. as a consequence, this enables an efficient mcmc simulation of the posterior distribution of the model, and also makes the model more robust to the outliers. secondly, it models the neighboring dependence of probes by introducing a latent indicator vector. a truncated poisson prior distribution is assumed for the latent indicator variable, with the rationale being justified at length.  the bayesian latent method is successfully applied to real and ten simulated datasets, with comparisons with some of the existing bayesian methods, hidden markov model methods, and sliding window methods. the numerical results indicate that the bayesian latent method can outperform other methods, especially when the data contain outliers. 