METHOD AND SYSTEM FOR SELECTIVE ACCESS OF STORED OR TRANSMITTED BIOINFORMATICS DATA

Description

TECHNICAL FIELD

The present application provides new methods for the efficient storage, transmission and multiplexing of bioinformatics data, and in particular genomic sequencing data, in compressed form that enable efficient selective access and selective protection of the different data categories composing the genomic datasets.

BACKGROUND

An appropriate representation of genome sequencing data is fundamental to enable efficient processing, storage and transmission of genomic data to make possible and facilitate analysis applications such as genome variants calling and all analysis performed, with various purposes, by processing the sequencing data and metadata. Today, genome sequencing information is generated by High Throughput Sequencing (HTS) machines in the form of sequences of nucleotides (a. k. a. bases) represented by strings of letters from a defined vocabulary.

These sequencing machines do not read out an entire genomes or genes, but they produce short random fragments of nucleotide sequences known as sequence reads. A quality score is associated to each nucleotide in a sequence read. Such number represents the confidence level given by the machine to the read of a specific nucleotide at a specific location in the nucleotide sequence. This raw sequencing data generated by NGS machines are commonly stored in FASTQ files (see also FIG. 1).

The smallest vocabulary to represent sequences of nucleotides obtained by a sequencing process is composed by five symbols: {A, C, G, T, N} representing the four types of nucleotides present in DNA namely Adenine, Cytosine, Guanine, and Thymine plus the symbol N to indicate that the sequencing machine was not able to call any base with a sufficient level of confidence, so the type of base in such position remains undetermined in the reading process. In RNA Thymine is replaced by Uracil (U). The nucleotides sequences produced by sequencing machines are called “reads”. In case of paired reads the term “template” is used to designate the original sequence from which the read pair has been extracted. Sequence reads can be composed by a number of nucleotides in a range from a few dozen up to several thousand. Some technologies produce sequence reads in pairs where each read can be originated from one of the two DNA strands.

In the genome sequencing field the term “coverage” is used to express the level of redundancy of the sequence data with respect to a reference genome. For example, to reach a coverage of 30× on a human genome (3.2 billion bases long) a sequencing machine shall produce a total of about 30×3.2 billion bases so that in average each position in the reference is “covered” 30 times.

State of the Art Solutions

The most used genome information representations of sequencing data are based on FASTQ and SAM file formats which are commonly made available in zipped form in the attempt of reducing the original size. The traditional file formats, respectively FASTQ and SAM for non-aligned and aligned sequencing data, are constituted by plain text characters and are thus compressed by using general purpose approaches such as LZ (from Lempel and Ziv) schemes (the well-known zip, gzip etc). When general purpose compressors such as gzip are used, the result of the compression is usually a single blob of binary data. The information in such monolithic form results quite difficult to archive, transfer and elaborate particularly in the case of high throughput sequencing when the volumes of data are extremely large.

After sequencing, each stage of a genomic information processing pipeline produces data represented by a completely new data structure (file format) despite the fact that in reality only a small fraction of the generated data is new with respect to the previous stage.

FIG. 1 shows the main stages of a typical genomic information processing pipeline with the indication of the associated file format representation.

Commonly used solutions presents several drawbacks: data archival is inefficient for the fact that a different file format is used at each stage of the genomic information processing pipelines which implies the multiple replication of data, with the consequent rapid increase of the required storage space. This is inefficient and unnecessary and it is also becoming not sustainable for the increase of the data volume generated by HTS machines. This has in fact consequences in terms of available storage space and generated costs, and it is also hindering the benefits of genomic analysis in healthcare from reaching a larger portion of the population. The impact of the IT costs generated by the exponential growth of sequence data to be stored and analysed is currently one of the main challenges the scientific community and that the healthcare industry have to face (see Scott D. Kahn “On the future of genomic data”—Science 331, 728 (2011) and Pavlichin, D. S., Weissman, T., and G. Yona. 2013. “The human genome contracts again” Bioinformatics 29(17): 2199-2202). At the same time several are the initiatives attempting to scale genome sequencing from a few selected individuals to large populations (see Josh P. Roberts “Million Veterans Sequenced”—Nature Biotechnology 31, 470 (2013))

The transfer of genomic data is slow and inefficient because the currently used data formats are organized into monolithic files of up to several hundred Gigabytes of size which need to be entirely transferred at the receiving end in order to be processed. This implies that the analysis of a small segment of the data requires the transfer of the entire file with significant costs in terms of consumed bandwidth and waiting time. Often online transfer is prohibitive for the large volumes of the data to be transferred, and the transport of the data is performed by physically moving storage media such as hard disk drives or storage servers from one location to another.

These limitations occurring when employing state of the art approaches are overcome by the present invention.

Processing the data is slow and inefficient for to the fact that the information is not structured in such a way that the portions of the different classes of data and metadata required by commonly used analysis applications cannot be retrieved without the need of accessing the data in its totality. This fact implies that common analysis pipelines can require to run for days or weeks wasting precious and costly processing resources because of the need, at each stage of accessing, of parsing and filtering large volumes of data even if the portions of data relevant for the specific analysis purpose is much smaller.

These limitations are preventing health care professionals from timely obtaining genomic analysis reports and promptly reacting to diseases outbreaks. The present invention provides a solution to this need.

There is another technical limitation that is overcome by the present invention.

In fact the invention aims at providing an appropriate genomic sequencing data and metadata representation by organizing and partitioning the data so that the compression of data and metadata is maximized and several functionality such as selective access and support for incremental updates are efficiently enabled.

A key aspect of the invention is a specific definition of classes of data and metadata to be represented by an appropriate source model, coded (i.e. compressed) separately by being structured in specific layers. The most important achievements of this invention with respect to existing state of the art methods consist in:

• • the increase of compression performance due to the reduction of the information source entropy constituted by providing an efficient model for each class of data or metadata;
  • the possibility of performing selective accesses to portions of the compressed data and metadata for any further processing purpose directly in the compressed domain;
  • the possibility of defining user specified “labels” identifying genomic regions or sub-regions or aggregations of regions or sub-regions to enable efficient selective access to the compressed data by means of parsing a “labels list” contained in the genomic file header;
  • the possibility of implementing access control and protection to the different genomic regions or sub-regions identified by a label;
  • the possibility of incrementally (without the need of re-encoding) updating and adding encoded data and metadata with new sequencing data and/or metadata and/or new analysis results;
  • the possibility of efficiently processing data as soon as they are produced by the sequencing machine or alignment tools without the need of waiting the end of the sequencing or alignment process.

The present application discloses a method and system addressing the problem of efficient manipulation, storage and transmission of very large amounts of genomic sequencing data, by employing a structured access units approach combined with multiplexing techniques.

The present application overcomes all the limitations of the prior art approaches related to the functionality of genomic data accessibility, selective data protection, efficient processing of data subsets, transmission and streaming functionality combined with an efficient compression.

Today the most used representation format for genomic data is the Sequence Alignment Mapping (SAM) textual format and its binary correspondent BAM. SAM files are human readable ASCII text files whereas BAM adopts a block based variant of gzip. BAM files can be indexed to enable a limited modality of random access. This is supported by the creation of a separate index file.

The BAM format is characterized by poor compression performance for the following reasons:

• 1. It focuses on compressing the inefficient and redundant SAM file format rather than on extracting the actual genomic information conveyed by SAM files and using appropriate models for compressing it.
• 2. It employs a general purpose text compression algorithm such as gzip rather than exploiting the specific nature of each data source (the genomic information itself).
• 3. It lacks any concept and does not support any functionality related to data classification that would enable the implementation of mechanisms providing selective access to specific classes of genomic data.

