The first method for phasing polyploid genomes was HapCompass [ 2 ], which uses a graphical approach. Popular methods that followed include HapTree [ 4,5 ], H-PoP [ 36 ] and SDhaP [ 11 ]. HapTree and H-PoP heuristically maximize a likelihood function and an objective function based on the MEC model respectively while SDhaP takes a semi-definite programming approach. HapTree-X [ 5 ] additionally incorporates long-range expression correlations to allow phasing even of pairs of variants that cannot be covered by a single read, overcoming some of the problems with short-read phasing.

Due to the increased prevalence of long-read data from Oxford Nanopore or PacBio, newer methods taking advantage of the longer-range correlations accessible through long-read data have been proposed [ 1,33 ]. Unfortunately, because the error profiles of long-read technologies differ considerably from Illumina short-reads (e.g. a higher prevalence of indel errors compared to SNPs), methods tailored to short-reads such as [ 25,34 ] may be ineffective, so altogether new paradigms are required.

At a more theoretical level, in the diploid setting, the standard minimum error correction (MEC) model [ 6 ] has proven to be quite powerful. It is known to be APX-Hard and NP-Hard but heuristically solved in practice with relatively good results. Unfortunately, a good MEC score may not imply a good phasing when errors are present [ 21 ]. This shortcoming is further exacerbated in the polyploid setting because similar haplotypes may be clustered together since the MEC model does not consider coverage; this phenomenon is known as genome collapsing [ 33 ]. Thus, although the MEC model can be applied to the polyploid setting, it may be suboptimal; however, there is yet to be an alternative commonly agreed upon formulation of the polyploid phasing problem. Indeed, this is reflected in the literature: the mathematical underpinnings of the various polyploid phasing algorithms are very diverse. 1.2

In this paper, we first address the theoretical shortcomings highlighted above by giving two new mathematical formulations of polyploid phasing. We adopt a

Probabilistic and graphical haplotype phasing 3 probabilistic framework that allows us to (1) give a better notion of haplotype similarity between reads and (2) define a new objective function, the uniform probabilistic error minimization (UPEM) score. Furthermore, we introduce the idea of framing the polyploid phasing problem as one of partitioning a graph to minimize the sum of the max spanning trees within each cluster, which we show is related to the MEC formulation in a specific case.

We argue that these formulations are better suited for polyploid haplotype phasing using long-reads. In addition to our theoretical justifications, we also implemented a method we call flopp (fast local polyploid phaser). flopp optimizes the UPEM score and builds up local haplotypes through the graph partitioning procedure described. When tested on simulated data sets, flopp produces much more accurate local haplotype blocks than other methods and also frequently produces the most accurate global phasing. flopp’s runtime is additionally comparable to, and often much faster than, its competitors.

The code for flopp is available at https://github.com/bluenote-1577/flopp. flopp utilizes Rust-Bio [ 19 ], is written entirely in the rust programming language and is fully parallelizable. flopp takes as input either BAM + VCF files, or the same fragment file format used by AltHap [ 16 ] and H-PoP. 2 2.1

We represent every read as a sequence of variants (rather than as a string of nucleotides, which is commonly used for mapping/assembly tasks). Let R be the set of all reads that align to a chromosome and m be the number of variants in our chromosome. Assuming that tetra-allelic SNPs are allowed, every read ri is considered as an element of the space ri ∈ {−, 0, 1, 2, 3}m. A read in this variant space is sometimes called a fragment. Denoting ri[j] as the jth coordinate of ri, ri[j] ∈ {0, 1, 2, 3} if the jth variant is contained in the read ri where 0 represents the reference allele, 1 represents the first alternative allele, and so forth. ri[j] = − if ri does not contain the jth allele.

We note that flopp by default only uses SNP information, but the user may generate their own fragments, permitting indels and other types of variants to be used even if there are more than four possible alleles. The formalism is the same regardless of the types of variants used or the number of alternative alleles.

For any two reads r1, r2, let d(r1, r2) = |{k : r1[k] 6= r2[k], (r1[k] 6= −) ∧ (r2[k] 6= −)}| s(r1, r2) = |{k : r1[k] = r2[k], (r1[k] 6= −) ∧ (r2[k] 6= −)}|. d and s stand for different and same, representing the number of different and same variants respectively between two reads.

We use k to denote the ploidy. Given a k-ploid organism, a natural approach to phasing is to partition R into k subsets where the cluster membership of a read represents which haplotype the read belongs to. Let R1, . . . , Rk be a partition of R. Given a partition P = {R1, . . . , Rk}, we denote P [i] = Ri. 4

J. Shaw and Y.W. Yu.

