RU2807474C1 - Method for compressing genome sequence data - Google Patents

Abstract

FIELD: computer technology.

SUBSTANCE: invention relates to the field of computer technology, namely to a method for compressing genome sequence data. The method of compressing genomic sequence data includes obtaining a reading record by one or more computers, determining by one or more computers the correspondence of the reading record to the reading that is accurately mapped to the reference sequence or inaccurately mapped to the reference sequence, based on the determination by one or more computers that the reading record corresponds to the reading that is inaccurately mapped to the reference sequence, determination by one or more computers of how much the number of inconsistencies of an inaccurately matched reading satisfies a given threshold the number of nonconformities, and based on the determination that the number of nonconformities satisfies a given threshold number of nonconformities (i) obtaining by one or more computers an offset relative to the previous nonconformity that is less than the maximum encoded offset value, and (ii) encoding by one or more computers of each nonconformity of an inaccurately matched reading and an offset relative to the previous nonconformity in a record size 1 byte.

EFFECT: technical result is aimed at reducing data compression losses.

31 cl, 4 dwg,

Description

Field of technology

The technical field generally relates to methods for presenting genome sequencing data obtained by a sequencer, and in particular to computer-implemented methods for compressing such genome sequencing data. The description proposes a reference sequence-based compression method that provides fast compression and decompression while eliminating information loss, and which has a high compression ratio.

State of the art

Next-generation sequencers are now generating large volumes of sequencing data at an affordable price. Modern systems generate more than 6 billion 150-nucleotide sequences in a single 36-hour cycle, which is enough to sequence 20 complete human genomes. As a result, many new prospects are opening up for the diagnosis of genetic diseases or the development of personalized medicine, the goal of which is to tailor treatment taking into account the specifics of the human genome.

However, new challenges also arise, particularly related to the cost of storing large volumes of data. The most common file format for raw (unaligned) sequence data is the FASTQ format, which stores sequence data (a string of nucleotides A, C, T, G, also called a read), quality values (the sequencing platform's error probabilities in the sequence for each nucleotide) and sequence names. It is a plain ASCII text file that is often compressed using the standard LZ compression algorithm (the Lempel-Ziv algorithm implemented in gzip software). However, the use of such compression methods is associated with a number of problems:

- low compression ratio, since data repeatability is not fully used;

- slow compression and decompression.

There are also compression methods specifically designed for FASTQ encoding, which are divided into those that use a reference sequence and those that do not. However, none of them fully satisfies the requirements, since a) methods based on reference sequences provide optimal compression ratios, but they are slow, b) methods without reference sequences are faster, but have lower compression ratios. An example of such a reference-free method is the SPRING software, which is a reference-free compression program for FASTQ files (web address: github.com/shubhamchandak94/SPRING). However, the compression method offered by SPRING software has a low compression ratio.

A number of methods have been proposed, related to reference sequence compression methods, which use sequence alignments and should work faster with good compression ratios. However, such methods have a number of problems, with a particularly serious problem being that they do not completely eliminate losses. Such a known compression method based on a reference sequence is, for example, described in patent document WO 2018/068829 A1. In the described method, after alignment to one or more reference sequences, nucleotide sequences are classified according to degrees of match accuracy (thus creating classes of aligned reads), and then encoded into multiple layers of syntax elements using different source models and entropy encoders for each layer. into which data is divided. Thus, data classes are encoded separately and structured into different layers of syntax elements, with each layer including descriptors that uniquely reflect the classified and aligned reads of that layer. The method is designed to obtain clear information sources with reduced information entropy, thereby providing increased compression rates as well as selective access to specific classes of compressed data. However, this compression method results in the reordering of reads in an order that differs from that obtained at the end of the read alignment stage (i.e., the order of reads is changed in accordance with their classes). Therefore, during the compression process, certain information is lost, especially the order of the original sequence. Thus, this may affect the reproducibility of some analysis results, since some downstream analysis programs may be dependent on the order of the reads. In addition, unpacking data in an order that differs from the original reading order makes it very difficult to verify that the unpacked file is identical to the original file. Moreover, this compression method is comparatively slow, especially compared to existing compression methods in the art without a reference sequence.

Brief description of the invention

The following features of the independent claims solve the problem of existing prior art solutions by providing a method for compressing genome sequence data. In one aspect, a computer-implemented method for compressing genome sequence data obtained by a sequencer, wherein said genome sequence data includes nucleotide or base sequence reads that have been aligned to a reference sequence, thereby generating aligned reads, wherein said aligned reads are stored in form of a list of reads in the source file, and includes the stages at which:

- determine for each aligned read how accurately or inaccurately the specified read maps to the specified reference sequence, or whether the specified read does not map to the specified reference sequence,

- encoding the reads in accordance with the specified definition, wherein reads that are determined to be exactly matched are coded in accordance with the first coding process, and reads that are determined to be unmatched are coded in accordance with the second coding process,

- wherein the determination stage for each inaccurately mapped read includes a comparison of the number of mismatches between the specified read and the specified reference sequence, taking into account a threshold value,

- wherein at the coding stage, reads that are determined to be mismatched are coded in accordance with a second coding process or a third coding process, wherein mismatched reads are coded in accordance with the second coding process if the specified number of mismatches exceeds a threshold, and if the specified number of mismatches is less than the threshold value, then such inaccurately matched reads are encoded according to the third encoding process,

- wherein in said second encoding process, each nucleotide or reading base is encoded separately,

- wherein said first and third encoding processes include well-defined sets of descriptors, and each set of descriptors uniquely represents readings associated with the corresponding encoding process, and each of said first and third encoding processes represents an entropy-reducing encoding process of the information source.

