Gibbs-sampling

where X=(x₁,x₂,…,x _d ) is a d-dimensional vector and Z is a normalization factor. The function ϕ specifies the edge potential for two variables with an edge between them. We would like to collect samples of X based on this distribution P(X).

Readers can refer to [21] for a more detailed description of Monte Carlo Markov Chain.

Polyploid haplotype phasing

The inputs for the polyploid haplotype phasing problem is the ploidy k (the number of haplotypes to be phased), the set of aligned sequencing reads X (We assume the raw reads have been aligned to a reference sequence and thus the SNPs spanned by the reads are identified already), a sequencing error rate ε. The VCF (Variant Call Format) file containing the SNP positions and dosages could be optional. The dosage information gives the number of reference alleles and alternative alleles for a given SNP position and therefore can be used to reduce the phasing search space and to improve the phasing accuracy. Notice for some programs such as HapCompass, the dosage information is mandatory.

The output of the phasing algorithms is the k phased haplotypes. There are a few popular metrics to evaluate the performance of the phasing algorithms, such as MEC (minimum error correction) [22], minimum fragment removal, MSR (minimum single nucleotide polymorphism (SNP) removal) [23], and two recent models MFC (maximum fragments cut) [24] and BOP (balanced optimal partition) [25]. In this work, we focused on MEC as our metric.

MEC

where X is a set of m sequencing reads, H is the set of n haplotypes, x _j is the j-th read, h _k is the k-th haplotype, D(x _j ,h _k ) is the number of mismatches between x _j and h _k ,r _jk is 1 if the j-the read is assigned to the k-the haplotype and 0 vice versa. A read is assigned to a haplotype which minimizes its number of mismatches. Notice mismatches only occur at SNP positions as all other positions are homozygous. Thus the polyploid haplotype phasing problem is to phase k haplotypes H given the set of reads X and the SNP positions so that the objective function MEC(X,H) is minimized. It is known [8] that minimizing MEC is NP-hard even for diploid haplotype phasing when the length of the reads is greater than one.

Recently various methods have been proposed to the polyploid haplotype phasing problem. HapCompass [16] builds a compass graph from the sequencing reads, which is an undirected weighted graph. In the compass graph, the vertices are the SNPs and an edge between a pair of SNPs indicates that at least one read spans the two SNPs. There is an integer weight associated with each edge. It is shown that a compass graph has a unique phasing if it has no conflicting cycles, which is a simple cycle that contains either an odd number of negative edges or at least one 0-weight edge or both. Haplotype phasings correspond to spanning trees in the graph. The phasing problem is converted to a minimum weighted edge removal optimization on the graph and an algorithm based on cycle basis local optimizations for resolving conflicting cycles is proposed.

HapTree [17] aims to trim down the search space for all the possible haplotypes to a much smaller set of more likely solutions. It takes an inductive approach, generating a collection of likely phasing solutions for the first two SNPs in the genome, and then extending those to phasing solutions of the first three SNPs, and those to the first four SNPs, and so on. When extending any particular solution, HapTree chooses (based on computing likelihoods) how the alleles of the newly added SNP may be assigned to chromosomes; it includes only those assignments that are sufficiently likely. Upon including all SNPs to be phased, HapTree randomly chooses a solution of maximum likelihood from amongst the solutions it has found. It is shown [18] that the trimming process might be time consuming for some cases.

where k _sim denotes the number of overlapping positions where the reads have an identical alleles and k _dissim is the number of positions where they are different. Then giving the graph and the ploidy, SDHaP aims to find k−1 cuts such that the sum of intra-partition edge weights is maximized and inter-partition edge weights is minimized and the problem is solved via correlation clustering. SDhaP formulates the problem as a semi-definite program (SDP), and employs a low-rank Lagrangian scheme followed by randomized projections and a greedy refinement of the k-ploid haplotypes to solve the SDP.

H-PoP and H-PoPG [19] try to partition the DNA reads sequenced from a k-ploid organism into k groups such that the reads of the same group share the same alleles on as many SNP loci as possible and the reads from different groups are different on as many loci as possible. Heuristic strategies are proposed by limiting the number of intermediate solutions at each iteration of a dynamic programming algorithm. Notice H-PoP assumes no VCF file, which contains the dosage information of the variants. is provided while H-PoPG accepts VCF file as an input.

