EA043338B1 - METHOD AND DEVICE FOR COMPACT REPRESENTATION OF BIOINFORMATION DATA USING SEVERAL GENOMIC DESCRIPTORS - Google Patents

Description

Cross reference to related applications

This application claims priority and benefit to the following patent applications:

PCT/US2017/017842 filed February 14, 2017 and PCT/US2017/041591 filed July 11, 2017.

Field of technology

The present invention provides a new method for presenting genome sequencing data that reduces memory footprint and improves access characteristics by providing new functionality that is not available in prior art presentation methods.

State of the art

Proper representation of genome sequencing data is fundamental to enable efficient genomic analysis applications such as genomic variant detection and all other analyzes performed for various purposes by processing sequencing data and metadata.

Sequencing of the human genome has become available with the advent of low-cost, high-throughput sequencing technologies. This capability opens new perspectives in several areas, from cancer diagnosis and treatment to the detection of genetic diseases, from pathogenic surveillance to the identification of antibodies to the creation of new vaccines, drugs and the adaptation of personalized treatment procedures.

Hospitals, genomic data analytics providers, bioinformaticians and large biological data centers are seeking low-cost, fast, reliable and interconnected genomic solutions to enable the advancement of genomic medicine on a global scale. As data storage has become one of the bottlenecks in the sequencing process, increasing attention is being paid to ways to present genome sequencing data in a compressed form.

The most commonly used representations of genomic information from sequencing data are based on the FASTQ and SAM archiving formats. The goal is to compress traditionally used file formats (FASTQ and SAM for unaligned and aligned data, respectively). Such files are composed of plain text characters and compressed, as mentioned above, using general purpose approaches such as LZ schemes (from Lempel and Ziv, the names of the authors who published the first versions), such as the well-known zip, gzip, etc. When using general-purpose compression tools such as gzip, the compression result is typically a blob. The resulting information in such a monolithic form is quite difficult to archive, transmit and process, particularly as in the case of high-throughput sequencing where the volume of data is extremely large. The BAM format has poor compression performance because it focuses on compressing the inefficient and redundant SAM format rather than extracting the actual genomic information contained in the SAM files, and adapts general purpose text compression algorithms such as gzip rather than taking advantage of the specific nature of each source data (genomic data itself).

A more complex approach to genomic data compression, which is used less frequently but is more efficient than BAM, is CRAM. CRAM provides more efficient compression by implementing differential reference encoding (specifically, exploiting data source redundancy), but it still lacks features such as incremental updates, support for streaming data processing, and selective access to specific classes of compressed data. These approaches provide poor compression ratios and data structures that are difficult to navigate and manipulate once compressed. Subsequent analysis can be severely slowed by the need to process large and hard-coded data structures to even perform a simple operation or access selected regions of a genomic dataset. CRAM builds on the concept of CRAM writing. Each CRAM entry represents one mapped or unmapped read by encoding all the elements needed to reconstruct it.

CRAM exhibits the following limitations and disadvantages, which are addressed and overcome by the invention described herein.

1. CRAM does not support data indexing and random access to data subsets with common specific characteristics. Data indexing is outside the scope of this specification (see section 12 of the CRAM v3.0 Specification) and is implemented as a separate file. In contrast, the inventive approach described herein uses a data indexing method that is integrated with the encoding process, and the indices are embedded in the encoded (ie, compressed) binary stream.

2. CRAM is built from core data blocks that can contain mapped reads of any type (perfectly matched reads, reads with only substitutions, reads with insertions or deletions (also called indels)). There is no understanding of the classification and grouping of reads into classes according to the mapping result relative to the reference sequence. Those. You need to check all the data even when searching for reads with certain characteristics. Such a og-

- 1 043338 limitation is removed with the help of this invention by classifying and dividing the data into classes before encoding.

3. CRAM is based on the principle of encapsulating each read in a CRAM record. This means that when searching for reads that differ in certain biological characteristics (for example, reads with substitutions but no indels, or perfectly mapped reads), you need to check every complete record.

In contrast, in the present invention there is a concept of data classes encoded in separate information blocks, and there is no concept of a record encapsulating each read. This allows you to more efficiently access sets of reads with specific biological characteristics (for example, reads with substitutions but no indels, or perfectly mapped reads) without having to decode each read (block of reads) to check its characteristics.

4. In a CRAM record, each field of the record is associated with a specific flag, and each flag must always have the same value due to the absence of the concept of context, since each CRAM record may contain some other type of data. This encoding mechanism introduces redundant information and prevents the use of efficient contextual entropy encoding.

Instead, in the present invention there is no concept of a flag indicating data, since this is essentially determined by the information block to which the data belongs. This results in a significant reduction in the symbols that need to be used and, as a consequence, a reduction in the entropy of the information source, which leads to more efficient compression. This improvement is possible through the use of different blocks, allowing the encoder to reuse the same character within each block with different meanings according to the context. In CRAM, each flag must always have the same value because there is no concept of contexts, and each CRAM entry can contain any type of data.

5. In CRAM, substitutions, insertions and deletions are represented using different descriptors, and it is these choices that increase the size of the alphabet of the information source and give increased source entropy. In contrast, the approach of the invention described here uses one alphabet and one encoding for substitutions, insertions and deletions. This simplifies the encoding and decoding process and creates a low-entropy source model whose encoding produces highly compressible binary data streams.

The purpose of the present invention is to compress genomic sequences by classifying and partitioning sequencing data so that redundant information to be encoded is minimized and functions such as selective access and support for incremental updates operate directly on the compressed region.

One aspect of the presented approach is the definition of data and metadata classes, structured in different blocks and encoded separately. More relevant improvements to this approach compared to existing methods are as follows:

1) improving compression performance by reducing the entropy of the information source by providing an efficient source model for each class of data or metadata;

2) the ability to selectively access parts of the compressed data and metadata for the purpose of any further processing directly in the compressed area;

3) the ability to incrementally (i.e., without the need for decoding and re-encoding) compressed data and metadata with new sequencing data and/or metadata and/or new analysis results associated with specific sets of sequencing reads.

Brief description of drawings

In fig. Figure 1 shows how the position of mapped read pairs is encoded in the pos block as a difference relative to the absolute position of the first mapped read.

In fig. Figure 2 shows how two reads in a pair can come from two strands of DNA.

In fig. Figure 3 shows how the reverse complement of read 2 is encoded if strand 1 is used as a reference.

In fig. Figure 4 shows the four possible combinations of reads that make up a read pair and the corresponding encoding in the rcomp block.

In fig. Figure 5 for three pairs of reads shows how to calculate the pairing distance in the case of reads of constant length.

In fig. 6 shows how pairing errors encoded in the pair block allow the decoder to recover the correct pairing using the encoded MPPPD.

In fig. Figure 7 shows the coding of the pairing distance when a read is mapped to a different reference than its paired read. In this case, additional descriptors are added to the mapping distance. One is the signaling flag, the second is the reference ID and another is the pairing distance.

In fig. Figure 8 shows the encoding of an n-type mismatch in the nmis block.

In fig. Figure 9 shows a mapped read pair that contains substitutions compared to the reference sequence.

- 2 043338

In fig. Figure 10 shows how to calculate substitution positions as absolute or difference values.

In fig. Figure 11 shows how to evaluate symbols encoding substitution types when codes are used

IUPAC. The symbols represent the distance, in the circular substitution vector, between a molecule present in the read and a molecule present on the reference at that position.

In fig. Figure 12 shows how to encode substitutions in a snpt block.

In fig. Figure 13 shows how to calculate replacement codes when IUPAC ambiguity codes are used.

In fig. Figure 14 shows how to encode a snpt block when IUPAC ambiguity codes are used.

In fig. Figure 15 shows how for class I reads a replacement vector is used, the same as for class M, with the addition of special codes for inserting characters A, C, G, T, N.

In fig. Figure 16 shows some examples of coding mismatches and indels for IUPAC ambiguity codes. In this case, the replacement vector is much longer, and therefore there are more possible characters to be calculated than in the five-character case.

In fig. Figure 17 shows another source model for mismatches and indels, where each block contains the position of mismatches or insertions of the same type. In this case, no characters are encoded for the mismatch or indel type.

In fig. Figure 18 shows an example of coding mismatches and indels. When a read contains no mismatches or indels of a given type, a 0 is encoded in the corresponding block. 0 Acts as a read delimiter and terminator in each block.

In fig. Figure 19 shows how a change in the reference sequence can transform M reads into P reads. This operation can reduce the information entropy of the data structure, especially in the case of high coverage data.

In fig. 20 shows a genomic encoder 2010 in accordance with one embodiment of the present invention.

In fig. 21 shows a genomic decoder 218 in accordance with one embodiment of the present invention.

In fig. Figure 22 shows how an internal reference can be constructed by clustering reads and assembling segments taken from each cluster.

In fig. Figure 23 shows how the strategy for building a reference is to save the most recent reads After a certain sorting has been applied to the reads (for example, in lexicographical order).

In fig. Figure 24 shows how a read belonging to the class of unmapped reads (class U) can be encoded using six descriptors stored or contained in corresponding blocks.

In fig. 25 shows how alternative encoding of reads belonging to class U, where the descriptor pos has a sign, is used to encode the mapping position of the read on the constructed reference.

In fig. Figure 26 shows how reference transformations can be applied to remove mismatches from reads. In some cases, reference transformations may result in the creation of new mismatches or a change in the type of mismatches found when compared to the reference before the transformation was performed.

In fig. Figure 27 shows how reference transformations can change the class membership of reads by eliminating all or a subset of mismatches (i.e., a read belonging to class M before the transformation is assigned to class P after the reference transformation has been applied).

In fig. Figure 28 shows how semi-mapped read pairs (HM class) can be used to fill unknown regions of a reference sequence by assembling longer contigs using unmapped reads.

In fig. 29 shows how N, M, and I data class encoders are configured using threshold vectors and form separate subclasses of N, M, and I data classes.

In fig. 30 shows how all data classes can use the same transformed reference for re-encoding, or each class N, M and I, or any combination thereof, can use a different transform.

In fig. 31 shows the header structure of a genomic dataset.

In fig. 32 shows the general structure of a master index table (MIT), where each row contains genomic intervals of several data classes P, N, M, I, U, HM and, in addition, pointers to metadata and annotations. Columns refer to specific positions on reference sequences associated with the encoded genomic data.

In fig. Figure 33 shows an example of one MIT line containing genomic intervals related to reads of class P. Genomic regions related to different reference sequences are separated by a special flag (S in this example).

In fig. 34 shows the general structure of a local index table (LIT) and how it is used to store pointers to the physical location of encoded genomic information in stored or transmitted data.

In fig. 35 shows an example of a LIT used to access access blocks No. 7 and 8 in the block payload.

In fig. 36 shows the functional relationship between several MIT and LIT lines contained in the genomic block headers.

In fig. 37 shows how an access block is composed of several blocks of genomic data carried by different genomic streams containing data belonging to different classes. Each block further contains data packets used as data transmission blocks.

In fig. 38 shows how access blocks are composed of a header and multiplexed blocks belonging to one or more homogeneous data blocks. Each block may consist of one or more packets containing actual descriptors of genomic information.

In fig. Figure 39 shows multiple alignments without splicing. The leftmost read has N alignments. N is the first mmap descriptor value to be decoded and signals the number of alignments of the first read. The next N values of the mmap descriptor are decoded and used to calculate P, which is the number of alignments of the second read.

In fig. Figure 40 shows how the pos, pair, and mmap descriptors are used to encode multiple alignments without splices. The leftmost read has N alignments.

In fig. Figure 41 shows multiple splice alignments.

In fig. Figure 42 shows the use of pos, pair, mmap, and msar descriptors to represent multiple splice alignments.

The essence of the invention