The invention overcomes the disadvantages of prior art compression methods by allowing fast compression and decompression while eliminating information loss and providing a high compression ratio. More specifically, the invention focuses on encoding the most frequently occurring cases in the most compact manner, even if this means using reduced-performance coding modes for the rare, least frequently occurring cases. This provides significant improvement in compression performance. Moreover, due to the genomic information representation format used in the invention, the compression performed by the method according to the invention is faster. Last but not least, the method according to the invention allows the original order of the reads to be preserved as such and does not lead to a reordering of the reads according to their classes. Consequently, there is no loss of information during the process, facilitating subsequent analysis as well as efficient conformation checks after the decompression step.

These and other features and advantages of the present invention will become more apparent from the accompanying drawings and the following detailed description. In addition, while reference may be made herein to thresholds exceeded or not exceeded, it should be understood that such thresholds may be conceptually used to determine whether such threshold is satisfied, complied with, or otherwise recorded, regardless of whether the quantities or values used to evaluate such thresholds are described using positive or negative values.

In accordance with one innovative aspect of the present disclosure, a method for compressing genomic sequence data is disclosed. In one aspect, a method may include performing one or more operations by executing software instructions by one or more computers, the operations including obtaining a read record by one or more computers; determining by one or more computers that a read record matches a read that exactly maps to the reference sequence or does not exactly map to the reference sequence; based on the determination by one or more computers that a read entry matches such a read that does not accurately map to the reference sequence, determining by one or more computers how much the number of mismatches of the imprecisely matched read satisfies a predetermined threshold number of mismatches; and based on the determination that the number of mismatches satisfies a predetermined threshold number of mismatches, encoding by one or more computers each mismatch of the imprecisely mapped read into a 1-byte record.

Other aspects include appropriate systems, apparatus, and computer programs for performing the operations of the methods described herein, as defined in instructions encoded on computer-readable storage devices.

These and other versions may optionally include one or more of the following features. For example, in some embodiments, determining by one or more computers whether the number of mismatches of the loosely mapped read satisfies a predetermined threshold number of mismatches may include determining by one or more computers whether the number of mismatches of the loosely mapped reads exceeds the predetermined threshold number of mismatches.

In some embodiments, each read record may include data indicating the absolute position of the start of the aligned read relative to a reference sequence, data indicating an accurate or inaccurate read mapping, data indicating the length of the read, data indicating an accurate or inaccurate read mapping , data indicating the number of inconsistencies identified in the reading, and data indicating the relative position of each of the identified possible inconsistencies in the reading.

In some embodiments, encoding each mismatch of an imprecisely mapped read into a 1-byte record for each particular mismatch includes: encoding by one or more computers the first two bits of the byte to include data representing an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base in the reference sequence; and encoding by one or more computers the remaining six bits of the byte to include data representing the position of the mismatch in the reference sequence, said position being calculated as an offset relative to the previous mismatch read.

In some embodiments, the method may further include determining by one or more computers whether the offset exceeds a maximum encoded value, and based on the determination that the offset exceeds the maximum encoded value, inserting by one or more computers at least one dummy discrepancy between the particular discrepancy and previous discrepancy.

In some embodiments, based on a determination that the number of mismatches does not satisfy a predetermined threshold number of mismatches, one or more computers encode a list of reference sequence positions corresponding to the position of each of the mismatches relative to the reference sequence using an information entropy reduction encoding process.

In some embodiments, based on determining that the read record corresponds to a read that exactly maps to the reference sequence, the method may further include encoding by one or more computers at least a portion of the read record using information-reducing entropy coding.

In some embodiments, one or more computers may include one or more hardware processors.

In some embodiments, one or more hardware processors may include one or more field programmable gate arrays (FPGAs).

In some embodiments, a method for compressing genomic sequence data may be implemented using one or more hardware processors. In such embodiments, hardware processors may include a hardware processing system that is configured to perform one or more operations. In one aspect, the operations may include receiving a read record by the hardware processing system; determining by the hardware processing system whether the read record corresponds to a read that exactly maps to the reference sequence or does not exactly map to the reference sequence; based on the hardware processing system's determination of whether a read record corresponds to a read that does not accurately map to the reference sequence, determining by the hardware processing system how much the number of mismatches of the inaccurately matched read satisfies a predetermined threshold number of mismatches; and, based on the determination that the number of mismatches satisfies a predetermined threshold number of mismatches, encoding the hardware processing system of each mismatch of the imprecisely mapped read into a 1-byte record.

These and other features and advantages of the present invention will become more apparent from the accompanying drawings and the following detailed description.

In some embodiments, each read record may contain data indicating the absolute position of the start of the aligned read relative to a reference sequence, data indicating an accurate or imprecise read mapping, data indicating the length of the read, data indicating an accurate or inaccurate read mapping , data indicating the number of inconsistencies identified in the reading, and data indicating the relative position of these possible inconsistencies in the reading.

In some embodiments, determining by the hardware processing system whether the number of mismatches of the imprecisely mapped read satisfies a predetermined threshold number of mismatches may include determining by the hardware processing system whether the number of mismatches of the imprecisely mapped reads exceeds the predetermined threshold number of mismatches.

In some embodiments, encoding each mismatch of an imprecisely mapped read into a 1-byte record for each particular mismatch may include the hardware processing system encoding the first two bits of the byte to include data about an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base of the reference sequence ; and encoding the processing hardware of the six remaining bits of the byte to include information about the position of the mismatch in the reference sequence, wherein said position is calculated as an offset relative to the previous mismatch read.

In some embodiments, the hardware processing system is further configured to perform operations that include the hardware processing system determining whether the offset exceeds a maximum encoded value, and based on the determination that the offset exceeds the maximum encoded value, inserting at least one dummy discrepancies between a specific nonconformity and a previous nonconformity.