A more sophisticated approach to genomic data compression that is less commonly used, but more efficient than BAM is CRAM (CRAM specification: https://samtools.github.io/hts-specs/CRAMv3.pdf). CRAM provides a more efficient compression for the adoption of differential encoding with respect to an existing reference (it partially exploits the data source redundancy), but it still lacks features such as incremental updates, support for streaming and selective access to specific classes of compressed data.

CRAM relies on the concept of the CRAM record. Each CRAM record encodes a single mapped or unmapped reads by encoding all the elements necessary to reconstruct it.

CRAM presents the following drawbacks and limitations that are solved and removed by the invention described in this document:

• 1. CRAM does not support data indexing and random access to data subsets sharing specific features. Data indexing is out of the scope of the specification (see section 12 of CRAM specification v 3.0) and it is implemented as a separate file. Conversely the approach of the invention described in this document employs a data indexing method that is integrated with the encoding process and indexes are embedded in the encoded (i.e. compressed) bit stream.
• 2. CRAM does not support the aggregation of the data related to several sequencing runs so that selective access is efficient and segregation of runs (i.e. the process of extracting the genomic information from the actual organic sample) is preserved. CRAM does provide the possibility to label reads as belonging to different groups, but this is provided on a read by read base and reads from different groups are then mixed in the file structure. In the present invention a method is described to structure the data so as to keep segregation among different sequencing runs so that efficient selective access is available.
• 3. CRAM is built by core data blocks that can contain any type of mapped reads (perfectly matching reads, reads with substitutions only, reads with insertions or deletions (also referred to as “indels”)). There is no notion of data classification and grouping of reads in classes according to the result of mapping with respect to a reference sequence. This means that all data need to be inspected even if only reads with specific features are searched. Such limitation is solved by the invention by classifying and partitioning data in classes before coding.
• 4. CRAM is based on the concept of encapsulating each read into a “CRAM record”. This implies the need to inspect each complete “record” when reads characterized by specific biological features (e.g. reads with substitutions, but without “indels”, or perfectly mapped reads) are searched. Conversely, in the present invention there is the notion of data classes coded separately in separate information layers and there is no notion of record encapsulating each read. This enables more efficient access to set of reads with specific biological characteristics (e.g. reads with substitutions, but without “indels”, or perfectly mapped reads) without the need of decoding each (block of) read(s) to inspect its features.
• 5. In a CRAM record each field in a record is associated to a specific flag and each flag must always have the same meaning as there is no notion of context since each CRAM record can contain any different type of data. This coding mechanism introduces redundant information and prevents the usage of efficient context based entropy coding.
  • Conversely in the present invention there is no notion of flag denoting data because this is intrinsically defined by the information “layer” the data belongs to. This implies a largely reduced number of symbols to be used and a consequent reduction of the information source entropy which results into a more efficient compression. Such improvement is possible because the use of different “layers” enables the encoder to reuse the same symbol across each layer with different meanings according to the context. In CRAM each flag must always have the same meaning as there is no notion of contexts and each CRAM record can contain any type of data.
• 6. In CRAM, substitutions, insertions and deletions are represented by using different syntax elements, option that increases the size of the information source alphabet and yields a higher source entropy. Conversely the approach of the disclosed invention uses a single alphabet and encoding for substitutions, insertions and deletions. This makes the encoding and decoding process simpler and produces a lower entropy source model which coding yields bitstreams characterized by high compression performance.
• 7. CRAM does not provide any mechanism to uniquely identify specific regions or sub regions of the genomic data or aggregations thereof. Apart from the definition of loci in terms of start and end positions on the reference sequence, according to the CRAM specification there is no way to:
  • label a region and access it using the defined label instead of the genomic start and end position. Start and end positions of the same genomic region may change if a new reference sequence is published, while a defined label would hide such change to any end user. The encoding and decoding system would take care of adapting the actual region identified by the label to the newly published reference sequence
  • aggregate several regions or sub-regions under the same label so that any end user would be able to select the required data via a single query not involving complex nested queries. The entire aggregation mechanism would be embedded in the encoding and decoding system as described in this document.
• 8. CRAM does not provide or support any mechanism to implement selective protection and access control relative to specific regions or sub regions of the genomic data or aggregations thereof, neither when such regions are pre-defined nor when they are specified by the user inserting appropriate “Labels”.

Beside CRAM also the other approaches to genomic data compression and processing present strong limitations to most of the desired functionality and do not support features that are provided by this invention disclosure as described and specified in the following of the document.

Genomic compression algorithms used in the state of the art can be classified into these categories:

• • Transform-based
    • LZ-based
    • Read reordering
  • Assembly-based
  • Statistical modeling

The first two categories share the disadvantage of not exploiting the specific characteristics of the data source (genomic sequence reads) and process the genomic data as string of text to be compressed without taking into account the specific properties of such kind of information (e.g. redundancy among reads, reference to an existing sample). Two of the most advanced toolkits for genomic data compression, namely CRAM and Goby (“Compression of structured high-throughput sequencing data”, F. Campagne, K. C. Dorff, N. Chambwe, J. T. Robinson, J. P. Mesirov, T. D. Wu), make a poor use of arithmetic coding as they implicitly model data as independent and identically distributed by a Geometric distribution. Goby is slightly more sophisticated since it converts all the fields to a list of integers and each list is encoded independently using arithmetic coding without using any context. In the most efficient mode of operation, Goby is able to perform some inter-list modeling over the integer lists to improve compression. These prior art solutions yield poor compression ratios and data structures that are difficult if not impossible to selectively access and manipulate once compressed. Downstream analysis stages can result to be inefficient and very slow due to the necessity of handling large and rigid data structures even to perform simple operation or to access selected regions of the genomic dataset.

A simplified vision of the relation among the file formats used in genome processing pipelines is depicted in FIG. 1. In this diagram file inclusion does not imply the existence of a nested file structure, but it only represents the type and amount of information that can be encoded for each format (i.e. SAM contains all information in FASTQ, but organized in a different file structure). CRAM contains the same genomic information as SAM/BAM, but it has more flexibility in the type of compression that can be used, therefore it is represented as a superset of SAM/BAM.

The use of multiple file formats for the storage of genomic information is highly inefficient and costly. Having different file formats at different stages of the genomic information life cycle implies a linear growth of utilized storage space even if the incremental information is minimal. Further disadvantages of prior art solutions are listed below.

• 1. Accessing, analysing or adding annotations (metadata) to raw data stored in compressed FastQ files or any combination thereof requires the decompression and recompression of the entire file with extensive usage of computational resources and time.
• 2. Retrieving specific subsets of information such as read mapping position, read variant position and type, indels position and types, or any other metadata and annotation contained in aligned data stored in BAM files requires to access the whole data volume associated to each read. Selective access to a single class of metadata is not possible with prior art solutions.
• 3. Prior art file formats require that the whole file is received at the end user before processing can start. For example the alignment of reads could start before the sequencing process has been completed relying on an appropriate data representation. Sequencing, alignment and analysis could proceed and run in parallel.
• 4. Prior art solution do not support structuring and are not able of distinguishing genomic data obtained by different sequencing processes according to their specific generation semantic (e.g. sequencing obtained at different time of the life of the same individual). The same limitation occurs for sequencing obtained by different types of biological samples of the same individual.
• 5. The protection by means of access control mechanisms (e.g. encryption, watermarking, digital signature, hashing) of entire or selected portions of the data is not supported by prior art solutions. For example the protection of:
  • a. selected DNA regions
  • b. only those sequences containing variants
  • c. chimeric sequences only
  • d. unmapped sequences only
  • e. regions or sub-regions or aggregations of regions or sub-regions identified by user defined Labels
  • f. specific metadata (e.g. origin of the sequenced sample, identity of sequenced individual, type of sample)
  • is not supported in files and data formats of prior art solutions.
• 6. The transcoding from sequencing data aligned to a given reference (i.e. a SAM/BAM file) to a new reference requires to process the entire data volume even if the new reference differs only by a single nucleotide position from the previous reference.

Therefore there is the clear need of an appropriate Genomic Information Storage Format (Genomic File Format) and Transport Mechanism that enable efficient compression, support selective access and protection functionality in the compressed domain, of local and remotely stored data and support the incremental addition of heterogeneous metadata in the compressed domain at all levels of the different stages of the genomic data processing.

The present invention provides a solution to the limitations of the state of the art by employing the method, devices and computer programs as claimed in the accompanying set of claims.

LIST OF FIGURES

FIG. 1 shows the main steps of a typical genomic pipeline and the related file formats.

FIG. 2 shows the mutual relationship among the most used genomic file formats

FIG. 3 shows how genomic sequence reads are assembled in an entire or partial genome via de-novo assembly or reference based alignment.

FIG. 4 shows how reads mapping positions on the reference sequence are calculated.

FIG. 5 shows how reads pairing distances are calculated.

FIG. 6 shows how pairing errors are calculated.

FIG. 7 shows how the pairing distance is encoded when a read mate pair is mapped on a different chromosome.

FIG. 8 shows how sequence reads can be generated from the first or second DNA strand of a genome.

FIG. 9 shows how a read mapped on strand 2 has a corresponding reverse complemented read on strand 1.

FIG. 10 shows the four possible combinations of reads composing a reads pair and the respective encoding in the rcomp layer.

FIG. 11 shows how “n type” mismatches are encoded in a nmis layer.

FIG. 12 shows an example of substitutions in a mapped read pair.

FIG. 13 shows how substitutions positions can be calculated either as absolute or differential values.

FIG. 14 shows how symbols encoding substitutions without IUPAC codes are calculated.

FIG. 15 shows how substitution types are encoded in the snpt layer.

FIG. 16 shows how symbols encoding substitutions with IUPAC codes are calculated.

FIG. 17 shows an alternative source model for substitution where only positions are encoded, but one layer per substitution type is used.

FIG. 18 shows how to encode substitutions, insertions and deletions in a reads pair of class I when IUPAC codes are not used.

FIG. 19 shows how to encode substitutions, insertions and deletions in a reads pair of class I when IUPAC codes are used.

FIG. 20 shows the structure of the Genomic Dataset Header of the genomic information data structure disclosed by this invention.

FIG. 21 shows how the Master Index Table contains the positions on the reference sequences of the first read in each Access Unit.

FIG. 22 shows an example of partial MIT showing the mapping positions of the first read in each pos AU of class P.

FIG. 23 shows how the Local Index Table in the layer header is a vector of pointers to the AUs in the payload.

FIG. 24 shows an example of Local Index Table.

FIG. 25 shows the functional relation between Master Index Table and Local Index Tables

FIG. 26 shows how Access Units are composed by blocks of data belonging to several layers. Layers are composed by Blocks subdivided in Packets.

FIG. 27 shows how a Genomic Access Unit of type 1 (containing positional, pairing, reverse complement and read length information) is packetized and encapsulated in a Genomic Data Multiplex.

FIG. 28 shows how Access Units are composed by a header and multiplexed blocks belonging to one or more layers of homogeneous data. Each block can be composed by one or more packets containing the actual descriptors of the genomic information.

FIG. 29 shows the structure of Access Units of type 0 which do not need to refer to any information coming from other access units to be accessed or decoded and accessed.

FIG. 30 shows the structure of Access Units of type 1.

FIG. 31 shows the structure of Access Units of type 2 which contain data that refer to an access unit of type 1. These are the positions of N bases in the encoded reads.

FIG. 32 shows the structure of Access Units of type 3 which contain data that refer to an access unit of type 1. These are the positions and types of mismatches in the encoded reads.

FIG. 33 shows the structure of Access Units of type 4 which contain data that refer to an access unit of type 1. These are the positions and types of mismatches in the encoded reads.

FIG. 34 shows the first five type of Access Units.

FIG. 35 shows that Access Units of type 1 refer to Access Units of type 0 to be decoded.

FIG. 36 shows that Access Units of type 2 refer to Access Units of type 0 and 1 to be decoded.

FIG. 37 shows that Access Units of type 3 refer to Access Units of type 0 and 1 to be decoded.

FIG. 38 shows that Access Units of type 4 refer to Access Units of type 0 and 1 to be decoded.

FIG. 39 shows the Access Units required to decode sequence reads with mismatches mapped on the second segment of the reference sequence (AU 0-2).

FIG. 40 shows how raw genomic sequence data that becomes available can be incrementally added to pre-encoded genomic data.

FIG. 41 shows how a data structure based on Access Units enables genomic data analysis to start before the sequencing process is completed.

FIG. 42 shows how new analysis performed on existing data can imply that reads are moved from AUs of type 4 to one of type 3.

FIG. 43 shows how newly generated analysis data are encapsulated in a new AU of type 8 and a corresponding index is created in the MIT.

FIG. 44 shows how to transcode data due to the publication of a new reference sequence (genome).

FIG. 45 shows how reads mapped to a new genomic region with better quality (e.g. no indels) are moved from AU of type 4 to AU of type 3

FIG. 46 shows how, in case new mapping location is found, (e.g. with less mismatches) the related reads can be moved from one AU to another of the same type.

FIG. 47 shows how selective encryption can be applied on Access Units of Type 4 only as they contain the sensible information to be protected.

FIG. 48 shows the data encapsulation in a genomic multiplex where one or more genomic datasets 482-483 contain Genomic streams 484 and streams of Genomic Datasets Mapping Table Lists 481, Genomic Dataset Mapping Tables 485, and Reference Identifiers Mapping Tables 487. Each genomic stream is composed by a Header 488 and Access Units 486. Access Units encapsulate Blocks 489 which are composed by Packets 4810.

FIG. 49 shows how raw genomic sequence data (499) or aligned genomic data (produced by element 491) are processed to be encapsulated in a Genomic Multiplex. The alignment (491) and reference genome construction (492) stages can be necessary to prepare the data for encoding. Data classes (498) generated by a data classification unit (494) can be further classified with respect to one or more transformed reference generated by a reference transformation unit (4919). The transformed classes (4918) are then sent to layers encoders (495-497). The generated layers (4911) are encoded by entropy coders (4912-4914) which generate Genomic Streams of Access Units (4915) fed to the Genomic Multiplexer (4916).

FIG. 50 shows how a genomic demultiplexer (500) extracts Genomic Streams (501) from the Genomic Multiplex (5010), one decoder per AU type (502-504) extracts the genomic layers which are then decoded (506-507) into various data classes (5011) which are used by class decoders (509) to reconstruct genomic formats such as for example FASTQ and SAM/BAM. When present in the multiplexed bitstream (5010) a genomic stream containing one or more reference transformations is decoded by an entropy decoder (504) to produce reference transformation descriptors (5012). Reference transformation descriptors are processed by a reference transformation unit (5013) to transform one or more “external” references to generate one or more transformed references (5014) to be used by the class decoders (509).

FIG. 51 shows the process of encoding sequence reads belonging to class U using a self-generated reference sequence using six layers of descriptors. Four layers are the same used for other classes P, N, M, I while two layers are specific to class U reads.

FIG. 52 shows how a label is built to aggregate genomic regions belonging to two different references.

FIG. 53 shows how an existing label can be updated in case new results of analysis require to add an additional region R4 to the existing ones (R1, R2 and R3).

FIG. 54 shows how the labeling mechanism can be used to implement access control and data protection on specific genomic regions or sub regions. The simple case uses one access control rule (AC) and one protection mechanism (e. g. encryption) for all genomic regions identified by one label.

FIG. 55 shows how the different genomic regions identified by the same label can be protected by several different access control rules (AC) and several different encryption keys.

FIG. 56 shows how an alternative encoding of reads of class U where a signed POS descriptor is used to encode the mapping position of a read on the computed reference FIG. 57 shows how half mapped read pairs can help in filling unknown regions of the reference sequence by assembling longer contigs with unmapped reads.

FIG. 58 shows the hierarchical structure of headers for genomic data stored following the structure described in this invention.

FIG. 59 shows how a device implementing the labeling mechanism described by this invention enables concurrent access to data related to several genomic regions when they are stored in different records of a database. This can happen either in presence of controlled access or not.

FIG. 60 shows how vectors of thresholds are used in encoders of classes N, M and I to generate separated subclasses of data

FIG. 61 provides an example of how reference transformations can change the class reads belong to when all or a subset of mismatches are removed (i.e. the read belonging to class M before transformation is assigned to class P after the transformation of the reference has been applied).

FIG. 62 shows how reference transformations can be applied to remove mismatches (MMs) from reads. In some cases reference transformations may generate new mismatches or change the type of mismatches found when referring to the reference before the transformation has been applied.

FIG. 63 The same reference transformation A0 can be used for all classes of data or different transformations AN, A_M, A_I are used for each class N, M, I

SUMMARY

The features of the claims below solve the problem of existing prior art solutions by providing

a method for selective access of regions of genomic data by employing labels, said labels comprising: an identifier of a reference genomic sequence (521), an identifier of said genomic regions (522), and an identifier of the data class (523) of said genomic data

In another aspect of the method said genomic data are sequences of genomic reads.

In another aspect of the method data classes can be of the following type or a subset of them:

• • “Class P” comprising genomic reads which do not present any mismatch with respect to a reference sequence
  • “Class N” comprising genomic reads including only mismatches in positions where the sequencing machine was not able to call any “base” and the number of said mismatches does not exceed a given threshold
  • “Class M” comprising genomic reads in which mismatches are constituted by positions where the sequencing machine was not able to call any base, named “n type” mismatches, and/or it called a different base than the reference sequence, named “s type” mismatches, and said numbers of mismatches do not exceed given thresholds for the number of mismatches of “n type”, of “s type” and a threshold obtained from a given function (f(n,s))
  • “Class I” when the genomic reads can possibly have the same type of mismatches of “Class M”, and in addition at least one mismatch of type: “insertion” (“i type”), “deletion” (“d type”), soft clips (“c type”), and wherein the numbers of mismatches for each type does not exceed the corresponding given thresholds and a threshold provided by a given function (w(n,s,i,d,c))
  • “Class U” comprising all reads that do not find any classification in the classes P, N, M, I

In another aspect of the method said genomic data are paired sequences of genomic reads.

In another aspect of the method said data class of paired reads can be of the following types or a subset of them:

• • “Class P” comprising genomic read pairs which do not present any mismatch with respect to a reference sequence
  • “Class N” comprising genomic reads pairs including only mismatches in positions where the sequencing machine was not able to call any “base” and said numbers of mismatches for each read do not exceed a given threshold
  • “Class M” comprising genomic read pairs including only mismatches in positions where the sequencing machine was not able to call any “base” and said numbers of mismatches for each read do not exceed a given threshold, named “n type” mismatches, and/or it called a different base than the reference sequence, named “s type” mismatches, and said numbers of mismatches does not exceed a given thresholds for the number of mismatches of “n type”, of “s type” and a threshold obtained from a given function (f(n,s))
  • “Class I” comprising read pairs which can possibly have the same type of mismatches of “Class M” pairs, and in addition at least one mismatch of type: “insertion” (“i type”) “deletion” (“d type”) soft clips (“c type”), and wherein the number of mismatches for each type does not exceed the corresponding given threshold and a threshold provided by a given function (w(n,s,i,d,c))
  • “Class HM” comprising read pairs for which only one read mate does not satisfy the matching rules for being classified in any of the classes P, N, M, I
  • Class “U” comprising all reads pairs for which both reads do not satisfy the matching rules for being classified in the classes P, N, M, I

In another aspect of the method said identifier of said genomic regions is comprised in a master index table.

In another aspect of the method said genomic data and said labels are entropy coded.

In another aspect of the method said master index table (4812) is comprised in a genomic dataset header (4813).

In another aspect of the method said regions of genomic data are dispersed among separate Access Units (524, 486).

In another aspect of the method the location of said regions of genomic data, in a file, is indicated in a local index table (525).

In another aspect of the method said labels are user specified.

In another aspect of the method said regions are protected and/or encrypted in a separate manner, without encrypting the whole genomic file.

In another aspect of the method said labels are stored in a genomic label list (GLL)

In another aspect the method further comprises encoding genomic data with selective access to regions of genomic data as previously defined.

In another aspect of the method said genomic label list is periodically retransmitted or updated in order to enable multiple synchronization points

In another aspect the method further comprises decoding a stream or a file of genomic data with selective access to regions of genomic data as previously defined.

The present invention further provides an apparatus for encoding genomic data as previously defined.

The present invention further provides an apparatus for decoding genomic data as previously defined.

The present invention further provides a storing mean for storing genomic data encoded as previously defined.

The present invention further provides a computer-readable medium comprising instructions that when executed cause at least one processor to perform the encoding method previously defined.

The present invention further provides a computer-readable medium comprising instructions that when executed cause at least one processor to perform the decoding method previously defined.

DETAILED DESCRIPTION

The present invention describes a labelling mechanism providing selective access and selective access control to genomic regions or sub-regions or aggregations of regions or sub-regions of compressed genomic data stored in a file format and/or the relevant access units to be used to store, transport, access and process genomic or proteomic information in the form of sequences of symbols representing molecules.

These molecules include, for example, nucleotides, amino acids and proteins. One of the most important pieces of information represented as sequence of symbols are the data generated by high-throughput genome sequencing devices.

The genome of any living organism is usually represented as a string of symbols expressing the chain of nucleic acids (bases) characterizing that organism. Current state of the art genome sequencing technology is able to produce only a fragmented representation of the genome in the form of several (up to billions) strings of nucleic acids associated to metadata (identifiers, level of accuracy etc.). Such strings are usually called “sequence reads” or “reads”.

The typical steps of the genomic information life cycle comprise Sequence reads extraction, Mapping and Alignment, Variant detection, Variant annotation and Functional and Structural Analysis (see FIG. 1).

Sequence reads extraction is the process —performed by either a human operator or a machine—of representation of fragments of genetic information in the form of sequences of symbols representing the molecules composing a biological sample. In the case of nucleic acids such molecules are called “nucleotides”. The sequences of symbols produced by the extraction are commonly referred to as “reads”. This information is usually encoded in prior art as FASTA files including a textual header and a sequence of symbols representing the sequenced molecules.

When the biological sample is sequenced to extract DNA of a living organism the alphabet is composed by the symbols (A,C,G,T,N).

When the biological sample is sequenced to extract RNA of a living organism the alphabet is composed by the symbols (A,C,G,U,N).

In case the IUPAC extended set of symbols, so called “ambiguity codes” are also generated by the sequencing machine, the alphabet used for the symbols composing the reads are (A, C, G, T, U, W, S, M, K, R, Y, B, D, H, V, N or −).

When the IUPAC ambiguity codes are not used, a sequence of quality score can be associated to each sequence read. In such case prior art solutions encode the resulting information as a FASTQ file. Sequencing devices can introduce errors in the sequence reads such as:

• 1. identification of a wrong symbol (i.e. representing a different nucleic acid) to represent the nucleic acid actually present in the sequenced sample; this is usually called “substitution error” (mismatch);
• 2. insertion in one sequence read of additional symbols that do not refer to any actually present nucleic acid; this is usually called “insertion error”;
• 3. deletion from one sequence read of symbols that represent nucleic acids that are actually present in the sequenced sample; this is usually called “deletion error”;
• 4. recombination of one or more fragments into a single fragment which does not reflect the reality of the originating sequence.

The term “coverage” is used in literature to quantify the extent to which a reference genome or part thereof can be covered by the available sequence reads. Coverage is said to be:

• • partial (less than 1×) when some parts of the reference genome are not mapped by any available sequence read
  • single (1×) when all nucleotides of the reference genome are mapped by one and only one symbol present in the sequence reads
  • multiple (2×, 3×, N×) when each nucleotide of the reference genome is mapped multiple times.

Sequence alignment refers to the process of arranging sequence reads by finding regions of similarity that may be a consequence of functional, structural, or evolutionary relationships among the sequences. When the alignment is performed with reference to a pre-existing nucleotides sequence referred to as “reference genome”, the process is called “mapping”. Sequence alignment can also be performed without a pre-existing sequence (i.e. reference genome) in such cases the process is known in prior art as “de novo” alignment. Prior art solutions store this information in SAM, BAM or CRAM files. The concept of aligning sequences to reconstruct a partial or complete genome is depicted in FIG. 3.

Variant detection (a.k.a. variant calling) is the process of translating the aligned output of genome sequencing machines, (sequence reads generated by NGS devices and aligned), to a summary of the unique characteristics of the organism being sequenced that cannot be found in other pre-existing sequence or can be found in a few pre-existing sequences only. These characteristics are called “variants” because they are expressed as differences between the genome of the organism under study and a reference genome. Prior art solutions store this information in a specific file format called VCF file.

Variant annotation is the process of assigning functional information to the genomic variants identified by the process of variant calling. This implies the classification of variants according to their relationship to coding sequences in the genome and according to their impact on the coding sequence and the gene product. This is in prior art usually stored in a MAF file.

The process of analysis of DNA (variant, CNV=copy number variation, methylation etc) strands to define their relationship with genes (and proteins) functions and structure is called functional or structural analysis. Several different solutions exist in the prior art for the storage of this data.

Genomic File Format

The invention disclosed in this document consists in the definition of a selective and controlled data access applied to a compressed data structure for representing, processing manipulating and transmitting genome sequencing data that differs from prior art solutions for at least the following aspects:

• • It does not rely on any prior art representation formats of genomic information (i.e. FASTQ, SAM).
  • It supports efficient handling and selective random access to data produced by multiple sequencing runs structured as multiple genomic datasets. Partitioning data from different sequencing runs into the same data structure enables analysts to simultaneously perform queries on them with great advantage for population genetics studies.
  • It implements a new original classification of the genomic data and metadata according to their specific characteristics. Sequence reads are mapped to a reference sequence and grouped in distinct classes according to the results of the alignment process. This results in data classes with lower information entropy that can be more efficiently encoded applying different specific compression algorithms such as Huffman coding, arithmetic coding (CABAC, CAVLAC), Asymmetric Numerical Systems, Lempel-Ziv and its derivations.
  • It implements a new method to associate data classes or subsets of data classes to specific genomic regions, or sub-regions or aggregations of regions or sub-regions, by means of user defined Labels that enable the selective access and protection of said compressed data classes corresponding to specific genomic regions or sub-regions or aggregations of regions or sub-regions.
  • It defines syntax elements and the related encoding/decoding process conveying the sequence reads and the alignment information into a representation which is more efficient to be processed for downstream analysis applications.

Classifying the reads according to the result of mapping and coding them using descriptors to be stored in layers (position layer, mate distance layer, mismatch type layer etc, etc, . . . ) present the following advantages:

• • A reduction of the information entropy when the different syntax elements are modelled by a specific source model which yields higher compression performance.
  • A more efficient access to data that are already organized in groups/layers that have a specific meaning for the downstream analysis stages and that can be accesses separately and independently directly in the compressed domain.
  • The presence of a modular data structure that can be updated incrementally by accessing only the required information without the need of decoding (i.e. decompressing) the whole data content.
  • The genomic information produced by sequencing machines is intrinsically highly redundant due to the nature of the information itself and to the need of mitigating the errors intrinsic in the sequencing process. This implies that the relevant genetic information which needs to be identified and analyzed (the variations with respect to a reference) is only a small fraction of the produced data. Prior art genomic data representation formats are not conceived to “isolate” the meaningful information at a given analysis stage from the rest of the information so as to make it promptly available and understandable by the analysis applications.
  • The solution brought by the disclosed invention is to represent genomic data in such a way that any relevant portion of data is readily available to the analysis applications without the need of accessing and decompressing the entirety of data and the redundancy of the data is efficiently reduced by efficient compression to minimize the required storage space and transmission bandwidth.

The key elements of the invention are:

• 1. The specification of a file format that “contains” structured and user-defined selectively accessible data elements called Access Units (AU) in compressed form. Such approach can be seen as the opposite of prior art approaches, SAM and BAM for instance, in which data are structured in non-compressed form and then the entire file is compressed. A first clear advantage of the approach is to be able to efficiently and naturally provide various forms of user-defined structured selective access to the data elements in the compressed domain which is impossible or extremely awkward in prior art approaches.
• 2. The structuring of the genomic information into specific “layers” of homogeneous data and metadata presents the considerable advantage of enabling the definition of different models of the information sources characterized by low entropy. Such models not only can differ from layer to layer, but can also differ inside each layer when the compressed data within layers are partitioned into Data Blocks included into Access Units. This structuring enables the use of the most appropriate compression for each class of data or metadata and portion of them with significant gains in coding efficiency versus prior art approaches.
• 3. The information is structured in Access Units (AU) so that any relevant subset of data used by genomic analysis applications is efficiently and selectively accessible by means of appropriate interfaces. These features enable faster access to data and yield a more efficient processing.
• 4. The definition of a Master Index Table and Local Index Tables enabling selective access to the information carried by the layers of encoded (i.e. compressed) data without the need to decode the entire volume of compressed data.
• 5. The possibility of accessing only the AUs that correspond to specific user defined genomic regions or sub-regions or aggregations of regions or sub-regions and data classes of interest by parsing a “Label List” present in the file header.
• 6. The possibility of providing different types of access control to different AUs and portions of data contained into the AU according to the user defined “Labels” identifying associated genomic regions.
• 7. The possibility of performing realignment of already aligned and compressed genomic data sets when they need to be re-aligned versus newly published reference genomes by performing an efficient transcoding of selected data portions in the compressed domain. The frequent release of new reference genomes currently requires resource consuming and time for the transcoding processes to re-align already compressed and stored genomic data with respect to the newly published references because all data volume need to be processed.

The method described in this document aims at exploiting the available a-priori knowledge on genomic data to define an alphabet of syntax elements with reduced entropy. In genomics the available knowledge is represented by an existing genomic sequence usually —but not necessarily —of the same species as the one to be processed. As an example, human genomes of different individuals differ only of a fraction of 1%. However, such small amount of data contain enough information to enable early diagnosis, personalized medicine, customized drugs synthesis etc. This invention aims at defining a genomic information representation format where the relevant information is efficiently accessible, access can be selectively controlled and data protected, the information is efficiently transportable and all such processing is performed handling compressed data structures.

The technical features used in the present invention are:

• 1. Partitioning genomic information generated by different sequencing runs into different genomic datasets in order to enable efficient data retrieval and processing when querying one or more of the available datasets.
• 2. Partition of the genome sequence data and metadata in “classes” sharing common features;
• 3. Definition of the structure of the genomic information carried by each data classes in which the genomic data is partitioned, into a sets of “layers” of descriptors in order to reduce the information entropy as much as possible;
• 4. Definition of a Master Index Table and Local Index Tables to enable selective access to the data classes and associated information by accessing only the desired layers of coded information (i.e. compressed) without the need to decode the entire coded genomic information;
• 5. Usage of different source models and entropy coders to code the syntax elements belonging to different layers of the data classes defined as specified in point 2;
• 6. Definition of specific mechanisms establishing a correspondence among dependent layers to enable selective access to the data without the need to decode all the layers if not necessary or desired;
• 7. Definition of a mechanism for labelling compressed data corresponding to specific genomic regions or sub-regions or aggregations of regions or sub-regions and corresponding data “classes” or subsets of data classes by “Labels” enabling efficient selective access;
• 8. Definition of mechanisms for the selective protection of specific genomic regions or sub-regions or aggregations of regions or sub-regions and corresponding data “classes” or subsets of data classes and any combination thereof.
• 9. Coding of the datasets or data “classes” with respect to one or more pre-existing or constructed reference sequences that can be further transformed to reduce the entropy of the sequence data representation.

In order to solve all the mentioned problems of the prior art in terms of efficient selective access and selective access control to specific data “classes”, specific genomic regions or sub-regions or aggregations of regions or sub-regions, while preserving efficient transmission and storing by means of an efficient compressed representation, the present invention application provides a specific data structure specification that implements appropriate data reordering into accessible units of homogeneous and/or semantically significant data enabling seamless access and processing required by state of the art genome data analysis applications.

In particular the present invention adopts a data structure based on the concept of Access Unit, “Labels” and the multiplexing of the relevant data, concepts which are absent from all state of the art genomic data formats.

Genomic data are structured and encoded into different Access Units. Hereafter follows a description of the genomic data that are contained into different Access Units and can be identified by “Labels” associating genomic data to specific genomic regions or sub-regions or aggregations of regions or sub-regions versus reference genomes.

Genomic Data Classification According to Matching Rules

The sequence reads generated by sequencing machines are classified by the disclosed invention into five different “classes” according to the matching results of the alignment with respect to one or more pre-existing reference sequences.

When aligning a DNA sequence of nucleotides with respect to a reference sequence the following cases can be identified:

• 1. A region in the reference sequence is found to match the sequence read without any error (i.e. perfect mapping). Such sequence of nucleotides is referenced to as “perfectly matching read” or denoted as “Class P”.
• 2. A region in the reference sequence is found to match the sequence read with a type and a number of mismatches determined only by the number of positions in which the sequencing machine generating the read was not able to call any base (or nucleotide). Such type of mismatches are denoted by an “N” the letter used to indicate an undefined nucleotide base. In this document this type of mismatch referred to as “n type” mismatch. Such sequences is referenced to as “N mismatching reads” or “Class N”. Once the read is classified to belong to “Class N” it is useful to limit the degree of matching inaccuracy to a given upper bound and set a boundary between what is considered a valid matching and what it is not. Therefore, the reads assigned to Class N are also constrained by setting a threshold (MAXN) that defines the maximum number of undefined bases (i.e. bases called as “N”) that a read can contain. Such classification implicitly defines the required minimum matching accuracy (or maximum degree of mismatch) that all reads belonging to Class N shares when referred to the corresponding reference sequence, which constitute an useful criterion for applying selective data searches to the compressed data.
• 3. A region in the reference sequence is found to match the sequence read with types and number of mismatches determined by the number of positions in which the sequencing machine generating the read was not able to call any nucleotide base, if present (i.e. “n type” mismatches), plus the number of mismatches in which a different base, than the one present in the reference, has been called. Such type of mismatch denoted as “substitution” is also called Single Nucleotide Variation (SNV) or Single Nucleotide Polymorphism (SNP). In this document this type of mismatch is also referred to as “s type” mismatch. The sequence read is then referenced to as “M mismatching reads” and assigned to “Class M”. Like in the case of “Class N”, also for all reads belonging to “Class M” it is useful to limit the degree of matching inaccuracy to a given upper bound, and set a boundary between what is considered a valid matching and what it is not. Therefore, the reads assigned to Class M are also constrained by defining a set of thresholds, one for the number “n” of mismatches of “n type” (MAXN) if present, and another for the number of substitutions “s” (MAXS). A third constraint is a threshold defined by any function of both numbers “n” and “s”, f(n,s). Such third constraint enable to generate classes with an upper bound of matching inaccuracy according to any meaningful selective access criterion. For instance, and not as a limitation, f(n,s) can be (n+s)1/2 or (n+s) or any linear or non-linear expression that sets a boundary to the maximum matching inaccuracy level that is admitted for a read belonging to “Class M”. Such boundary constitutes a very useful criterion for applying the desired selective data searches to the compressed data when analyzing sequence reads for various purposes because it makes possible to set a further boundary to any possible combination of the numbers “n” of “n type” mismatches and “s” of “s type” mismatches (substitutions) beyond the simple threshold applied to the one type or to the other.
• 4. A fourth class is constituted by sequencing reads presenting at least one mismatch of any type among “insertion”, “deletion” (a.k.a. indels) and “clipped”, plus, if present, any mismatches type belonging to class N or M. Such sequence is referenced to as “I mismatching reads” and assigned to “Class I”. Insertions are constituted by an additional sequence of one or more nucleotides not present in the reference, but present in the read sequence. In this document this type of mismatch is referred to as “i type” mismatch. In literature when the inserted sequence is at the edges of the sequence it is also referred to as “soft clipped” (i.e. the nucleotides are not matching the reference but are kept in the aligned reads contrarily to “hard clipped” nucleotides which are discarded). In this document this type of mismatch is referred to as “c type” mismatch. Keeping or discarding nucleotides is a decisions taken by the aligner stage and not by the classifier of reads disclosed in this invention which receives and processes the reads as they are determined by the sequencing machine or by the following alignment stage. Deletion are “holes” (missing nucleotides) in the read with respect to the reference. In this document this type of mismatch is referred to as “d type” mismatch. Like in the case of classes “N” and “M” it is possible and appropriate to define a limit to the matching inaccuracy. The definition of the set of constraints for “Class I” is based on the same principles used for “Class M” and is reported in Table 1 in the last table lines. Beside a threshold for each type of mismatch admissible for class I data, a further constraint is defined by a threshold determined by any function of the number of the mismatches “n”, “s”, “d”, “i” and “c”, w(n,s,d,i,c). Such additional constraint make possible to generate classes with an upper bound of matching inaccuracy according to any meaningful user defined selective access criterion. For instance, and not as a limitation, w(n,s,d,i,c) can be (n+s+d+i+c)1/5 or (n+s+d+i+c) or any linear or non-linear expression that sets a boundary to the maximum matching inaccuracy level that is admitted for a read belonging to “Class I”. Such boundary constitutes a very useful criterion for applying the desired selective data searches to the compressed data when analyzing sequence reads for various purposes because it enables to set a further boundary to any possible combination of the number of mismatches admissible in “Class I” reads beyond the simple threshold applied to each type of admissible mismatch.
• 5. A fifth class includes all reads that do now find any matching considered valid (i.e not satisfying the set of matching rules defining an upper bound to the maximum matching inaccuracy as specified in Table 1) for each data class when referring to the reference sequence. Such sequences are said to be “Unmapped” when referring to the reference sequences and are classified as belonging to the “Class U”.

Classification of Read Pairs According to Matching Rules

The classification specified in the previous section concerns single sequence reads. In the case of sequencing technologies that generates read in pairs (i.e. Illumina Inc.) in which two reads are known to be separated by an unknown sequence of variable length, it is appropriate to consider the classification of the entire pair to a single data class. A read that is coupled with another is said to be its “mate”.

If both paired reads belong to the same class the assignment to a class of the entire pair is obvious, the entire pair is assigned to the same class for any class (i.e. P, N, M, I, U). In the case the two reads belong to a different class, but none of them belongs to the “Class U”, then the entire pair is assigned to the class with the highest priority defined according to the following expression:

P<N<M<I

in which “Class P” has the lowest priority and “Class I” has the highest priority.

In case only one of the reads belongs to “Class U” and its mate to any of the Classes P, N, M, I a sixth class is defined as “Class HM” which stands for “Half Mapped”.

The definition of such specific class of reads is motivated by the fact that it is used for attempting to determine gaps or unknown regions existing in reference genomes (a.k.a. little known or unknown regions). Such regions are reconstructed by mapping pairs at the edges using the pair read that can be mapped on the known regions. The unmapped mate is then used to build the so called “contigs” of the unknown region as it is shown in FIG. 57. Therefore providing a selective access to only such type of read pairs greatly reduces the associated computation burden enabling much efficient processing of such data originated by large amounts of data sets that using the state of the art solutions would require to be entirely inspected.

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

TABLE 1
Type of mismatches and set of constrains that each sequence reads must satisfy to be
classified in the data classes defined in this invention disclosure.

Number and types of mismatches found when
matching a read with a reference sequence

Number of				Number of
unknown	Number of	Number of	Number of	clipped	Set of matching	Assignement
bases (“N”)	substitutions	deletions	Insertions	bases	accuracy constraints	Class
0	0	0	0	0	0	P
n > 0	0	0	0	0	n ≤ MAXN	N
					n > MAXN	U
n ≥ 0	s > 0	0	0	0	n ≤ MAXN and	M
					s ≤ MAXS and
					f(n, s) ≤ MAXM
					n > MAXN or	U
					s > MAXS or
					f(n, s) > MAXM
n ≥ 0	s ≥ 0	d ≥ 0*	i ≥ 0*	c ≥ 0*	n ≤ MAXN and	I

	*At least one mismatch	s ≤ MAXS and
	of type d, i, c must be resent	d ≤ MAXD and
	(i.e. d > 0 or i > 0 or > 0)	i ≤ MAXI and
		c ≤ MAXC
		w(n, s, d, i, c) ≤
		MAXTOT

	d ≥ 0	i ≥ 0	c ≥ 0	n > MAXN or	U
				s > MAXS or
				d > MAXD or
				i > MAXI or
				c > MAXC
				w(n, s, d, i, c) >
				MAXTOT

Matching Rules Partition of Sequence Read Data Classes N, M and I into Subclasses with Different Degrees of Matching Accuracy

The data classes of type N, M and I as defined in the previous sections can be further decomposed into an arbitrary number of distinct sub-classes with different degrees of matching accuracy. Such option is an important technical advantage in providing a finer granularity and as consequence a much more efficient selective access to each data class. As an example and not as a limitation, to partition the Class N into a number k of subclasses (Sub-Class N₁, . . . , Sub-Class N_k) it is necessary to define a vector with the corresponding components MAXN₁, MAXN₂, MAXN_(k-1), MAXN_(k), with the condition that MAXN₁<MAXN₂< . . . <MAXN_(k-1)<MAXN and assign each read to the lowest ranked sub-class that satisfy the constrains specified in Table 1 when evaluated for each element of the vector. This is shown in FIG. 60 where a data classification unit 601 contains Class P, N, M, I U, HM encoder and encoders for annotations and metadata. Class N encoder is configured with a vector of thresholds, MAXN₁ to MAXN_k 602 which generates k subclasses of N data (606).

In the case of the classes of type M and I the same principle is applied by defining a vector with the same properties for MAXM and MAXTOT respectively and use each vector components as threshold for checking if the functions f(n,s) and w(n,s,d,i,c) satisfy the constraint. Like in the case of sub-classes of type N, the assignment is given to the lowest sub-class for which the constraint is satisfied. The number of sub-classes for each class type is independent and any combination of subdivisions is admissible. This is shown in FIG. 60 where a Class M encoder and a Class I encoder are configured respectively with a vector of thresholds MAXM₁ to MAXM_j (603) and MAXTOT₁ to MAXTOT_h (604). The two encoders generate respectively j subclasses of M data (607) and h subclasses of I data (608). When two reads in a pair are classified in the same sub-class, then the pair belongs to the same sub-class.

When two reads in a pair are classified into sub-classes of different classes, then the pair belongs to the sub-class of the class of higher priority according to the following expression:

N<M<I

where N has the lowest priority and I has the highest priority.

When two reads belong to different sub-classes of one of classes N or M or I, then the pair belongs to the sub-class with the highest priority according to the following expressions:

N₁<N₂< . . . <N_k

M₁<M₂< . . . <M_j

I₁<I₂< . . . <I_h

where the highest index has the highest priority.

Transformations of the “External” Reference Sequences

The mismatches found for the reads classified in the classes N, M and I can be used to create “transformed references” to be used to compress more efficiently the read representation. Reads classified as belonging to the Classes N, M or I (with respect to the pre-existing (i.e. “external”) reference sequence denoted as RS₀) can be coded with respect to the “transformed” reference sequence RS₁ according to the occurrence of the actual mismatches with the transformed reference. For example if read^M_in belonging to Class M (denoted as the i^th read of class M) containing mismatches with respect to the reference sequence RS_n, then after “transformation” read^M_in=read^P_i(n+1) can be obtained with A(Ref_n)=Ref_n+1 where A is the transformation from reference sequence RS_n to reference sequence RS_n+1.

FIG. 61 shows an example on how reads containing mismatches (belonging to Class M) with respect to reference sequence 1 (RS₁) can be transformed into perfectly matching reads with respect to the reference sequence 2 (RS₂) obtained from RS₁ by modifying the bases corresponding to the mismatch positions. They remain classified and they are coded together the other reads in the same data class access unit, but the coding is done using only the descriptors and descriptor values needed for a Class P read. This transformation can be denoted as:

RS₂=A(RS₁)

When the representation of the transformation A which generates RS₂ when applied to RS₁ plus the representation of the reads versus RS₂ corresponds to a lower entropy than the representation of the reads of class M versus RS₁, it is advantageous to transmit the representation of the transformation A and the corresponding representation of the read versus RS₂ because an higher compression of the data representation is achieved.

The coding of the transformation A for transmission in the compressed bitstream requires the definition of two additional syntax elements as defined in the table below.

Syntax
elements	Semantic	Comments
rftp	Reference	position of difference between
	transformation	reference and contig
	position	used for prediction
rftt	Reference	type of difference between reference and
	transformation	contig used for prediction. Same syntax
	type	described for the snpt descriptor defined
		below.

FIG. 62 shows an example on how a reference transformation is applied to reduce the number of mismatches to be coded on the mapped reads.

It has to be observed that, in some cases the transformation applied to the reference:

• • May introduce mismatches in the representations of the reads that were not present when referring to the reference before applying the transformation.
  • May modify the types of mismatches, a read may contain A instead of G while all other reads contain C instead of G), but mismatches remain in the same position.
  • Different data classes and subsets of data of each data class may refer to the same transformed reference sequence or to reference sequences obtained by applying different transformations to the same pre-existing reference sequence.