The features of the following claims solve the problem of prior art solutions by providing a method for encoding genomic sequence data, wherein the genomic sequence data comprises nucleotide sequence reads, the method including the steps of aligning the reads to one or more reference sequences thereby creating aligned reads;

classifying the aligned sequences according to specified rules of comparison with one or more reference sequences, thereby creating classes of aligned reads;

encoding the classified aligned reads as a plurality of descriptor blocks, wherein encoding the classified aligned reads as a plurality of descriptor blocks includes selecting descriptors according to classes of the aligned reads;

structuring descriptor blocks using header information, thereby creating successive access blocks.

In accordance with yet another aspect, the encoding method further includes further classifying reads that do not satisfy the given matching rules into a class of unmapped reads;

constructing a set of reference sequences using some unmapped reads;

alignment of a class of unmapped reads to a set of constructed reference sequences;

encoding classified aligned reads as multiple descriptor blocks;

encoding a set of constructed reference sequences;

structuring descriptor blocks and encoded reference sequences using header information, thereby creating successive access blocks.

In accordance with yet another aspect, the encoding method further includes identifying genomic reads without mismatches to a reference sequence as a first class P.

In another aspect, the encoding method further includes identifying genomic reads as a second class N when mismatches are found only at positions where the sequencer fails to detect any base and the number of mismatches in each read does not exceed a given threshold.

In yet another aspect, the encoding method further includes identifying genomic reads as a third class M when mismatches are found at positions where the sequencer failed to detect any base (referred to as n-type mismatches) and/or it identified a different base compared to the reference sequence (called s-type mismatches), and the number of mismatches does not exceed the given thresholds for the number of n-type, s-type, and the threshold derived from a given function (f(n,s)) calculated on the number of mismatches n-type and s-type.

In accordance with yet another aspect, the encoding method further includes identifying genomic reads as a fourth class I when they may have mismatches of the same type as class M mismatches, and, additionally, at least one mismatch of the following type: insertion (i-type), deletion (d-type), soft cuts (c-type), and the number of mismatches of each type does not exceed the corresponding given threshold and the threshold provided by this function (w(n,s,i ,d,c)), calculated on the number of n-type, s-type, i-type, d-type and c-type mismatches.

In accordance with yet another aspect, the encoding method further includes identifying genomic reads as a fifth class U, which contains all reads not assigned to any of the classes P, N, M, I defined previously.

According to yet another aspect, the encoding method is further characterized in that the genomic sequence reads to be encoded are paired.

In accordance with yet another aspect, the encoding method is further characterized in that the classification further includes identifying the genomic sequences as a sixth HM class, as containing all pairs of reads where one read belongs to class P, N, M, or I, and the other read belongs to class U.

In accordance with yet another aspect, the encoding method further includes the steps of determining whether two paired reads are classified into the same class (each of P, N, M, I, U), then assigning the pair to the same identified Class.

Determining whether two paired reads belong to different classes, and if neither of them belongs to class U, then assigning the pair of reads to the class with the highest priority, determined according to the following priority:

P<N<M<I, where class P has the lowest priority, and class I has the highest priority;

determining whether only one of the paired reads is classified as belonging to class U, and classifying the pair of reads as belonging to class HM sequences.

In yet another aspect, the encoding method is further characterized in that each read class N, M, I is further divided into two or more subclasses (296, 297, 298) according to a threshold vector (292, 293, 294) defined by respectively, for each class N, M and I by the number of n-type mismatches (292), the function f(n,s) (293) and the function w(n,s,i,d,c) (294).

Determining whether two paired reads are classified into the same subclass, then assigning the pair to the same subclass;

determining whether two paired reads are classified into subclasses of different classes, then assigning the pair to the subclass belonging to the higher priority class according to the following expression:

N<M<I, where N has the lowest priority and I has the highest priority;

determining whether two paired reads are classified into the same class, which class is N, or M, or I, or into different subclasses, in which case assigning the pair to the highest priority subclass according to the following expressions:

Ni<N2<.. .<N _k ,

M ₁ <M ₂ <_<Mj, ^I 1< ^I ₂ <^< ^I h, where the highest index has the highest priority.

In accordance with yet another aspect, the mapping position information of each read is encoded by a pos descriptor block.

In yet another aspect, the stranding information (eg, which strand of DNA the read was sequenced from) of each read is encoded by an rcomp descriptor block.

In accordance with yet another aspect, information about paired-end reads is encoded by a pair descriptor block.

In another aspect, additional information about the alignment, such as whether the read is mapped to the correct pair, whether it has passed platform/vendor quality checks, whether it is a PCR duplicate or an optical duplicate, whether it is an auxiliary alignment, is encoded by a descriptor block flags.

In accordance with yet another aspect, information about unknown bases is encoded by a block of nmis descriptors.

In accordance with yet another aspect, information about the position of the substitutions is encoded by a snpp descriptor block.

In accordance with yet another aspect, information about the type of replacements is encoded by a block of special descriptors snpt.

In accordance with yet another aspect, information about the position of mismatches such as substitutions, insertions or deletions is encoded by an indp descriptor block.

In accordance with yet another aspect, information about the type of mismatches, such as substitutions, insertions or deletions, is encoded by an indt descriptor block.

In accordance with yet another aspect, the cut base information of the mapped read is encoded by an indc descriptor block.

In accordance with yet another aspect, information about unmapped reads is encoded by a ureads descriptor block.

In accordance with yet another aspect, information about the type of a reference sequence used for encoding is encoded by an rtype descriptor block.

In accordance with yet another aspect, information about multiple alignments of mapped reads is encoded by an mmap descriptor block.

In accordance with yet another aspect, information about spliced alignments and multiple alignments of the same read is encoded by an msar descriptor block and an mmap descriptor block.

In accordance with yet another aspect, information about read alignment scores is encoded by an mscore descriptor block.

In accordance with yet another aspect, information about the groups to which the reads belong is encoded by a block of special rgroup descriptors.

In accordance with yet another aspect, the encoding method is further characterized in that the descriptor blocks contain a main index table containing one partition for each class and subclass of aligned reads, wherein the partition contains mapping positions to one or more reference sequences of the first read of each access block of each class or data subclass; wherein the main index table and the access block data are encoded together.

In accordance with yet another aspect, the encoding method is further characterized in that the descriptor blocks also contain information related to the type of reference used (pre-existing or constructed) and segments of reads that do not match the reference sequence.

In yet another aspect, the encoding method is further characterized in that reference sequences are first converted into other reference sequences by applying substitutions, insertions, deletions, and truncations, and then the classified aligned reads are encoded as a plurality of descriptor blocks related to the converted reference sequences.

In accordance with yet another aspect, the encoding method is further characterized in that the same transformation is applied to the reference sequences for all data classes.

In accordance with yet another aspect, the encoding method is further characterized in that different transformations are applied to the reference sequences of each data class.

According to yet another aspect, the encoding method is further characterized in that the reference sequence transformations are encoded as descriptor blocks and structured by header information, thereby creating sequential access blocks.

In accordance with yet another aspect, the encoding method is further characterized in that encoding the classified aligned reads and associated reference sequence transformations into a plurality of descriptor blocks includes the step of associating a specific source model and a specific entropy encoder with each descriptor block.

According to yet another aspect, the encoding method is further characterized in that the entropy encoder is a context-adaptive arithmetic encoder, a variable length encoder, or a Golomb encoder.

The present invention also provides a method for decoding encoded genomic data, including the following steps:

parsing access blocks containing encoded genomic data to extract a plurality of descriptor blocks by using header information;

decoding multiple descriptor blocks to extract aligned reads in accordance with specified matching rules that determine their classification relative to one or more reference reads.

In accordance with yet another aspect, the decoding method further includes decoding unmapped genomic reads.

In accordance with yet another aspect, the decoding method further includes decoding classified genomic reads.

In accordance with yet another aspect, the decoding method further includes decoding a main index table containing one partition for each class of reads and associated corresponding mapping positions.

In accordance with yet another aspect, the decoding method further includes decoding information related to the type of reference being used—preexisting, transformed, or constructed.

In accordance with yet another aspect, the decoding method further includes decoding information related to one or more transformations to be applied to pre-existing reference sequences.

In another aspect, the decoding method is characterized in that the genomic reads are paired.

According to another aspect, the decoding method further includes a case where genomic data is entropy decoded.

The present invention also provides a genomic encoder (2010) for compressing genomic sequence data (209), wherein the genomic sequence data (209) comprises nucleotide sequence reads, wherein the genomic encoder (2010) includes an aligner assembly (201) configured to align the reads to one or more reference sequences, thereby creating aligned reads;

a constructed sequence generator unit (202), configured to create constructed reference sequences;

a data classification node (204), configured to classify aligned reads in accordance with specified matching rules with one or more pre-existing reference sequences or constructed reference sequences, thereby creating classes of aligned reads (208);

one or more block encoding nodes (205-207) configured to encode the classified aligned reads as descriptor blocks by selecting descriptors according to the classes of the aligned reads;

multiplexer (2016) for multiplexing compressed genomic data and metadata.

In accordance with yet another aspect, the genomic encoder further comprises a reference sequence transformation unit (2019) configured to transform pre-existing references and data classes (208) into transformed data classes (2018).

In accordance with yet another aspect, the genomic encoder is further characterized in that the data classification node (204) contains N, M, and I data class encoders configurable using threshold vectors, forming subclasses of N, M, and I data classes.

In accordance with yet another aspect, the genomic encoder is further characterized in that the reference transform node (2019) applies the same reference transform (300) to all classes and subclasses of the data.

In yet another aspect, the genomic encoder is further distinguished in that the reference transform decoder (2019) applies different transforms (301, 302, 303) to different classes and subclasses of data.

In accordance with yet another aspect, the genomic encoder further comprises functions suitable for implementing all aspects of the previously mentioned encoding methods.

The present invention also provides a genomic decoder (218) for decompressing a compressed genomic stream (211), wherein the genomic decoder (218) includes a demultiplexer (210) for demultiplexing the compressed genomic data and metadata;

parsing means (212-214) configured to parse the compressed genomic stream to produce genomic descriptor blocks (215);

one or more block decoders (216-217) configured to decode genomic blocks into classified nucleotide sequence reads (2111);

decoders (219) of classes of genomic data configured to selectively decode classified nucleotide sequence reads relative to one or more reference sequences to obtain uncompressed nucleotide sequence reads.

In accordance with yet another aspect, the genomic decoder further comprises a reference transform decoder (2113) configured to decode reference transform descriptors (2112) and create a transformed reference (2114) for use by genomic data class decoders (219).

In accordance with yet another aspect, the genomic decoder is further characterized in that one or more reference sequences are stored in a compressed genomic data stream (211).

In accordance with yet another aspect, the genomic decoder is further characterized in that one or more reference sequences are provided to the decoder via an out-of-band mechanism.

In accordance with yet another aspect, the genomic decoder is further characterized in that one or more reference sequences are constructed in the decoder.

In accordance with yet another aspect, the genomic decoder is further characterized in that one or more reference sequences are transformed at the decoder by a reference transform decoder (2113).

The present invention also provides a computer-readable storage medium comprising instructions whose execution causes at least one processor to perform all aspects of the previously mentioned encoding methods.

The present invention also provides a computer-readable storage medium comprising

- 7 043338 instructions, the execution of which causes at least one processor to perform all aspects of the previously mentioned decoding methods.

The present invention also provides support for storing genomic data encoded in accordance with all aspects of the previously mentioned encoding methods.

Detailed description

Genomic or proteomic sequences mentioned in this invention include, for example, but are not limited to, nucleotide sequences, deoxyribonucleic acid (DNA), ribonucleic acid (RNA) and amino acid sequences. Although genomic information in the form of a nucleotide sequence is described in some detail herein, it is clear that compression methods and systems can also be implemented for other genomic or proteomic sequences, albeit with some modifications understood by one skilled in the art.