In some embodiments, based on a determination that the number of mismatches does not satisfy a predetermined threshold number of mismatches, the hardware processing system is further configured to perform operations that include encoding the hardware processing system of a list of reference sequence positions corresponding to the position of each of the mismatches relative to the reference sequence, with using a coding process that reduces the entropy of information.

In some embodiments, based on determining that the read record corresponds to a read that exactly maps to the reference sequence, the hardware processing system is further configured to perform operations that include the hardware processing system encoding at least a portion of the read record using reduced entropy coding information.

In some embodiments, the hardware processing system includes one or more field programmable gate arrays (FPGAs).

According to another innovative aspect of the present disclosure, a computer-implemented method of compressing genome sequence data obtained by a sequencer, said genome sequence data comprising nucleotide or base sequence reads that have been aligned to a reference sequence, thereby generating aligned reads, wherein the specified aligned reads are stored as a list of reads in the source file. In one aspect, the method may include acting on each aligned read, determining whether said read maps accurately or poorly to said reference sequence, or determining that said read does not map to said reference sequence; coding a read in accordance with said definition, wherein those reads that are classified as accurately mapped are coded in accordance with a first coding process, and those reads that are classified as inaccurately mapped are coded in accordance with a second coding process, wherein the determination step for of each inaccurately mapped read includes comparing the number of mismatches between said read and a specified reference sequence based on a threshold value, wherein at the encoding stage, those reads that are determined to be inaccurately mapped are encoded in accordance with a second encoding process or with a third encoding process, wherein inaccurately mapped reads are encoded according to a second encoding process if said number of mismatches is greater than a threshold, and if said number of mismatches is less than a threshold, then such loosely mapped reads are encoded according to a third encoding process, wherein in said second encoding process each nucleotide or the basis of the readings are encoded separately, wherein said first and third encoding processes include distinct sets of descriptors, and each set of descriptors uniquely represents the readings associated with the corresponding encoding process, and each of the first and third encoding processes is a down-coding process entropy of the information source.

Other aspects include appropriate systems, apparatus, and computer programs for performing the operations of the methods described herein, as defined in instructions encoded on computer-readable storage devices.

These and other versions may optionally include one or more of the following features. For example, in some embodiments, if a read is classified as inaccurately mapped to a reference sequence and its number of mismatches is less than a threshold, then the determination step may include determining how globally or locally the read maps to the specified reference sequence, and wherein the third encoding process includes a first additional encoding process and a second additional encoding process, wherein those reads that are determined to be globally mapped are encoded in accordance with the first additional encoding process, and those reads that are determined to be locally mapped are encoded in accordance with the second additional process, wherein said first and second additional encoding processes include well-defined sets of descriptors, with each set of descriptors uniquely representing the reads associated with the corresponding additional encoding process.

In some embodiments, said descriptors of said first additional encoding process may include an alignment start position in a reference sequence, a read length, and a list of character substitution mismatches, and wherein said descriptors of said second additional encoding process include a local alignment start position in a reference sequence, length readings, a list of inconsistencies in the form of character replacements, as well as the length of clipped sections that are not part of the alignment.

In some embodiments, at the encoding stage, the read cuts to be encoded in accordance with the second additional encoding process are concatenated, with each nucleotide or base of the read cuts being individually encoded.

In some implementations, during the encoding stage, each mismatch of an imprecisely mapped read is encoded as 1 byte.

In some embodiments, at the encoding stage, each mismatch of an imprecisely mapped read is encoded as follows: the first two bits of the byte are used to encode an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base of the reference sequence; and the last six bits of the byte are used to encode the position of the mismatch in the reference sequence, with the position calculated as an offset relative to the previous mismatch read.

In some embodiments, at the encoding stage, if the offset calculated between a particular disparity and a previous disparity is greater than the maximum encoded value, then at the encoding stage at least one dummy disparity is inserted between the two disparities until each offset between each of the specified mismatches and said at least one dummy mismatch will not be less than the specified maximum encodable value, wherein a dummy mismatch is defined as a mismatch for which the byte bits are used to encode the mismatch or to encode a nucleotide or base that is identical to the corresponding reference nucleotide or base in the reference sequence .

In some embodiments, the initial step of dividing the list of reads into blocks of reads, each block starting with a header containing information necessary to decode the block, the compression method being implemented on a block-by-block basis.

In some embodiments, read blocks are characterized by the same block size.

In some embodiments, the final step is to generate a compressed file containing a list of encoded reads, wherein said encoded reads are stored in the compressed file in the same order as the reads stored in the original file.

In some embodiments, this threshold value is 31.

In some embodiments, for each aligned read, the step of determining whether the read includes at least one mismatch corresponding to the case where the sequencer could not recognize any base or nucleotide is performed.

In some embodiments, for each read containing at least one mismatch corresponding to a case where the sequencer could not recognize any base or nucleotide, the step of determining the number of such mismatches and the step of comparing said number to a reference threshold is performed.

In some embodiments, at the encoding stage, if the number of such mismatches is greater than a reference threshold, each nucleotide or read base to be encoded according to the second encoding process is individually encoded into 4 bits, and if the number of such mismatches is less than the reference threshold, each nucleotide or base read to be encoded according to the second encoding process is individually encoded in 2 bits, and the encoding step further includes encoding a list of positions along the reference sequence, said positions corresponding to the positions of such mismatches in the reference sequence.

Brief description of drawings

In FIG. 1 is a flowchart illustrating the steps of a compression method in accordance with the invention.

In FIG. 2 is a diagram showing apparatus for implementing the steps of the compression method in accordance with the invention.

In FIG. Figure 3 shows the first example of a read that is globally mapped to a reference sequence.

