CN111640467B - DNA sequencing quality fraction lossless compression method based on self-adaptive coding sequence - Google Patents

Abstract

The invention provides a DNA sequencing quality score lossless compression method based on a self-adaptive coding sequence, which mainly solves the problem that the compression ratio is low due to the fact that a prediction model of the existing quality score compression method is not accurate enough. The implementation scheme is as follows: 1) Compressing a block P by two encodings ₁ And P ₂ Extracting mass fraction data and base number data in the FASTQ file; 2) Calculating a first encoded compressed block P ₁ The mean value of the mass fraction of each row in the extracted file is quantized to obtain a row mean value matrix F of Mx 1; 3) Counting context information, base information and line mean information of the coded characters, 4) setting two identifiers C and D, and uniformly quantizing the information counted in step 3) to construct a coding model; 5) Driving an adaptive arithmetic coder by a coding model, and compressing a block P by a first code along the direction with the strongest correlation by a snake-shaped coding sequence ₁ And performing traversal compression. The invention improves the compression efficiency and can be used for storing and transmitting gene data.

Description

A Lossless Compression Method for DNA Sequencing Quality Score Based on Adaptive Coding Order

technical field

The invention belongs to the technical field of data compression, and in particular relates to a DNA sequencing quality score lossless compression method, which can be used for the compression of biological gene sequencing data.

Background technique

Sequencing has gradually become a widely used technology in biological research. Obtaining genetic information of different organisms can help us improve our understanding of the organic world. With the rapid development of the new generation of high-throughput gene sequencing technology NGS, sequencing companies represented by Illumina have continuously introduced new sequencing technologies, which has led to a rapid decline in sequencing costs. The price of human whole genome sequencing WGS has dropped to US$1,000 or even lower. And still decline at a rate higher than Moore's Law. In this case, the amount of next-generation sequencing data generated will exceed that of astronomical data, and the corresponding overhead of storing and transmitting these data will also increase. Therefore, it is of great significance to reduce the size of gene sequencing data through data compression, thereby reducing storage and transmission costs. At present, many achievements have been made in the research of gene compression tools, but there is no plan to reduce the code stream from the encoding sequence, so there is still room for improvement in compression efficiency.

Next-generation sequencing generates tens of thousands of short reads, which are typically stored in the widely accepted text-based FASTQ format, containing all the information generated by the sequencing. Each short read contains three parts: one is metadata, which is used to describe the sequencing platform and other information; the second is the DNA base sequence, which is used to record the DNA fragments obtained in the current short read; the third is the quality score, which is used to The degree of confidence in the determination of each symbol in the corresponding DNA base sequence. The quality score data in the FASTQ format has high randomness and noise, which is related to factors such as sequencing instruments and sequencing methods. It usually contains dozens of different characters and is difficult to compress. It usually accounts for about 70% of the compressed file. Therefore, the compression result of the quality score data plays a key role in the compression effect of the entire FASTQ format data.

At present, typical methods for lossless compression of quality scores in gene sequencing data mainly include the following:

The first is to use existing text compression tools as the most commonly used compression methods for FASTQ files, such as Gzip and 7z. These methods are mainly designed to process ordinary character sequences, and do not consider the unique characteristics of quality score data. Therefore, when compressing Doesn't work well with genetic sequencing data.

The second is the improved run-length method and dictionary method for genetic data compression. In most cases, the compression effect of these methods is worse than that of the entropy coding method, and the purpose of reducing the compression rate to the greatest extent cannot be achieved.

The third is some compression algorithms for quality scores, such as Quip, etc. This method uses a high-order Markov model to predictively encode it. The model is too complex and does not take into account the impact of encoding order on compression, resulting in long compression time and poor robustness of the algorithm.

Contents of the invention

The purpose of the present invention is to overcome the shortcomings of the above-mentioned prior art, and propose a DNA sequencing quality score lossless compression method based on adaptive coding sequence, so as to maximize the compression effect without increasing the compression time.

The technical solution of the present invention is: first, extract the base sequence and quality score sequence in the FASTQ file; then calculate its mean value for each line of quality score data and quantify it, and build a prediction model according to context information, mean value information, and base information; Finally, the arithmetic coder is driven by the serpentine coding sequence to encode the sequence to achieve the purpose of compressing the quality score. The specific implementation is as follows:

(1) Extract the quality score data and base data in the FASTQ file:

