Phasing Data from Genotype Queries via the $\mu$-PBWT

Within the field of haplotype analysis, the Positional Burrows-Wheeler Transform (PBWT) [9] has emerged as a fundamental data structure for efficiently encoding and searching large-scale haplotype data. Its ability to rapidly detect long shared substrings between a query sequence and a reference panel makes it a crucial tool for various genomic applications, including haplotype imputation [19, 18], ancestral inference [13, 25], and genotype phasing [15]. One of the key advantages of PBWT is its efficient representation of haplotype panels, allowing fast and memory-efficient computations, particularly compared to traditional sequence alignment-based approaches. Over the years, PBWT has been extended and refined to support more complex operations [24, 26, 17, 21], making it an essential component of modern haplotype analysis pipelines.

Sanaullah et al. [22], who formulated the haplotype threading problem, the task of representing a query haplotype using a minimal set of substrings derived from a reference panel, introduced a significant development in PBWT-based haplotype analysis. To address this problem, they introduced the Minimal Positional Substring Cover (MPSC) problem, which seeks to identify the smallest collection of substrings that fully cover a query while maintaining positional consistency with the reference haplotypes. Their work proposed an efficient PBWT-based solution by modeling the solution space as a graph, enabling computationally efficient haplotype threading. However, their approach focused primarily on haplotypes rather than genotypes, leaving an open question as to how these methods could be extended to genotype phasing, which requires handling heterozygous loci and potential recombination events.

Building on this work, in 2024, Bonizzoni et al. [3] demonstrated that the MPSC problem could be solved in sublinear space, significantly improving its scalability for large reference panels. That approach was based on efficiently computing $k$-Set Maximal Exact Matches (SMEMs) via $k$-Matching Statistics ($k$-MS), a technique that enabled the identification of minimal positional substring covers with reduced memory requirements. This advancement allowed for an efficient solution to the MPSC problem and its variations proposed by Sanaullah et al. The ability to solve these problems in sublinear space is particularly valuable for large-scale genomic datasets, where traditional PBWT-based methods may become memory-intensive.

In this paper, we extend the results of Bonizzoni et al. [3] to the problem of genotype phasing, adapting the MPSC framework to work with genotype queries instead of haplotypes. The main idea that inspires this work is the observation that large reference panels likely include a haplotype pair that explains a genotype query. However, we need methods that can scale to such large panels. Here, we introduce a novel approach using $\mu$-PBWT [8], a refined version of run-length encoded PBWT (RLPBWT) [2], to efficiently phase genotypes while optimizing for minimal recombination events. Our method is designed to reconstruct haplotype pairs that explain a given genotype by leveraging positional substrings (PSs) and complementary positional substrings (CPSs) extracted from a reference panel. Unlike previously proposed combinatorial methods that assume a perfect haplotype match [14, 4, 1], our approach can handle cases where recombinations are required to fully explain the genotype, making it more adaptable to real-world data. In addition, current state-of-the-art tools such as Beagle [5, 6] are based on the Li and Stephens Hidden Markov Model [16], while our approach is deterministic in the absence of mutation and recombination.

Our study focuses on two primary scenarios: (1) an ideal case where the reference panel contains two haplotypes that fully explain the genotype without requiring recombination, and (2) a more realistic setting where recombinations across haplotype segments of the panel are needed due to the absence of a single haplotype pair that explains the genotype query. To address these challenges, we introduce a combinatorial approach that processes 2-positional substrings (2PSs) and identifies valid haplotype pairs that minimize recombination events. Furthermore, we implement and evaluate our method against Beagle [6, 5], a state-of-the-art statistical phasing tool. Our experimental results demonstrate that $\mu$-PBWT achieves competitive accuracy while significantly reducing computational costs in scenarios with low mutation rates. On the other hand, introducing target mutations (and consequently recombinations not having the target parents in the reference panel) reduces the performance of our approach, in terms of switch error rate/mismatch rate. Our method requires little memory and can scale on very large reference panels, thanks to the $\mu$-PBWT index. This makes $\mu$-PBWT a promising alternative for applications where memory efficiency and processing speed are critical.