In FIG. Figure 4 shows a second example of a read that is globally mapped to a reference sequence in a case where dummy mismatch insertion is required.

Detailed Description of the Invention

Genomic sequences provided herein include, by way of example and without limitation, nucleotide sequences, deoxyribonucleic acid (DNA) sequences, ribonucleic acid (RNA) sequences, and amino acid sequences. Although genomic information in the form of a nucleotide sequence has been described in some detail herein, it should be understood that the compression method of the invention can also be implemented for other genomic sequences, albeit with a number of modifications that will be apparent to one skilled in the art.

Genome sequencing information is obtained by sequencers in the form of nucleotide sequences (or, more generally, bases) represented as a string of letters from a defined set of terms. The smallest set of terms is represented by five characters: {A, C, G, T, N}, which represent the 4 types of nucleotides present in DNA, namely adenine, cytosine, guanine and thymine. In RNA, thymine is replaced by uracil (U). N indicates that the sequencer could not recognize any base, so the true nature of the position remains uncertain.

The nucleotide sequences produced by sequencers are called “reads.” Sequence reads can range in length from a few tens to several thousand nucleotides. Some technologies allow sequence reads to be obtained in pairs, with one read from one DNA strand and a second read from a different strand. Throughout this specification, the term “reference sequence” means any sequence against which the nucleotide or base sequences produced by the sequencer are aligned. One example of such a reference sequence could actually be a reference genome, i.e. a sequence constructed by scientists as a representative example of the set of genes of a particular species. However, the reference sequence may also consist of a synthetic sequence designed only to improve the compression ratio of the reads due to their subsequent processing. Sequencers may introduce errors in sequence reads, and in particular use the wrong symbol (ie, representing a different nucleic acid) to represent the nucleic acid or base actually present in the sample being sequenced; this is usually called a substitution error, or "mismatch."

The subject of the invention is a reference sequence-based compression method, the input data of which are nucleotide or base reads, and such reads have been previously aligned with respect to the reference sequence, thereby forming aligned reads. The aligned reads are then stored as a read list in the source file. The method of aligning reads and storing them in the source file after alignment is not critical to the invention and is not the purpose of this description. Each read is then encoded as a position on the reference sequence and as a list of differences with the specified reference sequence. Each read can then be reconstructed based on the encoded alignment information and the reference sequence using appropriate decompression software made in accordance with the present invention.

It is preferable that the alignment software that processes the reads and aligns them to a reference sequence before passing them as input to the compression software and apparatus does not take into account certain types of errors introduced into sequence reads, such as, for example, insertion errors or deletion errors. An insertion error occurs when one or more additional characters are inserted into a single sequence read that do not refer to any nucleic acid actually present. A deletion error occurs when one or more characters belonging to nucleic acids that are actually present in the sequenced sample are removed from a single sequence read. More precisely, in the case of an insertion error or deletion error in a given sequence read, the alignment software will treat the resulting erroneous nucleic acids as substitution errors, also called “mismatches.” This preferred choice for alignment software configuration allows for faster subsequent encoding, providing a significantly better trade-off between speed and compression ratio.

For each read, the alignment software transmits the corresponding read record to the compression software and machine. Each read record contains at least the following information: the absolute position of the start of the aligned read relative to the reference sequence, the length of the read, the type of read alignment, the number of mismatches identified in the read, and the relative position of these possible mismatches in the read (as applicable).

The compression method according to the present invention will be described with reference to FIG. 1. For example, the method may be performed by the apparatus 20 shown in FIG. 2. The apparatus includes at least one processor 22 and one storage device 24 operably coupled to the processor 22 to form a computing device. Storage device 24 may store computer program code or software 26 containing computer executable instructions that, when executed by processor 22, cause processor 22 to perform operations including steps of the compression method in accordance with the invention.

A source file containing aligned reads in the form of a read list, for example, is stored in the apparatus 20. As shown in FIG. 1, the method preferably includes an initial step 2, wherein the initial list of aligned reads is divided into blocks of reads. Typically, the list of aligned reads is divided into blocks of 50,000 reads, and this particular value should not be construed as limiting the scope of the present invention, which can be applied in the same way to other values. Preferably, read blocks are characterized by the same block size. Each block of reads begins with a header containing information necessary to decode the block, such as, for example, the size of the block's contents in bytes, and/or the identifier of the block or its contents, and/or the number of reads contained in the block. This provides support for compressed file concatenation as well as transfer capabilities (each block of reads contains all the information needed to decode the block's reads). In addition, since the compression method can be carried out block by block, it also allows for multi-threaded processing of blocks of reads, thus providing parallelization and, as a result, a certain gain in processing time. If all reads of a given block have the same length, the read length is also stored in the header; otherwise, when performing the compression method, a list of the length values of each read is explicitly stored.

Each read record contains information about the read alignment type. In general, there are two main types of alignment: precise alignment and imprecise alignment, with an additional type corresponding to "unmatched" alignment. "Imprecise alignment" means that the alignment contains at least one mismatch other than N, while at least a region of the read matches a region of the reference sequence (according to this definition, an imprecisely mapped read may contain one or more Ns, provided that it also contains one or more other inconsistencies). In one embodiment, each read record begins with the following bit flags, with each bit flag taking one of two possible values:

- the first bit flag indicates forward or reverse orientation relative to the reference sequence,

- the second bit flag indicates whether the alignment is accurate,

- the third bit flag indicates whether the read contains at least one N,

- the fourth bit flag indicates the bit depth of position data encoding - 16 bits or 32 bits.

Steps 4-12 below perform sequential processing of blocks of reads and sequential processing of reads within blocks.

