# Discreteness of Space from Anisotropic Spin-Orbit Interaction

Ahmed Farag Ali 1<EMAIL_ADDRESS>Barun Majumder 2
<EMAIL_ADDRESS>1Dept. of Physics, Faculty of Sciences, Benha
University, Benha 13518, Egypt.
2 University of Tennessee, Knoxville, USA

###### Abstract

Various approaches to Quantum Gravity suggest an existence of a minimal
measurable length. The cost to have such minimal length could be modified
uncertainty principle, modified dispersion relation, non-commutative geometry
or breaking of continuous Lorentz symmetry. In this paper, we propose that
minimal length can be obtained naturally through spin-orbit interaction. We
consider Dresselhaus anisotropic spin-orbit interaction as the perturbative
Hamiltonian. When applied to a particle, it implies that the space, which
seizes this particle, should be quantized in terms of units that depend on
particle’s mass. This suggests that all measurable lengths in the space are
quantized in units depending on existent mass and the Dresselhaus coupling
constant. On one side, this indicates a breakdown of the space continuum
picture near the scale of tabletop experiments, and on the other side, it
proposes that spin-orbit interaction is a possible quantum gravity effect at
low energy scale that leads naturally to space quantization.

Physicists have developed various approaches to quantum gravity, however still
today we do not have a conclusive evidence of quantum gravitational effects
from experiments. The quantization of space from various theories predicts the
existence of minimal length and its manifestation is widely predicted to exist
in the form of a generalized uncertainty principle a1 ; a2 ; a3 ; a4 ; b1 . A
more detailed review on the topic can be found in arev . The effects of
quantum gravity have been investigated in various experiments and
observations. For example, observations from gamma-ray bursts are still an
open avenue for the investigations test1 ; test2  but far from any conclusive
evidence test3 ; test4 . Another promising path is to perform quantum
mechanical tabletop experiments to look for deviations from standard quantum
mechanical results as a consequence of modifying the Heisenberg uncertainty
principle based on the predictions from a theory of quantum gravity a4 ; test5
. Recently there has been a suggestion that experiments with a quantum optical
oscillatory system can provide a measurement of the canonical commutator by
optical interferometric techniques. The set up is highly effective in
investigating the bounds on quantum gravitational parameters in the
generalized uncertainty principle test6 . Bekenstein test7 ; test8  proposed a
tabletop experiment arguing that the present ultrahigh vacuum and cryogenic
technology may have enough sensitivity to detect signatures of quantum
gravity. The main objective of his table top experiment was to prove that at
length scales of the order of Planck length the structure of space is not a
smooth manifold. In other words, insights on the discreteness of space at
Planck scale may be drawn from state-of-art tabletop experimental setups.

If we only consider elementary particles as probes for quantum gravity we
would need to accelerate the particles which corresponds to very high energies
which is not feasible even within the realm of future particle accelerators.
As an alternative approach to investigate the nature of quantum gravity we
must try to relate our understanding to the present quantum laboratory
experiments. In this paper we will discuss a spin-orbit coupling perturbative
Hamiltonian that is reponsible for splitting of energy bands in condense
matter systems with broken inversion symmetry commonly known as the
Dresselhaus effect d1 . Later we show that if we consider such an Hamiltonian
as a perturbation (whose effect has already been observed in condensed matter
systems) and quantize the system with respect to Schr$\ddot{o}$dinger
equation, we essentially arrive at the _discreteness of space_. A recent study
has shown that a similar spin-orbit coupling whose strength is of the same
order as that of Dresselhaus coupling can generate a Lorentz violation in a
non-commutative theory epl1 .

Spin-orbit coupling is an artifact of special theory of relativity where the
electric and magnetic fields of a moving electron in its own reference frame
can undergo transformations and interact with it’s own spin to break spin
state degeneracies. Rashba d2  and Dresselhaus d1  spin-orbit coupling are
very important in the field of spintronics which in terms of strength are
similar and they are realized in systems with broken inversion symmetry. In
semiconductor quantum well systems it has been found that if the strengths of
the Rashba and Dresselhaus interactions are equal then SU(2) symmetry is
realized, however the SU(2) is broken in the presence of a cubic Dresselhaus
coupling well1 ; well2 .

Let us consider the minimal Hamiltonian ($H_{RD}$) with Rashba and Dresselhaus
coupling up to third-order in momentum with all the symmetry allowed terms in
the $p_{z}=0$ plane m1

$H_{RD}=H_{0}+\frac{1}{\hbar^{2}}\alpha(p)p_{y}\sigma_{z}~{}~{},$ (1)

where $p_{y}$, $p_{z}$ are the cartesian components of the momentum and
${\bf\sigma}$ are the Pauli matrices for the spin degree of freedom. The
Hamiltonian can be derived from a ${\bf k\cdot p}$ perturbation theory. The
free particle parabolic part of the Hamiltonian is

$H_{0}=\frac{p_{y}^{2}}{2m_{y}^{*}}~{}~{},$ (2)

where $m_{y}^{*}$ is the effective mass along orthogonal direction $p_{y}$.
The pure linear and cubic Rashba-Dresselhaus limit is recovered when

$\alpha(p)=\alpha_{l}+\frac{\alpha_{c}}{\hbar^{2}}p_{y}^{2}~{}~{},$ (3)

where $\alpha_{l}$ is the linear spin-orbit coupling constant and $\alpha_{c}$
account for p-cubic anisotropic interactions. The experimental values of
$\alpha_{l}$ and $\alpha_{c}$ can be found in TABLE I of m1 . For a purely
cubic Dresselhaus spin-polarization pattern let us consider only the
anisotropic coupling term ($\alpha_{l}=0$). The Hamiltonian of Eq.(1) can be
written as

$H_{RD}=\frac{p_{y}^{2}}{2m_{y}^{*}}+\frac{\alpha_{c}}{\hbar^{4}}p_{y}^{3}\sigma_{z}~{}~{}.$
(4)

From the study in b1  we know that a cubic order term in the Hamiltonian is
crucial for achieving discreteness of space however we need to consider
corrections to the uncertainty principle. In this present study we consider an
interaction found naturally in many condense matter systems with a cubic term
in momentum to get discreteness of space without any other assumptions. Let us
study a quantum particle in a box problem with such an Hamiltonian of Eq.(4)
considering the cubic Dresselhaus anisotropic term as a perturbation. With the
diagonal Pauli matrix $\sigma_{z}$ we have to solve a Schr$\ddot{o}$dinger
equation of the form

$\left[-\frac{\hbar^{2}}{2m_{p}}\frac{\partial^{2}}{\partial
x^{2}}+\frac{i\alpha_{c}}{\hbar}\frac{\partial^{3}}{\partial
x^{3}}\right]\psi=E\psi~{}~{},$ (5)

from $x=0$ to $L$ with boundary conditions $\psi(0)=\psi(L)=0$. A general
solution to the equation in leading order in $\alpha_{c}$ can be written as

$\psi=Ae^{ik(1+k\alpha_{c}m_{p}/\hbar^{3})x}+Be^{-ik(1-k\alpha_{c}m_{p}/\hbar^{3})x}+Ce^{\frac{-ix\hbar^{3}}{2\alpha_{c}m_{p}}}~{}~{},$
(6)

where $k=\sqrt{2m_{p}E/\hbar^{2}}$ and $A$, $B$ and $C$ are arbitrary
integration constants to be fixed by boundary conditions. With $\alpha_{c}=0$
we must get the energy quantization in the form
$E_{n}=\frac{n^{2}\pi^{2}\hbar^{2}}{2m_{p}L^{2}}$. This requires for
$\alpha_{c}\rightarrow 0$, $|C|\rightarrow 0$. The boundary condition
$\psi(0)=0$ gives the condition $B=-A-C$. This along with the condition
$\psi(L)=0$ gives an equation

$\displaystyle
2A\sin(kL)\sin\left(\frac{\alpha_{c}m_{p}k^{2}L}{\hbar^{3}}\right)$
$\displaystyle=2iA\sin(kL)\cos\left(\frac{\alpha_{c}m_{p}k^{2}L}{\hbar^{3}}\right)$
$\displaystyle-Ce^{-ikL(1-k\alpha_{c}m_{p}/\hbar^{3})}$
$\displaystyle+Ce^{-iL\hbar^{3}/2\alpha_{c}m_{p}}$ (7)

If we separate the real and imaginary parts then we get the following
equations:

$\displaystyle
2A\sin(kL)\sin\left(\frac{\alpha_{c}m_{p}k^{2}L}{\hbar^{3}}\right)$
$\displaystyle+C\cos\left(kL(1-k\alpha_{c}m_{p}/\hbar^{3})\right)$
$\displaystyle=C\cos\left(\frac{L\hbar^{3}}{2\alpha_{c}m_{p}}\right)$ (8)

and

$\displaystyle
2A\sin(kL)\cos\left(\frac{\alpha_{c}m_{p}k^{2}L}{\hbar^{3}}\right)$
$\displaystyle+C\sin\left(kL(1-k\alpha_{c}m_{p}/\hbar^{3})\right)$
$\displaystyle=C\sin\left(\frac{L\hbar^{3}}{2\alpha_{c}m_{p}}\right)$ (9)

As mentioned earlier when $\alpha_{c}\rightarrow 0$, $|C|\rightarrow 0$ which
is the usual energy quantization for a quantum particle in a box and we get
the condition $kL=q\pi~{}(q\in\mathbb{N})$. But when $\alpha_{c}\neq 0$ (its
experimental value can be found in m1 ), $kL=q\pi$ must still hold. As
$|C|\neq 0$ when $\alpha_{c}\neq 0$ the following additional relation must be
satisfied to satisfy the boundary conditions given by Eqns.(8) and (9):

$L=q\pi\frac{2\alpha_{c}m_{p}}{\hbar^{3}}+{\cal
O}(\alpha_{c}^{2})~{}~{}~{}~{}~{}\text{where}~{}~{}~{}~{}~{}q\in\mathbb{N}.$
(10)

So in order to contain a particle in the box, length must be quantized
according to Eq.(10). Observe that the length is quantized in terms of the
mass of the particle ($m_{p}$). This indicates a dynamical relation between
particle’s mass and space quantization. If the space has a discrete nature,
its discreteness unit will depend on the mass contained in this space. The
discreteness of space in that sense behaves like curvature of space. Both are
varying with mass. Similar results of length quantization have been obtained
with the generalized uncertainty principle b1  and are also found in Loop
Quantum Gravity lqg1 . However in those approaches the length quantization
comes in units of Plank length ($l_{PL}$). Here we show that a perturbative
spin-orbit interaction with broken SU(2) symmetry can automatically give
discreteness of space in terms of the mass of the particle without the need of
any extra assumptions. This indicates that Spin-Orbit interaction which is a
natural interaction in many condense matter systems is a possible quantum
gravity effect. Even though studied in a different context, our result is
analogous to the area quantization which is evident in loop quantum gravity
where $\Delta A\propto M^{2}$, where $A$ is the area of event horizon and $M$
is the black hole mass hod . Also the mass of a Schwarzschild black hole can
be shown to be quantized ma1 ; ma2 ; ma3  which agrees with Bekenstein’s
argument bek1 . The constant associated in loop quantum gravity is the Immirzi
parameter ash1  and is our case a natural anisotropic coupling $\alpha_{c}$.
More careful investigation is needed on this point. It will be also
interesting to investigate such discreteness in two and three dimensions to
look at the discreteness of area and volume. It will be interesting to study
the phenomenological implications of such discreteness of space in tabletop
experiments. We hope to report on these in the future.