Genome sequencing information is generated by high-throughput sequencing (HTS) instruments in the form of sequences of nucleotides (also called bases) represented by strings of alphabetic characters from a defined dictionary. The minimal vocabulary is represented by five characters: {A, C, G, T, N}, which represent the 4 types of nucleotides present in DNA, namely adenine, cytosine, guanine and thymine. In ribonucleic acid (RNA), thymine is replaced by uracil (U). N indicates that the sequencer was unable to detect any base, leaving the actual nature of the nucleotide at that position undetermined. If the sequencer accepts IUPAC ambiguity codes, the alphabet used for the characters is (A, C, G, T, U, W, S, M, K, R, Y, B, D, H, V, N or -). The nucleotide sequences produced by sequencers are called reads. Sequence reads can range in length from several tens to several thousand nucleotides. Some sequencing technologies generate sequence reads in pairs, with one read coming from one strand of DNA and the other from a different strand. In genome sequencing, the term coverage is used to express the level of redundancy of sequence data relative to a reference sequence. For example, to achieve 30x coverage of the human genome (3.2 billion bases long), the sequencer must generate a total of 30x3.2 billion bases so that on average each position in the reference is covered 30 times.

As used herein, a reference sequence is any sequence to which nucleotide sequences generated by sequencers are aligned/mapped. One example of a sequence would be a reference genome—a sequence collected by scientists as a representative example of a species' gene sets. For example, GRC37, the Human Genome Reference Consortium genome (build 37), was derived from samples from thirteen anonymous volunteers in Buffalo, New York. However, a reference sequence may also consist of a synthetic sequence designed and simply constructed to improve the compressibility of the reads for further processing. This is described in more detail in the section U class descriptors and construction of internal references for unmapped U class and HM class reads and is depicted in Fig. 22 and 23.

Sequencers can introduce errors into sequence reads, such as the following.

1. The decision to omit the definition of a basis due to lack of confidence in the definition of any specific basis. This is called an unknown base and is labeled N (denoted as an n-type mismatch).

2. Using an incorrect symbol (i.e., representing a different nucleic acid) to represent the nucleic acid actually present in the sequenced sample; this is usually called a substitution error (denoted as an s-type mismatch).

3. Insertion of additional characters into one sequence read that are not related to any actually present nucleic acid; this is usually called an insertion error (denoted as an i-type mismatch).

4. Removing from one read a sequence of characters that represent nucleic acids actually present in the sequenced sample; this is usually called a deletion error (denoted as a d-type mismatch).

5. Recombination of one or more fragments into one fragment that does not reflect the true essence of the original sequence; this usually results in the aligner deciding to truncate the bases (denoted as a c-type mismatch).

The term coverage is used in the literature to quantify the extent to which a reference genome or part of it can be covered by available sequence reads. Coverage is called partial (less than 1X) when not a single available sequence read is mapped to some parts of the reference genome;

single-fold (1X), when one and only one symbol represented in the sequence reads is mapped to all nucleotides of the reference genome;

- 8 043338 multiple (2X, 3X, NX), when each nucleotide of the reference genome is mapped several times.

The purpose of the present invention is to define a format for representing genomic information in which relevant information can be efficiently accessed and transported and in which the weight of redundant information is reduced. The disclosed invention has the following novel elements.

1. Sequence reads are classified and divided into data classes according to the alignment results relative to reference sequences. Such classification and separation allows selective access to the encoded data according to criteria related to the alignment results and the accuracy of the match.

2. Classified sequence reads and associated metadata are represented as homogeneous descriptor blocks to obtain clearly distinguishable sources of information with low information entropy.

3. The ability to model each individual source of information using a clearly distinguishable source model adapted to the statistical characteristics of each class, and the ability to change the source model within each class and within a descriptor block for each individually accessible data block (access block). Implementation of suitable context-adaptive probabilistic models and corresponding entropy encoders according to the statistical properties of each source model.

4. Define correspondences and dependencies between descriptor blocks to allow selective access to sequencing data and associated metadata without having to decode all descriptor blocks if not all information is required.

5. Encoding each sequence data class block and associated metadata relative to pre-existing (also called external) reference sequences or relative to transformed reference sequences obtained by applying appropriate transformations to pre-existing reference sequences to reduce the entropy of the descriptor block information sources. Descriptors represent reads divided into different data classes. After any encoding of reads using appropriate descriptors relative to a pre-existing reference or a transformed pre-existing reference sequence, the occurrence of various mismatches can be used to determine suitable transformations into reference sequences in order to find the final low-entropy encoded representation and achieve higher compression efficiency.

6. Construction of one or more reference sequences (also called internal references to distinguish them from pre-existing reference sequences, also called external reference sequences) used to encode a class of reads that demonstrate a degree of match accuracy relative to pre-existing reference sequences that does not satisfy specified restrictions. Such restrictions are established so that the cost of encoding a representation in compressed form of a class of reads aligned to internal reference sequences, and the cost of representing the internal reference sequences themselves, is lower than when encoding an unaligned class of representations literally or using external reference sequences without transformations or with transformations. In the following, each of the above-mentioned aspects will be further described in detail.

Classification of sequence reads according to matching rules.

Sequence reads generated by sequencers in accordance with the present invention are divided into six different classes depending on the results of the alignment with one or more pre-existing reference sequences.

When aligning a DNA nucleotide sequence relative to a reference sequence, the following cases can be distinguished.

1. The region in the reference sequence matches the sequence reads without a single error (i.e. perfect mapping). This nucleotide sequence is called a perfect match read or designated as class P.

2. The region in the reference sequence coincides with the sequence read, and the type and number of mismatches are determined only by the number of positions in which the sequencer that generated this read was unable to determine any base (or nucleotide). Mismatches of this type are designated by the letter N, indicating the unidentified nucleotide base. This type of mismatch is referred to herein as an n-type mismatch. Such sequences belong to reads of class N. Once a read has been classified as belonging to class N, it is useful to limit the degree of inaccuracy of the match to a given upper bound and to distinguish between matches that are considered valid and those that are not considered valid. Therefore, reads assigned to class N are also limited by setting a threshold (MAXN), which determines the maximum allowed number of reads.

- 9 043338 The number of undefined bases (i.e. bases defined as N) that a read may contain. This classification implicitly defines the required minimum match accuracy (or maximum mismatch rate) common to all reads belonging to class N when compared with the corresponding reference sequence, which provides a useful criterion for applying selective data retrieval to compressed data.

3. The region in the reference sequence matches the sequence read, with the types and numbers of mismatches determined by the number of positions at which the sequencer generating the read was unable to identify any nucleotide base, if any (i.e. n-type mismatch), plus the number of positions in which a base was defined that was different from that present in the reference. This type of mismatch, referred to as a substitution, is also called a single nucleotide variation (SNV) or single nucleotide polymorphism (SNP). This type of mismatch is also referred to herein as an s-type mismatch. In such a case, the sequence read is called an M-mismatch read and is assigned to class M. As with class N, for all reads belonging to class M, it is useful to limit the degree of match uncertainty to a given upper bound and to distinguish between matches that are considered valid and those that are not considered valid. Therefore, reads assigned to class M are also constrained by setting a set of thresholds, one for the number n of n-type mismatches (MAXN), if any, and another for the number of substitutions s (MAXS). The third constraint is a threshold defined by any function of both n and s, f(n,s). This third constraint allows classes to be formed with an upper bound on the inaccuracy of the match in accordance with any significant criterion for selective access. For example, without limitation, f(n,s) may be (n+s)1/2, or (n+s), or a linear or nonlinear expression that sets a bound for the maximum level of match imprecision assumed for a read belonging to the class M. Such a bound is a very useful criterion for applying the required selective data search to compressed data when analyzing sequence reads for various purposes, since it allows an additional bound to be set for any possible combination of numbers of n-type mismatches and s-type mismatches (substitutions). beyond the simple threshold applied to one type or another type.

4. The fourth class consists of sequence reads that represent at least one mismatch of any type among insertions, deletions (otherwise known as indels), and truncations, plus (if any) any mismatches of a type belonging to class N or M. Such sequences are called reads with I-mismatches and are assigned to class I. Insertions are formed by an additional sequence of one or more nucleotides not present in the reference but present in the read sequence. This type of mismatch is referred to herein as an i-type mismatch. In the literature, when the inserted sequence is at the edges of the sequence, it is also called soft-clipped (i.e., the nucleotides do not match the reference sequence, but remain in the aligned reads, unlike hard-clipped nucleotides, which are discarded). This type of mismatch is referred to herein as a c-type mismatch. Decisions to retain or discard nucleotides are made at the aligner stage rather than the read classifier disclosed herein, which accepts and processes the reads as determined by the sequencer or at the next alignment stage. Deletions are holes (missing nucleotides) in the read compared to the reference sequence. In this document, mismatches of this type are referred to as d-type mismatches. As in the case of classes N and M, it is possible and necessary to determine the limit of mismatch inaccuracy. The definition of a set of restrictions for class I is based on the same principles that are used for class M and are indicated in the last rows of the table. 1. In addition to the threshold for each type of mismatch allowed for class I data, an additional constraint is defined using a threshold defined by any function of the numbers of mismatches n, s, d, i and c - w(n,s,d,i,c ). This additional constraint allows classes to be formed with an upper bound on the imprecision of the match according to any meaningful user-defined selective access criterion. By way of example, and not limitation, w(n,s,d,i,c) may be (n+s+d+i+c)1/5, or (n+s+d+i+c), or any linear or non-linear expression that sets a bound on the maximum level of match uncertainty that is accepted for a class I read. Such a bound provides a very useful criterion for applying the required selective data search to compressed data when analyzing sequence reads for various purposes, since it allows an additional bound to be set for any possible combination of numbers of mismatches allowed in Class I reads, beyond the simple threshold applied to each type of mismatch allowed.

5. The fifth class includes all reads for which no mappings are found that are considered valid (i.e., do not satisfy the set of rules that define the upper bound for the maximum imprecision of the match, as indicated in Table 1) for each data class at comparison with the reference sequence. Such sequences are called unmapped to reference sequences and are classified as belonging to class U.

- 10 043338

Classification of pairs of reads in accordance with matching rules.

The classification defined in the previous section applies to single sequence reads. For sequencing technologies that generate reads in pairs (i.e., Illumina Inc. technologies) where the two component reads are known to be separated by an unknown sequence of variable length, it is appropriate to consider assigning the entire pair to a single data class. A read that is paired with another read is called a paired read.

If both paired reads belong to the same class, the assignment of the pair to the class is obvious: the entire pair is assigned to the same class as any of the parts of the pair (i.e. P, N, M, I, U). If two reads belong to different classes, but neither belongs to class U, then the entire pair is assigned to the class with the highest priority, determined according to the following expression:

P<N<M<I, where class P has the lowest priority and class I has the highest priority.

In the case when one of the reads belongs to class U, and its paired read belongs to any of the classes P, N, M, I, define the sixth class as class HM, which means half mapped.

The definition of such a specific class of reads is motivated by the fact that it is used in an attempt to identify gaps or unknown regions that exist in reference genomes (also called obscure or unknown regions). Such regions are recovered by mapping pairs at the edges using pairwise reads, which can be mapped to unknown regions. The unmapped read pairs are then used to construct so-called unknown region contigs, as shown in FIG. 28. Therefore, providing selective access to only pairs of reads of this type greatly reduces the associated computational load, allowing such data, originating from large numbers of data sets, to be processed much more efficiently than using prior art solutions that would require full verification.

The table below summarizes the matching rules applied to reads to determine the data class to which each read belongs. The rules are defined in the first five columns of the table in the form of the presence or absence of mismatch types (mismatch types n, s, d, i, and c). The sixth column provides rules in the form of a maximum threshold for each type of mismatch and any function (if any) f(n,s) and w(n,s,d,i,c) of the possible types of mismatch.