The method includes the following step 4 of determining, for each aligned read, how closely or inaccurately the specified read maps to the reference sequence, or whether the specified read does not map to the reference sequence. Such determination step 4 for each mismatched read involves comparing 4A the number of mismatches between the specified read and the reference sequence based on a threshold value. In a preferred embodiment, which should not be construed as limiting the scope of the present invention, the specified threshold value is 31. This particular value was chosen specifically to provide the best possible compromise for storing the number of nonconformities in a fairly compact form, which will become more clear below in relation to the stage 12. Indeed, observational statistics show that in the vast majority of cases, inaccurately mapped reads involve fewer than 31 mismatches. The principle behind this choice is to code in the most compact form the most frequently occurring cases, leaving room for a number of very few low quality cases. If a read is determined to be inaccurately mapped with less than a threshold number of mismatches, determination step 4 involves further determining how globally or locally the read maps to the reference sequence. A “globally mapped read” is an imprecisely mapped read whose entire sequence, including the start and end of the read, is not accurately mapped to the reference sequence. A “locally mapped read” is an inaccurately mapped read containing a segment of nucleotides or bases that are inaccurately mapped to the reference sequence. Thus, the specified segment of nucleotides or bases corresponds to the region of the original read.

Preferably, for each aligned read, the method further includes the step 6 of determining whether said read includes at least one N, i.e. whether the specified read includes at least one mismatch corresponding to a case where the sequencer could not recognize any base or nucleotide. Then, for each read containing at least one N, the method includes the step 8 of determining the number of such N mismatches and the step 10 of comparing said number N of mismatches with a reference threshold value. In a preferred embodiment, which should not be construed as limiting the scope of the present invention, said reference threshold value is 31.

Regardless of the results of the definition step 4, the method includes the following step 12 of encoding the reads at least in accordance with the specified definition. More precisely, reads that are classified as exactly mapped to the reference sequence, regardless of whether they contain no Ns or the number of Ns in them is less than the reference threshold, are encoded according to the first encoding process. Reads that are classified as unmapped, or reads that are classified as exactly mapped but in which the number N is greater than the reference threshold, are encoded according to a second coding process, where each nucleotide or base is encoded separately, regardless of whether the specified one has been aligned nucleotide or base or not. Reads that are classified as inaccurately matched are coded according to the second coding process or the third coding process. More precisely, reads that are classified as inaccurately mapped with a number of mismatches greater than a threshold are coded according to the second coding process. If a read is determined to be inaccurately mapped to a number of mismatches less than a threshold, and if said read does not contain any Ns or the number of Ns in such a read is less than the reference threshold, then said read is encoded according to the third encoding process. Otherwise, i.e. if the number N in such a read is greater than the reference threshold, then the specified read is encoded in accordance with the second encoding process.

Regardless of whether a particular read is classified as accurately mapped, inaccurately mapped, or unmapped, if the specified read contains at least one N but the number of Ns is less than the reference threshold, then encoding step 12 includes encoding a list of positions along the reference sequence, wherein positions correspond to N positions in the reference sequence. The list of positions is then stored in a storage device of the computing device, wherein said device implements the compression method. If a read contains at least one N value, but the number of N is less than the reference threshold, and it is to be encoded by a second encoding process, then each nucleotide or base of the read is individually encoded with 2 bits.

If a read contains at least one N value, but the number of N is greater than the reference threshold, then the said read is in all cases encoded by the second encoding process, and each nucleotide or base of the read is individually encoded with 4 bits. In this case, the encoding step 12 does not include encoding and storing the list of positions N in the reference sequence. Indeed, each N mismatch is then directly encoded according to a second encoding process in exactly the same way as other nucleotides or base reads.

The first and third encoding processes involve clearly defined sets of descriptors. Each set of descriptors uniquely represents the readings associated with the corresponding encoding process, and each of the first and third encoding processes represents an information entropy decreasing encoding process. More precisely, the third encoding process includes a first additional encoding process and a second additional encoding process. Inaccurately mapped reads, which are classified as globally mapped in step 4, are encoded according to the first additional encoding process. Inaccurately mapped reads that are classified as locally mapped in step 4 are encoded according to a second additional encoding process. The first and second additional encoding processes include well-defined sets of descriptors that uniquely represent the reads associated with the corresponding additional encoding process.

Thus, the alignment information that is encoded for each read and which allows the complete sequence of the read to be recovered during data decompression depends on the corresponding process or additional encoding process that was used for said read. For example, descriptors used in the first encoding process could be:

the absolute position of the start of the exactly mapped read relative to the reference sequence (in 16- or 32-bit encoding), and

read length (coded using differential encoding relative to the length of the previous read, with a variable length code ranging from 2 bits to 34 bits).

The descriptors used in the first additional encoding process may be:

absolute position of the start of the loosely mapped read relative to the reference sequence (in 16- or 32-bit encoding),

read length (coded using differential encoding relative to the length of the previous read with a variable length code ranging from 2 bits to 34 bits) and

list of reading inconsistencies.

Descriptors used in the second additional encoding process may be:

the absolute position of the start of the loosely mapped read region relative to the reference sequence, also called the start position of the local alignment (in 16- or 32-bit encoding),

read length (encoded using differential encoding relative to the length of the previous read with a variable length code ranging from 2 bits to 34 bits),

list of reading inconsistencies and

length of read trims that are not part of the alignment (8-bit encoding for each trim).

