# Efficient and effective methods for variant selection

Srinivas Aluru, Anne Benoit, Bora Uçar

► To cite this version:

Srinivas Aluru, Anne Benoit, Bora Uçar. Efficient and effective methods for variant selection. 16th ACM Conference on Bioinformatics, Computational Biology, and Health Informatics (BCB 2025), Oct 2025, Philadelphie, PA, United States. hal-05237949

HAL Id: hal-05237949

https://hal.science/hal-05237949v1

Submitted on 11 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Efficient and effective methods for
variant selection

Srinivas Aluru¹, Anne Benoit², and Bora Uçar³

¹IDEaS, Georgia Institute of Technology, Atlanta, USA, aluru@cc.gatech.edu

²LIP, ENS Lyon & Institut Universitaire de France & IDEaS, Georgia Institute of Technology,
Atlanta, USA, anne.benoit@ens-lyon.fr

³CNRS, LIP (UMR5668, CNRS - ENS Lyon - UCB Lyon 1 - Inria), and IDEaS, GaTech,
Atlanta, USA, bora.ucar@ens-lyon.fr

October 11, 2025

Abstract

Variation graphs succinctly capture genetic variations among individuals within a species or a target population. The use of variation graphs instead of a single reference genome is credited with reducing bias and increasing the accuracy of sequence mapping algorithms. However, complete variation graphs that comprehensively incorporate all genetic variations are often found to be ineffective and inaccurate in practice due to the presence of a combinatorially explosive number of paths in the graph that do not correspond to any observed genome. Thus, a balance is struck in carefully selecting a subset of variants to be incorporated, for which mathematical frameworks have recently been developed. We advance the mathematical framework proposed by Jain et al.[Bioinformatics, 37 (2021)], where integer linear programming formulations were developed for optimal variant selection. We propose novel graph-based formulations and develop exact and fast algorithms for certain cases, approximation methods for some others, and empirically close to optimal results in all cases. The primary advantage of the algorithms designed here is that they provide near-optimal results at orders of magnitude faster run time of an ILP solver.

# 1 Introduction

Rapid decline in sequencing costs caused a paradigm shift away from merely sequencing a representative genome per species to the current practice of sequencing numerous individuals to catalog the underlying genetic diversity. Such efforts are needed to understand the complex influence of genetic traits in determining phenotypic characteristics, and are vital to diverse applications from engineering plants for yield or tolerance of specific environmental conditions to understanding disease propensity and delivering the promise of precision medicine for humans.

Graphs are a natural and compact representation for multiple related genomes due to the overwhelming identity among them, with variations captured as short branching paths. The term *genome graph* is used when representing any set of genomes/haplotypes (some organisms have multiple haplotypes, the variants among which confer a degree of robustness) not necessarily from the same species, although such a graph is typically used with genomes exhibiting high levels of similarity. *Variation graphs* are a specific subclass of genome graphs, when the graphs are used to capture variations among individuals within a species or a target population. In what follows, we use the term variation graphs to more explicitly call out the focus of this paper, *variant selection*, though everything in our work, as well as other works cited here, is broadly applicable.

In a graph, each constituent haplotype is present as a path. Unavoidably, there are graph paths that represent components from different haplotypes artificially stitched together, and are not part of any observed haplotypes. The combinatorial explosion in the number of such paths with an increase in the number of variants incorporated is considered the primary reason for diminishing returns and eventual decline in accuracy as more and more variants are incorporated. In addition, the run time performance of sequence-to-graph algorithms
degrades with an increase in graph complexity. Thus, variation graphs are built taking only a subset of known variations into account, the selection of which is termed the *variant selection problem*.

Initial approaches to Variant selection used ad hoc methods, such as incorporating variants based on their frequency of occurrence [9], or targeting variants relevant to a particular population or the context of a biological study [6, 7]. More recently, mathematical frameworks with provable guarantees that perform well across a range of scenarios have been developed [14, 17]. In this work, we consider the variant selection framework proposed by Jain et al. [14]. Given two parameters $\alpha$ and $\delta$, this framework seeks to minimize the number of variants selected subject to retaining substrings of length up to $\alpha$ with at most $\delta$ differences. The types of variants considered (SNPs, indels), the objective function (minimizing the number of variants or the number of variant loci), and whether the substrings to be preserved are from the haplotypes or more inclusively paths from the complete variation graph, result in a rich set of combinatorial problems. Most of these problems are NP-hard and are solved using Integer Linear Programming (ILP) approaches [14, 21], although certain problems only involving SNPs are polynomially solvable [14].

We revisit four variant selection problems under Jain et al.'s framework in which the substrings to be preserved are from the variation graph. Two of these problems are solvable in polynomial time, and the other two are NP-hard. Our main contributions are as follows. We provide a novel graph-theoretic formulation using the well-known class of interval graphs [2], that is applicable to all four variant selection problems. We develop efficient and effective heuristics for one of the polynomial-time solvable problems; Our heuristics run in linear or near-linear time, while existing algorithms have super-quadratic complexity. We show that our heuristics solve some of the sub-cases exactly. One of the heuristics is shown to have an approximation guarantee of 1/2, and others behave better empirically in terms of run time or variation in the solution quality. We also adapt two of our heuristics to the two NP-hard problems. Their performance is analyzed empirically and is demonstrated to be close to optimal results. While the problems are NP-hard, a commercial-strength ILP solver solves most of the problem instances in under four hours. In contrast, our best performing heuristic runs in under 19 seconds for all instances of the two NP-hard problems in our experiments.

The rest of the paper is organized as follows: We first provide essential background on Jain et al.'s variant selection framework [14] and the graph fundamentals we rely on (Section 2). Our core algorithmic contributions are described in Section 3, where we propose graph-theoretic modeling, and subsequently design several algorithms to solve the different variant selection problem scenarios. Extensive experimental results are presented in Section 4, both for the simpler case with only SNP variations, and when including insertions and deletions (indels). Finally, we conclude and discuss future research directions in Section 5.

# 2 Background and problems studied

## 2.1 Variant reduction framework

We summarize the variant reduction framework by Jain et al. [14]. Using their notation, let $R_1, R_2, \dots, R_m$ be $m$ input haplotype sequences, among which $R_1$ is designated as the reference haplotype and all others are described as variations from $R_1$. The whole input is represented by an edge-labeled directed multigraph $G = (V, E, \sigma)$, with the vertex set $V$, edge set $E$, and edge-labels $\sigma$. The edge labels are in $\Sigma \cup \{\varepsilon\}$, where $\Sigma = \{A, C, G, T\}$ is the alphabet of DNA bases and $\varepsilon$ is the empty character. First, $R_1$ is represented as a directed chain on $|R_1| + 1$ vertices, whose edge-labels read as the sequence $R_1$. This chain is called coordinate axis for the whole graph. Then, this chain is augmented with variants present in other $R_i$'s. An SNP variant adds a labeled edge between two adjacent coordinates. A deletion variant adds an edge with the label $\varepsilon$ between two coordinates separated by the deletion length. An insertion variant is represented by a chain of one or more labeled edges starting and ending at the same coordinate. The number of variants at a coordinate $i$ is equal to the out-degree of the vertex $v_i$ minus one.

