Compression and decompression method for DNA sequencing data

Abstract

The invention discloses a compression and decompression method for DNA sequencing data, and belongs to the technical field of biological information. The technical problems that in the prior art, when third-generation sequencing data are compressed, flexibility is poor, and the compression rate is low are solved. According to the compression method for the DNA sequencing data, global compression or block compression can be dynamically selected before compression, different compression methods are adopted for different types of data, and the flexibility of the compression method is improved; when the base sequence is compressed, the adopted self-indexing structure is a lossless compression method, the initial positions of all the sequences can be compressed and restored by using smaller data volume, and the compression rate is improved. Due to the fact that a special design is adopted in a run length coding structure, a character part is removed, an obtained run length segment is shorter, lossless compression of the position and the sequence is guaranteed, and meanwhile a better compression effect is provided. The method is mainly used for compression and decompression of DNA sequencing data.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to a method for compressing and decompressing DNA sequencing data. Background Technology

[0002] DNA sequencing converts the sequence of the four bases "A," "T," "C," and "G" on a DNA strand into a digital format that can be processed and stored by a computer. Third-generation sequencing (TGS) is a new sequencing technology that differs from traditional first-generation sequencing and second-generation high-throughput short-read sequencing. Third-generation sequencing produces sequences with long reads, non-uniform read lengths, high throughput but high redundancy, leading to a sharp increase in storage costs. Traditional compression tools are primarily designed for second-generation short-read data, and their compression strategies are inflexible and have low compression rates when processing third-generation sequencing data. Summary of the Invention

[0003] To overcome the technical problems of poor flexibility and low compression ratio in the compression of third-generation sequencing data in existing technologies, this invention provides a method for compressing and decompressing DNA sequencing data.

[0004] This invention is achieved through the following technical solution:

[0005] A method for compressing DNA sequencing data includes the following steps:

[0006] The raw sequencing data was aligned with the reference genome to generate a standardized file, which was then divided into N independent files. Each independent file was further divided into three parts, and each part was compressed as follows:

[0007] Header file and sequence compression: Preprocess each sequence by converting the base characters in each sequence into the first and second characters and recording the generated abnormal information; construct an ascending set X from all sequences; write the data in set X into matrix Y; vertically read matrix Y to obtain the character stream file Z; apply a customized run-length encoding algorithm to process the character stream file Z; then use Huffman coding to process the generated abnormal information and character stream file Z respectively.

[0008] Quality value compression: Converts the quality value to a WebP format image;

[0009] Apart from the header file, sequence, and quality value data mentioned above, all other data is retained directly.

[0010] Furthermore, the preprocessing involves determining the match between each base sequence recorded in the CIGAR field and the reference genome, converting the base sequence into a sequence containing only the first and second characters, and generating anomaly information. Specific steps include:

[0011] When the CIGAR field contains "*", the entire original sequence is deleted from the SEQ field, and the entire original base sequence is written into the error message.

[0012] When the CIGAR field contains “uM” or “u=", replace the corresponding u bases in the SEQ field with u first characters;

[0013] When the CIGAR field contains “uD”, replace the corresponding u bases in the SEQ field with u second characters, and append u “D” characters to the exception information;

[0014] When the CIGAR field contains “uI”, replace the corresponding u bases and the immediately upstream base in the SEQ field with a second character, and append an “I” character and the u original bases in the SEQ field at the current position to the error information;

[0015] When the CIGAR field contains “uS” and appears at the beginning of the sequence, replace the corresponding u bases in the SEQ field and its immediately downstream base with a second character, and append an “S” character and the original u bases at the current position to the abnormal information.

[0016] When the CIGAR field contains “uS” and appears at the end of the sequence, replace the corresponding u bases in the SEQ field and its immediate upstream base with a second character, and append an “S” character and the original u bases at the current position to the abnormal information;

[0017] When the CIGAR field contains “uH”, the corresponding u original bases are deleted from the SEQ field, and one “H” character is appended to the error message.

[0018] When the CIGAR field contains “uX”, replace the corresponding u bases in the SEQ field with u second characters, and append one “X” character and the u original bases at the current position to the exception information;

[0019] When the CIGAR field contains “uP”, append u “P” characters to the exception information.

[0020] Furthermore, writing the data from set X into matrix Y means writing the first and second characters of each sequence of data in set X into matrix Y. The specific operation includes the following steps:

[0021] Create an ordered set X a Iterate through each site k of the reference genome sequentially, performing the following operations during each iteration:

[0022] (a) Dynamically maintain the ordered set X a Add the sequence with starting position k from set X to the ordered set X. a At the end of set X, remove the sequence with ending position k from the ordered set X. a Delete;