FIG. 61 shows an example on how reads can change the type of coding from a data class to another by means of the appropriate set of descriptors (e.g. using the descriptors of a Class P to code a read from Class M) after a reference transformation is applied and the read is represented using the transformed reference. This occurs for example when the transformation changes all bases corresponding to the mismatches of a read in the bases actually present in the read, thus virtually transforming a read belonging to Class M (when referring to the original non transformed reference sequence) into a virtual read of Class P (when referring to the transformed reference). The definition of the set of descriptors used for each class of data is provided in the following sections.

FIG. 63 shows how the different classes of data can use the same “transformed” reference R₁=A₀(R₀) (630) to re-encode the reads or different transformations A_N (631), A_M (632), A_I (633) can be separately applied to each class of data

Genomic Data Headers for Global Parameters

The data structure of said genomic data requires the storage of global parameters and metadata to be used by the decoding engine. These data are organized in the following structures: For file based storage:

• • Datasets Multiplex Header
  • Dataset Header
  • Descriptors Layer Header
  • Block Header

The hierarchical relationship among these headers is shown in FIG. 58.

For transport in a streaming scenario:

• • Datasets Mapping Tables List
  • Datasets Mapping Table
  • Transport Block Header
  • Packet Header

A dataset is defined as the ensemble of coding elements needed to reconstruct the genomic information related to a single genomic sequencing run and all the following analysis. If the same genomic sample is sequenced twice in two distinct runs, the obtained data will be encoded in two distinct datasets.