The polyploid haplotype phasing method we proposed in this work is an extension of Harsh [14], which applies the Gibbs Sampling algorithm to phase diploid haplotypes (namely two haplotypes). For diploid scenario, the two haplotypes are complimentary. Therefore indeed only one haplotype needs to be phased. For polyploid scenario, the haplotypes are more than two and they can share common segments and are not necessarily complimentary. Thus the problem is much more challenging. Harsh can not be applied to polyploid haplotype phasing in that the conditional probability estimation for the Gibbs Sampling process needs to be completely re-invented.

Poly-Harsh

Gibbs sampling

In this work, we developed a Gibbs Sampling based method Poly-Harsh for polyploid haplotype phasing. Our algorithm consists of two major steps: fix the haplotypes, compute read assignments; then fix the read assignments, compute the haplotypes. Here when we compare a haplotype to the reference haplotype, if the allele is the same as the reference, the genotype value is 0. If the allele is alternative, the genotype value is 1. In cases where the allele is neither the reference nor the alternative, we simply assign the genotype value as 1. For each SNP position i, we define a genotype value vector h _i = [ g_1,i,g_2,i,…,g_k,i], where g_j,i is the genotype value of the j-th haplotype at the i-th SNP. For a given ploidy k (for illustration purpose and without losing generality, assuming k = 4), the genotype value vector could be one of 2 ^k binary vectors [1, 0, 0, 0], [1, 1, 0, 0], …, [0, 1, 1, 1]. Notice we ignore vectors [0,0,0,0] and [1,1,1,1] as they indicate the allele is homozygous. We also define a read assignment vector r as a k-dimensional vector r= [ a₁,a₂,…,a _k ], where a _i = 1 if the read is assigned to the i-th haplotype and a _i = 0 if the read is not assigned to the i-th haplotype. For example, for k = 4, if the read is assigned to the first haplotype, r= [ 1,0,0,0]. If the read is assigned to the second haplotype, r= [ 0,1,0,0] and so on so forth. Notice a read can not be assigned to more than one haplotype and on the other hand it has to be assigned to one haplotype. Given the genotype value vector h _i for the i-th SNP and a read assignment vector r _j for the j-th read x _j , we define a function as below:

$$\begin{array}{@{}rcl@{}} \theta(h_{i}, r_{j}, x_{j}) & = & ln(1-\epsilon)^{t} + ln(\epsilon)^{k-t} \\ t & = & match(h_{i}, r_{j} \times x_{j,i}) \end{array} $$

(5)

where ε is the sequencing error rate, x_j,i is the i-th value of the x _j ,match(A,B) is the vector-wise matches between two vectors A and B where the number of matches is increased by 1 if two vector elements are identical (either both 1 or both 0). For example, for h _i = [ 1,1,0,0], r _j =[1,0,0,0],x_j,i=1, we have r _j ×x_j,i= [ 1,0,0,0] and thus t =match(h _i ,r _j ×x_j,i)= 3. For h _i = [ 1,1,0,0],r _j =[1,0,0,0],x_j,i=0, we have r _j ×x_j,i=[0,0,0,0] and thus t =match(h _i ,r _j ×x_j,i)= 2. The θ function essentially models the probability of the correct read assignment given the matches between the read and the haplotypes. Notice the function considers the mismatches due to sequencing error.

Given ploidy as k, the set of genotype value vectors \(H =\, [\!h_{1}, h_{2}, \dots, h_{n}]\) where \(h_{i} =\, [g_{1,i}, g_{2,i}, \dots, g_{k_{i}}]\) as we have defined and n is the number of SNPs, the set of read assignment vectors \(R = \,[\!r_{1}, r_{2}, \dots, r_{n}]\) where \(r_{i} = [a_{1}, a_{2}, \dots, a_{k}]\) is the assignment vector for the i-th read as we have defined, the sampling process proceeds as follows: We first randomly initiate H, then we fix H and compute the conditional probability P(R|H). We sample the reads assignment R based on P(R|H). Next we fix R and compute the conditional probability P(H|R). We sample the genotype value vectors H based on P(H|R). For ploidy k, we have 2 ^k haplotype values for a specific SNP. Assuming the genotype value vector is h and the set of reads X= [ x₁,x₂,…,x _n ], we could compute its probability as