[0023] (b) Reverse prefix sorting and record index mapping: Calculate the ordered set X a The offset value between the start position and k of each sequence is used to extract set X. a The character whose index in each sequence is equal to the offset value is used, and then the ordered set X is sorted according to the lexicographical order of the characters. a Update the sorting to get X b And record the same sequence in set X a With X b Index mapping relationships in the database;

[0024] (c) Matrix write operation: Traverse and update the sorted set X b During each traversal, the offset value of the next position after the starting position of the unterminated sequence is calculated: if the offset value is less than the length of the current sequence, the character with the offset value in the current sequence is written into matrix Y, and the row index is the current sequence in set X. a The index value in the sequence is k+1; if the offset value is greater than or equal to the current sequence length, the sequence is marked as terminated.

[0025] Furthermore, the specific process of vertically reading matrix Y to obtain the character stream file Z includes:

[0026] Traverse each column of matrix Y according to the order of the starting alignment sites of each sequence, and read each row in each column from top to bottom:

[0027] a: If the character read is the starting character of the sequence, then write the third character into the character stream Z and then write the currently read character;

[0028] b: If the character read is the last character of the sequence, calculate the difference x between the row index of the current character and the row index of the previous ending sequence. If the difference is positive, write the currently read character into the character stream Z, then write x fourth characters and one sixth character; if the difference is negative, write the currently read character into the character stream Z, then write |x| fifth characters and one sixth character.

[0029] c: If the character read is not the start or end character, then write the character directly to the character stream Z.

[0030] Furthermore, the specific execution steps of the customized run-length encoding algorithm include:

[0031] The character stream Z is fused and encoded into a run-length encoded string in pure integer sequence format according to the encoding rules, and then the alignment point of the first sequence in set X is added to the beginning of the run-length encoded string.

[0032] The encoding rules include:

[0033] If the third character is followed by the first character, the third character is encoded as "0"; otherwise, the third character is encoded as "-1".

[0034] If the fourth character is followed by the first character, the fourth character is encoded as "-2"; otherwise, the fourth character is encoded as "-3".

[0035] If the fifth character is followed by the first character, the fifth character is encoded as "-4"; otherwise, the fifth character is encoded as "-5".

[0036] The first or second character is encoded as the number of times the character itself is repeated.

[0037] Furthermore, the specific steps for converting the quality value to WebP format include:

[0038] The quality values ​​are converted into integer data and then allocated to the RGB channels in sequence. Each set of integer data in the R, G, and B channels is encoded into an image pixel, and then an image in WebP format is generated.

[0039] A method for decompressing compressed files of third-generation DNA sequencing data, wherein the decompression method is used to decompress compressed files generated by any of the compression methods described above;

[0040] If the file to be decompressed is a single file, decompress it directly; if the file to be decompressed consists of multiple block files, filter the block files in the target range and decompress each block file separately.

[0041] Perform Huffman decoding on the header file and the compressed sequence data to obtain an integer string; decode the integer string using a custom run-length decoding method to obtain a character stream Z′; construct a matrix Y′ based on the character stream Z′; use reverse inverse prefix matching to operate on the matrix Y′ to generate a set X′ containing only the first and second character sequences;

[0042] Perform Huffman decoding on the abnormal information data, and combine the set X′ with the decoded abnormal information to restore the original base sequence;

[0043] Compressing quality values ​​involves parsing the RGB channel values ​​of each pixel in a WebP image into unsigned integers and then encoding them into a quality value string.

[0044] Furthermore, the specific steps for decoding the integer string into the character stream Z′ include:

[0045] Extract the first integer from the integer string as the starting position pos′ and write it to a temporary file; create a dynamically convertible character pointer P, and mark the character pointed to as the first character or the second character by controlling the parity of P;

[0046] Read the integer string sequentially from the second element to the end, performing the following operations during the reading process:

[0047] When a non-positive integer A is read: reset pointer P = A, then read the next integer B, and write B special characters corresponding to A according to the decoding rules;

[0048] When a positive integer C is read: write the first or second character of the currently marked character to C pointers P, and increment pointer P;

[0049] The decoding rules include:

[0050] If the element is 0 or -1, decode the element into the third character;

[0051] If the element is -2 or -3, decode the element into the fourth and sixth characters, respectively.

[0052] If the element is -4 or -5, decode the element into the fifth and sixth characters.

[0053] Furthermore, the step of constructing matrix Y' based on character stream Z' involves vertically filling matrix Y' based on character stream Z' and the starting position pos' recorded in the temporary file. The specific steps for filling each column include:

[0054] Initialize dynamic state variables: fill the state array full; counter for rows to be filled y; index of the end row δ; position pointer p′ = pos′;

[0055] Parse the character stream Z' character by character and update the state variables:

[0056] When encountering 'a' third characters, perform the following operations: Query the first 'y' elements in the array 'full' that have a value of false, and take the index L of the last element among the 'y'; get the number N of elements in the row corresponding to the index in matrix Y'; update the position pointer p' = p' + N, and record p' in a temporary file; update the value of 'a' to a + y, and the element at index L in the array 'full' has a value of true;