## References

  * (1) D. Amati, M. Ciafaloni, G. Veneziano, Phys. Lett. B 216 (1989) 41.
  * (2) M. Maggiore, Phys. Lett. B 304 (1993) 65.
  * (3) J. Magueijo, L. Smolin, Phys. Rev. Lett. 88 (2002) 190403.
  * (4) S. Das, E. C. Vagenas, Phys. Rev. Lett. 101 (2008) 221301.
  * (5) A. F. Ali, S. Das and E. C. Vagenas, Phys. Lett. B 678 (2009) 497.
  * (6) S. Hossenfelder, Living Rev. Relativity 16 (2013) 2.
  * (7) Amelino-Camelia, G., Ellis, J., Mavromatos, N. E., Nanopoulos, D. V. and Sarkar, S., Nature 393 (1998) 763-765.
  * (8) Jacob, U. and Piran, T., Nature Phys. 7 (2007) 87-90.
  * (9) Abdo, A. A. et al., Nature 462 (2009) 331-334.
  * (10) Tamburini, F., Cuofano, C., Della Valle, M. and Gilmozzi, R., Astron. Astrophys. 533 A71 (2011).
  * (11) Ali, A. F., Das, S. and Vagenas, E. C., Phys. Rev. D 84 (2011) 044013.
  * (12) Pikovski, I., Vanner, M. R., Aspelmeyer, M., Kim, M.S. and Brukner, C., Nature Phys. 8 (2012) 393-397.
  * (13) J. D. Bekenstein, Phys. Rev. D 86 (2012) 124040.
  * (14) J. D. Bekenstein, Found Phys 44 (2014) 452-462.
  * (15) G. Dresselhaus, Phys. Rev. 100 (1955) 580.
  * (16) S. Aghababaei and G. Rezaei, Euro. Phys. Lett. 132 (2020) 11002.
  * (17) E. I. Rashba, Sov. Phys. Solid State 2 (1960) 1109.
  * (18) J. D. Koralek et. al., Nature 458 (2009) 610-613.
  * (19) Lorenz Meier et. al., Nature Physics 3 (2007) 650-654.
  * (20) C. Autieri, P. Barone, J. Slawinska and S. Picozzi, Phys. Rev. Materials 3 (2019) 084416.
  * (21) T. Thiemann, J. Math. Phys. 39 (1998) 3372.
  * (22) S. Hod, Phys. Rev. Lett. 81 (1998), 4293; S. Hod, Gen. Rel. Grav. 31 (1999), 1639.
  * (23) J. Makela, P. Repo, Phys. Rev. D 57 (1998) 4899.
  * (24) O. Obregon, M. Sabido, V.I. Tkach, Gen. Relativ. Gravit. 33 (2001) 913.
  * (25) B. Majumder, Phys. Lett. B 701 (2011) 384.
  * (26) J.D. Bekenstein, Phys. Rev. D 7 (1973) 2333; J.D. Bekenstein, Phys. Rev. D 9 (1974) 3292;
  * (27) A. Ashtekar, J. Baez, A. Corichi and K. Krasnov, Phys. Rev. Lett. 80 (1998) 904.


# Distributed-Memory Parallel Contig Generation for _De Novo_ Long-Read Genome
Assembly

Giulia Guidi  University of California, BerkeleyLawrence Berkeley National
LaboratoryBerkeleyCaliforniaUSA ,  Gabriel Raulet  Lawrence Berkeley National
LaboratoryBerkeleyCaliforniaUSA ,  Daniel Rokhsar  University of California,
BerkeleyDOE Joint Genome InstituteBerkeleyCaliforniaUSA ,  Leonid Oliker
Lawrence Berkeley National LaboratoryBerkeleyCaliforniaUSA ,  Katherine Yelick
University of California, BerkeleyLawrence Berkeley National
LaboratoryBerkeleyCaliforniaUSA and  Aydın Buluç  University of California,
BerkeleyLawrence Berkeley National LaboratoryBerkeleyCaliforniaUSA

(2022)

###### Abstract.

_De novo_ genome assembly, i.e., rebuilding the sequence of an unknown genome
from redundant and erroneous short sequences, is a key but computationally
intensive step in many genomics pipelines. The exponential growth of genomic
data is increasing the computational demand and requires scalable, high-
performance approaches.

In this work, we present a novel distributed memory algorithm that, from a
string graph representation of the genome and using sparse matrices, generates
the contig set, i.e., overlapping sequences that form a map representing a
region of a chromosome.

Using matrix abstraction, we mask branches in the string graph, and compute
the connected component to group genomic sequences that belong to the same
linear chain (i.e., contig). Then, we perform multiway number partitioning to
minimize the load imbalance in local assembly, i.e., concatenation of
sequences from a given contig. Based on the assignment obtained by
partitioning, we compute the induce subgraph function to redistribute
sequences between processes, resulting in a set of local sparse matrices.
Finally, we traverse each matrix using depth-first search to concatenate
sequences.

Our algorithm shows good scaling with parallel efficiency up to 80% on 128
nodes, resulting in uniform genome coverage and showing promising results in
terms of assembly quality. Our contig generation algorithm localizes the
assembly process to significantly reduce the amount of computation spent on
this step. Our work is a step forward for efficient _de novo_ long read
assembly of large genomes in a distributed memory.

††journalyear: 2022††copyright: rightsretained††conference: 51st International
Conference on Parallel Processing; August 29-September 1, 2022; Bordeaux,
France††booktitle: 51st International Conference on Parallel Processing (ICPP
’22), August 29-September 1, 2022, Bordeaux, France††doi:
10.1145/3545008.3545050††isbn: 978-1-4503-9733-9/22/08

##  1\. Introduction

Genomics is facing exponential growth in data coming from improved and cheaper
instrumentation that is overwhelming conventional analytical infrastructures.
It is especially true for _de novo_ genome assembly (Zhang et al., 2011), that
is decoding the sequence of an unknown genome from redundant and erroneous
short sequences. A key but computationally intensive step in many genomics
pipelines.

On the one hand, large-scale parallel systems promise to overcome the
computational and memory limitation of conventional shared-memory computing.
On the other hand, programming such large-scale systems is challenging,
especially when computation is highly irregular, as is often the case in
genomics.

The advent of long–read sequencing technologies (Eid et al., 2009; Goodwin et
al., 2015), which yield sequences with an average length of more than 10,000
base pairs (bp), improved our ability to assemble complex genomic repeats and
obtain more accurate assemblies that were not possible with short–read
technologies (Phillippy et al., 2008; Nagarajan and Pop, 2009). In the past,
longer sequences were associated with higher error rates (Hon et al., 2020),
resulting in significantly higher algorithm complexity and computational
costs. Recent advances in sequencing technologies have lowered error rates,
reducing the need for error correction and polishing phases, but computational
costs remain high as we generate more data and reconstruct higher quality
assemblies.

It is common practice to assemble long-read data according to the
Overlap–Layout–Consensus (OLC) paradigm (Berlin et al., 2015). The first step
(O) is to identify overlaps between reads to create an _overlap_ graph. The
second step (L) simplifies the overlap graph by removing redundant information
and transforming it into an _string_ graph. The third step (C) takes the
string graph and processes it to obtain a first draft of the unknown genome,
returning a contig set as output. A _contig_ is a group of overlapping DNA
sequences that together form a region of a chromosome. A common approach to
contig generation takes the string graph as input, branching nodes in this
graph are masked, and the set of linear, unbranched paths in the graph is
extracted to form the contig set. This step can be followed by further
polishing phases to merge and correct the contig set and to separate
haplotypes.

In this paper, we build on previous work in the literature, diBELLA 2D (Guidi
et al., 2021c, b), and present a novel distributed-memory algorithm that
generates the contig set starting from a string graph representation of the
genome and using a sparse matrix abstraction. The proposed approaches are
implemented in a software package called ELBA, which also includes the overlap
and layout phases of diBELLA 2D, which can scale to hundreds of nodes. Each
stage of the computation in ELBA is implemented for distributed memory.

The algorithm first determines the contig set starting from a string graph by
masking out the branches and extracting the non-branched paths by computing
connected components over the graph (i.e., a sparse matrix). Using a greedy
multiway number partitioning algorithm, ELBA then determines how to
redistribute sequences (i.e., rows and columns of the sparse matrix) between
processes so that sequences belonging to the same contig are stored on the
same processor. Once the sequences-to-processor assignment is computed, ELBA
uses this information to implement the induced subgraph function that
redistributes the sequences using sparse matrix abstraction so that each
processor has a local sparse matrix representing the non-branching contig(s)
assigned to the processor. ELBA then computes a local assembly step on each
processor independently and in parallel. In doing so, it uses depth-first
search to concatenate the sequences belonging to a particular contig and
return the contig set. ELBA is the first distributed-memory implementation of
the OLC assembly paradigm.

Our algorithm shows good scaling with parallel efficiency up to 80% on 128
nodes on a representative dataset. Our approach leads to uniform coverage of
the genome and shows promising results in terms of assembly quality. Our
contig generation algorithm is fast and efficient because it localizes the
graph traversal problem, i.e., the practical concatenation of sequences, so
that most of the computation can be performed locally without requiring fine-
grained communication. Compared to related work, this approach significantly
reduces the fraction of the runtime required by the entire assembly pipeline
for contig generation. Our implementation is publicly available at
https://github.com/PASSIONLab/ELBA.

##  2\. Background

A genome consists of one or more linear chains of DNA molecules, i.e.,
nucleotides or bases, represented by the alphabet $\Sigma=\\{\rm A,C,G,$
$T\\}$. Each linear chain has a direction and is paired with a second chain
called its _reverse complement_ , which has the opposite direction. Together
they are called _strands_ and wind around each other to form a double helix;
this double helix is organized in three-dimensional space and forms a
chromosome. Different species may have different numbers of chromosomes. On
the opposite strands, A always pairs with T and C pairs with G. Given a string
$v=\rm ATTCG$, its reserve complement is $v^{\prime}=\rm CGAAT$.

The goal of _de novo_ long-read genome assembly is to reconstruct an unknown
genome sequence, possibly consisting of multiple chromosomes, starting from a
collection of redundant, erroneous, and significantly shorter sequences called
_reads_. For long-read sequencing data, this is commonly done following the
OLC assembly paradigm introduced in the previous section. In this paper, we
take the overlap (O) and layout (L) phases as implemented in diBELLA 2D (Guidi
et al., 2021c, b) and extend them to include the third (C) phase, that is
contig generation. diBELLA 2D uses a $k$-mer–based approach and relies on
parallel sparse matrix multiplication for overlap detection and transitive
reduction to obtain a string graph (i.e., the input to the contig generation
algorithm described in Section 4).

