MixSIH: a mixture model for single individual haplotyping

Notation

Throughout the paper, we denote the number of elements of any set A by |A|, and the direct product set $\underset{n}{\underbrace{A\times \cdots \times A}}$ by A^⊗n. Let X = {1, 2, . . . , M} be the set of SNP loci, and $\mathcal{H}=\left\{0,1\right\}$ be the two haplotypes. It is convenient to introduce a phase vector Φ = φ₁ ... φ _M . The pair φ _j = (φ _j0 , φ _{j 1}) is referred to as phase, and represents the two alleles of haplotype 0 and 1 at site j, respectively. Because the haplotype assembly problem is trivial for homozygous sites, and because it is usually much easier to determine the genotype than to determine the haplotypes, it is often convenient to restrict the SNP loci X to heterozygous sites. Furthermore, if sequence-specific sequencing errors are not considered, it is convenient to use a simple binary representation of alleles; we randomly assign 0 to one of the two alleles at each heterozygous site j, and 1 to the other allele. In this case, the set of alleles is denoted by Σ = {0, 1}, and the set of possible phases is denoted by Δ = {(0, 1), (1, 0)}. We assume this binary representation throughout the paper.

Let F = {f _i |i = 1, . . . , N} be the set of input fragments which are supposed to be aligned to the reference genome, and each fragment f _i takes value f _ij ∈ Σ at locus j ∈ X if a nucleotide is aligned and equal to one of two alleles, and f _ij = ∅ if fragment f _i is unaligned, gapped, ambiguous, or a base different from the two alleles, at site j. For any subset X' ⊆ X, we say fragment f _i spans the sites X' if f _ij ≠ ∅ for all j ∈ X'. We refer to the subset of X spanned by fragment f as X(f). We say fragment f _i covers site j if there exists a pair of spanning two different (possible non consecutive) SNP sites j₁, j₂ ∈ X(f _i ) such that j₁ < j ≤ j₂. The set of fragments that cover site j is denoted by F ^c (j). Further, we refer to the set of all the possible haplotypes for sites X(f _i ) as $\mathrm{\Delta }\left(f_i\right)={\mathrm{\Delta }}^{\otimes \left|X\left(f_i\right)\right|}$.

The SIH problem takes a set of aligned SNP fragments F as input and outputs a hidden phase vector Φ (Figure 1). Because the SIH problem does not associate the inferred haplotypes $\mathcal{H}$with the real paternal and maternal chromosomes, the switched configuration $\bar{\mathrm{\Phi }}={\bar{\phi }}_1\cdots {\bar{\phi }}_M,{\bar{\phi }}_j=\left({\bar{\phi }}_{j\bar{0}},{\bar{\phi }}_{j\bar{1}}\right)$ with $\bar{0}=1$ and $\bar{1}=0$, must be regarded as a completely equivalent prediction. Therefore, SIH has no meaning if there is only one heterozygous site, and it is only meaningful if one considers co-occurrences of alleles on the same haplotype for two or more heterozygous sites.

Mixture model

is the probability that we observe σ ∈ Σ when the true allele is σ' ∈ Σ and α represents the sequence error rate which we assume is independent of fragments and positions.

We take α as a fixed constant because it is better estimated from other resources rather than from only the bases at the SNP sites. For example, we may estimate α by using the all the read sequences or by using information from other dedicated studies about sequencing and mapping errors. In the following, we use α = 0.1 unless otherwise mentioned and the dependency of the α is described in Additional file 1. We further assume the mixture probabilities are equal, p ^m (0) = p ^m (1) = 0.5, as they often converge to around 0.5. Therefore, the parameter set Θ that needs to be optimized consists only of the set of phase probabilities: $\mathrm{\Theta }=\left\{{\theta }_{j\nu }\right\}=\left\{p_j^{\mathrm{\Phi }}\left(\nu \right)\right\}$.

We explain the difference between our model and the models of Kim [14] and Li [16] in Additional file 1.

The minimum connectivity score

