DiVerG: Scalable Distance Index for Validation of Paired-End Alignments in Sequence Graphs

Sequence graphs provide compressed representations of sets of similar – evolutionary or environmentally related – genomic sequences [25]. They are typically defined as node-labelled bidirected graphs where each node contains a string label. Bidirectionality allows traversal in both forward and reverse directions, which accounts for the complementary structure of genomes. For theoretical analysis, we use an equivalent representation called a character graph, defined as $G(V,E,\lambda )$ where $V$ and $E$ are the node and edge sets, respectively, and $\lambda :V\to \mathrm{\Sigma }$ assigns a single character from alphabet $\mathrm{\Sigma }$ to each node. Character graphs and general sequence graphs can be converted to each other in time linear in the size of the character graph [27]. By definition, in a character graph, the length of a sequence spelt out by a walk of length $k$ is exactly $k$.

Lastly, a chain graph is a character graph that represents a single sequence which is equivalent to the sequential chaining of all (single-letter) nodes. A sequence graph $G$ is called sparse if the average node degree in $G$ is close to 1.

Fragments in sequencing libraries typically vary in size depending both on the sequencing technology in use and library preparation procedures. In our study, the distance between paired reads in a library is modelled by an interval indicating the expected lower ($d_1$) and upper bounds ($d_2$). This interval, which is referred to as distance constraints and denoted as $(d_1,d_2)$, are assumed to be provided as input parameters. In Appendix A, we discussed how the distance constraints can be determined for a library.

As stated before, aligning short reads to a reference genome, whether linear or graph, can be ambiguous due to several factors such as repetitive regions in the genome, sequencing errors, and the genetic differences between the donor genome and the reference. In such cases, reads might have multiple candidate alignments. Utilizing the distance information in paired-end sequencing data can help identify the correct alignments (Figure 6(b) in Appendix A).

Two alignments are considered as paired if the the distance between them is consistent with the inferred fragment model. In addition to the distance, the orientation of two paired reads should also be as expected depending on the sequencing technology in use; e.g. one of the pair should be aligned on the forward strand while the other is on the reverse strand.

Due to the sheer number of queries in a typical read mapping experiment, the algorithm solving Section 2.2 needs to be efficient.

The distance index $𝒯$ computed by Equation 1 quickly becomes infeasible for large graphs due to its space requirements. For example, consider a human pangenome graph constructed using the autosomes of the GRCh37 reference genome and incorporating variants from the 1000 Genome Project. Constructing the distance index $𝒯$ for this graph with distance criteria $d_1=150$ and $d_2=450$ results in a matrix of order 2.8B with approximately 870B non-zero values. Assuming that each non-zero value consumes 4 bytes, the total space required for the final distance matrix in sorted CRS would be about 3.5 TB.

Even for smaller graphs, storage requirements are impractical for GPU architectures as they have much smaller memory than the accessible main memory in a CPU. Such many-core architectures are beneficial for faster computation of SpGEMM, which is the computational bottleneck for constructing the distance index.

A key observation for identifying the compression potential is the structure of the distance matrix. Genome graphs usually consist of multiple connected components corresponding to each genomic region or chromosome. If the nodes are indexed such that their indices are localized for each region, the distance index $𝒯$ is a block-diagonal matrix, where each block corresponds to a different genomic region or chromosome.

Further examinations reveal that in each row, non-zero values are also localized within a certain range – often grouped into a few clusters of consecutive columns – if nodes are indexed in a “near-topological order”. This is because almost all edges in sequence graphs are local, and it is reflected in the matrix as long as the ordering that defines node indices mostly preserves this locality – i.e. adjacent nodes appear together in the ordering. This observation can be seen clearly in the distance index constructed for a chain graph with distance constraints $(d_1,d_2)$ as shown in Figure 3.

The nodes of a graph are indexed along the rows of its adjacency matrix with the same order governing the columns as well. This induces a total order on the set of nodes in the graph, assigning each node a unique index. Different total orders yield different adjacency matrices. As will be shown in Section 4.2, our method greatly benefits if the total order on the node set maximizes the locality of the indices for adjacent nodes; i.e. adjacent nodes are assigned nearby indices by such ordering. It can be shown that this problem is equivalent to finding a sparse matrix with minimum bandwidth by permuting its rows and columns. The problem of finding such ordering has been proven to be NP-Complete [31].