A path in $G = (V, E, \sigma)$ with $\alpha$ edges denotes a string of length $\alpha$. A path is restricted to visit a vertex at most once, unless there is an insertion variant at a visited vertex, in which case that vertex can be visited twice if the path includes the variant. Therefore, a path is allowed only one insertion variation at a vertex. The graph $G = (V, E, \sigma)$ contains paths that are substrings of the original haplotypes, and other genotypes that are not necessarily in the input. The variant reduction framework [14, 21] aims to reduce the graph $G = (V, E, \sigma)$ to $G' = (V', E', \sigma')$ by selectively removing edges and keeping the labels of the retained edges intact. The constraint is to ensure that all $\alpha$-length strings corresponding to a path in $G$ are mapped to the same starting vertex in $G'$ with at most a user-specified error threshold $\delta$. The error threshold is defined by Hamming or edit distance, which apply, respectively, to the SNP-only and SNP-and-indels cases. For a given variation graph $G$, a reduced
graph $\mathcal{G}'$ is $(\alpha, \delta)_h$-compatible if all $\alpha$-long strings in $\mathcal{G}$ can be mapped to their corresponding paths in graph $\mathcal{G}'$ with Hamming distance at most $\delta$—that is, the number of positions that differ from an $\alpha$-length path in $\mathcal{G}$ and its corresponding path in $\mathcal{G}'$ is less than $\delta$. Similarly, $\mathcal{G}'$ is $(\alpha, \delta)_e$-compatible if all $\alpha$-length strings in $\mathcal{G}$ can be mapped to the corresponding paths in $\mathcal{G}'$ with an edit distance at most $\delta$—that is, the number of insertions, deletions, and substitutions required to change an $\alpha$-long path in $\mathcal{G}$ to its corresponding path in $\mathcal{G}'$ is less than $\delta$.

Jain et al. define the following four problems in their variant reduction formulation:

**Problem 1 :** Given a variation graph $\mathcal{G}$ where all variations are SNPs, find an $(\alpha, \delta)_h$-compatible reduced graph $\mathcal{G}'$ with the minimum number of coordinates containing one or more variants.

**Problem 2 :** Given a variation graph $\mathcal{G}$ where all variations are SNPs, find an $(\alpha, \delta)_h$-compatible reduced graph $\mathcal{G}'$ with the minimum number of variants.

**Problem 3 :** Given a variation graph $\mathcal{G}$, find an $(\alpha, \delta)_e$-compatible reduced graph $\mathcal{G}'$ with the minimum number of coordinates containing one or more variants.

**Problem 4 :** Given a variation graph $\mathcal{G}$, find an $(\alpha, \delta)_e$-compatible reduced graph $\mathcal{G}'$ with the minimum number of variants.

The first two problems are polynomial-time solvable, whereas the decision versions of the last two are NP-hard, because of the added complexity of considering insertion and deletion variants [21]. For the first problem, Jain et al. give a linear time exact algorithm. For the second one, they formulate an integer linear program (ILP), and recognize that the constraint matrix is totally unimodular, and hence an LP solver finds an optimal solution. For the third and fourth problems, they propose a greedy algorithm for a heuristic solution and an ILP formulation for an exact solution. Our aim in this paper is to develop fast algorithms that deliver good results for all four problems.

There are also haplotype-aware versions of the four problems [21], in which the compatibility is defined with respect to the given haplotypes, not all possible strings from the variation graph. Haplotype-aware versions are beyond the scope of our study.

## 2.2 Graph preliminaries

We use undirected graphs to tackle the four problems, with a total ordering on the vertices. For a graph $G = (V, E)$ with $n = |V|$ vertices, we use $\pi(\cdot)$ to denote the ordering, where $\pi(v)$ gives the order of the vertex $v$. We use $\text{adj}(v)$ to denote the neighbors of the vertex $v$. We use $\text{hadj}(v)$ to denote the neighbors of the vertex $v$ that are numbered higher than $v$ in the ordering, i.e., $\text{hadj}(v) = \{u : u \in \text{adj}(v) \text{ and } \pi(v) < \pi(u)\}$. We also use $\text{ladj}(v)$ to denote the neighbors of the vertex $v$ that are numbered lower than $v$ in the ordering. We use $\text{hadj}^+(v)$ to denote the set $\{v\} \cup \text{hadj}(v)$, and similarly $\text{ladj}^+(v) = \{v\} \cup \text{ladj}(v)$.

A chordal graph is a graph $G = (V, E)$, where for every cycle of length at least four, $G$ has an edge joining two nonconsecutive vertices of the cycle [11, Ch. 4]. Chordal graphs have a special type of vertex orderings, called perfect elimination orderings. A perfect elimination ordering for a graph with $v$ vertices and $m$ edges can be found in $\Theta(n+m)$ time (see [11, Ch. 4] or [20]). An ordering $\pi(\cdot)$ of the vertices of a chordal graph $G = (V, E)$ is a perfect elimination ordering whenever $\text{hadj}^+(v)$ is a clique for any vertex $v \in V$. All maximal cliques of a chordal graph are of the form $\text{hadj}^+(v)$ for a vertex $v$. These two properties are the key to solve the maximum clique and the maximum independent set problems for chordal graphs in linear time.

An undirected graph $G = (V, E)$ is called an interval graph if its vertices can be represented by intervals in the real line in such a way that $(u, v) \in E$ whenever the corresponding intervals of $u$ and $v$ overlap. Interval graphs form a proper subset of chordal graphs.

The vertices in our graph models have weights and another attribute that we call charge. We use $G = (V, E, w, c)$ to denote a graph with vertex weights $w$ and charges $c$. We use $w(v)$ and $c(v)$ to denote the weight and charge of a vertex $v$, and this notation is extended to a set of $S$ vertices as $w(S) = \sum_{v \in S} w(s)$, or $c(S) = \sum_{v \in S} c(s)$. We abuse this notation in algorithms, where $w(S) \leftarrow 0$ means “set $w(v) = 0$ for all $v \in S$”.

Matrices are represented with bold capital letters, as in $\mathbf{A}$, whose entries are represented with the corresponding lower case letters, as in $a_{ij}$. Vectors are represented with bold lower case letters, as in $\mathbf{x}$, whose entries are represented with the corresponding lower case letters, as in $x_i$.
# 3 Proposed Algorithms

Jain et al. use the following definitions to obtain an integer linear programming formulation to tackle Problems 2, 3 and 4. This formulation also covers Problem 1.

Let $p_1, \dots, p_n$ be the coordinates of the sites at which there are variants, in the ascending order. Let $\mathbf{A}$ be an $n \times n$ matrix such that $a_{ij} \neq 0$ iff the coordinate $p_j$ is in the range $(p_i - \alpha, p_i]$. Observe that $\mathbf{A}$ is a lower triangular matrix. In Problems 1 and 2, all nonzeros of $\mathbf{A}$ are equal to 1. In the other two problems, a nonzero $a_{ij} = c_j$, where $c_j$ is an estimated upper bound on the edit distance if the insertion and deletion variants are removed at site $j$ and does not depend on $p_i$. That is, in Problems 3 and 4, all nonzeros in the $j$th column of $\mathbf{A}$ have the same value $c_j$. The value $c_j$ cannot be determined prior to the construction of the reduced variation graph, as $c_j$ measures the edit distance between a path in $\mathcal{G}$ and its corresponding path $\mathcal{G}'$. That is why Jain et al. estimate $c_j$; we use their framework to use the same estimates. When there are only SNPs at the coordinate $p_j$ in Problems 3 and 4, we have $c_j = 1$, and the nonzeros in the $j$th column of $\mathbf{A}$ all are equal to one. Let us collect all $n$ estimated edit distances in a vector $\mathbf{c}$ for later use. Let us define $n$ boolean variables, $x_1, \dots, x_n$ associated with the coordinates $p_1, \dots, p_n$, where $x_i = 1$ if the variants at coordinate $p_i$ are removed. The constraints of all four problems can therefore be written as $\mathbf{Ax} \le \vec{\delta}$, where $\vec{\delta}$ is an $n$-vector whose components are $\delta$.

The objective function is $\max \mathbf{w}^T \mathbf{x}$, where $w_i$ is equal to the reduction if we remove variants at $p_i$. In Problems 1 and 3, $w_i = 1$ for all $i$. In Problems 2 and 4, $w_i$ is equal to the number of variants at $p_i$, which is equal to one less than the out-degree of the vertex $p_i$ in the variation graph $\mathcal{G}$. Jain et al. [14] recognize the constraint matrix to be totally unimodular for problems 1 and 2. That is why an LP solver, or a specific method based on minimum-cost flow solvers [8] are usable for solving these two problems. These are however costly methods—in the matrix multiplication time [4] for LP solvers or $O(mn \log n + m^2 \log n)$ time for the minimum-cost flow solvers [1, Ch. 10] on a graph with $m$ edges and $n$ vertices; while there are near linear time solvers for the minimum cost flow problem [3], these are theoretical tools (rather than practical methods).