diBELLA 2D sees the overlap detection step as a sparse matrix multiplication
between a
$\lvert\textit{sequences}\rvert$-by-$\lvert\textit{{k}-mer}{s}\rvert$ matrix
$\mathbf{A}$ and its transpose $\mathbf{A}^{\mbox{\scriptsize\sf T}}$. The
resulting matrix $\mathbf{C}=\mathbf{A}\mathbf{A}^{\mbox{\scriptsize\sf T}}$
is a sparse
$\lvert\textit{sequences}\rvert$-by-$\lvert\textit{sequences}\rvert$ that
stores information about which sequences share at least one $k$-mer. diBELLA
2D then computes an element-wise pairwise alignment step for each non-zero in
the $\mathbf{C}$ matrix, removing non-zeros that are below a certain alignment
threshold. The result of this operation is the result matrix $\mathbf{R}$,
which in turn is the input for the transitive reduction phase, whose goal is
to remove redundancies from the overlap graph and return a string graph in the
form of a sparse matrix $\mathbf{S}$.

A string graph (or matrix) is a graph $G=(V,E)$, where $V$ is the set of
sequences and $E$ is the set of overlap _suffixes_ between any two vertices. A
string graph contains neither redundant edges nor redundant vertices. A
redundant vertex is a read that is completely contained (i.e., aligned) within
another read and whose edges therefore contain information that the containing
read already has. A redundant edge or transitive edge is an edge that carries
less or the same information as a parallel path in the graph and therefore can
be removed without loss of information during transitive reduction.

In diBELLA 2D, so also in ELBA, the string graph is encoded as a bidirected
graph, so we can encode information about both DNA strands. A bidirected graph
(Edmonds and Johnson, 2003) is one in which each edge at each end has an
independent orientation (or arrow). There are three types of bidirected edges:
(a) those with arrows at both ends pointing outward, (b) those with both
arrows pointing inward, and (c) those with both arrows pointing in the same
direction, one away from its vertex and the other toward its vertex.

Given two vertices $v_{i}$ and $v_{j}$, the edge $e_{i,j}$ between them
defines an overlap between the two vertices (i.e., sequences). The weight
$w_{i,j}$ of the edge $e_{i,j}$ stores the _overhang length_ or suffix length,
i.e., the part of $v_{j}$ that goes beyond the overlap between $v_{i}$ and
$v_{j}$, and its direction, i.e., how such sequences overlap with each other (
for example, $v_{i}$ and $v_{j}$ could have the same strand orientation and
the final portion of $v_{i}$ could overlap with the initial portion of
$v_{j}$).

If $G=(V,E)$ is a bidirected string graph, then a valid walk in $G$ is a
continuous sequence of edges where each vertex is entered by an inward head
and exited by an outward head (unless it is the end of the path) or vice
versa.

A contig is a collection of overlapping DNA sequences that together form a
consensus region of DNA, such as a chromosome. To be considered a contig, a
set of overlapping sequences must form a linear chain that is a valid walk.

##  3\. Related Work

Contig generation is an important step in any _de novo_ genome assembly
pipeline, regardless of the type of sequencing technology (long-read or short-
read). In this section, we summarize the literature on this kernel and explain
how it relates to our work.

The paradigm for assembling long-read sequencing data is composed of three
main stages: finding overlapping sequences to create an overlap graph,
removing redundant information to simplify the graph, and generating the
contig set.

HiCanu (Nurk et al., 2020) and Falcon-Unzip (Chin et al., 2016) implement
similar approaches inspired by the Bogart algorithm presented by Canu (Koren
et al., 2017) to generate the contig set in shared memory starting from a
sparse long-read overlap graph or a string graph. The Bogart module creates an
assembly graph using a variant of the _best overlap graph_ strategy of Miller
et al (Miller et al., 2008), where a _best_ overlap is the longest overlap to
a given read end excluding _contained_ sequences (i.e., when all bases in one
sequence are aligned to another sequence). The Bogart algorithm removes
overlapping sequences from the overlap graph to include only those that are
within some tolerance of the global median error rate, and recalculates the
longest overlapping sequences using only that subset (i.e., a sparse overlap
graph). Bogart generates the initial contig set from the maximum non-branching
paths in that graph.

The Shasta assembler (Shafin et al., 2020) also uses a similar procedure by
creating an undirected graph where each vertex is an oriented read (i.e., each
read contributes two vertices to the read graph, one in its original
orientation and one in the reverse complement orientation) and an undirected
edge is created between two vertices when we find an alignment between the
corresponding oriented sequences. To reduce the high connectivity in the
repeat regions, the Shasta assembler preserves only the $k$-nearest neighbor
subset of the edges. The contig set is created from the linear structures of
the graph.

In contrast, Hifiasm (Cheng et al., 2021) generates a primary assembly based
on the topological structures of the graph and the phasing relationship
between the different haplotypes using a bubble-popping procedure (Li, 2016).
Finally, a best overlap graph is used to deal with remaining unresolved
substructures in the assembly graph.

A De Bruijn graph-based approach is common (Medvedev et al., 2007) for genome
assembly of short-read sequencing data, but it has not been suitable for long-
read data in the past because of higher error rates.

PaKman (Ghosh et al., 2020) is a parallel distributed memory contig generation
algorithm for short-read sequencing technology. PaKman introduces a new
compact data representation of the De Bruijn graph, named PaK-Graph, which
uses iterative compression to fit the graph into the memory available on each
node. This compression enables low-cost replication of the graph across nodes,
reducing the need for communication and creating an embarrassingly parallel
procedure. PaKman also introduces an algorithm to perform non-redundant contig
generation, i.e., to avoid two processes traversing the same path in the graph
and generating the same contig.

MetaHipMer (Georganas et al., 2018; Hofmeyr et al., 2020) and the earlier
HipMer (Georganas et al., 2015) are distributed-memory _de novo_ metagenome
and genome assembly pipelines, respectively, designed for short-read data and
thus also use the De Bruijn-based approach. Both are implemented using a
partitioned global address space model in either UPC (Carlson et al., 1999) or
UPC++ (Bachan et al., 2019). Contig generation is performed after the
construction of a distributed hash table of $k$-mers with the left and right
base extension from the input data stored with each $k$-mer. Each process then
creates contigs by starting at a $k$-mer, walking left and right, appending
the extension, and looking up the resulting $k$-mer in the hash table. These
lookups often take place on remote nodes and are performed with either a
remote memory operation or a remote procedure call. If a previously visited
$k$-mer is reached during this process, the two contigs are merged. Fine-
grained synchronization prevents a data race that can occur when two processes
attempt to merge at the same time.

As error rates decrease for long-read data, we find in the state of the art
some attempt to use a De Bruijn graph approach for this type of sequencing
technology, such as the shared-memory assembler Flye (Kolmogorov et al.,
2019). Instead of generating a contig set, Flye first generates a _disjointig_
set, i.e. concatenating multiple disjoint genomic sequences, and then
concatenates these error-prone disjointigs into a single string (in any
order), constructs an assembly graph from the resulting concatenation, uses
sequences to disentangle this graph, and resolves bridged repetitive areas
(which are bridged by some sequences in the repeat graph). It then uses the
repeat graph to resolve unbridged repetitive areas ( that are not bridged by
sequences) based on the differences between repeat copies. The output is a
contig set generated from the paths in this graph. Using a De Brujin approach
in conjunction with long read sequencing data leads to more complex algorithms
than those commonly used in an OLC-based assembler. This makes an approach
such as Flye’s unsuitable for scalable implementation on distributed memory
architectures, as the access patterns are extremely fine-grained.

Figure 1. A high-level representation of ELBA built on diBELLA 2D (Guidi et
al., 2021c) and extended to include contig generation, which is critical to
the assembly process.

Regarding the use of sparse matrices and matrices in general for genomic
computations, there is related work for both shared and distributed memory
machines. The shared-memory software BELLA (Guidi et al., 2021b) is the first
work to propose the use of sparse matrices in the context of _de novo_ genome
composition, focusing on the overlap detection phase. A sparse matrix
$\mathbf{A}$ is used to indicate the presence of $k$-mers in sequences, and by
multiplying by their transpose, i.e.
$\mathbf{A}\mathbf{A}^{\mbox{\scriptsize\sf T}}$, BELLA identifies overlapping
sequences.

diBELLA 2D (Guidi et al., 2021c) resembles BELLA, in that it computes both
overlap detection and transitive reduction over an overlap graph as
distributed SpGEMM and sparse computation. PASTIS (Selvitopi et al., 2020),
similarly, computes protein homology search as distributed SpGEMM.

Besta et al. (Besta et al., 2020) present an approach, similar to BELLA, for
calculating Jaccard similarity between sequences of different genomes using
distributed SpGEMM. The main difference is that their software is optimized
for the case where the output matrix is dense.

##  4\. Algorithm

In this section, we introduce ELBA, a parallel distributed memory algorithm
designed to generate a contig set from a string graph representation of long
read sequencing data.

ELBA uses sparse matrices as the core data structure throughout the
computation. This abstraction leads to a flexible and modular software design,
as well as high parallelization that performs well in distributed memory. To
implement parallel computation over sparse matrices, we use the Combinatorial
BLAS (CombBLAS) library (Azad et al., 2022), a well-known framework for
implementing graph algorithms in the language of linear algebra, and use a
semiring abstraction to overload the classical multiplication and addition
operation as needed. ELBA is inspired by Canu’s best-overlap graph strategy
(Koren et al., 2017) and improves upon the transitive reduction algorithm
presented in diBELLA 2D (Guidi et al., 2021c) by extracting a set of high-
quality, non-branching paths from the string graph.

Unlike MetaHipMer, which uses fine-grained communication at nearly every step
of contig construction, ELBA localizes the graph traversal problem so that the
sequences that make up a contig are stored locally at each rank, avoiding
communication.

###  4.1. Overview

Algorithm 1 Parallel matrix-based computation in ELBA.

1:procedure ELBA

2: $\textit{sequences}\leftarrow\textsc{FastaReader}($)

3: $\textit{{k}-mer}{s}\leftarrow\textsc{KmerCounter}($)

4: $\mathbf{A}\hskip
4.29993pt\leftarrow\textsc{GenerateA}(\textit{sequences},\hskip
0.82504pt\textit{{k}-mer}{s})$

5: $\mathbf{A}^{\mbox{\scriptsize\sf
T}}\leftarrow\textsc{Transpose}(\mathbf{A})$

6: $\mathbf{C}\leftarrow\mathbf{A}\mathbf{A}^{\mbox{\scriptsize\sf
T}}$$\triangleright$ Candidate overlap matrix

7: $\mathbf{C}\leftarrow\textsc{Apply}(\mathbf{C},\rm
Alignment())$$\triangleright$ Run alignment

8: $\mathbf{R}\leftarrow\textsc{Prune}(\mathbf{C},\rm
AlignmentScoreLessThan(\textit{t}))$

9: $\mathbf{R}\leftarrow\textsc{Prune}(\mathbf{R},\rm
IsContainedRead(\textit{}))$

10: $\mathbf{S}\leftarrow\textsc{TransitiveReduction}(\mathbf{R})$

11:
$\textit{cset}\leftarrow\textsc{ContigGeneration}(\mathbf{S})$$\triangleright$
Algorithm 2

12: return cset

Algorithm 2 Parallel contig generation on $\mathbf{S}$.

1:procedure ContigGeneration($\mathbf{S},sequences$)

2: $\mathbf{L}\leftarrow\textsc{BranchRemoval}(\mathbf{S})$

3: $\mathbf{v}\leftarrow\textsc{ConnectedComponent}(\mathbf{L})$