[0057] When the fourth or fifth character is encountered, the value of 'a' is reduced by the number of consecutive occurrences of the current fourth or fifth character; each time a fourth character is encountered, the value of the integer δ is incremented by 1, and the value of the element with index δ in the array 'full' is set to true; each time a fifth character is encountered, the value of the integer δ is decremented by 1, and the value of the element with index δ in the array 'full' is set to true.

[0058] When the first or second character is encountered, vertically fill the active row of the current column of matrix Y′.

[0059] Furthermore, the specific steps for generating a set X′ containing only the first and second character sequences using the reverse inverse prefix matching operation matrix Y′ include:

[0060] Create an intermediate matrix I with the same dimensions as Y′; read the starting positions of all sequences from the temporary file and sort them in ascending order;

[0061] Sequentially traverse each site k of the reference genome, performing the following operations during each traversal:

[0062] (a) Dynamically constructing an ordered set X a ′: Add the sequence with starting position k to the ordered set X. a At the end of ', remove the sequence with ending position equal to k from the ordered set X. a Delete in '

[0063] (b) For set X a Perform reverse prefix sorting and record the index mapping: Calculate the ordered set X a The offset value between the starting position and k of each sequence in ' is used to extract set X. a In the sequence X, the character whose index is equal to the offset value is used. Then, based on the lexicographical order of these characters, the ordered set X is... a Update the sorting to get set X b ′, and record the same sequence in set X a ′ and X b The index mapping relationship in ';

[0064] (c) Reverse write operation: Traverse and update the sorted set X b During each iteration, the index mapping relationship and the index n of the current sequence are used to determine the position of the current sequence in set X. a The index m in matrix Y' is used to retrieve the character at column index k and row index m in matrix Y', and write it to the position at column index k and row index n in matrix I.

[0065] After the site traversal of the reference genome is completed, matrix I is read row by row to obtain set X′.

[0066] The beneficial effects of this invention are:

[0067] This invention provides a compression method for DNA sequencing data. Before compression, it allows for dynamic selection of global compression or block compression, employing different compression methods for different data types to enhance the flexibility of the compression approach. When compressing base sequences, a self-indexed structure is used, a lossless compression method that can compress and restore the start positions of all sequences with a smaller data volume, improving the compression ratio. Due to the special design of the run-length encoding structure, which removes the character portion, the resulting run fragments are shorter, providing better compression performance while ensuring lossless compression of position and sequence. Furthermore, the compression of the quality fraction portion in this invention is also lossless.

[0068] Furthermore, this invention also provides a decompression method for DNA sequencing data, primarily used for decoding sequencing data compressed using this invention. This decompression method selects the appropriate decompression method based on the compression method and employs different decompression methods for different data types. If block compression was used during compression, the decompression interval can be selected as needed during decompression, offering high flexibility. Attached Figure Description

[0069] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0070] Figure 1 This is a schematic diagram of the pre-processed base sequence during compression in one embodiment of the present invention;

[0071] Figure 2 This is a schematic diagram of run-length encoding obtained during the compression process according to an embodiment of the present invention. Detailed Implementation

[0072] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0073] The terms “first,” “second,” “third,” “fourth,” etc. (if present) in the specification, claims, and accompanying drawings of this application are used to distinguish similar objects and are not necessarily used to describe a particular order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this application described herein can be implemented, for example, in a sequence other than those illustrated or described herein. Specific Implementation Method 1

[0075] The following is combined with Figure 1 , Figure 2 A specific embodiment of a method for compressing DNA sequencing data according to the present invention is described below.

[0076] In this embodiment, for ease of understanding, “E” represents the first character, “W” represents the second character, “#” represents the third character, “$” represents the fourth character, “%” represents the fifth character, and “|” represents the sixth character. This does not represent a limitation on the method of the present invention.

[0077] The raw sequencing data is aligned with the reference genome to generate a standardized file, which is then divided into N independent files. When N is 1, it means that the standardized file is not divided into blocks, which is a global compression mode. When N is greater than 1, it means that the standardized file is divided into N blocks according to the preset block size.

[0078] The standardized file format can be SAM or BAM; in this embodiment, SAM format is used. The SAM file, separated by tabs, records each DNA read fragment and its alignment with the reference genome. It is divided into two parts: a header and alignment information. The header provides annotation information such as the SAM format version used, the aligned reference genome information, and sequencing information. The alignment information contains 11 required and optional parameters. The required parameter fields used in this embodiment include: POS, CIGIA, SEQ, and QUAL. The POS field stores the starting coordinates of the sequence alignment to the reference genome; CIGIA stores the alignment information string, reflecting the sequence differences between the sequence and the reference genome starting at the POS position; SEQ stores the sequence base string; and QUAL stores the ASCII quality string corresponding to the sequence bases in the SEQ field, reflecting the sequencing accuracy of each base.