Due to the diploid nature of humans, chromosomes come in two copies, called haplotypes, each copy inherited from one of the two parents (see [4] for a survey of the main concepts). Sequencing an individual is a process that does not distinguish between the two haplotypes, as reads are fragments uniformly distributed over both chromosomal copies. Distinguishing the haplotypes from sequencing reads is biologically expensive, and thus computational methods are employed for this purpose.

A genotype is simply the conflated information of a pair of haplotypes and can be easily determined by sequencing. More precisely, a genotype is specified by a vector of allele pairs at each locus, where the pairs consist of homologous or heterozygous alleles. In a genotype, information on which heterozygous allele belongs to which haplotype is not available and must be determined.

Observe that each haplotype inherited from a parent is itself the result of a mosaic of the two parental haplotypes, due to recombination events, which are a typical evolutionary phenomenon occurring in human populations. Recombination events regulate Mendelian inheritance laws within a family trio and have been extensively studied in the context of combinatorial methods for the phasing or haplotyping of input genotype matrices [4].

From a combinatorial perspective, a haplotype at a set of biallelic loci, the genomic positions in which only two distinct alleles are observed in a population, is formalized as a binary string. In this representation, the $i$-th position of the string is assigned a 0 if the individual carries the reference allele at that locus (i.e. the allele matches the reference in that locus), and a 1 if the individual carries the alternative allele.

The genotype of an individual is the combined representation of the two haplotypes inherited from the two parents. Formally, a genotype is represented as a string on the alphabet $\{0,1,2\}$, where the $i$-th position is $0$ if the individual has inherited the reference allele from both parents at the $i$-th locus, it is $1$ if the individual has inherited the alternative allele from both parents at the $i$-th locus, or it is $2$ if the individual, at the $i$-th locus, has inherited the reference allele from a parent and the alternative allele from the other parent. The positions where the genotype is $0$ or $1$ are called homozygous, while the others are called heterozygous. We say that a pair $⟨h_1,h_2⟩$ of $w$-length haplotypes explains a $w$-length genotype $Q$ if and only if, for all $i\in [1,w]$ we have that $h_1[i]=h_2[i]=Q[i]$ if $Q[i]\ne 2$, or $h_1[i]\ne h_2[i]$ if $Q[i]=2$. In simpler terms, a pair of haplotypes explains a genotype if they match the genotype at homozygous positions and differ at heterozygous positions.

Then, the genotype phasing problem aims at reconstructing an individual’s haplotypes from their genotype. More specifically, determining the exact pair of haplotypes that correspond to the observed genotype. The problem is straightforward for homozygous positions, where both haplotypes match the genotype; however, at heterozygous positions, the phase is ambiguous because multiple pairs of haplotypes, differing only in the assignment of alleles to each parent, can equally explain the genotype. There is a vast literature on combinatorial methods for genotype phasing, mainly distinguished by how the problem is modeled and by the type of input data considered. For example, Gusfield proposed solving genotype phasing using the perfect phylogeny model [14, 1] when the input matrix consists of a collection of genotype vectors.

The availability of large panels of population haplotype data has shifted the focus of genotype phasing toward the use of statistical methods which, based on information from an input panel, can effectively explain genotype queries. An example of this is Beagle [5], which uses phased input data, that is, a known haplotype panel, to build a hidden Markov model specifically tailored to diploid data.

In this paper, genotype phasing is addressed through the so-called haplotype threading, which models the process of reconstructing a haplotype as a mosaic of haplotypes from a reference panel.

