Fast MCMC sampling for hidden markov models to determine copy number variations

Background Hidden Markov Models (HMM) are often used for analyzing Comparative Genomic Hybridization (CGH) data to identify chromosomal aberrations or copy number variations by segmenting observation sequences. For efficiency reasons the parameters of a HMM are often estimated with maximum likelihood and a segmentation is obtained with the Viterbi algorithm. This introduces considerable uncertainty in the segmentation, which can be avoided with Bayesian approaches integrating out parameters using Markov Chain Monte Carlo (MCMC) sampling. While the advantages of Bayesian approaches have been clearly demonstrated, the likelihood based approaches are still preferred in practice for their lower running times; datasets coming from high-density arrays and next generation sequencing amplify these problems. Results We propose an approximate sampling technique, inspired by compression of discrete sequences in HMM computations and by kd-trees to leverage spatial relations between data points in typical data sets, to speed up the MCMC sampling. Conclusions We test our approximate sampling method on simulated and biological ArrayCGH datasets and high-density SNP arrays, and demonstrate a speed-up of 10 to 60 respectively 90 while achieving competitive results with the state-of-the art Bayesian approaches. Availability: An implementation of our method will be made available as part of the open source GHMM library from http://ghmm.org.


Background
The Sirens' call of Bayesian methods is that we can do without the parameters of models and, instead, compute probabilities of interest directly, indicating for example how likely a biological fact is given our data. The price one pays for that convenience is on one hand the conundrum of which prior distributions to choose and how to set their hyper-parameters; the frequent use of uniform priors is evidence that this is indeed non-trivial. On the other hand, unless the choice of likelihood and prior yields simple posteriors which we can integrate symbolically, we have to resort to sampling for example with Markov Chain Monte Carlo (MCMC) [1].
In the following we will concentrate on HMMs, stochastic functions of Markov Chains, which have not only been used extensively for discrete sequences-pairwise-sequence alignments with pair-HMMs [2], gene finding with labeled HMMs [3], and detecting remote homologs using profile HMMs [4]-but also for continuous-valued observations, such as gene expression time-courses [5]. Continuous observation sequences from either DNA microarrays or next generation sequencing experiments, note that the proportion of mapped reads in an interval is frequently used as a continuous measure of copy number, to detect chromosomal aberrations or copy number variations lead to the same fundamental computational problem and share characteristics of the data. The goal is to segment an observation sequence into regions in which there is little variation around a common mean. In other words, the assumption is that the data can be approximately described by piece-wise constant functions. Indeed, if hybridization intensity was an exact, un-biased measurement of DNA concentration before amplification, the sequence of hybridization intensities of probes along a chromosome would yield a piece-wise constant function in ArrayCGH experiments. This assumption holds true for a mixture of different cell populations because a finite sum of piece-wise constant functions is again a piece-wise constant function.
A wide range of methods for copy number detection in ArrayCGH data have been developed in recent years, including change-point detection based methods [6,7], smoothing based methods [8,9], and hierarchical clustering [10]. Here, we concentrate on HMM-based approaches which have been proposed for segmenting sequences of continuous-valued observations and shown to match or improve upon the state-of-the-art [11][12][13]. Once a model is trained from the data, either using maximum likelihood (ML) or maximum a posteriori (MAP), the segmentation is given by the most likely state sequence obtained with the Viterbi algorithm [14]. The segmentation, however, is highly dependent on the model parameters. A small change in the computed parameter values can adversely affect the recovered segmentation. A full Bayesian approach alleviates this dependence by integrating out the model parameters. As analytic integration of a complex high dimensional model is impossible for most distributions, the Bayesian approach requires the use of numerical integration techniques like MCMC [15], for example by direct Gibbs sampling [16] of model parameters and state paths. Scott [17] reports that combining the forward-backward recursions [18] and Gibbs sampling improves the converge rate and consequently the running time. Nevertheless, MCMC remains substantially slower than training one model and running Viterbi once and the loss in reliability introduced by relying on one ML or MAP model is ignored in practice. For discrete emissions, compressing sequences and computing forward and backward variables and Viterbi paths on the compressed sequences yields impressive speed-ups [19]. However, discretization of continuous emissions, similar to vector quantization used in speech recognition [18], is not viable as the separation between the different classes of observations is low since the observations are one-dimensional. Moreover, maximal compression is to be expected for small number of discrete symbols and clearly compression ratio conflicts with fidelity in the analysis.
For a different task, arguments about spatial relations between groups of multi-variate data points were used to achieve considerable speed-up. Moore and colleagues used modified kd-trees, a data structure to efficiently execute spatial queries such as determining the nearest neighbor of a given point, to accelerate k-means [20]. The fundamental insight is illustrated in Figure 1 (left). In the reassignment step of k-means one has to find the nearest centroid for every data point. Due to the kd-tree, there are groups of points contained in a node of the tree for which this decision about the nearest centroid can be made simultaneously by a geometrical argument about the vertices of the hyperrectangle defined by this node. A similar kd-tree based approach was used in speech recognition [21,22] to quickly find the most important components in a mixture of large number of Gaussians and thus approximate the full observation density in one individual HMM state with multi-variate emissions.
At the core of our approach is a similar geometrical argument about several uni-variate data points based on a modified kd-tree. We adaptively identify blocks of observations, cf. Figure 1 (middle). For all observations in a block we now estimate, at least conceptually, the most likely state simultaneously depending on the means of the Gaussians in each state to gain a considerable speed-up proportional to the average block length. Similarly, we can avoid sampling states for each individual observation in a block if we can bound the posterior. Considerable care has to be taken for combining Figure 1 Fundamental insight. When reassigning two-dimensional data points to cluster centroids c and d in k-means clustering (left) the hyperrectangles obtained from kd-trees reduce the computational effort by making an argument about all points in an hyperrectangle based on their vertices; consider for example the rightmost hyperrectangle. For sequences (middle) there is no overlap in y-direction and decisions about the most likely state can be made per block considering the means of the Gaussians of a three-state HMM (right), μ -, μ = and μ + . Note that at every given block only a decision between the two states with closest mean is necessary, if one assume comparable variances. Decision boundaries are displayed dashed. blocks and to bound the errors introduced by the approximations based on geometric arguments.
In summary, we • propose the first use of spatial index structures for several consecutive observations and approximate computations based on geometric arguments to substantially speed-up the problem of segmenting sequences of continuous observations using HMM, • demonstrate that very frequently used biological datasets of high relevance measuring chromosomal aberration and copy number variations are consistent with our assumptions of piece-wise constant observations with added noise, and • achieve speed-ups between 10 and 90 on those biological datasets while maintaining competitive performance compared to the state-of-the-art.