Preferably, the mismatch list that is encoded in the first and second additional processes includes a header (a bit flag encoded by 1 byte). The first five bits of the byte are used to encode the number of mismatches contained in the read (in the preferred embodiment where the threshold is 31, the number is in the range [0-31]). One bit is then used to encode information about whether the loosely mapped read is globally or locally mapped. Another bit is used to encode information about whether the 2-bit mode is enabled for the second encoding process. The last bit is used to encode information about whether 4-bit mode is enabled for the second encoding process. Preferably, for each read encoded in accordance with the second additional encoding process in encoding step 12, concatenation of the truncated regions of said read (i.e., those regions that are not part of the local alignment) is performed and each nucleotide or base of said truncated regions is encoded separately . In a preferred embodiment, each nucleotide or base of such truncated read regions is encoded separately by 2 bits.

In a preferred embodiment, each mismatch encoded in the mismatch list of the imprecisely mapped read (ie, encoded according to the first or second additional encoding process) is encoded as 1 byte. More specifically, each inaccurately mapped read mismatch that is to be coded according to the first or second additional coding process may be coded as follows:

the first two bits of the byte are used to encode an alternative nucleotide or base present in the read instead of the corresponding reference nucleotide or base of the reference sequence,

the last six bits are used to encode the position of the mismatch in the reference sequence, with the specified position being calculated as an offset relative to the previous read mismatch (the relative position of the mismatch, excluding the first read mismatch for which the absolute position is encoded). Thus, the range of such an offset, which is encoded by 6 bits, is [0-63].

In FIG. 3 shows an example of coding read discrepancies using the first additional coding process. A read is a loosely mapped read that is globally mapped to a reference sequence. The reading contains two inconsistencies:

the first mismatch, located at position 12 of the read, which is a replacement of nucleotide A in the reference sequence by nucleotide T in the read, and

the second mismatch, located at position 21 of the read, which is a replacement of nucleotide C in the reference sequence with nucleotide G in the read.

The list of nonconformities is then encoded as follows:

<12, T>, the value "12" corresponding to the absolute position of the first mismatch in the reading, and

<9, G>, the value "9" corresponding to the relative position of the second mismatch in the reading, i.e. the offset between the second discrepancy and the first discrepancy.

<12, T>can, for example, be converted to the value "51" (encoded by 1 byte), and <9, G>can be converted to the value "38" (encoded by 1 byte). This byte encoding is obtained as follows:

offset position x 4+ nucleotide code (at A=0, C=1, G=2, T=3)

Preferably, for each loosely mapped read that is to be encoded in accordance with the first or second additional encoding process, if the offset calculated between a given read mismatch and a previous read mismatch is greater than the maximum encoded value, then at least one “dummy” mismatch is inserted between those two discrepancies until each offset between each of the specified discrepancies and at least one “dummy” discrepancy is less than the specified maximum encoded value. A “dummy” mismatch is defined as a mismatch for which the byte bits used to encode the mismatch encode a nucleotide or base that is equivalent to the corresponding reference nucleotide or base in the reference sequence. In a preferred embodiment, which should not be construed as limiting the scope of the present invention, the maximum encoded value is 63, which corresponds to the maximum value that can be encoded in 6 bits.

In FIG. 4 shows an example of encoding read mismatches using a first additional encoding process in the case where insertion of a “dummy” mismatch is required. A read is a loosely mapped read that is globally mapped to a reference sequence. The reading contains two inconsistencies:

the first mismatch, located at position 22 of the read, which is a replacement of nucleotide A in the reference sequence by nucleotide T in the read, and

the second mismatch, located at position 134 of the read, which is a replacement of nucleotide C in the reference sequence with nucleotide G in the read.

The offset position between the second and first mismatch is 112, which is greater than the maximum coded value of 63. Thus, a "dummy" mismatch must be inserted between the two mismatches such that each offset between each of the mismatches and the "dummy" mismatch is less than that specified maximum encoded value. A “fictitious” T nucleotide mismatch (corresponding to a “real” T nucleotide in the reference sequence), for example, is introduced at read position 85. The offset position calculated between the "dummy" disparity and the first disparity is 63, which is the same as the maximum encoded value. The offset position calculated between the second disparity and the "dummy" disparity is 49, which is less than 63.

The list of nonconformities is then encoded as follows:

<22, T>, the value "22" corresponding to the absolute position of the first mismatch in the reading,

<63, T>, the value “63”, corresponding to the relative position of the “fictitious” discrepancy in the reading, i.e. the offset between the "dummy" mismatch and the first mismatch, and

<49, T>, the value "49" corresponding to the relative position of the second mismatch in the reading, i.e. offset between the second discrepancy and the “dummy” discrepancy.

For example, <22, T> can be converted to the value "91" (encoded by 1 byte), <63, T> can be converted to the value "255" (encoded by 1 byte), and <49, G> can be converted to value "198" (encoded by 1 byte). This byte encoding is obtained as follows:

offset position x 4+ nucleotide code (at A=0, C=1, G=2, T=3)

The method includes a final step 14 of generating a compressed file containing a list of encoded reads. The encoded reads are stored in the compressed file in the same order as the reads stored in the uncompressed original file. Each read can then be reconstructed based on the encoded alignment information and the reference sequence using appropriate decompression software and/or a method implemented in accordance with the present invention.

Although the inventive techniques disclosed herein are described with reference to an example architecture of computing device 20 (shown in FIG. 2 for illustrative purposes), they may be implemented in hardware, software, firmware, or any their combinations. When implemented in software, the computer program code may be stored on a computer-readable medium or executed by a hardware processing unit including one or more processors, as is the case with device 20 in FIG. 2. It should be understood that the term "processor" as used herein includes one or more processing devices, including a signal processor, microprocessor, microcontroller, application specific integrated circuit (ASIC), field programmable gate array (FPGA), or other type of processing system, as well as sections or combinations of such system elements. Also, the term “storage device” as used herein includes an electronic storage device associated with a processor, such as random access memory (RAM), read-only memory (ROM), or other types of storage devices, in any combination.