Datasets Multiplex Header

This is the data structure prepended to one or more datasets aggregated in a “Multiplex”.

Syntax	Description
Datasets_ Multiplex _Header {
Multiplex_id	Label to identify this Datasets Multiplex
	from any other Datasets Multiplex.
Version_number	Version number of the Dataset
	Multiplex. The version number shall be
	incremented by a unit whenever the
	definition of the Datasets Multiplex
	changes.
List_number	Number of the current datasets list.
gd_number	Number of datasets composing
	the datasets Multiplex.
for (i=0; i<gd_number;i++) {
genomic_dataset_ID	Field identifying the dataset. This field
	shall not take any single value more
	than once within one version of the
	Dataset List
}
Metadata	Data structure carrying metadata to
	be used for application-specific
	processing such as data analysis and
	content protection.
}

This is the data structure prepended to an encoded dataset.

TABLE 2
Genomic Dataset Header structure.

Element	Type	Description
Dataset_ID	Byte array	Unique identifier
		for the encoded
		content
Major_Brand	Byte array	Major + Minor
Minor_Version	Byte array	version of the
		encoding
		algorithm
Header Size	Integer	Size in bytes of the
		entire encoded
		content
Reads Length	Integer	Size of reads in
		case of constant
		reads length. A
		special value (e.g.
		0) is reserved for
		variable reads
		length
Ref count	Integer	Number of
		reference
		sequences used
Access Units counters	Byte array	Total Number of
	(e.g.	encoded Access
	integers)	Units per
		reference
		sequence
Ref ids	Byte array	Unique identifiers
		for reference
		sequences
Ref_count	Integer	number of
		references
for (i=0; i<Ref_count; i++) {
Reference_genome:Ref_ID	string:string	Unambiguous ID,
		as a characters
		string, identifying
		the reference
		sequence(s) used
		in this Dataset
}
for (i=0; i<Ref_count; i++) {
Ref blocks	Byte array	Number of
		encoded blocks
		per each reference
}
Dataset label size	Integer	The size of the
		following element
Dataset label	String	A string of
		character used to
		identify the
		dataset
Dataset type	Integer	The type of data
		encoded in the
		dataset (e.g.
		aligned, not
		aligned)
Master index table	Byte array	This is a
Alignment positions of first read in each block (Access Unit).		multidimensional
l.e. smaller position of the first read on the reference genome		array supporting
per each block of the six classes		random access to
1 per pos class (six) per reference		Access Units.
Label List	Byte array	This is a list of
Sub-part of the Genomic Dataset Header indicating	(e.g.	Labels, each one
number of Labels	integers)	represented as a
for each Label:		multidimensional
the Label ID		array in order to
the number of reference sequences concerned		support selective
by the label		access to specific
for each reference sequence		genomic regions
the reference identifier		or sub-regions or
the number of regions covered by the		aggregations of
label,		regions or sub-
for each region:		regions.
the class ID
the start position in the
genomic range
the end position in the
genomic range
Start position and end position can be replaced by “block
numbers”, composing, together with reference sequence ID and
class ID, a three dimensional vector addressing the coordinates
of the Master Index Table.
Parameters set	Byte array	Encoding
		parameters used
		to configure the
		encoding process
		and sent to the
		decoder.

Descriptors Layer Header

Descriptors (a.k.a. syntax elements) are described in the following sections of this document and are the building blocks of the genomic information representation described by this invention. They are organized in layers (a.k.a. descriptors streams) of homogeneous elements partitioned according to the specific statistical properties of each descriptor. This has the advantage of reducing the entropy of each layer and improving compression efficiency.

Each layer is prepended by the Descriptors Layer Header described below.

Syntax	Description
Descriptors_Layer_Header {
Descriptors_Layer_ID	Descriptors layer ID, table specified in this specification
Num_Of_Blocks	Number of Blocks in the Descriptors Layer
Label size	Size of the human readable label
Label	(Human-Readable) Label
Flag	Flag used to interpret the following metadata
Local Index Table	The Local Index Table structure as described in this invention
Metadata	Data structure carrying metadata to be used for application-
	specific processing such as data analysis and content protection.
}

Block Header

Every Descriptors Layer is composed by one or multiple Genomic Data Blocks. One or more Blocks from different Layers compose an Access Unit, depending on the Class of data.

An Access Unit is a set of Genomic Blocks that can be decoded either independently from other Access Units by using only globally available data (e.g. decoder configuration) or by using information contained in other Access Units.

Syntax	Semantic
Block_Header {
Descriptors_Layer_ID	Unambiguously identifies the descriptors stream. Same as
	Descriptors_Layer_ID in Descriptor Layer Header
Block size (BS)	Number of bytes composing Block, including this header and
	payload, and excluding padding (total Block size will be BS +
	padding size).
}

Definition of the Information Necessary to Represent Sequence Reads into Layers of Descriptors

Once the classification of reads is completed with the definition of the Classes, further processing consists in defining a set of distinct syntax elements which represent the remaining information enabling the reconstruction of the DNA read sequence when represented as being mapped on a given reference sequence.

A sequence read (e.g. a DNA segment) referred to a given reference sequence can be fully expressed by:

• • The starting position on the reference sequence pos (292).
  • A flag signaling if the read has to be considered as a reverse complement versus the reference rcomp (293).
  • A distance to the mate pair in case of paired reads pair (294).
  • The value of the read length (295) in case of the sequencing technology produces variable length reads. In case of constant reads length the read length associated to each reads can obviously be omitted and can be stored in the Genomic Dataset Header.
  • For each mismatch:
    • Mismatch position nmis (300) for class N, snpp (311) for class M, and indp (321) for class I)
    • Mismatch type (not present in class N, snpt (312) in class M, indt (322) in class I)
  • Flags (296) indicating specific characteristics of the sequence read such as:
    • template having multiple segments in sequencing
    • each segment properly aligned according to the aligner
    • unmapped segment
    • next segment in the template unmapped
    • signalization of first or last segment
    • quality control failure
    • PCR or optical duplicate
    • secondary alignment
    • supplementary alignment
  • Soft clipped nucleotides string (323) when present for class I
  • Flag indicating the reference used for alignment and compression (e.g. internal reference for class U) if applicable (descriptor rtype).
  • For class U, descriptor indc identifies those parts of the reads (typically the edges) that do not match, with a specified set of matching accuracy constraints, with the “internal” reference sequences.
  • Descriptor ureads is used to encode verbatim the reads that cannot be mapped on any available reference being it “external” (i.e pre-existing like an actual reference genome) or a “internal” reference sequence.

This classification creates groups of descriptors (syntax elements) that can be used to univocally represent genome sequence reads. The table below summarizes the syntax elements needed for each class of reads aligned with “pre-existing” (i.e. “external”) or “constructed” (i.e. “internal”) references.

TABLE 3
Defined layers per class of data.

	P	N	M	I	U	HM

	pos	X	X	X	X	X	X
	pair	X	X	X	X	X
	rcomp	X	X	X	X		X
	flags	X	X	X	X		X
	rlen	X	X	X	X		X
	nmis		X
	snpp			X		X
	snpt			X		X
	indp				X		X
	indt				X		X
	indc				X		X
	ureads					X	X
	rtype					X

Reads belonging to class P are characterized and can be perfectly reconstructed by only a position, a reverse complement information and an offset between mates in case they have been obtained by a sequencing technology yielding mated pairs, some flags and a read length.

The next section details how these descriptors are defined for classes P, N, M and I while for class U they are described in a later section.

Class HM is applied to read pairs only and it is a special case where one read belongs to class P, N, M or I and the other to class U.

Position Descriptors Layer

In each Access Unit, only the mapping position of the first encoded read is stored in the AU header as absolute position on the reference genome. All the other positions are expressed as a difference with respect to the previous position and are stored in a specific layer. This modeling of the information source, defined by the sequence of read positions, is in general characterized by a reduced entropy particularly for sequencing processes generating high coverage results. Once the absolute position of the first alignment has been stored, all positions of other reads are expressed as difference (distance) with respect to the first one.

For example FIG. 4 shows how after encoding the starting position of the first alignment as position “10000” on the reference sequence, the position of the second read starting at position 10180 is coded as “180”. With high coverage data (>50×) most of the descriptors of the position vector will show very high occurrences of low values such as 0 and 1 and other small integers. FIG. 10 shows how the positions of three read pairs are encoded in a pos Layer.

The same source model is used for the positions of reads belonging to classes N, M, P and I. In order to enable any combination of selective access to the data, the positions of reads belonging to the four classes are encoded in separate layers as depicted in Table I.

Reverse Complement Descriptor Layer

Each read of the read pairs produced by sequencing technologies can be originated from either genome strands of the sequenced organic sample. However, only one of the two strands is used as reference sequence. FIG. 8 shows how in a reads pair one read (read 1) can be originated from one strand and the other (read 2) can be originated from the other strand.

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

In case of coupled reads, four are the possible combinations of direct and reverse complement mate pairs. This is shown in FIG. 10. The rcomp layer codes the four possible combinations.