More precisely, haplotype threading aims to cover a haplotype query using segments from haplotypes in an input panel. The main idea is that a haplotype results from the accumulation of recombination events, combining chromosomal segments from various haplotypes in the panel. The Li and Stephens (LS) hidden Markov model (HMM) [16] produces a haplotype threading for a query haplotype using a panel of known haplotypes, requiring time that is linear in the size of the input panel. However, this running time becomes infeasible when the haplotype panel is large, and thus sublinear-time solutions have recently been investigated, mainly using the PBWT (Positional Burrows-Wheeler Transform) framework.

Moreover, modern haplotype panels typically exhibit much lower mismatch rates and switch error rates, two key measures commonly used in the probabilistic modeling of the LS framework. Observe that a zero mismatch rate between the query haplotype and the input panel corresponds to the situation where the query can be entirely covered by segments from the panel, with the number of such segments directly corresponding to the switch error. In this case, combinatorial approaches become particularly attractive, especially when they offer faster solutions than LS-based methods.

In particular, the solution proposed by Sanaullah et al. [20], the Minimum Positional Substring Cover (MPSC) with non-overlapping segments, corresponds exactly to haplotype threading under a zero mismatch error assumption. When overlapping segments are permitted, identifying the breakpoint between two matching positional substrings corresponds to determining the transition point between two distinct haplotypes in the panel.

Now, the input panel is not guaranteed to include the two haplotypes from which a query genotype may have been generated, nor do we know whether the two haplotypes of a parent are present in the panel when explaining a haplotype as the result of recombinations from that pair. This is why the mosaic of segments may involve multiple haplotypes of the panel.

The introduction of the PBWT has significantly improved the time efficiency of computational approaches to haplotype threading. The first major improvement concerns the time required to find coverage of a haplotype, which has been reduced from $O(N\cdot M)$ to $O(N)$, where $N$ is the number of sites in a haplotype sequence [20]. This reduction in time complexity allows the methods to scale to a large number $M$ of haplotype sequences, which was not feasible with the original LS-based method. Finally, since the memory required by these approaches can be extremely prohibitive, a further improvement was proposed in [3], where a compressed data structure is introduced to store and manage the input panel using the $\mu$-PBWT [8]. In particular, space efficiency is crucial for handling large cohorts of haplotype data, helping to reduce mismatch errors and switch errors. In this paper, we will assume an ideal future framework where the large amount of data allows for the assumption that the two haplotypes that explain a genotype are part of the input data.

In this section, we extend the definition of positional substring cover to the case of a genotype and a panel $S$ of strings. To accomplish this, we first define the concept of a complementary positional substring (CPS) for patterns $P$ and $P^{\prime }$. A CPS consists of two positional substrings $(i,j,P)$ and $(i,j,P^{\prime })$, which are bitwise complementary. For simplicity, a complementary positional substring will be simply identified by a triple $(i,j,P)$.

Building on this concept, we let $Q$ be a $w$-long genotype query and $S$ a set of $w$-long strings, then a PSC_g of genotype $Q$ consists of a set $C$ of positional substrings $(i,j,Q)$ that occur in at least two rows of panel $S$, i.e., 2PS, and a CPS $(k,l,S)$ of $S$, such that: (i) $(k,l,Q)$ is a substring of $Q$ of only 2’s, and (ii) all positions $t$ of $Q$ are covered by $C$.

Next, we formally state the problem investigated in this paper. The core objective is to phase the query genotype $Q$ by covering it with positional substrings (PSs) or complementary positional substrings (CPSs) derived from a panel of haplotypes $S$.

In Problem 1, a natural optimization goal is to minimize $\mathrm{\ell }$, the size of the PSC_g. However, to further reduce the number of recombinations, we focus on minimizing the number of distinct pairs of haplotype sequences from $S$ that contribute to the positional substrings used to cover $Q$. In other words, our goal is to minimize the number of 2-PSs and CPSs that originate from different haplotype pairs in $S$. We refer to this variant of Problem 1 as minimum recombination genotype-positional phasing. The problem has no recombinations if all 2-PSs and CPSs are derived from the same pair of haplotypes. This criterion will be incorporated into the heuristic method we propose to solve the genotype-positional phasing problem. The primary case we investigate in this paper is when the two haplotypes that resolve $Q$ are already present in the panel $S$, which means that no recombinations occur. As we will discuss in Section 5, the case of no recombination has a straightforward and efficient solution. Finally, we recall that each 2PS is associated with a set of rows in the reference panel, and we will explicitly use this set to solve Problem 1.

