Compression of next-generation sequencing quality scores using memetic algorithm

Abstract

Background

The exponential growth of next-generation sequencing (NGS) derived DNA data poses great challenges to data storage and transmission. Although many compression algorithms have been proposed for DNA reads in NGS data, few methods are designed specifically to handle the quality scores.

Results

In this paper we present a memetic algorithm (MA) based NGS quality score data compressor, namely MMQSC. The algorithm extracts raw quality score sequences from FASTQ formatted files, and designs compression codebook using MA based multimodal optimization. The input data is then compressed in a substitutional manner. Experimental results on five representative NGS data sets show that MMQSC obtains higher compression ratio than the other state-of-the-art methods. Particularly, MMQSC is a lossless reference-free compression algorithm, yet obtains an average compression ratio of 22.82% on the experimental data sets.

Conclusions

The proposed MMQSC compresses NGS quality score data effectively. It can be utilized to improve the overall compression ratio on FASTQ formatted files.

Background

DNA sequencing provides fundamental data for many research areas e.g. genomics, bioinformatics, and biology [1]. Rapid progress has been made for DNA sequencing technologies, especially the high-throughput next-generation sequencing (NGS), in the last few years. Newly proposed high efficiency methods significantly stimulate the production and usage of NGS data [2]. However the exponential growth of NGS data poses huge challenge to data storage and transmission [3]. Thereby efficient compression algorithms are required.

General-purpose compression algorithms e.g. gzip and bzip2 usually fail to obtain satisfactory results on NGS data, because they are designed for ordinary plain text or binary files. To achieve higher compression ratio, many specialized methods are proposed. For instance, Cox et al. [4] proposed a Burrws-Wheeler transform based compression algorithm for large scale DNA sequence databases. Jones et al. [5] presented Quip, a high efficient reference-based NGS data compression tool relying on external reference genomes or light-weight de novo assembly to generate reference sequences from the target data. Popitsch et al. [6] proposed the NGC tool for SAM format files compression. Hach et al. [7] proposed SCALCE by introducing locally consistent parsing in data encoding. More comprehensive review on NGS data compression can be found in [8, 9].

Typically, raw NGS data includes a series of sequencing records (called reads). Each record consists of three major components: a metadata containing the read name, platform, and other useful information; a DNA sequence read obtained from one fold of the oversampling; and a quality score sequence estimating accuracies of the corresponding DNA bases. Existing algorithms usually focus on the compression of DNA sequence reads, and utilize conventional methods e.g. Huffman coding and run-length encoding (RLE) to compress quality scores [10]. Quality score data is considered more important than the metadata, and usually occupies similar or even more space than the DNA sequences. It may pose bigger challenges for compression than DNA sequence reads, due to the larger alphabet size. By introducing compressor specific for NGS quality scores, the overall compression ratio of NGS data can be further improved.

In this paper, we propose a memetic algorithm (MA) based NGS quality scores compression algorithm, namely MMQSC. The algorithm consists of three major parts: a Huffman coding based preprocessing is conducted in the first place, followed by MA based encoding codebook design. Finally, quality score data is compressed by using the codebook. MA is widely known as a synergy of population-based evolutionary algorithm and local search or individual learning procedures. MAs are capable of solving various complex optimization problems more effectively than their conventional counterparts [11]. In this work, the self-adaptive differential evolution combining with neighborhood search (SaNSDE) [12], and Davies, Swann, and Campey with Gram-Schmidt orthogonalization (DSCG) [13], or SaNSDE-DSCG for short, are introduced to MMQSC, to optimize the NGS quality scores compression codebook, with which most repetitive short score segments are identified and represented with much shorter encoding.

In conventional MAs, each individual represents a candidate solution of the entire problem, i.e., compression codebook. Its optimization is highly complex, because the codebook consists of hundreds of quality score symbols. Multimodal optimization tries to find all or most of the multiple solutions within a population in a single simulation run [14]. Based on multimodal optimization, the MMQSC uses an individual to represent only a single specific code vector, and composes codebook with the entire evolution population. Thereby computational complexity distributed to each individual's fitness evaluation is significantly reduced.

The proposed MMQSC obtains promising performance on NGS quality scores stored in the widely used FASTQ format [15]. Experimental results on five representative NGS data sets show that MMQSC obtains better compression ratio than other state-of-the-art methods. Particularly, MMQSC is a reference-free algorithm for lossless compression, yet obtains an average compression ratio of 22.82%, i.e., 1.81 bits per quality value (BPQ) on the experimental data sets.

The remainder of this paper is organized as follows: Section II describes details of the proposed MMQSC compression algorithm. Section III presents the experimental results of MMQSC on the five real-world NGS data sets. Finally, a conclusion is provided in Section IV.