[0079] For each individual file, it is split into three parts: header and sequence, quality values, and directly retained portion. The header and sequence portion includes the header of the SAM file and the contents of the POS, CIGIA, and SEQ fields; the quality values ​​portion includes the contents of the QUAL field; the directly retained portion refers to the contents of other necessary fields in the SAM file besides the header and sequence and quality values ​​portions. Since this part of the data accounts for a very small proportion of the SAM file, the raw data is saved directly without compression.

[0080] Compress each part after splitting:

[0081] (I) Compression of header files and sequences:

[0082] (1) Preprocess each sequence by converting the base characters in each sequence into {E,W} characters and recording the generated abnormal information; construct all sequences into an ascending set X.

[0083] In this field, the data in a SAM file is collectively referred to as reads, and each piece of data separated by tabs is called a read. Each sequence mentioned here refers to each read data.

[0084] The preprocessing involves determining the match between each base sequence recorded in the CIGAR field and the reference genome, converting the base sequences into {E,W} sequences, and generating anomaly information. In the SAM file used in this embodiment, the value range of CIGAR is *|([0-9]+[MIDNSHPX=])+. Assuming u and r are integers, the specific preprocessing steps include:

[0085] When the CIGAR field contains "*", it indicates that the sequence has a very low match with the reference genome. The entire original sequence will be deleted from the SEQ field and written into the abnormal information.

[0086] When the CIGAR field contains “uM” or “u=", it means that there are u consecutive base characters matching the reference genome at the current position, and the u bases in the SEQ field are replaced with u “E” characters.

[0087] When the CIGAR field contains “uD”, it indicates that there is a missing u bases at the current position. The u bases in the SEQ field are replaced with u “W” characters, and u “D” characters are appended to the exception information.

[0088] When the CIGAR field contains “uI”, it indicates that there is an insertion of u consecutive bases at the current position. Replace the corresponding u bases in the SEQ field and the one base immediately upstream of them with a “W” character. At the same time, add an “I” character and the u original bases in the SEQ field at the current position to the error information.

[0089] When the CIGAR field contains “uS” and appears at the beginning of the sequence, it indicates that there are u consecutive bases that cannot be matched at the current position. Replace the corresponding u bases in the SEQ field and the next downstream base with a “W” character, and append a “S” character and the original u bases at the current position to the error information.

[0090] When the CIGAR field contains “uS” and appears at the end of the sequence, it indicates that there are u consecutive bases that cannot be matched at the current position. Replace the corresponding u bases in the SEQ field and the one base immediately upstream of them with a “W” character. At the same time, add a “S” character and the u original bases at the current position to the error information.

[0091] When the CIGAR field contains “uH”, it means that there are u consecutive bases that cannot be matched at the current position. The corresponding u original bases will be deleted from the SEQ field, and an “H” character will be added to the error information.

[0092] When the CIGAR field contains “uX”, it indicates that there is a replacement anomaly at the current position where there are u consecutive bases. The u bases in the SEQ field are replaced with u “W” characters, and one “X” character and the u original bases at the current position are appended to the anomaly information.

[0093] When the CIGAR field contains "uP", it usually appears together with the insertion exception "rI", indicating that there are u empty paddings before the r bases inserted at the current position. In actual processing, the read is treated as an insertion exception, and u "P" characters are appended to the exception information.

[0094] The above preprocessing operation is based on the storage structure characteristics of SAM files. It transforms the long read sequence and the alignment operation information contained in the CIGAR field into sequences containing only two characters. Through the above preprocessing operation, the base sequence of the SEQ field in each read is transformed into the {E,W} character sequence, which can make all E and W sequences accurately align to the alignment site of the reference genome within their original matching intervals without generating gaps or insertions, thus facilitating subsequent maximum inverse prefix matching change operations.

[0095] After each read data in reads has been preprocessed, a set X is constructed using all the read data. Then, the data in set X is sorted in ascending order according to the POS value of each read data, i.e., the starting alignment point, to obtain the ascending set X.

[0096] After the set X is constructed, create a matrix Y with the same size as set X.

[0097] (2) Writing data from set X into matrix Y means writing the {E,W} characters from each sequence of data in set X into matrix Y. The specific operation includes the following steps:

[0098] Create an ordered set X a The target matrix Y is obtained by sequentially traversing each site k of the reference genome and completing the traversal.

[0099] On each iteration, perform the following operations:

[0100] (a) Dynamically maintain the ordered set X a Add the sequence in set X whose starting alignment point is k to the ordered set X. a At the end of set X, remove the sequence with ending position k from the ordered set X. a Delete;

