QUARTIC: QUick pArallel algoRithms for high-Throughput sequencIng data proCessing

Alignment

For the alignment, the software consists of a MPI wrapper to the BWA-MEM software such that the original code remains unchanged. The wrapper uses parallel IO and shared memory (Figure 1). It first parses the input FASTQ files in order to define chunks of equal size. The different chunks are actually represented as two offsets (start and end) from the original FASTQ file thus reducing the amount of information to store. The reference genome is loaded into the shared memory on each computer node. Then, each chunk is processed in parallel threads by the original BWA-MEM algorithm. The results is finally written thanks to a shared file pointer in a unique SAM file. For a human genome, the size of the reference genome file in the shared memory takes around 5 GB and the total memory used by a BWA-MEM thread is around 300 MB. The original BWA-MEM processes a chunk of reads at a time that contains the same number of nucleotides which is also the case with our parallel implementation in order to allow the reproducibility with the original algorithm. We name our code alignment mpiBWA (Hupé & Jarlier, 2020a).

Figure 1. Alignment with mpiBWA using two computing nodes and four cores per node.

Each core is in charge of aligning a chunk of the FASTQ file. The reference genome is stored in shared memory of each node of the computing cluster. Once aligned, the results are written into a SAM file. To ensure scalability, the read and write operation require a parallel filesystem.

Sorting

For the sorting, the software implements a parallel version of the bitonic sort proposed by Batcher (1968) for sorting any sequence of elements of size n = 2^k being of power of 2. Its complexity is 𝒪 (n log²n) that is higher than 𝒪 (n log n) from the popular merge sort algorithm. However, the bitonic sort is very suitable for parallel implementation since it always compares elements in predefined sorting network that is independent of the input data. To understand how the algorithm works, we first defined in what follows some concepts (see Grama et al. (2003) for details).

The algorithm relies on a bitonic sequence that is a sequence of values 〈a₀, a₁, … , a_n−1〉 with the property that 1) there is an index i, 0 ≤ i ≤ n−1 such that 〈a₀, a₁, … , a_i〉 is monotonically increasing and 〈a_i+1, … , a_n−1〉 is monotonically decreasing, or 2) there exists a cyclic shift of indices so that the condition 1) is satisfied. From any bitonic sequence s = 〈a₀, a₁, … , a_n−1〉, the bitonic split operation consists in transforming the input sequence s into these two subsequences:

Batcher (1968) proved that both s₁ and s₂ are bitonic sequences and the elements of s₁ are smaller than the elements of s₂. Thus, a recursive bitonic split from a bitonic sequence of size n = 2^k , until the sequence obtained are of size one allows the sorting of the input bitonic sequence in k splits as shown in Figure 2. The procedure of sorting a bitonic sequence from bitonic split operation is called bitonic merge. For each split of a bitonic merge, n/2 comparisons are performed during which the two numbers are exchanged if not in the right order using a compare-exchange operation. A bitonic merge procedure of a sequence of size n is noted $\text{B}{\text{M}}_n^{\oplus }$ if the comparisons sort the number in monotonically increasing order or $\text{B}{\text{M}}_n^{\ominus }$ for monotonically decreasing order.

Figure 2. Sorting network for the bitonic sort.

Each horizontal line represents an input and each arrow represents a two-input/two-output comparator which performs compare-exchange operation. Depending on its direction, the outputs are sorted in increasing (decreasing) order with an downwards (upwards) arrow. First, the comparator network transforms an input sequence of 16 unordered numbers into a bitonic sequence alternating $\text{B}{\text{M}}_2^{\oplus }/\text{B}{\text{M}}_2^{\ominus },$ $\text{B}{\text{M}}_4^{\oplus }/\text{B}{\text{M}}_4^{\ominus }\text{and}$ $\text{B}{\text{M}}_8^{\oplus }/\text{B}{\text{M}}_8^{\ominus }$ bitonic merges. Then, the obtained bitonic sequence is sorted with a recursive bitonic merge procedure $\text{B}{\text{M}}_{16}^{\oplus }$ (adapted from Grama et al. (2003)).

Sorting any sequence of unordered number thus requires to convert the input sequence into a bitonic sequence. This is obtained using a bitonic sorting network (see Figure 2) that sorts n = 2^k numbers using k − 1 bitonic sorting stages, where the i-th stage is composed of n/2ⁱ alternating increasing $\text{B}{\text{M}}_{2^i}^{\oplus }$ and decreasing $\text{B}{\text{M}}_{2^j}^{\ominus }$ bitonic merges. The final stage being $\text{B}{\text{M}}_n^{\oplus }$ to obtain the sorted sequence.