HMM
We only consider HMMs with Gaussian emission distributions; see Figure 1 (right) for an example and [18] for an introduction. We will use the following notation: N denotes the number of states, S = {s 1 , s 2 , ..., s N } the set of states, T the length of an observation ,j≤N the transition matrix, π = (π 1 , π 2 , ..., π N ) the initial distribution over states, B = μ 1 , σ 2 1 , . . . , μ N , σ 2 N with μ 1 ≤ ... ≤ μ N are parameters of the emission distributions, and Q = {q 1 , q 2 , ..., q T } the hidden state sequences with q t S. From an observation sequence O we can obtain a maximum likelihood point estimate of the parameters (A, B, π) using the Expectation-Maximization (EM) or Baum-Welch [23] algorithm and compute the most likely path using the Viterbi algorithm [14].

MCMC Sampling
Bayesian analysis requires choosing prior distributions on parameters and we use standard conjugate prior distributions following [1]. Specifically, we choose and π~Dirichlet(λ π ), whereμ i ,σ i , a i , b i , λ A i , and λ π are the hyperparameters of the model.
It is only possible in some trivial cases to compute posterior distribution in closed form using analytic integration. In most applications, specially for high dimensional distributions, Monte Carlo integration techniques, like MCMC sampling by Gibbs sampling or Metropolis-Hastings, are easier to compute and generally produce more accurate results than numerical integration [15]. Scott [17] compares various MCMC approaches and strongly argues in favor of forward-backward Gibbs sampling (FBG-sampling), which has been successfully used by others [24,25], for it's excellent convergence characteristics. We briefly summarize FBG-sampling for an HMM λ = (A, B, π); see [17,26] for details: 1. Choose initial parameters θ 0 = (A 0 , B 0 , π 0 ) 2. Perform the following steps for iteration 0 ≤ m <M.
(b) Sample q m T ∼ P(q T , O|θ m ) . (c) Sample the state sequence Q m in a backward fashion for T >t ≥ 1.
Despite its fast convergence, FBG-sampling takes O (MTN 2 ) time for M iterations. For long observation sequences with millions of observations, as they are common in biological applications, realistic values for M and N make FGB-sampling intractable. In the next section we discuss how FBG-sampling can be approximated to improve the running time to O(γMTN 2 ), where γ is the compression ratio of the observation sequence, while maintaining good statistical properties. We refer to our sampling technique as approximate sampling.