The same coding is used for the reverse complement information of reads belonging to classes P, N, M, I. In order to enable enhanced selective access to the data, the reverse complement information of reads belonging to the four classes are coded in different layers as depicted in Table 3.

Pairing Descriptors Layer

The pairing descriptor is stored in the pair layer. Such layer stores descriptors encoding the information needed to reconstruct the originating reads pairs, when the employed sequencing technology produces reads by pairs. Although at the date of the disclosure of the invention the vast majority of sequencing data is generated by using a technology generating paired reads, it is not the case of all technologies. This is the reason for which the presence of this layer is not necessary to reconstruct all sequencing data information if the sequencing technology of the genomic data considered does not generate paired reads information.

Definitions

• • mate pair: read associated to another read in a read pair (e.g. Read 2 is the mate pair of Read 1 in the example of FIG. 4)
  • pairing distance: number of nucleotide positions on the reference sequence which separate one position in the first read (pairing anchor, e.g. last nucleotide of first read) from one position of 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 number of nucleotide positions.
  • position pairing distance (PPD): the PPD is a way to express a pairing distance in terms of the number of reads separating one read from its respective mate present in a specific position descriptor layer.
  • most probable position pairing distance (MPPPD): is the most probable number of reads separating one read from its mate pair present in a specific position descriptor layer.
  • position pairing error (PPE): is defined as the difference between the MPPD or MPPPD and the actual position of the mate.
  • pairing anchor: position of first read last nucleotide in a pair used as reference to calculate the distance of the mate pair in terms of number of nucleotide positions or number of read positions.

FIG. 5 shows how the pairing distance among read pairs is calculated.

The pair descriptor layer is the vector of pairing errors calculated as number of reads to be skipped to reach the mate pair of the first read of a pair with respect to the defined decoding pairing distance.

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

The same descriptors are used for the pairing information of reads belonging to classes N, M, P and I. In order to enable the selective access to the different data classes, the pairing information of reads belonging to the four classes are encoded in different layer as depicted in.

Pairing Information in Case of Reads Mapped on Different References

In the process of mapping sequence reads on a reference sequence it is not uncommon to have the first read in a pair mapped on one reference (e.g. chromosome 1) and the second on a different reference (e.g. chromosome 4). In this case the pairing information described above has to be integrated by 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 mapped on two different sequences (different values indicate if read1 or read2 are mapped on the sequence that is not currently encoded)
2. a unique reference identifier referring to the reference identifiers encoded in the Genomic Dataset Header structure as described in Table 2.
3. a third element containing the mapping information on the reference identified at point 2 and expressed as offset with respect to the last encoded position.

FIG. 7 provides an example of this scenario.

In FIG. 7, since Read 4 is not mapped on the currently encoded reference sequence, the genomic encoder signals this information by crafting additional descriptors in the pair layer. In the example shown in FIG. 7 Read 4 of pair 2 is mapped on reference no. 4 while the currently encoded reference is no. 1. This information is encoded using 3 components:

1) One special reserved value is encoded as pairing distance (in this case 0xffffff)
2) A second descriptor provides a reference ID as listed in the Genomic Dataset Header (in this case 4)
3) The third element contains the mapping information on the concerned reference (170).

Mismatch Descriptors for Class N Reads

Class N includes all reads in which only “n type” mismatches are present, at the place of an A, C, G or T base a N is found as called base. All other bases of the read perfectly match the reference sequence.

FIG. 11 shows how:

• • the positions of “N” in read 1 are coded as
  • absolute position in read 1 or
  • as differential position with respect to the previous “N” in the same read. the positions of “N” in read 2 are coded as
  • absolute position in read 2+read 1 length or
  • differential position with respect to the previous N In the nmis layer, the coding of each reads pair is terminated by a special “separator” symbol.

Encoding Substitutions (Mismatches or SNPs)

A substitution is defined as the presence, in a mapped read, of a different nucleotide with respect to the one that is present in the reference sequence at the same position (see FIG. 12).

Each substitution can be encoded as

• • “position” (snpp layer) and “type” (snpt layer). See FIG. 13, FIG. 14, FIG. 16 and FIG. 15.
  • OR
  • “position” only but using one snpp layer per mismatch type. See FIG. 17

Substitutions Positions

A substitution position is calculated as for the values of the nmis layer, i.e.: In read 1 substitutions are encoded

• • as absolute position in read 1 OR
  • as differential position with respect to the previous substitution in the same read In read 2 substitutions are encoded

In read 1 substitutions are encoded

• • as absolute position in read 2+read 1 length OR
  • as differential position with respect to the previous substitution FIG. 13 shows how substitutions positions are encoded in layer snpp. Substitutions positions can be calculated either as absolute or as differential values.

In the snpp layer, the encoding of each reads pair is terminated by a special “separator” symbol.

Substitutions Types Descriptors

For class M (and I as described in the next sections), mismatches are coded by an index (moving from right to left) from the actual symbol present in the reference to the corresponding substitution symbol present in the read {A, C, G, T, N, Z}. For example if the aligned read presents a C instead of a T which is present at the same position in the reference, the mismatch index will be denoted as “4”. The decoding process reads the encoded syntax element, the nucleotide at the given position on the reference and moves from left to right to retrieve the decoded symbol. E.g. a “2” received for a position where a G is present in the reference will be decoded as “N”. FIG. 14 shows all the possible substitutions and the respective encoding symbols when IUPAC ambiguity codes are not used and FIG. 15 provides an example of encoding of substitutions types in the snpt layer.

In case of presence of IUPAC ambiguity codes, the substitution indexes change as shown in FIG. 16.

In case the encoding of substation types described above presents high information entropy, an alternative method of substitution encoding consists in storing only the mismatches positions in separate layers, one per nucleotide, as depicted in FIG. 17.

Encoding of Insertions and Deletions

For class I, mismatches and deletions are coded by an indexes (moving from right to left) from the actual symbol present in the reference to the corresponding substitution symbol present in the read: {A, C, G, T, N, Z}. For example if the aligned read presents a C instead of a T present at the same position in the reference, the mismatch index will be “4”. In case the read presents a deletion where a A is present in the reference, the coded symbol will be “5”. The decoding process reads the coded syntax element, the nucleotide at the given position on the reference and moves from left to right to retrieve the decoded symbol. E.g. a “3” received for a position where a G is present in the reference will be decoded as “Z” which indicates the presence of a deletion in the sequence read.

Inserts are coded as 6, 7, 8, 9, 10 respectively for inserted A, C, G, T, N.

In case of adoption of the IUPAC ambiguity codes the substitution mechanism results to be exactly the same however the substitution vector is extended as: S={A, C, G, T, N, Z, M, R, W, S, Y, K, V, H, D, B} and insertions use different codes: 16, 17, 18, 19, 20.

FIG. 18 and FIG. 19 show examples of how to encode substitutions, inserts and deletions in a reads pair of class I.

The following structures of file format, access units and multiplexing are described referring to the coding elements disclosed here above. However, the access units, the file format and the multiplexing produce the same technical advantage also with other and different algorithms of source modeling and genomic data compression.

Construction of “Internal” References for Unmapped Reads of “Class U” and “Class HM”

In the case of the reads belonging to Class U or the unmapped pair of “Class HM” since they cannot be mapped to any “external” reference sequence satisfying the specified set of matching accuracy constraints for belonging to any of the classes P, N, M, or I, one or more “internal” reference sequences are constructed and used for the compressed representation of the reads belonging to these data classes.

Several approaches are possible to construct appropriate “internal” references such as for instance and not as limitation:

• • the partitioning of the non-mapped reads into clusters containing reads that share a common contiguous genomic sequence of at least a minimal size (signature). Each cluster can be uniquely identified by its signature.
  • the sorting of reads in any meaningful order (e.g. lexicographic order) and the use of the last N reads as “internal” reference for the encoding of the N+1. This method is shown in FIG. 51.
  • performing a so called “de-novo assembly” on a subset of the reads of class U so as to be able to align and encode all or a relevant sub-set of the reads belonging to said class according to the specified matching accuracy constraints or a new set of constraints.

If the read being coded can be mapped on the “internal” reference satisfying the specified set of matching accuracy constraints, the information necessary to reconstruct the read after compression is coded using syntax elements that can be of the following types:

• • 1. Start position of the matching portion on the internal reference in terms of read number in the internal reference (pos layer). This position can be encoded either as absolute or differential value with respect to the previously encoded read.
  • 2. Offset of the start position from the beginning of the corresponding read in the internal reference (pair layer). E.g. in case of constant read length the actual position is pos*length+pair.
  • 3. Possibly present mismatches coded as mismatch position (snpp layer) and type (snpt layer)
  • 4. Those parts of the reads (typically the edges identified by pair) that do not match with the internal reference (or do so, but with a number of mismatches above a defined threshold) are encoded in the indc layer. A padding operation can be performed to the edges of the part of the internal reference used in order to reduce the entropy of the mismatches encoded in the indc layer, as shown in FIG. 51. The most appropriate padding strategy can be chosen by the encoder according to the statistical properties of the genomic data being processed. Possible padding strategies include:
    • a. No padding
    • b. Constant padding pattern chosen according to its frequency in the currently encoded data.
    • c. Variable padding pattern according to the statistical properties of the current context defined in terms of the latest N encoded reads

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

• • 5. A flag that indicates if the read has been encoded using an internal self-generated, external or no-reference (rtype layer)
  • 6. Reads which are encoded verbatim (ureads).

FIG. 51 provides an example of such encoding procedure.

FIG. 56 shows an alternative encoding of unmapped reads on the internal reference where pos+pair syntax elements are replaced by a signed pos. In this case pos would express the distance —in terms of positions on the reference sequence —of the left most nucleotide position of read n with respect of the position of the left most nucleotide of read n−1.

This coding approach can be extended to support N start positions per read so that reads can be split over two or more reference positions. This can be particularly useful to encode reads generated by those sequencing technology (e.g. from Pacific Bioscience) producing very long reads (50K+bases) which usually present repeated patterns generated by loops in the sequencing methodology. The same approach can be used as well to encode chimeric sequence reads defined as reads that align to two distinct portions of the genome with little or no overlap.

The approach described above can be clearly applied beyond the simple class U and could be applied to any layer containing syntax elements related to reads positions (pos layers).

File Format: Selective Access to Regions of Genomic Data by Using the Master Index Table

In order to support selective access to specific regions of the aligned data, the data structure described in this document implements an indexing tool called Master Index Table (MIT). This is a multi-dimensional array containing the loci at which specific reads map on the used reference sequences. The values contained in the MIT are the mapping positions of the first read in each pos layer so that non-sequential access to each Access Unit is supported. The MIT contains one section per each class of data (P, N, M, I, U and HM) and per each reference sequence. The MIT is contained in the Genomic Dataset Header of the encoded data. FIG. 20 shows the structure of the Genomic Dataset Header, FIG. 21 shows a generic visual representation of the MIT and FIG. 22 shows an example of MIT for the class P of encoded reads.

The values contained in the MIT depicted in FIG. 22 are used to directly access the region of interest (and the corresponding AU) in the compressed domain.

For example, with reference to FIG. 22, if it is required to access the region comprised between position 150,000 and 250,000 on reference 2, a decoding application would skip to the second reference in the MIT and would look for the two values k1 and k2 so that k1<150,000 and k2>250,000. Where k1 and k2 are 2 indexes read from the MIT. In the example of FIG. 22 this would result in positions 3 and 4 of the second vector of the MIT. These returned values will then be used by the decoding application to fetch the positions of the appropriate data from the pos layer Local Index Table as described in the next section.

Together with pointers to the layer containing the data belonging to the four classes of genomic data described above, the MIT can be uses as an index of additional metadata and/or annotations added to the genomic data during its life cycle.

Local Index Table

Each data layer described above is prefixed with a data structure referred to as local header. The local header contains a unique identifier of the layer, a vector of Access Units counters per each reference sequence, a Local Index Table (LIT) and optionally some layer specific metadata. The LIT is a vector of pointers to the physical position of the data belonging to each AU in the layer payload. FIG. 23 depicts the generic layer header and payload where the LIT is used to access specific regions of the encoded data in a non-sequential way.

In the previous example, in order to access region 150,000 to 250,000 of reads aligned on the reference sequence no. 2, the decoding application retrieved positions 3 and 4 from the MIT. These values shall be used by the decoding process to access the 3^rd and 4^th elements of the corresponding section of the LIT. In the example shown in FIG. 24, the Total Access Units counters contained in the layer header are used to skip the LIT indexes related to AUs related to reference 1 (5 in the example). The indexes containing the physical positions of the requested AUs in the encoded stream are therefore calculated as:

position of the data blocks belonging to the requested AU=data blocks belonging to AUs of reference 1 to be skipped+position retrieved using the MIT, i.e.
First block position: 5+3=8
Last block position: 5+4=9

The blocks of data retrieved using the indexing mechanism called Local Index Table, are part of the Access Units requested.

FIG. 26 shows how the data blocks retrieved using the MIT and the LIT compose one or more Access Units.

Access Units The genomic data classified in data classes and structured in compressed or uncompressed layers are organized into different Access Units.

Genomic Access Units (AU) are defined as sections of genome data (in a compressed or uncompressed form) that reconstructs nucleotide sequences and/or the relevant metadata, and/or sequence of DNA/RNA (e.g. the virtual reference) and/or annotation data generated by a genome sequencing machine and/or a genomic processing device or analysis application. An example of Access Unit is provided in FIG. 26.

An Access Unit is a block of data that can be decoded either independently from other Access Units by using only globally available data (e.g. decoder configuration) or by using information contained in other Access Units.

Access Units are differentiated by:

• • type, characterizing the nature of the genomic data and data sets they carry and the way they can be accessed,
  • order, providing a unique order to Access Units belonging to the same type.

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

Hereafter follows a non-exhaustive list of definition of different types of genomic Access Units:

• 1) Access units of type 0 do not need to refer to any information coming from other Access Units to be accessed or decoded and accessed. The entire information carried by the data or data sets they contain can be independently read and processed by a decoding device or processing application.
• 2) Access units of type 1 contain data that refer to data carried by Access Units of type 0. Reading or decoding and processing the data contained in Access Units of type 1 requires having access to one or more Access Units of type 0. Access unit of type 1 encode genomic data related to sequence reads of “Class P”
• 3) Access Units of type 2 contain data that refer to data carried by Access Units of type 0. Reading or decoding and processing the data contained in Access Units of type 2 requires having access to one or more Access Units of type 0. Access unit of type 2 encode genomic data related to sequence reads of “Class N”
• 4) Access Units of type 3 contain data that refer to data carried by Access Units of type 0. Reading or decoding and processing the data contained in Access Units of type 3 requires having access to one or more Access Units of type 0. Access unit of type 3 encode genomic data related to sequence reads of “Class M”
• 5) Access Units of type 4 contain data that refer to data carried by Access Units of type 0. Reading or decoding and processing the data contained in Access Units of type 4 requires having access to one or more Access Units of type 0. Access unit of type 4 encode genomic data related to sequence reads of “Class I”
• 6) Access Units of type 5 contain reads that cannot be mapped on any available reference sequence (“Class U”) and are encoded used an internally constructed reference sequence. Access Units of type 5 contain data that refer to data carried by Access Units of type 0. Reading or decoding and processing the data contained in Access Units of type 5 requires having access to one or more Access Units of type 0.
• 7) Access Units of type 6 contain read pairs where one read can belong to any of the four classes P, N, M, I and the other cannot be mapped on any available reference sequence (“Class HM”). Access Units of type 6 contain data that refer to data carried by Access Units of type 0. Reading or decoding and processing the data contained in Access Units of type 6 requires having access to one or more Access Units of type 0.
• 8) Access Units of type 7 contain metadata (e.g. quality scores) and/or annotation data associated to the data or data sets contained in the access unit of type 1. Access Units of type 7 may be classified and labelled in different layers.
• 9) Access Units of type 8 contain data or data sets classified as annotation data. Access Units of type 8 may be classified and labelled in layers.
• 10) Access Units of additional types can extend the structure and mechanisms described here. As an example, but not as a limitation, the results of genomic variant calling, structural and functional analysis can be encoded in Access Units of new types. The data organization in Access Units described herein does not prevent any type of data to be encapsulated in Access Units being the mechanism completely transparent with respect to the nature of encoded data.