Define the consensus haplotype H(Ri) ∈ {−, 0, 1, 2, 3}m associated to a subset of reads as follows. For all indices l = 1, . . . , m let H(Ri)[l] = arg maxa |{r ∈ Ri : r[l] = a}| and break ties according to some arbitrary order. If only − appear at position l over all reads, we take H(Ri)[l] = −. It is easy to check that H(Ri) is a sequence in {−, 0, 1, 2, 3}m such that H(Ri)[k] 6= − at indices for which some read overlaps, and Pr∈Ri d(H(Ri), r) is minimized.

In our formalism, we can phrase the MEC model of haplotype phasing as the k task of finding a partition {R1, . . . , Rk} of R such that Pi=1 Prj∈Ri d(rj , H(Ri)), which is called the MEC score, is minimized. For notational purposes, for a subset Ri ⊂ R, define S(Ri) = Pr∈Ri s(H(Ri), r) and D(Ri) = Pr∈Ri d(H(Ri), r).

k Pi=1 D(Ri) is just the MEC score for a particular partition. 2.2

Min-sum max tree partition (MSMTP) model. Let G(R) = (R, E, w) be an undirected graph where the vertices are R and edges E are present between two reads r1, r2 if r1, r2 overlap, i.e. d(r1, r2) + s(r1, r2) > 0. Let the weight of e = (r1, r2) be w(e) = w(r1, r2) for some weight function w. We call G(R) the read-graph, see Figure 1. A similar notion is found in [ 11,1,33,24,23 ].

For a partition of R into disjoint subsets {R1, . . . , Rk} we take G(Ri) as defined above. We only consider partitions of vertices such that all G(Ri) are connected, which we will denote as valid partitions. Let M ST (G) be the maximum spanning tree of a graph G. Define k SM T PRk (R1, . . . , Rk) = X X i=1 e∈MST (G(Ri)) w(e).

(1)

We formulate the min-sum max tree partition (MSMTP) problem as finding a valid partition {R1, . . . , Rk} of R such that SM T PRk (R1, . . . , Rk) is minimized. We refer to the objective function being minimized as the SMTP score and the computational problem of minimizing the SMTP score as the MSMTP problem. - - - - 0 0 0 0

1 1 1 1 - - - d = 1 s = 3 d = 4 s = 0 0 - 0 0 0 1 0 0 d = 5 s = 1 1 0 1 - 1 1 1 1 d = 3 s = 0 d = 1 s = 2

R1 

R2 H(R1) - 0 0 0 0 0 0/1 0 0 H(R2) - 1 0/1 1 1 1 1 1 1

MEC = 2 Fig. 1. An example of a read-graph (without the weighting w specified) along with corresponding d, s values between reads. The colors represent a partition of the readgraph into two subsets R1, R2. The consensus haplotypes H(R1), H(R2) are shown where the 0/1 indicates that either 1 or 0 are valid alleles for a consensus haplotype.

We now discuss the algorithms implemented in flopp. A high-level block diagram showing outlining flopp’s main processing stages are outlined in Figure 2. Importantly, flopp is a local clustering method, which means that we first attempt to find good partitions for smaller subsets of reads by optimizing the SMTP and UPEM functions and then joining haplotype blocks together afterwards. Choice of edge weight function Previous methods that use a read-graph formalism such as SDhaP [ 11 ], FastHap [ 24 ], WHP [ 33 ], and nPhase [ 1 ] all define different weightings between reads. Two previously used weight functions are wSDhaP (r1, r2) = s(r1, r2) − d(r1, r2) s(r1, r2) + d(r1, r2) (SDhaP) werror(r1, r2) =

s(r1, r2) s(r1, r2) + d(r1, r2) (nPhase; FastHap is similar).

These weightings are quite simple and have issues when the length of the reads have high variance, as is the case with long-reads. A more sophisticated approach is to use a probabilistic model. Let E(r1, r2) = 1 − werror(r1, r2) be

Suppose p < a < 1. Let H = a log ap + (1 − a) log( 11−−ap ) be the Kullback-Leibler divergence between an a−coin and a p−coin. Then

Pr(Binom(n, p) ≥ an) ≤ e−nH . as

. (3) 8

J. Shaw and Y.W. Yu.

In particular, this bound is asymptotically tight up to logarithms, i.e.

log Pr(Binom(n, p) ≥ an) lim = 1. (4) n→∞ −nH

The proof of the above theorem is simple; it follows by an application of a Chernoff bound and Sterling’s inequality. In [ 3 ], they show that this approximation is particularly good when a >> p. Theorem 2 gives an interpretation of w(r1, r2) as the negative log value of a 1-sided binomial test where the null hypothesis is that r1, r2 come from the same haplotype, the number of trials is s(r1, r2) + d(r1, r2) and the number of successes is d(r1, r2).