It is straightforward to generalize the bitonic sort algorithm on parallel architecture as shown in Figure 3. Each processor is in charge of an even number m = n/p elements (where p is the number of processors that is a power of 2) that are first sorted using a merge sort algorithm. Then, each comparison over the bitonic sorting network is replaced by a pair of processors with performs a compare-split operation. During the compare-split operation, the two sorted sequences from each processor are merged into one monotonic sorted list, then bisected into two (lower and higher) sequences. After the compare-split, one processor will keep the lower m elements from sequence and the other processor will keep the higher m elements according to the direction of the arrow in the step (see Kim et al. (2001) & Grama et al. (2003) for details). This parallelization allows the sorting of n elements in 𝒪 (log²n) time using a bitonic sorting network.

Figure 3. Parallel bitonic sort of 16 elements with four processors.

(a) Each horizontal line represents an input of elements attached to a processor. Each processor first locally sorts its list of 4 elements, then, the parallel sort performs three steps of compare-split operations over the bitonic sorting network. (b) During the compare-split operation, each pair of processors first exchanges their information: each processor sends its block to the other process, then it merges the two sorted blocks into one list and stores only the appropriate half of the merged block to be used in the next step (adapted from Kim et al. (2001) & Grama et al. (2003)).

The workflow for the sorting is described in Figure 4. The SAM file is read and split into p blocks that are dispatched across the p processors. Each block is parsed such that for each line of the SAM file, we extract the genomic coordinates of each read, its sequence (and all other information contained in the SAM file) and the line is indexed according to its offset in the SAM memory buffer of a processor. Then, the genomic coordinates are sorted with the bitonic sort as shown in Figure 3. During the sorting, a vector with five values follows the bitonic sorting network including the genomic coordinate c, the rank r_i of the processor that parsed the block of the input SAM file, the offset o_i of the line from the input SAM file, the rank r_d of the processor that will write the block of the sorted data in the destination SAM file and the offset o_d of the line in the sorted destination SAM file. Obviously, the values r_d and o_d are known when the bitonic sort is completed. It is important to highlight that the parallel computation is performed by different processors that are located on different compute nodes. This means that a block of data that has been read by a given processor is not accessible by another processor. Moreover, to optimize the writing of the sorted destination SAM file, it is essential to write the data in contiguous block. Thus, it is necessary that a processor in charge of the writing owns locally the data from a block of contiguous offsets o_d in the sorted destination file. This implies that all the data have to be shuffled across all the processors that have to exchange data in an all-to-all communication step. To optimize the communication between the processors during the shuffle, we have implemented a Bruck algorithm (Bruck et al., 1997) of time complexity 𝒪 (log² p). The Bruck algorithm is performed twice as shown in Figure 4 and can be seen as a joint procedure between two tables (the elements of the table being located on different processors). During the first Bruck phase, the r_d and o_d are sent back to the corresponding processor from where the data originated, and then, during the second Bruck phase, all the data (i.e. the appropriate lines of the SAM file) are sent to the processor that has been assigned to write the data. During the second phase the original data is copied into the memory buffer and then exchanged.

Figure 4. All-to-all communications with the Bruck algorithms for the sorting of 16 reads on four processors.

The SAM file is read and split over the 4 processors P₀ to P₃ that sort the data with the bitonic sort. Then, a first Bruck phase sent the r_d and o_d values back to processor from where the data originated. The data from the SAM file are sent by a second Bruck phase to the processor that has been assigned to write the contiguous blocks (during this step, the original data is copied into the memory buffer and then exchanged). Finally, the data can be written on a parallel filesystem. c: genomic coordinate. r_i: rank of the processor that parsed the block of the input SAM file. o_i: offset of the line from the input SAM file. r_d: rank of the processor that will write the block of the sorted data in the destination SAM file. o_d: offset of the line in the sorted destination SAM file.

Note that when the SAM file contains several chromosomes, each chromosome is sorted successively. Therefore, the upper memory bound depends on the size of the whole SAM file plus internal structures plus the size of the biggest chromosome. In this case, the total amount of memory used is around 1.5 the size of the original SAM file if the file contains read on all the chromosome of the human genome. If the SAM file contains only one chromosome, then the total amount of memory required is 2.5 the size of the original SAM file. To be efficient, the sorting requires a number of cores being a power of 2. We name our code for sorting mpiSORT (Hupé & Jarlier, 2020b).