Several heuristic algorithms have been proposed to minimize the bandwidth of sparse matrices by reordering rows and columns [10, 11, 8]. However, applying these methods to our problem presents some key limitations. First, although our method logically views sequence graphs as character graphs, in practice it works with string-labelled graphs for space efficiency. This requires preserving the sequential ordering of bases within nodes of a string-labelled graph, a constraint that existing heuristics do not necessarily meet when applied to the character graph representation. Breaking this intra-node ordering introduces memory overhead and fragmentation, impacting graph traversal performance in downstream tasks through reduced memory locality and increased cache misses. Second, when applied directly to string-labelled graphs while preserving node ordering, these methods are less effective at reducing bandwidth compared to their performance on character graphs.

To address these limitations, we propose a new heuristic algorithm inspired by the Cuthill–McKee (CM) algorithm [10, 11]. The new algorithm directly operates on string-labelled sequence graphs. In essence, the algorithm traverses the graph in a breadth-first search similar to the CM. However, instead of prioritizing visiting nodes by their degrees, it prioritizes nodes that have a predecessor with a lower topological sort order (Algorithm 1). Note that, in practice, sequence graphs, particularly variation graphs, are DAGs in the majority of relevant cases. If input sequence graphs are not acyclic, e.g. in de Bruijn graphs, we resort to a “semi-topological ordering” established by “dagifying” the graph, i.e. by running topological sort algorithm on the graph while the traversal algorithm ignores any edges forming a cycle.

Our experimental results (in Section 5) demonstrate that the ordering achieved by Algorithm 1 offers a very effective heuristic while assigning sequential indices to bases within each node. This can be justified by the fact that sequence graphs generally preserve the linear structure of the sequences from which they are constructed, and the variations in these sequences are often only local in the genomic coordinates, thereby affecting the graph topology locally. As a result, prioritizing sibling nodes for index assignment while following topological ordering tends to localize indices of adjacent nodes, thereby minimizing bandwidth.

DiVerG aims to exploit particular properties of the distance index matrix $𝒯$ mentioned in Section 3.1, as well as its sparsity, to reduce space and time requirements. Our method introduces a new sparse storage format, which will be explained subsequently. This format relies on transforming the standard CRS by replacing column index ranges with their lower and upper bounds to achieve high compression rate. Later, we will propose tailored algorithms to compute sparse matrix operations directly on this compressed data structure without incurring additional overhead for decompression.

It can be easily shown that rCRS can be constructed from sorted CRS in linear time with respect to the number of non-zero values. Accessing element $a_{ij}$ in $𝐀$ represented in rCRS reduces to determining whether $j$ lies within any range in row $i$. Since the ranges are disjoint and sorted, this operation can be performed using binary search which requires $O(\log (\widehat{z}))$ time in the worst case, where $\widehat{z}$ is the maximum size of $C_i^{\prime }$ across all rows $i$. Note that, unlike the CRS format, the size of $C^{\prime }$ is no longer equal to $\text{nnz}(𝐀)$.

The space complexity of the rCRS representation decisively depends on the distribution of non-zero values in the rows of the matrix. This representation can substantially compress the column indices array $C$ – as low as $O(n)$ compared to $O(nnz(A))$ in CRS – if the array contains long distinct ranges of consecutive integers. The limitation of this representation is that it can be twice as large when there is no stretch of consecutive indices in $C$. It is important to note that the compression rate can directly influence query time, as it is proportional to the compressed representation of non-zero values.

To mitigate this issue, we define a variant of rCRS, which is referred to as Asymmetrical Range CRS (aCRS). In the worst case scenario, aCRS format takes as much space as the standard CRS while keeping the same time complexity $O(\log (\widehat{z}))$ for the query operation. We elaborate on the aCRS definition and related analysis in Appendix B.

Although aCRS theoretically provides a more compact representation compared to rCRS, the experimental results (in Section 5) show that it occupies as much space as rCRS when storing distance index $𝒯$ for sequence graphs in practice. This is due to the fact that the ordering defined in Section 4.1, together with the topology of sequence graphs, rarely results in isolated column indices. Consequently, aCRS holds two integers per range, similar to rCRS. Considering this fact and the relative ease of implementing efficient sparse matrix algorithms with rCRS compared to aCRS, particularly on GPU, we base our implementation on the rCRS format. This can be observed in the index constructed for a chain graph (Figure 3), which often resembles the general structure of sequence graphs.

Given sparse matrices ${𝐀}_{m\times n}$ and ${𝐁}_{n\times s}$ where $m$, $n$, and $s$ are positive integers, SpGEMM is an algorithm that computes $𝐂=𝐀\cdot 𝐁$ using sparse representations of two input matrices. In this section, we introduce a new SpGEMM algorithm, called Range SpGEMM or rSpGEMM for short, computing matrix-matrix multiplication when input and output matrices are in the rCRS format. The algorithm is designed to scale well on both multi-core architectures (CPUs) as well as many-core architectures (GPUs).