(1a) Statistically analyze the characteristics of DNA sequencing data, and create two M×N coded compression blocks P ₁ and P ₂ , where M is the number of lines in the compression block, that is, the number of lines of quality score data processed at one time, and N is the compression block The number of columns, that is, the length of the mass fraction, N≤150;

(1b) extract the quality score data and base data stored in the FASTQ file through the first code compression block P ₁ and the second code compression block P ₂ respectively;

(2) calculate the mean value of each line quality score in the extracted FASTQ file of the first coded compression block _P1 and quantize, obtain the line mean value matrix F of M * 1;

(3) The context information, base information and row mean information of the encoded characters are counted and quantified uniformly, and the final encoding model is calculated:

(3a) Build a model for the current coded character q: count the first four characters q ₁ , q ₂ , q ₃ and q ₄ , and take the base information corresponding to the current character and the previous character in the second code compression block P ₂ as j ₁ and j ₂ , take the mean value _of the _row where the character q is located in the row mean matrix F and record it as f, which is _the quantized result; q ₄ takes the same sign or makes it equal to zero;

(3b) Reduce the model cost by quantifying the overall model, that is, take the larger value of the first two characters q ₁ and q ₂ as A, and the larger value of the last two characters q ₃ and q ₄ as B , creating two distinct identifiers C and D, computing the final encoding model for the current encoding symbol:

P _now ＝A·B·C·D·j ₁ ·j ₂ ·f

Wherein, when q ₁ =q ₂ , the identifier C=1, otherwise C=0; when q ₃ =q ₄ , D=1, otherwise D=0; P _now is the probability estimate of the current coding symbol;

(4) Utilize the designed final coding model to drive the self-adaptive arithmetic coder, and traverse and compress the first coded compression block P ₁ along the direction with the strongest correlation by adopting the serpentine coding sequence.

Compared with the prior art, the present invention has the following advantages:

1. Since the present invention makes full use of the probability update mechanism of the arithmetic coder, the compression rate of the quality score data in the equal-length FASTQ file is superior to all current algorithms.

2. Since the present invention compresses the mean value of each line of quality score while compressing the quality score data, it facilitates the statistics and access to the mean value in the downstream processing process.

3. Due to the simple structure of the coding model designed by the present invention, it has strong portability, and is convenient for re-optimization and integration into the compression of the entire FASTQ file. It can be widely used in various compression schemes using this module, and has good scalability.

Description of drawings

Fig. 1 is the realization flowchart of the present invention;

Fig. 2 is a schematic diagram of quantifying the quality score row mean value in the present invention;

FIG. 3 is a schematic diagram of a serpentine scanning sequence used in the present invention.

Detailed ways

The present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

With reference to Fig. 1, the realization steps of the present invention are as follows:

Step 1, extract the quality score data and base data in the FASTQ file.

Genetic sequencing generates thousands of short reads, which are typically stored in the widely accepted text-based FASTQ format, containing all the information generated by the sequencing. In the FASTQ file format, each short read consists of four lines, each line separated by a newline, where:

The first line starts with the '@' character, followed by a unique sequence ID identifier and optional sequence description content, and the identifier and description characters are separated by spaces;

The second line is the nucleotide sequence, representing the base data, consisting of a sequence consisting of only five characters {'A', 'T', 'C', 'G', 'N'}, where the character 'N' Indicates an ambiguous base;

The third line starts with the character '+', followed by the sequence logo and description information, or no information, as a separator;

The last line is the quality score line, each character corresponds to the quality of the base at the corresponding position in the second line, and the quality score corresponds to the number Q = -10log 10P, where P indicates the probability that the corresponding nucleotide in the read is wrong. Quality scores are usually denoted using ASCII letters [33:73] or [64:104] and are used both for quality control of raw data and for downstream processing.

The specific implementation of this step is as follows:

1.1) Statistically analyze the characteristics of DNA sequencing data, and create two coded compression blocks P ₁ and P ₂ of M×N size, where M is the number of lines in the compressed block, that is, the number of lines for processing quality score data at one time, and N is the compression The number of columns of the block, that is, the length of the quality fraction, N≤150;

1.2) Extract the quality score data and base data stored in the FASTQ file through the first code compression block _P1 and the second code compression block _P2 respectively;

