Method for compressing genomic sequence data

Abstract

A method for compressing genomic sequence data is disclosed. The present invention relates to a reference-based method for compressing genomic sequence data produced by a sequencing machine. Determining whether a sequence of nucleotides or bases previously aligned with the reference sequence is completely mapped, incompletely mapped or not mapped with the reference sequence; coding is then performed according to the determination. The determining step comprises: for each incompletely mapped sequence, comparing the number of mismatches between said sequence and a reference sequence with a reference threshold, and comparing the number of mismatches with the reference threshold as a result of said comparison method for compressing genomic sequence data generated by a sequencing machine. And coding the incompletely mapped sequence according to different coding processes.

Description

[0001] This application is a divisional application of the patent application with application number 202080062683.7 filed on March 7, 2022 and invention name “Method for compressing genomic sequence data”. Technical Field

[0002] The present invention generally relates to methods for representing genomic sequencing data generated by sequencing machines, and more particularly to computer-implemented methods for compressing such genomic sequencing data. The present disclosure provides a reference-based compression method that allows for fast compression and decompression without causing information loss and with a high compression ratio. Background Art

[0003] Next-generation sequencing machines now generate massive amounts of sequencing data at an affordable price. The latest systems produce over 6 billion 150-nucleotide sequences in a single 36-hour run, enough to sequence 20 complete human genomes. This opens up many new perspectives for the diagnosis of genetic diseases and the development of personalized medicine, designed to adapt treatments based on the specificity of the human genome.

[0004] However, this also brings new challenges, particularly the costs associated with storing massive amounts of data. The most commonly used file format for raw (unaligned) sequence data is the FASTQ format, which stores sequence data (strings of A, C, T, and G nucleotides, also known as reads), quality values (the probability that the sequencing platform will make a sequencing error for each nucleotide), and sequence names. This is a plain ASCII text file, typically compressed using the common text compression scheme LZ (Lempel-Ziv scheme, implemented in gzip software). However, using this type of compression method introduces several problems:

[0005] -Low compression ratio due to incomplete use of data redundancy

[0006] - Slow compression and decompression

[0007] There are also compression methods specifically designed for FASTQ encoding, which are categorized as reference-based or non-reference-based. However, none of these methods are entirely satisfactory because: a) reference-based methods have good compression ratios but are slow, and b) non-reference-based methods are faster but have lower compression ratios. An example of such a non-reference-based method is provided by the software SPRING, a reference-free compressor for FASTQ files (World Wide Web address: github.com / shubhamchandak94 / SPRING). However, the compression method provided by the SPRING software has a low compression ratio.

[0008] Among reference-based compression methods, some have been proposed that use sequence alignments and aim to be faster and have good compression ratios. However, such methods encounter several problems, a major one being that they are not completely lossless. This known reference-based compression method is described, for example, in patent document WO 2018 / 068829 A1. In the described method, after alignment with one or more reference sequences, nucleotide sequences are classified according to the accuracy of the match (thereby creating categories of aligned reads), and then these nucleotide sequences are encoded into multi-layer syntax elements using different source models and entropy encoders for each layer in which the data is partitioned. Thus, the categories of data are encoded separately and constructed in different layers of syntax elements, each layer containing descriptors that univocally represent the classified and aligned reads of the layer. This method aims to obtain different information sources with simplified information entropy, thereby allowing improved compression performance and selective access to compressed data of specific categories. However, this compression method reorders the reads in an order different from the order obtained at the end of the read alignment step (i.e., the reads are reordered according to their categories). Some information is lost during the compression process (especially the initial sequence ordering). Consequently, the reproducibility of some analytical results may be affected, as some downstream analysis software may rely on the order of the reads. Furthermore, decompressing the data in an order different from the original order of the reads makes it more difficult to verify that the uncompressed file is identical to the original file. Furthermore, this compression method is relatively slow, especially when compared to state-of-the-art non-reference-based compression methods. Summary of the Invention

[0009] The features of the following independent claims solve the problems of the prior art solutions by providing a method for compressing genomic sequence data. In one aspect, a computer-implemented method for compressing genomic sequence data generated by a sequencing machine is disclosed, the genomic sequence data comprising reads of a sequence of nucleotides or bases that have been aligned to a reference sequence, thereby generating aligned reads, the aligned reads being stored as a read list in an initial file, the method comprising the following steps:

[0010] - for each aligned read, determining whether the read is completely or incompletely mapped to the reference sequence, or whether the read is unmapped to the reference sequence,

[0011] - encoding the reads according to the determination, wherein the reads determined to be fully mapped are encoded according to a first encoding process and the reads determined to be unmapped are encoded according to a second encoding process,

[0012] - wherein the determining step comprises, for each incompletely mapped read, comparing the number of mismatches between the read and the reference sequence to a threshold value,

[0013] - wherein, in the encoding step, the read segments determined to be incompletely mapped are encoded according to the second encoding process or the third encoding process; when the number of mismatches is greater than the threshold, the incompletely mapped read segments are encoded according to the second encoding process; and when the number of mismatches is less than the threshold, the incompletely mapped read segments are encoded according to the third encoding process;