Table 1

The type of mismatch and set of constraints that must be satisfied by each sequence read assigned to the data classes defined in this description

Number and types of mismatches detected when comparing a read to a reference sequence A set of restrictions on the accuracy of the match Destination class Number of Unknown bases (“Ν”) Number of replacements Number of deletions Qty insertions Number of bases cut 0 0 0 0 0 0 R η > 0 0 0 0 0 η < ΜΑΧΝ N η > ΜΑΧΝ And η > 0 s>0 0 0 0 η < ΜΑΧΝ and s < MAXS and f(n,s) < МАХМ M η > ΜΑΧΝ or s > MAXS or f(n,s) > МАХМ And

- 11 043338

η > 0 s>0 d > 0* i > 0* c>0* η < ΜΑΧΝ and s < MAXS and d < MAXD and i < MAXI and c < MAXC w(n,s,d,i,c) < MAXTOT 1 *There must be at least one mismatch of type d, i, c (i.e. d > 0, or i > 0, or c > 0) d>0 i>0 c > 0 η > MAXN or s > MAXS or d > MAXD or i > MAXI or And c > MAXC w(n,s,d,i,c) > MAXTOT

Matching rules for dividing N, M, and I read data classes into subclasses with varying degrees of matching accuracy.

The data classes of types N, M and I, which are defined in the previous sections, can be further divided into an arbitrary number of different subclasses with different degrees of matching accuracy. This capability is an important technical advantage in providing partitioning into smaller subsets and, as a result, much more efficient selective access to each class of data. As an example, but not a limitation, to partition a class N into k subclasses (subclass N1, ..., subclass N _k ), it is necessary to define a vector with the corresponding components MAXN1, MAXN2, ..., MAXN(k-1), MAXN _W at provided that MAXN1<MAXN2<...<MAXN(k_ _D <MAXN, and assign each read to the lowest-ranking subclass that satisfies the constraints specified in Table 1 while performing an evaluation for each element of the vector. This is shown in FIG. 29, where the data classification block (291) contains a P, N, M, I, U, HM class encoder and encoders for annotations and metadata.The N class encoder is configured using a threshold vector 292 MAXN1-MAXN _k that generates k data subclasses N (296).

In the case of classes of type M and I, apply the same principle by defining a vector with the same properties for MAHM and MAHTOT, respectively, and use the components of each vector as a threshold to check whether the functions f(n,s) and w( n,s,d,i,c). As with subclasses of type N, assign to the lowest-ranking subclass that satisfies the constraint. The number of subclasses for each class type is independent and any combination of subdivisions is permitted. This is shown in Fig. 29, where the M class encoder 293 and the I class encoder 294 are configured, respectively, using the threshold vector MAXM ₁ -MAXM _j and MAXTOT1 -MAXTOT _h . These two encoders respectively generate j subclasses (297) of M data and h subclasses (298) of I data.

If two reads in a pair are assigned to the same subclass, then the pair belongs to the same subclass.

If two reads in a pair are classified as subclasses of different classes, then the pair belongs to the subclass of the class with higher priority according to the following expression:

N<M<I, where N has the lowest priority and I has the highest priority.

If two reads belong to different subclasses of one of the classes N, M or I, then the pair belongs to the subclass with the highest priority according to the following expressions:

N1<N2<_<Nk,

M1<M ₂ <_<Mj, l|<I2<—lh, where the highest index has the highest priority.

Transformations of external reference sequences.

Mismatches found for reads classified as N, M, and I can be used to create transformed references that will be used to more efficiently compress the read representation.

Reads classified as being in classes N, M, or I (relative to a preexisting (i.e., external) reference sequence designated RS0) may be coded relative to the converted reference sequence RS1 according to the occurrence of actual mismatches relative to the converted reference. For example, if a read ^M in belonging to class M (designated as the i-th read of class M) contains mismatches relative to the reference sequence RSn, then after converting read ^M in = read ^P _i ( _n+1 ) can be obtained

- 12 043338 using A(Ref _n )=Ref _n+1 , where A is the conversion from reference sequence RS _n to reference sequence RS _n+1 . In fig. Figure 19 shows an example of how reads containing mismatches (belonging to class M) relative to reference sequence 1 (RS1) can be converted into reads that perfectly match reference sequence 2 ( _RS2 ) derived from RS1 by changing the bases corresponding to the positions of the mismatch . They remain classified and are encoded along with other reads in the same access block, but the encoding is done using only the descriptors and descriptor values needed for the class P read. This transformation can be represented as the following expression:

RS2=A(RS1)

When the transform representation A that generates RS2, when applied to RS1, and the read representation relative to RS2 correspond to lower entropy than the class M read representation relative to RS1, it is appropriate to transmit the transform representation A and the corresponding read representation relative to RS2, since higher representation compression is achieved data.

Encoding Transform A for transmission in a compressed binary data stream requires defining two additional syntax syntax elements, as specified in the table below.

Descriptors Semantics Comments rftp Reference conversion position Positions of difference between the reference and the contig used for prediction rftt Conversion type reference Type of difference between reference and contig used for prediction The same syntax is described for the snpt descriptor defined later in this document.

In fig. Figure 26 shows an example of how to apply reference transformation to reduce the number of mismatches to be coded on mapped reads.

It is important to note the consequences of applying a transformation to a reference in some cases:

it is possible to introduce inconsistencies into the read representations that would not be present when compared with the reference before applying the transformation;

it is possible to change the types of nonconformities; a read may contain A instead of G while all other reads contain C instead of G, but the mismatches remain in the same position;

different data classes and data subsets of each data class may refer to the same transformed reference sequences or to reference sequences obtained by applying different transformations to the same pre-existing reference sequence.

In fig. Figure 27 also shows an example of how reads can change encoding type from one data class to another using the appropriate set of descriptors (for example, using class P descriptors to encode a read from class M) after applying a transformation and representing the read using the transformed reference. This occurs, for example, when a transformation changes all bases corresponding to read mismatches to bases actually present in the read, thereby effectively converting a read belonging to class M (if compared to the original untransformed reference sequence) into a virtual read of class P (if compare with the transformed reference). The definition of the set of descriptors used for each data class is given in the following sections.

In fig. Figure 30 shows how different data classes can use the same transformed reference R1=A ₀ (R0) (300) to re-encode reads, or different transformations A _N (301), AM ( can be applied to each data class separately 302), AI (303).

Determine the information needed to represent sequence reads as descriptor blocks.

Once the classification of reads is completed by defining classes, further processing consists of defining a set of different descriptors that represent the remaining information allowing the sequence of the read to be reconstructed when presented as mapping to a given reference sequence. The data structure of these descriptors requires storing global parameters and metadata that will be used by the decoding engine. This data is structured in a genomic dataset header, described in the table below. A dataset is defined as the collection of coding elements required to recover the genomic information associated with one genome sequencing run and all subsequent analysis. If the same sample for a genome is sequenced twice in two different runs, the resulting data will be encoded in two different data sets.

- 13 043338

Genomic dataset header structure

Table 1a

Element Type Description Unique ID Array Unique identificator bytes for coded content Main brand Array Main + bytes additional version Additional version Array bytes coding algorithm Header size Integer Size in bytes of total encoded content Read length Integer Size of reads in case of constant read length. Special meaning (eg 0) reserved for variable-length reads Ref. counter Integer Number of reference sequences used Access block counters Array Total bytes coded blocks (For example, access to reference whole numbers) subsequence Reference IDs Array Unique bytes identifiers for reference sequences for (i=0; i<Ref_count; i++) { Reference_genome: RefJD String: string A unique ID as a character string, defining the reference sequences used in a given data set ί for (i=0; i<Ref_count; i++) { Ref blocks Array Number of coded bytes blocks for each reference } Data set label size Integer Next size element

- 14 043338

Data set label Line String of characters used to identify a data set Data set type Integer Type of data encoded in the dataset (e.g., aligned, unaligned) Main index table The alignment positions of the first read in each block (access block). That is, the lowest position of the first read on the reference genome for each block of six classes. 1 per pos class (six) per reference sequence Array bytes This is a multi-dimensional array that supports random access to access blocks Label list A subset of the main title, indicating • number of tags; • for each tag: o Tag ID; o number of reference sequences related to the label; o for each reference sequences: reference identifier; number of areas, covered by the label; for each area: • Class ID; • starting position in the genomic range; • terminal position in the genomic range. Array bytes This is a list of tags, each of which is represented as a multidimensional array to support selective access to specific genomic regions or subregions, or unions of areas or subareas The start position and end position can be replaced by "block numbers", together with the reference sequence ID and class ID, forming a three-dimensional vector to address the coordinates of the main index table Parameter set Array bytes Encoding options used for configuring the encoding process and sent to the decoder

A sequence read (i.e., a segment of DNA) related to a given reference sequence can be expressed using the following elements.

Starting position on the reference sequence (pos).

A flag indicating whether the sequence should be treated as reverse complement to the reference (rcomp).

Distance to a paired read in the case of paired reads.

- 15 043338

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

For each mismatch, position (nmis for class N, snpp for class M, and indp for class I);

type of mismatch (absent in class N, snpt in class M, indt in class I).

Flags indicating certain characteristics of sequence reads, such as a pattern having multiple segments when sequenced;

each segment properly aligned according to the alignment means; unmapped segment;

next unmapped segment in template;

first or last segment signaling;

failure to pass quality control;

PCR duplicate or optical duplicate;

assisted alignment;

additional alignment.

A string of soft-clipped nucleotides, if present (indc in class I).

A flag indicating the reference used for alignment and compression (e.g. internal reference for class U), if applicable (rtype descriptor).

For class U, the indc descriptor indicates those parts of the reads (usually edges) that do not meet a given set of restrictions on the accuracy of matches with internal references.

The ureads descriptor is used to literally encode reads that cannot be mapped to any available reference, be it a pre-existing (i.e. external, such as the actual reference genome) or internal reference sequence.

This classification creates groups of descriptors (descriptors) that can be used to unambiguously represent genomic sequence reads. The table below summarizes the descriptors required for each class of reads aligned to external (i.e. pre-existing) or internal (i.e. constructed) references.

table 2

Defined descriptor blocks for each data class

P N M 1 And H.M. pos X X X X X X pair X X X X rcomp X X X X X X flags X X X X X X rlen X X X X X X nmis X snpp X X snpt X X indp X X indt X X indc X X X ureads X X rtype X rgroup X X X X X X mmap X X X X X msar X X X X X mscore X X X X X

Reads belonging to the P class are distinguished by the fact that they can be perfectly reconstructed using only position, reverse complement information and offset between paired reads in the case when they were obtained using sequencing technology that provides paired-end reads with a long insert, some flags and read length.

The next section details the definition of these descriptors for classes P, N, M and I, and for class U they are described in a subsequent section.

The HM class applies only to pairs of reads and represents a special case in which one read belongs to class P, N, M, or I, and the other belongs to class U.

Position descriptor.

- 16 043338

The position block (pos) stores only the mapping position of the first mapped read as an absolute value on the reference sequence.

All other position descriptors take a value expressing the difference relative to the previous position. This modeling of the source of information, determined by the sequence of read position descriptors, typically has reduced entropy, especially for sequencing processes that generate high-coverage results.

For example, in FIG. Figure 1 shows how, after describing the starting position of the first alignment as position 10,000 on the reference sequence, the position of the second read, starting at position 10,180, is described as 180. At high coverages (>50x), most position vector descriptors show a very high occurrence of low values, such like 0 and 1 and other small integers. In fig. Figure 1 shows how the positions of three paired reads are described in the pos block.

Descriptor for reverse complement.

Each read of a pair of reads generated by sequencing technologies can come from either of the two strands of the genome of the organic sample being sequenced.

However, only one of the two chains 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 strand, and the other read (read 2) can come from a different strand.

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

In the case of paired reads, four combinations of pairs are possible from forward and reverse pairwise complement pairs. This is shown in Fig. 4. The rcomp block encodes the four possible combinations.

The same encoding is used for the reverse complement information of reads belonging to classes N, M, P and I. To provide selective access to different classes of data, the reverse complement information of reads belonging to the four classes is encoded in different blocks, as shown in Table . 2.

Mating information descriptor.

The pairing descriptor is stored in a pair block. Such a block stores descriptors that encode the information necessary to reconstruct the resulting pairs of reads when the sequencing technology used creates reads in pairs. Although at the time of disclosure the vast majority of sequencing data are generated using paired-end read technology, not all technologies are. That is why the presence of this block is not necessary to restore all the information from the sequencing data if the genomic data sequencing technology in question does not provide information about paired reads.

Definitions.

Paired read: a read that is associated with another read in a read pair (for example, in the previous example, read 2 is a paired read of read 1);

pairing distance: the number of nucleotide positions of the reference sequence that separate one position in the first read (pairing anchor, e.g., the last nucleotide of the first read) from one position in the second read (e.g., the first nucleotide of the second read);

most probable pairing distance (MPPD): this is the most probable pairing distance expressed in the number of nucleotide positions;

position pairing distance (PPD): PPD is a way of expressing the pairing distance in terms of the number of reads separating one read from its corresponding paired read present in a particular position descriptor block;

Most Probable Position Pairing Distance (MPPPD): This is the most likely number of reads separating one read from its paired read present in a particular position descriptor block;

position pairing error (PPE): defined as the difference between MPPD or MPPPD and the actual position of the paired read;

pairing anchor: The position of the last nucleotide of the first read in a pair, used as a reference point for calculating the distance to the paired read in number of nucleotide positions or number of read positions.

In fig. Figure 5 shows how the pairing distance among pairs of reads is calculated.

A block of pair descriptors is a vector of pairing errors, calculated as the number of reads that must be skipped to reach the paired read of the first read of the pair relative to a certain decoding pairing distance.

In fig. Figure 6 shows an example of how pairing errors are calculated, both as an absolute value and as a difference vector (characterized by low entropy at high coverages).

The same descriptors are used for the pairing information of reads belonging to classes N, M, P and I. To provide selective access to different classes of data, the pairing information of reads belonging to these four classes is encoded in different blocks, as depicted in FIG. 8

- 17 043338 (class N), fig. 10, 12 and 14 (class M) and fig. 15 and 16 (Class I).

Pairing information for reads mapped to different reference sequences.

In the process of mapping sequence reads to a reference sequence, often the first read in a pair is mapped to one reference sequence (for example, chromosome 1), and the second is mapped to another reference sequence (for example, chromosome 4). In this case, the pairing information described above must be combined with additional information related to the reference sequence used to map one of the reads. This is achieved by coding

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

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

3) the third element contains information about mapping to the reference specified in paragraph 2, and is expressed as an offset relative to the last encoded position.

An example of this scenario is shown in Fig. 7.

Since in FIG. 7 read 4 is not mapped to the current encoded reference sequence, the genomic encoder signals this information by creating additional descriptors in the pair block. In the example shown below, read 4 of pair 2 is mapped to reference no. 4, while reference no. 1 is currently being mapped. This information is encoded using 3 components:

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

2) the second descriptor provides the reference ID, which is listed in the main header (in this case 4);