Accordingly, software instructions or code for implementing the methodologies and protocols described herein may be stored in one or more associated storage devices, such as RAM, embedded or removable storage, and then, when ready for use, loaded into RAM and executed by the processor.

The techniques of the present disclosure may be implemented in a wide variety of devices or apparatus, including, for example, mobile phones, computers, servers, tablet computers, and similar devices.

Although illustrative embodiments of the present invention have been described herein with reference to the accompanying drawings, it should be understood that the invention is not limited to these embodiments and that various other changes and modifications can be made by one skilled in the art without departing beyond the scope or spirit of the invention.

Statistical and numerical examples of the compression method according to the invention

The comparison example below was implemented for an uncompressed data file that contained 48 million reads or nucleotide sequences:

uncompressed data file size: 35,770 MB (megabytes)

file size compressed using gzip software: 6649 MB

file size compressed using software without SPRING reference sequence: 1402 MB

file size compressed using the reference sequence-based compression method according to the present invention: 1179 MB

compression time using software without SPRING reference sequence: 1722 s

compression time using the reference sequence-based compression method according to the present invention: 181 s

average size of bits/nucleotides in a data file without compression (ASCII encoding): 8 bits/nucleotide; average size of bits/nucleotides in a file that has been compressed with an encoding adapted for the 4 possible characters A, T, C, G: 2 bits/ nucleotide

average bits/nucleotide size in a file that has been compressed using the reference sequence compression method of the present invention: 0.33 bits/nucleotide

The above numerical examples show that the present invention achieves fast compression and decompression, providing a high compression ratio.

Claims (94)