Approximate sampling
Through application of a modified kd-tree algorithm (details below), we compress the observation sequence Figure 1 (middle), and store precomputed first moment, second moment, and the number of observations compressed into block o t for all i ≤ T'. In subsequent MCMC iterations we assume that observations compressed in a block o t arise from the same underlying state. In other words we ignore the contribution of the state paths that do not go through the same state for o t . By ignoring those state paths, we refer to them as weak state paths, when computing forward-variables and by reusing pre-computed statistics we are able to accelerate sampling.
At first ignoring weak state paths may sound like a very crude approximation for computing forward variables. But in many applications this is certainly not true. We demonstrate with a symmetric Gaussian HMM that the weak state path assumption is a fairly realistic approximation and leads to faster sampling. We define a symmetric HMM λ = (A, B, π) with N states s 1 , ..., s N , where we set self-transition probability a ii = t and non-self-transition Remark 1 For realistic values of τ, t, and n, the contribution from ignored weak state paths, which we call , can be very small. If ≪ 1, ignoring weak state paths will not introduce large errors in the computation. For the 2-state example in Section: Results, where t = 0.9, d = 1, and s 2 = 0.1, is at most 1 3 for a block length n ≤ 10 if we assume τ > 0.25 and a = 1. If τ is much larger and consequently c 2n N is much smaller, we can roughly say that n can be as large as 1 + log 1+rc (1 + ) in a symmetric Gaussian HMM.
Remark 2 We often encounter situations where P(q i = s x |O i-1 ) ≫ P(q i ≠ s x |O i-1 ). Even though it is not exploited in the lemma (a being greater than or equal to 1), as a consequence of this, the observation sequence can be compressed into larger blocks keeping small in practice.
From the above lemma and the remarks we see that the longer blocks created by an algorithm should concentrate around the means of the Gaussian distributions. While a brute force algorithm looks at local information, a kd-tree like algorithm alternately looks at both dimensions and utilizes global information like the density of data points (around the means data concentration should be high) to create better quality blocks. We use a modified kd-tree based algorithm to find such blocks and discuss the details below.

kd-tree Based Sequence Compression
Given a starting width parameter w we create a list of nodes from the observation sequence O = o 1 , ..., o T using the following steps.

If | max
w δ L or |O| = 1, create a node storing the first and second moment of the observations in O, append it to O', and then go to the end step. Otherwise, go to the next step.
and O r , go to step 2 keeping the level value L unchanged, and setting d = 1. 5. End step.
In the above recursive algorithm, w states the initial width, δ controls the rate of width shrinking in successive levels of the iterations, and O' accumulates the compressed blocks of observations. The current iteration level L, the current dimension d, and the current input set O are local variables in the recursive algorithm. Notice that we start with an empty list O' and at the end of the recursive procedure O' contains an ordered list of compressed observations. To gain further compression of the sequence, we sequentially go through the blocks of O' and combine consecutive blocks if the distance between their means is less than w. We also combine three consecutive blocks if the outer blocks satisfy this condition and the inner block has only one observation. In step 3 of the above algorithm, the input set is divided into two subsets and each subset contains half of the elements from the original set. Consequently, the height of the recursion tree is at most 2logT and the running time of the above algorithm is O(T log T). This overhead is negligible compared to the time that it takes to run M iterations of MCMC sampling.