4: $\mathbf{p}\leftarrow\textsc{GreedyPartitioning}(\mathbf{v},P)$

5: $\mathbf{P}\leftarrow\textsc{InducedSubgraph}(\mathbf{L},\mathbf{p})$

6: $\textit{cset}\leftarrow\textsc{LocalAssembly}(\mathbf{P},sequences)$

7: return cset

ELBA uses a highly parallelizable strategy to generate the contig set in a
distributed memory environment. This means that each processor may not have
access to the entire set of sequences it needs for contig generation because
they are distributed among processes. Read sequences must therefore be
communicated across the processor grid before contigs are generated via depth-
first search.

To efficiently communicate sequences where they are needed, ELBA uses a
parallel algorithm, implemented as a sparse matrix-based computation over the
string matrix, to extract information about which sequences belong to the same
contig (i.e., a _linear component_) and therefore need to be sent to the same
processor. A linear component of a graph is the maximal subgraph where each
vertex is connected to two adjacent vertices, except for two vertices that
mark the extremities of the linear chain, which are connected to a vertex.
Once the membership of each sequence is known (i.e., to which contig it
belongs), ELBA uses this information to achieve load balancing that ensures
that each processor has similar amount of work during the local assembly
process.

This is achieved by first estimating the contig sizes, i.e., the number of
reads belonging to a particular contig, in parallel and then assigning
approximately equal collections of contig sets to each process. Then, the
sequences are redistributed among the processes using a newly implemented
function that is also based on sparse matrix computation. This function
generates local matrices from a distributed matrix (in our case, the string
matrix) on each process according to the load-balanced assignment. ELBA then
assembles the contigs whose sequence connectivity is stored in the local
matrix in parallel on each processor without requiring any further
communication.

Algorithm 1 summarizes the entire computation, including overlap detection and
transitive reduction—lines 11-12 represent the main contribution of this work,
better represented in Algorithm 2: (a) how the contig set is determined (line
2), (b) how the workload is distributed among the processes (lines 3-4), and
(c) how the local assembly is performed (line 5) and the contig set is output.

###  4.2. Contig Set Determination

The first step in our contig generation algorithm is to determine which
sequences belong to which contig based on the connectivity of the string
graph.

Let us define a linear component of the graph as the maximal subgraph
$L\subseteq R=(V,E)$ where (a) each vertex has either degree one or two and
(b) the subgraph is connected. In a linear component there are exactly two
vertices which can have degree one due to the connectedness rule. These
vertices are the endpoint vertices.

Then, we define the _contig set_ of $R$ as the set of all linear chains which
are not themselves subgraphs of another linear chain in $R$. Thus, the contig
set consists of several independent collections of nodes in the bidirected
string graph, where each collection is a linear sequence of unique overlapping
sequences. The determination of the contig set can be done in two stages.

First we have to identify _branching vertices_ that are vertices whose degree
is $\geq$ 3\. These vertices would make it impossible to determine a unique
linear chain. Therefore, our goal is to mask out such branching vertices to
obtain only linear sequences of vertices. For example, consider a string graph
consisting of (a) $v_{1}\to v_{2}\to v_{3}$, (b) $v_{3}\to v_{4}\to v_{5}\to
v_{6}$, and (c) $v_{3}\to v_{7}\to v_{8}$, the vertex $v_{3}$ is a branching
vertex since its degree is equal to three, and it would not allow us to
extract a linear chain, since at this point in the computation we cannot know
whether the correct linear path is (i) $v_{1}\to v_{2}\to v_{3}\to v_{4}\to
v_{5}\to v_{6}$ or (ii) $v_{1}\to v_{2}\to v_{3}\to v_{7}\to v_{8}$. Thus, we
mask out $v_{3}$ to obtain three linear sequences of vertices: (i) $v_{1}\to
v_{2}$, (ii) $v_{4}\to v_{5}\to v_{6}$, and (iii) $v_{7}\to v_{8}$.

To identify branching vertices, we perform a summation reduction over the row
dimension of the string matrix (i.e., an adjacency matrix) and return a
distributed vector $\mathbf{d}$ whose values represent the degree of the
corresponding row sequence (i.e., the index of the vector). Then we perform an
element-wise selection operation on the degree vector $\mathbf{d}$ to extract
the indices (sequences) whose value is greater than or equal to 3. The result
is a distributed vector $\mathbf{b}$ whose values are the indices of the
branching vertices. The branch vector $\mathbf{b}$ is used to remove the
corresponding rows and columns from the string matrix $\mathbf{S}$ and create
a linear chain version of $\mathbf{S}$ that we name $\mathbf{L}$ (line 2,
Algorithm 2). From a graph-theoretic point of view, this operation removes
(i.e., sets to zero) the edges adjacent to branching vertices. The indexing of
the matrix does not change, but its nonzeros do. For example, if row 10 is a
branching vertex, the entire row—and column, since $\mathbf{S}$ is
symmetric—is cleared, but row 10 is still a row in the matrix. This way we can
avoid re-indexing sequences during the computation. The result is an even
sparser string graph with nodes of degree $0$, $1$ or $2$. For each contig,
there are exactly two vertices of degree $1$, which we call root nodes (or
root vertices) and use as starting points for depth-first search and local
assembly of the contig in the final stage of contig generation.

Once the string matrix is in its unbranched form $\mathbf{L}$, we want to
decompose it into its linear components to produce manageable independent
subproblems that we can work on in parallel (i.e., local assemblies).
Therefore, we use the sparse matrix based connected components (LACC)
algorithm presented by Azad et al. (Azad and Buluç, 2019) to determine the
contig set. LACC is a distributed-memory implementation of the Awerbuch-
Shiloach algorithm (Awerbuch and Shiloach, 1987) using the CombBLAS library.
LACC takes advantage of the sparsity of the vector to avoid processing
inactive vertices. It takes as input the unbranched sparse matrix $\mathbf{L}$
and returns a distributed vector $\mathbf{v}$ (line 3, Algorithm 2), which is
a mapping from global sequence indices (i.e., rows and columns of
$\mathbf{S}$) to contig indices $C_{i}$. In the previous graph example, LACC
would return ${v_{1},v_{2}}\in C_{1}$, ${v_{4},v_{5},v_{6}}\in C_{2}$, and
${v_{7},v_{8}}\in C_{3}$.

Once we have determined which sequences belong to which contig, we estimate
the size of each contig by counting how many sequences belong to it. However,
since $\mathbf{v}$ is a distributed vector, a processor may not know the full
set of sequences belonging to the contig it owns. Therefore, we need to
communicate this information across the processor grid. Each processor
computes a local size estimate based on the vertices it owns locally for each
local contig. An MPI Reduce-scatter collective operation is used to determine
the global size of each contig and redistribute the sizes across the processor
grid to create a distributed mapping of contig indices to their associated
global sizes.

###  4.3. Contig Load Balancing and Communication

Previously, we determined which sequences belong to which contig and estimated
the size of each contig based on the number of sequences associated with it.
Using the contig size as an optimization parameter, we want to distribute the
workload as evenly as possible across the processor grid.

#### Load Balancing Algorithm

Given a vector of contig sizes of length $n$ and $P$ processes, we want to
find the near-optimal contig-to-processor assignment such that the amount of
work estimated by the contig size that each processor has to do is similar.
Our objective and problem definition is similar to the _multiway number
partitioning_ optimization problem first introduced by Ronald Graham in the
context of the identical-machines scheduling problem (Graham, 1966). The
classical application is to schedule a set of $m$ jobs with different runtimes
on $k$ identical machines in such a way that the makespan, i.e., the elapsed
time until the schedule is completed, is minimized.

In the context of contig generation, this means that given a multiset of $S$
instances (in our case, the contig sizes), we want to partition this multiset
into $P$ subsets (the number of processes in our processor grid and must be a
positive integer) such that the sums of the subsets (i.e., the sums of the
contig sizes) are as similar as possible. The partitioning results in a
balanced distribution of the workload for the next and final phase of the
computation.

The multiway number partitioning problem is NP-hard. To overcome this
limitation, in ELBA we use an approximation algorithm known in the scheduling
literature as the Longest Processing Time (LPT) algorithm, whose goal is to
minimize the largest subset, which belongs to a class of algorithms known as
_greedy number partitioning_. The algorithm loops over the contig sizes and
inserts each number into the set whose current sum of sizes is smallest. The
result is a partitioning that minimizes the time processes spend waiting for
the most heavily loaded process to finish assembling its contig subset. If the
contig sizes were not sorted, then the runtime would be $O(n)$ and the
approximation ratio would be at most $2-1/P$. It is possible to improve the
approximation ratio to $(4P-1)/3P$ by sorting the input vector of contig sizes
(Graham, 1969).

The improvement in the approximation using LPT is accompanied by an increase
in the runtime to $O(n\textrm{log}n)$ due to sorting. However, the number of
contigs $n$ is smaller than the number of sequences by at least an order of
magnitude, and the increase in runtime does not create a computational
bottleneck. For the same reason, we collect the global information about
contig lengths in a single processor and run the partitioning algorithm on it
to avoid the unnecessary communication of small messages. The partitioning
algorithm returns the vector $\mathbf{p}$ specifying the assignment of contigs
to processes (line 4, Algorithm 2)—$\mathbf{p}$ is broadcasted to the entire
processor grid so that each local process can determine where to send its
local sequences and associated information.

As mentioned earlier, the problem size at this stage of the computation is
often smaller than the problem size at the initial stage (i.e., overlap
detection), so for some species it is possible that $n<P$. In this case, some
of the processes are idle for the final phase of the computation. For the two
species we use in the experimental evaluation in Section 6.2, $n$ is equal to
$6411$ and $4287$ and $P$ varies from $18$ to $128$. In the next section, we
explain how contigs are redistributed among processes based on their size.

#### Induced Subgraph Algorithm

Once we have determined where a contig and its sequences should be stored, we
must perform the actual communication step to send the linear component
information and associated sequences to the owner processor.

Communicating both the linear component information and the sequences involves
the same high-level procedure of reassigning vertices representing a linear
chain to their owner processor. However, the underlying data structures are
fundamentally different, namely that the overlap graph is stored as a sparse
matrix while the sequences are stored as distributed char arrays, and
therefore require a different implementation.

Let us first focus on the communication of linear component information, i.e.
the graph-like structure that stores connectivity information. We have as
input the sequence-by-sequence matrix $\mathbf{L}$, the resulting matrix after
the transitive reduction step and from which we have cut out vertices of
degree $\geq 3$. $\mathbf{L}$ is distributed over $P$ processes, which are
logically organized in a $\sqrt{P}\times\sqrt{P}$ grid. Let $n$ be the number
of vertices in $\mathbf{L}$, where $\mathbf{L}=(V,E)$ in the graph
interpretation. From the load balancing algorithm, we have a distributed
vector $\mathbf{v}:[n]\to[P]$ such that $\mathbf{v}[u]=P_{i}$ means that
vertex $u$ should belong to processor $P_{i}$.

The goal is to create an induced subgraph—or submatrix—$\mathbf{L}^{(P_{i})}$
locally on each processor $P_{i}$, i.e., a graph formed by a subset of the
vertices of the original graph $\mathbf{L}$ and the edges connecting the
vertices in this subset. Formally, we define an induced subgraph
$\mathbf{L}^{(P_{i})}$ such that
$\mathbf{L}^{(P_{i})}=(V^{(P_{i})},E^{(P_{i})})$ with $V^{(P_{i})}=\\{v\in
V:P_{i}=\mathbf{v}[v]\\}$ and $E^{(P_{i})}=\\{(u,v)\in E:u,v\in
V^{(P_{i})}\\}.$