We model the four problems as graphs with weights and charges assigned to vertices, aiming to develop fast solvers with both provable guarantees for certain cases and practical effectiveness in general. Let $G_\alpha = (V, E, w, c)$ be a graph with $n$ vertices such that $v_i$ corresponds to $p_i$. There is an edge $(v_i, v_j) \in E$, whenever $p_j$ is in $(p_i - \alpha, p_i]$, and $w(v_i)$ and $c(v_i)$ are the weight and charge of the vertex $v_i$. $G_\alpha$ is an interval graph as $v_i$ corresponds to the interval $(p_i - \alpha, p_i]$. Note that this graph is different from the variation graph, and depends on the variant coordinates and $\alpha$. For brevity, we drop $\alpha$ from $G_\alpha$. Note also that by associating the $i$th row of $\mathbf{A}$ with $v_i$, and interpreting the consecutive nonzeros in the $i$th row as the interval of $v_i$, one obtains the same graph.

**Problems 1'–4':** Given $G = (V, E, w, c)$ and $\delta$ with vertex weights $w$ and vertex charges $c$, find a set of vertices $\mathcal{I}$ with maximum weight $w(\mathcal{I})$, while the total charge of $\mathcal{I}$ around each $v \in V$, defined as $c(\mathrm{adj}^+(v) \cap \mathcal{I})$, is no larger than $\delta$.

For Problem 1, $w(v_i) = c(v_i) = 1$ for all $v_i$; for Problem 2, $w(v_i)$ is the out-degree of $p_i$ in $\mathcal{G}$ minus one, and $c(v_i) = 1$ for a vertex $v_i$; for Problems 3 and 4, $w(v_i)$ and $c(v_i)$ are equal to the corresponding entries in the vectors $\mathbf{w}$ and $\mathbf{c}$ defined above.

We observe that when $\delta = 1$, in Problems 1' and 2', a vertex $v \in \mathcal{I}$ cannot have any of its neighbors in $\mathcal{I}$, and $c(\cdot) = 1$ for all vertices. Therefore, for $\delta = 1$, these two problems are equivalent to the well-known maximum independent set and the maximum weight independent set problems. The case $\delta > 1$ for Problems 1' and 2' therefore is a generalization of the maximum independent set problem. In this generalization, we are asked to find a set of vertices $\mathcal{I}$ such that each vertex in $\mathcal{I}$ can be a neighbor of at most $\delta - 1$ other vertices in $\mathcal{I}$, and all other vertices are allowed to be a neighbor of at most $\delta$ vertices in $\mathcal{I}$. Problems 3' and 4' add another extension in which a vertex $v$ included in $\mathcal{I}$ charges all vertices in its neighborhood by $c(v)$ and a vertex is allowed to have been charged by at most $\delta$. For these reasons, we call a solution $\mathcal{I} \subseteq V$ to Problems 1'–4' a $\delta$-independent set.

The generalized independent set problems are interesting on their own. A related problem is called generalized vertex packing problem [19], which does not have the generalization we need for Problems 3' and 4'. The problems concerning generalized cliques are also relevant [16]—to the best of our knowledge, the generalizations do not cover Problems 3' and 4'. Our algorithms can also be considered as tackling generalized independent set problems on interval graphs (and a vertex is charged by other vertices in its lower-adjacency).

We present algorithms for Problems 1' and 2' and show their properties. We then adapt them for Problems 3' and 4'.
## 3.1 A greedy approach

We first translate an approach proposed by Jain et al. [14] to the graph point of view proposed here. This algorithm, called Greedy in the rest of the paper, visits vertices in the order of their intervals' starting position, and includes a vertex in the $\delta$-independent set $\mathcal{I}$ whenever doing so is feasible. A vertex $v$ can be added in $\mathcal{I}$ if $c(\mathrm{ladj}(v) \cap \mathcal{I}) < \delta$. For Problem 1', this is an exact algorithm. We re-write this as a lemma, whose proof is in [14, Section 3.1.1].

**Lemma 1.** Greedy is exact for Problem 1', in which $w(i) = c(i) = 1$ for all vertices $v_i$.

## 3.2 A sorting-based approach

Sorting-based approaches have been used in finding maximum weighted and unweighted independent sets in general graphs, see for example [13, 18]. We discuss their adaptation to our case.

We propose GreedySort, a sorting-based method built upon Greedy. GreedySort visits the vertices in a nonincreasing order of weights, breaking ties by the order of the vertices (smaller first), and includes a vertex into the $\delta$-independent set, if doing so is feasible. This greedy approach performs badly in general graphs; consider a star of $n+1$ vertices, where the vertex in the middle has a weight $w + \varepsilon$, and all other vertices $w$: GreedySort will pick the middle vertex and none of the others, while the optimal weight is $n \cdot w$. In our case, GreedySort performs much better, as shown in the following lemma, because of the particular form of the graphs that are considered.

**Lemma 2.** GreedySort is $1/2$ approximate for $\delta = 1$ for Problem 2'.

*Proof.* Let $w^*$ be the weight of an optimal solution. Consider a vertex $v$ in an optimal solution, but not chosen by GreedySort. Since $v$ is not chosen by GreedySort, there must be another vertex $u$ in $\mathrm{ladj}(v) \cup \mathrm{hadj}(v)$ with $w(u) \ge w(v)$ chosen by it. If we add these inequalities up over all such $v$, and find the largest number $\ell$ of times the inequalities of the form $w(u) \ge w(\cdot)$ appear for a vertex $u$ chosen by GreedySort, then we can bound GreedySort's approximation as $w^*/\ell$.

Observe that for any vertex $u$, the sets $\mathrm{ladj}^+(u)$ and $\mathrm{hadj}^+(u)$ are both cliques, and hence any feasible solution can have at most two vertices in $\mathrm{ladj}^+(u) \cup \mathrm{hadj}^+(u)$. In other words, a vertex $u$ chosen by GreedySort can appear in at most two inequalities of the form $w(u) \ge w(\cdot)$, concluding the proof. $\square$

In the general case of $\delta \ge 2$, one expects better performance from GreedySort, by following the reasoning in the proof above. This time, for a vertex $v$ not to be chosen by GreedySort, there must be more than $\delta$ vertices with a heavier weight in $\mathrm{ladj}(v)$ or $\mathrm{hadj}(v)$. However this is not so, as shown by the following lemma.

**Lemma 3.** GreedySort is $1/2$ approximate for $\delta \ge 2$ for Problem 2', and this bound is tight.

*Proof.* The proof follows along the same lines as before. Let $\mathcal{I}_O$ be the set of vertices in an optimal solution, and $\mathcal{I}_G$ be the set of vertices picked by GreedySort.

For each $v$ in $\mathcal{I}_O \setminus \mathcal{I}_G$, we can write $\delta$ inequalities of the form $w(u_i) \ge w(v)$, for $i = 1, \dots, \delta$ for vertices $u_i \in \mathrm{ladj}(v) \cap \mathcal{I}_G$; or there is a vertex $z \in \mathrm{hadj}(v)$ where $u_i \in \mathrm{ladj}(z) \cap \mathcal{I}_G$. Each vertex $u_i$ can appear for $\delta$ vertices in $\mathrm{hadj}(u_i)$ and another $\delta$ vertices in $\mathrm{ladj}(u_i)$. Since $\mathrm{ladj}(v) \cup \mathrm{hadj}(v) \cap \mathcal{I}_G$ contains at least $\delta$ vertices for a $v \in \mathcal{I}_O \setminus \mathcal{I}_G$, the vertex $v$ appears in $\delta$ inequalities. This implies that $2\delta w(\mathcal{I}_G) \ge \delta w(\mathcal{I}_O)$, hence the proof.

