EA040022B1 - METHOD AND DEVICE FOR COMPACT REPRESENTATION OF BIOINFORMATICS DATA - Google Patents

Description

The invention provides a new way to represent genome sequencing data that reduces memory usage and improves access performance by providing new functionality that is not available with prior art presentation methods.

State of the art

Proper representation of genome sequencing data is fundamental to enable efficient applications for genomic analysis, such as the identification of genome variants and any other analysis performed for various purposes by processing sequencing data and metadata.

Human genome sequencing has become accessible due to the advent of high-throughput, low-cost sequencing technologies. This capability opens up new vistas in several areas, from cancer diagnosis and treatment to detection of genetic diseases, from pathogen surveillance to detect antibodies, to the development of new vaccines, drugs and personalization of individualized treatment.

Hospitals, genome analysis service providers, bioinformaticians, and large biological data warehouses are looking for low-cost, fast, reliable, and interconnected genomic information processing solutions to scale up genomic medicine to the global scale. Since data storage has become one of the bottlenecks in the sequencing process, ways of presenting genome sequencing data in a compressed form are being increasingly explored.

The most commonly used representations of genome information based on sequencing data are based on the FASTQ and SAM formats. The challenge is to compress traditionally used file formats (respectively FASTQ and SAM for unaligned and aligned data). Such files consist of plain text characters and are compressed, as mentioned above, using general purpose approaches such as LZ schemes (named after Lempel and Ziv, the authors who published the first versions) (well known zip, gzip, etc.). With general purpose compressors such as gzip, the result of compression is usually a single blob of binary data. Information in such a monolithic form is very difficult to archive, transfer and process, especially when, as in the case of high-throughput sequencing, the amount of data is extremely large. The BAM format suffers from poor compression performance due to its focus on compressing the inefficient and redundant SAM format rather than extracting the actual genomic information carried by SAM files, and due to the use of general purpose text compression algorithms such as gzip instead of in order to exploit the specific nature of each data source (the genomic data itself).

CRAM is a more sophisticated approach to genomic data compression that is less used but more efficient than BAM. CRAM provides more efficient compression to apply differential encoding to the existing reference (it partially exploits data source redundancy), but it still lacks features such as incremental updates, streaming support, and selective access to certain classes of compressed data. These approaches produce unsatisfactory compression ratios and data structures that are difficult to navigate and manipulate after compression. Subsequent analysis can be very slow due to the need to process large and rigid data structures even to perform a simple operation or to access selected regions of a genomic data set. CRAM relies on the concept of CRAM recording. Each CRAM entry encodes one mapped or unmapped read, encoding all the elements needed to reconstruct it.

CRAM has the following disadvantages.

1) For CRAM, data indexing is outside the scope of the specification (see section 12 of the CRAM v 3.0 specification) and is implemented as a separate file. In contrast, the approach of the invention described in this document employs a data indexing technique that is integrated with the encoding process and the indices are embedded in the encoded bitstream.

2) In CRAM, all basic data blocks can contain any type of mapped reads (perfect match reads, substitution-only reads, insertion or deletion reads (also called indels)). There is no concept of classification and grouping of reads in classes according to the result of mapping relative to the reference sequence.

3) In the invention, there is no concept of a record encapsulating each read, because the data needed to restore each read is distributed over many data containers called layers. This provides more efficient access to a set of reads with specific biological characteristics (eg, reads with substitutions but no indels, or perfectly mapped reads) without having to decode each (block) of the read(s) to check for its features.

4) In a CRAM record, each data type is identified by a specific flag. In the present invention, there is no concept of a flag denoting data, because this is essentially determined by the layer to which the data belongs. This implies a significantly reduced number of symbols used.

- 1 040022 oxen and, as a result, a decrease in the entropy of the information source, which leads to more efficient compression. This is due to the fact that using different layers allows the encoder to reuse the same symbol in each layer with different values. In CRAM, each flag must always have the same value, because there is no concept of contexts, and each CRAM entry can contain data of any type.

5) In CRAM substitutions, insertions and deletions are expressed according to different syntaxes, while the proposed approach uses the same alphabet and coding for substitutions, insertions and deletions. This simplifies the process of encoding and decoding and creates a lower entropy source model whose encoding yields highly compressed bitstreams.

The invention is directed to compressing genomic sequences by organizing and partitioning data in such a way as to minimize redundant information to be encoded and to implement features such as selective access and support for incremental updates. One aspect of the presented approach is the definition of data and metadata classes that should be encoded separately and that should be structured in different layers. The most important improvements of this approach over existing methods are as follows:

1) an increase in compression efficiency due to a decrease in the entropy of the information source created by providing an efficient model for each class of data or metadata;

2) the ability to perform selective access to portions of compressed data and metadata for further processing;

3) the ability to incrementally (without the need for recoding) update the encoded data and metadata with new data and/or sequence metadata and/or new analysis results.

Brief description of the drawings

In FIG. 1 shows that the position of pairs of mapped reads is encoded in the pos layer as the difference with the absolute position of the first mapped read;

fig. 2 - two reads in a pair can come from two strands of DNA;

fig. 3 - the reverse complement of read 2 will be encoded in this way if strand 1 is used as a reference;

fig. 4 - four possible combinations of reads that make up a pair of reads and the corresponding encoding in the rcomp layer;

fig. 5 - scheme for calculating the mating distance at a constant length of reads for three pairs of reads;

fig. 6 - pairing errors encoded in the pair layer allow the decoder to recover the correct pairing of reads using the encoded MPPPD;

fig. 7 - coding of mating distance when a reed is mapped on a different reference than its partner in a pair. In this case, additional descriptors are added to the pairing distance. One is the signal flag, the second is the reference ID, and then the mating distance;