The vector $\mathbf{v}$ is also distributed across the
$\sqrt{P}\times\sqrt{P}$ processor grid and is therefore divided into $P$
subvectors, each of size $\approx n/P$; note that we write
$\mathbf{v}_{(i,j)}$ to denote the subvector on process $P(i,j)$. Each process
is only aware of which vertices it stores to send to other processes.
Therefore, we need a communication step to make each process aware of which
vertices to receive. Given the way $\mathbf{v}$ is distributed, we can avoid
an MPI_Allgather operation spanning the entire grid and instead use the square
process grid to communicate $\mathbf{v}$ in a scalable way. That is, we
perform an allgather operation over the $\mathbf{Row}$ dimension followed by
point-to-point communication to access the information stored on the
$\mathbf{Column}$ dimension.

Figure 2. An example of how the induced subgraph algorithm communicates
vertices across the processor grid.

More precisely, let $(u,v)\in\mathbf{L}_{i,j}$ be any nonzero value, where
$\mathbf{L}_{i,j}$ is the local submatrix stored on processor $P(i,j)$. The
goal is to determine $\mathbf{v}[u]$ and $\mathbf{v}[v]$. Given we know
$(u,v)$ is stored in the submatrix $\mathbf{L}_{i,j}$, we know that
$\mathbf{v}[u]$ is stored in the row $P(i,:)$. Thus, we can access
$\mathbf{v}[u]$ by computing an MPI_Allgather operation over the processor row
$P(i,:)$ and, more generally, over the $\mathbf{Row}$ dimension of
$\mathbf{S}$. To determine $\mathbf{v}[v]$, we first note that if we are a
diagonal processor, i.e. $i=j$, we already have access to $\mathbf{v}[v]$ via
the previous MPI_Allgather operation. If $i\neq j$, then we know that
$\mathbf{v}[v]$ is stored on the transposed processor $P(j,i)$ because we
performed the allgather operation over the $\mathbf{Row}$ dimension, i.e.
$P(j,:)$. Therefore, we perform a point-to-point communication between
$P(i,j)$ and $P(j,i)$ to exchange the subvector $\mathbf{v}_{i,*}$ and the
subvector $\mathbf{v}_{j,*}$. This gives the processor $P(i,j)$ access to
$(u,v)$. It follows that each processor $P(i,j)$ now has access to every entry
in $\mathbf{v}$ corresponding to a nonzero in $\mathbf{L}_{i,j}$. A similar
procedure is previously implemented for distributed-memory breadth-first
search (Buluç and Madduri, 2011).

Figure 2 illustrates an example with 1-based indexing, where $u=1$ (red square
in the leftmost matrix) is stored on process $P(1,1)$ and $v=12$ (green
square) is stored on process $P(3,4)$. If we assume that $P(1,3)$ requires
access to $\mathbf{v}[1]$ and $\mathbf{v}[12]$, an MPI_Allgather operation on
the $P(1,:)$ row would provide access to $\mathbf{v}[1]$, while the same
operation on the $P(3,:)$ row would provide access to $\mathbf{v}[12]$ on the
transposed processor $P(3,1)$. Point-to-point communication between $P(1,3)$
and $P(3,1)$ would then allow $P(1,3)$ to access $\mathbf{v}[12]$.

Once we have access to the processor assignment stored in $\mathbf{v}$, we
communicate edges to their target processor and build the local induced
subgraph. To do this, we loop through each nonzero $(u,v)\in\mathbf{L}_{i,j}$
stored on processor $P(i,j)$, and for each $(u,v)$ where $\mathbf{v}[u]$ and
$\mathbf{v}[v]$ are destined for the same destination processor, we construct
a triple $(u,v,\mathbf{S}(u,v))$ and place it on an outgoing buffer to the
destination processor. A custom all-to-all ensures the required re-
distribution of non-zeros. Once each processor has access to its edge set
$E^{(P_{i})}$, the local adjacency matrix of the induced subgraph is
constructed; while we re-index the local matrix to fit its new, smaller size,
we also keep a map of the original global vertex indices, since it is needed
in the final phase.

#### Read Sequence Communication

The communication of read sequences is implemented separately, since the read
sequences are not stored as nonzeros in the sparse matrix $\mathbf{L}$, but in
a distributed auxiliary data structure. Furthermore, a sequence is represented
and stored as a char array. A large data set could exceed the MPI limit of
$2^{31}-1$. Therefore, we need to treat sequence communication separately and
consider this potential limitation.

Read sequences that need to be sent to a processor other than the one they are
currently on are packed into a char buffer and communicated as a sequence of
non-blocking point-to-point messages in an all-to-all fashion. To deal with
the MPI $2^{31}-1$ count limit, we check the length of each message to be
communicated and its receive buffer. If it goes beyond the limit, we
communicate the sequences using a user-defined contiguous MPI data type whose
size is equal to the buffer length. This way we can send and receive each
character buffer in a single call.

###  4.4. Local Contig Assembly

Let $\mathbf{L}_{i,j}$ be the local graph (or matrix), composed of one or more
linear components, stored on the processor $P_{i,j}$ obtained via the induced
subgraph algorithm. There are $P$ such graphs, one for each processor, but we
can assume without loss of generality that we are dealing only with
$\mathbf{L}(P_{i,j})$ or $\mathbf{L}$ for short.

Suppose $\mathbf{L}$ has $n$ vertices and $m$ edges. Each vertex of
$\mathbf{L}$ is a read sequence $l$ assigned by the partitioning algorithm to
$P_{i,j}$. Then, we define $\Sigma=\\{A,C,T,G\\}$ as the alphabet of DNA
nucleotides. Consequently, we represent the $n$ sequences of $\mathbf{L}$ by
$l_{0},l_{1},\ldots,l_{n-1}\in\Sigma^{*}$, which means that each sequence is a
combination of $\\{A,C,T,G\\}$. Given any read sequence $l$, we express the
nucleotides of that sequence by $(l[0],l[1],\ldots,l[|l|-1])$, where $|l|$
denotes the length of the sequence. If $l[i]$ is a base or nucleotide in
$\Sigma$, we denote the Watson-Crick complement base of $l[i]$ by $l[i]^{c}$.
This allows us to generalize the notion of sequence to include its reverse
complement, as defined in Section 2: if $i<j$, then we write $l[i:j]$ for the
substring $(l[i],l[i+1],\ldots l[j])$ and $l[j:i]$ for its reverse complement
substring $(l[j]^{c},[j-1]^{c},\ldots,l[i]^{c})$.

$\mathbf{L}$ can also be viewed from the point of view of its matrix
representation, such that $\mathbf{L}$ is a sparse $(n\times n)$ matrix with
$m$ nonzeros. In the earlier stages of the pipeline, we use the doubly
compressed sparse column (DCSC) format (Buluç and Gilbert, 2008) to store our
matrices for scalability. However, for the local traversal algorithm described
in this section, we converted these matrices to compressed sparse column (CSC)
format for simplicity and faster vertex (column) indexing. The storage space
required to store the local matrices is an order of magnitude smaller than
before and hence CSC does not introduce memory scalability issues despite
hypersparsity. This conversion takes linear time in the number of local
vertices, as only column pointers needs to be uncompressed and row indices
array stays intact, and has negligible effect to overall runtime.

As a sparse matrix, $\mathbf{L}(u,v)=e$ is a nonzero provided that the source
sequence $l_{u}$ and the destination sequence $l_{v}$ overlap. In this nonzero
$e$ we store two values computed from the original string graph $\mathbf{S}$,
which is all we need for assembly: $\mathbf{pre}_{u}(e)$ and
$\mathbf{post}_{v}(e)$. To define $\mathbf{pre}_{u}(e)$ and
$\mathbf{post}_{v}(e)$, we use inclusive indexing, resulting in the following
asymmetric definition. The first value $\mathbf{pre}_{u}(e)$ stores the index
$i$ of $l_{u}$, which is the last nucleotide in $l_{u}$ that does not overlap
with $l_{v}$, i.e., the nucleotide on $l_{u}$ that precedes the overlap with
$l_{v}$. The second value $\mathbf{post}_{v}(e)$ stores the index $j$ of
$l_{v}$, which is the first nucleotide in $l_{v}$ that overlaps with $l_{u}$,
i.e., the index of the beginning of the overlap between $l_{u}$ and $l_{v}$.
The definition is asymmetric because $l_{u}[i]$ and $l_{v}[j]$ do not store
the same nucleotide. This asymmetry is necessary to compute the non-
overlapping prefixes of each read in a contig, which, when joined together,
produce a contig sequence.

Figure 3. An example of how the local assembly algorithm concatenates
sequences.

For example, if $l_{u}=l_{0}=AGAACT$ and $l_{v}=l_{1}=AACTGAAG$ as shown in
Figure 3 and we consider 0-based indexing, then $\mathbf{pre}(e)=1$ and
$\mathbf{post}(e)=0$ because $l_{u}[2:5]$ and $l_{v}[0:3]$ overlap. If instead
we consider $l_{u}=l_{1}=AACTGAAG$ and $l_{v}=l_{2}=TGAAGAA$, then
$\mathbf{pre}(e)=4$ and $\mathbf{post}(e)=2$ because $l_{u}[5:7]$ and
$l_{v}[2:4]$ overlap. It is also possible that $l_{u}=l_{0}$ would overlap
with $l_{v}^{c}=l_{1}^{c}=CTTCAGTT$ (the reverse complement of $l_{1}$). In
this case, the same procedure would apply to compute $\mathbf{pre}(e)=1$ and
$\mathbf{post}(e)=4$, where the overlapping subsequences would be $l_{u}[2:5]$
and $l_{v}^{c}[7:4]$, since the $l_{u}$ subsequence $AACT$ would overlap with
the $l_{v}^{c}$ reserve complement subsequence $AGTT$.

Theoretically, the beginning of the overlapping substring of $l_{v}$ should
always be $0$ or $|l_{v}|-1$, depending on the orientation. Yet, we still
store $\mathbf{post}(e)$ because the diBELLA 2D (Guidi et al., 2021c)
pipeline, which we extend in this work, uses an x-drop seed-and-extend
approach to align overlapping sequences. With x-drop, the alignment between
two sequences can potentially end early (i.e., before extending to the end of
a read), leaving a short overhang in the alignment coordinates at the end of
the sequence. If $l_{u}=l_{1}=AACTGAAG$ and $l_{v}=l_{2}=TGAAGAA$ are defined
as in Figure 3, and the alignment told us that $l_{u}[5:7]$ and $l_{v}[2:4]$
were the overlapping region, it would be incorrect to match
$\mathbf{pre}(e)=4$ with $\mathbf{post}(e)=0$, since we need
$\mathbf{post}(e)=2$ to correctly concatenate the subsequences that form a
contig. If we join the three partial sequences marked in red on
$l_{0},l_{1},l_{2}$ in Figure 3, as explained below, we get a contig.

Table 1. Details of the machines used for evaluation; the frequency reported
is the processor max turbo frequency.