The bound is tight for graphs of Figure 1, where we have three $\delta$-cliques. All vertices of the first $\delta$-clique are connected to all vertices of the second $\delta$-clique (constituting a 2$\delta$-clique), and all vertices of the second $\delta$-clique are connected to all vertices of the third $\delta$-clique. Assuming a weight of $w + \varepsilon$ for all vertices in the second $\delta$-clique for an infinitesimal $\varepsilon$, and $w$ for all other vertices, GreedySort chooses all $\delta$ vertices in the middle clique, while the first and the third cliques form the optimal solution, completing the proof. $\square$

Figure 1, when $\delta = 1$, shows also a tight example for Lemma 2 which is a three vertex path. Note that the worst case bound of $1/2$ is attained for this specific graphs and hence we expect better performance in practice.

## 3.3 A linear-time algorithm

We propose hMWISd (a heuristic for a $\delta$-independent set) shown in Algorithm 1, a linear-time algorithm for the maximum weighted $\delta$-independent set in an interval graph $G$. While we present this algorithm for interval graphs,
<figure><img src="image_2.png" /><figcaption>Figure 1: An example in which GreedySort's performance is 1/2 for any $\delta \ge 1$. The first and the second cliques, and the second and the third cliques are fully connected.</figcaption></figure>

**Algorithm 1: $\mathcal{I} = hMWISd(G, w, \delta)$**

**Input:** $G$, an interval graph, vertices are ordered
**Input:** $w$, the weight of vertices of $G$
**Output:** $\mathcal{I}$, a maximum weighted $\delta$-independent set in $G$

$\omega \leftarrow w$ \quad /* copy $w$ to $\omega$, which is modified */
$\mathcal{I} \leftarrow \emptyset$
mark($v$) $\leftarrow 0$ for each vertex $v$
**for** $v = 1$ **to** $n$ **do** \quad /* in the given order */
  **if** $\omega(v) > 0$ **then**
    mark($v$) $\leftarrow 1$
    **for** $u \in \text{hadj}(v)$ **do**
      $\lfloor \omega(u) \leftarrow \omega(u) - \omega(v) / \delta \rfloor$
  neighsInI($v$) $\leftarrow 0$ for all $v$
  **for** $v = n$ **downto** 1 **do** \quad /* in the reverse order */
    **if** mark($v$) = 1 **then**
      **if** $\max_{u \in \text{hadj}(v)}(\text{neighsInI}(u)) < \delta$ **then**
        $\mathcal{I} \leftarrow \mathcal{I} \cup \{v\}$
        **for** $u \in \{v\} \cup \text{hadj}(v)$ **do**
          $\lfloor \text{neighsInI}(u) \leftarrow \text{neighsInI}(u) + 1 \rfloor$
it can handle any chordal graph assuming that the vertices are numbered according to a perfect elimination ordering.

The algorithm hMWISd is adapted from an algorithm by Frank [10], which solves the maximum weight independent set problem in chordal graphs. This algorithm performs two passes over the graph. In the first pass, the vertices are visited in the non-decreasing order of their ids, and a set of vertices are marked as candidates. These candidate vertices bear somewhat heavier weight than the vertices in their lower adjacency. In the second pass, the candidate vertices are visited in the reverse order. A candidate vertex $v$ is added to the solution $\mathcal{I}$, if doing so is feasible for the vertices in $\text{hadj}(v)$. For this check, we keep a counter $\text{neighsInI}(\cdot)$ per vertex. Every time we add a vertex $v$ to $\mathcal{I}$, we increment the counter of the vertices in $\text{hadj}(v)$. Adding a vertex $v$ to $\mathcal{I}$ is feasible, if the counters of the vertices in $\text{hadj}(v)$ all are strictly less than $\delta$. Finding the maximum value $\max_{u \in \text{hadj}(v)}(\text{neighsInI}(u))$ can be done in linear time in total. This can be partly avoided by visiting the vertices in $\text{ladj}(u)$ of $u$ to un-mark them, when $\text{neighsInI}(u)$ becomes $\delta$.

We note the following lemma, whose proof is immediate as for the $\delta = 1$ case Algorithm 1 is a version of Frank's algorithm [10].

**Lemma 4.** *Algorithm 1 is exact for $\delta = 1$ for Problem 2'.*

When there are unit weights, the candidates generated by Algorithm 1 at the first for-loop are exactly the output of Greedy or GreedySort. Therefore, Algorithm 1 is exact for Problem 1' which we state as a lemma without proof.

**Lemma 5.** *Algorithm 1 is exact for Problem 1'.*

### 3.4 An algorithm with a stable performance

We propose hMWISdDec (a heuristic for $\delta$-independent set by decomposing $\delta$) whose run time is $O(\delta(m + n))$, and hence suitable for small $\delta$. The proposed algorithm is shown in Algorithm 2. It calls $\delta$ times Algorithm 1 to find a maximum weighted independent set, where the weights change between two calls. When a set of vertices is found at a step $d$ for $d = 1, \dots, \delta$ they are added to the output, and their weight is set to zero. This algorithm makes use of the fact that the exact algorithm hMWISd ($G, \omega, 1$) does not include any vertex $v$ with a zero weight in the solution.

**Algorithm 2:** $\mathcal{I} = \text{hMWISdDec}(G, w, \delta)$

**Input:** $G$, an interval graph, vertices are ordered

**Input:** $w$, the weight of vertices of $G$

**Input:** $\delta$, the constraint

**Output:** $\mathcal{I}$, a maximum weighted $\delta$-independent set in $G$

$\omega \leftarrow w$

/* the algorithm modifies the weights $\omega$ */

$\mathcal{I} \leftarrow \emptyset$

for $d = 1$ to $\delta$ do

$\mathcal{I}_d \leftarrow \text{hMWISd}(G, \omega, 1)$

$\mathcal{I} \leftarrow \mathcal{I} \cup \mathcal{I}_d$

$\omega(\mathcal{I}_d) \leftarrow 0$

/* for each $v \in \mathcal{I}_d$, set $\omega(v) = 0$ */

The algorithm does not keep track of the number of neighbors of a vertex included in the output, yet it obtains a $\delta$-independent set as we show below.

**Lemma 6.** *Algorithm 2 obtains a $\delta$-independent set.*

*Proof.* We have to show that for a vertex $v$, the set $\text{ladj}^+(v)$ contains at most $\delta$ vertices in $\mathcal{I}$. To do so, we show that for a vertex $v$ and any $d = 1, \dots, \delta$, we have $|\text{ladj}^+(v)| \cap \mathcal{I}_d| \le 1$.

If a vertex $v$ is selected to be in the independent set $\mathcal{I}_d$, then none of its lower numbered neighbors can be in $\mathcal{I}_d$.

If $v$ is not in $\mathcal{I}_d$, for any $u_1, u_2 \in \text{ladj}(v)$ only one of them can be in $\mathcal{I}_d$. This is so, as $\text{ladj}(v)$ is a clique in an interval graph, and hence two vertices in $\text{ladj}(v)$ cannot be in an independent set. $\square$
<figure><img src="image_3.png" /><figcaption>(a) A graph, with a maximum weight independent set of weight 9. The sets {1, 5} and {3, 5} are both optimal.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) After finding {1, 5} as a maximum weight independent set, {3, 6} becomes the maximum weight independent set.</figcaption></figure>

Figure 2: Vertices are labeled 1–6, and their weights are given either above or below them.

The algorithm above is not exact for $\delta \ge 2$. An example is shown in Figure 2. For this graph and $\delta = 2$, the optimal solution is formed by vertices {1, 3, 4, 5} and has a weight of 17, as can be verified in Figure 2a. If Algorithm 2 finds {1, 5} as a maximum weighted independent set in the first step, then it finds {3, 6} as the maximum weighted independent set in the second step as shown in Figure 2b.

## 3.5 A dynamic programming formulation

We propose a dynamic programming method to exactly solve Problem 2'. Our algorithm has a run time of $O(n^\delta)$. Despite the large run time, this exact solver could be useful in a toolbox, as it is simpler to use than an LP solver. It can also be used to solve the variant reduction problem for a shorter sub-chain of the input.