fig. 8 - coding of N-mismatches in the nmis layer;

fig. 9 is a mapped pair of reads that represents substitutions relative to the reference sequence;

fig. 10 - calculation of the positions of replacements both in absolute and differential values;

fig. 11 - calculation of characters encoding substitution types when IUPAC codes are not used. Symbols represent the distance - in a circular substitution vector - between a molecule present in the read and a molecule present in the reference at that position;

fig. 12 - coding of substitutions in the snpt layer;

fig. 13 - calculation of substitution codes using IUPAC ambiguity codes;

fig. 14 - snpt layer coding using IUPAC codes;

fig. 15 - for reads of class I, the substitution vector used is the same as for class M with the addition of special codes for insertions of symbols A, C, G, T, N;

fig. 16 are some examples of mismatch and indel coding when using IUPAC ambiguity codes. The substitution vector in this case is much longer, and therefore there are more possible calculated characters than in the case of five characters;

fig. 17 is another source model for mismatches and indels, where each layer contains the position of mismatches or insertions of the same type. In this case, the characters for the mismatch or indel type are not encoded;

fig. 18 is an example of mismatch and indel coding. If there are no mismatches or indels of that type for a read, the value 0 is encoded in the corresponding layer. The character 0 acts as the delimiter and terminator of the reads in each layer;

fig. 19 - Modification in the reference sequence can convert M-read to R-read. This operation can reduce the informational entropy of the data structure, especially in the case of high overlap;

fig. 20 is a genomic encoder 2010 according to one embodiment of this invention;

fig. 21 is a genomic decoder 218 according to one embodiment of this invention.

- 2 040022

Short description

The invention relates in particular to the following objects.

1) A computer-implemented method for compressing genomic sequence data, wherein said genomic sequence data contains reads of nucleotide sequences, comprising the following steps:

aligning said reads with one or more reference sequences, thereby obtaining aligned reads, classifying said aligned reads into different classes, including at least the first class: if said aligned reads match said one or more reference sequences without any errors ;

second class: if the specified aligned reads match an area in the specified one or more reference sequences with a mismatch count consisting of the number of positions where the sequencer could not determine any base;

third class: if said aligned reads match a region in said one or more reference sequences with a number of mismatches consisting of positions where the sequencer could not detect any bases and the presence of insertions or deletions or truncated nucleotides;

fourth class: if said aligned reads do not find any valid mapping on said one or more reference sequences according to said alignment constraints, thus obtaining aligned read classes; and encoding said classified and aligned reads as a plurality of layers of syntactic elements containing descriptors, wherein for said first class, said descriptors include at least an initial position in the reference sequence, a flag signaling that the read should be considered as an inverse complement to the reference, distance to the partner of the pair in the case of paired reads, length value in the case when the sequencing technology produces reads of variable length; for said second class, the descriptors include at least the first class descriptors and a mismatch position for each mismatch; for said third class, the descriptors include the descriptors of said second class and the mismatch position and mismatch type for each mismatch; for the specified fourth class, the descriptors include the descriptors of the specified first class and the mismatch type for each mismatch, where the specified encoding of the specified classified aligned reads as a plurality of layers of syntax elements includes the selection of the specified syntax elements containing the specified descriptors, in accordance with the specified classes of aligned reads, and the encoding of the specified classified aligned reads in the form of a plurality of layers of syntactic elements is carried out using a specific entropy encoder 2012-2014.

2) The method according to claim 1, characterized in that the layers of syntax elements additionally contain the position of the variant relative to the reference sequence, the type of the variant, the position of the deletion relative to the reference sequence, the position of one or more characters that are absent in the reference sequence but are present in aligned reads, the type inserts in this position.

3) The method of claim 1, wherein said entropy encoder is a context-adaptive arithmetic encoder.

4) A method for decompressing a genomic stream compressed by the method of claim 1, said method comprising the following steps:

parsing and decoding 212-214 of the compressed genomic stream into genomic layers of syntactic elements 215, decoding said genomic layers into data classes 216-217, unfolding said genomic layers into classified nucleotide sequence reads, selective decoding using decoders 219 classes of said classified nucleotide sequence reads and combining the results on one or more reference sequences to obtain uncompressed reads of nucleotide sequences.

5) A genomic encoder 2010 for compressing genomic sequence data 209, said genomic sequence data 209 comprising reads of nucleotide sequences, said genomic encoder 2010 comprising an alignment module 201 configured to align said reads with one or more reference sequences, thereby obtaining aligned reads, a data classification module 204 configured to classify said aligned reads into different classes, including at least the first class: if said aligned reads match said one or more reference sequences without any errors;

second class: if the specified aligned reads match the region in the specified one or more reference sequences with the number of mismatches consisting of the number of positions in which the sequencer could not determine any base;

third class: if said aligned reads match a region in said one or more reference sequences with a number of mismatches consisting of positions where the sequencer could not detect any bases and the presence of insertions or deletions or truncated nucleotides;

fourth class: if said aligned reads do not find any valid mapping on said one or more reference sequences according to said alignment constraints, thus obtaining aligned read classes;

one or more layer-encoding modules 205-207 configured to encode said classified aligned reads as syntax element layers containing descriptors by selecting said syntax elements according to said aligned read classes, where for said first class, said descriptors include at least initial position in the reference sequence, a flag indicating that the read should be considered as an inverse complement to the reference, the distance to the partner of the pair in the case of paired reads, the length value in the case when the sequencing technology produces reads of variable length; for said second class, the descriptors include at least the first class descriptors and a mismatch position for each mismatch; for said third class, the descriptors include the descriptors of said second class and the mismatch position and mismatch type for each mismatch; for the specified fourth class, the descriptors include the descriptors of the specified first class and the type of mismatch for each mismatch, entropy encoder 2012-2014 for entropy coding of the specified layers of syntactic elements.