Since the quality score characters in most FASTQ files are less than 40 and the jumps are not large, a good prediction model can be designed according to the correlation between data to improve the compression effect. At the same time, considering that too many related characters will not only lead to increased time and calculation complexity, but also sometimes cause model cost problems, it is necessary to use appropriate compression blocks to count the correlation between quality scores. When computing resources allow Within the range of , the larger the compression block is designed, the better the compression effect, but in order not to exceed the maximum memory, in this embodiment, a compression block of 2000000×160 is used. The total number of models is set to 40×40×40×16. In the actual compression process, data of one compressed block size is processed each time until the end of the file.

Step 2: Calculate and quantize the mean value of the quality score of each row in the FASTQ file extracted from the first coded compression block P ₁ to obtain an M×1 row mean value matrix F.

2.1) Perform an averaging operation on each line of the first coded compressed block P ₁ whose size is M×N, add and divide the N quality score values of each line by the total number N to obtain the mean value of the quality score of each line;

2.2) Quantify and store the obtained quality score values of each row:

Referring to Figure 2, after counting the mean value of the quality score of each row, clustering is performed according to the distribution of the mean value, the mean value with a large number is subdivided, and the quality values with a small number and low value are merged to facilitate Improve coding efficiency. For a specific compressed file, a unique quantization method can be designed according to the mean distribution to achieve the optimal effect, but this not only increases the calculation amount but also increases a lot of calculation time. Therefore, this example uses a quantization method with strong scalability and easy implementation, that is, two adjacent mean values are regarded as the same situation, and the part with a small quality value and a low quantity is regarded as a part as a whole. Summarizing the quantitative experience, the quantitative results are obtained as follows:

If f _i <(num-15), then f _i =(num-15);

If (num-15)≤f _i <(num-13), then f _i = (num-13);

If (num-13)≤f _i <(num-11), then f _i = (num-11);

If (num-11)≤f _i <(num-9), then f _i = (num-9);

If (num-9)≤f _i <(num-7), then f _i =(num-7);

If (num-7)≤f _i , then f _i =(num-6);

Among them, num is the total number of encoding symbols 40, _fi is the mean value of the current line, and the value of i is [1,M];

Merge the quantized mean values f _i of each row in a column-arranged manner to obtain an M×1 row mean value matrix F.

Step 3, the context information, base information and row mean information of the coded characters are counted and quantified uniformly, and the final coding model is calculated.

3.1) Build a model for the current coded character q: count the first four characters q ₁ , q ₂ , q ₃ and q ₄ , and take the base information corresponding to the current character and the previous character in P ₂ as j ₁ and j ₂ , The mean value of the row where the character q is located is recorded as f, where f is the quantized result; for marginal characters lacking the above information, its q ₁ , q ₂ , q ₃ and q ₄ can take the same sign or make them equal to zero.

For example: given the specific content of the first encoded compressed block _P1 is: E, F, G, H, I; the specific content of the second encoded compressed block _P2 is: A, T, C, G, G;

When the coding model is established for the third character "G" in the first coded compression block _P1 , the values of the first four characters are: q ₁ =F, q ₂ =E, q ₃ =0, q ₄ =0; the value of the base information corresponding to the current character and the previous character is: j ₁ =C, j ₂ =T; the mean value of the row where it is located is: f=mean(ASCII(E)+ASCII(F )+ASCII(G)+ASCII(H)+ASCII(I));

When the coding model is established for the fifth character "I" in P ₁ , the values of the first four characters are respectively: q ₁ =H, q ₂ =G, q ₃ =F, q ₄ =E; its current The value of the base information corresponding to the character and the previous character is: j ₁ =G, j ₂ =G; the mean value of the row where it is located is: f=mean(ASCII(E)+ASCII(F)+ASCII(G )+ASCII(H)+ASCII(I));

It can be seen that the same method can be used to build a model for marginal characters lacking the above information.

3.2) Considering the actual situation that the total number of models is limited, reduce the model cost by quantifying the overall model, that is, take the larger value of q ₁ and q ₂ as A, and the larger value of q ₃ and q ₄ as B, Create two different identifiers C and D, C is used to judge whether q ₁ and q ₂ are equal, and D is used to judge whether q ₃ and q ₄ are equal. Therefore, the final selected model of the current coding symbol is: P _now =A·B·C·D·j ₁ ·j ₂ ·f.

Among them, P _now is the probability estimation value of the current coding symbol.

Step 4, use the designed final coding model to drive the adaptive arithmetic coder, and use the serpentine coding order to traverse and compress the first coding compression block P ₁ along the direction with the strongest correlation.

4.1) Obtain a more accurate probability estimate P _now of the current coded character through the above-mentioned final coding model, and send it to the adaptive arithmetic coder as the optimal prediction result;