$T(i, n_1, n_2, \dots, n_\delta)$ is the maximum possible weight $\sum_{i \in I} w_i$ that can be obtained by selecting a subset of vertices $I \subseteq [v_i, v_n]$, knowing that the choice of $I$ might be constrained by other vertices from $[v_1, v_{i-1}]$ that will be part of the final set of selected vertices. The goal is to compute $T(1, 0, \dots, 0)$, which is the total weight of the selected vertices, with no constraints. The numbers $n_1, n_2, \dots, n_\delta$ tell us until which index the previously placed vertices have an impact in a solution. Hence, if $n_1 \le i$, it means that there are no more than $\delta$ chosen vertices in $[v_{i-\delta+1}, v_i]$, and it is therefore possible to add $v_i$ to the set of selected variants $I$. Otherwise, the first vertices that might be selected in $I$, given the current decisions, is $v_{n_1}$.

The recurrence is as follows:

$$T(i, n_1, n_2, \dots, n_\delta) = \max \begin{cases} T(i+1, n_1, n_2, \dots, n_\delta) \\ w(i) + T(i+1, n_2, n_3, \dots, n_{\delta-1}, i+\alpha) \\ \text{if } n_1 \le i \end{cases}$$

The first line corresponds to not choosing $v_i$, hence the next spots of available positions where we can choose new vertices remain the same. The second line corresponds to choosing $v_i$, which is possible only if $n_1 \le i$. In that case, we shift all the $n_i$'s (corresponding to the previous $\delta-1$ vertices that were chosen before $v_i$), and add the $n_\delta$ when $v_i$ will not be within distance $\alpha$ of it, hence not conflicted. The goal is to compute $T(1, 0, \dots, 0)$, with the initialization conditions $T(n+1, n_1, n_2, \dots, n_\delta) = 0$ for any values $n_1, n_2, \dots, n_\delta$.

We start by computing all values $T(n, n_1, \dots, n_\delta)$ for all values of $n_1, \dots, n_\delta$ by setting $T(n, n_1, \dots, n_\delta) = w_n$ if $n_1 \le n$ (we can select a vertex $n$ so we add it to our set, given previous choices), and 0 otherwise (there are already $\delta$ vertices in an interval of length $\alpha$ including vertex $n$). Subsequently, all values of $T(n-1, n_1, \dots, n_\delta)$ (case $i = n-1$) are computed, which reuse the values with $i = n$ that were previously computed. This iterative process is continued until $i = 1$ is reached, obtaining the desired result.

The complexity of this algorithm is $O(n^{\delta+1})$, as we compute all values of $T(i, n_1, \dots, n_\delta)$ with each index ranging from 1 to $n$.

## 3.6 Addressing Problems 3 and 4

For these two problems, GreedySort visits the vertices in a non-increasing order of $w(i)/c(i)$, breaking ties by id (orders the smaller id first), and includes a vertex into the solution if feasibility is met.

The adaptation of hMWISd keeps the first pass the same. In the second pass (which is in reverse order), vertices are included in the solution if feasibility is met. This is done by first ensuring that

$$\max_{u \in \text{hadj}(v)} (\text{neighsInI}(u)) \le \delta - c(v),$$

and then by incrementing the weights of vertices in hadj$^+$($v$) by $c(v)$ when $v$ is included in the solution. Since the charges of vertices vary, the solution at the end is not necessarily maximal, and we can possibly include
<table><caption>Table 1: Properties of graphs obtained for chromosomes (Chr) 1–22 with varying $\alpha$. The maximum weight of a vertex ($w_{max}$) and the number of vertices ($|V|$) do not vary with $\alpha$. The average number of vertices in an $\text{ladj}^+(\cdot)$ is given for varying $\alpha$.</caption><thead><tr><th rowspan="2">Chr</th><th rowspan="2">w<sub>max</sub></th><th rowspan="2">|V|</th><th colspan="4">Average |ladj<sup>+</sup>(·)| for α =</th></tr><tr><th>150</th><th>1,000</th><th>5,000</th><th>10,000</th></tr></thead><tbody><tr><td>1</td><td>3</td><td>6,215,031</td><td>5.5</td><td>30.0</td><td>144.2</td><td>286.2</td></tr><tr><td>2</td><td>3</td><td>6,807,707</td><td>5.6</td><td>31.0</td><td>148.8</td><td>295.2</td></tr><tr><td>3</td><td>3</td><td>5,602,388</td><td>5.6</td><td>30.8</td><td>148.6</td><td>295.1</td></tr><tr><td>4</td><td>3</td><td>5,499,200</td><td>5.7</td><td>31.6</td><td>152.2</td><td>302.0</td></tr><tr><td>5</td><td>3</td><td>5,054,708</td><td>5.5</td><td>30.8</td><td>148.4</td><td>294.7</td></tr><tr><td>6</td><td>3</td><td>4,816,035</td><td>5.6</td><td>31.1</td><td>149.5</td><td>296.5</td></tr><tr><td>7</td><td>3</td><td>4,532,453</td><td>5.8</td><td>32.1</td><td>154.5</td><td>306.5</td></tr><tr><td>8</td><td>3</td><td>4,433,731</td><td>6.2</td><td>35.1</td><td>169.6</td><td>337.1</td></tr><tr><td>9</td><td>3</td><td>3,426,740</td><td>5.9</td><td>33.0</td><td>158.9</td><td>315.7</td></tr><tr><td>10</td><td>3</td><td>3,836,552</td><td>5.7</td><td>32.0</td><td>154.2</td><td>306.2</td></tr><tr><td>11</td><td>3</td><td>3,890,938</td><td>5.7</td><td>32.1</td><td>154.8</td><td>307.4</td></tr><tr><td>12</td><td>3</td><td>3,709,628</td><td>5.5</td><td>30.6</td><td>147.2</td><td>292.3</td></tr><tr><td>13</td><td>3</td><td>2,736,545</td><td>5.5</td><td>30.5</td><td>146.8</td><td>291.5</td></tr><tr><td>14</td><td>3</td><td>2,547,626</td><td>5.7</td><td>31.3</td><td>149.7</td><td>296.6</td></tr><tr><td>15</td><td>3</td><td>2,328,194</td><td>5.7</td><td>31.7</td><td>152.1</td><td>301.9</td></tr><tr><td>16</td><td>3</td><td>2,606,617</td><td>6.7</td><td>38.4</td><td>185.8</td><td>369.2</td></tr><tr><td>17</td><td>3</td><td>2,234,274</td><td>5.7</td><td>31.3</td><td>150.5</td><td>298.7</td></tr><tr><td>18</td><td>3</td><td>2,178,376</td><td>5.6</td><td>31.0</td><td>149.5</td><td>297.1</td></tr><tr><td>19</td><td>3</td><td>1,758,134</td><td>6.2</td><td>34.5</td><td>165.9</td><td>328.9</td></tr><tr><td>20</td><td>3</td><td>1,744,899</td><td>5.6</td><td>31.3</td><td>151.0</td><td>300.1</td></tr><tr><td>21</td><td>3</td><td>1,058,338</td><td>5.9</td><td>32.8</td><td>158.1</td><td>313.8</td></tr><tr><td>22</td><td>3</td><td>1,059,516</td><td>6.2</td><td>34.3</td><td>163.4</td><td>323.7</td></tr></tbody></table>

more vertices. We do so by visiting vertices in the increasing order of ids and include additional vertices in the solution, as long as it is feasible to do so.

# 4 Experimental results

We use the framework by Jain et al. [15], which uses vcftools [5] to parse SNPs and indels from the 1000 Genomes Project data [22]. Our codes are available at https://gitlab.inria.fr/bora-ucar/vf-graphs. We tested the algorithms with $\alpha$ values of 150, 1,000, 5,000, and 10,000. For each $\alpha$, we experimented with three values of $\delta$, which are 0.01, 0.05, and 0.1 of $\alpha$. These choices are found useful for short and long-read sequencing and low-to-high error tolerance [15].