WHP [ 33 ] uses a log ratio of binomial probability densities as their edge weight. Mathematically, the resulting formula is almost the same as Equation 2, except they fix E(r1, r2) to be the average error rate between reads from different haplotypes, a constant. They end up with both negative and positive edges with large magnitude and Theorem 2 does not apply to their case. On the other hand, our edge weights are rigorously justified as negative log p-values for a binomial test, making it more interpretable.

MSMTP Algorithm We devised a heuristic algorithm for local clustering inspired by the MSMTP problem. From now on, let S ⊂ R. The pseudocode is shown in Algorithm 1. In the diploid setting, the main idea behind this algorithm is similar to that of FastHap [ 24 ], but the algorithm is quite different after generalizing to higher ploidies.

Algorithm 1: Greedy min-max read partitioning

Input : Read-graph G(S), ploidy k, iterations n

Output: A partition {S1, . . . , Sk} of S 1 {v1, . . . , vk} ←FindMaxClique(G(S), k) 2 for i = 1 to k do 3 Si ←{vi} 4 end 5 for i = 1 to n do 6 7

k V ← G(S) \ Si=1 Si Reverse-sort V by assigning to v ∈ V the value v → minSi∈{S1,...,Sk} maxr∈Si s(r, v) + d(r, v) V ← V [: d |Vn| e] for v in V do

S0 ← arg minSi maxr∈Si w(v, r)

S0 ← S0 ∪ {v} 8 9 10 11 12 end 13 end 14 Return {S1, . . . , Sk}

For the FindMaxClique method mentioned in Algorithm 1, we use a greedy clique-finding method which takes O(k|S| log |S| + k2|S|) time. First, we take the

clustering procedure by optimizing the UPEM score using a method similar to the Kernighan-Lin algorithm [ 18 ]. Pseudocode can be found in the appendix as Algorithm 2.

The algorithm takes in a partition P = {S1, . . . , Sk} of a subset of reads S ⊂ R as well as a parameter n which indicates the maximum number of iterations. The algorithm checks how moving reads from one partition to another partition changes the UPEM score for every read. We store the best moves and execute a fraction of them. Proceed for n iterations or until the UPEM score does not improve anymore. In practice, we set the parameter n = 10. The time complexity of Algorithm 2 is O(n|S|k log(|S|k)).

Local phasing procedure Note that Algorithms 1 and 2 work on subsets of reads or subgraphs of the underlying read-graph. Let b ∈ N be a constant representing the length of a local block. We consider subsets B1, . . . , Bl ⊂ R where

Bi = {r ∈ R : ∃j, b · (i − 1) ≤ j ≤ b · i, r[j] 6= −}.

The subsets are just all reads that overlap a shifted interval of size b, similar to the work done in [ 31 ]. After choosing a suitable b, we run the readpartitioning and iterative refinement on all B1, . . . , Bl to generate a set of partitions P1, . . . , Pl. We found that a suitable value of b is the 31 −quantile value of read lengths. By read length we mean the last non 0−0 position minus the first non 0−0 position of r ∈ {0, 1, 2, 3, −}.

It is important to note that computationally, the local clustering procedure is easily parallelizable. The local clustering step has therefore a linear speed-up in the number of threads. 10 2.4

J. Shaw and Y.W. Yu.

Filling in erroneous blocks Once we obtain a set of partitions P1, . . . , Pl of B1, . . . , Bl according to the above local clustering procedure, we can identify partitions with low UPEM score and correct them. We describe this procedure in Appendix C.

Local cluster merging Let P represent the final partition of all reads R. We build P given P1, . . . , Pl as follows. Start with P = P1. Let Sk be the symmetric group on k elements, i.e. the set of all permutations on k−elements. At the ith step of the procedure, let

k σi = arg max X |P [j] ∩ Pi+1[σ(j)]|.

σ∈Sk j=1 Then let P [j] = P [j] ∪ Pi+1[σi(j)] for all j = 1, . . . , k. Repeat this procedure for i = 1, . . . , l − 1.

We experimented with more sophisticated merging procedures such as a beam-search approach but did not observe a significantly better phasing on any of our data sets.

Phasing output Once a final partition P has been found, we build the final phased haplotypes. The output is k sequences in {0, 1, 2, 3, −}m where m is the number of variants. flopp can take fragment files as input, in which each line of the file describes a single read fragment in {0, 1, 2, 3, −}m, or it can take a VCF file and a BAM file of aligned reads. In the former case, without a VCF file, we do not have genotyping information, so we simply output {H1, . . . , Hk} where Hi = H(P [i]) is the consensus haplotype. If a VCF file is present, flopp constrains the final haplotype by the genotypes, that is, for some output variant, the number of reference and alternate alleles is the same as in the VCF file. We describe this procedure in Appendix D.