$$\begin{array}{@{}rcl@{}} P(h | R) = \frac{exp\left(\sum_{j=1}^{n} \theta(h, r_{j}, x_{j})\right)}{exp\left(\sum_{i=1, j=1}^{i=2^{k}, j=n} \theta(h_{i}, r_{j}, x_{j})\right)} \end{array} $$

(6)

where the function θ is defined in Eq. 5, h _i is the i-th genotype value vector, r _j is the assignment vector for the j-th read x _j . Then we apply sampling to update the genotype values of the SNP. We do a similar Gibbs sampling step for read origin given fixed haplotypes. Again, we conduct sampling to update the read assignment vector on the fixed haplotypes H for a given read x as below:

$$\begin{array}{@{}rcl@{}} P(r | H) = \frac{exp\left(\sum_{j=1}^{2^{k}} \theta(h_{j}, r, x)\right)}{exp\left(\sum_{i=1, j=1}^{i=k, j=2^{k}} \theta(h_{j}, r_{i}, x)\right)} \end{array} $$

(7)

where the function θ is defined in Eq. 5, h _j is the j-th genotype value vector, r _i is the assignment vector when the read x is assigned to the i-th haplotype.

Given H and R, we can easily compute the MEC of the phasing. Notice the k haplotypes can be constructed from the genotype value vectors H, by concatenating all the genotype values from the same haplotype. For example, haplotype one can be constructed as [g_1,1,g_1,2,…,g_1,n]. Therefore, the haplotype phasing problem becomes identifying the optimal H that minimizes MEC.

Notice for diploid scenario and polyploid scenario, the computations for P(r|H) and P(h|R) are significantly different, as for polyploid scenario, r and h are k-dimensional vectors. The probabilities not only depends on the haplotype the read is assigned to, but also depends on the remaining haplotypes: if one haplotype has a mismatch at a SNP position to a read and the read is not assigned to the haplotype, it is considered as a correct operation to not assign the read to the haplotype.

We repeat the above two steps iteratively. For each iteration, we compute the MEC score of the reads. The process converges when the MEC does not improve or we have reached certain number of iterations. As Gibbs Sampling is sampling based and its performance is affected by the initial random simulation of H, it may fall into local optimum solutions. Therefore we re-run the program multiple times, each time started from a different random seed H. This helps the program to escape from local optimum. A pseudocode of the algorithm is shown in Algorithm 1. Notice the algorithm shows only one run of the procedure. If we would like to run the algorithm multiple times, we need to repeat steps 1 to 5 multiple times, each time with a different randomly initialized H.

Generate candidates

At first we recursively construct all possible candidate haplotypes for each sample. Let the number of samples and length of each candidate be n and t, respectively. Without loss of generality, let each sample consists of k blocks of broken haplotypes. If each block contains 4 haplotypes, the number of possible candidates will be 4 ^k as stated above. We recursively construct candidates for each sample. The time complexity for constructing all possible candidate haplotypes for m samples is \(O\left (4^{k}mt\right)\) which is exponential with respect to k. It is computationally very expensive task when k is large and the execution time can be very large. However in reality k is very small and we can safely assume that 1≤k≤8.

Detect frequent candidates

Frequent candidates are those candidate haplotypes which occurs multiple times (i.e., at least twice) across m samples. This set of candidates contains highly accurate haplotypes because of their multiple occurrences. Exact algorithm is quadratic in the number of candidate haplotypes \(\left (\text {i.e}., O\left (\left (4^{k}m\right)^{2}t\right)\right)\). To reduce the runtime we are proposing a randomized algorithm. The expected runtime is sub-quadratic in the number of candidate haplotypes and at the same time the algorithm is also highly accurate in computing the set of frequent candidates. At first we illustrate how to find most similar pair of candidates. By naturally extending it, we will describe the process of finding the set of frequent candidates. We describe our randomized algorithm next.

Suppose we are given n vectors \(\hat {b}_{1},\hat {b}_{2},\ldots,\hat {b}_{n}\) each of length t. The problem is to find the pair of vectors that are the most similar (i.e., the Hamming distance between them is the smallest). Note that, given two vectors, we can find the Hamming distance between them in O(t) time. A straight forward algorithm to identify the most correlated pair of vectors takes \(O\left (n^{2}t\right)\) time. This algorithm computes the Hamming distance between every pair of vectors. We can achieve a better run time using randomization. We say that the correlation between a pair of strings is p if the Hamming distance between them is t(1−p). Let p₁ be the correlation between the most correlated pair of strings and p₂ be the correlation between the second most correlated pair of strings.