Similar to most parallel SpGEMM algorithms, our rSpGEMM follows Gustavson’s algorithm [23]. Algorithm 2 illustrates its general structure, in which $A(i,:)$ refers to row $i$ of matrix $𝐀$. This notation can be generalized to $B(j,:)$ and $C(i,:)$ in Algorithm 2 specifying row $j$ and $i$ in $B$ and $C$, respectively. Note that Algorithm 2 only iterates over non-zero entries of $𝐀$ and 2 and 3]$𝐁$.

The algorithm computes $𝐂$ one row at a time. To compute $C(i,:)$, it iterates over non-zero values in row $A(i,:)$ and computes its contribution to $C(i,:)$ by multiplying it to non-zero values in row $B_j$; i.e. the term $a_{ij}\cdot B(j,:)$ in Line 3. This contribution is then accumulated with the partial results computed from previous iterations. The row accumulation can be simplified to $C(i,:)\leftarrow C(i,:)+B(j,:)$ in Boolean matrices as $a_{ij}$ is 1.

Since computing each row of the final matrix is independent of the others, $𝐂=𝐀\cdot 𝐁$ can be seen as multiple independent vector-matrix multiplication $C(i,:)=A(i,:)\cdot 𝐁$. For simplicity, we will only focus on the equivalent vector-matrix multiplication.

Given that the output matrix is also in the sparse format, knowing the number of non-zero values in each row of $𝐂$ is essential before storing the calculated $C_i$ in the final matrix. This issue is usually tackled using a two-phase approach: firstly, the symbolic phase, delineates the structure of $𝐂$, primarily corresponding to the computation of its row map array. Secondly, the numeric phase, carries out the actual computation of $𝐂$.

Similar to SpGEMM algorithms based on Gustavson’s, the fundamental aspect of rSpGEMM resides in three anchor points: the data structure employed for row accumulation, memory access pattern to minimize data transfer latency, and distribution of work among threads in hierarchical parallelism. In the following, we address each of these points.

In this section, we shift our focus to calculating matrix-matrix addition in rCRS format which is required for computing the distance index. We propose the Range Boolean Sparse Add (or rSpAdd) algorithm which computes matrix $𝐂=𝐀\vee 𝐁$ where all matrices are in rCRS format and $\vee$ is Boolean matrix addition defined in Definition 1.

Considering that non-zero values in each row of rCRS matrices are represented by ranges, each with two integers, and these integers are sorted and disjoint, row $C(i,:)$ is equivalent to the sorted combination of ranges in both rows $A(i,:)$ and $B(i,:)$, where overlapping ranges are collapsed into one. For each row in $𝐂$, this is achieved by using two pointers, each initially pointing to the first element of the respective rows. The algorithm peeks at the pairs in $A(i,:)$ and $B(i,:)$ indicated by the pointers. If two ranges are disjoint, it inserts the one with smaller bounds in $C(i,:)$ and increments the pointer indicating the inserted pair by two. If two ranges overlap, they will be merged by inserting the minimum of lower bounds and the maximum of upper bounds of the two pointed elements into $C(i,:)$ in that order. Then, both pointers are incremented to indicate the next pairs. In case either pointer reaches the end of the row, the remaining elements from the other row are directly appended to $C(i,:)$.

The time complexity of this algorithm is bound by compressed representation of input matrices in rCRS format. Since the ranges in rCRS are disjoint and sorted, merging two rows $A(i,:)$ and $B(i,:)$ can be done in $\mathrm{\Theta }(|A(i,:)|+|B(i,:)|)$; in which $|A(i,:)|$ and $|B(i,:)|$ are the sizes of $A(i,:)$ and $B(i,:)$ in rCRS format, respectively. rSpAdd not only requires smaller working space compared to SpAdd but is arguably faster due to compressed representation of the matrices.

As stated in Section 3.1, the distance index constructed by Equation 1 is a block-diagonal matrix. For this reason, DiVerG builds the distance index incrementally for each block; i.e. computing the distance matrix for each component of the sequence graph individually. The adjacency matrix for each component can be constructed in rCRS format one row at a time. Therefore, it is not required to store the whole matrix in sorted CRS format to be able to convert it to rCRS. Once matrices $𝐀$ and $𝐈$ are in rCRS format, $𝒯$ can be computed using rSpAdd and rSpGEMM. The powers of $𝐀$ and $𝐀\vee 𝐈$ are computed using “exponentiation by squaring” [28]. Ultimately, matrix $𝒯$ is further compressed by encoding both the column indices and row map arrays using Elias-$\delta$ encoding [17].