[0014] - wherein, during the second encoding process, each nucleotide or base of the read is encoded individually,

[0015] - wherein the first encoding process and the third encoding process include different sets of descriptors, each descriptor set univocally representing a read segment associated with the corresponding encoding process, each of the first encoding process and the third encoding process being a simplified information source entropy encoding process.

[0016] The present invention overcomes the shortcomings of existing compression methods by allowing fast compression and decompression without causing information loss and providing a high compression ratio. More specifically, the focus of the present invention is to encode the most frequently occurring cases in the most compact way, even if this means using a degraded encoding mode for the rare least frequently occurring cases. This leads to a huge improvement in compression performance. In addition, due to the genomic information representation format used in the present invention, the compression performed by the method according to the present invention is faster. Last but not least, the method according to the present invention maintains the original order of the reads as such and does not reorder the reads according to their category. Therefore, no information is lost during the process, which makes it easier to perform downstream analysis and to perform effective consistency checks after the decompression step.

[0017] These and other features and advantages of the present invention will become more apparent with reference to the accompanying drawings and the following detailed description. In addition, although threshold values may be referred to herein as being exceeded or not exceeded, it should be understood that such threshold values may be conceptually employed such that it is determined whether such threshold values are met, complied with, or otherwise detected, regardless of whether the numbers or values used to achieve those threshold evaluations are described using positive or negative values.

[0018] According to an innovative aspect of the present disclosure, a method for compressing genomic sequence data is disclosed. In one aspect, the method may include performing one or more operations by executing software instructions through one or more computers, wherein the operations include: obtaining a read record by the one or more computers; determining by the one or more computers whether the read record corresponds to a read that is completely mapped to a reference sequence or a read that is incompletely mapped to the reference sequence; based on the determination by the one or more computers that the read record corresponds to a read that is incompletely mapped to the reference sequence, determining by the one or more computers whether the number of mismatches of the incompletely mapped reads meets a predetermined mismatch threshold number, and based on determining that the number of mismatches meets the predetermined mismatch threshold number, encoding by the one or more computers each mismatch of the incompletely mapped reads as a record having a size of 1 byte.

[0019] Other aspects include corresponding systems, apparatus, and computer programs that perform the actions of the methods as disclosed herein, as defined by instructions encoded on a computer-readable storage device.

[0020] These and other versions may optionally include one or more of the following features. For example, in some implementations, determining, by the one or more computers, whether the number of mismatches of the incompletely mapped reads satisfies a predetermined threshold number of mismatches may include determining, by the one or more computers, whether the number of mismatches of the incompletely mapped reads is greater than the predetermined threshold number of mismatches.

[0021] In some embodiments, each read record may include: data indicating the absolute starting position of the aligned read relative to the reference sequence; data indicating the length of the read; data indicating whether the read is fully mapped or incompletely mapped; data indicating the number of mismatches identified in the read; and data indicating the relative position of each of the possible mismatches in the read.

[0022] In some embodiments, encoding each mismatch of the incompletely mapped read segment as a record having a size of 1 byte includes: for each specific mismatch, encoding, by the one or more computers, the first two bits of the byte to include data representing an alternative nucleotide or base present in the read segment instead of the corresponding reference nucleotide or base in the reference sequence; and encoding, by the one or more computers, the remaining six bits of the byte to include data representing the position of the mismatch in the reference sequence, the position being calculated as an offset relative to a previous mismatch of the read segment.

[0023] In some specific implementations, the method may also include: determining, by one or more computers, whether the offset is greater than a maximum encodable value; and inserting, by one or more computers, at least one false mismatch between the particular mismatch and the previous mismatch based on determining that the offset is greater than the maximum encodable value.

[0024] In some specific implementations, the method may also include: based on determining that the number of mismatches does not meet the predetermined mismatch threshold number, one or more computers use a simplified information entropy coding process to encode a list of positions of the reference sequence corresponding to the position of each of the mismatches into the reference sequence.

[0025] In some implementations, the method may further include encoding, by one or more computers, at least a portion of the read record using simplified entropy coding based on determining that the read record corresponds to a read that is fully mapped to the reference sequence.

[0026] In some implementations, the one or more computers can include one or more hardware processors.

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

[0028] In some embodiments, the method for compressing genomic sequence data may be performed by one or more hardware processors. In such embodiments, the hardware processor may include a hardware processing circuit system configured to perform one or more operations. In one aspect, the operations may include: obtaining a read record by the hardware processing circuit system; determining by the hardware processing circuit system whether the read record corresponds to a read that is completely mapped to a reference sequence or a read that is incompletely mapped to the reference sequence; based on the determination by the hardware processing circuit system that the read record corresponds to a read that is incompletely mapped to the reference sequence, determining by the one or more computers whether the number of mismatches of the incompletely mapped reads meets a predetermined mismatch threshold number, and based on determining that the number of mismatches meets the predetermined mismatch threshold number, encoding by the hardware processing circuit system each mismatch of the incompletely mapped reads as a record having a size of 1 byte.

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

[0030] In some embodiments, each read record may include: data indicating the absolute starting position of the aligned read relative to the reference sequence; data indicating the length of the read; data indicating whether the read is fully mapped or incompletely mapped; data indicating the number of mismatches identified in the read; and data indicating the relative position of the possible mismatches in the read.