The idea of our algorithm is to iteratively collect pairs of strings that are candidates to be the most correlated. Once we collect enough pairs, we compute the distance between each pair in this collection and output the closest. In each iteration we pick q columns randomly. For any vector (or string), the values in these columns can be concatenated to get a q-bit integer. We hash the vectors based this integer value. Subsequently, we generate pairs as follows: Consider any bucket in the hash table. If there are m vectors in this bucket, then each pair of vectors in this bucket is added as a candidate to a list C. There are \(O\left (n^{\frac {\text {log} p_{1}}{\text {log} p_{2}}}\text {log}\ n\right)\) iterations in the algorithm. We can show that after \(O\left (n^{\frac {\text {log} p_{1}}{\text {log} p_{2}}}\text {log}\ n\right)\) iterations, C will have the most correlated pair of bulbs with a high probability (i.e., with a probability of 1−n^−Ω(1)). For details readers are refereed to [26].

We can naturally extend the above algorithm to get the frequent candidate haplotypes. In each iteration we compute similarity coefficients for all the the pairs in each hash bucket. The similarity coefficient is defined as \(SC = \frac {t - d}{t}\) where t is the length of the haplotype and d is the Hamming distance between a pair of interest. The more the similarity coefficient, the more similar will be the pair. We retain that pair which has the similarity coefficient ≥ a threshold, V. The individual candidate haplotype belonging in each pair is then hashed into a hash map H. The keys and values of hash map H is the candidate strings and frequency of occurrences, respectively. Finally we sort H with respect to its values (i.e., the frequency of occurrences) and output candidates occurring at least twice across m samples. For each iteration the expected number of pairs generated is O(n) as described in [26]. The expected time to build the hash map H will be also O(n) for a particular iteration. As the expected time to resolve a collision (i.e., look up and update) for H is O(1), the total time spent for all the iterations will be \(O\left (n^{1+\frac {\log p_{1}}{\log p_{2}}}\right)\). Since we can sort the map by using any integer sorting algorithm, the expected runtime of this step is \(O\left (n^{1+ \epsilon }t\log n\right)\) where n=4 ^k m and 0<ε<1. We have used fixed number of iterations (i.e., 50) in the experimentations. Let the number of frequent candidates be s.

Detect true candidates

Frequent set of candidates contains all the haplotypes shared twice across the samples. Each sample has 4 ^k candidate haplotypes. We need to find 4 haplotypes for each sample. Let S^′ and S^″ be the sets of 4 ^k candidate haplotypes of a particular sample and frequent candidates, respectively. Intuitively if we find the closest (i.e., most similar) pair of haplotypes (c^′∈S^′,c^″∈S^″), c^′ will be one of the 4 haplotypes of a particular sample. Haplotypes c^′ and c^″ may be identical but it is not guaranteed. c^′ may not be in the set of frequent candidates due to large number of errors. In this case the closest one (i.e., c^″) from s gives us the best possible information to find c^′. Next we illustrate this step algorithmically.

Each sample has 4 ^k candidates and the size of pruned set is s as described in 3. For each candidate c^′∈S^′ and c^″∈S^″, we compute pair-wise similarity coefficient. We then sort all the pairs (c^′,c^″) with respect to SC in non-increasing order and then rank of frequent candidates in non-decreasing order. We traverse the sorted list of pairs and collect first 4 candidate haplotypes c^′ with the following restrictions: (1) No candidate haplotype c^′ will be chosen more than once; and (2) Each subsequence from each block will be used exactly once. If the size of each candidate is t, we need \(O\left (4^{k}st\right)\) time to compute pair-wise similarity coefficients (since the number such pairs is 4 ^k s). Sorting the list of pairs with respect to similarity score and then by rank of frequent candidate haplotypes will take no more than \(O\left (4^{k}s(k+\log s)\right)\) time. After sorting detecting first 4 haplotypes could take \(O\left (4^{k}st\right)\) time. In total the time complexity to find 4 candidate haplotypes for m samples will be \(O\left (4^{k}ms(k+t+\text {log}\ s)\right)\).

Polyploid haplotype phasing