Width Parameter Selection
For increasing values of w the average block size increases exponentially in the above kd-tree based compression. As a result, the compression ratio γ = T T plotted as a function of w, has a knee which can inform the choice of w. Moreover, methods originally developed to find the optimal numbers of clusters in clustering can be used to find the knee of such a curve automatically. In particular, we use the L-method [27] which finds the knee as the intersection of two straight lines fitted to the compression curve.

Fast Approximate Sampling
Given the compressed input sequence O = o 1 , o 2 , . . . , o T computing forward variables and subsequent sampling is a straightforward modification of the uncompressed case. In particular we make the following two changes to the FBG-sampling algorithm.

Modified forward variables at positions
are computed using the following formula, constant time computation using precomputed statistics 2. After sampling the state sequence, parameters are updated ignoring non-self transitions between consecutive observations in o t .
Clearly, each iteration of approximate sampling takes

Results
We evaluate FBG-sampling and approximate sampling in three different settings. First, its effectiveness is verified for a simple two state model. Then, we test on simulated ArrayCGH data which is the accepted standard for method evaluation [28]. Finally, we report findings from an analysis of Mantle Cell Lymphoma (MCL) cell lines [29], Corriel cell lines [30], GBM datasets [31], and high resolution SNP arrays [13,32]. For biological data, if multiple chromosomes are present, we use pooling [25] across chromosomes, which does not allow transition between different chromosomes but assumes model parameters to be identical across chromosomes. Throughout this section we define s D to be the standard deviation of all observations in the dataset. We compress the dataset with increasing values of w = 0.25s D , 0.5s D , 0.75s D , .... For evaluation we consider the experiments as two class problems: aberrant clones belong to the positive class and normal clones belong to the negative class. When ground truth labels of a dataset are available we report F1-measure, recall, and precision for the experiment. With tp, fp, tn, fn we denote the number of true and false positives and true and false nega- 2 . From λ 2ST we sample an observation sequence O = o 1 , ..., o 10,000 , and run MCMC for M = 100 steps with hyperparameter valuesμ 1:2 = 0, 1 for the prior mean on μ, σ 1:2 = 0.5, 0.5 for the prior variance on μ, a 1:2 = 4, 4 for the shape of Gamma prior on σ -2 , b 1:2 = 1, 1 for the rate of Gamma prior on s -2 , δ π = 1, 1 for the Dirichlet prior on the initial distribution π, and δ A i 1:2 = 1, 1 for the Dirichlet prior on row i of transition matrix A.
After M iterations, we compare the posterior prob- FBG and θ M A are M-th parameter samples of FBGsampling and approximate sampling. Figure 2 shows that the posterior probability of being in state 1 for each position can be approximated fairly well even for large values of w. The average posterior error Table 1. Similarly, we compute the Viterbi paths and report total number of mismatches between them along with the likelihoods in Table 1.