[0101] (b) Reverse prefix sorting and record index mapping: Calculate the ordered set X a The offset value between the starting alignment site and k of each sequence is used to extract set X. a The character whose index in each sequence is equal to the offset value is used, and then the ordered set X is sorted according to the lexicographical order of the characters. a Update the sorting to get X b And record the same sequence in set X a With X b The index mapping relationship in the data; the set X at this time. b In the sequence, all E characters precede all W characters.

[0102] (c) Matrix write operation: Traverse and update the sorted set X b During each traversal, the offset between the starting alignment site of the unterminated sequence and the (k+1)th site of the reference genome is calculated. If the offset is less than the current sequence length, the character indexed by the offset is written into matrix Y, and the row index is the current sequence in set X. a The index value in the sequence is k+1; if the offset value is greater than or equal to the current sequence length, the sequence is marked as terminated.

[0103] During the traversal of reference genome locus k, when locus k satisfies k = max(pos)i +length i When -1), the traversal of point k ends. This operation can terminate the current traversal early, reducing time overhead. In the above formula, pos i Indicates the starting alignment position of the i-th sequence, length i This represents the length of the i-th sequence.

[0104] The algorithm is based on the similarity characteristics of adjacent sequences in sequencing data: for two adjacent sequences, the higher the similarity in their prefix sequences, the greater the probability that they will have the same bases at subsequent base positions. Therefore, this invention continuously maintains the maximum common reverse substring among sequence groups at each position k in the reference genome to maximize the probability that adjacent sequences at position k+1 have the same bases, thereby generating longer and fewer run-length encoded blocks.

[0105] (3) The specific process of vertically reading matrix Y to obtain character stream file Z includes:

[0106] Traverse each column of matrix Y according to the starting alignment site of each sequence, and read each row in each column from top to bottom:

[0107] a: If the character read is the starting character of the sequence, write the "#" character into the character stream Z and then write the currently read character;

[0108] b: If the character read is the last character of the sequence, calculate the difference x between the row index of the current character and the row index of the previous ending sequence. If the difference is positive, write the currently read character into the character stream Z, followed by x "$" characters and one "|" character; if the difference is negative, write the currently read character into the character stream Z, followed by |x| "%" characters and one "|".

[0109] c: If the character read is not the start or end character, then write the character directly to the character stream Z.

[0110] (4) Apply a customized run-length encoding algorithm to process the character stream file Z. The specific process includes:

[0111] The character stream Z is fused and encoded into a run-length encoded string in pure integer sequence format according to the encoding rules. Then, the alignment point of the first sequence in set X is added to the beginning of the run-length encoded string.

[0112] The encoding rules include:

[0113] If "#" is followed by "E", "#" is encoded as "0"; if "#" is followed by something other than "E", "#" is encoded as "-1".

[0114] If "$" is followed by "E", then "$" is encoded as "-2". If there is no character following "$", then the "$" character is encoded as "-3".

[0115] If "%" is followed by "E", "%" is encoded as "-4"; if "#" is followed by something other than "E", "%" is encoded as "-5".

[0116] "E" or "W" is encoded as the number of times the character itself is repeated.

[0117] (5) Use Huffman coding to process the generated exception information and character stream file Z respectively.

[0118] After the above operations (1)-(5), the compression of header files, sequences, and exception information is completed.

[0119] (ii) Compression of the mass value portion:

[0120] The quality values ​​are converted into integer data and then allocated to the RGB channels in sequence. Each set of integer data in the R, G, and B channels is encoded into an image pixel, and then a WebP format image is generated. This image is the compressed quality value file.

[0121] (iii) In the SAM file, apart from the header file and sequence and quality value data mentioned above, other data is directly retained. Specific Implementation Method Two

[0123] The DNA sequencing data decompression method provided in this embodiment is used to decompress the compressed file generated by the DNA sequencing data compression method provided in this invention.

[0124] In this embodiment, for ease of understanding, “E” represents the first character, “W” represents the second character, “#” represents the third character, “$” represents the fourth character, “%” represents the fifth character, and “|” represents the sixth character. This does not represent a limitation on the method of the present invention.

[0125] If the file to be decompressed is a single file, decompress it directly. If the file consists of multiple chunks, select the chunks within the target range and decompress each chunk separately. The decompression method is the same for both single-file global decompression and individual chunk decompression. During compression, the SAM data was split into streams; correspondingly, the data to be decompressed also needs to be split during decompression. Different data types use different compression methods, and different decompression methods are used for different compression methods.

[0126] The decompression method includes the following steps:

[0127] Decompression of header files and sequence compressed data:

[0128] (1) Perform Huffman decoding on the header file and the compressed sequence data to obtain an integer string; then decode the integer string using a custom run-length decoding method to obtain the character stream Z′. Specific steps include:

[0129] Extract the first integer from the integer string as the starting position pos′ and write it to a temporary file; create a dynamically convertible character pointer P, and mark the pointed-to character as "E" or "W" by controlling the parity of P. In this embodiment, when P is even, the pointed-to character is marked as "E"; when P is odd, the pointed-to character is marked as "W".

[0130] Read the integer string sequentially from the second element to the end, performing the following operations during the reading process:

[0131] When a non-positive integer A is read: reset pointer P = A, then read the next integer B, and write B special characters corresponding to A according to the decoding rules;

[0132] When a positive integer C is read: write the character "E" or "W" currently marked on the pointer P, and increment the pointer P.

[0133] The decoding rules used in this embodiment include:

[0134] If the element is 0 or -1, decode the element as "#";

[0135] If the element is -2 or -3, decode the element as "$|";

[0136] If the element is -4 or -5, decode the element as "%|".

[0137] (2) Constructing matrix Y′ based on character stream Z′ involves filling matrix Y′ vertically based on character stream Z′ and the starting position pos′ recorded in the temporary file. The specific steps for filling each column include:

[0138] Initialize dynamic state variables: fill the state array full, with all default values ​​set to false; the counter for rows to be filled, y, is initialized to 0; the index of the end row, δ, is initialized to -1; the position pointer p′ = pos′;

[0139] Parse the character stream Z' character by character and update the state variables:

[0140] When encountering a "#" characters, perform the following operations: Query the first y elements of the array full that have a value of false, and take the index L of the last element among the y elements; get the number N of elements in the row corresponding to index L in matrix Y′; update the position pointer p′ = p′ + N, and record p′ in a temporary file; update the value of a to a + y, and the element with index L in the array full has a value of true.

[0141] In the above operation, when y = 0, it means that the matrix Y′ has not yet been filled.

[0142] When a "$" or "%" character is encountered, the value of a is reduced by the number of consecutive occurrences of the current "$" or "%" character; for each "$" character encountered, the value of the integer δ is increased by 1, and the value of the element with index δ in the array full is set to true; for each "%" character encountered, the value of the integer δ is decreased by 1, and the value of the element with index δ in the array full is set to true.

[0143] When the character "E" or "W" is encountered, vertical filling is performed into the active row of the current column of matrix Y'. The active row refers to the row containing the sequence of characters that has not yet ended and is yet to be filled in the current column of matrix Y'.

[0144] (3) Use the reverse inverse prefix matching operation matrix Y′ to generate a set X' containing only the {E,W} sequences. Specific steps include:

[0145] Create an intermediate matrix I with the same dimensions as Y′; read the starting positions of all sequences from the temporary file and sort them in ascending order;

[0146] Sequentially traverse each site k of the reference genome, performing the following operations during each traversal:

[0147] (a) Dynamically constructing an ordered set X a ′: Add the sequence with starting position k to the ordered set X. a At the end of ', remove the sequence with ending position equal to k from the ordered set X. a Delete in '

[0148] (b) For set X a Perform reverse prefix sorting and record the index mapping: Calculate the ordered set X a The offset value between the starting position and k of each sequence in ' is used to extract set X. a In the sequence X, the character whose index is equal to the offset value is used. Then, based on the lexicographical order of these characters, the ordered set X is... a The updated sorting yields set X. b ′, and record the same sequence in set X a ′ and X b The index mapping relationship in ';

[0149] (c) Reverse write operation: Traverse and update the sorted set X b During each iteration, the index mapping relationship and the index n of the current sequence are used to determine the position of the current sequence in set X. aThe index m in matrix Y' is used to retrieve the character at column index k and row index m in matrix Y', and write it to the position at column index k and row index n in matrix I.

[0150] After the site traversal of the reference genome is completed, matrix I is read row by row to obtain the {E, W} sequence set X'.

[0151] Perform Huffman decoding on the abnormal information data, and combine the {E,W} sequence set X′ with the decoded abnormal information to reconstruct the original base sequence. The specific steps include:

[0152] When the character is "E", it is restored to the original bases of the corresponding reference genome;

[0153] When the character is "W", it is restored to the corresponding character in the corresponding error message.

[0154] After decompressing the compressed header and sequence files, delete the temporary files.

[0155] For decompressing the quality score compressed data, the pixel data of the WebP image is parsed into an unsigned integer sequence in the order of R->G->B channels, with each integer corresponding to an ASCII code value; finally, these ASCII code values ​​are decoded and merged to restore the original quality score string.

[0156] To verify the beneficial effects of the present invention, the following experiments were conducted:

[0157] The 10X–50X sequencing data of the whole genome of *E. coli* NC_08253 generated by PBSIM was compressed using the compression method provided in this invention. The compression results are shown in the table below, in MB:

[0158]