Platform | Cores/Node | Frequency (GHz) | Processor | Memory (GB) | Network and Topology | L1 | L2 | L3
---|---|---|---|---|---|---|---|---
Cori Haswell | 32 | 3.6 | Intel Xeon E5-2698V3 | 128 | Aries Dragonfly | 64KB | 256KB | 40MB
Summit CPU | 42 | 4.0 | IBM POWER9 | 512 | InfiniBand Non-Blocking Fat Tree | 32KB | 512KB | 10MB

Table 2.  Data sets used during evaluation: name, depth, number of sequences
in the input, average read length, input size, genome size, and error rate.

Label | Depth | Reads (K) | Length | Input (GB) | Size (Mb) |  Error (%)
---|---|---|---|---|---|---
O. sativa |  30 | 638.2 | 19,695 | 12.2 | 500 |  0.5
C. elegans |  40 | 420.7 | 14,550 |  3.8 | 100 |  0.5
H. sapiens |  10 | 4,421.6 | 7,401 | 31.1 | 3,200 | 15.0

Because the local assembly operates directly on the CSC format, we describe it
briefly. $\mathbf{L}$.JC is the column pointer array of length $n+1$,
$\mathbf{L}.IR$ is the row index array of length $m$, and
$\mathbf{L}.\textit{VAL}$ is the array of tuples
$(\mathbf{pre}(e),\mathbf{post}(e))$, also of length $m$.

The local assembly algorithm is a variant of depth-first search, simplified by
the fact that the maximum vertex degree of $\mathbf{L}$ is $2$ by
construction. For this reason, there is always only one vertex in the
frontier, and the search is thus a linear walk. Each local matrix
$\mathbf{L}_{i,j}$ is assigned a contig set by multiway number partitioning,
which is the set of connected components of $\mathbf{L}_{i,j}$. Each vertex of
$\mathbf{L}$ has degree $1$ or $2$. Moreover, we define contigs as linear
chains of at least two sequences, where $q$ is the number of vertices in a
given connected component. It follows that any connected component consisting
of $q\geq 2$ vertices must have exactly $2$ vertices of degree $1$ (i.e.,
_root_ vertices) and $q-2$ vertices of degree $2$ (i.e., _intermediate_
vertices). In short, the idea is to scan for root vertices in $\mathbf{L}$
and, if we find them, to take a walk from that root vertex (via as many
intermediate vertices as necessary) until another root vertex is found. As we
proceed, we perform the concatenation of subsequences described below. This
search for the root vertex is performed over all $n$ vertices. Therefore, we
must mark each final root vertex (found by the linear traversal) as visited to
avoid accidentally composing the same contig twice.

More precisely, we loop over the $n$ sequence vertices in each subgraph
$\mathbf{L}_{i,j}$ and compute
$\mathbf{L}_{i,j}.JC[i+1]-\mathbf{L}_{i,j}.JC[i]$ from the column pointer
array representing the degree of vertex $i$. Each time we find a vertex of
degree $1$, i.e., a root vertex, that has not been visited before, we perform
a traversal starting from this that root vertex $r$. Given a generic vertex in
the chain $c$ (which can be either a root vertex or an intermediate vertex),
we find the successor vertex $c^{\prime}$ by examining the vertices in the
slice $\mathbf{L}_{i,j}.IR[\mathbf{L}_{i,j}.JC[c]:\mathbf{L}_{i,j}.JC[c+1]]$
of the row index array where the edges are stored. Each vertex $c$ has at most
two successor vertices, and we select the unvisited one. This process
continues until we reach the vertex $r^{\prime}$ for which
$\mathbf{L}_{i,j}.JC[r^{\prime}+1]-\mathbf{L}_{i,j}.JC[r^{\prime}]=1$, i.e.,
the second root vertex of a contig. The result is a chain of $q$ vertices
$r,c_{1},\ldots,c_{q-2},r^{\prime}$ for each contig stored in the submatrix
$\mathbf{L}_{i,j}$. As we proceed, we collect the edges between them
$e_{0},e_{1},\ldots,e_{q-2}$. Because the values for $\mathbf{pre}(e_{i})$ and
$\mathbf{post}(e_{i})$ have already been computed and stored in each edge
$e_{i}$, we know exactly which subsequences of which read sequence to look up
and join to form the contig. Namely, $l_{r}[\alpha:\mathbf{pre}(e_{0})]\oplus
l_{c_{1}}[\mathbf{post}(e_{0}):\mathbf{pre}(e_{1})]\oplus\ldots\oplus
l_{c_{q-2}}[\mathbf{post}(e_{q-3}):\mathbf{pre}(e_{q-2})]\oplus
l_{c_{r^{\prime}}}[\mathbf{post}(e_{q-2}):\beta],$ where $\alpha=0$ or
$|l_{r}|-1$ depending on the orientation of $l_{r}$, and $\beta$ is similarly
defined for $l_{r^{\prime}}$.

This algorithm is performed by each process in parallel on its own induced
subgraph $\mathbf{L}_{i,j}$. For each read in the contig, we either look for
the subsequence in the locally stored char array that we started with, or in
the char array obtained by communicating the read sequence. Because we store
read sequences in large packed char arrays, we do not need to copy the entire
read sequence to find the correct subsequence. Instead, we can simply use the
offsets already computed, which tell us where each read is in the buffer, and
then read the subsequence directly from the buffer. The algorithm is $O(q)$,
where $q$ is the number of vertices in a connected component as previously
defined (i.e., the number of sequences in a contig), since the search for the
root vertex takes linear time and the traversal is linear in the number of
edges, which is $2(q-1)$.

##  5\. Experimental Setup

To evaluate our contig generation algorithms and the ELBA long-read assembly
pipeline we used two machines: the Haswell partition of the Cray XC40
supercomputer Cori at NERSC and the IBM supercomputer Summit at Oak Ridge
National Laboratory, on which we used only IBM POWER9 CPUs (Table 1).

Each node on the Cori Haswell partition is a dual-socket 2.3GHz Intel Xeon
with 16 cores each and 128GB of total memory, while each Summit node is a
dual-socket IBM POWER9 with 22 cores each and 512GB of DDR4 from RAM, but only
42 cores are available because 2 of them are reserved for the operating
system. On Cori Haswell, we used gcc-11.2 and the O3 flag to compile C/C++
codes, while on Summit we used gcc-9.2. On both machines, we used the default
MPI implementation.

To evaluate our algorithm, we use three different species: Oryza sativa (O.
sativa), Caenorhabditis elegans (C. elegans), and Homo sapiens (H. sapiens)
whose characteristics are summarized in Table 2. The evaluation is divided
into two categories: (a) runtime and scalability of the ELBA pipeline
including our novel contig generation algorithm on two machines, and (b)
runtime compared to state-of-the-art shared memory software. For low error
rate sequences (O. sativa and C. elegans), we also report assembly quality.

The state-of-the-art software we consider are Hifiasm (Cheng et al., 2021) and
HiCanu (Nurk et al., 2020), because of their speed and popularity,
respectively, and both are written for shared-memory parallelism. HiCanu can
optionally run on grid computing, but is not implemented for distributed
memory parallelism. It is worth noting that Hifiasm and HiCanu include
additional polishing stages that make their overall assembly quality higher
than ELBA’s. For the H. sapiens dataset, we consider Miniasm (Li, 2016) and
Canu (Koren et al., 2017), as this dataset has a much higher error rate that
is not suitable for Hifiasm and HiCanu. Our goal is to show the
competitiveness and potential of ELBA in terms of assembly quality,
demonstrating in particular clear advantages in terms of runtime. ELBA was run
with the k-mer length parameter $k=31$ and the x-drop threshold $x=15$ for the
low error rate data and with $k=17$ and $x=7$ for H. sapiens in Table 2.
Hifiasm and HiCanu were run with their default setting.

To show performance, we run ELBA on both Cori Haswell and Summit (except for
H. sapiens, where we only use Summit), while for comparison with the state of
the art, we only use Cori Haswell, since Hifiasm and HiCanu use SSE and AVX2
intrinsics, which are not supported on the IBM POWER9 processor on Summit.
Hifiasm and HiCanu were developed for shared memory, so we only give runtimes
for a single Cori Haswell node using multi-threading. The lack of support for
AVX2 intrinsics is also the reason why ELBA’s alignment is slower on Summit
than on Cori. For this reason, we want to emphasize that the goal of using two
machines is to show performance and scaling on different systems, not to
directly compare the two machines, since ELBA is optimized for a general HPC
system, i.e., our code is general and no architecture-specific optimizations
have been made.

To assess the quality of the contig set of ELBA, Hifiasm, and HiCanu, we use
QUAST (Gurevich et al., 2013) and report the following metrics: Completeness,
longest contig size, number of contigs, and misassembled contigs. Completeness
measures the percentage of the reference genome to which at least one contig
has been aligned. This is calculated by counting the number of nucleotides
aligned to the reference genome and dividing by the total length. The number
of misassembled contigs is defined as the number of contigs that contain
incorrect assemblies, e.g., a contig consisting of sequences originating from
different regions of the reference genome.

##  6\. Results

In this section, we evaluate the performance and quality of the overall ELBA
pipeline and our novel contig generation algorithm, both individually and
compared to the state of the art.

###  6.1. Performance and Scalability

Figure 4. ELBA strong scaling on Cori Haswell and Summit CPU using 32 MPI
rank/node on C. elegans (left) and on O. sativa (right). Figure 5. ELBA
runtime breakdown of the main stages of the pipeline on Cori Haswell and
Summit for C. elegans on the left and for O. sativa on the right.

Figure 4 illustrates the strong scaling of the entire ELBA pipeline for C.
elegans on the left and for O. sativa on the right. The C. elegans dataset was
run on $P=$ {18, 32, 50, 72, 128} nodes using 32 MPI ranks/nodes on both
machines, while the O. sativa was run on $P=$ {18, 32, 50, 72, 128} nodes
using 32 MPI ranks/nodes on Summit CPU and on $P=$ {50, 72, 128} nodes Cori
Haswell because the algorithm for $P=$ {18, 32} ran out of memory since Cori
Haswell has a smaller memory per node than Summit.

ELBA achieves a parallel efficiency of 75% on Cori Haswell and 69% Summit CPU
for C. elegans, while it achieves a parallel efficiency of 80% and 64% for O.
sativa. The parallel efficiency of O. sativa between 72 and 128 nodes on
Summit is 83%, which is similar to the parallel efficiency on Cori between 50
and 128 nodes. For the H. sapiens dataset, the parallel efficiency on Summit
between 200 and 392 nodes is close to 90%, as shown in Figure 6 on the left.
These results show a good scaling behavior of the whole ELBA pipeline with a
large input on two different architectures.

Figure 5 shows the runtime breakdown of ELBA on the two machines for C.
elegans on the left and for O. sativa on the right, while Figure 6 on the
right shows the runtime breakdown of ELBA on Summit for H. sapiens. To
highlight the impact of the major stages, we omit I/O and other minor
computation from the breakdown. The overall impact of the omitted computation
is negligible.

