The genetic composition of individuals within the same species may difer due to mutations in their genomes. Variants calling procedures are pivotal to detect these polymorphisms. Typically, variants calling algorithms use the De Bruijn graph data structure for storing the k-nucleotide long strings sequenced from the genomes of multiple individuals. Subsequently, mutations in an individual can be identified by searching for divergent paths (known as bubbles) in the corresponding De Bruijn graph.
CEUR
ceur-ws.org SNPs.
However, experimental results have shown that scalability is limited due to the overhead of traversing paths scattered across diferent computers.
In this work, we address one of the main problems that prevent distributed De Bruijn graphs from performing well in practice. Namely, we analyze how a bad distributed partitioning of a De Bruijn graph may hinder the eficiency of the
SNPs detection process. Then, we introduce a distributed partitioning strategy based on a hierarchical clustering method able to overcome this problem. We also provide an experimental analysis showing how this approach allows for a relevant performance boost of the considered pipeline and ensures for better scalability.
Advances in New Generation Sequencing (NGS) technologies have led to an explosion in the amount of genomics data at afordable costs. Because of this massive production of biological datasets, pivotal tasks in Bioinformatics have become computationally expensive.
One of the most relevant tasks is variants calling. We use this term to denote a procedure aimed at detecting variants between the genomes of two or more diferent individuals. Referencebased algorithms are commonly used to detect the presence of mutations with respect to the average genome of the corresponding species. However, high-quality reference genomes are not available for all species and for all applications. For these reasons, it is useful to consider CEUR Workshop Proceedings reference-free algorithms. These allow comparison of individuals by evaluating reads sequenced from the corresponding genomes, without relying on an average species sequence.
The state-of-the-art reference-free methods are usually based on the exploration of a graphbased data structure, called De Bruijn (DB) graph, that is built from an input collection of reads. The presence of a polymorphism is then detected by finding some specific types of divergent paths in the corresponding DB graph. However, the computational and memory requirements of these algorithms can be prohibitive, making them unsuitable for processing large DB on a single workstation. Some alternative methods that have emerged recently therefore use distributed computing to develop reference-free pipelines for variants calling, able to take advantage of the large availability of memory and computational resources of a computer cluster.
Despite the theoretical advantages, analyzing a DB graph on a distributed architecture poses significant challenges. For example, the data trafic required to explore a region of the graph where each vertex of a path is on a diferent node of the underlying distributed system can create excessive overhead that negatively afects algorithm performance and scalability.
In this work, we focus on the specific problem of reducing the communication overhead to
pay for exploring distributed DB graphs by considering a recent distributed pipeline presented
in [
We also present the results of an experimental analysis performed on a real dataset to compare
our improved partitioning scheme with the default partitioning strategy used by [
The DNA structure consists of two antiparallel strands comprising nucleotide bases, commonly represented by the alphabet Σ = {, , , } . The strands run in opposite directions, and the bases in the corresponding positions are paired according to specific rules (A-T and C-G). Consequently, a genomic sequence is usually encoded as a single string of characters drawn from the alphabet Σ, as the reverse strand is easily obtainable. However, modern sequencing technologies do not return the entire genome in one continuous sequence. Instead, they generate huge amounts of nucleotide fragments, called reads. These reads are randomly sequenced from both strands, and the coverage level determines the number of reads associated with a particular genomic region.
Consequently, reconstructing a genome becomes a complex task due to the chaotic reads samples that need to be assembled. The assembly phase is extremely expensive from the computational viewpoint since it need to face the time complexity of computing the matching score for each pair of reads.
Assembly techniques based on k-mers and the DB graph have revolutionized genomic analysis. The k-mer-based approach involves splitting the genomic reads into smaller overlapping subsequences of length k (called k-mers) and building the corresponding DB graph. The vertices of a DB graph represent the distinct k-mers extracted from the reads sample, while the edges connect the k-mers that overlap for − 1 characters. The DB graph provides a compact and eficient representation of the genome, and traversing its path of overlapping sequences enhances an eficient reconstruction of the original genomic sequence.
DB graphs play also crucial role in variants calling. Variants calling refers to the process of identifying genetic variations between diferent individuals, such as single nucleotide polymorphisms (SNPs) and insertions/deletions (indels), in a given genomic dataset. The structure of the DB graph enables eficient and accurate detection of variants. Indeed, when dealing with sequencing data referring to two or more individuals, the presence of polymorphisms generate divergent paths, called bubbles. There exist diferent type of bubbles according to corresponding polymorphisms. In the case of isolated SNPs, we have simple bubble. Figure 1 shows an example. In real world scenario the topology of graph is extremely complex and the detection of the divergent is not straightforward.
DB graphs tend to become very large and, consequently, dificult to store, when based on
genomic datasets returned by modern sequencing technologies. A possible solution to this
problem consists of introducing compressed data structures for storing a DB graph. This problem
is addressed in [
Further solutions take into account a probabilistic representation of the DB graph. Among
the proposed techniques, [
Despite the potential advantages, the usage of distributed computing for large-scale variants
calling analysis still remains a relatively unexplored field. To the best of our knowledge, [
Computations on distributed systems enable the processing of huge collections of data. In distributed computing architectures, a number of computers (called nodes) work together to solve a single problem. The nodes coordinate their hardware resources using messages exchanged over a network.
The MapReduce (MR) paradigm [
Apache Spark [
Apache Spark also provides libraries (built on top of RDDs) dedicated to specific tasks. The
GraphX API [
The distributed pipeline we consider in this work has been recently proposed in [
In a distributed system, the amount of available memory is usually much larger, since we can
leverage on the computational resources of all nodes of the system. Therefore, in the approach
presented in [
Specifically, the proposed algorithm is formulated in three steps, each made of a sequence of distributed transformations. In the first step, a DB graph is created from from an initial distributed collection of reads. In the second step, is processed by a distributed algorithm that searches for all simple bubbles. In the third step, all isolated bubbles, i.e., the SNPs, are returned in the output. More details on the algorithm follow.
Given two FASTA/FASTQ files, each containing a set of reads extracted from two diferent
individuals, their contents are loaded in memory using the FASTdoop library [
A map transformation is used to filter out ( +1) -mers with low coverage. Then, after merging all samples, a DB distributed graph is created using the GraphX library (see Section 2.3.2). Next, we extract consecutive k-mers from the ( + 1) -mers to feed the distributed collection of vertices of the DB graph. Then, ( + 1) -mers are added to the distributed collection storing the edges of the DB graph. The edges connect the vertices corresponding to the two consecutive k-mers. This process is also performed for the backward strand of each ( + 1) -mer.
In this step, two centrality indices telling the number of incoming and outgoing edges for each vertex , respectively denoted inDegree and outDegree, are evaluated.
To do this, for each vertex , the number of triples containing in the distributed representation of is counted, distinguishing the cases where is the source vertex from those where is the destination vertex. Then, for each vertex, these counts are computed by a distributed reduce transformation to obtain the corresponding centrality indices.
In this step, a distributed algorithm is executed over to find all paths that represent simple bubbles.
Initially, vertices with an outDegree greater than 1 are considered as starting vertices and are marked as enabled. In the first iteration, each enabled triplet sends a path-building message to its adjacent vertices, including information about the current path (which initially contains only the starting point) and the iteration number.
In subsequent iterations, the triplets containing vertices that received messages in the previous iteration are enabled. When a vertex receives a path-building message from a neighbor, it updates the received message by adding its own identity and increasing the iteration number. The vertex stores a copy of the iteration number in its state and broadcasts the updated message to its adjacent vertices.
To handle branching bubbles, the algorithm modifies the transmission of path-building messages. In the first -1 iterations, the message is only transmitted to destination vertices with one incoming edge and one outgoing edge. Once the number of iterations reaches , the messages are forwarded without checking the degrees of the destination vertex. This change is because at iteration , the destination vertex is expected to be the end vertex of the bubble, allowing multiple incoming and outgoing edges.
At the end of the execution, after running iterations, the algorithm identifies vertices with iteration number + 1 that have received two messages. These messages contain paths starting from the same vertex and with the same length. These vertices represent terminal nodes of sequence pairs containing simple bubbles and are considered as isolated SNPs. They are encoded as a distributed collection of string pairs.
When using a distributed approach, one expects that the solution time of a given problem could
be reduced by just employing more computational resources. However, the ability to eficiently
exploit the computational capability of a distributed system is not always a given. This is either
due to the inherent non-parallelism of the problem being solved, or to issues involving the
algorithm itself, its implementation and the way it interacts with the underlying distributed
system. This is the case of the pipeline presented in [
A more detailed investigation revealed that this lack of scalability is mainly due to the suboptimal strategy GraphX uses for partitioning the vertices of the DB graph across diferent nodes of a distributed system, resulting in a very high number of x-cross paths, i.e. paths of a distributed graph where vertices are placed over two or more nodes of the distributed system.
In this scenario, the advantages of adding further computing nodes may be burdened by the communication overhead introduced by the message-passing mechanism required to explore the distributed graph. This overhead includes both the exchange of data between computing units within a node and across diferent nodes. The latter case is significantly more expensive in terms of network communication times and scalability level.
Our work focuses on assessing the positive impact on the performance of the considered algorithm coming from the adoption of a partitioning strategy diferent from the standard one available with GraphX. Our expectation is that the usage of an improved scheduling strategy would significantly reduce, when not eliminating at all, the number of x-cross paths. In turn, this would lead to a significant reduction in the communication overhead between diferent computing nodes and/or computing units, thus allowing for shorter overall execution times. To achieve this, our proposal leverages a hierarchical clustering approach to group reads referring to the same genomic region, enabling an intra-cluster solution for bubbles detection.
The strategy we propose to improve the partitioning of a DB graph over a distributed system is based on the assumption that clustering the reads according to an appropriate similarity metric isolates specific genomic regions, allowing to perform the bubble detection procedure mostly within them. Consequently, assigning these clusters to diferent independent computing units for parallel processing is expected to require less communication overhead than using the standard partitioning strategy, thus enhancing the scalability of the DB graph exploration in a distributed environment.
It requires two steps. In the Hierarchical Clustering step, the collection of input reads about
two distinct individuals are clustered using a bottom-up approach and a greedy stop criterion.
In the Clusters Binning step, the clusters returned in the previous step are grouped into distinct
bins scattered across the computational units of the distributed architecture. The generation of
the bins is driven by the application of the longest processing time first rule (i.e., rule)[
Given an input set of reads sequenced from two individuals, we use a hierarchical bottom-up algorithm to group those reads that are likely to cover the same genomic region.
An algorithm of this type implements an iterative procedure that progressively merges similar elements into larger clusters to eventually obtain a dendrogram, i.e., a hierarchical tree-like structure representing the clustering relationships. In the initialization step, each read is considered as a cluster with a single element. Then, in the -th iteration, for each cluster, the read closest to the others in the same cluster is selected as representative. Then, the pairwise similarity matrix between all clusters is calculated and the two closest clusters are merged. This process continues until obtaining a single large cluster that contains all the reads in the input dataset.
Creating the entire dendrogram can be very computationally intensive when working with large datasets. Moreover, solutions with a small number of clusters should be avoided, as our strategy aims to detect clusters consisting of very similar reads that refer to the same genomic region. To achieve this, we incorporated a greedy stopping criterion into the algorithm. Namely, our hierarchical iterative process stops when the similarity between the two clusters to be merged at the -th iteration is less than a certain threshold . In this way, the choice of the number of clusters is also automated without the need for a global evaluation of the dendrogram. If the threshold is suficiently large, hierarchical clustering yields a large number of clusters containing a relatively small number of reads. However, if is too large, reads related to the same genomic section may be assigned to diferent clusters, causing many bubbles to be missed during the detection step.
We recall that a genomic read is simply a string of characters from the alphabet Σ = {, , , } . Our clustering framework, therefore, requires a similarity measure that can capture the underlying similarity patterns between strings and is not extremely sensitive to slight variations in characters. This helps to ensure that reads sequenced from multiple individuals covering the same genomic regions are clustered together, despite the presence of SNPs in the same cluster.
In this work, we consider the Jaccard index, also known as Jaccard similarity coeficient . This
measure is widely used in computational genomics for assessing similarity between nucleotide
strings and for clustering [
The clusters returned in the previous step need to be distributed among the computational units of the distributed system to perform bubble detection in parallel. In this way, each computing unit processes a bin of clusters according to the available resources. However, the number of reads changes within the cluster, and some clusters may be overloaded due to the presence of highly repetitive genomic regions. For this reason, the binning strategy is controlled by an appropriate scheduling model, to ensure a balanced workload on the computational units.
We consider this problem as a scheduling problem with identical parallel machines and no preemption. The objective is to minimize the maximum completion time (makespan) by assigning each of the jobs, characterized by their respective processing times { , = 1, … , } , to one of the parallel identical machines. In our context, these jobs correspond to the tasks associated with building and exploring the local DB graph belonging to a cluster of reads. The processing time of each job is estimated based on the number of k-mers extracted from the reads within the same cluster. The machines, on the other hand, represent the available computational units. This problem, referred to as || mathematically as follows: , is known to be NP-hard and can be formulated (1) (2) (3) (4) min
≥ ∑=1
∀ ∈ {1..} ∑=1 = 1 ∀ ∈ {1.. } ∈ {0, 1} ∀ ∈ {1..} ∈ {1.. } where is the
, is the number of jobs, is the number of machines, is the processing time of job , is the decision variable. The value is 1 if the job is assigned to machine , otherwise it is 0. Given the complexity of this problem, many heuristic solutions have been proposed in literature, as the
rule. It involves sorting jobs in a non-increasing order based on their processing times and assigning each job iteratively to the machine that currently has the minimum completion time. The completion time of a machine is determined by summing the processing times of all the jobs assigned to that particular machine.
We conducted an experimental analysis to assess the quality of the proposed clustering-driven partitioning strategy and evaluate its impact on the scalability.
We downloaded the FASTA file corresponding to chromosome 7 of the average human genome
assembly GRCh37.p13 [
The individual NA12878 is accompanied by a VCF file containing the set of high-confidence variants with respect to the average GRCh37.p13 genome. From this VCF file, we selected 10, 000 SNPs for our analysis. Thereafter, we extracted the portions of the DB graph corresponding to these SNPs, to compare the GraphX default partitioning strategy with the one outcoming from our cluster-driven partitioning strategy. From now on, we refer to the filtered DB graph used in our analysis as DB1e4.
We conducted our experiments on a high-performance computing system that consists of eight computing nodes. Each node runs on the Linux operating system and is equipped with two AMD Epyc 7452 processors, 64 compute cores, and 256 GB of RAM. For our experiments, we employed the 3.1.3 Apache Spark version and the 2.7 Apache Hadoop version.
First, we evaluated the GraphX default partitioning strategy on DB1e4. Thanks to the information provided by the VCF file, containing mutations about the individual NA12878, with respect to the GRCh37.p13 version of the average human genome, we easily identified the bubbles occurring in the graph. For each bubble, we determined the number of computing units that need to communicate in order to traverse each edge of the bubble. Since we set = 31 , each bubble may be split into up to 31 diferent computing units.
In this case, the histogram shows that there is no bubble that can be traversed by staying on a single computing unit. Instead, it is required to involve at least two computing units while, the most frequent case, requires the involvement of three computing units.
We compared the time-performance of the distributed pipeline described in Section 3, when using the GraphX default partitioning strategy against version based on our cluster-driven partitioning scheme, in terms of scalability.
First, we performed the SNPs calling procedure on DB1e4 using the aforementioned pipeline and the default GraphX partitioning strategy, while increasing the number of computing units. Then, we repeated the same experiment, but using our hierarchical clustering framework and binning strategy. In this last case, we considered for each run a number of bins that was 4 times the number of allocated computing units. The resulting bins were then loaded into memory on the distributed computing architecture, to perform the bubble detection procedure.
Figure 4 compares the scalability deriving from the application of the two partitioning approaches, highlighting how our proposed strategy significantly improves the ability to leverage on an increasing number of computing resources.
In this work, we focused on a performance bottleneck afecting Spark-based distributed pipelines used for SNPs detection and due to the bad performance of the default partitioning strategy employed by the Spark framework when creating distributed graphs. Thus, we have proposed an improved partitioning strategy, based on a hierarchical clustering algorithm. According to our experimental results, the proposed strategy is able to significantly improve to enhance the scalability of the considered distributed computing pipeline for SNPs calling. In future directions, it will be pivotal to explore diferent similarity measures and algorithms for genomic reads clustering, as the Jaccard index and hierarchical clustering algorithms can be computationally expensive when dealing with huge-scale analyses. Furthermore, ongoing research is evaluating strategies for detecting more complex mutations and dealing with genomic regions characterized by a high mutation intensity, whose corresponding reads may not fall within the same cluster.
The authors would like to thank the Department of Statistical Sciences of University of Rome - La Sapienza for computing time on the TeraStat cluster. This work was partially supported by Università di Roma - La Sapienza Research Project 2021 “Caratterizzazione, sviluppo e sperimentazione di algoritmi eficienti”. It was also supported in part by INdAM – GNCS Project 2023 “Approcci computazionali per il supporto alle decisioni nella Medicina di Precisione”.