Integer programming framework for pangenome-based genome inference

Abstract

Affordable genotyping methods are essential in genomics. Commonly used genotyping methods primarily support single nucleotide variants and short indels but neglect structural variants. Additionally, accuracy of read alignments to a reference genome is unreliable in highly polymorphic and repetitive regions, further impacting genotyping performance. Recent works highlight the advantage of haplotype-resolved pangenome graphs in addressing these challenges. Building on these developments, we propose a rigorous alignment-free genotyping framework. Our formulation seeks a path through the pangenome graph that maximizes the matches between the path and substrings of sequencing reads (e.g., k-mers) while minimizing recombination events (haplotype switches) along the path. We prove that this problem is NP-Hard and develop efficient integer-programming solutions. We benchmarked the algorithm using downsampled short-read datasets from homozygous human cell lines with coverage ranging from 0.1× to 10×. Our algorithm accurately estimates complete major histocompatibility complex (MHC) haplotype sequences with small edit distances from the ground-truth sequences, providing a significant advantage over existing methods on low-coverage inputs. Although our algorithm is designed for haploid samples, we discuss future extensions to diploid samples.

Fig. 1:

A small example of our reduction from Hamiltonian Path Problem to Problem 1 (Theorem 1). (Top) The starting instance of $G$ of Hamiltonian Path Problem. (Bottom) The vertex labeled graph $G^{\prime }$ constructed from $G$. Here, $n=4$ and we assume $c=2$, making $b=\langle {\log }_2(n+2(c(n+1)+\text{1))}\rangle +1=6$. Each edge is supported by a unique haplotype (not shown). The string set is $𝓢=\{0000010000001,0000100000001,\dots ,0110100000001\}$.

Fig. 2:

Evaluation of the performance of the IQP method with and without relaxation of the binary edge variables $x_{uv}$. We compared runtime using various short-read datasets.

Fig. 3:

Evaluation of the performance of the ILP method with and without relaxation of the binary edge variables $x_{uv}$. We compared runtime using various short-read datasets.

Fig. 1:

A simple illustration of an haplotype-resolved pangenome graph with two haplotype paths highlighted in pink and blue colors. An inferred path with a single recombination is shown as a dashed line.

Fig. 2:

(A) A pangenome graph with four haplotype paths $h_1,h_2,h_3$ and $h_4$. Set of haplotype paths passing through a vertex is listed below each vertex. (B) The corresponding expanded graph which includes four disjoint paths, one for each haplotype path. The recombination edges are shown in purple, these edges have a weight of $c$. We consider only the useful recombinations (Lemma 3). The edges which are not recombination edges in the expanded graph have a weight of 0. (C) The corresponding optimized expanded graph.

Fig. 3:

Accuracy of haplotype sequences estimated by PHI, VG and PanGenie using short reads from MHC sequences of five haplotypes (APD, DBB, MANN, QBL, SSTO). The x-axes indicate the coverage of short-read data. The y-axes indicate the edit distance between the estimate haplotype sequence and the ground-truth sequence on a logarithmic scale.

Fig. 4:

Performance comparison between the ILP and IQP solutions implemented in PHI. We compared their runtime and memory-usage using short-read sequencing datasets sampled from five haplotypes.

Fig. 5:

Assessement of PHI’s performance with the increasing number of genomes in pangenome graph. The left figure shows the accuracy in terms of edit distance between the output sequences and ground-truth sequences. The middle and right figure show the runtime and memory-usage respectively.