Methods

SaNSDE and DSCG based memetic algorithm for multimodal optimization

MA is introduced whereby the concept of "meme", which was coined by Dawkins [16], is employed within an evolutionary computation framework to improve search efficiency. Typically, MA utilizes a population-based global optimization as fundamental framework, and introduces separate local searches or 'memes' embedded in each generation of the global evolution to refine the population [17]. The procedure of a canonical MA framework is illustrated in Algorithm 1 [18].

In MAs both global search and local search strategies can be selected flexibly according to the target problem. Typically, NGS data consists of thousands or even millions of read entries, wherein each of them contains hundreds of quality score symbols. Finding a codebook for compressing such data is naturally a high-dimensional optimization problem. Differential evolution (DE) [19] is capable of solving large scale problems effectively. In this paper, a high performance DE variant namely the SaNSDE is utilized as the global optimizer. Particularly, SaNSDE uses three self-adaptive mechanisms to select mutation strategies and control parameter values. By introducing neighborhood search in the optimization process, SaNSDE obtains higher performance than conventional algorithms. Moreover, the widely-used local search strategy DSCG is introduced to increase convergence speed. DSCG is a gradient-based optimizer that searches solution space in a greedy manner. Combining the exploration of SaNSDE and exploitation of DSCG, the proposed MA obtains promising performance on quality scores compression codebook design.

As shown in Figure 1, multimodal optimization searches for not only the global best solution gbest (as single-objective optimization), but also all the local optimal pbest _i , i = 1, 2, 3... Multimodal optimization has been used in a wide range of applications, because it can locate all or most of the optimal solutions in a single simulation run. Fitness sharing [20] is introduced in the proposed SaNSDE-DSCG to conduct multimodal optimization. Given a raw fitness value f _R ( x _i ), wherein x _i is the candidate solution of individual ps _i . Fitness sharing transforms it into shared fitness f _S ( x _i ) using following equation:

Figure 1

Multimodal optimization.

$$f_S\left({\mathit{x}}_i\right)=f_R\left({\mathit{x}}_i\right)\times {\tau }_i$$

(1)

in which:

$${\tau }_i=\sum _{j\in \left|ps\right|,j\ne i}{\left(1-\frac{d_{i,j}}{\epsilon }\right)}^{\alpha }$$

(2)

where ε is the niching radius, parameter α controls the shape of the sharing function, distance d _i,j is defined as:

$$d_{i,j}=\left\{\begin{array}{c}\hfill dist\left({\mathit{x}}_i,{\mathit{x}}_j\right)ifdist\left({\mathit{x}}_i,{\mathit{x}}_j\right)\le \epsilon \hfill \\ \hfill 0otherwise\hfill \end{array}\right.$$

(3)

where dist( x _i , x _j ) is the Manhattan distance between x _i and x _j . If evolution population is gathering in the same optimal region, its shared fitness values will deteriorate significantly to disperse the individuals. By utilizing f _S ( x _i ) to guide the search process, SaNSDE-DSCG is capable of finding all optimums effectively.

Compression codebook design using SaNSDE-DSCG

As shown in Figure 2, a NGS quality score sequence is compressed by substituting original scores with the index of its most similar code vector in the codebook, and their corresponding symbol differences.

Figure 2

Codebook based NGS quality score sequences compression.

Given a quality score sequence Q= "CCCGFF' and code vector C= "CCGHFFC". Sequence Q is encoded as {i, Q^∗}, where i is in the index of C, and Q^∗ records the symbol differences as:

$$\begin{array}{c}\hfill \mathit{Q}=\mathsf{\text{C}}\mathsf{\text{C}}\mathsf{\text{C}}\mathsf{\text{G}}-\mathsf{\text{F}}\mathsf{\text{F}}-\hfill \\ \hfill \mathit{C}=\mathsf{\text{C}}\mathsf{\text{C}}\wedge \mathsf{\text{G}}\mathsf{\text{H}}\mathsf{\text{F}}\mathsf{\text{F}}\mathsf{\text{C}}\hfill \\ \hfill {\mathit{Q}}^{*}=UU\left(I,\mathsf{\text{"C"}}\right)UDUUD\hfill \end{array}$$

(4)

in which U denotes symbol matching (unchanged), I stands for insertion (marked with "∧"), D for deletion (marked with "−"), and S for substitution. For insertions and substitutions, the original symbol should also be recorded (for instance "C" on the third quality score). This matching process is conducted by using dynamic programming (DP).

Data size of the original P-dimensional sequence is L _O = 8 × P, because raw quality scores are stored in 8 bits ASCII format. On the other hand, each difference type in {U, I, D, S} takes 2 bits to represent. Thereby the encoded data size is:

$$L_C=\lceil lo{g}_2\left(M\right)\rceil +2\times P^{*}+8\times R$$

(5)

where M is the number of code vectors in a codebook, P^∗ denotes length of the symbol differences sequence (i.e. Q^∗), and R is the number of all original symbols recorded for insertions and substitutions. Given the example above, original quality score sequence takes L _O = 48 bits of storage, while the encoded data uses only about 24 bits. Usually the encoding process makes L _C smaller than L _O , i.e., conducts compression. The more the code vector is similar to the original quality score sequence, the higher compression ratio is achieved. That is, quality of the codebook decides the overall compression performance.

In this paper, we utilize the proposed SaNSDE-DSCG to optimize compression codebook design. During the initialization, code vectors in the codebook are generated randomly, and encoded as individuals to form an evolution population. In each fitness evaluation, input candidate solution x _i = [x _i,1 , x _i,2 ,..., x _i,N ] is mapped into code vector C _i = "s₁ s₂ ... s _N " using the following equation:

$$s_n=\mathit{S}[\left.⌈x_n\rceil \right],n=1,2,\dots ,N$$

(6)

where characters set S includes all unique symbols in the original quality score sequences, variable N is the predefined code vector length. This mapping is conducted because candidate solutions consist of continuous values, while code vectors are discrete symbol strings. The C _i is then matched to all quality score sequences, and calculates the corresponding encoded data size. Raw fitness value of x _i is defined as:

$$f_R\left({\mathit{x}}_i\right)=\sum _{k\in K}{L}_C\left({\mathit{C}}_i,{\mathit{Q}}_k\right)$$

(7)

in which K is the total number of quality score sequences, L _C ( C _i , Q _k ) denotes encoded data size on the k th sequence Q _k using code vector C _i . Small f _R ( x _i ) indicates that C _i is more similar to the original score data, i.e., obtains higher compression ratio. Shared fitness f _S ( x _i ) is then calculated accordingly.

It is noted that accurate symbol differences information, e.g. mismatched symbol positions, is not necessary for fitness values calculation. In most cases the approximate Levenshtein distance [21] is good enough to guild the codebook optimization process. Moreover, calculation of Levenshtein distance requires much less computational resources than the DP based matching algorithm. By utilizing Levenshtein distance in fitness evaluation, we can achieve similar optimization performance in a relatively high speed. Approximate size of encoded data can be calculated as:

$$L_C^{\prime }=\lceil lo{g}_2\left(M\right)\rceil +2\times P_k+4\times lev\left({\mathit{C}}_i,{\mathit{Q}}_k\right)$$

(8)

in which P _k is the length of Q _k , and lev(·) denotes Levenshtein distance between the two input sequences. It's value is multiplied by 4, because there is a half chance (I and S) in {U, I, D, S} needs to record the original quality score. Accurate matching information is obtained only after the codebook design process for actual sequences compression.

Procedure of the SaNSDE-DSCG based compression codebook design algorithm is illustrated in Figure 3. In conventional single-objective optimization based design algorithms, the entire compression codebook is encoded in each individual in the evolution population. Typically, an individual in such methods is constructed by connecting all code vectors in the codebook from end to end, i.e. x _i ^' = {C₁, C₂,..., C _M } [22]. Optimal codebook is obtained by searching the global best solution. Dimensionality of solution space is M × N, and the algorithm calculates encoded data size (Equ. (8)) for M × N × |ps| time in each generation of the evolution optimization. Its computational complexity is too high to be applied on large NGS data. In contrary, MMQSC searches the solution space by using multimodal optimization, wherein each individual is utilized to represent one specific code vector, and the entire evolution population is utilized to compose the optimal compression codebook. MMQSC evolves each individual to make the code vector more representative to original quality score sequences, and accordingly the optimal codebook as a whole can compress input data more effectively. Moreover, dimensionality of solution space in MMQSC is reduced to N, because individual x _i and corresponding code vector C _i have the same length. In each generation of the evolution optimization, the encoded data size is calculated for only M × N times.

Figure 3

Compression codebook design using SaNSDE-DSCG based multimodal optimization.

The MMQSC algorithm

The proposed MMQSC algorithm consists of three major parts: Huffman coding based preprocessing, SaNSDE-DSCG based codebook design, and quality score data compression. Details of SaNSDE-DSCG optimization algorithm and its application on compression codebook design have been discussed in the previous sections.

During the preprocessing, raw quality score sequences are extracted from target FASTQ file, and undergo a Huffman coding. Each sequence is then converted into a symbol string by mapping every 6 bits of the encoded binary data to readable ASCII character:

$$h_t=chr\left(33+int\left(huff\left({\mathit{Q}}_k\right)\left[t\times 6~\left(t+1\right)\times 6\right]\right)\right)$$

(9)

in which h _t is the t th symbol in converted string H _k = "h₁ h₂ ... h _T ", function huff(·) denotes Huffman coding, int(·) converts binary data to corresponding integer value, and chr(·) maps integer number into ASCII character. Offset value 33 is added to the input number in chr(·) to ensure exported character is readable. In the majority of cases H _k consists of fewer symbols than original sequence Q _k . Thereby dimensionality of code vectors is reduced, and the codebook design problem is simplified.

The SaNSDE-DSCG is conducted afterward on the encoded sequences. After optimization, MMQSC maps individuals in the evolution population into code vectors using Equ. (6) to construct a compression codebook. This codebook is then utilized to compress the input data, wherein accurate symbol differences information is obtained by using DP based matching algorithm.

Procedure of the MMQSC algorithm is demonstrated in Algorithm 2. It is worth noting that the codebook design process can also be conducted in an offline manner. That is, a universal compression codebook obtained from representative NGS quality score data sets is utilized to encode all input sequences. The time-consuming optimization process is performed for only one time. Successive compressions, which are usually conducted repeatedly, require much less computational time.

Results and discussion

Five representative NGS data obtained from various species, and also of different read numbers and file sizes, are selected to evaluate the overall performance of MMQSC. Details of the data sets are summarized in Table 1. All data are downloaded in FASTQ format from the National Center for Biotechnology Information - Sequence Read Archive (NCBI-SRA) database [23].

Table 1 NGS data sets for MMQSC performance evaluation.

In SaNSDE-DSCG optimization, the compression codebook size M is used as the number of individuals, i.e., |ps|. The value is decided as:

$$M={\left(lo{g}_2\left(\lceil \frac{K}{10}\rceil \right)\right)}^2$$

(10)

The length of code vectors, i.e. dimensionality of solution space, is calculated using the following equation:

$$N=\frac{1}{2}\times \left({\mathsf{\text{min}}}_{k\in K}{P}_k+{\mathsf{\text{max}}}_{k\in K}{P}_k\right)$$

(11)

Parameters setting for SaNSDE-DSCG based multimodal optimization is listed in Table 2 in which |S| is the number of unique symbols in the original quality scores, and FEs denotes the maximum fitness evaluation calls of the optimization.

Table 2 Parameters setting for SaNSDE-DSCG optimization.

Five widely used compression algorithms including the RLE, Huffman coding, gzip, bzip2, and Lempel-Ziv-Markov chain algorithm (LZMA) are utilized for comparison with the proposed algorithm. All algorithms are compared in terms of compression ratios (CR) and bits per quality value (BPQ). The BPQ is defined as follows:

$$BPQ=\frac{{{\sum }_{k\in K}\mathsf{\text{min}}}_{i\in M}{L}_C\left({\mathit{C}}_i,{\mathit{Q}}_k\right)}{{\sum }_{k\in K}{P}_k}$$

(12)

Compression results of all algorithms on the NGS data sets are summarized in Table 3.

Table 3 Compression performance on experimental NGS data sets.

Results in Table 3 show that, the proposed MMQSC obtains better performance than the counterpart representative algorithms. Particularly, MMQSC obtains average compression ratio of 22.82%, resulting in an over 77.18% size reduction in the quality score data. The average 1.81 BPQ result is much smaller than the original 8 BPQ in ASCII format. Moreover, performance of MMQSC remains stable in the experimental data sets, indicating that the algorithm has good robustness on different types of NGS data.

Convergence trace of codebook optimization processes on experimental data sets is illustrated in Figure 4, in which y-axis, labeled as function value, is the optimal shared fitness value in SaNSDE-DSCG optimization. The figure shows that by combining SaNSDE and DSCG in an MA framework to conduct multimodal search, compression codebook is optimized effectively. Particularly, DSCG increases convergence speed in the early stage of optimization. Premature convergence is successfully prevented by using the high performance SaNSDE algorithm.

Figure 4

Convergence trace of compression codebook optimization on experimental data sets.

Conclusions

This paper presents MMQSC, a MA based NGS quality scores compression algorithm. The MMQSC utilizes Huffman coding to preprocess raw quality score sequences stored in FASTQ format. To obtain higher performance, a SaNSDE and DSCG based MA is proposed to optimize the compression codebook design. The Levenshtein distance is used in fitness evaluations to estimate encoded data size, and improves computation speed. After the codebook optimization, a DP based matching algorithm is conducted to obtain accurate symbol differences information. This information, as well as the optimized codebook, is utilized to compress input quality score data. Experimental results on five NGS data show that the proposed MMQSC obtains higher compression ratio than counterpart state-of-the-art algorithms. Particularly, MMQSC reduces about 77% of the storage space on the experimental data sets.