As described above, the two haplotypes $\mathcal{H}$ in the SIH problem have no particular identity and it is not possible to predict which of them converges to the actual paternal or maternal chromosome. In relation to this, the likelihood function P (F, H, Ψ| Θ) has a symmetry between the switched configurations: $P\left(F,\bar{H},\bar{\mathrm{\Psi }}|\bar{\mathrm{\Theta }}\right)=P\left(F,H,\mathrm{\Psi }|\mathrm{\Theta }\right)$, where $\bar{H}=\left\{{\bar{h}}_i|i=1,\dots ,N\right)$ and $\bar{\mathrm{\Psi }}=\left\{{\bar{\mathrm{\Phi }}}^{\left(i\right)}|i=1,\dots ,N\right\}$ represent the configuration that all the haplotype origins of the fragments are exchanged, and $\bar{\mathrm{\Theta }}=\left\{{\bar{\theta }}_{jv}\right\},{\bar{\theta }}_{jv}={\theta }_{j\bar{v}}$ are the switched phase probabilities. Therefore, the marginal likelihood $P\left(F|\mathrm{\Theta }\right)={\sum }_{H,\mathrm{\Psi }}P\left(F,H,\mathrm{\Psi }|\mathrm{\Theta }\right)$ is symmetric for the two parameter sets: $P\left(F|\bar{\mathrm{\Theta }}\right)=P\left(F|\mathrm{\Theta }\right)$.

Where ${\mathrm{\Theta }}^{\prime }=\left\{{\theta }_{jv}^{\prime }\right\}$ with ${\theta }_{jv}^{\prime }={\theta }_{jv}$ for j <j₀ and ${\theta }_{jv}^{\prime }={\bar{\theta }}_{jv}$ for j ≥ j₀. The second equality follows from the symmetry of P (F| Θ) described above, and shows that only the fragments covering site j₀ are necessary to compute the connectivity of site j₀. The connectivity measures the resilience of the assembly result against swapping the two haplotypes 0 and 1 in the right part j = j₀, . . . , M of the sites. We refer to this change of parameters Θ → Θ' as twisting the parameters at site j₀.

We extract confidently assembled regions by selecting the pairs (j₁, j₂) with high MC values. From the above definition, it is obvious that if the MC value is higher than a given threshold for some pair (j₁, j₂), then all the pairs inside range [j₁, j₂] have MC values higher than the threshold. In this sense, MC(j₁, j₂) can be considered as defined on the range [j₁, j₂].

Variational bayesian inference

where Z ^{H Ψ} is a normalization constant, β _ihjν and λ _jν represent the hyperparameters that specify the posterior distributions, and Dir(θ _j |λ _j ) is the Dirichlet probability distribution of | Δ| parameters. Because Q ^{H Ψ}(H, Ψ) and Q^Θ(Θ) are connected through the dependencies among the hyperparameters, they cannot be found simultaneously. Therefore, we optimize β _ihjν and λ _jν by an iterative method.

Here, the twist of hyperparameters Λ = {λ _jν } is defined similarly to that of parameters Θ = {θ _jν }. We describe the details of this procedure in Additional file 1.

Inferring haplotypes

We select the phase ν at site j for which this $p_j^{\mathrm{\Phi }}\left(v\right)$ is the highest. We limit the predicted haplotype segments to the regions with high MC values.

Possible extensions of the model

In this paper, we consider only the binary representation of heterozygous sites. We also constrain the error rate to be constant throughout the sequence. However, some of these constraints are easily removed. We can include homozygous sites and four nucleotide alleles by expanding the phase set Δ. For example, the phase set of a multi-allelic variant is represented like Δ = {(A,C),(A,G),(C,A),(C,G),(G,A),(G,C)}. We can even include small structural variations if they can be represented by additional allele symbols and the phase set of a structural variant is represented such as Δ₁ = {(A,-),(-,A)} for indel and Δ₂ = {("AC","ACAC"),("ACAC","AC")} for short tandem repeats. With these extensions, the accuracy of genotype calling of multi-allelic variants from sequencing data might be improved by considering haplotypes simultaneously [30] and the accuracy and the recall of the haplotype region might be improved because all variant sites add information to infer the derivation of the fragments. Furthermore, we can make the error probability matrix p ^e (σ|σ') dependent on the alleles of each fragment, which may be useful for incorporating the quality scores of sequenced reads.