6) A genomic decoder 218 for decompressing a genomic stream 211 compressed by the genomic encoder of claim 6, said genomic decoder 218 comprising parsing and decoding means 210, 212-214 configured to parse said compressed genomic stream into genomic layers of syntax elements 215 , one or more layer decoders 216-217 configured to decode genomic layers into data classes, and further configured to process said genomic layers into classified nucleotide sequence reads 2111, genomic data class decoders 213 configured to selectively decode said classified nucleotide sequence reads, and combining the result for one or more reference sequences to obtain uncompressed reads of nucleotide sequences.

7) Genomic decoder according to claim 6, characterized in that one or more reference sequences are stored in the compressed genome stream 211.

8) The genomic decoder according to claim 6, characterized in that one or more reference sequences are fed into the decoder by an out-of-band mechanism.

9) Genomic decoder according to claim 6, characterized in that one or more reference sequences are built in the specified decoder.

10) A computer-readable medium containing instructions that, when executed, cause at least one processor to perform the method of any one of claims 1-4.

The features of the independent claims below and above solve the problem of the existing prior art solutions by providing a method for classifying genome sequences and a method for compressing using said classification. In one aspect, a method is provided for classifying genomic sequence data generated by a sequencer, said genomic sequence data comprising nucleotide base sequences, said classification being performed according to a reference sequence, said method comprising the following steps:

identifying class P sequences containing matching (matching) regions in the reference sequence without mismatches; identifying class N sequences containing matching regions in the reference sequence, with a number of mismatches represented by positions in which the sequencer could not determine any base; identifying sequences of class M containing matching regions in the reference sequence, with a number of mismatches represented by positions in which the sequencer could not determine any base or determined a base different from the base in the reference sequence; identifying class I sequences containing class M mismatches plus the presence of insertions or deletions;

identification of class U sequences containing all reads that do not find a valid mapping in the reference sequence.

In another aspect, a method is provided for compressing genomic sequence data generated by a sequencer, said genomic sequence data comprising the sequences

- 4 040022 nucleotides, and this method includes the following steps:

aligning said reads with a reference sequence to obtain aligned reads;

classifying said aligned reads according to a multiplicity of degrees of matching with the reference sequence, thereby creating aligned read classes;

encoding said aligned reads as layers of syntactic elements;

and the specified syntax elements are selected in accordance with the specified classes of aligned reads.

In another aspect, a method is provided for decompressing a compressed genomic stream, said method comprising the following steps:

parsing said compressed genomic stream into genomic layers of syntactic elements, unfolding said genomic layers into classified nucleotide sequence reads, selectively decoding said classified nucleotide sequence reads based on one or more reference sequences to obtain uncompressed nucleotide sequence reads.

In yet another aspect, a genomic encoder 2010 for compressing genomic sequence data 209, said genomic sequence data 209 comprising reads of nucleotide sequences, said genomic encoder 2010 comprising an alignment module 201 configured to align said reads with one or more reference sequences to obtain such thus, aligned reads, a data classification module 204 configured to classify said aligned reads according to degrees of accuracy of matching with one or more reference sequences, thereby obtaining aligned read classes;

one or more layer encoding modules 205-207 configured to encode said classified aligned reads as layers of syntax elements by selecting said syntax elements according to said aligned read classes.

In another aspect, a genomic decoder 218 for decompressing a compressed genomic stream 211, said genomic decoder 218 comprising:

parsers 210, 212-214 configured to parse said compressed genomic stream into genomic layers of syntactic elements 215, one or more layer decoders 216-217 configured to decode genomic layers into classified nucleotide sequence reads 2111, genomic class decoders data 213 configured to selectively decode said classified nucleotide sequence reads at one or more reference sequences to produce uncompressed nucleotide sequence reads.

Detailed description

Genomic or proteomic sequences referred to in this invention include, for example, but not limited to, nucleotide sequences, deoxyribonucleic acid (DNA), ribonucleic acid (RNA), and amino acid sequences. Although the description in this document is quite detailed with respect to genomic information in the form of a nucleotide sequence, it should be understood that these storage methods and systems can be implemented for other genomic or proteomic sequences as well, albeit with minor variations, as will be understood by one skilled in the art. areas.

Genome sequencing information is generated by high-throughput sequencers (HTS) in the form of nucleotide (base) sequences, represented by strings of letters from a specific dictionary. The smallest vocabulary is represented by five characters: {A, C, G, T, N} representing 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 was unable to identify any base, and therefore the actual nature of this provision is undetermined. In case the sequencer works with IUPAC disambiguation codes, the alphabet used for characters is (A, C, G, T, U, W, S, M, K, R, Y, B, D, H, V, N or -).

Nucleotide sequences obtained by sequencers are called reads. The read sequence can be from several tens to several thousand nucleotides. Some technologies produce a sequence of reads in pairs, where one read comes from one strand of DNA, and the second from another strand. In genome sequencing, the term overlap is used to express the level of redundancy of sequence data relative to a reference sequence. For example, to achieve 30-fold overlap of the human genome (3.2 billion bases long), the sequencer must produce a total of 30x3.2 billion bases, so that on average each position in the reference is overlapped 30 times.