[0031] In some specific implementations, determining by the hardware processing circuit system whether the number of mismatches of the incompletely mapped read segments meets a predetermined mismatch threshold number may include determining by the hardware processing circuit system whether the number of mismatches of the incompletely mapped read segments is greater than the predetermined mismatch threshold number.

[0032] In some embodiments, encoding each mismatch of the incompletely mapped read segment as a record having a size of 1 byte may include: for each specific mismatch, encoding by the hardware processing circuit system the first two bits of the byte to include data representing an alternative nucleotide or base present in the read segment instead of the corresponding reference nucleotide or base in the reference sequence; and encoding by the hardware processing circuit system the remaining six bits of the byte to include data representing the position of the mismatch in the reference sequence, the position being calculated as an offset relative to a previous mismatch of the read segment.

[0033] In some specific implementations, the hardware processor circuit system is further configured to perform the following operations, which include: determining by the hardware processing circuit system whether the offset is greater than a maximum encodable value; and based on determining that the offset is greater than the maximum encodable value, inserting by the hardware processing circuit system at least one false mismatch between the specific mismatch and the previous mismatch.

[0034] In some specific implementations, the hardware processor circuit system is further configured to perform the following operations, which include: based on determining that the number of mismatches does not meet the predetermined mismatch threshold number, encoding, by the hardware processing circuit system, a list of positions of the reference sequence corresponding to the position of each of the mismatches into the reference sequence using a simplified information entropy encoding process.

[0035] In some specific implementations, the hardware processor circuit system is further configured to perform the following operations, which include: based on determining that the read segment record corresponds to a read segment that is completely mapped to the reference sequence, encoding at least a portion of the read segment record using simplified information entropy coding by the hardware processing circuit system.

[0036] In some implementations, the hardware processing circuitry includes one or more field programmable gate arrays (FPGAs).

[0037] According to another innovative aspect of the present disclosure, a computer-implemented method for compressing genomic sequence data generated by a sequencing machine is disclosed, the genomic sequence data including reads of sequences of nucleotides or bases that have been aligned to a reference sequence, thereby generating aligned reads, the aligned reads being stored as a read list in an initial file. In one aspect, the method may include the following actions for each aligned read: determining whether the read is fully mapped or incompletely mapped to the reference sequence, or whether the read is unmapped to the reference sequence; encoding the reads based on the determination, wherein the reads determined to be fully mapped are encoded according to a first encoding process, and the reads determined to be unmapped are encoded according to a second encoding process, wherein the determining step includes, for each incompletely mapped read, comparing the number of mismatches between the read and the reference sequence with a threshold value, wherein, in the encoding step, the reads are encoded according to the second encoding process or the third encoding process. The read segments determined to be incompletely mapped are encoded, when the number of mismatches is greater than the threshold, the incompletely mapped read segments are encoded according to the second encoding process, and when the number of mismatches is less than the threshold, the incompletely mapped read segments are encoded according to the third encoding process, wherein in the second encoding process, each nucleotide or base of the read segment is individually encoded, wherein the first encoding process and the third encoding process include different descriptor sets, each descriptor set univocally represents the read segment associated with the corresponding encoding process, and each of the first encoding process and the third encoding process is a simplified information source entropy encoding process.

[0038] Other aspects include corresponding systems, apparatus, and computer programs that perform the actions of the methods as disclosed herein, as defined by instructions encoded on a computer-readable storage device.

[0039] These and other versions may optionally include one or more of the following features. For example, in some implementations, the determining step may include a further determination, upon determining that the read is not completely mapped to the reference sequence and has a number of mismatches below the threshold, of whether the read is globally mapped or locally mapped to the reference sequence, and wherein the third encoding process includes a first encoding sub-process, encoding the reads determined to be globally mapped according to the first encoding sub-process, and encoding the reads determined to be locally mapped according to the second encoding sub-process, and a second encoding sub-process, the first encoding sub-process and the second encoding sub-process including different sets of descriptors, each descriptor set univocally representing the read associated with the corresponding encoding sub-process.

[0040] In some specific implementations, the descriptor of the first encoding sub-process may include the alignment start position in the reference sequence, the read length, and a mismatch list represented by symbolic replacements, and wherein the descriptor of the second encoding sub-process includes the local alignment start position in the reference sequence, the read length, a mismatch list represented by symbolic replacements, and the length of the cropped portion of the read that is not part of the alignment.

[0041] In some implementations, in the encoding step, the clipped portions of the reads to be encoded according to the second encoding sub-process are concatenated, and each nucleotide or base of the clipped portions is encoded separately.

[0042] In some implementations, in the encoding step, each mismatch of an incompletely mapped read segment is encoded on 1 byte.

[0043] In some embodiments, in the encoding step, each mismatch of the incompletely mapped read is encoded as follows: the first two bits of the byte are used to encode the alternative nucleotide or base present in the read instead of the corresponding reference nucleotide or base in the reference sequence, and the last six bits of the byte are used to encode the position of the mismatch in the reference sequence, which is calculated as an offset relative to the previous mismatch of the read.