Simulation from Genetic Template
We use 500 simulated datasets published in [28]. Each dataset has 20 chromosomes and 100 clones per chromosome for a total of 2,000 clones per dataset. A fourstate HMM predicts the aberrant regions-loss defined as state S 1 and gain defined as state S 3 or S 4 . The neutral region is modeled as state S 2 . We put an ordering constraint on the means, μ 1 <μ 2 <μ 3 <μ 4 , to prevent label switching of the states [17]. Hyperparameter choices follow [25] and areμ 1:4 = −0.5, 0, 0.58, 1 for the prior mean on μ,σ 1:4 = 0.5, 0.001, 1.0, 1.0 for the prior variance on μ, a 1:4 = 10,100, 5, 5 for the shape of gamma prior on s -2 , and b 1:4 = δ π = δ A i 1:4 = 1, 1, 1, 1 for the rate of gamma prior on σ -2 , the Dirichlet prior on initial distribution π, and the Dirichlet prior on row i of transition matrix A, respectively. Table 2 shows the mean and standard deviation of F1measure, recall, and precision over the 500 datasets for FBG-sampling, approximate sampling, and Expectation Maximization (EM) with the ground truth provided by [28]. Even for this collection of relatively small datasets we see a 10-fold speed up. For each dataset we run FBG and approximate sampling for M = 100 steps (we have visually monitored the parameters and noticed convergence within 50 steps, see Figure 3 for a representative example). The last 10 samples are used to compute 10 samples of the posteriors for each state and for each position in the observation sequence. Subsequently, aberrant regions are predicted based on the average of those distributions. We report the speed-up of approximate vs. FBG sampling based on the time it takes to compress the sequence and run M steps of MCMC. For one individual dataset EM requires 58 seconds on average, which allows for a total of 800-1000 repetitions from randomized points sampled from the prior distributions in the time needed for FBG sampling. Each run continues until the likelihood converges and the best model based on likelihood is selected. Aberrant regions are predicted and compared against the ground truth based on the Viterbi path. We report the mean and standard deviation of F1-measure, recall, and precision over the results of EM on 500 datasets.

Biological Data Mantle Cell Lymphoma (MCL)
De Leeuw and colleagues identified recurrent variations across cell lines using ArrayCGH data of MCL cell lines [29]. Out of the eight cell lines [29] HBL-2 was fully annotated with marked gain and loss regions in the autosomes. This dataset contains about 30,000 data points (combining all the autosomes). We have used a four-state HMM for predicting aberrant regions. State 1 represents copy number loss, state 2 represents normal copy number, state 3 represents single copy gain, and state 4 multiple gain. For HBL-2 we report the F1-measure, recall, precision and speed-up. Similar to the synthetic case we put an ordering constraint on the means, μ 1 <μ 2 <μ 3 <μ 4 . Hyperparameter choices follow [25] and  for each simulated dataset sampling methods run once and we report the average and standard deviation over 500 datasets, but for HBL-2 we let them run 10 times and report the average and standard deviation of these 10 F1-measures, recalls, and precisions in Table  2. Each EM run starts with the initial parameter values sampled either from the prior distributions, or from uniform distributions, and continues until the likelihood value converges. We report the performance of the most likely model (which is the preferred criteria to select a model), the likelihood of the best model based on F1-measure, and the average and standard deviation of F1-measures, recalls, and precisions of all the models generated by EM.
As representative examples, we also plot the segmentation of chromosome 1 and 9 computed by FBG-sampling and approximate sampling along with the ground truth labels in Figure 4.

Corriel
Corriel cell lines were used by Snijders et al. [30] and are widely considered a gold standard in ArrayCGH data analysis. This dataset is smaller and, in fact, fairly easy compared to the MCL cell lines. For the Corriel cell line we use a 4-state HMM and report the results for GM05296 and GM00143 in Table 2. Again, approximate sampling performs competitively while achieving more than a 10-fold speed-up. Hyperparameter choices follow [24].

GBM
The glioma data from Bredel et al. [31] has previously been used to analyze the performance of CNV detection methods [9,33]. According to [33], GBM datasets are noisy Approximate sampling results are reported for different choices of w. The w value which is closest to the one estimated by the L-method is shown in italic. Width w is reported in σ D of the corresponding dataset, time is reported in seconds, and compression is defined as T T . For HBL-2, the initial parameter values for EM are sampled from the prior or uniform distributions, and the average (mean), most likely (ML), and best (in terms of F1-measure) performances along with likelihoods are reported.
but contains a mixture of aberrant regions with different width and amplitude. In particular, chromosome 13 of GBM31 is reported to have low amplitude loss in p-arm and chromosome 7 of GBM29 is reported to have high amplitude gains near the EGFR locus by previous studies [9,33]. The segmentation of these two chromosomes are shown in Figure 5. Although [33] reports that EM based HMM failed to detect these aberrations we see that Bayesian HMM has successfully detected both the gain in show the outcomes of approximate sampling. Note that the reduction in standard deviation suggests that approximate sampling converges quicker than FBG sampling for these probes. and the Dirichlet prior on row i of transition matrix A, respectively, and at the recommended w value we see a 10 fold speed-up.