4.2) The encoder performs traversal encoding and compression:

When encoding, it is necessary to encode and scan characters one by one. The traditional scanning method defaults to traversing line by line. After traversing a whole line, continue scanning from the starting position of the second line. This example uses column-by-column scanning, and after one column is encoded, the end of the next column is used as the starting point, and the reverse upward traverse is performed, and the whole process looks like a serpentine scan, as shown in Figure 3. The ultimate lossless compression is achieved by traversing all characters.

The above description is only a specific example of the present invention, and does not constitute any limitation to the present invention. Obviously, for those skilled in the art, after understanding the content and principles of the present invention, they can make various modifications and changes in form and details without departing from the principles and structures of the present invention, but these are based on the present invention. The modification and change of the inventive concept are still within the protection scope of the claims of the present invention.

Claims (5)

1. A DNA sequencing quality score lossless compression method based on an adaptive coding sequence, characterized in that, comprising the following: (1) Extract the quality score data and base data in the FASTQ file: (1a) Statistically analyze the characteristics of DNA sequencing data, and create two M×N coded compression blocks P ₁ and P ₂ , where M is the number of lines in the compression block, that is, the number of lines of quality score data processed at one time, and N is the compression block The number of columns, that is, the length of the mass fraction, N≤150; (1b) extract the quality score data and base data stored in the FASTQ file through the first code compression block P ₁ and the second code compression block P ₂ respectively; (2) calculate the mean value of each line quality score in the extracted FASTQ file of the first coded compression block _P1 and quantize, obtain the line mean value matrix F of M * 1; (3) The context information, base information and row mean information of the encoded characters are counted and quantified uniformly, and the final encoding model is calculated: (3a) Build a model for the current coded character q: count the first four characters q ₁ , q ₂ , q ₃ and q ₄ , and take the base information corresponding to the current character and the previous character in the second code compression block P ₂ as j ₁ and j ₂ , take the mean value _of the _row where the character q is located in the row mean matrix F and record it as f, which is _the quantized result; q ₄ takes the same sign or makes it equal to zero; (3b) Reduce the model cost by quantifying the overall model, that is, take the larger value of the first two characters q ₁ and q ₂ as A, and the larger value of the last two characters q ₃ and q ₄ as B , creating two distinct identifiers C and D, computing the final encoding model for the current encoding symbol: P _now ＝A·B·C·D·j ₁ ·j ₂ ·f Wherein, when q ₁ =q ₂ , the identifier C=1, otherwise C=0; when q ₃ =q ₄ , D=1, otherwise D=0; P _now is the probability estimate of the current coding symbol; (4) Utilize the designed final coding model to drive the self-adaptive arithmetic coder, and traverse and compress the first coded compression block P ₁ along the direction with the strongest correlation by adopting the serpentine coding order. 2. The method according to claim 1, wherein the DNA sequencing data feature in (1a) means that the DNA sequencing data contains thousands of reads, each read has four lines, and the second line is quality Score data, the fourth line is base data, the overall DNA sequencing quality score data is encoded in ASCII code, and the number of encoding symbol types = maximum value - minimum value + 1 quality score data. 3. method according to claim 1, it is characterized in that, in (2), the mean value of every row's quality fraction extracted to the first encoding compression block _P1 is quantized, and its realization is as follows: If f _i <(num-15), then f _i =(num-15); If (num-15)≤f _i <(num-13), then f _i = (num-13); If (num-13)≤f _i <(num-11), then f _i = (num-11); If (num-11)≤f _i <(num-9), then f _i = (num-9); If (num-9)≤f _i <(num-7), then f _i =(num-7); If (num-7)≤f _i , then f _i =(num-6); Among them, f _i is the mean value of the quality score of each row, num is the total number of coded symbols whose size is 40, and the value of i is [1,M], and each row f _i is quantized to obtain an M×1 row mean matrix F. 4. method according to claim 1, is characterized in that, the final encoding model of utilization design in (4) drives adaptive arithmetic coder, refers to the P _now that current symbol is carried out probability estimation as optimal prediction result into the adaptive arithmetic coder. 5. The method according to claim 1, characterized in that, in (4), the serpentine coding order is adopted to traverse and compress the first coding compression block _P1 along the direction with the strongest correlation, which means traversing the first coding and compressing Block P ₁ is scanned column by column from top to bottom, and then reversely traverses from bottom to top after scanning a column, and reciprocates in turn until the entire compressed block is traversed.

Images