The legend is arranged in reverse order with respect to the layers of the bar,
i.e. the first entry of the legend at the top is associated with the first
layer in the stacked bar from the bottom. CountKmer corresponds to the k-mer
counting step, then follows the step DetectOverlap, which represents the time
to create and compute the candidate overlap matrix $\mathbf{C}$, i.e., the
first step shown in Figure 1. Then comes the Alignment step, i.e., the time
required to perform pairwise alignment on each nonzero in the candidate
overlap matrix $\mathbf{C}$, followed by TrReduction, i.e., the transitive
reduction time and the second step shown in Figure 1. Finally, ExtractContig
is the time we spend extracting the contig set from the sparse matrix
$\mathbf{S}$ and is shown as the third step in Figure 1. The ExtractContig
step is the core contribution of this work, but its implementation is an
essential part of making the entire ELBA pipeline work and for this reason it
is key to showing the scalability of the pipeline.

Figure 6. ELBA strong scaling (on the left) and runtime breakdown (on the
right) of the main stages of the pipeline on Summit for H. sapiens.

Figure 5 shows the breakdown of ELBA performance for the main phases of
computation for C. elegans and O. sativa using $P{=}\\{50,72,128\\}$ nodes on
each machine while Figure 6 shows the breakdown of ELBA performance breakdown
for H. sapiens using $P{=}\\{200,288,338,392\\}$ Summit nodes. ELBA is faster
overall on Cori Haswell than on Summit CPU. The CountKmer, DetectOverlap, and
Alignment phases show nearly linear scaling on both machines, with the
exception of CountKmer for H. sapiens on Summit, which scales only
sublinearly. The relative contribution of pairwise alignment to the overall
runtime increases on Summit CPU compared to Cori Haswell and largely
contributes to the discrepancy in runtime because the alignment library used
in ELBA is not optimized for the IBM processor. The TrReduction and
ExtractContig phases are also significantly faster on Cori Haswell than on
Summit CPU. On Summit, the computation for these two phases does not scale and
consumes a higher percentage of runtime. This is because the amount of work is
smaller in these two phases and the algorithms are latency-bound. Further,
Summit’s network has lower performance compared to Cori Haswell’s. Summit CPU
has lower network bandwidth per core, and we also only use 32 of the 42
available cores on Summit to make the comparison with Cori fair, which does
not saturate Summit’s bandwidth.

In each dataset, 65-85% of the runtime of contig generation on both machines
is taken by the induced subgraph function described in Section 2, which mainly
involves communication. The ExtractContig never requires more than 5% of the
computation for each species and each machine, demonstrating the efficiency of
our contig generation algorithm.

ELBA’s contig generation focuses on localizing the graph traversal problem so
that the sequences that make up a contig are organized into localized matrices
on each processor. The working set at this stage is smaller than at the
beginning of the computation, so the result is a fast and efficient, but
latency-bound algorithm.

###  6.2. Comparison with the State-of-the-Art

To demonstrate the competitiveness and advantages of ELBA over state-of-the-
art long read assembly software, we compare the runtime and scaling of ELBA
with Hifiasm (Cheng et al., 2021) and HiCanu (Nurk et al., 2020) for O. sativa
and C. elegans and with Mifiasm (Li, 2016) and Canu (Koren et al., 2017) for
H. sapiens. The competing software are designed for shared memory parallelism.
Therefore, we ran them on a single Cori Haswell node with multithreading,
while we ran ELBA on $P=\\{18,50,128\\}$ with 32 MPI ranks/nodes. Due to the
hard-to-separate differences between these software, we decided to make a
comparison based on total runtime. Hifiasm, HiCanu, Mifiasm, Canu, and ELBA
perform very different computation, but they ultimately aim to solve the same
problem. For the sake of completeness, we also compare the quality of
assembly. It is worth noting that both Hifiasm and HiCanu perform a polishing
phase. This can lead to a slight disadvantage in runtime, but ensures a better
assembly quality.

Table 3 summarizes the runtime performance of Hifiasm and HiCanu, in the
rightmost column the speedup of ELBA over those software for $P=\\{18,128\\}$
and $P=\\{50,128\\}$ for C. elegans and O. sativa, respectively. ELBA is up to
$15\times$ faster than Hifiasm and up to $58\times$ faster than HiCanu for the
C. elegans dataset, while the speedup for O. sativa, which is a larger genome,
reaches $36\times$ over Hifiasm and up to $159\times$ over HiCanu. On Cori,
runtimes for Hifiasm and HiCanu for O. sativa are 17 and 64 minutes,
respectively, while ELBA takes less than 2 minutes on 128 nodes. For C.
elegans, runtimes of Hifiasm and HiCanu are approximately 1 hour and 5 hours,
while ELBA requires 1 minute on 128 nodes. For H. sapiens, ELBA on 392 nodes
on Summit takes 17 minutes, while Mifiasm on Cori Haswell takes 1 hour and 30
minutes with 32 threads and Canu ran out of time after more than 64 hours of
runtime with 32 threads on Cori Haswell. However, these times are from
different architectures and are not directly comparable, so we do not give a
speedup for H. sapiens. The results from O. sativa and C. elegans show that
performance on Summit was worse than on Cori Haswell due to the less powerful
processor on Summit.

Finally, we compare the quality of the assemblies of ELBA with that of Hifiasm
and HiCanu. Table 4 summarizes metrics obtained from QUAST (Gurevich et al.,
2013). Specifically, we show the completeness, i.e., the proportion of the
reference genome that was covered by at least one contig (the higher, the
better), then the length of the longest contig (the higher, the better), the
size of the contig set (the lower, the better if associated with high
completeness), and the number of misassemblies (the lower, the better). In
both species, ELBA has competitive genome completeness, especially in C.
elegans, where it is higher than both Hifiasm and HiCanu, and misassemblies.
In ELBA, the contigs are significantly shorter than in the two competing
software. This is understandable, since ELBA does not currently perform a
polishing step, which is reserved as future work.

Table 3. ELBA’s speedup over state-of-the-art software. Hifiasm and HiCanu are
designed for shared memory parallelism and were run on a single Cori node with
multithreading.

Tool | Organism | Runtime (s) | Nodes (32 MPI Rank/Node) | ELBA Speed-Up
---|---|---|---|---
Hifiasm | C. elegans |  1,015.2 | 18–128 |  3–15$\times$
HiCanu | C. elegans |  3,819.0 | 18–128 |  11–58$\times$
Hifiasm | O. sativa |  4,131.9 | 50–128 |  18–36$\times$
HiCanu | O. sativa | 18,131.0 | 50–128 | 78–159$\times$

Table 4. Comparison of assembler quality for O. sativa (top) and C. elegans
(bottom). Hifiasm and HiCanu implement additional polishing stages to improve
their metrics.

Tool | Completeness (%) | Longest Contig (Mb) | Contigs | Misassembled Contigs
---|---|---|---|---
ELBA | 37.09 |  0.172 | 6411 | 2
Hifiasm | 26.94 |  7.083 | 1661 | 1
HiCanu | 25.94 | 37.523 |  168 | 2
ELBA | 98.93 |  0.313 | 4287 | 5
Hifiasm | 99.96 |  6.438 |  133 | 0
HiCanu | 99.90 | 18.332 |  32 | 2

##  7\. Conclusions

Recent advances in sequencing technologies have increased the need for high-
performance approaches, as we can now generate more data at lower cost, which
in turn requires higher computational resources to reconstruct high-quality
assemblies in a timely manner. In this paper, we presented the contig
generation phase of the distributed-memory long-read assembler ELBA.

Contig generation is critical to assembly functionality, as it enables the
construction of longer genomic sequences that represent a physical map of a
region of a chromosome. ELBA has achieved speed increases of up to $36\times$
and $159\times$ compared to two state-of-the-art software for the O. sativa
dataset, opening the door for high-performance genome assembly.

ELBA and its contig generation step rely on distributed sparse matrices to
first determine the contig set starting from a string graph. Then, using a
greedy multiway number partitioning algorithm, it determines how the rows and
columns of the sparse matrix representing DNA sequences are redistributed
between processes so that sequences belonging to the same contig are stored
locally on the same processor. ELBA then uses such partitioning and the sparse
matrix abstraction to implement the induced subgraph function and redistribute
the sequences among the processes. Finally, ELBA’s contig generation step
computes a local assembly step, i.e., the actual concatenation of sequences
into a contig, on each processor independently and in parallel.

Future work includes developing a polishing or scaffolding phase to further
improve the quality of ELBA assembly. One possibility is to once again use the
sparse matrix abstraction to find similarities within the contig set and
obtain even longer sequences. In addition, we plan to reduce the memory
consumption of ELBA so that we can assemble large genomes at low concurrency,
and optimize ELBA for running in a cloud environment as high-performance
scientific computing in the cloud becomes more popular (Guidi et al., 2021a;
Reed et al., 2022; Mesnard and Barba, 2019). Finally, we plan to use GPUs in
several stages, such as alignment (Zeni et al., 2020; Awan et al., 2020) and
k-mer counting (Erbert et al., 2017; Nisa et al., 2021), to improve our
performance and take full advantage of today’s heterogeneous systems.

###### Acknowledgements.

This work is supported by the Advanced Scientific Computing Research program
within the Office of Science of the DOE under contract number DE-
AC02-05CH11231, and by the Exascale Computing Project (17-SC-20-SC), a
collaborative effort of the DOE Office of Science and the NNSA. We used
resources of the NERSC supported by the Office of Science of the DOE under
Contract No. DEAC02-05CH11231 and of the OLCF supported by the Office of
Science of the DOE under Contract No. DE-AC05-00OR22725.