SNP Array
High-resolution Single Nucleotide Polymorphism (SNP) arrays are capable of detecting smaller CNVs than ArrayCGH. To demonstrate the computational advantage of approximate sampling on SNP arrays we have chosen publicly available Affymetrix 100 k pancreatic cancer datasets from [32] and Illumina HumanHap550 arrays of HapMap individuals which are provided as examples in PennCNV [13]. An Affymetrix 100 k dataset consists of two different arrays each with ≈ 60, 000 SNP markers and, in total, 10 5 data points per sample. On the other hand, the Illumina HumanHap550 array has around half a million SNP markers. We have applied FBG-sampling and approximate sampling with w = 1.8s D , the recommended value by the L-method, to the sample datasets from Harada et al. [32] and found that the computational speed-up is 22-fold (100 runs of FBG-sampling takes 1844 seconds). Both sampling approaches mostly agree on their predictions but they, specially FBG-sampling, detect several more CNVs than previously identified [32]. For example, the high amplification in chromosome 11 (sample 33) is successfully identified by all methods but in chromosome 18 (sample 16) the sampling algorithms find a few normal regions previously predicted [32] as loss using the CNAG tool [34] (see Figure 6). One possible reason for these differences is that while we use 269 HapMap samples as reference they use 12 unpublished normal samples as reference. Similarly, we have tested our method with 2.0s D ≤ w ≤ 3.0s D against Illumina HumanHap samples and observed 70 to 90 fold speed-up in computations (100 runs of FBG-sampling takes 9693 seconds). These samples are taken from apparently healthy individuals and contain very few CNVs. As expected, both sampling algorithms' predictions are nearly identical and they seem to predict extreme outliers as aberrant markers. In contrast, PennCNV [13] does not report CNVs which are covered by less than 3 SNPs, thus suppressing the outliers as normal. We plot a typical region (from 1.4e + 08bp to 1.7e + 08bp) of chromosome 6 from sample 3 (ID 99HI0700A) in Figure 7.
To set hyperparameters we follow the default parameters of the HMM used in PennCNV [13]. We have observed that HMMs for large arrays are particularly sensitive to the self-transition probabilities (which is also reflected in the default parameter values of the HMM used in PennCNV). Hence, hyperparameters were set to Other hyperparameters for the 3-state HMM were set such that the expected values of prior distributions match the default values for PennCNV. In particular, they wereμ 1:3 = 10.66, 0, 0.54 for the prior mean on μ, σ 1:3 = 0.001, 0.001, 0.001 for the prior variance on μ, a 1:3 = 12, 30, 25 for the shape of gamma prior on s -2 , b 1:3 = 1, 1, 1 for the rate of gamma prior on s -2 , and δ π = 1, 9, 1 for the Dirichlet prior on initial distribution π, respectively.