3) the third element contains information about mapping to the corresponding reference (170).

Mismatch descriptors for class N reads.

Class N includes all reads in which there are only n-type mismatches - instead of a base A, C, G or T, N was found as the specific base. All other bases of the reads perfectly match the reference sequence.

In fig. 8 shows how the position of N in read 1 is encoded as an absolute position in read 1, or as a differential position relative to the previous N in the same read;

the position of N in read 2 is encoded as the absolute position in read 2 + the length of read 1, or the difference position relative to the previous N.

In the nmis block, the encoding of each pair of reads is terminated with a special delimiter character.

Descriptors coding for substitutions (mismatches and SNPs), insertions and deletions.

A substitution is defined as the presence of a different nucleotide base in the mapped read compared to 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 substitution is coded as position (snpp block) and type (snpt block). Depending on the statistical occurrence of substitutions, insertions or deletions, different models of the source of the corresponding descriptors can be determined and the characters encoded in the corresponding block can be generated.

Source model 1: replacements in the form of type provisions.

Substitution position descriptors.

The replacement position is calculated similarly to the values of the nmis block, namely the following.

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

In read 2, replacements are encoded as the absolute position in read 2 + the length of read 1, or as a difference 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 coding);

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

In the snpp block, the encoding of each pair of reads is terminated with a special delimiter character.

- 18 043338

Substitution type descriptors.

For class M (and for class I, as described in the following sections), mismatches are encoded by index (moving from right to left) from the actual character present in the reference sequence to the corresponding replacement character present in the read (A, C, G, T, N, Z). For example, if an aligned read contains a C instead of a T that is present at the same position in the reference sequence, the mismatch index will be designated as 4. The decoding process reads the encoded descriptor, the nucleotide at a given position in the reference, and moves from left to right to extract the decoded character. For example, 2 received for a position where G is present in the reference will be decoded as N. In FIG. 11 shows all possible situations and the corresponding encoding symbols. To minimize the entropy of descriptors, each index can be assigned distinctly distinct and context-adaptive probabilistic models according to the statistical properties of each replacement type for each data class.

If IUPAC ambiguity codes are accepted, the results of the replacement mechanism will be the same, however 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. Figure 12 shows an example of encoding types of replacements in the snpt block. Some examples of encoding substitutions when implementing IUPAC ambiguity codes are shown in FIG. 13. Another example of substitution indices is shown in FIG. 14.

Coding of insertions and deletions.

For class I, mismatches and deletions are encoded by indices (moving from right to left) from the actual character present in the reference sequence to the corresponding substitution character present in the read {A, C, G, T, N, Z}. For example, if an aligned read contains a C instead of a T at the same position in the reference, the mismatch index will be 4. In the case where the read contains a deletion at a position where an A is present in the reference, the encoded symbol will be 5. The decoding process reads the encoded descriptor, nucleotide at a given position on the reference and moves from left to right to extract the decoded symbol. For example, 3 taken for a position where G is present in the reference will be decoded as Z.

For inserted A, C, G, T, N insertions are coded as 6, 7, 8, 9, 10, respectively. In fig. Figure 15 shows an example of how substitutions, insertions, and deletions are encoded in class I read pairs. To support the full set of IUPAC ambiguity codes, the substitution vector

S={A, C, G, T, N, Z} should be replaced with

S={A, 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, if the replacement vector has 16 elements, the insertion codes should have different values, namely 16, 17, 18, 19, 20. This mechanism is illustrated in FIG. 16.

Source model 2: one block per replacement type and indels.

For some data statistics, a different model for encoding substitutions and indels than described in the previous section may be developed, resulting in a source with lower entropy. This coding model is an alternative to the methods described above for mismatches only and for mismatches and indels. In this case, one data block is defined for each possible substitution character (5 without IUPAC codes, 16 with IUPAC codes) plus one block for deletions and another 4 blocks for insertions. For ease of explanation, but 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 block contains the position of mismatches or insertions of the same type. If there are no mismatches or insertions of this type in the read pair being encoded, a 0 is encoded in the corresponding block. To enable the decoder to begin the decoding process for the blocks described in this section, the header of each access block contains a flag indicating by signaling the first block to be encoded. In the example in FIG. 18, the first element to be decoded is at position 2 in block C. When there are no mismatches or indels of this type in a read pair, 0 is added to the corresponding blocks. When on the decode side, the decode pointer for each block points to the value 0, the decoding process proceeds to next pair of reads.

Coding of additional alarm flags.

For each data class introduced above (P, M, N, I), additional information about the nature of the reads being encoded may need to be encoded. This information may relate, for example, to a sequencing experiment (eg, indicating the probability of a duplicate read) or may express some characteristics of the read mapping (eg, first or second in a pair). In the context of the present invention, this information is encoded in a separate block 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 region of the reference sequence. Other examples of using such flags:

- 19 043338 paired reads;

read mapped to the correct pair;

unmapped read or paired read;

read or paired read from the reverse chain;

first/second read in a pair;

not the main read;

The read does not pass the platform/vendor quality check;

the read is a PCR duplicate or an optical duplicate;

additional alignment.

Descriptors for class U and construction of internal references for unmapped reads of class U and class HM

In the case of reads belonging to class U or an unmapped pair of class HM, since they cannot be mapped to any external reference sequence that satisfies a given set of match accuracy constraints for belonging to class P, N, M or I, for a compressed representation of reads belonging to these data classes, construct and use one or more internal reference sequences.

Several approaches are possible to constructing suitable internal references, such as, for example, but not limited to the following.

Dividing unmapped reads into clusters containing reads that share a common contiguous genomic sequence of at least a minimum size (signature). Each cluster can be uniquely identified by its signature, as shown in FIG. 22.

Sorting the reads in any meaningful order (for example, lexicographic order) and using the last N reads as the internal reference for (Ν+1)-γο encoding. This method is shown in Fig. 23.

Performing a so-called de-assembly on a subset of reads of class U to be able to align and encode all or a corresponding subset of the reads belonging to that class according to specified match accuracy constraints or a new set of constraints.

If the read being encoded can be mapped to an internal reference that satisfies a given set of match accuracy constraints, the information needed to reconstruct the read after compression is encoded using descriptors, which can be one of the following types.

1. The initial position of the matching part on the internal reference in the form of a read number in the internal reference (pos block). This position can be encoded as an absolute or differential value relative to a previously encoded read.

2. Offset of the starting position from the beginning of the corresponding read in the internal reference (pair block). For example, in the case of a constant read length, the actual position would be pos *length+pair.

3. Mismatches may be present, coded as position (snpp block) and type (snpt block).

4. Those parts of the reads (usually edges indicated by pair) that do not match the internal reference (or match, but with a number of mismatches above a given threshold) are encoded in the indc block. To reduce the entropy of the mismatches encoded in the indc block, a padding operation can be applied to the edges of the portion of the internal reference being used, as shown in FIG. 24. The most appropriate padding strategy can be selected by the encoder according to the statistical properties of the genomic data being processed. Possible filling strategies are:

a) without filling;

b) a constant padding pattern selected according to its frequency in the current data being encoded;

c) a variable padding pattern according to the statistical properties of the current context determined by the last N coded reads.

The specific type of padding strategy will be signaled by special values in the indc block header.

5. A flag that indicates whether the read is encoded using an internal self-generated reference, an external reference, or without a reference (rtype block).

6. Reads, which encode literally (reads). An example of such an encoding procedure is shown in Fig. 24.

In fig. Figure 25 shows an alternative encoding of mapped reads to an internal reference, where the pos+pair descriptors are replaced by a signed pos descriptor. In this case, pos expresses the distance, in units of positions, on the reference sequence - the position of the leftmost nucleotide of read n relative to the position of the leftmost nucleotide of read n-1.

In the case of variable-length reads of class U, an additional rlen descriptor is used to store the length of each read.

This coding approach can be extended to support N starting positions for each

- 20 043338 reads, so that reads can be split into two or more reference positions. This can be particularly useful for encoding reads generated by sequencing technologies (such as those from Pacific Bioscience) that produce very long reads (over 50K bases) that typically contain repetitive structures generated by cycling in the sequencing technology. The same approach can also be used to encode chimeric sequence reads, defined as reads that are aligned to two different positions in the genome with little or no overlap.

It is clear that the approach described above can be applied beyond the simple class U and could be applied to any block containing descriptors related to read positions (pos blocks).

Alignment evaluation descriptor.

The mscore descriptor provides a score for each alignment. In the context of the present invention, it is used to represent the mapping/alignment score for each read generated by the genomic sequence read alignment tools.