Parameter estimation We already mentioned in previous sections how we set all parameters for algorithms except for and the parameter σ in the UPEM score. We set σ to be the median length of the reads divided by 25, which we found empirically to be a good balance between the binomial and chi-squared terms.

To estimate , we start with an initial guess of = 0.03. We then select 10 subsets Bi ⊂ R at random and perform local clustering and refinement. We estimate from the 10 · k total clusters by choosing = the 110 th quantile error, D(C) where for each cluster C ⊂ Bi the error is S(C)+D(C) . We pick a bottom quantile because we assume that there is some error in our method, so to get the true error rate we must underestimate the computed errors. 3 Results and Discussion 3.1

There are a plethora of phasing metrics developed for diploid and polyploid phasing [ 26 ]. We use three different metrics of accuracy, but argue that each

We used the v4.04 assembly of Solanum tuberosum produced by the Potato Genome Sequencing Consortium [ 15 ] as a reference. We took the first 3.5 Mb of chromosome 1 and removed all “N”s, leaving about 3.02 Mb and simulated bi-allelic SNPs using the haplogenerator software [ 26 ], a script made specifically for realistic simulation of variants in polyploid genomes.

We generated SNPs with mean distance between SNPs of 45 bp. This is in line with the 42.5 bp average distance between variants as seen in [ 35 ] for Solanum tuberosum; in that study they observe that > 80% of variants are biallelic SNPs. This is also in line with the 58 bp mean distance between variants 12

J. Shaw and Y.W. Yu. seen for the hexaploid sweet potato (Ipomoea batatas) observed in [ 38 ]. The dosage probabilities provided to haplogenerator are the same parameters as used in [ 26 ] for the tetraploid case. When simulating triploid genomes, we use the same dosage probabilities but disallow the case for three alternate alleles. We repeated our test for mean distance between SNPs of 90 and 135 bp in Appendix G.

We used two different software packages for simulating reads. We used PaSS [ 39 ] with the provided default error profiles for simulating PacBio RS reads. PaSS has higher error rates than other methods like PBSIM [ 29 ] which tends to underestimate error [ 39 ]. We used NanoSim [ 37 ] for simulating nanopore reads using a default pre-trained error model based on real human dataset provided by the software.

After generating the haplotypes and simulating the reads from the haplotypes, we obtain a truth VCF from the haplotypes. We map the reads using minimap2 [ 20 ] to the reference genome. The scripts for the simulation pipeline can be found at https://github.com/bluenote-1577/flopp test. 3.3

We primarily test against H-PoPG [ 36 ], the genotype constrained version of HPoP. We tried testing against WHP, but found that the results for WHP were extremely discontiguous and the accuracy was relatively poor across our data sets. We show an example of this in Appendix F. Other methods such as HapTree and AltHap were tested, but we ran into issues with either computing time or poor accuracy due to the methods not being suited for long-read data. We did not test against nPhase because the output of nPhase does not have a fixed number of haplotypes.

The switch error and Hamming error rates are shown in Figure 3. For each test, we ran the entire pipeline three times; each run at high ploidies takes on the timescale of days to complete. The run times on PacBio reads for H-PoPG, flopp, as well as one instance of WHP on 3x ploidy data are shown in Figure 4.

For the nanopore data set simulated from NanoSim, the results for H-PoPG and flopp are very similar. The PaSS PacBio simulator outputs reads that are more erroneous, and we can see that flopp generally performs better than HPoPG across ploidies except for the Hamming error rate when the coverage is relatively low; interestingly, the switch error rate is still lower in this case. flopp’s switch error rate is consistently 1.5-2 times lower than H-PoPG for the simulated PacBio reads. On the low coverage data sets, H-PoPG’s global optimization strategy leads to a better global phasing than flopp’s local methodology.

Note that in these tests we phase a contig of length 3.02 Mb, whereas most genes of interest are < 50 kb. In Figure 5 we plot the mean q-block error rates of the 5x and 6x ploidy phasings at 10x coverage; these runs have higher Hamming error rates for flopp than H-PoPG. For blocks of up to around 850 · 45 = 38250 bases, flopp outputs phasings with lower q-block error rates than H-PoPG despite a larger global error rate. Although flopp may sometimes give worse global phasings than H-PoPG, flopp can give extremely accurate local phasings.

Computationally, Figure 4 shows that flopp is at least 3 times faster than H-PoPG and up to 30 times faster than H-PoPG for 6x ploidy on a single thread.