EM vs. MCMC
As already a 4-state Gaussian HMM has 23 free parameters applying EM is often difficult due to the existence of multiple local optima and the local convergence of EM. Consequently, the estimate has to be repeated many times with randomly initialized parameter values to find the most likely model. It should also be noted that not necessarily the model maximizing the likelihood exhibits the best performance in classifying aberrations 2. Additionally, applying any constraint in an EM settings, i.e. ordering of the state means, is harder than in MCMC. EM also produces inferior results on datasets exhibiting class imbalance, for example where one type of aberrations (observations for one HMM state) are rare or missing, while MCMC can overcome this problem using informative priors. In Table 2 we see that MCMC sampling performs better than EM on real biological data which is consistent with prior reports from Guha [24] and Shah [25] who also report difficulties with EM and better MCMC performances. In particular, for HBL-2 we observe that the best model in terms of F1-measure-which is not available in de novo analysis-is not the most likely model and the most likely model has very low precision and, consequently, worse F1-measure than MCMC sampling. On the simulated datasets, EM performs poorly if the dataset is imbalanced. As a result we observe slightly worse standard deviation for the precisions and F1-measures computed by EM in Table 2. We also notice from the confusion matrix of three classes-loss, neutral, and gain-that even though the mean F1-measure, recall, and precision (defined on two classes, neutral and aberrant) are high, due to label switching [17], EM does not distinguish gain from loss, and vice versa, very well (see Table 3). By re-ordering the already learned state means the label switching problem can be addressed, but that increases misclassification rate due to state collapsing as the parameter values, specially means of the Gaussians, become almost identical [17]. In contrast, Bayesian methods cope with class imbalance problem by applying order constraints. Moreover, using McNemar's test [35] on the combined result of the 500 datasets we have verified that the predictions based on EM are significantly different from the predictions made by Bayesian methods with p-values being less than 0.001 in both cases.

FBG vs. Approximate Sampling
In an ideal setting, like the 2-state HMM example, approximate sampling closely mimics the performance of FBG sampling up to moderate compression level. For the simulated and real dataset approximate sampling's performance is comparable to FBG's while achieving a speed-up of 10 or larger. We also observe that for higher compression levels approximate sampling reports small number of aberrant clones, which results in small tp and fp values, but large fn value. As a result, we see low recall and high precision rate when the compression  level is too high for a particular dataset (see the rows with w ≥ 4.0s D for HBL-2 in Table 2). From Figures 4,5,6, and 7 we observe that segmentations by both sampling methods are almost identical at the recommended width w value. In case of HBL-2, they differ from the ground truth in some places. They predict a few extra gain regions and outliers are generally predicted as gains. We, as well as Shah et al. [25], have noticed that the HBL-2 dataset has many outliers, and the variance of emission distribution of gain state 4 converges to a high value which tries to explain the outliers. In contrast, GBM has fewer outliers (see Figure 5) and approximate sampling seems robust to those outliers. As the compression algorithm forces possible outliers to be included in a compressed block, it is robust to moderate frequencies of outliers.

Width Parameter
By varying the width parameter w we can control the compression ratio γ and the speed-up achieved by approximate sampling. But from Table 1 and 2, and Lemma 1 it is also clear that by setting a large value one can get unfavorable results. We have analyzed the effect of different w values using a synthetic dataset with a controlled level of noise following [36]. Each dataset has three chromosomes with total probe counts 500, 750, and 1000. Ten aberrant regions of size 11-20 probes, randomly assigned gain or loss, are inserted in random positions of the 500 probe chromosome. Similarly, 15 aberrant regions of size 11-25 probes, randomly assigned gain or loss, are inserted into larger chromosomes. A noise component N(0, s) is added to the theoretical log-ratios -1,0,0.58 (loss, neutral, and gain respectively) to model the data. For a set of noise levels, s ranging from 0.1 to 0.5, 50 synthetic datasets are generated. We use a 3-state HMM with width parameter values in the range 0s D , ..., 4s D (where s D is the standard deviation of the dataset). Signal-to-noise ratio (SNR) is defined as 0.58 σ . In Figure 8 we plot the mean compression ratio, F1-measure, recall, and precision for 50 synthetic datasets and the real biological data HBL-2. For all noise levels the compression ratio drops exponentially with increasing values of w. Ideally, one would like to compress as much as possible without affecting the quality of the predictions. We visually identified a best value for width as the point after which the F1-measure drops substantially. Comparing the knee of the curve with the best value, we notice that while using the knee computed by Lmethod [27] is a conservative choice for width, in most cases we can still obtain a considerable speed-up.