The evaluation is expressed using a power and a mantissa. The number of bits used to represent the exponent and mantissa are passed as configuration parameters. As an example, but not a limitation, the table. Figure 2a shows how it is defined in IEEE RFC 754 for an 11-bit exponent and a 52-bit mantissa.

The score of each alignment can be represented using a single sign bit (S);

bits for degree (E);

bits for mantissa (M).

Table 2a

Alignment estimates can be expressed as 64-bit double precision floating point values

The radix used to calculate scores is 10, so score = -1 ^s x 10 ^e x M

Read groups.

During the sequencing process, different types of sequenced reads can be obtained. By way of example, and not limitation, types may apply to different sequenced samples, different experiments, or different sequencer configurations.

Once sequenced and aligned, this information is stored in accordance with the present invention using a dedicated descriptor called rgroup. rgroup is a label associated with each encoded read that allows the decoder to separate the decoded reads into groups after decoding.

Descriptors for multiple alignments

The following descriptors are defined to support multiple alignments. For the case of spliced reads, the present invention defines a global flag spliced_reads_flag, which should be set to 1.

mmap handle.

The mmap descriptor is used to signal how many read positions or the leftmost read of a pair have been aligned. A genomic record containing multiple alignments is associated with a single multi-byte mmap descriptor.

The first two bytes of the mmap descriptor represent an unsigned integer N that refers to the read as a single segment (if there are no splices in the encoded dataset) or, alternatively, to all segments into which the read was spliced for multiple possible alignments (if there are no splices in the encoded dataset) data contains splices).

The value N indicates how many pos descriptor values are encoded for the pattern in a given entry. N is followed by one or more unsigned integers Mi, as described below.

Chaining of multiple alignments.

The rcomp descriptor described in the present invention is used to indicate the chainness of each read alignment using the syntax defined in the present invention.

Multiple alignment estimates.

In the case of multiple alignments, each alignment is assigned one mscore descriptor as defined in the present invention.

Multiple alignments without splicing.

If there are no splices in the access block, spliced_reads_flag is not set.

For paired-end sequencing, the mmap descriptor consists of a 16-bit integer without

- 21 043338 sign N followed by one or more 8-bit unsigned integers Mi, with i taking values from 1 to the number of complete alignments of the first (in this case leftmost) read. For each alignment of the first read, spliced or not, Mi is used to signal how many segments are used to align the second read (in this case without splicing it is equal to the number of alignments) and, furthermore, how many pair descriptor values are encoded for that alignment of the first read .

Σ™

The Mi values are used to calculate P - ' ¹ , which indicates the number of alignments of the second read.

The special value Mi (Mi=0) indicates that the i-e alignment of the leftmost read is paired with the alignment of the rightmost read, which is already paired with the k-th alignment of the leftmost read, with k<i (in which case no new alignment is found, which is consistent with the above equation). As an example, in the simplest cases, if there is one alignment for the leftmost read and two alternative alignments for the rightmost read, N will be equal to 1 and M1 will be equal to 2;

if two alternative alignments are found for the leftmost read, but only one for the rightmost read, N will be equal to 2, M1 will be equal to 1, and _M2 will be equal to 0;

when Mi is 0, the corresponding pair value should be bound to the existing alignment of the second read;

otherwise, a syntax error will occur and the alignment will be considered broken.

Example. If the first read has two mapping positions and the second read only one, N is 2, M ₁ is 1, and M ₂ is 0, as above. If this is followed by another alternative auxiliary mapping for the entire pattern, N will be equal to 3 and M ₃ will be equal to 1.

Fig. 39 illustrates the meaning of N, P and Mi in the case of multiple alignments without splices, and in Error! Reference source not found. In fig. Figure 10 shows how the pos, pair, and mmap descriptors are used to encode multiple alignment information.

Referring to FIG. 40 the following is true:

The rightmost read has P = ¹ - ¹ alignments.

Some values of Mi can be =0 when the i-e alignment of the leftmost read is paired with the alignment of the rightmost read, which is already paired with the k-th alignment of the leftmost read, and k<i.

One reserved pair descriptor value may be present to signal alignments that belong to other access block ranges. If present, it is always the first pair descriptor for the current entry.

Multiple alignments with splices.

When encoding a data set with spliced reads, the msar descriptor allows the length and stringiness of the splices to be represented.

After decoding the mmap and msar descriptors, the decoder knows how many reads or pairs of reads were encoded to represent multiple mappings, and how many segments constitute each mapping of each read or pair of reads. This is shown in Fig. 41 and 42.

Referring to FIG. 41 the following is true.

The leftmost read has N1 alignments with N splices (N1<N).

N represents the number of splices present in all leftmost read alignments and is encoded as the first value of the mmap descriptor.

,νι

Σ ^Μί

The rightmost read has P = ¹¹ splices, where Mi is the number of splices of the rightmost read that is paired with the ith alignment of the leftmost read (1<i<N1). In other words, P represents the number of splices of the rightmost read and is calculated using the N values following the first value of the mmap descriptor.

N1 and N2 represent the number of alignments of the first and second read and are calculated using the N + P values of the msar descriptor.

Referring to FIG. 42 the following is true.

The leftmost read has N1 alignments with N splices (N1<N). If N ₁ =N and N ₂ =P, there will be no splices.

JV1

Σ ^Μί

The rightmost read has P - ¹ - ¹ splices (tj 1<]<P) = and N2 (N ₂ <P) alignments.

The number of pair descriptors can be calculated as

- 22 043338

NP=Max(N ₁ , P)+M ₀ , where M ₀ is the number of Mi with value 0;

NP must be incremented by 1, in which case one special pair descriptor indicates the presence of alignments in other databases.

Alignment assessment.

The mscore descriptor allows the alignment score to be indicated by signaling. With single-end sequencing it will have N1 values per template; with paired-end sequencing it will have a value for each alignment of the entire template (the number of different alignments of the first read and possibly plus the number of additional alignments of the second read, i.e. when Mi-1>0).

Number of evaluations= MAX(N1, N ₂ )+M _o , where M ₀ represents the total number Mi=0.

In the present invention, more than one score value may be associated with each alignment. The number of alignments is signaled by the as_depth configuration parameter.

Descriptors for multiple alignments without splicing.

Table 3

Determining the number of descriptors required to represent multiple alignments in a single genomic record in the case of multiple alignments without splices

Semantics mmap mscore Result Reed (or paired-terminal read) with several mappings, Single read: N Pair of reads: N, Mi where Single read: N values Pair of reads: ΜΑΧ(Ν, Σ(Μ,)) values + Σ(Μ, = 0) Introducing a separator allows you to have any Single read: a read has multiple mappings and its encoded as sequences of N Not spliced 1 <i<N number of ratings. IN otherwise, if it is not present, you should use 0. These are floating point values as defined in the present invention. consecutive segments belonging to the class with the highest ID. N pos descriptors are used. Read Pair: A read pair has multiple mappings and is encoded as a sequence of: N segments for the first read; Ρ =Σ(Μι) pairings for second-read alignments. N pos descriptors are used. N x P pair descriptors are used, and N Mind/ NP = ί = ι + (number Mi = 0) + (optional)1 An optional pairing descriptor is used, when leveling are present on other reference sequences than the one currently being encoded. N+1 mmap descriptors

- 23 043338

Reed (couple), uniquely mapped 0 0

Descriptors for multiple splice alignments.

In table Figure 4 shows how to determine the number of descriptors needed to represent multiple alignments in a single genomic record in the case of multiple alignments without splices.

Table 4

Descriptors used to represent multiple alignments and associated scores

Semantics mmap mscore Result Reed (or paired-terminal read) with several mappings, with splices Single read: N Pair of reads: Ν, Μι where 1 <i < N Single read: Ni values Pair of reads: Max(Ni, N ₂ ) + Σ(Μι== 0) values Ni and Ν ₂ using Ν + P descriptors msar W1 _Ρ= Σ- These are floating point values as defined herein invention. Single read: a read has multiple mappings and its encoded as sequences of N consecutive segments belonging to the class with the highest identifier. N pos descriptors are used. Read Pair: A read pair has multiple mappings and is encoded as a sequence of: N segments for the first read; P =Σ(Μι) pairings for second-read alignments. N pos descriptors are used. Reed (couple), uniquely mapped 0 0

Multiple alignments on different sequences.

It may happen that the alignment process reveals alternative mappings to a different reference sequence than the one on which the main mapping is located.

For pairs of reads that are uniquely aligned, the pair descriptor should be used to represent the absolute positions of the read when there is, for example, a chimeric alignment with a paired read on another chromosome. The pair descriptor should be used to indicate, by signaling, the reference and position of the next entry containing additional alignments for the same pattern. The last record (for example, the third, if alternative mappings are encoded in 3 different databases) must contain the reference and position of the first record. In the case where one or more alignments for the leftmost read in a pair are not on the reference sequence that belongs to the current database being encoded, I use a reserved value for the pair descriptor. The reserved value is followed by the reference sequence identifier and the leftmost alignment position among all those contained in the next database (i.e., the first decoded pos descriptor value for that entry).

Multiple alignments with insertions, deletions, unmapped parts.

When alternative auxiliary mapping does not preserve the continuity of the reference region where the sequence is being aligned, reconstruction of the accurate mapping generated by the aligner may not be possible because the actual sequence (and thus descriptors associated with mismatches such as substitutions or indels) is encoded only

- 24 043338 for basic alignment. To represent how auxiliary alignments are mapped to a reference sequence, the msar descriptor should be used when they contain indels and/or soft cuts. If msar is represented by the special character * for auxiliary alignment, the decoder will reconstruct the auxiliary alignment from the main alignment and the auxiliary alignment mapping position.

msar handle.

The msar (Multiple Segments Alignment Record) descriptor supports spliced reads and alternative auxiliary alignments that contain indels or soft cuts.

msar is intended to convey information about the length of the mapped segment;

other mapping continuity (i.e. presence of insertions, deletions or cut bases) for assisted alignment and/or spliced read.

msar uses the CIGAR extended notation syntax described below, plus the additional character described in Table. 5.

Table 5

A special character used for the msar descriptor in addition to the syntax described in Table. 6

Symbol Semantics Description * Auxiliary alignment is not contains indels or soft cuts This is used when reconstructing the auxiliary alignment does not require any additional information other than the alignment position and the primary alignment

Extended cigar syntax.

This section defines the extended CIGAR (E-CIGAR) syntax for strings that must be associated with sequences and corresponding mismatches, indels, trim bases, and information about multiple alignments and spliced reads.

The editing operations described in the present invention are listed in table. 6.

Table 6

E-CIGAR string syntax in MPEG-G format

Operation Semantics E-CIGAR introduction CIGAR representation equivalent to SAM Increment of the pointer to reference R and the pointer to read g by n positions (match) η matching reasons n= PM in more old versions (not equivalent), = in recent versions Replacement of a nucleotide in a read with base C from reference, increment Replacement of the symbol C (C is present in the read and absent in reference) WITH M in older versions,

- 25 043338

pointer to reference R and pointer to read g by 1 X in recent versions (not equivalent) Increment of the pointer on the read by n positions (insertion from the read) η bases inserted in read (not present in reference) n+ ηΙ Increment of the pointer to the reference I by l positions (deleting sequences S in the read) η bases are removed from the read (but present in reference) P- nD Increment of the pointer to the reference I by l positions (insertion into the read). Can only take place at the beginning or end of a read η soft cutoffs (P) nS Hard cutting. Can only take place at the beginning or end of the read η hard cutoffs [P] nH Increment of the pointer to the reference R by l positions in compliance with the splice consensus (splice in read) Non-directional splice from η reasons P* nN Increment of the pointer to the reference R by l positions, compliance with the splice consensus on the direct chain (direct splice in the read) Direct splicing of η bases P/ Does not exist Increment of the pointer to the reference R by l positions, compliance with the splice consensus on the reverse chain (reverse splice in read) Reverse splicing from η bases P% Does not exist