Benchmark

We benchmarked mpiBWA and mpiSORT on the HG001 NA12878 sample from GIAB (Zook et al., 2014). This sample is a whole genome with 300X depth of coverage composed of 2.13 billions 2x250 pair-ended reads. The BAM file has been downloaded from:

url: ftp://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/NA12878/

folder: NIST_NA12878_HG001_HiSeq_300x/NHGRI_Illumina300X_novoalign_bams/

The BAM file has been converted into FASTQ files using bedtools bamtofastq. In order to test the scalability of our software with respect to the size of the data, we downsampled the original 300X sample to obtain FASTQ files corresponding to depths of coverage ranging from 28X to 300X with a geometric progression with a common ratio of 1.6. The sequences have been aligned on the human genome GRCh38.

We ran the benchmark on a computing cluster equipped with Intel Xeon^® Gold 6148 CPU @ 2.40GHz (Skylake architecture). The nodes are interconnected with Intel^® Omni Path Architecture (OPA) of 100 Gbps speed. The parallel file system is BeegFS with two servers. We compile the programs with GCC 8.3 and use Open MPI 3.1.4. Jobs have been submitted using slurm scheduler.

Implementation

The code has been written in C programming language using MPI directives. mpiBWA encapsulates the original BWA-MEM version 7.15. Two implementations for mpiBWA exist: the first outputs one SAM file with all the chromosomes, the second (named mpiBWAByChr) outputs one SAM file per chromosome. Moreover, mpiBWA comes along with mpiBWAIdx, which is responsible for creating a binary image of the reference genome. This image is subsequently loaded in the shared memory by mpiBWA. Then, every mpiBWA process on the same computing node will share the same genome reference in order to save memory usage. mpiBWAIdx does not need MPI to run.

Operation

As our software rely on the MPI standard, mpirun must be available to run the program. Several MPI implementations exist such as mpich, open-mpi or Intel^® MPI Library.

mpiBWA

The first step consist in building the binary image of the reference genome that will be used in the shared memory by mpiBWA:

mpiBWAIdx hg19.small.fa

This creates the file myReferenceGenome.fa.map.

Importantly, mpiBWAIdx requires the following files that are generated by BWA index:

  •  myReferenceGenome.fa.sa
  •  myReferenceGenome.fa.bwt
  •  myReferenceGenome.fa.ann
  •  myReferenceGenome.fa.pac
  •  myReferenceGenome.fa.amb

The .map file needs only to be generated once. Then, mpiBWA is executed with the mpirun program, for example:

mpirun -n 2 mpiBWA mem -t 8 -o ./HCC1187C.sam hg19.small.fa \
    HCC1187C_R1_10K.fastq \
    HCC1187C_R2_10K.fastq

This command launches two processes MPI with eight BWA threads each (16 cores will be used). If you want to split the results of the alignment by chromosome, use mpiBWAByChr, for example:

mpirun -n 2 mpiBWAByChr mem -t 8 -o ./HCC1187C.sam hg19.small.fa \
    HCC1187C_R1_10K.fastq \
    HCC1187C_R2_10K.fastq

The total memory used during the alignment if approximately the size of the .map file plus the size of the SAM chunks loaded by each BWA tasks. A BWA thread takes around 300 MB of memory.

mpiSORT

From the results obtained from mpiBWAByChr, each SAM files can be sorted as follows:

mpirun -n 4 mpiSORT -u chr1.sam ${HOME}/mpiSORTExample

This command launches four processes MPI.

As mpiSORT requires all the input SAM file to be loaded into the memory, the program is memory bounded. Therefore, a lot of attention has to be paid to define how many cores are needed to process the data. For example, let’s assume that the computing cluster consists of nodes with 190 GB of RAM memory with 40 cores each (thus 4.75 GB per core is available but this value can be rounded to the lower limit of 4.5 GB to leave some free memory on the node). To choose the number of cores, the following rule has to be applied: the total memory for sorting a SAM file that contains only one chromosome is around 2.5 times the size of the SAM file. For instance, sorting a chr1.sam file of size 209 GB with 4.5 GB per core requires 128 cores, for a of size 110 GB it requires 64 cores. We remind that the bitonic sort algorithm requires a number of cores that is a power of 2, this is the reason why the number of cores has to be set to the upper bound to the closest power of 2.