In order to show that re-running Poly-Harsh multiple times in usual leads to better performance as the program could escape from the local optimum, we show in Fig. 1 the minimum, the mean and the standard deviation of 8 rounds of running of Poly-Harsh with respect to each parameter settings. As we can see that for all parameter settings, with 8 rounds of running, the minimum MEC we obtained is in usual much better than the mean MEC, especially when the MEC is relatively large. This indicates that running the program multiple times could in general improve the performance. In our experiments, we see that in general the performance converges in 8 rounds of running. Also we can observe that larger coverage and longer reads lead to larger MEC and larger standard deviation.

Next we compare the phasing performance of Poly-Harsh with Hapcompass [16] and H-PoP/H-PoPG [19]. We again take the parameter settings specified in Table 1. For each parameter setting, we randomly simulate 10 data sets and we show the average performance. For Poly-Harsh, we conducted 8 rounds of running, each with 100 iterations. Notice that HapCompass does require the dosage information while both Poly-Harsh and H-PoPG (the version of H-PoP that takes VCF) take the dosage information as optional. Therefore, for a fair comparison, we first feed the dosage information to all three methods. We show the MEC of the three methods in Fig. 2. We can see that both Poly-Harsh and H-PoPG achieved better results compared with HapCompass. Poly-Harsh achieved the best performance for all experiments.

Next we feed only the aligned sequencing reads information to H-PoP and Poly-Harsh for comparison purpose. HapCompass is excluded in the experiments as it requires dosage information. We show that the MEC of both methods in Fig. 3. We can see again Poly-Harsh achieved better results compared with H-PoP. Also the MEC is in general larger when dosage is not provided, indicating that without dosage, the phasing becomes less accurate.

Finally, the execution time of Poly-Harsh depends on the number of rounds of running. For one round of running, on our data sets, the execution time for all three methods are comparable, all less than 1 second. The execution time of Poly-Harsh increases linearly with respect to the number of rounds of running.

Contiguous haplotype reconstruction

We created 15 simulated datasets each having 30 samples. As described above each sample contains 4 haplotypes. Each haplotype is a binary string of length 1,300 bps. Binary string is created by concatenating randomly chosen binary values (i.e., either 0 or 1) under uniform distribution. As the phased haplotypes usually contain errors, we introduce errors in the haplotypes by flipping bases according to given error rates. For example, if an error rate is 1%, the probability that a base in a haplotype will be flipped is 0.01. Each sample is randomly partitioned into a number of blocks. The number of such blocks will be 1 through 8. As described we are able to construct a haplotype correctly, if it occurs at least twice across the dataset of interest. To make the datasets more realistic, we introduce the concept of Shared. If a dataset is told to be X% shared, it means X haplotypes occurs 2× across all the 30 samples.

Now consider datasets D6-D10. 80% haplotypes in each dataset occur twice across the samples. Thus the sensitivity should be at least 80%. But the correctness of our algorithm is above 90%. It is due to the fact that if we can construct 3 haplotypes of a sample correctly, the rest will be formed accurately. The same observation applies to datasets D11-D15. We visualize the results in Figs. 4 and 5.

As our method is the very first algorithm that tries to construct contiguous haplotypes from the phased haplotype blocks, we compared our method against a simple baseline method: concatenate the blocks randomly. The baseline method has a sensitivity of around 25% among all scenarios, as the random concatenation doesn’t rely on any information from the dataset. Thus our method is much more accurate than the random concatenation algorithm.

The experimental evaluations show that our algorithm is indeed effective and efficient in terms of both accuracy and runtime. For synthetic dataset, our algorithm achieves nearly 100% accuracy where the datasets have less errors and modest number of shared haplotypes. Here accuracy is defined as the fraction of haplotypes constructed correctly. In some datsets the median of accuracy is 70−80%. This is due to the fact that the haplotype construction may be affected if the dataset of interest has very low discriminative power. In this case, there is no sufficient information to distinguish true and false haplotypes accurately. It has cumulative effects also. Each sample has 4 haplotypes as stated above. Suppose we are not able to correctly construct the first haplotype. Then there is a high chance that the rest will be constructed incorrectly. This case arises when (1) dataset contains large number of errors and/or (2) there is little shared information across the samples. In case 2 it is impossible to construct all the haplotypes accurately. In order to construct haplotypes correctly, 3 out of 4 haplotypes of each sample must be occurred at least twice across the samples. In case 1 shared information may be lost because of noise (such as, missing values, phasing errors) in the dataset.