DiVerG provides an exact solution to the DVP through its compressed format for storing both final and intermediate matrices when computing the distance index $𝒯$. The compression ratio directly impacts the efficiency of the matrix multiplication operations, with higher compression resulting in faster construction. Since compression efficiency depends on the topology of the sequence graph and the structure of its adjacency matrix, establishing tight theoretical bounds on space and time complexities is challenging. However, we can derive meaningful lower bounds to facilitate comparison with previous approaches.

It can be demonstrated that for a sequence graph $G=(V,E)$, indexing cost is lower bounded by that of the chain graph $G_c=(V^{\prime },E^{\prime })$ representing the longest path in $G$ [27]. This analysis applies to our rSpGEMM algorithm as well. For a chain graph, the index $𝒯$ has a particularly simple structure (Figure 3): in each row $i$, all columns from $d_1$ to $d_2$ is 1. This structure offers two advantages in the rCRS format:

• $\mathrm{■}$

  Space complexity: Each row requires only two integers regardless of the distance range, resulting in $\mathrm{\Omega }(|V|)$ space complexity – independent of distance parameters $d_1$ and $d_2$, unlike PairG’s $\mathrm{\Omega }(|V|\cdot d)$ where $d=d_2-d_1$. Note that Elias encoding of the row map and column indices arrays compresses the final index even further in practice.

• $\mathrm{■}$

  Time complexity: While rSpGEMM performs the same operations as SpGEMM, bit-parallelism reduces the computation time by a factor of the word size $W$. Jain et al. [27] showed that computing $𝒯$ for $G$ using SpGEMM requires $\mathrm{\Omega }\left(|V^{\prime }|\left(d^2+\log d_1\right)\right)$; therefore, DiVerG’s rSpGEMM requires $\mathrm{\Omega }\left(|V^{\prime }|\left(d^2/W+\log d_1\right)\right)$. Note that, this improvement is particularly significant when computing ${\left(𝐀\vee 𝐈\right)}^{d_2-d_1}$, where multiple non-zero values per row enable effective bit-parallelism.

Given that sequence graphs often closely resemble their linear structure, these bounds provide a reasonable approximation of actual performance while highlighting the fundamental advantages of our approach.

DiVerG is implemented in C++17. While we develop our theory using character graphs, our implementation handles general string-labelled sequence graphs without loss of generality. The matrix $𝒯$, in the final stage, is converted to encoded rCRS format (described in Section 4.5) whose size is reported as index size in Section 5.3. Our implementation relies on Kokkos [13] for performance portability across different architectures, and we specifically focused on two execution spaces: OpenMP (CPU), and CUDA (GPU).

We compared the performance of our algorithm with PairG which uses the kkSpGEMM meta-algorithm implemented in Kokkos [13] for computing sparse matrix multiplication. Kokkos offers several algorithms for computing SpGEMM. The kkSpGEMM algorithm attempts to choose the best algorithm based on the inputs. As the PairG implementation was in practice not runnable on GPU, we re-implemented the corresponding part of PairG to support GPU and used this implementation for comparison.

All CPU experiments were conducted on an instance on de.NBI cloud with a 28-core (one thread per core) Intel® Xeon® Processor running Ubuntu 22.4. Distance index construction on GPU were carried out on an instance with an NVIDIA A40 GPU running Ubuntu 22.4.

DiVerG is assessed using seven different graphs. Four of these graphs were previously employed by PairG in their evaluation, while three additional graphs have been introduced to address significantly larger and more complex scenarios. Together, these graphs span over a spectrum of scales, complexities, and construction methods.

Three of these graphs are variation graphs of different regions of the human genome, constructed from small variants in the 1000 Genomes Project [1] based on GRCh37 reference assembly: the mitochondrial DNA (mtDNA) graph, the BRCA1 gene graph (BRC), and the killer cell immunoglobulin-like receptors (LRC_KIR) graph. One graph is a de Bruijn graph constructed using whole-genome sequences of 20 strains of B. anthracis. We added three new variation graphs to this ensemble: the MHC region graph constructed from alternate loci released with GRCh38 reference assembly [34], and the HPRC CHM13 graph [29] constructed by minigraph. Table 1 presents some statistics of the sequence graphs used in our evaluation.

We evaluate the index size (size), construction time (ctime), and high water mark memory footprint during index construction (mem) for different distance constraints using the graphs explained in Section 5.2. Moreover, the query time (qtime) is measured for the distance index by averaging the time spend for querying 1M randomly sampled paired positions throughout each graph. Our method is evaluated using two separate sets of distance constraints.