[0044] In some implementations, in the encoding step, if the offset calculated between a given mismatch and the previous mismatch is greater than a maximum encodable value, at least one false mismatch is inserted between the two mismatches until each offset between each of the mismatches and the at least one false mismatch is below the maximum encodable value, a false mismatch being defined as a mismatch for which bits of the byte are used to encode the mismatch or to encode a nucleotide or base that is equal to a corresponding reference nucleotide or base in the reference sequence.

[0045] In some implementations, an initial step is to divide the read list into read blocks, where each block begins with a header containing information needed to decode the block, and where the compression method proceeds block by block.

[0046] In some implementations, the read chunks have the same chunk size.

[0047] In some embodiments, the final step is to provide a compressed file containing a list of encoded reads stored in the compressed file in the same order as the reads were stored in the initial file.

[0048] In some implementations, the threshold is equal to 31.

[0049] In some implementations, for each aligned read, a step is provided of determining whether the read contains at least one mismatch corresponding to a case where the sequencing machine was unable to call any base or nucleotide.

[0050] In some implementations, for reads that each contain at least one mismatch corresponding to an instance where the sequencing machine was unable to call any base or nucleotide, there is provided a step of determining the number of such mismatches and a step of comparing the number to a reference threshold.

[0051] In some specific implementations, in the encoding step, if the number of such mismatches is greater than the reference threshold, each nucleotide or base of the read to be encoded according to the second encoding process is individually encoded with 4 bits, and if the number of such mismatches is less than the reference threshold, each nucleotide or base of the read to be encoded according to the second encoding process is individually encoded with 2 bits, and the encoding step also includes encoding a list of positions along the reference sequence that correspond to the positions of such mismatches in the reference sequence. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 is a flow chart showing the steps of a compression method according to the present invention.

[0053] Figure 2 is a schematic diagram illustrating an apparatus for implementing the steps of the compression method according to the present invention.

[0054] Figure 3 A first example of reads globally mapped to a reference sequence is shown.

[0055] Figure 4 A second example of reads globally mapped to a reference sequence is shown where false mismatches must be inserted. DETAILED DESCRIPTION

[0056] The genomic sequences mentioned in the present invention include, for example, but not limited to, nucleotide sequences, deoxyribonucleic acid (DNA) sequences, ribonucleic acid (RNA) sequences, and amino acid sequences. Although the specific embodiments herein describe genomic information in the form of nucleotide sequences in considerable detail, it should be understood that, as will be appreciated by those skilled in the art, the compression methods according to the present invention can also be implemented for other genomic sequences, albeit with some variations.

[0057] Genomic sequencing information is generated by a sequencing machine in the form of a sequence of nucleotides (or more generally, bases) represented by a string of letters from a limited vocabulary. The minimal vocabulary is represented by five symbols: {A, C, G, T, N}, representing the four types of nucleotides present in DNA: adenine, cytosine, guanine, and thymine. In RNA, thymine is replaced by uracil (U). An N indicates that the sequencing machine cannot detect any base, so the true nature of the position is uncertain.

[0058] The nucleotide sequences produced by a sequencing machine are called "reads". The length of a sequence read can range from tens of nucleotides to thousands of nucleotides. Some technologies produce sequence reads in pairs, where one read is from one DNA strand and the second read is from the other strand. Throughout this disclosure, a "reference sequence" is any sequence on which a nucleotide or base sequence produced by a sequencing machine is aligned / mapped. An example of such a reference sequence can actually be a reference genome, a sequence assembled by scientists as a representative example of a collection of genes for a species. However, a reference sequence can also consist of synthetic sequences that are conceived to simply improve the compressibility of the reads, given that they are to be further processed. Sequencing machines may introduce errors in sequence reads, and it is worth noting that the wrong symbols (i.e., representing different nucleic acids) may be used to represent nucleic acids or bases that are actually present in the sample being sequenced; this is commonly referred to as substitution errors or "mismatches".

[0059] The present invention is a reference-based compression method that receives as input reads of a nucleotide or base sequence, such reads having been previously aligned with a reference sequence, thereby generating aligned reads. The aligned reads are then stored as a list of reads in an initial file. The aligned reads and the manner in which these reads are stored once aligned in the initial file are not critical to the present invention and are not the purpose of this disclosure. Each read is then encoded as a list of its position on a reference sequence and its differences from the reference sequence. Each read can then be reconstructed from the alignment encoding information and the reference sequence by appropriate decompression software configured in accordance with the present invention.

[0060] Preferably, the alignment software that processes the reads and aligns them with the reference sequence before providing them as input to the compression software and device does not take into account certain types of errors introduced in the sequence reads, such as insertion errors or deletion errors. An insertion error consists in the insertion of one or more extra symbols in a sequence read that do not relate to any actually existing nucleic acid. A deletion error consists in the deletion of one or more symbols representing nucleic acids actually present in the sequencing sample from a sequence read. More precisely, in the case of an insertion error or deletion error in a given sequence read, the alignment software then treats the resulting erroneous nucleic acid as a substitution error, also referred to as a "mismatch". This preference for alignment software configuration allows for faster subsequent encoding, in particular providing a better compromise between speed and compression ratio.

