BigSeqKit: a parallel Big Data toolkit to process FASTA and FASTQ files at scale

GigaScience , 2023, 12 , 1–12 DOI: 10.1093/gigascience/giad062 Technical Note BigSeqKit : a parallel Big Data toolkit to process FASTA and FASTQ files at scale César Piñeiro and Juan C. Pichel * CiTIUS, Universidade de Santiago de Compostela, Santiago de Compostela 15782, Spain ∗Corr espondence addr ess . J uan C . Pichel, CiTIUS , Univ ersidade de Santia go de Compostela, Santia go de Compostela 15782, Spain. E-mail: juancarlos.pic [email protected] Abstract Bac kgr ound: High-throughput sequencing technologies have led to an unprecedented explosion in the amounts of sequencing data av aila b le , whic h ar e typicall y stor ed using F ASTA and F ASTQ files. W e can find in the literatur e sev eral tools to pr ocess and manipulate those type of files with the aim of transforming sequence data into biolog ical knowledg e . How ever, none of them are well fitted for processing efficiently very large files, likely in the order of terabytes in the following years, since they are based on sequential pr ocessing. Onl y some routines of the well-known seqkit tool are partly parallelized. In any case, its scalability is limited to use few threads on a single computing node. Results: Our approach, BigSeqKit , takes adv anta ge of a high-performance computing–Big Data framework to parallelize and optimize the commands included in seqkit with the aim of speeding up the manipulation of F ASTA/F ASTQ files. In this way, in most cases, it is from tens to hundreds of times faster than several state-of-the-art tools. At the same time, our toolkit is easy to use and install on any kind of hardw ar e platform (local server or cluster), and its routines can be used as a bioinformatics library or from the command line. Conclusions: BigSeqKit is a v er y complete and ultra-fast toolkit to process and manipulate large F ASTA and F ASTQ files. It is pub licl y av aila b le at https://github.com/citiususc/BigSeqKit . Ke yw ords: F ASTA/F ASTQ files, Performance, Parallelism, Big Data Introduction The history of modern DNA sequencing started se v er al decades ago and, since then, has seen astounding growth in sequencing capacity and speed. From the first genomes with a few thousand bases, DNA sequencing has advanced to sequence the human genome of 3 billion bases. In recent years, next-generation sequencing (NGS) technology, also known as massive parallel sequencing (MPS), has made it possible to expand the amount of sequencing data a vailable . For example , the Illumina No vaSeq 6000 [ 1 ] platform can generate a maximum output of 6 Tb of data and read about 20 billion sequences per run. Note that sequences, commonl y named reads , ar e composed of ASCII c har acters r epr esenting a nucleotide (base) from the sequence. In the DNA case, we can only find 4 possible bases (A—adenine , C—cytosine , G— guanine, and T—thymine). T he NGS ra w data are mainly stored in FASTA [ 2 ] and FASTQ [ 3 ] text-based file formats. In particular, nucleotide and protein sequences are typically stored in the FASTA file format, whereas FASTQ is the most widely used format for sequencing read data. An example of FASTA file is shown in Fig. 1 . A sequence in FASTA format begins with a single-line description about the sequence in the subsequent lines . T he description line is distinguished from the sequence data by a greater-than ( > ) symbol at the beginning. On the other hand, the FASTQ format was designed to handle the quality metrics of the sequences obtained from the sequencers. In FASTQ, e v ery 4 lines describe a sequence or read. An example is displa yed in Fig. 2 . T he information pro vided per read is as follows: identifier and an optional description (first line), sequence (second line), and the quality score of the read (fourth line). An extra field, r epr esented by symbol “+ ,” is used as separator between the data and the quality information (third line). Manipulating these files efficiently is essential to analyze and inter pr et data in any genomics pipeline. Common operations on F ASTA and F ASTQ files include searching, filtering, sampling, deduplication, and sorting, among others. We can find se v er al tools in the liter atur e for FASTA/Q file manipulation such as HTSeq [ 5 ], FASTX [ 6 ], fqtools [ 7 ], seqtk [ 8 ], Biopython [ 9 ], samtools [ 10 ], p yfadix [ 11 ], p yfastx [ 12 ], and seqkit [ 13 ]. These tools can be classified according to how the sequences are parsed [ 12 ]. In the first category, sequences ar e pr ocessed in order, whic h causes important overheads when extracting and randomly sampling sequences . T hat is the case of HTSeq , FASTX , fqtools , and seqtk . In the second category, we find tools that support random access to sequences by establishing an index file. Tools belonging to this category are more efficient in terms of performance and memory consumption. Ho w e v er, none of them are well fitted for processing v ery lar ge files of hundr eds of GB (likel y TBs in the near future) since they are based on sequential processing. The exception is seqkit that allows some routines to use a few threads, but in any case, its scalability is very limited. To deal with this issue, in this article, we introduce BigSeqKit , a parallel toolkit to manipulate FASTA and FASTQ files at scale with speed and scalability at its core. BigSeqKit takes adv anta ge of IgnisHPC [ 14 , 15 ], a computing engine that unifies the de v elopment, combination, and execution of high-performance computing (HPC) and Big Data parallel tasks using different languages Recei v ed: Mar c h 7, 2023. Revised: May 25, 2023. Accepted: July 10, 2023 ©The Author(s) 2023. Published by Oxford Uni v ersity Pr ess GigaScience. This is an Open Access article distributed under the terms of the Cr eati v e Commons Attribution License ( https://creativecommons.org/licenses/by/4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is pr operl y cited. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 2 | GigaScience , 2023, Vol. 12, No. 1 Figure 1: Example of FASTA file showing the first part of the PAX6 gene (obtained from [ 4 ]). Figure 2: Example of FASTQ file format (obtained from [ 4 ]). and pr ogr amming models. As it was demonstrated, IgnisHPC outperforms the state-of-the-art fr ame work Spark [ 16 ] in terms of performance and scalability running applications that represent the most typical algorithmic patterns in Big Data and scientific computing. BigSeqKit uses the seqkit routines as basis since that toolkit covers a wide range of utilities and is one of the most used by the bioinformatics r esearc h comm unity. As a consequence, BigSeqKit will offer the same functionalities and command interface [ 17 ]. BigSeqKit can be used from the command line, but it is at the same time a library, so its routines can also be called from a C/C ++ , Python, Go, or Java application. Another important c har acteristic of BigSeqKit is that it is fully containerized, which isolates the execution environment from the physical system and av oids dependenc y problems. As a consequence, BigSeqKit is very easy to install and can run on a local server or on any type of cluster since it supports some of the most important resource and scheduler managers (e.g., Mesos [ 18 ], Nomad [ 19 ], and Slurm [ 20 ]). Bac kgr ound IgnisHPC [ 14 , 15 ] unifies the execution of Big Data and HPC workloads in the same computing engine. Unlike other fr ame works such as Hadoop [ 21 ] and Spark [ 16 ], IgnisHPC has native support for m ultilangua ge a pplications using both JVM (Jav a Virtual Machine) and non-JVM-based languages. In this way, applications can be implemented using 1 or se v er al pr ogr amming langua ges following an API inspired by Spark’s one. The pr e vious v ersion of IgnisHPC supported nativ el y C, C ++ , J ava, and Python. Ho w ever, seqkit w as implemented using the Go pr ogr amming langua ge. Since BigSeqKit par allelizes and optimizes the seqkit routines using IgnisHPC, it was necessary to add support for this language in the fr ame work. Another solution would r equir e to port the complete toolkit to a different langua ge, whic h is a difficult task prone to errors. It is worth noting that, to the best of our kno wledge, no w adays IgnisHPC is the first parallel computing framework to include native support for this language. Considering Spark instead of IgnisHPC is not an option because, as it was demonstrated in [ 14 ], when using a nonnative language code, data transfers between the JVM and external pr ocesses degr ade noticeabl y Spark’s ov er all performance. Go is a pr ogr amming langua ge with a simple syntax that was designed to be easy to learn and use. With the release of Go v1.18, the language included support for Generics, which allows the creation of functions , types , and methods that can work with any data type . T his mak es Go an effecti v e and user-friendl y way to implement Big Data interfaces . T he implementation of Go in IgnisHPC is similar to that of C ++ , as both are compiled and statically typed languages. Ho w ever, Go replaces the concept of inheritance with composition, which does not change the philosophy of use in IgnisHPC. Big Data functions are still accessible through the IgnisHPC API, and users can create their own code by implementing the same interfaces. One of the k e y features of IgnisHPC is its use of containers to isolate and execute code. Containers are lightweight and portable, making it easy to run IgnisHPC on a variety of different clusters, including both HPC and Big Data. IgnisHPC is also tolerant to failures, as the containers or processes can be easil y r estarted if ther e are issues. In particular, if some data are lost, IgnisHPC has enough information about how it was derived. In this way, only those operations needed to recompute the corresponding portion of data are performed. We must highlight that although the IgnisHPC API [ 22 ] uses a sequential notation, operations on data are performed in parallel. As we pointed out, the IgnisHPC API was inspired by the Spark API in such a way that IgnisHPC codes are easily understandable by users who are familiar with Spark. Table 1 shows a list of some of the most important functions supported by IgnisHPC. In particular: r Map functions: The common c har acteristic to r outines belonging to this type is that they a ppl y the same function to each element in the data. As a result of the transformation, the output could be of different size with respect to the input. r Reduce functions: reduce and treeReduce methods a ggr egate all the elements in the input data using a function. aggregate and treeAggregate are a sort of reduction where the type of the input and output data is different. In this case, 2 functions are necessary; the first one is applied to each Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 BigSeqKit | 3 Table 1: Some of the most important IgnisHPC API functions Type Functions Map map, flatmap, mapWithIndex, filter, keyBy, keys, values, mapPartitions, mapValues, etc. Reduce reduce, treeReduce, aggregate, treeAggregate, reduceByKey, aggregateByKey, etc. Group groupBy, groupByKey Sort sort, sortBy, sortByKey I/O parallelize,collect, top, take, saveAsObjectFile, saveAsTextFile, saveAsJsonFile, etc SQL union, join, distinct Math sample, sampleByKey, take, takeSample, count, countByKey, countByValue, max, min, etc. Balancing repartition, partitionByHash, partitionByRandom, partitionBy Persistence persist, cache, unpersist, uncache element in a data partition, and the second one combines the partial results obtained for each partition. reduceByKey and aggregateByKey are variations where the operation is performed only among elements with the same k e y in such a way that the final result is a set of unique pairs with values calculated using reduce or aggregate operations, r espectiv el y. r Group functions: These methods group elements in a data frame according to their k e y value ( groupByKey ) or a userdefined function ( groupBy ). r Sort functions: In order to sort elements , IgnisHPC pro vides 3 functions: sort , sortByKey , and sortBy . The first method uses the natural order and does not need any additional function. sortByKey sorts the k e ys using their natural order. sortBy allows to use a user-defined function to specify the order of the elements. If the result of applying that function to 2 elements is true , then the first element should precede the second one. All methods support ascending and descending order. r SQL functions: These functions operate on data frames. union concatenates 2 data frames, join merges elements of 2 data frames whose k e ys match, and distinct returns a new data frame after removing the duplicate records . T hese methods are necessary, for example, in many graph processing problems. r Other functions: IgnisHPC implements se v er al oper ations that return a value to the driver code, but they do not modify or gener ate ne w stor ed data. Spark r efers to this type of operations as actions . For instance, IgnisHPC supports methods such as count , take , takeSample , and collect . The most basic operation is count that returns the number of elements of a stored data collection. collect returns a collection with all the elements stored in the executors of a task. take applies a collect operation but obtains only the first n elements, where n is chosen by the user. takeSample returns a random sample of n elements from the distributed data, with or without r eplacement. Finall y, another inter esting r outine is parallelize , which distributes the elements of a collection among the executors to form a distributed dataset. In this case, new stored data are created. It is worth noting that the IgnisHPC API functions allow users to parallelize a code with a high le v el of abstraction. In this way, it is only necessary to focus on data dependencies. Table 2: List of commands included in both BigSeqKit and seqkit . Those commands with an asterisk support new functionalities not included in seqkit Basic commands seq Transform sequences (extract ID, filter by length, r emov e ga ps, r e v erse complement, etc.) subseq Get subsequences by region/gtf/bed, including flanking sequences stats Simple statistics of FASTA/Q files: #seqs, min/max length, N50, Q20%, Q30%, etc. faidx ∗Create FASTA or FASTQ index file and extract subsequences Forma t con version fa2fq Retrie v e corr esponding FASTQ r ecor ds b y a FASTA file fq2fa Convert FASTQ file to FASTA format translate Translate DN A/RN A to protein sequence Searching grep Search sequences by ID/name/sequence/sequence motifs locate Locate subsequences/motifs Set oper a tions sample Sample sequences by number or proportion rmdup Remove duplicated sequences by ID/name/sequence common Find common sequences of multiple files by ID/name/sequence duplicate Duplicate sequences N times head Print first N FASTA/Q records head-genome Print sequences of the first genome with common prefixes in name pair Match up paired-end reads from 2 FASTQ files range Print FASTA/Q records in a range (start:end) Edit concat Concatenate sequences with the same ID from multiple files replace Replace name/sequence using a regular expression rename Rename duplicated IDs Ordering sort Sort sequences by ID/name/sequence/length shuffle Shuffle sequences Methods As we commented pr e viousl y, BigSeqKit ( RRID:SCR _ 023592 ) speeds up the seqkit routines through parallelization and optimization techniques. Table 2 shows the routines supported by the current version of BigSeqKit . Despite most of the commands in seqkit are sequential, we can classify each command implementation into 3 categories according to its inherent parallelism: r Independent: it is a embarr assingl y par allel workload. As a consequence, the computation could be applied to all sequences in parallel. An example is seq , a function that transforms sequences. In this case, the tr ansformation onl y affects each sequence individually. r P artiall y dependent: computations could be done in parallel, but the method r equir es some type of consensus to obtain the result. For instance, stats should merge the partial results computed for each sequence to calculate some statistics of the considered FASTA/Q file. r De pendent: de pendencies between sequences pr e v ent the method from being executed in parallel. As a consequence, Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 4 | GigaScience , 2023, Vol. 12, No. 1 BigSeqKit r equir es a complete ne w algorithm to perform the same command in parallel. rmdup is a good example because with the aim of removing duplicated sequences, it is necessary to read all of them before generating a result. The integr ation, par allelization, and optimization of each seqkit command in IgnisHPC will be different depending on its category. More details are provided below. Independent routines For these commands, the computation can be applied to all sequences in parallel because there are no dependencies (communication) among them. In other w or ds, routines belonging to this category can be processed using an embarrassingly parallel appr oac h. Considering the IgnisHPC (and Spark) API, it is only necessary to use map functions to parallelize the computations. As we pointed out, the common c har acteristic to these API functions is that they a ppl y the same operation to each element in the data. The following BigSeqKit commands belong to this category: seq , subseq , stats , fq2fa , fa2fq , translate , grep , locate , duplicate , and replace (see Table 2 for details). Partially dependent routines As we mentioned, this category includes commands in which computations can be done in parallel using map functions, but the methods r equir e some type of consensus to get the desired outcome . T his consensus can be easily implemented using the IgnisHPC API. The following BigSeqKit commands belong to this category: r stats : statistics can be generated in parallel but the final r esult m ust be unique, so all partial r esults m ust be merged using a reduction ( reduce operation in the IgnisHPC API). r head : sequences should know their position inside the file to c hec k if they ar e inside the head window . T o do that, it is necessary to use mapWithIndex , a special map operation included in the IgnisHPC API that allows each element to know its global index within a data structure. r head-genome : similar to head , but not all sequences are valid. In order to determine the window, invalid sequences must be r emov ed first. r range : also similar to head . Sequences should know their position inside the file to c hec k if they are within the range window. r grep : although this command was included in the pr e vious category, a command option ( --delete-matched ) limits the number of results to just 1 per search pattern. In such cases, it is necessary to r emov e the extra results. r faidx : also similar to head , sequences compute their offsets inside the input file using mapPartitionWithIndex and exchange the information between executors to perform a parallel indexing operation with a simple map . Dependent routines Commands belonging to this category have an implementation in seqkit that by its nature cannot be parallelized. Ho w ever, IgnisHPC allows us to define the implementation at a high le v el, whic h incr eases noticeabl y the pr oductivity. Behaviors and functionalities will be pr eserv ed in BigSeqKit but thr ough a complete ne w par allel implementation. In particular: r sample : a sequential sampling can be performed in parallel if we split the sequences and run a sample for each partition. It was mathematically proven that sampling without replacement follows a hypergeometric function [ 23 ]. In this way, we can calculate the proportion of the sample that corresponds to each partition. r rmdup : sequences ar e gr ouped ( groupBy API function) using a hash with the ID, name, or sequence. In those groups containing more than 1 element, a search for duplicates is carried out to r emov e them. r pair and concat : sequences of the input files generate k e y– v alue pairs wher e the k e y is the ID and the value is the sequence with its index file ( map ). P airs ar e unified by means of union and grouped using groupByKey . Afterw ar d, sequences in the same group are paired or concatenated if they belong to different files. r common : the first stage of the command is the same one explained above for pair and concat . Then if a sequence can be found in all files, we c hec k its index file, to avoid its deletion. r rename : sequences ar e gr ouped ( groupBy ) using their ID, and then IDs in the same group are renamed. r sort : the sequential sort algorithm implemented in seqkit is replaced by a sample MergeSort [ 24 ] algorithm that can be efficiently executed in parallel in a distributed environment. r shuffle : sequences shuffling can be implemented using the IgnisHPC API function partitionByRandom . Another implementation details In order to parallelize and integrate the seqkit routines into IgnisHPC, it was necessary to start considering the sequence parser. It takes a stream of c har acters in F ASTA and F ASTQ format and generates a data structure with the sequence r epr esentation. In seqkit , this stream can be r epr esented by a file or the standard input. In BigSeqKit , this stream is implemented using the IgnisHPC iter ators, whic h gr ant the users access to the data partitions. In this way, BigSeqKit will read the data from a file and split it in m ultiple partitions, whic h facilitates their par allel pr ocessing. In particular, each worker reads a portion of the input file, so the input/output (I/O) operation is performed in parallel. There is 1 worker per computing node. Within each worker, its portion of the file is further divided among the available threads , impro ving the ov er all I/O performance. As a result, the seqkit command arguments that affect file processing will have no effect in BigSeqKit . For example, the --two-pass option, which reads a file multiple times instead of storing all the sequences in memory, does not make sense in BigSeqKit . We must highlight that the fact of splitting the input files between se v er al computing nodes in BigSeqKit means that the memory consumed by node is also split, which allows our tool to work with larger datasets. In addition, BigSeqKit also reduces the memory footprint by only storing the IDs and indices of each sequence. Another important adv anta ge of using IgnisHPC is how memory is handled. Users can choose a type of stor a ge according to their particular case. For instance, if an input file is too large to be k e pt completely in the server memory, it could be stored compressed in memory or in disk. Performance would be lower, but it could be successfull y pr ocessed. That scenario is not considered by seqkit that simply would raise an “out of memory” error . In particular , BigSeqKit supports the following stor a ge options: Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 BigSeqKit | 5 r In-Memory: it is the best performer since all data ar e stor ed in memory. It is the default option. r Raw memory: data are stored in a memory buffer using a serialized binary format. Extra memory consumption is minimal and the buffer is compressed by Zlib. r Disk: similar to raw memory but the buffer is stored as a POSIX file. Although the performance is significantly worse, it enables working with vast amounts of data that cannot be entir el y k e pt in memory. On the other hand, rmdup , common , and pair commands in seqkit use hash functions to c hec k duplicates. It is well known that hash functions can produce the same result for different values . T his event is commonly known as a hash collision. Howe v er, seqkit does not c hec k for collisions, so it is possible to gener ate incorr ect r esults. BigSeqKit uses hashes to gr oup sequences but then c hec ks for collisions by comparing the real values. Finally, seqkit and other state-of-the-art tools build index files ( faidx routine) to speed up some other tasks (e.g., searches). Although BigSeqKit is also capable of creating those index files, it does not r equir e them to impr ov e its performance since data within IgnisHPC ar e alr eady indexed. In other w or ds, the index is created while reading the input file. New functionalities BigSeqKit not only enables the parallelization of seqkit functions but also impr ov es its algorithms to provide benefits e v en for sequential executions and includes additional functionalities. In particular, the faidx command in seqkit implements indexing of FASTA files using the samtools format, but FASTQ files are not supported. BigSeqKit adds support for this type of files and generates an index file using the samtools format as well. Note that this is the most widespread format and is also supported by other state-of-the-art tools . T her efor e, BigSeqKit allows indexing of both FASTA and FASTQ files using the same syntax than seqkit . How to Use BigSeqKit BigSeqKit can be used in 2 different wa ys . T he first one is by means of a command-line interface (CLI). This a ppr oac h is similar to the “command subcommand” structure adopted by seqkit [ 13 ]. In this way, it is only necessary to select a subcommand or routine (see a complete list in Table 2 ) and pass its ar guments thr ough command line. As we mentioned pr e viousl y, to impr ov e the usability and facilitate the adoption of BigSeqKit , it implements the same command interface as seqkit . Since BigSeqKit runs within the IgnisHPC fr ame work, it is necessary to send the BigSeqKit routine through the IgnisHPC submitter. For instance, if we are running BigSeqKit on a local server, the following expression uses the routine seq to print the name of the sequences included in a FASTA file to an output file: Ther efor e, the syntax should be: ignis-submit ignishpc/full bigseqkit < cmd > < arguments > . In addition, users can also specify through arguments the number of instances , cores , and memory (in GB) to be used in the execution. By default, those values are set to 1. For example, we can execute the pr e vious command using 2 cores: Figure 3: Example of Python code using the BigSeqKit routines. Unlike the other state-of-the-art tools, BigSeqKit can also be executed on a parallel cluster. Typical HPC clusters have Slurm [ 20 ] as the pr eferr ed r esource mana ger and Singularity [ 25 ] as a container-based technology. In this case, users will send the job using the ignis-slurm submitter instead of ignissubmit . On the other hand, BigSeqKit can also be used as a bioinformatics library. It is worth noting that BigSeqKit was implemented in Go language. Ho w ever, thanks to the multilanguage support provided by IgnisHPC, it is possible to call BigSeqKit routines from C/C ++ , Python, Java, and Go applications without additional overhead. An example of Python code is shown in Fig. 3 . This example is equivalent to the previous one used in the explanation of the CLI. Since BigSeqKit has been created as a library, it only needs to be imported to be used. Functions in BigSeqKit do not use files as input; they use DataFrames instead, an abstract representation of parallel data used by IgnisHPC (similar to RDDs in Spark). P ar ameters ar e gr ouped in a data structur e wher e eac h field r epr esents the long names of a parameter. We must highlight that BigSeqKit functions can be link ed (lik e system pipes using “| ”), so the DataFrame generated by one can be used as input to another. In this way, integrating BigSeqKit routines in a more complex code is really easy. The code starts initializing the IgnisHPC fr ame work (line 5 in Fig. 3 ). Next, a cluster of containers is configured and built (lines from 7 to 15). Multiple parameters can be used to configure the envir onment suc h as ima ge , number of containers , number of cores , and memory per container. In this example, we will use 1 node (instances) and 2 cores by node. After configuring the IgnisHPC execution environment, the BigSeqKit code actually starts. First, we read the input file (line 17). There is a different function for reading F ASTA and F ASTQ files. All the input sequences are stored as a single data structure . T he next stage consists of printing the name of the sequences included in the FASTA file (line 19). The function Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 6 | GigaScience , 2023, Vol. 12, No. 1 takes as parameters the sequences and the options that specify its behavior. Finally, the names of the sequences are written to disk. It is important to highlight that lazy e v aluation is performed, so functions ar e onl y executed when the result is r equir ed to be saved on disk. Experimental Results In this section, we analyze the performance results obtained by BigSeqKit with respect to other state-of-the-art tools. In particular, we hav e consider ed samtools , pyf astx , and seqkit for their performance and number of commands supported. Experiments were conducted using up to 8 computing nodes of the FinisT er - rae III [ 26 ] supercomputer installed at CESGA (Spain). Each node contains a 32-core Intel Xeon Ice Lake 8352Y @2.2 GHz processor and 256 GB of memory interconnected with Infiniband HDR 100. It is a Linux cluster running Rocky Linux v8.4 (kernel v4.18.0). We have used SingularityCE v3.9.7 (containers), IgnisHPC v2.2, pyfastx v0.8.4, samtools v1.16.1, and seqkit v2.3.1 (with Slurm as cluster manager and Lustre as distributed file system). The performance e v aluation was carried out using as input 6 different F ASTA/F ASTQ files that cov er a wide v ariety of c har acteristics and sizes . T he main features of these files are the following: r D 1 ( m64013e_210227_222017.hifi_reads —FASTA— 24 GB): Number of sequences: 1.2M, Minimum length: 85, Average length: 19.7K, Maximum length: 48.5K. r D 2 ( SRR642648_1.filt —FASTQ—24.1 GB): Number of sequences: 98.7M, Minimum length: 100, Av er a ge length: 100, Maximum length: 100. r D 3 ( Homo_sapiens.GRCh38.dna_sm.toplevel —FASTA— 59.7 GB): Number of sequences: 639, Minimum length: 970, Av er a ge length: 98.8M, Maximum length: 248.9M. r D 4 ( ERR4667750 —FASTQ—79.1 GB): Number of sequences: 318.1M, Minimum length: 101, Av er a ge length: 101, Maximum length: 101. r D 5 ( uniprot_trembl —FASTA—104 GB): Number of sequences: 229.9M, Minimum length: 7, Av er a ge length: 351.6, Maximum length: 45.3K. r D 6 ( DRR002180_2 —FASTQ—395 GB): Number of sequences: 1.625B, Minimum length: 101, Av er a ge length: 101, Maximum length: 101. As example to illustrate the benefits of our tool, we will e v aluate the following utilities (see Table 2 for a complete list of commands): faidx builds an index for F ASTA/F ASTQ files, locate locates sequences following some search pattern, replace replaces a name/sequence using a r egular expr ession, rmdup r emo ves duplicated sequences , sample selects sequences by number or proportion, seq transforms sequences (extract ID, filter by length, etc.) and r emov es ga ps, and sort sorts sequences by ID/name/sequence/length. We will also include the performance results of the corresponding utilities, if they exist, for samtools , pyfastx , and seqkit . Execution times for all the tools considered include the overhead of loading sequences into memory and the subsequent writing of results to disk. Note that the “2-pass” argument of seqkit was not used in the experiments. Eac h r esult was computed as the median of 5 experiments. For the sake of r epr oducibility, all the codes and scripts used for performing the benchmarks are freely available at the BigSeqKit repository. First, in order to provide an ov er all idea about the scalability and performance of BigSeqKit with respect to the other stateof-the-art tools, we will show the speedups obtained for the D 4 dataset using different number of cores . T he beha vior is very similar when considering the other datasets. Results in log scale are displayed in Fig. 4 . Speedups were calculated using as reference the sequential execution (1 core) of the corresponding BigSeqKit command. According to the r esults, se v er al conclusions can be made. It can be observed that the scalability of BigSeqKit is quite good, r eac hing speedups up to 27.7 ×and 95.7 ×( seq command) using 1 server (32 cores) and 8 computing nodes (256 cores), respectiv el y. Note that speedups of some routines are not higher when using 256 cores due to a small fraction of the code that should be executed sequentially (Amdahl’s law). While samtools and pyfastx routines are always processed sequentially, seqkit uses a multithreaded approach to (partly) parallelize some commands. Ho w e v er, its scalability is limited to use a few threads on a single server (computing node). This is the case of locate . Its best speedup onl y r eac hes 11.3 ×(32 cores) while this value increases until 19.6 ×with BigSeqKit . If 8 nodes are used, BigSeqKit is 49.9 ×faster than the sequential execution. For all the commands studied, BigSeqKit clearly outperforms samtools , pyfastx , and seqkit . Ther e ar e onl y a fe w cases using 1 core where the speedups of these tools are slightly greater than 1 (e.g., executing the faidx routine with samtools and pyfastx ). Ho w ever, other commands such as sort and sample are processed faster with BigSeqKit e v en using 1 core. Tables from 3 to 9 display, for all the datasets, the execution times of BigSeqKit and the other state-of-the-art tools when running faidx , locate , replace , rmdup , sample , seq , and sort utilities, r espectiv el y. Speedups with respect to the sequential execution of the corresponding BigSeqKit command are shown between br ac kets . Highlighted is the fastest time o v er all and the corresponding speedup. Note that BigSeqKit stores compressed in memory the largest dataset D 6 when using 1 computing node since it exceeds the memory capacity of an individual server (see the Raw memory storage option in the Background section). For all the experiments conducted, BigSeqKit is always the fastest tool both considering a single server (1 node) or se v er al computing nodes. In any case, let’s take a look in detail of the behavior for each command: r faidx (Table 3 ): BigSeqKit speedups range from 5.4 ×to 27.4 × considering a single server (32 cores) and from 7.2 ×to 144 × with 8 nodes. It means, for example, building the index file for our largest dataset D 6 (395 GB) in just 5.8 minutes (single server), while samtools and pyfastx require about 2.1 hours. This time decreases to 1 minute when BigSeqKit uses 8 nodes. As mentioned pr e viousl y, the faidx r outine in seqkit does not support FASTQ files (D 2 , D 4 , and D 6 ). r locate (Table 4 ): the searching routines, grep and locate , ar e v ery expensiv e in terms of computations. Note that considering sequential processing, locate takes more than 3 hours to process our smallest dataset D 1 independently of the tool considered. This time increases to more than 3 days of computation for D 6 . seqkit has a m ultithr ead v ersion of locate , which obtains speedups from 10.5 ×to 18.8 ×. These speedups are always lower than the ones obtained by BigSeqKit on a single server. It is important to highlight that seqkit raises an out-of-memory error when processing D 6 with 1, 2, and 4 cores. On the other hand, when using 8 nodes, BigSeDownloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 BigSeqKit | 7 Figure 4: Speedups (in log scale) obtained by BigSeqKit and other state-of-the-art tools with respect to the BigSeqKit sequential time when executing different commands using D 4 as input. Note that locate was parallelized in seqkit . qKit ac hie v es noticeable speedups up to 104.1 ×. In this way, it is able to reduce the time necessary to execute the locate command with our largest dataset D 6 from 3 days to only 0.8 hours. r replace (Table 5 ): this routine (or an equivalent) is not supported by samtools and pyfastx . In this case, BigSeqKit is from tens to hundreds of times faster than seqkit , r eac hing speedups up to 159.8 ×. r rmdup (Table 6 ): this routine is also not supported by samtools and pyfastx . In this case, BigSeqKit is tens of times faster than seqkit , ac hie ving a maxim um speedup of 74.7 ×when removing the duplicated sequences in D 5 . r sample (Table 7 ): operation not supported by samtools . BigSeqKit is again faster than the other tools, increasing the speedups as the input data size grows. It can be observed that BigSeqKit is able to sample sequences in seconds. For instance , p yfastx and seqkit take about 3 hours to process D 6 , while BigSeqKit r equir es just 2 minutes. r seq (Table 8 ): operation not supported by samtools . Performance results are similar to the sample ones in such a way that BigSeqKit filters sequences by ID in a few seconds, achieving a noticeable speedup of 169.7 ×. It should be noted that among the routines examined in this study, seq is the least computationally demanding. r sort (Table 9 ): this routine was not included in pyfastx . In general, the performance of samtools and seqkit is poor. And, most importantly, both tools produce memory errors when processing the largest dataset D 6 , so it cannot be sorted. Howe v er, BigSeqKit sorts D 6 21.8 ×and 131.1 ×faster than the sequential version using a single server and 8 computing nodes, r espectiv el y. It means that the time decr eases fr om 5 hours to bar el y 2 minutes. Finall y, we m ust highlight that one of the main reasons for the differences in the speedups between datasets running the same command with BigSeqKit is the load balance between threads. It will depend on the c har acteristics of the dataset: number of sequences and their length. Conclusions Current state-of-the-art tools such as seqkit , pyfastx , and samtools ar e not r eady for pr ocessing and manipulating v ery lar ge FASTA and FASTQ files because all of them ar e mainl y based on sequential processing. To that end, we have presented BigSeqKit , which parallelizes and optimizes the seqkit routines using the IgnisHPC computing fr ame w ork. Since seqkit w as pr ogr ammed in Go, IgnisHPC was extended to support that language. As a consequence, IgnisHPC is no w adays the first par allel computing fr ame work that supports Go. BigSeqKit can be easily installed on a local server or on a cluster. In addition, it can be used from the command line or as a library. Thanks to the multilanguage support of IgnisHPC, BigSeqKit routines can be called from C/C ++ , Python, Java, and Go codes. Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 8 | GigaScience , 2023, Vol. 12, No. 1 Table 3: Execution times (seconds) using different number of cores: faidx command. Highlighted are fastest time and number of times faster than sequential BigSeqKit 1 2 4 8 16 32 64 (2 nodes) 128 (4 nodes) 256 (8 nodes) D 1 samtools 86.2 [1.03 ×] – – – – – – – – pyfastx 109.2 [0.81 ×] – – – – – – – – seqkit 75.4 [1.17 ×] – – – – – – – – BigSeqKit 88.4 46.0 35.3 26.3 19.4 16.3 [5.4 ×] 13.6 12.3 [7.2 ×] 12.5 D 2 samtools 165.6 [1.06 ×] – – – – – – – – pyfastx 177.9 [0.99 ×] – – – – – – – – BigSeqKit 175.9 90.8 67.4 50.3 39.1 31.4 [5.6 ×] 23.4 19.1 15.5 [11.3 ×] D 3 samtools 210.0 [0.77 ×] – – – – – – – – pyfastx 131.2 [1.23 ×] – – – – – – – – seqkit 131.8 [1.23 ×] – – – – – – – – BigSeqKit 161.9 83.9 61.7 24.5 17.5 15.7 [10.3 ×] 13.6 13.4 [12.1 ×] 14.7 D 4 samtools 538.4 [1.27 ×] – – – – – – – – pyfastx 615.5 [1.11 ×] – – – – – – – – BigSeqKit 684.2 346.6 175.2 90.3 45.4 29.3 [23.3 ×] 19.6 15.3 12.5 [54.7 ×] D 5 samtools 771.0 [1.08 ×] – – – – – – – – pyfastx 634.3 [1.31 ×] – – – – – – – – seqkit 1,096.2 [0.76 ×] – – – – – – – – BigSeqKit 829.8 361.3 179.4 89.3 49.4 30.3 [27.4 ×] 23.6 19.3 16.5 [50.3 ×] D 6 samtools 7,651.6 [1.14 ×] – – – – – – – – pyfastx 7,712.5 [1.13 ×] – – – – – – – – BigSeqKit 8,712.3 4,423.3 2,282.2 1,191.9 640.2 350.4 [24.9 ×] 129.5 85.3 60.5 [144 ×] Table 4: Execution times (seconds) using different number of cores: locate command. Highlighted are fastest time and number of times faster than sequential BigSeqKit 1 2 4 8 16 32 64 (2 nodes) 128 (4 nodes) 256 (8 nodes) D 1 pyfastx 11,523.5 [1.0 ×] – – – – – – – – seqkit 12,822.9 6,385.0 3,210.9 1,731.4 940.5 612.4 [18.8 ×] – – – BigSeqKit 11,486.2 6,286.1 3,180.0 1,637.3 850.9 470.6 [24.4 ×] 264.6 156.9 110.3 [104.1 ×] D 2 pyfastx 8,841.2 [1.2 ×] – – – – – – – – seqkit 12,319.8 6,909.4 3,335.9 1,746.2 997.3 971.2 [10.5 ×] – – – BigSeqKit 10,168.6 5,264.5 2,711.5 1,412.2 814.6 545.4 [18.6 ×] 384.7 293.5 234.9 [43.3 ×] D 3 pyfastx 13,075.3 [1.1 ×] – – – – – – – – seqkit 14,281.6 8,161.7 5,009.6 3,184.1 1,832.4 1,054.9 [14.1 ×] – – – BigSeqKit 14,834.2 8,223.3 4,572.8 2,585.6 1,494.6 872.1 [17.0 ×] 532.8 365.9 262.5 [56.5 ×] D 4 pyfastx 30,028.3 [1.05 ×] – – – – – – – – seqkit 39,640.5 21,257.6 10,803.1 5,715.1 3,369.7 2,795.2 [11.3 ×] – – – BigSeqKit 31,615.2 16,832.1 8,531.9 4,433.3 2,466.8 1,609.9 [19.6 ×] 1,074.7 794.6 633.5 [49.9 ×] D 5 pyfastx 27,876.5 [1.06 ×] – – – – – – – – seqkit 31,301.8 16,884.7 9,141.1 4,698.4 2,971.8 2,802.9 [10.5 ×] – – – BigSeqKit 29,540.7 15,431.3 8,120.2 4,401.4 2,454.5 1,443.9 [20.5 ×] 908.1 599.5 440.9 [67 ×] D 6 pyfastx 270,214 [1.02 ×] – – – – – – – – seqkit Out of Mem. Out of Mem. Out of Mem. 40,122 23,075 18,309 [15.0 ×] – – – BigSeqKit 275,680 141,095 72,110 37,140 19,810 11,477 [24.0 ×] 7,003 4,422 3,080 [89.5 ×] Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023 BigSeqKit | 9 Table 5: Execution times (seconds) using different number of cores: replace command. Highlighted are fastest time and number of times faster than sequential BigSeqKit 1 2 4 8 16 32 64 (2 nodes) 128 (4 nodes) 256 (8 nodes) D 1 seqkit 132.4 [1.02 ×] – – – – – – – – BigSeqKit 134.5 69.5 36.1 25.0 18.7 12.7 [10.6 ×] 13.1 13.6 12.5 [10.8 ×] D 2 seqkit 395.7 [1.04 ×] – – – – – – – – BigSeqKit 410.6 213.5 110.1 74.5 56.9 29.7 [13.8 ×] 16.8 13.9 13.5 [30.4 ×] D 3 seqkit 410.5 [0.99 ×] – – – – – – – – BigSeqKit 406.7 209.5 109.4 74.0 56.1 29.5 [13.8 ×] 15.3 13.6 12.9 [31.5 ×] D 4 seqkit 543.7 [1.05 ×] – – – – – – – – BigSeqKit 570.3 293.5 109.4 74.0 55.1 29.4 [19.4 ×] 20.3 13.5 12.5 [45.6 ×] D 5 seqkit 1,572.1 [1.03 ×] – – – – – – – – BigSeqKit 1,621.7 819.9 420.1 217.2 115.1 62.9 [25.8 ×] 37.2 24.2 18.5 [87.7 ×] D 6 seqkit 8,980.8 [1.07 ×] – – – – – – – – BigSeqKit 9,620.8 5,000.3 2,605.2 1,364.2 717.7 387.5 [24.8 ×] 142.1 90.5 60.2 [159.8 ×] Table 6: Execution times (seconds) using different number of cores: rmdup command. Highlighted are fastest time and number of times faster than sequential BigSeqKit 1 2 4 8 16 32 64 (2 nodes) 128 (4 nodes) 256 (8 nodes) D 1 seqkit 178.9 [1.01 ×] – – – – – – – – BigSeqKit 180.5 94.3 50.2 35.1 27.1 15.8 [11.4 ×] 14.8 14.4 13.8 [13.1 ×] D 2 seqkit 320.6 [1.04 ×] – – – – – – – – BigSeqKit 333.3 174.7 93.5 65.9 49.9 26.5 [12.6 ×] 15.9 14.1 [23.6 ×] 15.0 D 3 seqkit 515.5 [0.91 ×] – – – – – – – – BigSeqKit 469.5 246.7 182.7 127.5 96.1 51.4 [9.1 ×] 27.4 20.9 20.6 [22.8 ×] D 4 seqkit 729.9 [0.99 ×] – – – – – – – – BigSeqKit 720.5 378.5 197.5 139.7 102.9 54.0 [13.3 ×] 30.5 16.4 14.1 [51.1 ×] D 5 seqkit 2,173.6 [0.97 ×] – – – – – – – – BigSeqKit 2,100.2 1,110.4 612.3 341.2 195.1 115.2 [18.2 ×] 70.5 43.2 28.1 [74.7 ×] D 6 seqkit 9,937.1 [1.11 ×] – – – – – – – – BigSeqKit 11,022.3 5,578.5 3,006.7 1,709.6 1,004.1 600.1 [18.4 ×] 275.2 241.6 228.8 [48.2 ×] Regarding the experimental results, BigSeqKit clearly outperforms seqkit , pyfastx , and samtools for all the tasks considered. On a single server, BigSeqKit is overall tens of times faster than those state-of-the-art tools, r eac hing speedups with respect to the BigSeqKit sequential time up to 27.7 ×. Considering an 8-node cluster, BigSeqKit is e v en faster, r eac hing speedups higher than 160 ×. It means that most of the tasks can be performed in just a few seconds . For instance , our toolkit effectiv el y r educes the execution time of the locate command on our largest dataset from 3 days to a mere 0.8 hours. It is important to highlight that seqkit and samtools were unable to process that dataset with some routines due to memory issues, which confirms that current state-of-theart tools are not well fitted for processing very large files. As future w ork, w e plan to add also the remainder seqkit commands not included in the current version of BigSeqKit : sliding , sana , fx2tab , tab2fx , convert , amplicon , fish , split , split2 , restart , and mutate . Note that all of them are independent routines, so their implementation using IgnisHPC will be straightforw ar d. Availability of Source Code and Requirements Project name: BigSeqKit r Pr oject homepa ge: https://github.com/citiususc/BigSeqKit Downloaded from https://academic.oup.com/gigascience/article/doi/10.1093/gigascience/giad062/7233988 by UDC - Universidade da Coruna user on 26 October 2023