Throughout this specification, a reference sequence is any sequence against which nucleotide sequences generated by sequencers are aligned/mapped. One example of a sequence would be a reference genome, a sequence compiled by scientists as a representative example of a set of genes in a species. For example, GRC37, the human genome of the Genome Reference Consortium (assembly 37), was obtained from thirteen anonymous volunteers from Buffalo, pc. NY. However, the reference sequence may also consist of a synthetic sequence designed simply to improve the compressibility of the reads in view of their further processing.

Sequencing devices can introduce errors into the read sequence, such as

1) using the wrong symbol (i.e., representing a different nucleic acid) to represent the nucleic acid actually present in the sequenced sample; this is usually called a replacement error (mismatch);

2) insertion in one read of a sequence of additional characters that do not refer to any nucleic acid actually present; this is usually called an insertion error;

3) deletion from one read of a sequence of characters representing nucleic acids that are actually present in the sequenced sample; this is usually called a deletion error;

4) recombination of one or more fragments into one fragment that does not reflect the reality of the original sequence.

The term overlap is used in the literature to quantify the extent to which a reference genome, or part of it, can be overlapped by available sequence reads. An overlap is called partial (less than 1-fold) when some parts of the reference genome are not mapped by any available sequence read;

once (1x) when all nucleotides of the reference genome are mapped with one and only one character present in the sequence read;

multiple (2x, 3x, Nx) when each nucleotide of the reference genome is mapped several times.

The present invention is directed to the development of a genomic information presentation format in which relevant information is efficiently accessible and portable, and the weight of redundant information is reduced.

The main aspects of the disclosed invention are as follows.

1) Classification of sequence reads into different classes based on alignment results relative to reference sequences in order to provide selective access to encoded data according to criteria related to alignment results and match accuracy.

2) Decomposition of a sequence of data reads and metadata into homogeneous layers in order to obtain separate sources of information with reduced information entropy.

3) The ability to model each individual source with a separate source model adapted to each statistical characteristic, including the ability to change the source model in each class of reads and layer for each available data block (access blocks). Adopt appropriate context-adaptive probability models and appropriate entropy encoders according to the statistical properties of each source model.

4) Definition of correspondences and dependencies between layers to provide selective access to data without the need to decode all layers if not all information is needed.

5) Encoding of each sequence data class and associated metadata layers based on a reference sequence, which can be modified to reduce the entropy of data classes and layer information sources. After the first encoding based on the reference sequence, the detected mismatches can be used to adapt/modify the reference sequence in order to further reduce the overall information entropy. This is a process that can be performed iteratively as long as the reduction in information entropy is significant.

Next, each of the above aspects will be described further.

The main header of the file.

Classification of sequence reads.

According to the results of alignment to one or more reference sequences, the read values generated by the sequencers are classified in accordance with the disclosed invention into five different classes.

When aligning a DNA nucleotide sequence with respect to a reference sequence, five outcomes are possible.

1) A region in a reference sequence is found to match a read of the sequence

- 6 040022 without any errors (perfect mapping). Such a sequence of nucleotides will be called a perfect match read or denoted as class P.

2) A region in the reference sequence is found to match a sequence read with multiple (multiple) mismatches consisting of a number of positions where the sequencer could not detect a single base (or nucleotide). Such mismatches are denoted by the letter N. Such sequences will be denoted as mismatched N-reads or class N.

3) A region in the reference sequence is found to match a read of the sequence with several mismatches consisting of a number of positions where the sequencer could not detect any base (or nucleotide), OR another base other than that specified in the reference genome was detected. This type of mismatch is called single nucleotide variation (SNV) or single nucleotide polymorphism (SNP). Such sequences will be referred to as mismatched M-reads or Class M.

4) The fourth class consists of reads representing a type of mismatch, which includes a class M mismatch plus the presence of an insertion or deletion (also called indels). Insertions are represented by a sequence of one or more nucleotides not present in the reference but present in the read sequence. In the literature, when an inserted sequence is at the edges of a sequence, it is said to be soft-truncated (i.e., nucleotides do not match the reference, but are retained in aligned reads as opposed to hard-truncated nucleotides, which are discarded). Deletions are holes (missing nucleotides) in an aligned read with respect to the reference. Such sequences will be called mismatched I-reads or class I.

5) The fifth class includes all reads that do not find any valid mapping on the reference sequence according to the specified alignment constraints. Such sequences are called unmapped and belong to the class U.

The remaining unmapped reads relative to the reference sequence can be assembled into one sequence using de-novo assembly algorithms. After creating a newly assembled reference sequence, non-mapped reads can be further mapped relative to it and classified into one of 4 classes P, N, M and I.

Decomposition of the information needed to represent sequence reads into layers of descriptors.

After the classification of the read is completed with the definition of classes, further processing consists in determining a set of various syntactic elements representing the remaining information, allowing the DNA read sequence to be reconstructed when it is presented as mapped to this reference sequence. The data structure of these syntax elements requires the storage of global parameters and metadata to be used by the decoding engine. This data is structured in the main header described in the table below.

- 7 040022

Table 1. Structure of the main heading

Element Type Description Unique identificator Byte array Unique identifier for encoded content Version Byte array Main + auxiliary encoding algorithm version Header Size Integer Size in bytes of all encoded content Read length Integer The size of the read at a constant length of the read. A special value (e.g. 0) is reserved for variable read length Number of reference sequences Integer Number of reference sequences used Access block counters Byte array (like integers) Total number of encoded access blocks per reference sequence Reference IDs sequences Byte array Unique identifiers for reference sequences Main index table Position Alignment the first read in each block (access block). That is, the smaller position of the first read of the reference genome for each block of 4 classes Byte array (like integers) This is a multi-dimensional array that supports random access to access blocks 1 per pos class (4) per reference

The DNA segment related to a given reference sequence can be fully expressed by the following parameters:

starting position in the reference sequence (pos);

a flag indicating that the read should be treated as a backcomplement to the reference (rsomp);

distance to the partner of the pair in case of paired reads (pair);

the value of the read length in the case when the sequencing technology gives a variable length (len). In the case of a constant read length, the read length associated with each read can obviously be omitted and can be stored in the main header of the file;

for each mismatch, the position of the mismatch (nmis for class N, snpp for class M, and indp for class I), type of mismatch (none in class N, snpt in class M, indt in class I);

flags (296) indicating specific characteristics of the sequence read, such as a template having multiple segments in sequencing;

each segment is correctly aligned according to the aligner;

unmapped segment;

the next segment in the template is not mapped;

-8040022 first or last segment signaling;

failure of quality control;

PCR or optical duplicate;

secondary alignment;

additional alignment;

an optional string of softly truncated nucleotides (indc in class I).

This classification creates groups of descriptors (syntax elements) that can be used to uniquely represent genomic sequence reads. The table below lists the syntax elements required for each aligned read class.

Table 2. Definition of layers for each data class

P Ν ΜI pos X X XX pair X X XX rcomp X X XX

Flags X X XX (flags) rlen X X XX nmisX snppX snptX indpX indtX indcX

Reads belonging to class P are characterized and can only be fully reconstructed by position, reverse complement and offset information between pair members if they were obtained using paired sequencing technology, by some flags and by the length of the read.

The next section details how these descriptors are defined.

Position descriptor layer.

In the position (pos) layer, only the mapping position of the first encoded read is stored in the AU header as an absolute position in the reference genome. All other position descriptors take on a value expressing the difference from the previous position. This modeling of the source of information, determined by the sequence of read positions, is generally characterized by reduced entropy, especially for sequencing processes that give results with high overlap.

For example, in FIG. 1 shows how, after describing the start position of the first alignment as position 10000 in the reference sequence, the position of the second read starting at position 10180 is described as 180. At high overlap values (>50x), most position vector descriptors will have a very high occurrence of low values, such as 0 and 1, and other small integers. In FIG. Figure 9 shows how the positions of three read pairs are described in the pos layer.

The reverse complement descriptor layer.

Each read of a pair of reads obtained using sequencing technologies can originate from any strand of the genome of the sequenced organic sample. However, only one of the two strands is used as a reference sequence.

In FIG. Figure 2 shows how, in a pair of reads, one read (read 1) can come from one thread, and the other (read 2) from the other.

When strand 1 is used as the reference sequence, reid 2 can be encoded as the reverse complement of the corresponding fragment on strand 1. This is shown in FIG. 3.

- 9 040022

In the case of linked reads, four combinations of direct and reverse complement pairs are possible.

This is shown in FIG. 4. The rcomp layer encodes these four possible combinations.

The same coding is used for reverse complement information for a read belonging to classes P, N, M, I. To provide enhanced selective access to data, the reverse complement information for a read belonging to these four classes is encoded in different layers, as shown. in table. 2.

Mating information descriptor layer.

The pairing descriptor is stored in the pair layer. Such a layer stores descriptors that encode the information needed to reconstruct the original pairs of reads when the sequencing technology used generates reads by pairs. While at the time of this disclosure, the vast majority of sequencing data is generated using paired read technology, this is not the case for all technologies. For this reason, the presence of this layer is not necessary to recover all information from the sequencing data if the sequencing technology of the genomic data in question does not generate pairwise read information.

Definitions.

Pair partner: A read that is related to another read in a pair of reads (for example, read 2 is the pair of read 1 in the previous example).

Pairing Distance: The number of nucleotide positions in the reference sequence that separate one position in the first read (pairing anchor, eg the last nucleotide of the first read) from one position in the second read (eg the first nucleotide of the second read).

Most Probable Pairing Distance (MPPD): The most probable pairing distance, expressed as the number of nucleotide positions.

Pairing Distance in Positions (PPD): PPD is a way of expressing the pairing distance in terms of the number of reads separating one read from the corresponding pair present in a specific position descriptor layer.

Most Probable Position Pairing Distance (MPPPD): The most likely number of reads separating one read from its pair present in the specific position descriptor layer.

Mating position error (PPE): is defined as the difference between the MPPD or MPPPD and the actual position of the partner in the pair.

Pairing anchor: the positions of the last nucleotide of the first read in a pair, used as a reference to calculate the distance of a pair of conjugated elements, expressed as the number of nucleotide positions or the number of positions read.

In FIG. Figure 5 shows how the mating distance between pairs of reads is calculated.

The descriptor layer of a pair is a pairing error vector calculated as the number of reads that must be skipped to reach the pair partner of the pair's first read, given a given pairing decoding distance.

In FIG. 6 shows an example of how the pairing errors are calculated, both as an absolute value and as a differential vector (characterized by lower entropy for high overlap values).

For information about the mating of reads belonging to classes N, M, P and I, the same descriptors are used. In order to implement selective access to different data classes, the pairing information for the reads belonging to the four classes is encoded in another layer, as shown in the figure.

Mating information in case of reads mapped by different references.

In the process of mapping a sequence read to a reference sequence, it is not uncommon for the first read in a pair to be mapped to one reference (eg, chromosome 1) and the second read to a different reference (eg, chromosome 4). In this case, the mating information described above should be combined with additional information related to the reference sequence used to map one of the reads. This is achieved by encoding the following parameters:

1) a reserved value (flag) indicating that the pair is mapped on two different sequences (different values indicate whether read 1 or read 2 is mapped to a sequence that is not currently encoded);

2) a unique reference identifier referring to the reference identifiers encoded in the main header structure, as described in Table 1. 1;