[0061] For each read, the alignment software provides a corresponding read record to the compression software and device. Each read record contains at least the following information: the absolute starting position of the aligned read relative to the reference sequence, the length of the read, the alignment type of the read, the number of mismatches identified in the read, and the relative position of the potential mismatch in the read (where appropriate).

[0062] Now refer to Figure 1 The compression method according to the present invention is described. Figure 2 20. The apparatus includes at least one processor 22 and a memory 24 operatively coupled to the processor 22 to form a computing device. The memory 24 may store computer program code or software 26 comprising computer executable instructions that, when executed by the processor 22, cause the processor 22 to perform operations including the steps of the compression method according to the present invention.

[0063] The initial file in which the aligned reads are stored as a read list is stored, for example, in the memory of the device 20. Figure 1, the method preferably includes an initial step 2 in which the initial list of aligned reads is divided into read blocks. Typically, the list of aligned reads is divided into 50,000 read blocks, this particular value not being construed as limiting the scope of the invention which may be applied in the same manner with other values. Preferably, the read blocks have the same block size. Each read block begins with a header containing information required to decode the block, such as the size in bytes of the contents of the block, and / or an identifier of the block or its contents, and / or the number of reads contained in the block. This allows for support of concatenation of compressed files, as well as streaming capabilities (each read block contains all the information required to decode the reads of the block). Furthermore, since the compression method can then proceed block by block, this also allows for multi-threaded processing of the read blocks, thereby allowing for parallelization and some resulting gains in processing time. If all reads for a given block have the same length, the read lengths are also stored in the header, otherwise a list of each read length is explicitly stored during the compression method.

[0064] Each read record contains information about the type of read alignment. In general, two main types of alignments can be identified: complete alignment and incomplete alignment, plus an additional type corresponding to "unmapped" reads. "Incomplete alignment" means that the read contains at least one non-N mismatch, while at least a portion of the read matches a portion of the reference sequence (according to this definition, an incompletely mapped read can contain one or more Ns if it also contains one or more other mismatches). In an exemplary embodiment, each read record begins with the following bit markers, each of which has a value between two possible values:

[0065] - a first marker indicating forward or reverse orientation relative to a reference sequence,

[0066] - a second flag indicating a complete or incomplete alignment,

[0067] - a third bit flag indicating whether the read segment contains at least one N,

[0068] - The fourth bit flag indicates whether the position information is encoded in 16 bits or 32 bits.

[0069] The following steps 4-12 are performed on a read-by-block basis, and within a block on a read-by-read basis.

[0070] The method comprises a next step 4 of determining, for each aligned read, whether the read is completely mapped or incompletely mapped to the reference sequence, or whether the read is unmapped to the reference sequence. This determination step 4 comprises, for each incompletely mapped read, comparing the number of mismatches between the read and the reference sequence to a threshold value 4A. In a preferred embodiment, although this should not be construed as limiting the scope of the invention, the threshold value is equal to 31. This particular value is deliberately chosen so as to provide the best possible compromise for storing the number of mismatches in a sufficiently compact manner, as will be better understood later with respect to step 12. In fact, it has been statistically observed that in the vast majority of cases, incompletely mapped reads have fewer than 31 mismatches. The principle behind this selection is to encode the most frequently occurring cases in the most compact manner, leaving behind some very few degraded cases. If the read is determined to be incompletely mapped with a number of mismatches below the threshold, the determination step 4 comprises a further determination as to whether the read is globally mapped or locally mapped to the reference sequence. A "globally mapped read" is an incompletely mapped read whose entire sequence (including the start and end points of the read) is incompletely mapped to the reference sequence. A "locally mapped read" is an incompletely mapped read that contains a segment of nucleotides or bases that is incompletely mapped to the reference sequence. Thus, the segment of nucleotides or bases corresponds to a portion of the original read.

[0071] Preferably, the method further comprises a step 6 of determining, for each aligned read, whether the read contains at least one N, i.e., whether the read contains at least one mismatch corresponding to a case where the sequencing machine could not call any base or nucleotide. For each read containing at least one N, the method then comprises a step 8 of determining the number of such N mismatches; and a step 10 of comparing the number of N mismatches to a reference threshold. In a preferred embodiment, although this should not be construed as limiting the scope of the invention, the reference threshold is equal to 31.

[0072] Regardless of the result of determining step 4, the method includes a next step 12, which encodes the reads at least according to the determination. More precisely, the reads determined to be completely mapped to the reference sequence are encoded according to the first encoding process, regardless of whether the reads do not contain N or have a number of N below the reference threshold. The reads determined to be unmapped or the reads determined to be completely mapped but have a number of N greater than the reference threshold are encoded according to the second encoding process, wherein each nucleotide or base is encoded separately, regardless of whether the nucleotide or base is aligned or not. The reads determined to be incompletely mapped are encoded according to the second encoding process or the third encoding process. More precisely, the reads determined to be incompletely mapped with a number of mismatches greater than the threshold are encoded according to the second encoding process. If the read is determined to be incompletely mapped with a number of mismatches below the threshold, the read is encoded according to the third encoding process if the read does not contain N or has a number of N below the reference threshold. If this is not the case, ie if the reads have a number N greater than a reference threshold, the reads are encoded according to a second encoding process.