Access Units of type 0 are ordered (e.g. numbered), but they do not need to be stored and/or transmitted in an ordered manner (technical advantage: parallel processing/parallel streaming, multiplexing)

Access Units 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 ordered manner (technical advantage: parallel processing/parallel streaming).

FIG. 26 shows how Access Units are composed by a header and one or more layers of homogeneous data. Each layer can be composed by one or more blocks. Each block contains several packets and the packets are a structured sequence of the descriptors introduced above to represent e.g. reads positions, pairing information, reverse complement information, mismatches positions and types etc.

Each Access unit can have a different number of packets in each block, but within an Access Unit all blocks have the same number of packets.

Each data packet can be identified by the combination of 3 identifiers X Y Z where:

• • X identifies the access unit it belongs to
  • Y identifies the block it belongs to (i.e. the data type it encapsulates)
  • Z is an identifier expressing the packet order with respect to other packets in the same block

FIG. 28 shows an example of Access Units and packets labelling where AU T N is an access unit of type T with identifier N which may or may not imply a notion of order according to the Access Unit Type. Identifiers are used to uniquely associate Access Units of one type with those of other types required to completely decode the carried genomic data.

Access Units of any type can be further classified and labelled in different “categories” according to different sequencing processes. For example, but not as a limitation, classification and labelling can take place when

• • 1. sequencing the same organism at different times (Access Units contain genomic information with a “temporal” connotation),
  • 2. sequencing organic samples of different nature of the same organisms (e.g. skin, blood, hair for human samples). These are Access Units with “biological” connotation.

The access units of type 1, 2, 3, 4, 5 and 6 are built according to the result of a matching function applied on genome sequence fragments (a.k.a. reads) with respect to the reference sequence encoded in Access Units of type 0 they refer to.

For example access units (AUs) of type 1 (see FIG. 30) may contain the positions and the reverse complement flags of those reads which result in a perfect match (or maximum possible score corresponding to the selected matching function) when a matching function is applied to specific regions of the reference sequence encoded in AUs of type 0. Together with the data contained in AUs of type 0, such matching function information is sufficient to completely reconstruct all genome sequence reads represented by the data set carried by the access units of type 1.

With reference to the genomic data classification previously described in this document, the Access Units of type 1 described above would contain information related to genomic sequence reads of class P (perfect matches).

In case of variable reads length and paired reads the data contained in AUs of type 1 mentioned in the previous example, have to be integrated with the data representing the information about reads pairing and reads length in order to be able to completely reconstruct the genomic data including the reads pairs association. With respect to the data classification previously introduced in the present document, pair and rlen layers would be encoded in AU of type 1.

The matching functions applied with respect to access units of type 1 to classify the content of AU for the type 2, 3 and 4 can provide results such as:

• • 1. each sequence contained in the AU of type 1 perfectly matches sequences contained in the AU of type 0 in correspondence to the specified position;
  • 2. each sequence contained in the AU of type 2 perfectly matches a sequence contained in the AU of type 0 in correspondence to the specified position, except for the “N” symbols present (base not called by the sequencing device) in the sequence in the AU of type 2;
  • 3. each sequence contained in the AU of type 3 includes variants in the form of substituted symbols (variants) with respect to the sequence contained in the AU of type 0 in correspondence to the specified position;
  • 4. each sequence contained in the AU of type 4 includes variants in the form of substituted symbols (variants), insertions and/or deletions with respect to the sequence contained in the AU of type 0 in correspondence to the specified position.
  • 5. each sequence contained in the AU of type 5 do not map any sequence contained in the AU of type 0.
  • 6. each sequence pair contained in the AU of type 6 presents one sequence that can belong to any class P, N, M and I (points 1 to 4 above) while the other sequence does not map any sequence contained in the AU of type 0.

Access units of type 0 are ordered (e.g. numbered), but they do not need to be stored and/or transmitted in an ordered manner (technical advantage: parallel processing/parallel streaming, multiplexing)

Access units 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 ordered manner (technical advantage: parallel processing/parallel streaming).

Identifying Access Units by Using “Labels” Associated to Specific Genomic Regions

An additional mechanism is provided by the disclosed invention enabling user-defined selective access to data classes referring to specific genomic regions or sub-regions or aggregations of regions or sub-regions.

A “Label” is an identifier which is assigned to a specific genomic region or sub-region or aggregations of regions or sub-regions. Labels identify genomic regions by specifying: the reference sequence id (“Ref ids”), the index of the MIT corresponding to the desired region of the reference sequence, and the data classes. An example is provided in FIG. 52.

A single, a subset, or all data classes can be referenced by a Label, enabling selective access to only a sub-set of the data associated to a specific genomic region or sub-regions or aggregations of regions or sub-regions.

A Label list should be created by a Genomic Labels Generator (4917 FIG. 49), in a storage scenario, and/or in a streaming scenario to make available the available Labels to the analysis applications applying a selective access to the stored or streamed data.

A Label List might include the following elements:

• • the number of Labels
  • for each Label in the list:
    • the Label ID
    • the number of reference sequences concerned by the label
    • for each reference sequence
      • the reference identifier
      • the number of regions covered by the label,
      • for each region:
        • the class ID
        • the start position in the genomic range
        • the end position in the genomic range

The table below reports a pseudo-syntax for a generic “Label List”.

TABLE 4
Syntax of the generic “Label List” data format.

Syntax	Description

Label_list( ) {

	num_Labels	total number of labels in the list
	for (i=0; i<num_Labels;i++) {

	Label_id	label identifier
	num_ref	number of references concerned by the

current label

for (j = 0; j < num_ref; j++) {

	ref_id	current reference
	num_regions	number of different regions of this

reference identified by the label

for (k = 0; k < num_regions; k++) {

	class_id	type of class, start and end position of
	start_pos	this region
	end_pos

}

}

}

}

In case Genomic Data are compressed and streamed, one or more Access Units can be identified using a specific “Label” by means of a Block Header field (“Label ID”), which serves as an identifier for the “Label” in the “Label List” which the current block belongs to. Such field enables a dynamic mapping of blocks to “Labels”, typical for streaming scenarios.

In the Genomic File Format, the “start_pos” and “end_pos” fields can be replaced by the block numbers referring to all “blocks” belonging to a specific “Label”, as follows:

TABLE 5
Efficient implementation of the “Label List” Syntax
data format in the case of a compressed file.

Syntax	Data type	Description

num_Labels

Bitstring

number of labels in

the genomic dataset

for (i=0; i<num_Labels;i++) {

	Label_id	Bitstring	label identifier
	Label_length_in_blocks	Bitstring	number of data

	blocks identified by
	one label

for (j = 0; j < Label_length_in_blocks;j++) {

ref_id

Bitstring

reference id for this

block

	class_id	Bitstring	class id for this block
	block_num	Bitstring	block number in the

Master Index Table

}

}

The use of block numbers instead of “start_pos” and “end_pos” presents a relevant technical advantage because it enables a direct access to the “Master Index Table” (MIT), considering that the three dimensional vector consisting of: “ref_num”, “class_id” and “block_num” can be used as coordinates to directly address the MIT itself.

In storage scenarios, the “Label List” is created by a Genomic Labels Generator (4917) and sent to the genomic multiplexer (see also FIG. 49). The demultiplexer parses the Label List syntax and exposes the available Labels to the data access application, which according to the specific data access required selects the Access Units corresponding to the subset of “Labels”.

The possibility of using “Labels” to identify Access Units associated to specific genomic regions does not prevent using the indexing tools such as MIT and LIT without “Labels” to achieve random data access functionality. Generic random access can be achieved by specifying a three dimensional vector determining the MIT and LIT coordinates of interest (reference id, position range and classes) and ignoring the information carried by the Label List.

FIG. 51 shows how labels are used to aggregate and uniquely identify several genomic regions by using indexes contained in the MIT.

FIG. 59 shows how a device (592) implementing the labelling mechanism disclosed by this invention can enable concurrent access to several records of data (596) stored in a database (595). Selective protection of one or more regions identified by the same label is supported as well by means of a dedicated module (591) in charge of parsing the query (591) and dispatching the required metadata to the security module (594) in charge of enforcing access control. The labels decoder (593) is in charge of translating the label syntax into object identifiers that can be protected (and therefore access is controlled by the security module 594) or not.

Technical Effects

The technical effect of structuring genomic information in Access Units or Access Units identified by Labels as described here is that the genomic data:

1. can be selectively queried in order to access:

• • specific “categories” of data (e.g. with a specific temporal or biological connotation) without having to decompress the entire genomic data or data sets and/or the related metadata.
  • specific regions of the genome for all “categories”, a subset of “categories”, a single “category” (with or without the associated metadata) without the need to decompress other regions of the genome
  • specific genomic regions or sub-regions or aggregations of regions or sub-regions identified by user defined “Labels” by only parsing the “Label List” main header and accessing (i.e. retrieving and decompressing) only the corresponding Access Units
    2. can be incrementally updated with new data that can be available when:
  • new analysis is performed on the genomic data or data sets
  • new genomic data or data sets are generated by sequencing the same organisms (different biological samples, different biological sample of the same type, e.g. blood sample, but acquired at a different time, etc.)
    3. can be efficiently transcoded to a new data format in case of
  • new genomic data or data sets to be used as new reference (e.g. new reference genome carried by AU of type 0)
  • update of the encoding format specification
    4. can be protected with different levels of granularity in terms of both access control (e.g. encryption) and permissions enforcement. For example these scenarios are enabled:
  • the same access control rule and encryption keys can be applied to all the genomic regions or sub-regions identified by one label (see FIG. 54 for an example);
  • different access control rules and different encryption keys can be used to protect each single region or sub-regions aggregated under the same label (see FIG. 55 for an example).

With respect to prior art solutions such as SAM/BAM, the described technical features address the issues of requiring data filtering to happen at the application level when the entire data has been retrieved and decompressed from the encoded format.

Hereafter follows examples of application scenario where the association of access unit structure, file format and Labelling mechanism becomes instrumental for a technological advantage.

Selective Access

In particular the disclosed data structure based on Access Units of different types possibly including user defined “Labels” enables to:

• • extract only the read information (data or data sets) of the whole sequencing of all “categories” or a subset (i.e. one or more layers) or a single “category” without having to decompress also the associated metadata information (limitation of current state of the art: SAM/BAM that cannot even support distinction between different categories or layers);
  • extract all the reads aligned on specific regions of the assumed reference sequence for all categories, subsets of the categories, a single category (with or without the associated metadata) without the need of decompressing also other regions of the genome (limitation of current state of the art: SAM/BAM);
  • extract all the reads belonging to a single, a subset or all data “classes” aligned on specific genomic regions or sub-regions or aggregations of regions or sub-regions identified by user specified “Labels” for all categories, subsets of the categories, a single category (with or without the associated metadata) without the need of decompressing also other data associated to other regions of the genome (limitation of current state of the art: SAM/BAM).

FIG. 39 shows how the access to the genomic information mapped on the second segment of the reference sequence (AU 0-2) with mismatches only requires the decoding of AUs 0-2, 1-2 and 3-2 only. This is an example of selective access according to both a criteria related to a mapping region (i.e. position on the reference sequence) and a criteria related to the matching function applied to the encoded sequence reads with respect to the reference sequence (e.g. mismatches only in this example).

A further technical advantage is that the querying on the data is much more efficient in terms of data accessibility and execution speed because it can be based on accessing and decoding only selected “categories”, specific regions of longer genomic sequences and only specific layers for access units of type 1, 2, 3, 4 that match the criteria of the applied queries and any combination thereof.

The organization of access units of type 1, 2, 3, 4 into layers allow for efficient extraction of nucleotides sequences

• • with specific variations (e.g. mismatches, insertions, deletions) with respect to one or more reference genomes;
  • that do not map to any of the considered reference genomes;
  • that perfectly map on one or more reference genomes;
  • that map with one or more accuracy levels.

FIG. 52 shows how the access to the genomic information associated only to specific genomic regions or sub-regions or aggregations of regions or sub-regions associated to user defined “Labels”. The syntax of a label is based on a three coordinates system where each region or sub-region associated to a label can be uniquely identified by:

• • 1. reference ID,
  • 2. data type (class)
  • 3. block number in the MIT (corresponding to a genomic region).

These three coordinates can be used to identify

• • the MIT location containing the genomic position of the region on the corresponding reference and
  • the LIT location containing the physical location of the data representing the respective genomic region or sub-region

Like in the case of accessing data related to a specified genomic region, a further technical advantage is that the querying on the data results to be much more efficient in terms of data accessibility and execution speed because it can be based on accessing and decoding only selected “categories”, of the labelled specific regions and only specific layers for access units of type 1, 2, 3, 4 that corresponds to the “Labels” of the applied queries and any combination thereof.

Another technical advantage of this labelling mechanism is the possibility of efficiently retrieving encoded genomic information that has been scattered among several Access Units due to its characteristics such as position on the reference genome, type of mismatches with respect to the reference (524).

Filtering genomic data according to the characteristics of the mapped reads (e.g. perfectly matching, substitutions only, etc.) today can take hours when using the traditional formats such as BAM and CRAM. This is due to the fact that the data are sparse within the compressed format and require decompression and filtering using pipelines of commands. The present invention describes a data structure that enables data filtering in a matter of seconds. Memory usage can be as well reduced by a factor that is proportional with the file size (from 10× to 100×) since the present invention does not require the decoding (i.e. memory allocation) of the entire file.

Selective Access to Specific Genomic Regions Identified by User Specified “Labels” in “Storage” and “Streaming” Scenarios.

For example let's suppose that sequencing data is compressed and selective access to “GeneXY” and “GeneWZ” is required. The two genomic regions corresponding to “GeneXY” and “GeneWZ” in the compressed file format or in the compressed stream must be labelled. Depending if a compressed data file is generated for storage or a compressed data stream is generated for streaming, two methods are used.

In the case of a compressed data file, the multiplexer creates a “Label List” which includes two Labels with: “Label_ID”=GeneXY and “Label_ID”=GeneWZ. The Label parameter “Label_lenght_in_blocks” and for each block the parameters: “ref_num”, “class_ID”, “block_num” are determined by the multiplexer based on the position on the reference of the “GeneXY” and “GeneWZ” regions and the class of data for which the selective access is desired. The complete syntax is reported in Table 5.

In the case of a compressed stream, the multiplexer creates a “Label List” which includes two Labels with: “Label_ID”=GeneXY and “Label_ID”=GeneWZ. The Label parameters “ref ID”, “class_ID”, “start_pos” and “end_pos” are determined by the multiplexer based on the position on the reference of the “GeneXY” and “GeneWZ” regions and the class of data for which the selective access is desired. The complete syntax is reported in Table 4.

The method used in the case of a compressed stream is generic and could be used also in the case of a compressed file for storage, but the corresponding implementation would result less efficient because the use of block numbers, as described in the case of compressed file, enables a direct access to the “Master Index Table” (MIT).

In both cases mentioned above (streaming and storage), the mechanism of retrieval of the genomic data identified by the labels follows is the same.

When parsing a label a decoding device will:

• • 1. Identify the reference sequence from the first element of the label
  • 2. Identify the class of data from the second element of the label
  • 3. Identify the block of the MIT (corresponding to a genomic region) from the third element of the label
  • 4. The two coordinates parsed in 1 and 2 enable the decoder to identify the required Genomic Streams (484),
  • 5. Each Genomic Stream starts with a header containing a LIT (525) containing pointers to the descriptors encoding data mapped to each genomic region. The third coordinate parsed in 3 is used to access the correct pointer in the LIT of each Genomic Stream.
  • 6. The decoder can efficiently retrieve all the descriptors to decode the genomic data identified by the decoded Genomic Label even if they are scattered among different Access Units (524).