Outliers
Gaussian HMMs are known to be sensitive to outliers which is evident from our results of HBL-2 and SNP arrays. Traditionally, outliers have been handled either by using a mixture distribution as the emission distribution or by preprocessing the data to remove possible outliers or impute more appropriate values. We have observed that a simple local median approach works very well to identify the outliers in a time series of log 2 -ratio values. Although using a mixture distribution or a distribution with fat tails, i.e. Student's-t distribution, is a better choice we lose a significant computational advantage in approximate sampling. Another approach, which we prefer in this context, is to use a HMM state with a very wide Gaussian and low selftransition probability to model the outliers. We have observed very good performance in this way. However, as our primary focus is to compare FBG-sampling with approximate sampling we choose to use a simple Gaussian model at the end.

Conclusions
Analyzing CGH data either from DNA microarrays or next generation sequencing to estimate chromosomal aberrations or investigate copy number variations (CNV), leads to the problem of segmenting sequences of observations which are essentially noisy versions of piecewise-constant functions. For reasons of efficiency, ML or MAP point estimates of HMM parameters combined with the Viterbi-algorithm to compute a most likely sequence of hidden states and thus a segmentation of the input are most popular in practice. This ignores research which clearly demonstrates that Bayesian approaches, where MCMC is used for sampling and thus integrating out model parameters, is more robust with higher recall and higher precision [24]. Additionally, our experiments show that likelihood is not informative with respect to the quality of CNV calls putting the use of ML into question even if the estimation problem could be solved. We propose a method using approximate sampling to accelerate MCMC for HMMs on such data. Our method constitutes the first use of ideas from spatial index structures for several consecutive observations and approximate computations based on geometric arguments for HMM; the effectiveness of this approach was previously demonstrated for k-means clustering, mixture estimation, and fast evaluation of a mixture of Gaussians.
We demonstrate that these very abundant biological CGH datasets, which measure chromosomal aberrations and copy number variations, are consistent with our assumptions of piece-wise constant ground truths, and we are able to achieve speed-ups between 10 and 60 respectively 90, on these biological datasets while maintaining competitive prediction accuracy compared to the state-of-the-art. As datasets with even higher resolution, both from higher density DNA microarrays and next generation sequencing, become available, we believe that the need for precise and efficient MCMC techniques will increase. The added precision over popular ML/MAP-based methods is of particular biological relevance as for complex neurodegenerative diseases such as Autism de-novo copy number variations have recently been shown to play a role [37]; a precise and quick analysis on large collectives of patients is desirable.
Applying approximate sampling to multi-dimensional observations-to jointly analyze data sets for recurrent CNVs [38] instead of analyzing individuals and postprocessing results-and considering more complicated HMM topologies and observation densities are directions for our future work.

Appendix
Proof. Using the assumption on τ, for any position i ≤ l ≤ i + n -1, we can argue that, For any partial state path q i , . . ., q i+n-1 , We partition S n , the set of all possible partial state paths of length n, into N subsets S s 1 . . . S s N such that, S sj = {S ∈ S n : (∀s l =sj C(S, sj) > C(S, sl)) ∨ ((∀s l =sj C(S, sj) ≥ C(S, sl)) ∧S1 = sj)} for 1 ≤ j ≤ N, where C(S, s) = The size of S n can be expressed in terms of total number of non-self-transitions present in a path, As the sets S s j are equal sized partitions of S n , Now we derive an upper bound of the contribution from state paths in S s x . In the following equations we make use of the fact that a state path with k non-selftransitions goes through at least k 2 non-s x states.
Similarly, we derive an upper bound of the contribution from state paths in S s y , where 1 ≤ y ≠ x ≤ N. Now we use the fact that, because of the pigeonhole principle any state path in S s y has to go through at least n N nons x states.
Applying (5) and (6) in (4)  Note: For simplicity of the notation, we follow the convention that μ x 0 = −∞ and μ x N+1 = ∞ so that the proof holds for x = 1 or x = N.