[0073] Regardless of whether a given read has been determined to be fully mapped, incompletely mapped, or unmapped, if the read contains at least one N but has a number of N below a reference threshold, the encoding step 12 includes encoding a list of positions along the reference sequence that correspond to the positions of N in the reference sequence. This list of positions is then stored in a memory of a computing device that implements the compression method. If the read contains at least one N but has a number of N below the reference threshold, it is encoded according to a second encoding process, and each nucleotide or base in the read is individually encoded using 2 bits.

[0074] If a read contains at least one N but has a number of N greater than a reference threshold, the read is encoded in any case according to the second encoding process, and each nucleotide or base of the read is encoded individually with 4 bits. In this case, the encoding step 12 does not include encoding and storing a list of N positions in the reference sequence. In fact, each N mismatch is then directly encoded according to the second encoding process in much the same way as the other nucleotides or bases of the read.

[0075] The first encoding process and the third encoding process include different descriptor sets. Each descriptor set univocally represents a read segment associated with the corresponding encoding process, and each of the first encoding process and the third encoding process is a simplified information entropy encoding process. More precisely, the third encoding process includes a first encoding sub-process and a second encoding sub-process. The read segments of the incomplete mapping determined to be globally mapped during step 4 are encoded according to the first encoding sub-process. The read segments of the incomplete mapping determined to be locally mapped during step 4 are encoded according to the second encoding sub-process. The first encoding sub-process and the second encoding sub-process include different descriptor sets, each descriptor set univocally represents a read segment associated with the corresponding encoding sub-process.

[0076] The alignment information encoded for each read and enabling reconstruction of the entire read sequence during data decompression then depends on the corresponding encoding process or sub-process used for the read. For example, the descriptor for the first encoding process could be:

[0077] ○ The absolute starting position of the fully mapped read relative to the reference sequence (encoded in 16 or 32 bits), and

[0078] o The length of the read (encoded relative to the length of the previous read using differential encoding, with variable length codes ranging from 2 to 34 bits).

[0079] The descriptor for the first encoding sub-process may be:

[0080] ○ The absolute starting position of the incompletely mapped read relative to the reference sequence (encoded in 16 or 32 bits),

[0081] o The length of the read (encoded relative to the length of the previous read using differential encoding, with variable length codes ranging from 2 to 34 bits), and

[0082] ○ List of mismatches for the reads.

[0083] The descriptor for the second encoding sub-process may be:

[0084] ○ The absolute starting position of the incompletely mapped portion of the read relative to the reference sequence – also called the local alignment start position (encoded in 16 or 32 bits),

[0085] o The length of the read (encoded relative to the length of the previous read using differential encoding, with variable length codes ranging from 2 to 34 bits),

[0086] ○ A list of mismatches for the reads, and

[0087] o The length of the clipped parts of the read that were not part of the alignment (encoded in 8 bits for each clipped part).

[0088] Preferably, the mismatch list encoded in the first sub-process and the second sub-process includes a header (a bit flag encoded on 1 byte). The first five bits of this byte are used to encode the number of mismatches contained in the read segment (in a preferred embodiment, where the threshold is equal to 31, the number is in the range [0-31]). Then one bit is used to encode whether the incompletely mapped read segment is a global map or a local map. Another bit is used to encode whether the 2-bit mode is activated for the second encoding process. The last bit is used to encode whether the 4-bit mode is activated for the second encoding process. Preferably, for each read segment encoded according to the second encoding sub-process during encoding step 12, the cut portion of the read segment (i.e., those parts that are not part of the local alignment) is concatenated, and each nucleotide or base of the cut portion is encoded separately. In a preferred embodiment, each nucleotide or base of such a cut portion of the read segment is encoded separately with 2 bits.

[0089] In a preferred embodiment, each mismatch encoded in the mismatch list of an incompletely mapped read (i.e., encoded according to the first encoding sub-pass or the second encoding sub-pass) is encoded on 1 byte. More precisely, each mismatch of an incompletely mapped read to be encoded according to the first encoding sub-pass or the second encoding sub-pass may be encoded as follows:

[0090] o The first two bits of the byte are used to encode the alternative nucleotide or base present in the read instead of the corresponding reference nucleotide or base in the reference sequence,

[0091] o The last six bits are used to encode the position of the mismatch in the reference sequence, which is calculated as an offset relative to the previous mismatch of the read (the relative position of the mismatch, except for the absolute position of the first mismatch of the read, which is encoded). Therefore, the range of this offset (encoded in 6 bits) is [0-63].

[0092] Figure 3 An example of encoding mismatches of a read according to a first encoding sub-process is provided. The read is an incompletely mapped read that is globally mapped to a reference sequence. The read has two mismatches:

[0093] o The first mismatch, located at position 12 in the read, consists in replacing an A nucleotide in the reference sequence with a T nucleotide in the read, and

[0094] o The second mismatch, located at position 21 in the read, consists in the replacement of a C nucleotide in the reference sequence with a G nucleotide in the read.

[0095] The mismatch list for this read is then encoded as:

[0096] ○ <12,T>, the value "12" corresponds to the absolute position of the first mismatch in the read segment, and

[0097] ○<9,G>, the value “9” corresponds to the relative position of the second mismatch in the read segment, ie, the offset between the second mismatch and the first mismatch.

[0098] For example, <12, T> can be converted to the value "51" (encoded on 1 byte), and <9, G> can be converted to the value "38" (encoded on 1 byte). This byte encoding is obtained in the following way:

[0099] Offset position x 4 + nucleotide value (where A=0, C=1, G=2, T=3)

[0100] Preferably, for each incompletely mapped read to be encoded according to the first encoding sub-process or the second encoding sub-process, if the offset calculated between a given mismatch in the read and the previous mismatch is greater than the maximum encodable value, at least one "false" mismatch is inserted between the two mismatches until each offset between each of the mismatches and the at least one "false" mismatch is below the maximum encodable value. A "false" mismatch is defined as a mismatch for which the bits of a byte are used to encode the mismatch or to encode a nucleotide or base that is equal to the corresponding reference nucleotide or base in the reference sequence. In a preferred embodiment, although this should not be construed as limiting the scope of the invention, the maximum encodable value is equal to 63, corresponding to the maximum value encoded in 6 bits.

[0101] Figure 4 An example of encoding mismatches of a read according to the first encoding sub-process is provided in the case where a "false" mismatch must be inserted. The read is an incompletely mapped read that is globally mapped to the reference sequence. The read has two mismatches:

[0102] o The first mismatch, located at position 22 in the read, consists in replacing an A nucleotide in the reference sequence with a T nucleotide in the read, and

[0103] o The second mismatch, located at position 134 in the read, consists in the replacement of a C nucleotide in the reference sequence with a G nucleotide in the read.

[0104] The position offset between the second mismatch and the first mismatch is 112, which is greater than the maximum codable value of 63. Therefore, a "false" mismatch must be inserted between the two mismatches so that each offset between each of the mismatches and the "false" mismatch is lower than the maximum codable value. A "false" mismatch with a T nucleotide (corresponding to a "true" T nucleotide in the reference sequence) is inserted, for example, at position 85 in the read. The calculated position offset between this "false" mismatch and the first mismatch is 63, which corresponds to the maximum codable value. The calculated position offset between the second mismatch and the "false" mismatch is 49, which is lower than 63.

[0105] The mismatch list for this read is then encoded as:

[0106] ○ <22,T>, the value "22" corresponds to the absolute position of the first mismatch in the read segment,

[0107] ○<63,T>, the value "63" corresponds to the relative position of the "false" mismatch in this read segment, i.e.

[0108] the offset between the “false” mismatch and the first mismatch, and

[0109] ○ <49,G>, the value “49” corresponds to the relative position of the second mismatch in the read segment, ie, the offset between the second mismatch and the “false” mismatch.

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

[0111] Offset position x 4 + nucleotide value (where A=0, C=1, G=2, T=3)

[0112] The method includes a final step 14 of providing 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 original uncompressed file. Each read can then be reconstructed from the alignment of the encoded information and the reference sequence using appropriate decompression software and / or methods configured according to the present invention.

[0113] Although reference is made to the exemplary architecture of computing device 20 (for illustrative purposes in Figure 2 ), but the technology of the present invention disclosed herein can be implemented in hardware, software, firmware, or any combination thereof. When implemented in software, the computer program code can be stored on a computer medium and executed by a hardware processing unit including one or more processors, which is similar to using Figure 2The same is true for the device 20. It should be understood that the term "processor" as used herein is intended to include one or more processing devices, including signal processors, microprocessors, microcontrollers, application specific integrated circuits (ASICs), field programmable gate arrays (FPGAs), or other types of processing circuit systems, as well as parts or combinations of such circuit system elements. In addition, the term "memory" as used herein is intended to include electronic memory associated with a processor, such as random access memory (RAM), read-only memory (ROM), or other types of memory used in any combination.

[0114] Thus, software instructions or code for executing the methods and protocols described herein may be stored in one or more of the associated memory devices (e.g., ROM, fixed memory, or removable memory) and, when ready for use, loaded into RAM and executed by the processor.

[0115] The techniques of this disclosure may be implemented in a wide variety of devices or apparatuses, including, for example, mobile phones, computers, servers, tablet computers, and similar devices.

[0116] Although illustrative embodiments of the present invention have been described herein with reference to the accompanying drawings, it will be understood that the invention is not limited to those precise embodiments and that various other changes and modifications may be made by those skilled in the art without departing from the scope or spirit of the invention.

[0117] Statistical and numerical examples of the compression method according to the present invention

[0118] The following comparative example was performed on an uncompressed data file containing 48 million nucleotide reads or sequences:

[0119] ○ Uncompressed data file size: 35,770MB (megabytes)

[0120] ○ File size after gzip compression: 6,649MB

[0121] ○ File size compressed with non-reference-based SPRING software: 1,402MB

[0122] ○ Size of the file compressed using the reference-based compression method according to the present invention: 1,179

[0123] MB

[0124] Compression time using non-reference-based SPRING software: 1,722 seconds