Source models, entropy encoders and encoding modes.

For each data class, subclass, and corresponding genomic data structure descriptor block in the present invention, different encoding algorithms may be implemented in accordance with the specific features of the data or metadata contained in each block and their statistical properties. The encoding algorithm should be understood as the relationship of a specific descriptor block source model with a specific entropy encoder. A specific source model can be defined and selected to obtain the most efficient data encoding in terms of minimizing source entropy. The selection of an entropy encoder may be guided by considerations of coding efficiency and/or probability distribution characteristics and associated implementation issues. Each choice of a specific encoding algorithm, also called an encoding mode, can be applied to the entire block of descriptors associated with a class or subclass of data for the entire data set, or a different encoding mode can be applied to each portion of the descriptors divided into access blocks. Each source model associated with an encoding mode is characterized by defining the descriptors produced by each source (i.e., the set of descriptors used to represent a class of data, such as read positions, read pairing information,

- 26 043338 mismatch with the reference sequence, as defined in table. 2);

defining a coupled probabilistic model;

Definition of the associated entropy encoder.

Further benefits.

Classifying sequence data into specific classes and subclasses of data allows for the implementation of efficient encoding modes that take advantage of the lower entropy of the information source, characterized by modeling descriptor sequences as single individual data sources (e.g., distance, position, etc.).

Another advantage of the invention is the ability to access only a subset of data of the type of interest. For example, one of the most important applications in genomics is the search for differences in a genomic sample compared to a reference (NRV) or population (RNP). Today, this type of analysis requires the processing of full sequence reads, whereas when implementing the data representation disclosed in the present invention, mismatches are already isolated into just three classes of data (n-type and i-type mismatches are also considered if of interest).

Another advantage is the ability to perform efficient transcoding of data and metadata compressed against a particular external reference sequence into compressed relative to another external reference sequence when publishing new reference sequences or when performing remapping on already mapped data (for example, using a different mapping algorithm). with obtaining new alignments.

In fig. 20 shows an encoder 207 in accordance with the principles of the present invention. Encoder 207 receives as input raw sequence data 209, such as generated by genome sequencing device 200. Devices for genome sequencing are known in the field, such as Illumina HiSeq 2500, Thermo-Fisher Ion Torrent devices. The raw sequence data 209 is provided to an aligner node 201, which prepares sequences for coding by aligning the reads to a reference sequence 2020. In an alternative embodiment, a dedicated module 202 may be used to generate a reference sequence from available reads using various strategies as described herein document in the section Construction of internal references for unmapped reads of class U and class HM. After processing by the reference generator 202, the reads can be mapped to the resulting longer sequence. The aligned sequences are then classified using a classification module 204. Thereafter, a reference transformation step is applied to the reference to reduce the entropy of the data generated by the data classification node 204. It involves processing the external reference 2020 in a reference transformation node 2019, which creates transformed data classes 2018 and reference transformation descriptors 2021. The transformed data classes 2018 are then supplied to block encoders 205-207 along with reference transform descriptors 2021. The genomic blocks 2011 are then provided to arithmetic encoders 2012-2014, which encode the blocks according to statistical properties of the data or metadata contained in the block. The result is Genome Stream 2015.

In fig. 21 shows a decoder 218 in accordance with the principles of the present invention. The decoder 218 receives the multiplexed genomic binary data stream 2110 from the network or storage element. The multiplexed genomic binary data stream 2110 is fed to a demultiplexer 210 to create individual streams 211, which are then fed to entropy decoders 212-214 to produce genomic blocks 215 and reference transform descriptors 2112. The extracted genomic blocks are supplied to block decoders 216-217 for further decoding into data classes, and the reference transformation descriptors are supplied to the reference transformation node 2113. Class decoders 219 further process the genomic descriptors 2111 and the transformed reference 2114 and combine the results to produce uncompressed sequence reads, which are then also stored in formats known in the art, such as a text file or a compressed zip file, or FASTQ files. or SAM/BAM.

Class decoders 219 can recover the original genomic sequences by using the original reference sequence information contained in one or more genomic streams and the reference transform descriptors 2112 in the encoded binary data stream. If reference sequences are not transported via genomic streams, they must be present on the decoding side and accessible to class decoders. The inventive methods described herein may be implemented in hardware, software, firmware, or any combination thereof. When implemented in software, they can be stored on a computer-readable storage medium and executed by a hardware processing device. A hardware processing device may include one or more processors, digital signal processors, general purpose microprocessors, application specific integrated circuits, or other discrete logic circuits.

- 27 043338

The methods of this description can be implemented in a wide variety of devices and instruments, including mobile phones, desktop computers, servers, tablets and similar devices.

File format: Selective access to genomic data regions using a master index table.

To support selective access to specific areas of aligned data, the data structure described in this document implements an indexing facility called a master index table (MIT). This is a multidimensional array containing the mapping locations of specific reads to associated reference sequences. The values contained in MIT represent the mapping positions of the first read in each pos block in order to support non-sequential access to each access block. MIT contains one section for each data class (P, N, M, I, U, and HM) and each reference sequence. MIT is contained in the header of the genomic dataset of encoded data. In fig. 31 shows the header structure of a genomic dataset, FIG. 32 shows an overall visual representation of MIT, and FIG. Figure 33 shows an MIT example for class P encoded reads.

Depicted in FIGS. The 33 values contained in MIT are used to directly access the region of interest (and the corresponding database) in the compressed region. For example, as shown in FIG. 33, if it is necessary to access the region contained between positions 150000 and 250000 in reference 2, the decoding application must go directly to the second reference in MIT and find two values of k1 and k2 such that k1<150000 and k2>250000. Where k1 and k2 are two indices read from MIT. In the example in FIG. 33 will result in the 3rd and 4th positions of the second MIT vector. These returned values will then be used by the decoding application to retrieve the corresponding data positions from the pos block's local index table, as described in the next section. Together with pointers to a block containing data belonging to the four classes of genomic data described above, MIT can be used as an index of additional metadata and/or annotations added to genomic data during its lifetime.

Local pointer table.

At the beginning of each genomic data block there is a data structure called a local header. The local header contains a unique block identifier, a vector of access block counters for each reference sequence, a local index table (LIT), and optionally some block-specific metadata. The LIT is a vector of pointers to the physical location of the data belonging to each access block in the block's payload. In fig. 34 depicts a common block header and payload where LIT is used to access specific areas of encoded data in a non-sequential manner. In the previous example, to access a region between 150,000 and 250,000 reads aligned to reference sequence #2, the decoding application retrieved positions 3 and 4 from MIT. These values must be used by the decoding process to access the 3rd and 4th elements of the corresponding LIT section. In the example shown in FIG. 35, the access block total counters contained in the block header are used to jump directly to the LIT indices associated with the DB related to reference 1 (5 in the example). Therefore, indices containing the physical positions of the requested DBs in the codestream are calculated as follows:

position of data blocks belonging to the requested DBs = data blocks belonging to reference DB 1 to be skipped + position retrieved using MIT, i.e.

first block position: 5+3=8;

position of the last block: 5+4=9.

Data blocks retrieved by an indexing mechanism called a local index table are part of the requested access blocks.

In fig. Figure 36 shows how the blocks contained in the MIT correspond to the LIT blocks for each data class or subclass.

In fig. 37 shows how data blocks retrieved by MIT and LIT form one or more access blocks as defined in the next section.

In an embodiment of the present invention, the LIT may be embedded as a substructure of the MIT. The advantage of this approach is the speed of access to indexed data in the case of sequential parsing of a compressed file. If LIT is embedded in MIT in the file header, then in the case of selective access, the decoder will only need to parse a small portion of the data to extract the requested compressed information. Another advantage that will be apparent to one skilled in the art is that, if streamed over a network, the indexed information contained in the MIT and LIT will be delivered among the first blocks of data, thereby allowing the receiving device to perform operations such as sorting and selective access until all data has been transferred.

Access block.

Genomic data, classified into data classes and structured into compressed and uncompressed blocks, is organized into different access blocks.

- 28 043338

Genomic access blocks (DBs) are defined as sections of genomic data (compressed and uncompressed) that recover nucleotide sequences, and/or corresponding metadata, and/or DNA/RNA sequences (for example, virtual reference), and/or annotation data, generated by a genome sequencer and/or genomic processing device or analysis application. An example of an access block is shown in Fig. 37.

An access block is a block of data that can be decoded either independently of other access blocks using only globally available data (eg, decoder configuration) or using information contained in other access blocks. Access blocks are distinguished by type, which characterizes the nature of the genomic data and data sets that contain these access blocks, and the possible way to access them;

an order that provides a unique order to access blocks belonging to the same type.

Access units of any type can be further classified into different categories.

The following is a non-exhaustive list of definitions of various types of genomic access blocks.

1) Accessing or decoding and accessing Type 0 access blocks does not require reference to any information coming from other access blocks. All information in the data or data sets they contain can be independently read and processed by a decoder or processing application.

2) Type 1 access blocks contain data that references data contained in type 0 access blocks. To read or decode and process data contained in type 1 access blocks, access to one or more type 0 access blocks is required. Access block type 1 encodes genomic data related to class P sequence reads.

3) Type 2 access blocks contain data that references data contained in type 0 access blocks. Reading or decoding and processing data contained in type 2 access blocks requires access to one or more type 0 access blocks. Access block type 2 encodes genomic data related to class N sequence reads.

4) Type 3 access blocks contain data that references data contained in type 0 access blocks. Reading or decoding and processing data contained in type 3 access blocks requires access to one or more type 0 access blocks. Access blocks type 3 encode genomic data related to class M sequence reads.

5) Type 4 access blocks contain data that references data contained in type 0 access blocks. Reading or decoding and processing data contained in type 4 access blocks requires access to one or more type 0 access blocks. Access block type 4 encodes genomic data related to class I sequence reads.

6) Type 5 access blocks contain reads that cannot be mapped to any of the available reference sequences (class U) and are encoded using an internally constructed reference sequence. Type 5 access blocks contain data that references data contained in type 0 access blocks. Reading or decoding and processing data contained in type 5 access blocks requires access to one or more type 0 access blocks.

7) Access blocks of type 6 contain pairs of reads, where one read can belong to any of the four classes P, N, M, I, and the other cannot be mapped to any of the available reference sequences (class HM). Type 6 access blocks contain data that references data contained in type 0 access blocks. Reading or decoding and processing data contained in type 6 access blocks requires access to one or more type 0 access blocks.

8) Type 7 access blocks contain metadata (eg, quality ratings) and/or annotation data associated with the data or data sets contained in a type 1 access block. Type 7 access blocks can be classified and labeled in different blocks.

9) Type 8 access blocks contain data or data sets classified as annotation data. Type 8 access blocks can be classified and labeled in blocks.

10) Access blocks of additional types can extend the structures and mechanisms described here. By way of example, and not limitation, the results of genomic variant determination, structural and functional analysis may be encoded in new types of access units. The organization of data into access blocks described herein does not prevent any type of data from being encapsulated into access blocks, which is a mechanism that is completely transparent about the nature of the encoded data.

Type 0 access blocks are ordered (eg, numbered), but do not need to be stored and/or transmitted in an orderly manner (technical advantage: parallel processing/parallel streaming, multiplexing).