[0159] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A method for compressing DNA sequencing data, characterized in that, Includes the following steps: The raw sequencing data was aligned with the reference genome to generate a standardized file, which was then divided into N independent files. Each independent file was further divided into three parts, and each part was compressed as follows: Header file and sequence compression: Preprocess each sequence by converting the base characters in each sequence into the first and second characters and recording the generated abnormal information; construct an ascending set X from all sequences; write the data in set X into matrix Y; vertically read matrix Y to obtain the character stream file Z; apply a customized run-length encoding algorithm to process the character stream file Z; then use Huffman coding to process the generated abnormal information and character stream file Z respectively. Quality value compression: Converts the quality value to a WebP format image; Apart from the header file, sequence, and quality value data mentioned above, all other data is retained directly.

2. The method for compressing DNA third-generation sequencing data according to claim 1, characterized in that, The preprocessing involves determining the match between each base sequence recorded in the CIGAR field and the reference genome, converting the base sequence into a sequence containing only the first and second characters, and generating anomaly information. Specific steps include: When the CIGAR field contains "*", the entire original sequence is deleted from the SEQ field, and the entire original base sequence is written into the error message. When the CIGAR field contains "uM" or "u=", replace the corresponding u bases in the SEQ field with u first characters; When the CIGAR field contains "uD", replace the corresponding u bases in the SEQ field with u second characters, and append u "D" characters to the exception information; When the CIGAR field contains "uI", replace the corresponding u bases and the immediately upstream base in the SEQ field with a second character, and append an "I" character and the u original bases in the SEQ field at the current position to the error information; When the CIGAR field contains "uS" and appears at the beginning of the sequence, replace the corresponding u bases and the immediately downstream base in the SEQ field with a second character, and append an "S" character and the original u bases at the current position to the abnormal information; When the CIGAR field contains "uS" and appears at the end of the sequence, replace the corresponding u bases in the SEQ field and its immediate upstream base with a second character, and append an "S" character and the original u bases at the current position to the abnormal information; When the CIGAR field contains "uH", the corresponding u original bases are deleted from the SEQ field, and one "H" character is appended to the error message. When the CIGAR field contains "uX", replace the corresponding u bases in the SEQ field with u second characters, and append one "X" character and the u original bases at the current position to the error information; When the CIGAR field contains "uP", append u "P" characters to the exception information.

3. The method for compressing DNA third-generation sequencing data according to claim 2, characterized in that, Writing data from set X into matrix Y means writing the first and second characters of each sequence of data from set X into matrix Y. The specific operation includes the following steps: Create an ordered set X a Iterate through each site k of the reference genome sequentially, performing the following operations during each iteration: (a) Dynamically maintain the ordered set X a Add the sequence with starting position k from set X to the ordered set X. a At the end of set X, remove the sequence with ending position k from the ordered set X. a Delete; (b) Reverse prefix sorting and record index mapping: Calculate the ordered set X a The offset value between the start position and k of each sequence is used to extract set X. a The character whose index in each sequence is equal to the offset value is used, and then the ordered set X is sorted according to the lexicographical order of the characters. a Update the sorting to get X b And record the same sequence in set X a With X b Index mapping relationships in the database; (c) Matrix write operation: Traverse and update the sorted set X b During each traversal, the offset value of the next position after the starting position of the unterminated sequence is calculated: if the offset value is less than the length of the current sequence, the character with the offset value in the current sequence is written into matrix Y, and the row index is the current sequence in set X. a The index value in the sequence is k+1; if the offset value is greater than or equal to the current sequence length, the sequence is marked as terminated.

4. The method for compressing DNA third-generation sequencing data according to claim 3, characterized in that, The specific process of vertically reading matrix Y to obtain character stream file Z includes: Traverse each column of matrix Y according to the order of the starting alignment sites of each sequence, and read each row in each column from top to bottom: a: If the character read is the starting character of the sequence, then write the third character into the character stream Z and then write the currently read character; b: If the character read is the last character of the sequence, calculate the difference x between the row index of the current character and the row index of the previous ending sequence. If the difference is positive, write the currently read character into the character stream Z, then write x fourth characters and one sixth character; if the difference is negative, write the currently read character into the character stream Z, then write |x| fifth characters and one sixth character. c: If the character read is not the start or end character, then write the character directly to the character stream Z.

5. The method for compressing DNA third-generation sequencing data according to claim 4, characterized in that, The specific execution steps of the customized run-length encoding algorithm include: The character stream Z is fused and encoded into a run-length encoded string in pure integer sequence format according to the encoding rules, and then the alignment point of the first sequence in set X is added to the beginning of the run-length encoded string. The encoding rules include: If the third character is followed by the first character, the third character is encoded as "0"; otherwise, the third character is encoded as "-1". If the fourth character is followed by the first character, the fourth character is encoded as "-2"; otherwise, the fourth character is encoded as "-3". If the fifth character is followed by the first character, the fifth character is encoded as "-4"; otherwise, the fifth character is encoded as "-5". The first or second character is encoded as the number of times the character itself is repeated.