Incremental Update

The access units of type 7 and 8 allow for easy insertion of annotations without the need of depacketizing/decoding/decompressing the whole file thereby adding to the efficient handling of the file which is a limitation of prior art approaches. Existing compression solutions may have to access and process a large amount of compressed data before the desired genomic data can be accessed. This will cause inefficient RAM bandwidth utilization and more power consumption also in hardware implementations. Power consumption and memory access issues may be alleviated by using the approach based on Access Units described here.

The data indexing mechanism described in the Master Index Table (see FIG. 21) together with the utilization of Access Unites and the possibility of identifying Access Units with user-defined “Labels” associated to specific genomic regions or sub-regions or aggregations of regions or sub-regions enables incremental update of the encoded content as described below. This mechanism is shown with an example in FIG. 53.

Insertion of Additional Data

New genomic information can be periodically added to existing genomic data for several reasons. For example when:

• • An organism is sequenced at different moments in time;
  • Several different samples of the same individual are sequenced at the same time;
  • New data generated by a sequencing process (streaming).

In the above mentioned situations, structuring data using the Access Units described here and the data structure described in the file format section enables the incremental integration of the newly generated data without the need to re-encode the existing data. The incremental update process can be implemented as follows:

• • 1. The newly generated AUs can simply be concatenated in the file with the pre-existing AUs and
  • 2. the indexing of the newly generated data or data sets are included in the Master Index Table described in the file format section of this document. One index shall position the newly generated AU on the existing reference sequence, other indexes consist in pointers of the newly generated AUs in the physical file to enable direct and selective access to them.
  • 3. The existing and/or newly generated AU can be identified with user defined “Labels” corresponding to specific genomic regions or sub-regions or aggregations of regions or sub-regions and a “Label List” can be included or updated.

This mechanism is illustrated in FIG. 40 where pre-existing data encoded in 3 AUs of type 1 and 4 AUs per each type from 2 to 4 are updated with 3 AUs per type with encoding data coming for example from a new sequence run for the same individual.

The mechanism of creating or updating “Labels” and the “Label List” are illustrated in FIG. 52 and FIG. 53.

In the specific use case of streaming genomic data and data sets in compressed form, the incremental update of a pre-existing data set may be useful when analyzing data as soon as they are generated by a sequencing machine and before the actual sequencing is completed. An encoding engine (compressor) can assemble several AUs in parallel by “clustering” sequence reads that map on the same region of the selected reference sequence. Once the first AU contains a number of reads above a pre-configured threshold/parameter, the AU is ready to be sent to the analysis application. Together with the newly encoded Access Unit, the encoding engine (the compressor) shall make sure that all Access Units the new AU depends on have already been sent to the receiving end or is sent together with it. For example an AU of type 3 will require the appropriate AU of type 0 and type 1 to be present at the receiving end in order to be properly decoded.

By means of the described mechanism, a receiving variant calling application would be able to start calling variants on the AU received before the sequencing process has been completed at the transmitting side. A schematic of this process is depicted in FIG. 41.

New Analysis of Results.

During the genome processing life cycle several iterations of genome analysis can be applied on the same data (e.g. different variant calling using different processing algorithm). The use of AUs as defined in this document and the data structure described in the file format section of this document enable incremental update of existing compressed data with the results of new analysis. For example, new analysis performed on existing compressed data can produce new data in these cases:

• • 1. A new analysis can modify existing results already associated with the encoded data. This use case is depicted in FIG. 42 and it is implemented by moving entirely or partially the content of one Access Unit from one type to another. In case new AUs need to be created (due to a pre-defined maximum size per AU), the related indexes in the Master Index Table must be created and the related vector is sorted when needed.
  • 2. New data are produced from new analysis and have to be associated to existing encoded data. In this case new AUs of type 7 can be produced and concatenated with the existing vector of AUs of the same type. This and the related update of the Master Index Table are depicted in FIG. 43

The use cases described above and depicted in FIG. 42 and FIG. 43 are enabled by:

• • 1. The possibility to have direct access only to data with poor mapping quality (e.g. AUs of type 4);
  • 2. The possibility to remap reads to a new genomic region by simply creating a new Access Unit possibly belonging to a new type (e.g. reads included in a Type 4 AU can be remapped to a new region with less (type 2-3) mismatches and included in a newly created AU);
  • 3. The possibility to create AU of type 8 (433) containing only the newly created analysis results and/or related annotations. In this case the newly created AUs only require to contain “pointers” to the existing AUs to which they refer to.
  • 4. The possibility of performing in a single run new analysis on several genomic regions and sub-regions identified by the same Label without the need to repeat the analysis on each single genomic region or sub-region. Labels as described in this document would enable users to manipulate non-contiguous genomic segments as if they were a single genomic sequence.
  • 5. The possibility of updating with new analysis results several genomic regions or sub regions identified by a single Label. The new results (usually expressed in the form of metadata) would be linked to the label identifying the aggregation of potentially several genomic regions and sub regions without the need of creating several links from the results to each genomic region or sub region.

Transcoding

Compressed genomic data can require transcoding, for example, in the following situations:

• • Publication of new reference sequences;
  • Use of a different mapping algorithm (re-mapping).

When genomic data is mapped on an existing public reference genome, the publication of a new version of said reference sequence or the desire to map the data using a different processing algorithm, today requires a process of re-mapping. When remapping compressed data using prior art file formats such as SAM or CRAM the entire compressed data has to be decompressed into its “raw” form to be mapped again with reference to the newly available reference sequence or using a different mapping algorithm. This is true even if the newly published reference is only slightly different from the previous or the different mapping algorithm used produces a mapping that is very close (or identical) to the previous mapping.

The advantage of transcoding genomic data structured using Access Units described here is that:

• 1. Mapping versus a new reference genome only requires re-encoding (decompressing and compressing) the data of AUs that map on the genome regions that have changes. Additionally the user may select those compressed reads that for any reason might need to be re-mapped even if they originally do not map on the changed region (this may happen if the user believes that the previous mapping is of poor quality). This use case is depicted in FIG. 44.
• 2. In case the newly published reference genome differs from the previous only in terms of entire regions shifted to different genomic locations (“loci”), the transcoding operation results particularly simple and efficient. In fact in order to move all the reads mapped to the “shifted” region it is sufficient to change only the value of the absolute position contained in the related (set of) AU(s) header. Each AU header contain the absolute position the first read contained in the AU is mapped to on the reference sequence, while all other reads positions are encoded differentially with respect to the first. Therefore, by simply updating the value of the absolute position of the first read, all the reads in the AU are moved accordingly. This mechanism cannot be implemented by state of the art approaches such as CRAM and BAM because genome data positions are encoded in the compressed payload, thus requiring complete decompression and re-compression of all genome data sets.
• 3. When a different mapping algorithm is used, it is possible to apply it only on a portion of compressed reads that was deemed mapped with poor quality. For example it can be appropriate to apply the new mapping algorithm only on reads which did not perfectly match on the reference genome. With existing formats today it is not possible (or it's only partially possible with some limitations) to extract reads according to their mapping quality (i.e. presence and number of mismatches). If new mapping results are returned by the new mapping tools the related reads can be transcoded from one AU from another of the same type (FIG. 46) or from one AU of one type to an AU of another type (FIG. 45).

Moreover, prior art compression solutions may have to access and process a large amount of compressed data before the desired genomic data can be accessed. This will cause inefficient RAM bandwidth utilization and more power consumption and in hardware implementations. Power consumption and memory access issues may be alleviated by using the approach based on Access Units described here.

A further advantage of the adoption of the genomic access units described here is the facilitation of parallel processing and suitability for hardware implementations. Current solutions such as SAM/BAM and CRAM are conceived for single-threaded software implementation.

Selective Protection

The approach based on Access Units organized in several types an layers as described in this document enables the implementation of content protection mechanisms otherwise not possible with state of the art monolithic solutions.

A person skilled in the art knows that the majority of genomic information related to an organism's genetic profile relies in the differences (variants) with respect to a known sequence (e.g. a reference genome or a population of genomes). An individual genetic profile to be protected from unauthorized access will therefore be encoded in Access Units of type 3 and 4 as described in this document. The implementation of controlled access to the most sensible genomic information produced by a sequencing and analysis process can therefore be realized by encrypting only the payload of AUs of type 3 and 4 (see FIG. 47 for an example). This will generate significant savings in terms of both processing power and bandwidth since the resources consuming encryption process shall be applied on a subset of data only.

Selective Protection of Specific Genomic Regions Identified by “Labels”

The labelling mechanism enables different mechanisms of data protection and access control. For example FIG. 54 shows how one protection mechanism (e.g. encryption) and one access control rule (AC) can be applied to several genomic regions identified by the same label. In a more sophisticated scenario, data protection can be implemented by applying a different access control rule and a different protection mechanism (encryption) to each region identified by a label. This is shown in FIG. 55.

Additionally, selective encryption of genomic regions or sub-regions or aggregations of regions or sub-regions identified by different “Labels” can be easily implemented by applying encryption only to compressed data corresponding to a “Label” for both file and streamed scenarios. For instance two genomic regions labelled as “GeneXY” and “GeneWZ” like in the example of section “Selective Access to Specific Genomic Regions identified by User Specified “Labels” in “storage” and “streaming” scenarios” can be differentiated by only encrypting data labelled by “GeneXY” and leaving in clear the compressed data labelled as “GeneWZ”. Encryption rules can be carried by the metadata fields (in both storage and streaming scenarios) and associated to each element of the “Label List”

Transport of Genomic Access Units

Genomic Data Multiplex

Genomic Access Units can be transported over a communication network within a Genomic Data Multiplex. A Genomic Data Multiplex is defined as a sequence of packetized genomic data and metadata represented according to the data classification disclosed as part of this invention, transmitted in network environments where errors, such as packet losses, may occur.

The Genomic Data Multiplex is conceived to ease and render more efficient the transport of genomic coded data over different environments (typically network environments) and has the following advantages not present in state of the art solutions:

• 1. it enables encapsulation of either a stream or a sequence of genomic data (described below) or Genomic File Format generated by an encoding tool into one or more Genomic Data Multiplex, in order to carry it over a network environment, and then recover a valid and identical stream or file format in order to render the transmission and access to information more efficient
• 2. It enables selective retrieval of encoded genomic data from the encapsulated Genomic Data Streams, for decoding and presentation.
• 3. It enables multiplexing several Genomic Datasets into a single container of information for transport and it enables de-multiplexing a subset of the carried information into a new Genomic Data Multiplex.
• 4. It enables the multiplexing of data and metadata produced by different sources (with the consequent separate access) and/or sequencing/analysis processes and transmit the resulting Genomic Data Multiplex over a network environment.
• 5. It supports identification of errors such as packet losses.
• 6. It supports proper reorder data which may arrive out of order due to network delays, therefore rendering more efficient the transmission of genomic data when compared with the state of the art solutions

An Example of Genomic Data Multiplexing is Shown in FIG. 49.

Genomic Dataset

In the context of the present invention a Genomic Dataset is defined as a structured set of Genomic Data including, for example, genomic data of a living organism, one or more sequences and metadata generated by several steps of genomic data processing, or the result of the genomic sequencing of a living organism. One Genomic Data Multiplex may include multiple Genomic Datasets (as in a multi-channel scenario) where each dataset refers to a different organism. The multiplexing mechanism of the several datasets into a single Genomic Data Multiplex is governed by information contained in data structures called Genomic Datasets List (GDL), Genomic Dataset Mapping Tables List (GDMTL) and Genomic Dataset Mapping Table (GDMT).

Genomic Dataset List

A Genomic Dataset List (GDL) is defined as a data structure listing all Genomic Datasets available in a Genomic Data Multiplex. Each of the listed Genomic Datasets is identified by a unique value called Genomic Dataset ID (GID).

Each Genomic Dataset listed in the GDL is associated to:

• • one Genomic Data Stream carrying one Genomic Dataset Mapping Table (GDMT) and identified by a specific value of Stream ID (genomic_dataset_map_SID);
  • one Genomic Data Stream carrying one Reference ID Mapping Table (RIDMT) and identified by a specific value of Stream ID (reference_id_map_SID).

The GDL is sent as payload of a single Transport Packet at the beginning of a Genomic Data Stream transmission; it can then be periodically re-transmitted in order to enable random access to the Stream.

The syntax of the GDL data structure is provided in the table below with an indication of the data type associated to each syntax element.

Syntax	Data type

genomic_dataset_list( ) {

	list_length	bitstring
	multiplex_id	bitstring
	version_number	bitstring
	applicable_section_flag	bit
	list_ID	bitstring
	for (i = 0; i < N; i++) {	N = number of Genomic

	Datasets in this
	Genomic Multiplex

	genomic_dataset_ID	bitstring
	genomic_dataset_map_SID	bitstring
	reference_id_map_SID	bitstring
	genomic_Label_list_SID	bitstring

	}
	Checksum	bitstring

}

The syntax elements composing the GDL described above have the following meaning and function.

section_length	bitstring field, specifying the number of bytes composing the
	section, starting immediately following the section_length field, and
	including the CRC.
multiplex_id	bitstring field which serves as a label to identify this multiplexed
	stream from any other multiplex within a network.
version_number	bitstring field indicating the version number of the whole Genomic
	Dataset List Section. The version number shall be incremented by 1
	whenever the definition of the Genomic Dataset Mapping Table
	changes. When the applicable_section_flag is set to ‘1’, then the
	version_number shall be that of the currently applicable Genomic
	Dataset List. When the applicable_section_flag is set to ‘0’, then the
	version_number shall be that of the next applicable Genomic
	Dataset List.
applicable_section_flag	A 1 bit indicator, which when set to ‘1’ indicates that the Genomic
	Dataset Mapping Table sent is currently applicable. When the bit is
	set to ‘0’, it indicates that the table sent is not yet applicable and
	shall be the next table to become valid.
list_ID	This is a bitstring field identifying the current genomic dataset list.
genomic_dataset_ID	genomic_dataset_ID is a bitstring field which specifies the genomic
	dataset to which the genomic_dataset_map_SID is applicable. This
	field shall not take any single value more than once within one
	version of the Genomic Dataset Mapping Table.
genomic_dataset_map_SID	genomic_dataset_map_SID is a bitstring field identifying the
	Genomic Data Stream carrying the Genomic Dataset Mapping Table
	(GDMT) associated to this Genomic Dataset. No
	genomic_dataset_ID shall have more than one
	genomic_dataset_map_SID associated. The value of the
	genomic_dataset_map_SID is defined by the user.
reference_id_map_SID	reference_id_map_SID is a bitstring field identifying the Genomic
	Data Stream carrying the Reference ID Mapping Table (RIDMT)
	associated to this Genomic Dataset. No genomic_dataset_ID shall
	have more than one reference_id_map_SID associated. The value of
	the reference_id_map_SID is defined by the user.
genomic_Label_list_SID	genomic_Label_list_SID is a bitstring field identifying the Genomic
	Data Stream carrying the Genomic Label List (GLL) associated to this
	Genomic Dataset. No genomic_dataset_ID shall have more than
	one genomic_Label_list_SID associated. The value of the
	genomic_Label_list_SID is defined by the user.
Chacksum	This is a bitstring field that contains an integrity check value for the
	entire GDL. One typical algorithm used for this purpose function is
	the CRC32 algorithm producing a 32 bit value other algorithms
	include the hashing functions MD5 and SHA-256.

Genomic Dataset Mapping Table

The Genomic Dataset Mapping Table (GDMT) is produced and transmitted at the beginning of a streaming process (and possibly periodically re-transmitted, updated or identical in order to enable the update of correspondence points and the relevant dependencies in the streamed data). The GDMT is carried by a single Packet following the Genomic Dataset List and lists the SIDs identifying the Genomic Data Streams composing one Genomic Dataset. The GDMT is the complete collection of all identifiers of Genomic Data Streams (e.g., the genomic sequence, reference genome, metadata, etc) composing one Genomic Dataset carried by a Genomic Multiplex. A genomic dataset mapping table is instrumental in enabling random access to genomic sequences by providing the identifier of the stream of genomic data associated to each genomic dataset.

The syntax of the GDMT data structure is provided in the table below with an indication of the data type associated to each syntax element.

genomic_dataset_mapping_table( ) {

	table_length	bitstring
	genomic_dataset_ID	bitstring
	version_number	bitstring
	applicable_section_flag	bit
	mapping_table_ID	bitstring
	genomic_dataset_ef_length	bitstring
	for (i=0; i<N; i++) {	N = number of extension fields

	associated to this
	Genomic Dataset

extension_field( )

data structure

	}
	for (i = 0;i < M ; i++) {	M = number of Genomic Data

	Streams associated to
	this specific Dataset

	data_type	bitstring
	genomic_data_SID	bitstring
	gd_component_ef_length	bitstring
	for (I = 0; I < K; i++) {	K = number of extension fields

	associated to each
	Genomic Data Stream

extension_field ( )

data structure

}

	}
	Chaecksum	bitstring

}

The syntax elements composing the GDMT described above have the following meaning and function.

version_number,	These elements have the same meaning as for the GDL
applicable_section_flag
table_length,	bitstring field specifying the number of bytes composing the table,
	starting after the table_length field, and including the Checksum field.
genomic_dataset_ID	bitstring field identifying a Genomic Dataset
mapping_table_ID	bitstring bit field identifying the current Genomic Dataset Mapping
	Table
genomic_dataset_ef_length	bitstring field specifying the number of bytes of the optional
	extension_field associated with this Genomic Dataset
data_type	bitstring field specifying the type of genomic data carried by the
	packets identified by the genomic_data_SID.
genomic_data_SID	bitstring bit field specifying the Stream ID of the packets carrying the
	encoded genomic data associated with one component of this
	Genomic Dataset (e.g. read p positions, read p pairing information
	etc. as defined in this invention)
gd_component_ef_length	bitstring field specifying the number of bytes of the optional
	extension_field associated with the genomic Stream identified by
	genomic_data_SID.
Checksum	This is a bitstring field that contains an integrity check value for the
	entire GDMT. One typical algorithm used for this purpose function is
	the CRC32 algorithm producing a 32 bit value or hashing functions
	such as MD5 and SHA-256.

extension_fields are optional descriptors that might be used to further describe either a Genomic Dataset or one Genomic Dataset component.

The data_type field can have the following values

	data_type	Description
	0	Dataset Header
	1	Layer Header
	2 to 15	User-defined extensions
	16 to N	16 + Descriptors_Layer_ID

Genomic Datasets Mapping Tables List

This structure carries information about all the datasets mapping tables related to a Genomic Datasets Multiplex.

Syntax	Description

Datasets_mapping_tables_list{

Multiplex_id

Datasets Multiplex ID, as in

Datasets Multiplex Header.

for (i=0; i<gd_number;i++) {

Note: gd_number as in Datasets

Multiplex Header.

dataset_mapping_table_SID

Stream ID of Dataset Mapping

Table of i-th Dataset.

}

}

Reference ID Mapping Table

The Reference ID Mapping Table (RIDMT) is produced and transmitted at the beginning of a streaming process. The RIDMT is carried by a single Packet following the Genomic Dataset List. The RIDMT specifies a mapping between the numeric identifiers of reference sequences (REFID) contained in the Block header of an access unit and the (typically literal) reference identifiers contained in the Genomic Dataset Header specified in Table 2.

The RIDMT can be periodically re-transmitted in order to:

• • enable the update of correspondence points and the relevant dependencies in the streamed data,
  • support the integration of new reference sequences added to the pre-existing ones (e.g. synthetic references created by de-novo assembly processes)

The syntax of the RIDMT data structure is provided in the table below with an indication of the data type associated to each syntax element.

Syntax	Data type

reference_id_mapping_table( ) {

	table_length	bitstring
	genomic_dataset_ID	bitstring
	version_number	bitstring
	applicable_section_flag	bit

reference_id_mapping_table_ ID

bitstring

for (i = 0; i < N; i++) {

N = number of reference

	sequences associated with
	the Genomic Dataset identified
	by genomic_dataset_ID

	ref_string_length	bitstring
	for

(i=0;i<ref_string_length;i++){

ref_string[i]

byte

	}
	REFID	bitstring

	}
	Checksum	bitstring (e.g. CRC-32 or

	MD5 hash)
}

The syntax elements composing the RIDMT described above have the following meaning and function.

table_length, genomic_dataset_ID,	These elements have the same meaning as for the
version_number, applicable_section_flag	GDMT
reference_id_mapping_table_ID	bitstring field identifying the current Reference ID
	Mapping Table
ref_string_length	bitstring field specifying the number of characters
	(bytes) composing ref_string, excluding the end of
	string (‘\0’) character.
ref_string[i]	byte field encoding each character of the string
	representation of a reference sequence (e.g. “chr1”
	for chromosome 1). The end of string (‘\0’) character
	is not necessary, as it is implicitly inferred from the
	ref_string_length field
REFID	This is a bitstring field uniquely identifying a reference
	sequence. This is encoded in the data Block header as
	REFID field.
Checksum	This is a bitstring field that contains an integrity check
	value for the entire RIDMT. One typical algorithm
	used for this purpose function is the CRC32 algorithm
	producing a 32 bit value or any hash function
	producing longer strings of bits.

Genomic Label List

As described above, a label is an identifier which is assigned to a specific genomic regions or sub-regions or aggregations of regions or sub-regions.

Labels identify genomic regions by specifying the reference sequence id, the position range with respect to the reference sequence and the data classes that they identify.

For such purpose, the Genomic Label List (GLL) is created during the packetization process by the multiplexer and transmitted.

The depacketizer of the demultiplexer parses the GLL syntax and exposes the available “Labels” to the data access application, which has the possibility to select and access the desired sub-set of data.

The GLL is (optionally) produced and transmittedat the beginning of a stream and typically transmitted periodically in order to enable multiple synchronization points (4811), and provides the list of “Labels” associated to the Multiplex and Dataset identified by the multiplex_id and dataset_id fields.

The syntax of the GLL data structure is provided in the table below with an indication of the data type associated to each syntax element.

TABLE 6
Complete syntax of “Label List” data format
for the streamed compressed data scenario.

Syntax	Description

genomic_label_list( ) {

	table_length
	multiplex_id
	dataset_id
	num_labels	total number of labels in the list
	for (i=0; i<num_labels;i++) {

	Label_id	label identifier
	num_ref	number of references concerned by the

current label

for (j = 0; j < num_ref; j++) {

	ref_id	current reference
	num_regions	number of different regions of this

reference identified by the label

for (k = 0; k < num_regions; k++) {

	class_id	type of class, start and end position of
	start_pos	this region
	end_pos

}

}

}

Checksum	e.g. CRC-32 or MD5 hash
}

The syntax elements composing the GLL described above have the following meaning and function.

TABLE 7
Description of syntax elements of Table 6.

table_length	Bitstring field specifying the number of bytes composing the list,
	starting after the table_length field, and including the Checksum field
multiplex_ID	Byte which serves as a label to identify the Genomic Multiplex from
	any other multiplex within a network
dataset_ID	Byte which serves as a label to identify the Genomic Dataset from any
	other dataset within the multiplex identified by multiplex_id
num_Labels	Bitstring representing the total number of Labels in this GLL
Label_id	Bitstring identifying the i-th Label
num_ref	Bitstring identifying the number of references concerned by the
	current label
ref_id	Bitstring identifying the j-th reference sequence the i-th Label refers
	to
num_regions	Bistring identifying the number of regions conveyed by the i-th Label
class_id	Bitstring identifying the class of the k-th region in the j-th reference in
	the i-th Label
start_pos	Bitstring indicating the position in the j-th reference sequence of the
	first read of the k-th region in the i-th Label
end_pos	indicating the position in the j-th reference sequence of the last read
	of the k-th region in the i-th Label
Checksum	Bitstring field that contains an integrity check value for the entire GLL.
	One typical algorithm used for this purpose function is the CRC32
	algorithm producing a 32 bit value or hashing functions producing
	longer strings of bits (e.g. MD5, SHA-256).

Genomic Data Stream

A Genomic Data Multiplex contains one or several Genomic Data Streams where each stream can transport

• • data structures containing transport information (e.g. Genomic Dataset List, Genomic Dataset Mapping Table etc.)
  • data belonging to one of the Genomic Data Layers described in this invention.
  • Metadata related to the genomic data
  • Any other data

A Genomic Data Stream containing genomic data is essentially a packetized version of a Genomic Data Layer where each packet is prepended with a header describing the packet content and how it is related to other elements of the Multiplex.

The Genomic Data Stream format described in this document and the File Format described in this document are mutually convertible. Whereas a full file format can be reconstructed in full only after all data have been received, in case of streaming a decoding tool can reconstruct and access, and start processing the partial data at any time.

A Genomic Data Stream is composed by several Genomic Data Blocks each containing one or more Genomic Data Packets. Genomic Data Blocks (GDBs) are containers of genomic information composing one genomic AU. GDB can be split into several Genomic Data Packets, according to the communication channel requirements.

Genomic access units are composed by one or more Genomic Data Blocks belonging to different Genomic Data Streams.

Genomic Data Packets (GDPs) are transmission units composing one GDB. Packet size is typically set according to the communication channel requirements.

FIG. 27 shows the relationship among Genomic Multiplex, Streams, Access Units, Blocks and Packets when encoding data belonging to the P class as defined in this invention. In this example three Genomic Streams encapsulate information on position, pairing and reverse complement of sequence reads.

Genomic Data Blocks are composed by a header, a payload of compressed data and padding information.

The table below provides an example of implementation of a GDB header with a description of each field and a typical data type.

TABLE 8
Description of Genomic Data Block syntax elements.

Data type	Description	Data type
Block Start Code	Reserved value used to unambiguously identify the beginning	bitstring
Prefix (BSCP)	of a Genomic Data Block.
Block Header	Block Header as defined in this document	bitstring
POS Flag (PSF)	If the POS Flag is set, the block contains the 40 bit POS field at	bit
	the end of the block header and before the optional fields.
Padding Flag (PDF)	If the Padding Flag is set, the block contains additional padding	bit
	bytes after the payload which are not part of the payload.
Block size (BS)	Number of bytes composing the block, including this header	bitstring
	and payload, and excluding padding (total block size will be
	BS + padding size).
Access Unit ID (AUID)	Unambiguous ID, linearly increasing (not necessarily by 1, even	bitstring
	though recommended). Needed to implement proper random
	access, as described in the Master Index Table defined in this
	invention.
Label ID (LID)	Unambiguous ID, linearly increasing by 1, identifying the	bitstring
	genomic region/classes (Label) this block belongs to. It
	corresponds to the i-th index in the main for loop in the
	Genomic Label List described above.
(Optional) Reference	Unambiguous ID, identifying the reference sequence the AU	bitstring
ID (REFID)	containing this block refers to. This is needed, along with POS
	field, to have proper random access, as described in the
	Master Index Table.
(Optional) POS (POS)	Present if PSF is 1. Position on the reference sequence of the	bitstring
	first read in the block.
(Extra optional fields)	Additional optional fields, presence signaled by BS.	bytestring
(Optional) Padding	(Optional, presence signaled by PDF) Fixed bitstring value that	bitstring
	can be inserted in order to meet the channel requirements. If
	present, the first byte indicates how many bytes compose the
	padding. It is discarded by the decoder.

The use of AUID, POS and BS enables the decoder to reconstruct the data indexing mechanisms referred to as Master Index Table (MIT) and Local Index Table (LIT) in this invention. In a data streaming scenario the use of AUID and BS enables the receiving end to dynamically re-create a LIT locally, without the need to send extra-data. The use of AUID, BS and POS will enable to recreate a MIT locally without the need to send additional data.

This has the technical advantage to

• • reduce the encoding overhead which might be large if the entire LIT is transmitted;
  • avoid the need of a complete mapping between genomic positions and Access Units which is not normally available in a streaming scenario

A Genomic Data Block can be split into one or more Genomic Data Packets, depending on network layer constraints such as maximum packet size, packet loss rate, etc. A Genomic Data Packet is composed by a header and a payload of encoded or encrypted genomic data as described in the table below.

TABLE 9
Description of Genomic Data Packet syntax elements.

Data type	Description	Data size
Stream ID (SID)	Unambiguously identifies data type carried by this	bitstring
	packet. A Genomic Dataset Mapping Table is needed at
	the beginning of the stream in order to map Stream IDs
	to data types. Used also for updating correspondence
	points and relevant dependencies.
Access Unit Marker Bit	Set for the last packet of the access unit. Allows to	bit
(MB)	identify the last packet of an AU.
Packet Counter	Counter associated to each Stream ID linearly increasing	bitstring
Number (SN)	by 1. Needed to identify gaps/packet losses. Wrap
	around at 255.
Packet Size (PS)	Number of bytes composing the packet, including	bitstring
	header, optional fields and payload.
Extension Flag (EF)	Set if extension fields are present.	bit
Extension Fields	Optional fields, presence signaled by PS.	bytestring
Payload	Block data (entire block or fragment)	bytestring

The Genomic Multiplex can be properly decoded only when at least one Genomic Dataset List, one Genomic Dataset Mapping Table and one Reference ID Mapping Table have been received, allowing to map every packet to a specific Genomic Dataset component.

Genomic Packet Header

Every Genomic Data Block may be split in fragments, which may be transmitted in the payload of Genomic Data Packets, depending on channel requirements, such as packet loss rate, protocol maximum packet size, etc.

A Genomic Data Packet is defined as follows.

Syntax	Description

Packet_header( ) {

Layer ID (LID)

Unambiguously identifies data type carried by this Packet.

	Unique for each sub-stream/data type. Mapping Table
	needed at beginning of stream in order to map Layer IDs to
	data types.

	Reserved	To maintain byte-alignment
	Access Unit Marker Bit (MB)	Set for the last Packet of the Access Unit. Allows identifying

the end of an AU as a set of Blocks.

Sequence Number (SN)

Packet counter, linearly increasing by 1. Needed to identify

	packet losses as gaps in SNs for each individual sub-stream.
	Associated to LID, i.e., different SN for every LID.

Packet Size (PS)

Number of bytes composing Packet, including header,

optional fields and payload.

	Extension Flag (EF)	Set if extension field is present.
	[optional] Extension field	Optional field, present if EF is set.

}

Multiplex Encoding Process

FIG. 49 shows how before being transformed in the data structures presented in this invention, raw genomic sequence data need to be mapped (491) on one or more reference sequence known a-priori (4920). In case a reference sequence is not available a “constructed” reference can be built from the raw sequence data (492). Already aligned data can be re-aligned in order to reduce the information entropy. After alignment, a genomic classifier (494) creates the data classes according to the matching functions described in Table land separates metadata (e.g. quality values) and annotation data from the genomic sequences. A reference transformation (4919) can be applied on the external reference (4920) in order to further reduce the entropy of the generated classes of data (498). The transformed data classes (4918) are fed to layers encoders (495-497) to produce genomic layers (491) which are then encoded by entropy encoders (4912-4914). The genomic streams generated by the entropy encoders are then sent to Genomic Multiplexer (4916) which generates the Genomic Multiplex. Genomic labels generated by a Genomic Labels Generator (4917) can be associated to the genomic streams (4915) by the Multiplexer (4916).