Access blocks of type 1, 2, 3, 4, 5 and 6 do not need to be ordered and do not need to be stored and/or transmitted in an orderly manner (technical advantage: parallel processing/parallel streaming

-

Claims (19)

broadcast). In fig. 37 shows how access blocks are composed using a header and one or more homogeneous data blocks. Each block can be composed of one or more blocks. Each block contains several packages, and the packages are a structured sequence of descriptors introduced above to represent, for example, positions of reads, pairing information, reverse complement information, positions and types of mismatches, etc. Each access block may have a different number of packets in each block, but within an access block, all blocks have the same number of packets. Each data packet can be identified by a combination of 3 identifiers X, Y and Z, where X indicates the access block to which it belongs; Y indicates the block it belongs to (i.e. the type of data it encapsulates); Z is an identifier expressing the order of a packet relative to other packets in the same block. In fig. 38 shows an example of marking access units and packets, where AU_T_N is an access unit of type T with an identifier N, which may or may not imply ordering according to the access unit type. Identifiers are used to uniquely associate access units of one type with access units of other types required to fully decode the contained genomic data. Access blocks of any type can be further classified and labeled into different categories according to different sequencing processes. By way of example, and not limitation, classification and labeling may occur when 1. A computer-implemented method for encoding genomic sequence data, where the genomic sequence data contains reads of nucleotide sequences, the method including the steps of aligning said reads to one or more first reference sequences, thereby creating aligned reads, classifying aligned reads according to specified matching rules with the specified one or more first reference sequences into different classes, thereby creating classes of aligned reads, wherein the classification includes identifying as the first class (class P) genomic reads without any mismatches with the reference sequence, where in the mapped read there are no mismatches with the reference sequence, used for mapping; identifying as a second class (class N) genomic reads where mismatches are found only at positions where the sequencer was unable to detect any base and the number of mismatches in each read does not exceed a given threshold; identification as a third class (class M) of genomic reads where mismatches are found at positions in which the sequencer failed to identify any n-type base and/or it identified a different base compared to the s-type reference sequence and the number of mismatches is not exceeds the given thresholds for the number of n-type, s-type mismatches and the threshold obtained from a given function (f(n,s)); identification as a fourth class (class I) of genomic reads that may have mismatches of the same type as mismatches of the third class (class M), and additionally at least one mismatch of the following type: insertion (i-type), deletion (d) -type), soft or hard cutoffs (c-type), and the number of mismatches for each type does not exceed the corresponding specified threshold and the threshold provided by the specified function (w(n,s,i,d,c)), coding classified aligned reads as a plurality of descriptor blocks, wherein encoding said classified aligned reads into a plurality of descriptor blocks includes selecting descriptors according to classes of the aligned reads, structuring said descriptor blocks using header information thereby creating successive access blocks, wherein the access blocks of the first class (class P) are built using descriptor blocks for information about the mapping position, descriptor blocks for information about which strand of DNA the read was sequenced from (chainness), and flags for information about certain characteristics of the sequence read, while in the specified access blocks class P mating information - 30 043338 paired-end reads are encoded using the appropriate descriptor block, with second-class (N-class) access blocks constructed using the same descriptor blocks as for first-class (P-class) access blocks, plus a descriptor block for information about the position of unknown bases , while access blocks of the third class (class M) are built using the same descriptor blocks as for access blocks of the first class (class P), plus descriptor blocks for information about the position and types of replacements, while access blocks of the fourth class (class I ) are constructed using the same descriptor blocks as for first class (P class) access blocks, plus descriptor blocks for information about the position and type of substitutions, insertions, deletions and base truncations, where the method further includes identifying genomic reads that are not assigned to none of the first to fourth classes (classes P, N, M, I) as a fifth class (class U), constructing a set of second reference sequences using at least some of the fifth class reads, aligning said fifth class reads on a set of second reference sequences, encoding the fifth class reads into corresponding descriptors based on certain constraints on the accuracy of the match with respect to the second reference sequences, structuring said corresponding descriptors using header information thereby creating access blocks of the fifth block. 1) sequencing the same organism at different times (access blocks contain genomic information on the genome with a temporal background); 2. The method according to claim 1, wherein the fifth class access blocks are constructed using one or more of the descriptor blocks for information about mapping positions; descriptor blocks for information about which DNA strand the read was sequenced from (chainness), and flags for information about certain characteristics of the sequence read, wherein information about the pairing of paired-end reads is encoded using the corresponding descriptor block; descriptor blocks for information about the position and type of substitutions; descriptor blocks for information about parts of reads that do not match the second reference sequences; descriptor blocks to literally encode reads that cannot be mapped to any reference sequence. 2) sequencing organic samples of various natures of some organisms (for example, skin, blood, hair for human samples). These access blocks are biologically inspired. CLAIM 3. The encoding method according to claim 1 or 2, in which the genomic sequence reads to be encoded are paired; and the classification further includes identifying the genomic reads as a sixth class (HM class) containing all pairs of reads where one read belongs to class P, N, M or I and the other read belongs to class U. 4. The encoding method according to claim 3, further including the following steps: determining whether two paired reads are classified into the same class (each of classes P, N, M, I, U) and then assigning the pair to the same identified class; determining whether two paired reads belong to different classes, and if neither of them belongs to class U, then assigning the pair of reads to the class with the highest priority, determined according to the following priority: p<n<m<i, where class P has the lowest priority, and class I has the highest priority; determining whether only one of the paired reads is classified as belonging to class U, and classifying the pair of reads as belonging to class HM sequences. 5. The method according to claim 4, wherein each read class N, M, I is further divided into two or more subclasses (296, 297, 298) in accordance with a threshold vector (292, 293, 294) defined respectively for each class N, M and I through the number of n-type mismatches (292), the function f(n,s) (293) and the function w(n,s,i,d,c) (294), further including the following steps: determining whether two paired reads are classified into the same subclass, then assigning the pair to the same subclass; determining whether two paired reads are classified into subclasses of different classes, then assigning the pair to the subclass belonging to the higher priority class according to the following expression: n<m<i, where N has the lowest priority and I has the highest priority; determining whether two paired reads are classified into the same class, which class is N, or M, or I, or into different subclasses, in which case assigning the pair to the highest priority subclass according to the following expressions: Ni<N ₂ <^<N _k , - 31 043338 M1<M2<,.<Mj, ¹ 1< ¹ 2< - < ¹ Н, where the highest index has the highest priority. 6. The method according to claim 5, wherein information about the mapping position of each read is encoded by a pos descriptor block; information about which DNA chain the read (chain) was sequenced from is encoded for each read using the rcomp descriptor block; and information about paired-end reads is encoded by a pair descriptor block. 7. The method of claim 6, wherein additional information about the alignment, such as whether the read is mapped to the correct pair, whether it has failed platform/vendor quality control, whether it is a PCR duplicate or an optical duplicate, whether it is ancillary alignment, encoded using a flags descriptor block. 8. The method according to claim 7, wherein information about unknown bases is encoded using a block of descriptors nmis. 9. The method of claim 8, wherein information about the position of the substitutions is encoded by a snpp descriptor block; information about the type of substitutions is encoded using a snpt descriptor block. 10. The method of claim 9, wherein information about the position of mismatches such as substitutions, insertions or deletions is encoded by an indp descriptor block; information about the type of mismatches, such as substitutions, insertions or deletions, is encoded by the indt descriptor block; The cut base information of the mapped read is encoded using the indc descriptor block. 11. The method according to claim 10, wherein information about unmapped reads is encoded using a ureads descriptor block; information about the type of the reference sequence used for encoding is encoded by an rtype descriptor block; information about multiple alignments of mapped reads is encoded using an mmap descriptor block; information about spliced alignments and multiple alignments of the same read is encoded by an msar descriptor block and an mmap descriptor block; information about read alignment scores is encoded using the mscore descriptor block; information about the groups to which the reads belong is encoded using the rgroup descriptor block. 12. The method of claim 11, wherein class HM access blocks are constructed using the same descriptor blocks as class I access blocks for mapped reads, and using ureads descriptor blocks for unmapped reads. 13. The method of claim 12, wherein class I access blocks are constructed using the same descriptor blocks as for class P access blocks, plus indp, indt and indc descriptor blocks for information about the position and type of substitutions, insertions, deletions and base trims, with information about multiple alignments being conveyed using mmap and msar descriptor blocks. 14. The method of claim 13, wherein information about the spliced alignments is conveyed using an extended cigar string containing the symbol = to indicate matching bases, the symbol + to indicate insertions, the symbol - to indicate deletions, the symbol / to indicate splicing on the straight strand, % symbol to indicate reverse strand splice, * symbol to indicate non-directional splice, text symbol from IUPAC codes for DNA to indicate substitution, (n) symbol to indicate n soft truncated bases where n is an integer, [n] symbol to indicate n hard-cut bases, where n is an integer. 15. The method of claim 14, wherein the descriptor blocks comprise a master index table containing one partition for each class and subclass of aligned reads, wherein the partition contains mapping positions to one or more reference sequences of the first read of each access block of each class or subclass of data ; said main index table and access block data are encoded together. 16. The method according to claim 15, in which the descriptor blocks also contain information about the type of reference used, selected - 32 043338 from a pre-existing or constructed one, and read segments that are not mapped to this reference sequence, wherein preferably the reference sequences are first converted to other reference sequences by applying substitutions, insertions, deletions and truncations, and then encoding the classified aligned reads as a set of descriptor blocks related to the transformed reference sequences. 17. The method of claim 16, wherein the same transformation is applied to the reference sequences used for all data classes; or different transformations are applied to the reference sequences used for each data class, the transformations of the reference sequences being encoded as descriptor blocks and structured using header information thereby creating successive access blocks. 18. The method of claim 14, wherein different transformations are applied to the reference sequences used for each data class, wherein the transformations of the reference sequences are encoded as descriptor blocks and structured using header information thereby creating successive access blocks, wherein the encoding of the classified aligned reads and associated transformations of reference sequences in the form of a plurality of descriptor blocks includes the step of associating a specific source model and a specific entropy encoder with each descriptor block, wherein the entropy encoder is a context-adaptive arithmetic encoder, a variable length encoder, or a Golomb encoder. 19. A computer-implemented method for decoding encoded genomic data, which includes the following steps: parsing access blocks containing encoded genomic data to extract a plurality of descriptor blocks by using header information; and decoding a plurality of descriptor blocks to select reads in accordance with specified matching rules that determine their classification relative to one or more reference sequences, wherein for an access block belonging to the first, second, third or fourth class, the descriptor blocks describe the match of said reads with the first reference sequences according to specified matching rules, where an access block belongs to the first class (class P) if it represents genomic reads without mismatches with the reference sequence used for mapping; where an access block belongs to the second class (class N) if it represents genomic reads for which mismatches are found only at positions where the sequencer failed to detect any base, and the number of mismatches in each read does not exceed a given threshold; where an access block belongs to the third class (class M) if it represents genomic reads for which mismatches are found at positions where the sequencer failed to identify any n-type base and/or it identified a different base compared to the s-type reference sequence type, and the number of mismatches does not exceed the specified thresholds for the number of mismatches of n-type, s-type and the threshold obtained from the given function (f(n,s)); where an access block belongs to the fourth class (class I) if it represents genomic reads that may have mismatches of the same type as mismatches of the third class (class M), and additionally at least one mismatch of the following type: insertion (i-type ), deletion (d-type), soft or hard cuts (c-type), and the number of mismatches for each type does not exceed the corresponding specified threshold and the threshold provided by the specified function (w(n,s,i,d,c )); where first class (P class) access blocks are constructed using descriptor blocks for information about the mapping position, descriptor blocks for information about which DNA strand the read was sequenced from, and flags for information about certain characteristics of the sequence read, wherein in the specified class P access blocks, information about the pairing of paired-end reads is encoded using the corresponding descriptor block; where the access blocks of the second class (class N) are constructed using the same blocks of descriptors as for the access blocks of the first class (class P), plus a block of descriptors for information about the positions of unknown bases; where third class (M class) access blocks are constructed using the same descriptor blocks as for first class (P class) access blocks, plus descriptor blocks for information about the position and types of replacements; where fourth class (Class I) access blocks are constructed using the same descriptor blocks as for first class (P Class) access blocks, plus descriptor blocks for information about -