3) a third element containing mapping information on the reference identified at point 2 and expressed as an offset from the last encoded position.

In FIG. 7 shows an example of this scenario.

In FIG. 7, since read 4 does not map to the currently encoded reference sequence, the genomic encoder conveys this information by creating additional descriptors in the pair layer. In the example below, read 4 of pair 2 is mapped to reference #4 while the currently encoded reference is #1. This information is encoded using 3 components:

- 10 040022

1) One special reserved value is encoded as the pairing distance (in this case 0xffffff)

2) The second descriptor contains the reference id as specified in the main header (in this case - 4)

3) The third element contains mapping information in the corresponding reference (170).

Mismatch descriptors for class N reads.

Class N includes all reads in which only the mismatches represented by N are present at the location of the definition of bases A, C, G or T. All other bases of the read perfectly match the reference sequence.

In FIG. Figure 8 shows how the positions of the N mismatches in read 1 are encoded as follows:

absolute position in read 1; or a differential position relative to the previous N in the same read.

The positions of the N mismatches in read 2 are coded as follows:

absolute position in read 2 + length of read 1 OR;

differential position relative to the previous N.

In the nmis layer, the encoding of each pair of reads ends with a special delimiter character.

In FIG. 8 shows how an N''-mismatch (where at a given mapping position N is present in the read instead of the actual base in the reference sequence) is encoded as a mismatch position only:

1) relative to the beginning of the read; or

2) relative to the previous mismatch (differential coding).

Descriptors encoding substitutions (mismatches or SNPs), insertions and deletions.

A substitution is defined as the presence in the mapped read of a nucleotide base other than that present in the reference sequence at the same position.

In FIG. Figure 9 shows examples of substitutions in a mapped pair of reads. Each replacement is encoded as a position (snpp layer) and a type (snpt layer). Depending on the statistical occurrence of substitutions, insertions or deletions, the given source models of the associated descriptors and the generated symbols encoded in the associated layer may differ.

Source model 1: replacements as positions and types.

Substitution position descriptors.

The replacement positions are calculated in the same way as the nmis layer values, i.e.

in read 1, replacements are encoded as an absolute position in read 1; or as a differential position relative to a previous replacement in the same read.

In read 2, the replacements are encoded as absolute position in read 2 + length of read 1; or as a differential position relative to the previous replacement

In FIG. 10 shows how substitutions (where, at a given mapping position, the character in the read differs from the character in the reference sequence) are encoded as:

1) the position of the mismatch relative to the beginning of the read; or relative to the previous mismatch (differential encoding),

2) the type of mismatch represented as a code calculated as described in FIG. 10.

In the snpp layer, the encoding of each pair of reads ends with a special delimiter character. Replacement type descriptors.

For class M (and I, as described in the following sections), mismatches are coded by an index (moving from right to left) from the actual character present in the reference to the corresponding replacement character present in the read {A, C, G, T, N, Z}. For example, if an aligned read shows a C instead of a T that is present at the same position in the reference, the mismatch index will be indicated as 4. The decoding process reads the encoded syntax element, the nucleotide at the given position on the reference, and moves from left to right, returning the decoded character. For example, a 2 obtained for a position where G is present in the reference would be decoded as N. In FIG. 11 shows all possible substitutions and corresponding coding symbols. Obviously, to minimize descriptor entropy, each replacement index can be assigned different and context-adaptive probability models according to the statistical properties of each replacement type for each data class.

In the case of IUPAC ambiguity codes, the replacement mechanism is exactly the same, but the replacement vector is expanded as follows: S = {A, C, G, T, N, Z, M, R, W, S, Y, K, V, H, D, B}.

In FIG. 12 shows an example of encoding replacement types in the snpt layer.

Some examples of substitution coding when IUPAC disambiguation codes are used are shown in FIG. 13. Another example of substitution indices is shown in FIG. 14.

Deletion insertion coding.

For class I, mismatches and deletions are encoded by (with moving from right to left encoding) a replacement from the actual character present in the reference to the corresponding replacement character present in the read: {A, C, G, T, N, Z}. For example, if an aligned read shows a C instead of a T present in the same position in the reference, the mismatch index will be 4. In the event that the read shows a deletion where an A is present in the reference, the encoded character will be 5. The decoding process reads the encoded syntax element, nucleotide at a given position on the reference and moves from left to right, returning the decoded symbol. For example, a 3 obtained for a position where G is present in the reference will be decoded as Z. Insertions are encoded as 6, 7, 8, 9,10 respectively for inserted A, C, G, T, N.

In FIG. Figure 15 shows an example of how to encode substitutions, insertions and deletions in a class I read pair. , C, G, T, N, Z, M, R, W, S, Y, K, V, H, D, B} as described in the previous paragraph for mismatches. In this case, the insertion codes must have different values, namely 16, 17, 18, 19, 20, if the replacement vector has 16 elements. This mechanism is shown in Fig. 16.

Source model 2: one layer for each type of substitutions and indels.

For some statistics, an encoding model for substitutions and indels can be developed that is different from that described in the previous section, resulting in a source with lower entropy. This coding model is an alternative to the techniques described above for mismatches only, and also for mismatches and indels.

In this case, for each possible replacement character, one data layer is defined (5 without IUPAC codes, 16 with IUPAC codes), plus one layer for deletions and 4 more layers for insertions. For ease of explanation, and not as a limitation on the application of the model, the following description will focus on the case where IUPAC codes are not supported.