[0125] ○ Compression time using the reference-based compression method according to the present invention: 181s

[0126] Average size of uncompressed data files in bits / nucleotide (ASCII encoding): 8 bits

[0127] / nucleotide

[0128] ○ Average size in bits / nucleotide of files compressed with an encoding suitable for the 4 possible characters A, T, C, G: 2 bits / nucleotide

[0129] o Average size in bits / nucleotide of files that have been compressed using the reference-based compression method according to the present invention: 0.33 bits / nucleotide

[0130] The numerical examples indicated above illustrate that the present invention allows fast compression and decompression while providing a high compression ratio.

Claims

1. A method for compressing genomic sequence data, the method comprising: Obtaining read records by one or more computers; determining, by the one or more computers, whether the read record corresponds to a read that is completely mapped to a reference sequence or a read that is incompletely mapped to the reference sequence; Based on a determination by the one or more computers that the read record corresponds to a read that is incompletely mapped to the reference sequence, determining by the one or more computers whether a number of mismatches for the incompletely mapped reads satisfies a predetermined threshold number of mismatches; and Based on determining that the number of mismatches satisfies the predetermined threshold number of mismatches: For a specific mismatch of the stated mismatch threshold number: obtaining, by the one or more computers, an offset relative to a previous mismatch, determining, by one or more computers, whether the offset is greater than a maximum encodable offset value, wherein the maximum encodable offset value is defined by a first number of bits of a predetermined record size minus a second number of bits encodable by a predetermined encoding of the particular mismatch, the first number of bits being greater than the second number of bits; Based on determining that the offset is greater than the maximum encodable offset value, at least one false mismatch is inserted by one or more computers, wherein the false mismatch is (i) inserted between the particular mismatch and the previous mismatch, and (ii) serving as a reference point between the particular mismatch and the previous mismatch; and determining, by one or more computers, (i) a first offset between the false mismatch and the previous mismatch, and (ii) a second offset between the false mismatch and the particular mismatch; and The one or more computers (i) encode the position of the previous mismatch and the previous mismatch into a first record having the predetermined record size, (ii) encode the first offset and the false mismatch into a second record having the predetermined record size, and (iii) encode the second offset and the specific mismatch into a third record having the predetermined record size. The method according to claim 1 , wherein the predetermined record size is 1 byte.

3. The method of claim 1 or 2, wherein determining, by the one or more computers, whether the number of mismatches of the incompletely mapped read segments satisfies a predetermined threshold number of mismatches comprises: A determination is made by the one or more computers whether the number of mismatches for the incompletely mapped reads is greater than the predetermined threshold number of mismatches.

4. The method according to claim 1 or 2, wherein each read record comprises: data indicating the absolute starting position of the aligned read relative to the reference sequence; data indicating the length of the read segment; data indicating whether the read segment is completely mapped or incompletely mapped; data indicating a number of mismatches identified in the read segment; as well as Data indicating a relative position of each of the mismatches in the read segment.

5. The method according to claim 1 or 2, further comprising: Based on determining that the number of mismatches does not satisfy the predetermined threshold number of mismatches, encoding, by one or more computers, a list of positions of the reference sequence corresponding to the position of each of the mismatches into the reference sequence using a simplified information entropy encoding process.

6. The method according to claim 1 or 2, further comprising: Based on determining that the read record corresponds to a read that is completely mapped to the reference sequence, encoding at least a portion of the read record using simplified information entropy coding by one or more computers.

7. The method of claim 1 or 2, wherein the one or more computers comprise one or more hardware processors.

8. The method of claim 7, wherein the one or more hardware processors comprise one or more field programmable gate arrays (FPGAs).

9. The method of claim 1 , wherein encoding the location of the previous mismatch and the previous mismatch into a first record having the predetermined record size comprises: encoding, by the one or more computers, two bits of the predetermined record size to include data representing an alternative nucleotide or base present in the read associated with the previous mismatch; as well as The other six bits of the predetermined record size are encoded by one or more computers to include data representing the location of the previous mismatch.

10. The method of claim 1 , wherein encoding the first offset and the false mismatch into a second record having the predetermined record size comprises: encoding, by the one or more computers, two bits of the predetermined record size to include data representing a nucleotide or base present in a read associated with the false mismatch; as well as The other six bits of the predetermined record size are encoded by one or more computers to include data representing the first offset.

11. The method of claim 1 , wherein encoding the second offset and the specific mismatch into a third record having the predetermined record size comprises: encoding, by the one or more computers, two bits of the predetermined record size to include data representing an alternative nucleotide or base present in the read associated with the particular mismatch; as well as The other six bits of the predetermined record size are encoded by one or more computers to include data representing the second offset.

12. A system for compressing genomic sequence data, the system comprising: One or more computers, and one or more storage devices storing instructions, the instructions being operable when executed by the one or more computers to cause the one or more computers to perform operations including the method of any one of the preceding claims.

13. A computer-readable storage device having instructions stored thereon that, when executed by a data processing apparatus, cause the data processing apparatus to perform operations including the method of any one of claims 1 to 11.

14. A hardware processor comprising hardware processing circuitry configured to perform operations including the method of any one of claims 1 to 11.

Images