## 4.1 SNP-only problems

We first present results for Problems 1 and 2, which only allow SNP variations. Table 1 presents the properties of graphs $G_\alpha$ obtained from variation graphs $G$ of the human chromosomes 1–22. In this table, the column “Chr” gives the chromosome number; the column “$w_{max}$” gives the maximum weight of a vertex, where the minimum was one for all graphs; the column “$|V|$” gives the number of vertices in $G_\alpha$; the last three columns give the average number of vertices in an $\text{ladj}^+(\cdot)$. We compared the proposed algorithms, the greedy algorithm from earlier work [14], and used Gurobi version 11.0.3 to solve the problems optimally using the ILP formulation.

Recall that all three algorithms GreedySort, hMWISd, hMWISdDec are exact for the case of unit weights and unit charges, thus solving Problem 1 optimally. Therefore, the purpose of the comparison is to assess the run time taken by the different algorithms. GreedySort performs a single pass over data assuming it does not need to sort, and is thus the most efficient one. hMWISd performs two passes, while hMWISdDec performs $2\delta$ passes. While hMWISd has the same worst-case run time complexity for the unit and non-unit weighted cases, GreedySort needs to know this beforehand, or check the weights beforehand to avoid sorting. We present run times in seconds with the unit weights, for graphs of chromosomes 1, 11, and 22, in Table 2. The associated
<table><caption>Table 2: Run time in seconds of various algorithms for $\alpha$ = 5,000 and $\delta$ = 100, on the variation graphs for chromosomes 1, 11, and 22.</caption><thead><tr><th>Chr</th><th>Gurobi</th><th>Greedy</th><th>hMWISd</th><th>hMWISdDec</th><th>GreedySort</th></tr></thead><tbody><tr><td>1</td><td>3,944.554</td><td>0.100</td><td>1.985</td><td>11.110</td><td>2.386</td></tr><tr><td>11</td><td>2,948.298</td><td>0.063</td><td>1.408</td><td>5.646</td><td>1.500</td></tr><tr><td>22</td><td>963.736</td><td>0.018</td><td>0.397</td><td>0.982</td><td>0.414</td></tr></tbody></table>

graphs are among those with the largest, medium, and smallest number of vertices. We used $\alpha$ = 5,000 and $\delta$ = 100 in these experiments.

The run times reported in Table 2 are computed on a system with Intel Xeon CPU E7-8890 v4 2.20GHz. When Gurobi is called to solve an ILP, it first solves a root relaxation (which is an LP). Since the LP formulation has integral solutions, calling Gurobi for solving the ILP has the same effect. In our experience, the ILP solver is much more memory efficient than the LP solver. We thus call Gurobi's ILP solver, which returns an optimal solution without having to perform branch-and-bound and related operations. We fed the solution found by GreedySort as an initial solution to Gurobi and allowed 4 hours of compute time. As seen from this table, Greedy is the fastest, followed by GreedySort, and then by hMWISdDec. GreedySort uses stable_sort from C++std, which could have some overheads. The run time of GreedySort is low enough that we did not try to improve its run time.

We next compare the algorithms GreedySort, Greedy, and hMWISdDec on Problem 2, where the weights are non-uniform. For this problem, the algorithms are not exact. Therefore, we investigate the quality of the solutions found by the proposed algorithms. A point in Figure 3 shows the ratio of the number of variants retained by an algorithm to the optimal solution (on the y-axis) for different $\alpha$ and $\delta$ configurations on all 22 chromosome graphs (on the x-axis). In these plots, the smaller ratios are better, as we want to minimize the number of variations to be retained while ensuring $(\alpha, \delta)_h$-compatibility. In eight cases with $\alpha$=10,000 and $\delta$ = 100, the ILP solver timed out due to the associated graphs having a large number of edges, and the pre-solve and pre-processing of the Gurobi solver took a long time. In those cases, we have computed the ratios using the best results at hand. In 10 cases where $\alpha$=10,000 and $\delta$ = 1,000, and in another five cases where $\alpha$=5,000 and $\delta$ = 500, the error tolerance was so high that an optimal solution discarded all the variations; we have chosen not to include them in the results.

As seen in these plots, all algorithms, even Greedy which does not consider the weights, are very close to the optimal. This is presumably due to the maximum and minimum weights being 3 and 1 (see Table 1), which is nearly uniform, resulting in problem instances close to those for which Greedy is exact. Even though the difference is small, GreedySort is the best-performing method—its maximum deviation from the optimal was for the graph of chromosome 14, $\alpha$=1,000, and $\delta$ = 50, in which case the optimal value is 24,653 and GreedySort found a solution with retaining 24,665 variations, just two more than the optimal.

In order to assess the performance of the four algorithms in more challenging cases, we have created a family of interval graphs, parametrized with $L$, the domain of the intervals of vertices; $n$, the number of vertices; and $\alpha$, the length of the intervals. We choose $n$ integer points, no repetitions, from $[1, L]$ as vertices and connect two vertices with an edge if the corresponding points are less than $\alpha$ away from each other. The expected degree of a vertex is $\Delta = \frac{\alpha}{L/n}$. This setting corresponds to a variation graph with $n$ variation sites. We then assign an integral weight, sampled uniformly at random between 2 and $\alpha$ to each vertex. In our experiments we use $L$=500,000, $n$=10,000 and $\alpha$= 5,000, in which case the expected degree $\Delta$ of a vertex is 100. We then choose a random $\delta$ in range $[2, \Delta]$; as $\delta$ is less than the average degree of a vertex, the decisions to be made at each vertex are non-trivial. A graph and a chosen $\delta$ define a problem instance. We created 2,048 problem instances with the above parameters and compared the four algorithms with an optimal solution in terms of the weight of the set to be removed (hence the larger the better). The optimal solution was found by Gurobi in all 2,048 instances. In Table 3, we report the minimum, the maximum, and the geometric mean of the ratio of each algorithm. As seen from this table, on average, GreedySort and hMWISdDec are the best methods, with solution qualities of 0.97 and 0.98 on average, where hMWISdDec's performance is much more stable. We also note that hMWISd found an optimal solution in 19 instances, while hMWISdDec solved 24 instances optimally—the other two methods did not find any optimal solutions. We thus confirm that all three algorithms proposed in this work deliver good results in general with at least 0.93 of the optimal solution, while hMWISdDec's performance is more stable.
<figure><img src="image_5.png" /><figcaption>Figure 3: Performance of different algorithms for Problem 2, for varying α and δ. x-axis: 22 chromosome graphs; y-axis: the ratio of the number of variants retained by an algorithm to that retained in an optimal solution—lower values indicate better performance. Nothing is plotted for problems in which the optimal solution retained only the reference haplotype.</figcaption></figure>