1. A method for compressing genomic sequence data, comprising: one or more computers receiving a read record; determining by one or more computers that the record of reads corresponds to a read that exactly maps to the reference sequence or does not exactly map to the reference sequence; based on the determination by one or more computers that the read record corresponds to a read that is not accurately mapped to the reference sequence, determining by one or more computers how much the number of mismatches of the inaccurately mapped read satisfies a predetermined threshold number of mismatches; And based on determining that the number of mismatches satisfies a predetermined threshold number of mismatches, (i) one or more computers obtaining an offset from a previous mismatch that is less than the maximum encoded offset value, and (ii) one or more computers encoding each mismatch of the imprecisely matched read and offset relative to a previous read-to-write mismatch of 1 byte. 2. The method of claim 1, wherein determining by one or more computers how much the number of mismatches of the imprecisely mapped read satisfies a predetermined threshold number of mismatches includes determining by one or more computers whether the number of mismatches of the imprecisely mapped read exceeds a predetermined threshold number of mismatches. 3. The method according to claim 1, in which each read record contains: data indicating the absolute position of the start of the aligned read relative to the reference sequence; data indicating read length; data indicating accurate or inaccurate read mapping; data indicating the number of inconsistencies identified in the reading; And data indicating the relative position of each of the indicated possible discrepancies in the reading. 4. The method of claim 1, wherein encoding each imprecisely mapped read mismatch into a 1-byte record for each specific mismatch includes: encoding by one or more computers the first two bits of a byte to include data reflecting an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base in the reference sequence; And the encoding by one or more computers of the last six bits of a byte to include data indicating the offset. 5. The method according to claim 4, additionally including: determining by one or more computers whether the offset exceeds a maximum encoded value; And based on the determination that the offset exceeds the maximum encoded value, insertion by one or more computers of at least one fictitious discrepancy between the particular discrepancy and the previous discrepancy. 6. The method according to claim 1, additionally including based on the determination that the number of mismatches does not satisfy a predetermined threshold number of mismatches, encoding by one or more computers a list of reference sequence positions corresponding to the position of each of the mismatches relative to the reference sequence using an information entropy reduction encoding process. 7. The method according to claim 1, further comprising based on the determination that the read record corresponds to a read that exactly maps to the reference sequence, encoding by one or more computers at least a portion of the read record using an information-reducing entropy coding process. 8. The method of claim 1, wherein the one or more computers comprise one or more hardware processors. 9. The method of claim 8, wherein the one or more hardware processors comprise one or more field programmable gate arrays (FPGAs). 10. A hardware processor that includes a hardware processing system that is configured to perform one or more operations, wherein the one or more operations include: receiving a reading record by a hardware processing system; determining by the hardware processing system that the read entry corresponds to a read that maps accurately to the reference sequence or that maps inaccurately to the reference sequence; based on the determination by the hardware processing system that the read record corresponds to a read that is not accurately mapped to the reference sequence, determining by one or more computers how much the number of mismatches of the inaccurately mapped read satisfies a predetermined threshold number of mismatches; And based on determining that the number of mismatches satisfies a predetermined threshold number of mismatches, (i) the hardware processing system obtains an offset from the previous mismatch that is less than the maximum encoded offset value, and (ii) the hardware processing system encodes each mismatch of the inaccurately matched read and the offset from the previous 1-byte read-to-write mismatch. 11. The hardware processor according to claim 10, wherein each read record contains: data indicating the absolute position of the start of the aligned read relative to the reference sequence; data indicating read length; data indicating accurate or inaccurate read mapping; data indicating the number of inconsistencies identified in the reading; And data indicating the relative position of specified possible discrepancies in reading. 12. The hardware processor of claim 10, wherein encoding each imprecisely matched read mismatch into a 1-byte record for each specific mismatch includes: encoding the hardware system to process the first two bits of a byte to include data reflecting an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base in the reference sequence; And the hardware system's encoding of the last six bits of a byte to include data that represents the offset. 13. The hardware processor of claim 12, wherein the hardware processing system is further configured to perform operations including: determining by the hardware processing system whether the offset exceeds a maximum encoded value; based on the determination that the offset exceeds the maximum encoded value, the hardware processing system inserting at least one dummy mismatch between the particular mismatch and the previous mismatch. 14. The hardware processor of claim 10, wherein the hardware processing system is further configured to perform operations including: based on determining that the number of discrepancies does not satisfy a predetermined threshold number of discrepancies, encoding by the hardware processing system a list of reference sequence positions corresponding to the position of each of the discrepancies relative to the reference sequence using an information entropy reduction encoding process. 15. The hardware processor of claim 10, wherein the hardware processing system is further configured to perform operations including based on the determination that the read record corresponds to a read that exactly maps to the reference sequence, encoding at least a portion of the read record by the hardware processing system using information entropy reduction coding. 16. The hardware processor of claim 10, wherein the hardware processing system comprises one or more field programmable gate arrays (FPGAs). 17. The hardware processor of claim 10, wherein determining by the hardware processing system how much the number of mismatches of the imprecisely mapped read satisfies a predetermined threshold number of mismatches includes determining by the hardware processing system whether the number of mismatches of the inaccurately mapped read exceeds a predetermined threshold number of mismatches. 18. Genome sequence data compression system containing: one or more computers and one or more storage devices that store instructions that, when executed by one or more computers, perform operations including: one or more computers receiving a read record; determining by one or more computers that the record of reads corresponds to a read that exactly maps to the reference sequence or does not exactly map to the reference sequence; based on the determination by one or more computers that the read entry corresponds to a read that is not accurately mapped to the reference sequence, determining by one or more computers how much the number of mismatches of the inaccurately mapped read satisfies a predetermined threshold number of mismatches; And based on determining that the number of mismatches satisfies a predetermined threshold number of mismatches, (i) one or more computers obtaining an offset from a previous mismatch that is less than the maximum encoded offset value, and (ii) one or more computers encoding each mismatch of the imprecisely matched read and offset relative to a previous read-to-write mismatch of 1 byte. 19. The system according to claim 18, in which each read record contains: data indicating the absolute position of the start of the aligned read relative to the reference sequence; data indicating read length; data indicating accurate or inaccurate read mapping; data indicating the number of inconsistencies identified in the reading; And data indicating the relative position of each of the indicated possible discrepancies in the reading. 20. The system of claim 18, wherein encoding each imprecisely mapped read mismatch into a 1-byte record includes, for each specific mismatch: encoding by one or more computers the first two bits of a byte to include data reflecting an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base in the reference sequence; And the encoding by one or more computers of the last six bits of a byte to include data indicating the offset. 21. The system according to claim 20, in which the operations additionally include: determining by one or more computers whether the offset exceeds a maximum encoded value; And based on the determination that the offset exceeds the maximum encoded value, insertion by one or more computers of at least one fictitious discrepancy between the particular discrepancy and the previous discrepancy. 22. The system according to claim 18, in which the operations further include based on the determination that the number of mismatches does not satisfy a predetermined threshold number of mismatches, encoding by one or more computers a list of reference sequence positions corresponding to the position of each of the mismatches relative to the reference sequence using an information entropy reduction encoding process. 23. The system according to claim 18, in which the operations further include based on the determination that the read record corresponds to a read that exactly maps to the reference sequence, encoding by one or more computers at least a portion of the read record using entropy-reducing coding of the information. 24. The system of claim 18, wherein the one or more computers comprise one or more hardware processors. 25. The system of claim 24, wherein the one or more hardware processors comprise one or more field programmable gate arrays (FPGAs). 26. A machine-readable storage device having instructions stored therein that, when executed by a data processing apparatus, cause the data processing apparatus to perform genomic sequence data compression operations comprising: receiving a read record; determining whether a read record matches a read that maps accurately to the reference sequence or that maps inaccurately to the reference sequence; based on determining that the read record corresponds to a read that is not accurately mapped to the reference sequence, determining how much the number of mismatches of the inaccurately mapped read satisfies a predetermined threshold number of mismatches; And based on determining that the number of mismatches satisfies a specified threshold number of mismatches, (i) obtaining an offset from the previous mismatch that is less than the maximum encoded offset value, and (ii) encoding each mismatch of the imprecisely mapped read and the offset from the previous mismatch into a record of size 1 byte. 27. The machine-readable storage device of claim 26, wherein each read record contains: data indicating the absolute position of the start of the aligned read relative to the reference sequence; data indicating read length; data indicating accurate or inaccurate read mapping; data indicating the number of inconsistencies identified in the reading; And data indicating the relative position of each of the indicated possible discrepancies in the reading. 28. The machine-readable storage device of claim 26, wherein encoding each imprecisely matched read mismatch into a 1-byte record for each particular mismatch includes: encoding by one or more computers the first two bits of a byte to include data reflecting an alternative nucleotide or base present in the read in place of the corresponding reference nucleotide or base in the reference sequence; And the encoding by one or more computers of the last six bits of a byte to include data indicating the offset. 29. The machine-readable storage device of claim 28, wherein the operations further include: determining by one or more computers whether the offset exceeds a maximum encoded value; And based on the determination that the offset exceeds the maximum encoded value, insertion by one or more computers of at least one fictitious discrepancy between the particular discrepancy and the previous discrepancy. 30. The machine-readable storage device of claim 26, wherein the operations further include: based on determining that the number of mismatches does not satisfy a predetermined threshold number of mismatches, encoding a list of reference sequence positions corresponding to the position of each of the mismatches relative to the reference sequence using an information entropy reduction encoding process. 31. The machine-readable storage device of claim 26, wherein the steps further include based on the determination that the read record corresponds to a read that exactly maps to the reference sequence, encoding at least a portion of the read record using entropy-reducing coding of the information.

Images