In FIG. 17 shows how each layer contains the position of mismatches or insertions of the same type. If there are no mismatches or insertions for that type in an encoded pair of reads, a 0 is encoded in the corresponding layer. To enable the decoder to start the decoding process for the layers described in this section, the header of each access block contains a flag signaling the first layer to be decoded. In the example in FIG. 18, the first decodable element is position 2 in layer C. When there are no mismatches or indel indices of this type in a pair of reads, 0 is added to the corresponding layers. On the decoding side, when the decoding pointer for each layer points to 0, the decoding process proceeds to the next a couple of reads.

Coding of additional signal flags.

Each data class presented above (P, M, N, I) may require encoding additional information about the nature of the encoded reads. This information may be associated, for example, with a sequencing experiment (eg, an indication of the probability of duplication of one read) or may express some characteristic of the mapping of the read (eg, the first or second in a pair). In the context of this invention, this information is encoded in a separate layer for each data class. The main advantage of this approach is the ability to selectively access this information only when necessary and only in the required area of the reference sequence. Other examples of the use of such flags are paired read;

reed mapped in the right pair;

rider or partner is not charted;

read or partner in a reverse chain pair;

first / second in a pair;

not primary alignment;

read does not pass platform/provider quality checks;

the read is a PCR or optical duplicate;

additional alignment.

Adaptation of reference sequences.

Mismatches coded for classes N, M, and I can be used to generate modified references that will be used to recode reads in the N, M, or I layer (relative to the first reference sequence, R0) as p-reads relative to the adapted R1 genome. For example, if we denote by r_in ^A M ie the class M read containing mismatches relative to the reference genome n, then after adaptation we can get r_in ^A M=r_(i(n+1)) ^A P with

A(Refn)=Refn+1, where A is the conversion from reference sequence n to reference sequence n+1.

In FIG. 19 shows how reads containing mismatches (M-reads) relative to reference sequence 1 (RS1) can be converted into perfectly matched reads (P-reads) relative to reference sequence 2 (RS2) obtained from RS1 by modification mismatched positions. This transformation can be expressed as

RS2 = A(RS1)

If the conversion expression A performed from RS1 to RS2 requires fewer mismatch expression bits present in M-reads, this encoding method results in less information entropy and hence better compression.

Source models, entropy encoders and coding modes.

For each layer of the genomic data structure disclosed in this invention, different coding algorithms can be adopted according to the specific characteristics of the data or metadata carried by the layer and its statistical properties. An encoding algorithm should be understood as the association of a particular descriptor source model with a particular entropy encoder. A particular source model can be defined and selected to obtain the most efficient data encoding in terms of minimizing source entropy. The choice of entropy encoder may be driven by coding efficiency considerations and/or by features of the probability distribution and associated implementation issues. Each choice of a particular encoding algorithm will be referred to as an encoding mode applied to the entire layer.

Each model associated with an encoding mode is characterized by a definition of the syntax elements generated by each source (eg, read positions, read pairing information, mismatches relative to a reference sequence, etc.);

defining an associated probabilistic model;

defining an associated entropy encoder.

Other benefits.

This classification makes it possible to implement efficient coding modes that use a lower entropy of the information source, characterized by the modeling of sequences of syntactic elements by separate independent data sources (for example, distance, position, etc.). Another advantage of the invention is the ability to access only a subset of the data types of interest. For example, one of the most important applications in genomics is to find differences in a genomic sample relative to a reference (SNV) or population (SNP). Today, this type of analysis requires the processing of full sequence reads, whereas when accepting the data representation disclosed by the invention, mismatches are already allocated to only one or three data classes (depending on the interest in considering N-codes and indels).

Another advantage is the ability to perform efficient transcoding from data and metadata compressed with reference to a specific reference sequence to a different reference sequence when a new reference sequence is published or when re-mapping is performed on previously mapped data (eg, using a different mapping algorithm).

In FIG. 20 shows an encoder 207 in accordance with the principles of this invention. Encoder 207 takes as input raw sequence data 209, such as those generated by genome sequencing machine 200. Encoder 207 takes as input raw sequence data 209, such as generated by genome sequencing machine 200. Genome sequencing machine 200 is known in the art, such as Illumina HiSeq 2500 or Thermo-Fisher Ion Torrent devices. The raw sequence data 209 is passed to an alignment module 201 which prepares the sequences for encoding by aligning the reads with the reference sequence. Alternatively, the de-novo assembler 202 can be used to generate a reference sequence from the available reads by looking for overlapping prefixes or suffixes so that longer segments (called contigs) can be assembled from the reads. After being processed by the de-novo 202 assembler, the reads can be mapped onto the resulting longer sequence. The aligned sequences are then classified by data classification module 204. The data classes 208 are then fed to layer encoders 205-207. The genomic layers 2011 are then fed into arithmetic encoders 2012-2014, which encode the layers according to the statistical properties of the data or metadata contained in that layer. As a result, the genomic stream 2015 is obtained.

In FIG. 21 shows a decoding device 218 in accordance with the principles of this invention. The decoder 218 receives the multiplexed genomic bitstream 2110 from the network or storage element. The multiplexed genomic bit stream 2110 is fed to a demultiplexer 210 to create separate streams 211, which are then fed to entropy decoders 212-214 to obtain genomic layers 215. The extracted genomic layers are fed to layer decoders 216-217 for further decoding of the layers into data classes. Class 219 decoders further process the genome descriptors and combine the results to produce uncompressed sequence reads, which can then be further stored in formats known in the art, such as a text file or a compressed zip file, or FASTQ files or SAM/BAM.

Class 219 decoders are capable of reconstructing original genomic sequences using information about the original reference sequences carried by one or more genomic streams. If reference sequences are not transported by genomic streams, they should be available on the decoding side and available to class decoders.