<table><caption>Table 3: The weight of the δ-independent sets found by different algorithms normalized to the optimal value.</caption><thead><tr><th></th><th>Greedy</th><th>hMWISd</th><th>hMWISdDec</th><th>GreedySort</th></tr></thead><tbody><tr><th>min</th><td>0.48</td><td>0.88</td><td>0.97</td><td>0.78</td></tr><tr><th>max</th><td>0.94</td><td>1.00</td><td>1.00</td><td>1.00</td></tr><tr><th>geomean</th><td>0.70</td><td>0.93</td><td>0.98</td><td>0.97</td></tr></tbody></table>
<table><caption>Table 4: Properties of graphs obtained for chromosomes (Chr) 1–22 with varying $\alpha$. The maximum weight of a vertex ($w_{max}$), the maximum charge of a vertex ($c_{max}$), and the number of vertices ($|V|$) do not vary with $\alpha$. The average number of vertices in an $\text{ladj}^+$ is given for different $\alpha$.</caption><thead><tr><th rowspan="2">Chr</th><th rowspan="2">w<sub>max</sub></th><th rowspan="2">c<sub>max</sub></th><th rowspan="2">|V|</th><th colspan="4">Average |ladj<sup>+</sup>(.)| for α =</th></tr><tr><th>150</th><th>1,000</th><th>5,000</th><th>10,000</th></tr></thead><tbody><tr><td>1</td><td>5</td><td>389</td><td>6,453,040</td><td>5.6</td><td>31.1</td><td>149.9</td><td>297.5</td></tr><tr><td>2</td><td>6</td><td>166</td><td>7,065,050</td><td>5.7</td><td>32.1</td><td>154.6</td><td>306.7</td></tr><tr><td>3</td><td>6</td><td>341</td><td>5,818,558</td><td>5.7</td><td>32.0</td><td>154.4</td><td>306.9</td></tr><tr><td>4</td><td>6</td><td>167</td><td>5,718,840</td><td>5.8</td><td>32.8</td><td>158.4</td><td>314.5</td></tr><tr><td>5</td><td>5</td><td>136</td><td>5,253,166</td><td>5.7</td><td>32.0</td><td>154.3</td><td>306.7</td></tr><tr><td>6</td><td>5</td><td>279</td><td>5,011,704</td><td>5.8</td><td>32.4</td><td>155.7</td><td>309.1</td></tr><tr><td>7</td><td>6</td><td>283</td><td>4,705,084</td><td>5.9</td><td>33.3</td><td>160.5</td><td>318.6</td></tr><tr><td>8</td><td>6</td><td>159</td><td>4,586,542</td><td>6.3</td><td>36.2</td><td>175.3</td><td>348.4</td></tr><tr><td>9</td><td>7</td><td>411</td><td>3,552,230</td><td>6.0</td><td>34.1</td><td>164.8</td><td>327.5</td></tr><tr><td>10</td><td>6</td><td>416</td><td>3,982,654</td><td>5.9</td><td>33.2</td><td>160.3</td><td>318.3</td></tr><tr><td>11</td><td>8</td><td>136</td><td>4,036,074</td><td>5.9</td><td>33.2</td><td>160.7</td><td>319.1</td></tr><tr><td>12</td><td>6</td><td>661</td><td>3,858,551</td><td>5.7</td><td>31.8</td><td>153.2</td><td>304.4</td></tr><tr><td>13</td><td>5</td><td>143</td><td>2,850,983</td><td>5.7</td><td>31.7</td><td>153.1</td><td>304.2</td></tr><tr><td>14</td><td>5</td><td>124</td><td>2,648,677</td><td>5.8</td><td>32.5</td><td>155.7</td><td>308.7</td></tr><tr><td>15</td><td>6</td><td>171</td><td>2,418,975</td><td>5.9</td><td>32.9</td><td>158.2</td><td>314.2</td></tr><tr><td>16</td><td>8</td><td>217</td><td>2,691,570</td><td>6.8</td><td>39.5</td><td>191.9</td><td>381.3</td></tr><tr><td>17</td><td>6</td><td>327</td><td>2,323,202</td><td>5.8</td><td>32.5</td><td>156.7</td><td>311.2</td></tr><tr><td>18</td><td>5</td><td>202</td><td>2,261,738</td><td>5.7</td><td>32.2</td><td>155.5</td><td>309.0</td></tr><tr><td>19</td><td>5</td><td>646</td><td>1,827,391</td><td>6.3</td><td>35.8</td><td>172.5</td><td>342.3</td></tr><tr><td>20</td><td>5</td><td>403</td><td>1,808,612</td><td>5.8</td><td>32.4</td><td>156.8</td><td>311.6</td></tr><tr><td>21</td><td>5</td><td>101</td><td>1,102,670</td><td>6.1</td><td>34.2</td><td>165.0</td><td>327.6</td></tr><tr><td>22</td><td>5</td><td>413</td><td>1,100,716</td><td>6.3</td><td>35.6</td><td>170.1</td><td>337.0</td></tr></tbody></table>

## 4.2 SNP and indels

We now present results for Problems 3 and 4, in which we deal with SNP and indel variations. Table 4 presents the properties of graphs $G_\alpha$ obtained from variation graphs $\mathcal{G}$ of the chromosomes 1–22. In this table, the column “Chr” gives the chromosome number; the column “$w_{max}$” gives the maximum weight of a vertex, where the minimum was one for all graphs; the column “$c_{max}$” gives the maximum charge of a vertex, where the minimum was one for all graphs; the column “$|V|$” gives the number of vertices in $G_\alpha$; the last three columns give the average number of vertices in an $\text{ladj}^+(.)$.

We compare the algorithms *GreedySort* and *hMWISd* on Problem 4, where the weights are non-uniform, and the charges are defined according to the variation graphs. A point in Figure 4 shows the ratio of the number of variants retained by an algorithm to the optimal solution (on the y-axis) for different $\alpha$ and $\delta$ configurations on all 22 chromosome graphs (on the x-axis). In these plots, the smaller ratios are better. As the performance of the algorithms varied significantly for different $\alpha$ and $\delta$ settings, the plots use different scales for the y-axis. In 36 cases, the ILP solver timed out, and we used the best solution at hand for normalization. In 5 cases where $\alpha=10,000$ and $\delta=1,000$, and in another case where $\alpha=5,000$ and $\delta=500$, an optimal solution without any of the variants was possible; we have discarded them. As seen in these plots, the best performing algorithm is *GreedySort* in all cases. *hMWISd* is usually better than *Greedy*; for the graph of chromosome 9, $\alpha=10,000$, and $\delta=1,000$, the performance of *Greedy* is significantly better than that of *hMWISd*. The maximum deviation of *GreedySort* from the optimum is 1.08, on the graph of chromosome 8, with $\alpha=10,000$ and $\delta=1,000$, where the optimal solution retains 202 variants and *GreedySort* found a solution that retains 219. Among other interesting data points in the figure where the numbers of retained variants are more sizable, we note the result on the graph of chromosome 9, with $\alpha=5,000$, $\delta=250$, where the optimum is 15,005, *GreedySort* is the closest at 15,293, and *Greedy* and *hMWISd* found solutions retaining 23,802 and 24,054 variants, respectively.

We also assessed the performance of the algorithms on Problem 3, where all vertex weights are unit (and other parameters are the same as those of Problem 4). Since weights were nearly uniform in Problem 4, one expects the performance of all algorithms on Problem 3 to be similar to their respective performance on Problem 4.
<figure><img src="image_6.png" /><figcaption>Figure 4: Performance of different algorithms for Problem 4, for varying α and δ. x-axis: 22 chromosome graphs; y-axis: the ratio of the number of variants retained by an algorithm to that retained in an optimal solution. In the first two column of plots, Greedy and hMWISd have similar performance, while GreedySort performs better than both globally.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Performance of different algorithms with respect to the optimal solution for Problem 4, for varying $\alpha$ and $\delta$, with random vertex weights between 2 and $\alpha$. x-axis: 22 chromosome graphs; y-axis: the ratio of the number of variants retained by an algorithm to that retained in an optimal solution—lower values indicate better performance.</figcaption></figure>

4. We have confirmed this by running the algorithms on a subset of problem instances, with $\alpha$ values of 150 and 1,000, and $\delta = \{0.01, 0.05\}$ of each $\alpha$. The geometric mean of the ratio of the values of the optimal solutions of Problems 4 and 3 is 1.02. The same ratio is observed for each algorithm (the ratio of the value of the solution obtained by an algorithm for Problems 4 and 3). That is, the plots in Figure 4 are also representative for Problem 3. As GreedySort is always better than the other two alternatives, and its results are always very close to the optimal solution, we conclude that it is the method of choice for Problems 3 and 4.