In this section, we present our approach to solving Problem 1 using $\mu$-PBWT. From this point forward, we make the simplifying assumption that each column in the reference panel contains at least two 0’s and two 1’s. This ensures that every homozygous position in the query genotype shares a 2PS with the reference panel, simplifying the presentation of our method. To effectively apply $\mu$-PBWT, we first refine how it identifies matches and 2-positional substrings (2PSs) shared between the target genotype and the haplotype panel. Specifically, we introduce a constraint that limits the length of these 2PSs to at most two. The rationale behind this constraint will be explained at the end of this section, where we demonstrate how the use of 2-MPSCs affects the number of phased heterozygous sites.

In this section, we present some preliminary experimental results of our approach. Since this approach is a proof-of-concept, the current implementation is inefficient for most subtasks, from the extraction of complementary rows to the greedy algorithm used to select the best pairs of haplotypes. All of these steps could be further optimized in the future, both theoretically and practically.

To assess the applicability of our theoretical results, we integrate the genotype phasing algorithm into the $\mu$-PBWT codebase and tested it in a simple genotype phasing scenario. For this experiment, we considered the HGSVC2 reference panel [7, 10] for human chromosome 2, with 68 haplotypes (34 samples) and $\mathrm{229\hspace{0.17em}300}$ biallelic variation sites. We randomly selected a sample (target sample) from this reference panel as a query, combining its haplotypes into an unphased genotype target. The target sample we selected has 60,604 heterozygous sites ($\sim 26\%$ of the total number of variants on chromosome 2).

We used VCF as input format where $0|0$ and $1|1$ represent phased homozygous sites, $0|1$ and $1|0$ phased heterozygous sites, $0/0$ and $1/1$ unphased homozygous sites, and $0/1$ unphased heterozygous sites.

We compared $\mu$-PBWT with Beagle [5, 6], a state-of-the-art phasing tool that is also based on the use of a reference panel. We ran Beagle with default parameters, using the HapMap [23] chromosome 2 genetic map for the GRCh38 reference genome.

We ran our experiments on a machine equipped with an Intel Xeon CPU E5-4610 v2 (2.30GHz), 256 GB of RAM, 8 GB of swap, and Ubuntu 20.04.6 LTS 64-bit with kernel 5.15.0. We measured execution times and peak memory usage using the Unix time command.

Moreover, we simulate mutations in the target genotype by replacing the genotype at each position with probability $\epsilon$. We considered five different values for $\epsilon$: 1%, 3%, 5%, 10%, and 20%. We refer to $\epsilon$ as mutation rate.

Table 1 summarizes the total number of heterozygous sites after the mutation step. As explained above, mutations cause small spurious recombinations at nearby sites. Thus, this small experiment allows us to test more realistic scenarios.

The accuracy of the results was evaluated by considering (i) the switch error rate, i.e., the rate of phase changes between the predicted haplotypes w.r.t. the ground truth haplotypes, (ii) the mismatch rate, i.e., the rate at which two switch errors occur at consecutive positions, and (iii) the percentage of heterozygous sites that have been phased. The three metrics were computed using the utility script publicly shared by the developers of HapCUT2 [11]¹¹1https://github.com/vibansal/HapCUT2/tree/master/utilities, considering the two haplotypes of the target sample as ground truth.