The methods of the invention disclosed herein may be implemented in hardware, software, firmware, or any combination thereof. When implemented in software, they may be stored on a computer medium and executed by a hardware processor. A hardware processor may comprise one or more processors, digital signal processors, general purpose microprocessors, application specific integrated circuits, or other discrete logic circuits.

The methods of the present invention may be implemented in a variety of devices or devices, including mobile phones, desktop computers, servers, tablets, and the like.

Claims (10)

CLAIM 1. A computer-implemented method for compressing genomic sequence data, wherein said genomic sequence data contains reads of nucleotide sequences, comprising the following steps: aligning said reads with one or more reference sequences, thereby obtaining aligned reads, classifying said aligned reads into different classes, including at least the first class: if said aligned reads match said one or more reference sequences without any errors ; second class: if the specified aligned reads match an area in the specified one or more reference sequences with a mismatch count consisting of the number of positions where the sequencer could not determine any base; third class: if said aligned reads match a region in said one or more reference sequences with a number of mismatches consisting of positions where the sequencer could not detect any bases and the presence of insertions or deletions or truncated nucleotides; fourth class: if said aligned reads do not find any valid mapping on said one or more reference sequences according to said alignment constraints, thus obtaining aligned read classes; and encoding said classified and aligned reads as a plurality of layers of syntactic elements containing descriptors, wherein for said first class, said descriptors include at least an initial position in the reference sequence, a flag signaling that the read should be considered as an inverse complement to the reference, distance to the partner of the pair in the case of paired reads, length value in the case when the sequencing technology produces reads of variable length; for said second class, the descriptors include at least the first class descriptors and a mismatch position for each mismatch; for said third class, the descriptors include the descriptors of said second class and the mismatch position and mismatch type for each mismatch; for the specified fourth class, the descriptors include the descriptors of the specified first class and the mismatch type for each mismatch, where the specified encoding of the specified classified aligned reads as a plurality of layers of syntax elements includes the selection of the specified syntax elements containing the specified descriptors, in accordance with the specified classes of aligned reads, and the encoding of the specified classified aligned reads in the form of a plurality of layers of syntactic elements is carried out using a specific entropy encoder 2012-2014. 2. The method according to claim 1, characterized in that the layers of syntax elements additionally contain the position of the variant relative to the reference sequence, the type of the variant, the position of the deletion relative to the reference sequence, the position of one or more characters that are absent in the reference sequence but are present in aligned reads, the type inserts in this position. 3. The method of claim 1, wherein said entropy encoder is a context-adaptive arithmetic encoder. 4. A method for decompressing a genomic stream compressed by the method of claim 1, said method comprising the following steps: parsing and decoding 212-214 of the compressed genomic stream into genomic layers of syntactic elements 215, decoding said genomic layers into data classes 216-217, unfolding said genomic layers into classified nucleotide sequence reads 14 040022 leotides, selective decoding using decoders 219 of classes of said classified reads of nucleotide sequences; and combining the results on one or more reference sequences to obtain uncompressed reads of nucleotide sequences. 5. Genomic encoder 2010 for compressing genomic sequence data 209, said genomic sequence data 209 comprising reads of nucleotide sequences, said genomic encoder 2010 comprising an alignment module 201 configured to align said reads with one or more reference sequences, thereby obtaining aligned reads, a data classification module 204 configured to classify said aligned reads into different classes, including at least the first class: if said aligned reads match said one or more reference sequences without any errors; second class: if the specified aligned reads match an area in the specified one or more reference sequences with a mismatch count consisting of the number of positions where the sequencer could not determine any base; third class: if said aligned reads match a region in said one or more reference sequences with a number of mismatches consisting of positions where the sequencer could not detect any bases and the presence of insertions or deletions or truncated nucleotides; fourth class: if said aligned reads do not find any valid mapping on said one or more reference sequences according to said alignment constraints, thus obtaining aligned read classes; one or more layer-encoding modules 205-207 configured to encode said classified aligned reads as syntax element layers containing descriptors by selecting said syntax elements according to said aligned read classes, where for said first class, said descriptors include at least initial position in the reference sequence, a flag indicating that the read should be considered as an inverse complement to the reference, the distance to the partner of the pair in the case of paired reads, the length value in the case when the sequencing technology produces reads of variable length; for said second class, the descriptors include at least the first class descriptors and a mismatch position for each mismatch; for said third class, the descriptors include the descriptors of said second class and the mismatch position and mismatch type for each mismatch; for the specified fourth class, the descriptors include the descriptors of the specified first class and the type of mismatch for each mismatch, entropy encoder 2012-2014 for entropy coding of the specified layers of syntactic elements. 6. A genomic decoder 218 for decompressing a genomic stream 211 compressed by the genomic encoder of claim 6, said genomic decoder 218 comprising parsing and decoding means 210, 212-214 configured to parse said compressed genomic stream into genomic layers of syntax elements 215 , one or more layer decoders 216-217 configured to decode genomic layers into data classes, and further configured to process said genomic layers into classified nucleotide sequence reads 2111, genomic data class decoders 213 configured to selectively decode said classified nucleotide sequence reads, and combining the result for one or more reference sequences to obtain uncompressed reads of nucleotide sequences. 7. Genomic decoder according to claim 6, characterized in that one or more reference sequences are stored in the compressed genome stream 211. 8. Genomic decoder according to claim 6, characterized in that one or more reference sequences are fed into the decoder by an out-of-band mechanism. 9. Genomic decoder according to claim 6, characterized in that one or more reference sequences are built in the specified decoder. 10. A computer-readable medium containing instructions that, when executed, cause at least one processor to perform the method of any one of claims 1-4.