In order to assess the performance of the three algorithms on more challenging cases, we have created instances of Problem 4' using the same graphs of chromosomes 1-22, except with vertex weights randomly chosen between 2 and $\alpha$. We did not change the charges of the vertices as that parameter affects feasibility with respect to a given $\delta$. We have used $\alpha$ values of 150, 1,000 and 5,000 and $\delta$ to be 0.01, 0.05 and 0.1 of a given $\alpha$ as before. The optimal solutions to all problem instances were found by Gurobi. The results are plotted in Figure 5, with the same signification of data points and the axes as in the previous plots. As seen in this figure, GreedySort is again the best solver, with performance always close to the optimal solution. The largest ratio of any solution found by GreedySort to the optimal is 1.88 and observed on the graph of chromosome 17, with $\alpha = 5,000$ and $\delta = 500$, where the optimal solution has the value 4,022 and GreedySort obtained 7,569; Greedy and hMWISd obtained solutions with values, respectively, 79,094 and 6,393. This was one of the extreme cases, as the initial variation graph is largely reduced (the number of variant containing positions was reduced from 2,323,202 to 25). The next largest ratio of any solution found by GreedySort was 1.15 on the graph of chromosome 9, with $\alpha = 1,000$, $\delta = 100$; the optimal, GreedySort', Greedy', and hMWISd's values were 129,031, 149,400, 612,339, and 325,244, respectively. In all other instances, the largest ratio of a solution found by GreedySort to the optimal value was 1.09. We thus conclude that GreedySort is the method of choice for Problem 4' for a wide range of parameters.

Recall that GreedySort's comparison function accesses two different arrays, and hMWISd visits vertices to ensure maximality—in most cases one pass was sufficient. These increase run time slightly, while allowing both heuristics to complete in a fraction of the run time of the Gurobi solver. For example, on the graph of chromosome 1, with $\alpha = 5,000$ and $\delta = 250$, GreedySort runs in 11.5 seconds (its result is 1.02 of the optimal) hMWISd runs in 23.5 seconds (its result is 1.12 of the optimal), Greedy completes in 0.45 seconds (is result is 1.42 of the optimal), and Gurobi runs in 8,498.1 seconds.
# 5 Conclusion and future work

We have revisited four variant selection problems in a mathematical framework that allows provable guarantees concerning the retained variants. Prior work in this area proposed a practical heuristic without any performance guarantee in general, and integer linear programming formulations to solve the problems exactly. We have improved upon this work by proposing three heuristics and a dynamic programming formulation for two of these problems which are solvable in polynomial time, and adapted two of these heuristics for the remaining two NP-hard problems. Our heuristics are designed to be fast, robust, and deliver optimal or near-optimal results for all problems, some of which with provable guarantees. In particular, the most versatile and effective heuristic (GreedySort) is optimal for the first problem, 1/2 approximate for the second one, and has yielded results that are within 1.08 of the optimal for the NP-hard problems. This method is also very fast: it takes only a few seconds on one of the largest problem instances, while obtaining a solution within a factor 1.02 of the optimal.

To the best of our knowledge, the graph theoretical formulation of Problems 3 and 4 gives rise to new independent set-like problems. There are multiple opportunities to explore this connection. One can devise data reduction methods for the ILP formulations [12]. Second, one can address these problems in a dynamic context, where the solution set is incrementally updated by the addition of new haplotypes. We have not addressed the haplotype-aware versions of the problems [21]. One way to address them is to consider a set of haplotypes without common variations (with respect to a reference haplotype) and solve the problem using the algorithms developed in our work. When haplotypes share variations, conflicting decisions arise; a good way to resolve them should be investigated. Incorporating complex genomic variations, such as inversions and tandem duplications, in the framework would certainly be of interest.

# References

[1] Ravindra K. Ahuja, Thomas L. Magnanti, and James B. Orlin. *Network flows: theory, algorithms, and applications*. Prentice-Hall, Inc., USA, 1993.

[2] Andreas Brandstädt, Van Bang Le, and Jeremy P. Spinrad. *Graph Classes: A Survey*. Society for Industrial and Applied Mathematics, 1999.

[3] Li Chen, Rasmus Kyng, Yang P. Liu, Richard Peng, Maximilian Probst Gutenberg, and Sushant Sachdeva. Almost-linear-time algorithms for maximum flow and minimum-cost flow. *Commun. ACM*, 66(12):85–92, November 2023.

[4] Michael B. Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time. In *Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing*, STOC 2019, page 938–942, New York, NY, USA, 2019. Association for Computing Machinery.

[5] Petr Danecek, Adam Auton, Goncalo Abecasis, Cornelis A. Albers, Eric Banks, Mark A. DePristo, Robert E. Handsaker, Gerton Lunter, Gabor T. Marth, Stephen T. Sherry, Gilean McVean, Richard Durbin, and 1000 Genomes Project Analysis Group. The variant call format and vcftools. *Bioinformatics*, 27(15):2156–2158, 06 2011.

[6] Agnieszka Danek, Sebastian Deorowicz, and Szymon Grabowski. Indexes of large genome collections on a PC. *PloS one*, 9(10):e109384, 2014.

[7] Alexander Dilthey, Charles Cox, Zamin Iqbal, Matthew R Nelson, and Gil McVean. Improved genome inference in the MHC using a population reference graph. *Nature genetics*, 47(6):682–688, 2015.

[8] Michael Dom. Algorithmic aspects of the consecutive-ones property. *Bull. EATCS*, 98:27–59, 2009.

[9] Hannes P. Eggertsson, Hakon Jonsson, Snaedis Kristmundsdottir, Eirikur Hjartarson, Birte Kehr, Gisli Masson, Florian Zink, Kristjan E. Hjorleifsson, Aslaug Jonasdottir, Adalbjorg Jonasdottir, et al. Graphtyper enables population-scale genotyping using pangenome graphs. *Nature genetics*, 49(11):1654, 2017.

[10] András Frank. Some polynomial algorithms for certain graphs and hypergraphs. In Nash-Williams, C. St. J. A. and Sheen, J., editors, *Proceedings of the 5th British Combinatorial Conference*, pages 211–226. Utilitas Mathematica, 1975.
[11] Martin Charles Golumbic. *Algorithmic Graph Theory and Perfect Graphs*, volume 57 of *Annals of Discrete Mathematics*. North-Holland Publishing Co., Amsterdam, The Netherlands, The Netherlands, 2004.

[12] Ernestine Großmann, Kenneth Langedal, and Christian Schulz. A comprehensive survey of data reduction rules for the maximum weighted independent set problem, 2024.

[13] Magnús Halldórsson and Jaikumar Radhakrishnan. Greed is good: approximating independent sets in sparse and bounded-degree graphs. In *Proceedings of the Twenty-Sixth Annual ACM Symposium on Theory of Computing*, STOC '94, page 439–448, New York, NY, USA, 1994. Association for Computing Machinery.

[14] Chirag Jain, Neda Tavakoli, and Srinivas Aluru. A variant selection framework for genome graphs. *Bioinformatics*, 37(Supplement 1):i460–i467, 07 2021.

[15] Chirag Jain, Neda Tavakoli, and Srinivas Aluru. Vf: Variant selection framework in genome graphs (v1.0), 2021.

[16] Jeffrey Pattillo, Nataly Youssef, and Sergiy Butenko. On clique relaxation models in network analysis. *European Journal of Operational Research*, 226(1):9–18, 2013.

[17] Jacob Pritt, Nae-Chyun Chen, and Ben Langmead. Forge: prioritizing variants for graph genomes. *Genome Biology*, 19:1–16, 2018.

[18] Shuichi Sakai, Mitsunori Togasaki, and Koichi Yamazaki. A note on greedy algorithms for the maximum weighted independent set problem. *Discrete Applied Mathematics*, 126(2):313–322, 2003.

[19] Hanif D. Sherali and J. Cole Smith. A polyhedral study of the generalized vertex packing problem. *Math. Program.*, 107(3):367–390, July 2006.

[20] Robert E. Tarjan and Mihalis Yannakakis. Simple linear-time algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs. *SIAM Journal on Computing*, 13(3):566–579, 1984.

[21] Neda Tavakoli, Daniel Gibney, and Srinivas Aluru. Graphslimmer: Preserving read mappability with the minimum number of variants. *Journal of Computational Biology*, 31(7):616–637, 2024.

[22] The 1000 Genomes Project Consortium. A global reference for human genetic variation. *Nature*, 526(7571):68–74, 2015.