Figure 3 presents the accuracy obtained by Beagle and $\mu$-PBWT at different mutation rates (detailed results are available in Table 2). The results obtained in the scenario without mutations (i.e., with a mutation rate equal to 0%) are consistent with how the two tools approach the solution of the phasing problem. Beagle is a statistical approach that is not expected to find the exact solution even if the target genotype can be perfectly explained by two haplotypes in the reference panel. In contrast, $\mu$-PBWT implements a combinatorial approach which is designed to always find a haplotype pair that explains the target, if such a pair exists in the reference panel. In fact, in this scenario $\mu$-PBWT correctly identifies the two haplotypes of the target sample, so both switch error rate and mismatch rate are equal to 0.0%, while Beagle exhibits some fluctuations, with both metrics greater than zero.

As the mutation rate increases, the accuracy of $\mu$-PBWT decreases dramatically. In fact, as we can see in Figure 3c, even with a mutation rate of 5%, $\mu$-PBWT is able to phase only 62.6% of heterozygous sites. We recall that $\mu$-PBWT does not phase heterozygous sites that are not explained by any haplotype pair that can be used as both a left and a right anchor. As a consequence, the higher the mutation rate, the larger the heterozygous genotype regions that do not have valid anchors and thus remain unphased. This fact further demonstrates the dependency of our approach on the data available in the reference panel, which benefits from panels including a large number of haplotypes. In addition to this strong dependency on the reference panel, our approach is unable to estimate haplotype pairs that differ from the ones obtained using the positional cover. Thus, around each mutation site, $\mu$-PBWT will report pairs of haplotypes that differ from the ones of the target sample. However, a manual inspection of the results revealed that one of the two reported haplotypes often coincides with a haplotype of the target sample. However, these cases are handled better by Beagle, again due to its probabilistic approach, which seems to be more robust against these mutations and successfully phases each genotype site.

Figure 4 presents the computational resources used by the two approaches. The running times of the two approaches are comparable. However, the running time of $\mu$-PBWT increases as the mutation rate increases, since, as expected, the number of haplotype switches that $\mu$-PBWT has to compute (instead of simply continuing the match) also increases. Each time the algorithm cannot proceed with any pair, it has to start again, considering larger sets of possible pairs whose analysis increases both the overall execution time and RAM usage. Note that enlarging the number of reference samples would worsen performance in terms of indexing and storing the reference panel, but it would likely reduce the number of possible pairs to consider at each iteration since longer matches would become more likely. On the other hand, the running time of Beagle remains constant as the mutation rate increases since it does not depend on the number of haplotype switches in the reported solution.

As expected, $\mu$-PBWT is very memory-efficient, requiring approximately 7 times less memory than Beagle. These results show that $\mu$-PBWT could be run on low-end hardware or, on the public cloud, at a fraction of the cost of running Beagle. We argue that further refinements and improvements to the model and algorithm would not significantly impact memory requirements, since memory usage in genotype phasing approaches is often dominated by the space required to store the reference panel and $\mu$-PBWT employs a very compact and efficient index representation of the panel.

In this paper, we present a proof of concept showing how $\mu$-PBWT can be used to solve the genotype phasing problem. First, we outline the theoretical framework underlying our method, based on the foundations presented in [3]. We then introduce our combinatorial approach to tackling the phasing problem, considering both scenarios: (i) when a haplotype pair fully explains the target genotype and (ii) when multiple haplotype pairs are required to cover different genotype regions. Experimental results indicate that our approach is sensitive to the size of the reference panel, particularly when the mutation rate increases.

Improving the resolution of haplotype reconstruction in complex genomic regions remains an important direction for future work. In particular, anchoring heterozygous regions using flanking homozygous segments may enhance stability during phasing, while incorporating methods to detect recombination events within heterozygous regions could improve accuracy in the presence of structural variation and mutation.

These enhancements build on the current strengths of the proposed method, which, in mutation-free settings, reliably identifies the correct haplotype pair. Furthermore, consistent with the evaluations of $\mu$-PBWT, the method achieves high efficiency with minimal memory requirements, making it suitable for execution on commodity hardware. Taken together, these characteristics provide a strong foundation for future comparisons with state-of-the-art phasing tools, particularly in real-world scenarios involving complex variant patterns and sequencing noise.