## References

  * (1)
  * Awan et al. (2020) Muaaz G Awan, Jack Deslippe, Aydın Buluç, Oguz Selvitopi, Steven Hofmeyr, Leonid Oliker, and Katherine Yelick. 2020\.  ADEPT: a domain independent sequence alignment strategy for GPU architectures.  _BMC Bioinformatics_ 21, 1 (2020), 1–29.
  * Awerbuch and Shiloach (1987) Baruch Awerbuch and Yossi Shiloach. 1987.  New connectivity and MSF algorithms for shuffle-exchange network and PRAM.  _IEEE Trans. Comput._ 36, 10 (1987), 1258–1263.
  * Azad and Buluç (2019) Ariful Azad and Aydın Buluç. 2019.  LACC: A Linear-Algebraic Algorithm for Finding Connected Components in Distributed Memory. In _IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. ACM New York, NY, USA.
  * Azad et al. (2022) Ariful Azad, Oguz Selvitopi, Md Taufique Hussain, John R. Gilbert, and Aydın Buluç. 2022\.  Combinatorial BLAS 2.0: Scaling Combinatorial Algorithms on Distributed-Memory Systems.  _IEEE Transactions on Parallel & Distributed Systems_ 33, 4 (2022), 989–1001.  https://doi.org/10.1109/TPDS.2021.3094091
  * Bachan et al. (2019) John Bachan, Scott B. Baden, Steven Hofmeyr, Mathias Jacquelin, Amir Kamil, Dan Bonachea, Paul H. Hargrove, and Hadia Ahmed. 2019\.  UPC++: A High-Performance Communication Framework for Asynchronous Computation. In _33rd IEEE International Parallel and Distributed Processing Symposium (IPDPS’19)_ (Rio de Janeiro, Brazil). IEEE.  https://doi.org/10.25344/S4V88H
  * Berlin et al. (2015) Konstantin Berlin, Sergey Koren, Chen-Shan Chin, James P Drake, Jane M Landolin, and Adam M Phillippy. 2015\.  Assembling large genomes with single-molecule sequencing and locality-sensitive hashing.  _Nature Biotechnology_ 33, 6 (2015), 623–630.
  * Besta et al. (2020) Maciej Besta, Raghavendra Kanakagiri, Harun Mustafa, Mikhail Karasikov, Gunnar Rätsch, Torsten Hoefler, and Edgar Solomonik. 2020\.  Communication-efficient jaccard similarity for high-performance distributed genome comparisons. In _International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, IEEE Computer Society, 1122–1132.
  * Buluç and Gilbert (2008) Aydın Buluç and John R Gilbert. 2008.  On the representation and multiplication of hypersparse matrices. In _2008 IEEE International Symposium on Parallel and Distributed Processing_. IEEE, 1–11.
  * Buluç and Madduri (2011) Aydın Buluç and Kamesh Madduri. 2011.  Parallel breadth-first search on distributed memory systems. In _Proceedings of International Conference for High Performance Computing, Networking, Storage and Analysis (SC)_. 1–12.
  * Carlson et al. (1999) William W Carlson, Jesse M Draper, David E Culler, Kathy Yelick, Eugene Brooks, and Karen Warren. 1999\.  _Introduction to UPC and language specification_.  Technical Report. Technical Report CCS-TR-99-157, IDA Center for Computing Sciences.
  * Cheng et al. (2021) Haoyu Cheng, Gregory T Concepcion, Xiaowen Feng, Haowen Zhang, and Heng Li. 2021.  Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm.  _Nature methods_ 18, 2 (2021), 170–175.
  * Chin et al. (2016) Chen-Shan Chin, Paul Peluso, Fritz J Sedlazeck, Maria Nattestad, Gregory T Concepcion, Alicia Clum, Christopher Dunn, Ronan O’Malley, Rosa Figueroa-Balderas, Abraham Morales-Cruz, et al. 2016\.  Phased diploid genome assembly with single-molecule real-time sequencing.  _Nature methods_ 13, 12 (2016), 1050–1054.
  * Edmonds and Johnson (2003) Jack Edmonds and Ellis L Johnson. 2003.  Matching: A well-solved class of integer linear programs.  In _Combinatorial Optimization—Eureka, You Shrink!_ Springer, 27–30.
  * Eid et al. (2009) John Eid, Adrian Fehr, Jeremy Gray, Khai Luong, John Lyle, Geoff Otto, Paul Peluso, David Rank, Primo Baybayan, Brad Bettman, et al. 2009\.  Real-time DNA sequencing from single polymerase molecules.  _Science_ 323, 5910 (2009), 133–138.
  * Erbert et al. (2017) Marius Erbert, Steffen Rechner, and Matthias Müller-Hannemann. 2017\.  Gerbil: a fast and memory-efficient k-mer counter with GPU-support.  _Algorithms for Molecular Biology_ 12, 1 (2017), 1–12.
  * Georganas et al. (2015) Evangelos Georganas, Aydın Buluç, Jarrod Chapman, Steven Hofmeyr, Chaitanya Aluru, Rob Egan, Leonid Oliker, Daniel Rokhsar, and Katherine Yelick. 2015\.  HipMer: an extreme-scale de novo genome assembler. In _The International Conference for High Performance Computing, Networking, Storage and Analysis (SC)_. ACM, ACM, New York, NY, 14.
  * Georganas et al. (2018) Evangelos Georganas, Rob Egan, Steven Hofmeyr, Eugene Goltsman, Bill Arndt, Andrew Tritt, Aydın Buluç, Leonid Oliker, and Katherine Yelick. 2018.  Extreme scale de novo metagenome assembly. In _SC18: International Conference for High Performance Computing, Networking, Storage and Analysis_. IEEE, IEEE Computer Society, 122–134.
  * Ghosh et al. (2020) Priyanka Ghosh, Sriram Krishnamoorthy, and Ananth Kalyanaraman. 2020\.  PaKman: A Scalable Algorithm for Generating Genomic Contigs on Distributed Memory Machines.  _IEEE Transactions on Parallel and Distributed Systems_ 32, 5 (2020), 1191–1209.
  * Goodwin et al. (2015) Sara Goodwin, James Gurtowski, Scott Ethe-Sayers, Panchajanya Deshpande, Michael C Schatz, and W Richard McCombie. 2015.  Oxford Nanopore sequencing, hybrid error correction, and de novo assembly of a eukaryotic genome.  _Genome Research_ 25, 11 (2015), 1750–1756.
  * Graham (1966) Ronald L Graham. 1966\.  Bounds for certain multiprocessing anomalies.  _Bell system technical journal_ 45, 9 (1966), 1563–1581.
  * Graham (1969) Ronald L. Graham. 1969\.  Bounds on multiprocessing timing anomalies.  _SIAM journal on Applied Mathematics_ 17, 2 (1969), 416–429.
  * Guidi et al. (2021a) Giulia Guidi, Marquita Ellis, Aydın Buluç, Katherine Yelick, and David Culler. 2021a.  10 years later: Cloud computing is closing the performance gap. In _Companion of the ACM/SPEC International Conference on Performance Engineering_. 41–48.
  * Guidi et al. (2021b) Giulia Guidi, Marquita Ellis, Daniel Rokhsar, Katherine Yelick, and Aydın Buluç. 2021b.  BELLA: Berkeley efficient long-read to long-read aligner and overlapper. In _SIAM Conference on Applied and Computational Discrete Algorithms (ACDA21)_. SIAM, SIAM, 123–134.
  * Guidi et al. (2021c) Giulia Guidi, Oguz Selvitopi, Marquita Ellis, Leonid Oliker, Katherine Yelick, and Aydın Buluç. 2021c.  Parallel string graph construction and transitive reduction for de novo genome assembly. In _IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, IEEE Computer Society, 517–526.
  * Gurevich et al. (2013) Alexey Gurevich, Vladislav Saveliev, Nikolay Vyahhi, and Glenn Tesler. 2013.  QUAST: quality assessment tool for genome assemblies.  _Bioinformatics_ 29, 8 (2013), 1072–1075.
  * Hofmeyr et al. (2020) Steven Hofmeyr, Rob Egan, Evangelos Georganas, Alex C Copeland, Robert Riley, Alicia Clum, Emiley Eloe-Fadrosh, Simon Roux, Eugene Goltsman, Aydın Buluç, et al. 2020\.  Terabase-scale metagenome coassembly with MetaHipMer.  _Scientific reports_ 10, 1 (2020), 1–11.
  * Hon et al. (2020) Ting Hon, Kristin Mars, Greg Young, Yu-Chih Tsai, Joseph W Karalius, Jane M Landolin, Nicholas Maurer, David Kudrna, Michael A Hardigan, Cynthia C Steiner, et al. 2020\.  Highly accurate long-read HiFi sequencing data for five complex genomes.  _Scientific data_ 7, 1 (2020), 1–11.
  * Kolmogorov et al. (2019) Mikhail Kolmogorov, Jeffrey Yuan, Yu Lin, and Pavel A Pevzner. 2019.  Assembly of long, error-prone reads using repeat graphs.  _Nature biotechnology_ 37, 5 (2019), 540–546.
  * Koren et al. (2017) Sergey Koren, Brian P Walenz, Konstantin Berlin, Jason R Miller, Nicholas H Bergman, and Adam M Phillippy. 2017.  Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation.  _Genome research_ 27, 5 (2017), 722–736.
  * Li (2016) Heng Li. 2016.  Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences.  _Bioinformatics_ 32, 14 (2016), 2103–2110.
  * Medvedev et al. (2007) Paul Medvedev, Konstantinos Georgiou, Gene Myers, and Michael Brudno. 2007.  Computability of models for sequence assembly. In _International Workshop on Algorithms in Bioinformatics_. Springer, 289–301.
  * Mesnard and Barba (2019) Olivier Mesnard and Lorena A Barba. 2019.  Reproducible workflow on a public cloud for computational fluid dynamics.  _Computing in Science & Engineering_ 22, 1 (2019), 102–116.
  * Miller et al. (2008) Jason R Miller, Arthur L Delcher, Sergey Koren, Eli Venter, Brian P Walenz, Anushka Brownley, Justin Johnson, Kelvin Li, Clark Mobarry, and Granger Sutton. 2008\.  Aggressive assembly of pyrosequencing reads with mates.  _Bioinformatics_ 24, 24 (2008), 2818–2824.
  * Nagarajan and Pop (2009) Niranjan Nagarajan and Mihai Pop. 2009.  Parametric complexity of sequence assembly: theory and applications to next generation sequencing.  _Journal of Computational Biology_ 16, 7 (2009), 897–908.
  * Nisa et al. (2021) Israt Nisa, Prashant Pandey, Marquita Ellis, Leonid Oliker, Aydın Buluç, and Katherine Yelick. 2021.  Distributed-Memory k-mer Counting on GPUs. In _IEEE International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 527–536.
  * Nurk et al. (2020) Sergey Nurk, Brian P Walenz, Arang Rhie, Mitchell R Vollger, Glennis A Logsdon, Robert Grothe, Karen H Miga, Evan E Eichler, Adam M Phillippy, and Sergey Koren. 2020\.  HiCanu: accurate assembly of segmental duplications, satellites, and allelic variants from high-fidelity long reads.  _Genome research_ 30, 9 (2020), 1291–1305.
  * Phillippy et al. (2008) Adam M Phillippy, Michael C Schatz, and Mihai Pop. 2008\.  Genome assembly forensics: finding the elusive mis-assembly.  _Genome Biology_ 9, 3 (2008), R55.
  * Reed et al. (2022) Daniel Reed, Dennis Gannon, and Jack Dongarra. 2022\.  Reinventing High Performance Computing: Challenges and Opportunities.  _arXiv preprint arXiv:2203.02544_ (2022).
  * Selvitopi et al. (2020) Oguz Selvitopi, Saliya Ekanayake, Giulia Guidi, Georgios Pavlopoulos, Ariful Azad, and Aydın Buluç. 2020.  Distributed Many-to-Many Protein Sequence Alignment Using Sparse Matrices.  _The International Conference for High Performance Computing, Networking, Storage and Analysis (SC)_ (2020).
  * Shafin et al. (2020) Kishwar Shafin, Trevor Pesout, Ryan Lorig-Roach, Marina Haukness, Hugh E Olsen, Colleen Bosworth, Joel Armstrong, Kristof Tigyi, Nicholas Maurer, Sergey Koren, et al. 2020\.  Nanopore sequencing and the Shasta toolkit enable efficient de novo assembly of eleven human genomes.  _Nature biotechnology_ 38, 9 (2020), 1044–1053.
  * Zeni et al. (2020) Alberto Zeni, Giulia Guidi, Marquita Ellis, Nan Ding, Marco D Santambrogio, Steven Hofmeyr, Aydın Buluç, Leonid Oliker, and Katherine Yelick. 2020.  LOGAN: High-performance GPU-based X-drop long-read alignment. In _International Parallel and Distributed Processing Symposium (IPDPS)_. IEEE, 462–471.
  * Zhang et al. (2011) Wenyu Zhang, Jiajia Chen, Yang Yang, Yifei Tang, Jing Shang, and Bairong Shen. 2011\.  A practical comparison of de novo genome assembly software tools for next-generation sequencing technologies.  _PloS One_ 6, 3 (2011), e17915.