6. The method for compressing DNA third-generation sequencing data according to claim 1, characterized in that, The specific steps for converting the quality value to the WebP format image include: The quality values ​​are converted into integer data and then allocated to the RGB channels in sequence. Each set of integer data in the R, G, and B channels is encoded into an image pixel, and then an image in WebP format is generated.

7. A method for decompressing compressed files of third-generation DNA sequencing data, characterized in that, The decompression method is used to decompress the compressed file generated by the compression method according to any one of claims 1 to 6; If the file to be decompressed is a single file, decompress it directly; if the file to be decompressed consists of multiple block files, filter the block files in the target range and decompress each block file separately. Perform Huffman decoding on the header file and the compressed sequence data to obtain an integer string; decode the integer string using a custom run-length decoding method to obtain a character stream Z′; construct a matrix Y′ based on the character stream Z′; use reverse inverse prefix matching to operate on the matrix Y′ to generate a set X′ containing only the first and second character sequences; Perform Huffman decoding on the abnormal information data, and combine the set X′ with the decoded abnormal information to restore the original base sequence; Compressing quality values ​​involves parsing the RGB channel values ​​of each pixel in a WebP image into unsigned integers and then encoding them into a quality value string.

8. The method for decompressing compressed files of third-generation DNA sequencing data according to claim 7, characterized in that, The specific steps for decoding an integer string into a character stream Z′ include: Extract the first integer from the integer string as the starting position pos′ and write it to a temporary file; create a dynamically convertible character pointer P, and mark the character pointed to as the first character or the second character by controlling the parity of P; Read the integer string sequentially from the second element to the end, performing the following operations during the reading process: When a non-positive integer A is read: reset pointer P = A, then read the next integer B, and write B special characters corresponding to A according to the decoding rules; When a positive integer C is read: write the first or second character of the currently marked character to C pointers P, and increment pointer P; The decoding rules include: If the element is 0 or -1, decode the element as the third character; If the element is -2 or -3, decode the element into the fourth and sixth characters, respectively. If the element is -4 or -5, decode the element into the fifth and sixth characters.

9. A method for decompressing compressed files of third-generation DNA sequencing data according to claim 8, characterized in that, The step of constructing matrix Y′ based on character stream Z′ involves filling matrix Y′ vertically with the character stream Z′ and the starting position pos′ recorded in the temporary file. The specific steps for filling each column include: Initialize dynamic state variables: fill the state array full; counter for rows to be filled y; index of the end row δ; position pointer p′ = pos′; Parse the character stream Z' character by character and update the state variables: When encountering 'a' third characters, perform the following operations: Query the first 'y' elements in the array 'full' that have a value of false, and take the index L of the last element among the 'y'; get the number N of elements in the row corresponding to the index in matrix Y'; update the position pointer p' = p' + N, and record p' in a temporary file; update the value of 'a' to a + y, and the element at index L in the array 'full' has a value of true; When the fourth or fifth character is encountered, the value of 'a' is reduced by the number of consecutive occurrences of the current fourth or fifth character; each time a fourth character is encountered, the value of the integer δ is incremented by 1, and the value of the element with index δ in the array 'full' is set to true; each time a fifth character is encountered, the value of the integer δ is decremented by 1, and the value of the element with index δ in the array 'full' is set to true. When the first or second character is encountered, vertically fill the active row of the current column of matrix Y′.

10. A method for decompressing compressed files of third-generation DNA sequencing data according to claim 9, characterized in that, The specific steps for generating a set X′ containing only the first and second character sequences using the reverse inverse prefix matching operation matrix Y′ include: Create an intermediate matrix I with the same dimensions as Y′; read the starting positions of all sequences from the temporary file and sort them in ascending order; Sequentially traverse each site k of the reference genome, performing the following operations during each traversal: (a) Dynamically constructing an ordered set X a ′: Add the sequence with starting position k to the ordered set X. a At the end of ', remove the sequence with ending position equal to k from the ordered set X. a Delete in ' (b) For set X a Perform reverse prefix sorting and record the index mapping: Calculate the ordered set X a The offset value between the starting position and k of each sequence in ' is used to extract set X. a In the sequence X, the character whose index is equal to the offset value is used. Then, based on the lexicographical order of these characters, the ordered set X is... a Update the sorting to get set X b ′, and record the same sequence in set X a ′ and X b The index mapping relationship in '; (c) Reverse write operation: Traverse and update the sorted set X b During each iteration, the index mapping relationship and the index n of the current sequence are used to determine the position of the current sequence in set X. a The index m in matrix Y' is used to retrieve the character at column index k and row index m in matrix Y', and write it to the position at column index k and row index n in matrix I. After the site traversal of the reference genome is completed, matrix I is read row by row to obtain set X′.

Images