使用参考序列表示和处理生物信息数据的方法和系统
相关申请的交叉引用
本申请要求专利申请PCT/EP2016/074311、PCT/EP2016/074301、PCT/EP2016/074307、PCT/EP2016/074297、PCT/US2017/17842、PCT/US2017/17841的优先权和权益。
技术领域
本公开提供了一种表示基因组测序数据的新方法，该基因组测序数据相对于具有多个比对坐标或比对点的一个或更多个参考序列进行了比对，该表示需要将读序分割成更小的片段(也称为“拼接读序”)。所公开的表示方式通过提供已知现有技术的表示方法所不具有的新功能而减少了所用存储空间并改善了访问性能。
背景技术
基因组测序数据的适当表示是实现高效的基因组分析应用(例如基因组变异检测(variants calling)以及通过处理基因组测序数据和元数据而出于各种目的进行的所有其它分析)的基础。
由于高产出低成本测序技术的出现，人类基因组测序变得可负担得起。这种机会在从癌症的诊断和治疗到遗传疾病的辨认、从抗体辨认的病原体监测到新疫苗、药物的开发以及个性化治疗的定制的若干领域中开辟了新的前景。
医院、基因组分析提供商、生物信息和大型生物数据存储中心正在寻找可负担得起的、快速的、可靠的和互连的基因组信息处理方案，该基因组信息处理方案应该能够使基因医学发扬到世界范围。由于测序处理中的瓶颈之一已经变成数据存储问题，因此越来越多地研究以压缩形式表示基因组测序数据的方法。
测序数据的最常用的基因组信息表示基于压缩FASTQ和SAM格式。目标是压缩传统上使用的文件格式(对于非比对数据和比对数据，分别是FASTQ和SAM)。这样的文件由纯文本字符构成，并且如上所述通过使用诸如LZ(来自Lempel和Ziv，发布第一版的作者)方案(公知的zip、gzip等)的通用方法来进行压缩。当使用诸如gzip的通用压缩器时，压缩的结果通常是一大团的二进制数据。这种单体形式的信息导致很难存档、传送和说明，特别是在高产出测序的情况下，数据量非常大时。BAM格式的特点是压缩性能较差，因为它侧重于压缩低效和冗余的SAM格式，而不是提取SAM文件传达的实际基因组信息，并且采用了通用的文本压缩算法(如gzip)而没有利用每个数据源(基因组数据本身)的特定性质。
SAM的另一个重要制约是不能恰当地支持与基因组读序或读序对相关联的多个比对(也称为多个映射)的表示。基因组读序比对包括根据下一代测序(Next GenerationSequencing)技术产生的读序来重建测序样本的基因组信息的过程。可以在先前完全不知道来源基因组的情况下，或者使用已有的基因组作为参考，来进行该重建。后一种方法在本领域中被称为“基于参考的比对”。在基于参考的比对中，将由测序样本产生的基因组读序与已有的参考序列进行比较，以找到相对于该读序表现出最小数量的差异(如果有的话)的参考序列的区域。这个过程被称为将读序与参考序列进行“比对”或“映射”。
由于一些基因组区域的重复特性，读序可能以相同的精度比对至多个位置。例如，同一读序可以完美地(即，没有任何失配)与参考序列上的相同长度的两个以上片段匹配。在这种情况下，这两个以上的比对被认为是等效的，并且该读序被认为具有“多个比对”。图15例示了这种情况。在其它情况下，不同的比对可以具有不同的精度，例如，一个比对可以不存在失配(完全匹配)，而另一个比对可以具有一个或更多个失配。在这种情况下，使用评分系统来对多个比对进行排名。
在一些情况下，为了找到满足预先建立的匹配标准(例如最大失配数目)的比对位置，需要把读序分成两个以上的子片段。在这种情况下，读序被称为“拼接读序”，并且各个子片段被称为“拼接段”。图16例示了这种情况。
当前的SAM规范不支持使用11个必填字段来表示多个比对和拼接段，而是需要用到通过用于读序比对的不同工具实现以不同且因此低效的方式使用的可选字段。本文中描述的发明对于表示多个比对和拼接读序以及保持压缩和访问效率的兼顾提供了解决方案。
一种使用较少但在压缩方面比BAM更高效的更复杂的基因组数据压缩方法是CRAM。CRAM因为采用了相对于参考的差分编码而提供了更高效的压缩(它部分利用了数据源冗余)，但仍然缺乏诸如增量更新、支持流、对特定类压缩数据的选择性访问以及多个比对和拼接读序的恰当表示的特性。
这些方法产生较差的压缩比和数据结构，一旦压缩就很难浏览和操作该数据结构。下游分析可能非常慢，因为即使执行简单的操作或者访问基因组数据集的选定区域，也需要处理大且刚性的数据结构。CRAM依赖于CRAM记录的概念。每个CRAM记录通过对重构它所需的所有元素进行编码来表示单个映射读序或非映射读序。
CRAM存在以下缺陷和限制，这些缺陷和限制由本文中描述的发明解决和消除：
1.CRAM不支持数据索引和对共有特定特征的数据子集的随机访问。数据索引超出了规范的范围(参见CRAM规范v 3.0的第12节)，并且该数据索引是作为单独的文件实现的。相反，本文中描述的发明的方法采用集成有编码处理并且索引被嵌入在编码(即，压缩)比特流中的数据索引方法。
2.CRAM是由核心数据块构建的，这些核心数据块可以包含任何类型的映射读序(完全匹配的读序、仅具有替换的读序、具有插入或缺失(也称为“indel”)的读序)。不存在根据相对于参考序列的映射结果对类中的读序进行数据分类和分组的概念。这意味着即使仅搜索具有特定特征的读序，也需要检查所有的数据。本发明通过在编码之前对类中的数据进行分类和划分来解决这种限制。
3.CRAM基于将各个读序封装到“CRAM”中的概念。这意味着当搜索以特定生物特征为特征的读序(例如具有替换但不具有“indel”的读序或者完全映射的读序)时，需要检查各个完整的“记录”。
相反，在本发明中，有着在单独的信息层中单独地编码的数据类的概念，并且没有封装了各个读序的记录的概念。这使得能够更高效地访问具有特定生物特征的读序(例如具有替换但不具有“indel”的读序或者完全映射的读序)的集合，而无需解码各个读序(块)以检查其特征。
4.在CRAM记录中，由于各个CRAM记录可以包含任何不同类型的数据，因而没有上下文的概念，所以各个记录字段与特定标志相关联，并且每个标志必须总是具有相同的含义。这种编码机制引入了冗余信息且使得无法利用高效的基于上下文的熵编码。
相反，在本发明中，不存在对数据进行标注的标志的概念，因为这本质上是由数据所属的信息“层”限定的。这意味着要使用的符号数量会大大减少，进而导致信源熵的减少，这带来了更高效的压缩。因为使用不同的“层”使得编码器能够根据上下文在具有不同含义的各个层中重复使用相同符号，所以可以实现这样的改进。在CRAM中，各个标志必须总是具有相同的含义，因为没有上下文的概念且各个CRAM记录可以包含任何类型的数据。
5.在CRAM替换中，利用不同的语法元素来表示插入和缺失，这种选项会增加信息源字母表的大小并产生更高的源熵。相反，所公开的发明的方法使用单个字母表以及针对替换、插入和缺失的编码。这使得编码和解码过程更简单，并且产生更低的熵源模型，其编码产生以高压缩性能为特征的比特流。
6.CRAM不能兼顾于单端读序(single reads)和双端读序(paired-end reads)而支持多重比对以及拼接读序的恰当表示。
一旦被映射到参考序列上的一个或更多个位置，基因组读序可以与它们所映射到的参考序列片段完全匹配或者呈现一些失配。失配可以是以下类型：
·替换：映射读序中的单个核苷酸不同于参考序列上的相应核苷酸。
·插入：映射读序包含参考序列中不存在的核苷酸。在比对过程之后，它们导致映射核苷酸之间的“插入”。
·缺失：参考序列中存在的一个或更多个核苷酸不存在于映射读序中的相应位置处。
·软剪辑：读序或读序对的边缘处的一个或更多个核苷酸不映射在参考序列上，而是由比对工具保留。
·硬剪辑：读序或读序对的边缘处的一个或更多个核苷酸不映射在参考序列上，并且不被比对工具保留作为映射读序的一部分。在这种情况下，映射读序比未映射读序短。
本发明旨在通过根据读序比对过程的结果对数据进行分类并分割测序数据来压缩基因组序列，以使得待编码的冗余信息最少，并且在压缩域中直接实现诸如选择性访问和支持增量更新的特征。
当根据比对过程的结果对基因组数据进行分类时，用本发明中公开的语法元素进行的表示能够实现更高效的熵编码、对数据的选择性访问以及增量更新。增强的压缩归因于数据被分割成具有同种统计特性的独立数据块。当数据被构造在这样的可以独立解压的块中时，选择性访问需要更少的计算能力和带宽，并且可以实现增量更新，可以无需对整个数据集进行重新编码而添加新的编码数据块。
所提出的方法的一个方面是定义构成在不同层中并且分别编码的数据和元数据的类。这种方法相对于现有方法的更相关的改进包括：
1.由于通过为每类数据或元数据提供高效的源模型而构成的信源熵减小，压缩性能得到提高；
2.可以直接在压缩域中对压缩数据和元数据的部分进行选择性访问，以进行任意的进一步处理；
3.可以用与特定测序读序集相关联的新测序数据和/或元数据和/或新分析结果来对压缩数据和元数据进行增量更新(即不需要解码和重新编码)。
本发明中公开的基因组序列数据的表示方式依赖于“描述符”的概念。描述符被定义为用于表示要使用熵编码器进行压缩的基因组序列数据的语法元素。通过不同描述符来表示原始基因组序列数据使得能够更高效地压缩数据并且能够更好地对数据进行选择性访问。通过对共有相同统计特性的每种类型的描述符或描述符子集采用不同的熵编码器而实现了更高效的压缩。通过定义能够根据不同的生物学意义将基因组信息分割成数据子集的描述符，能够实现更高效的选择性访问。由于各个数据子集可以独立于其它数据进行解码，因此所需的处理能力降低并且解码时间更短。
在本发明公开中定义的描述符被构造为多个“描述符流”，这些“描述符流”是相同类型的压缩语法元素的块。
附图说明
图1示出了参考序列上的坐标系以及读序和读序对在参考序列上的映射。
图2例示了如何将基因组记录长度定义为参考序列上的将读序或读序对的最左侧映射碱基与最右侧映射碱基分隔开的基因组位置的数目。在读序对的情况下，这是当两个读序都被映射到相同的参考序列上时，参考序列上的将读序1的最左侧碱基与其配对读序2的最右侧碱基分隔开的基因组位置的数目。
图3示出了在成对读序的情况下计算基因组记录长度的示例。
图4示出了比对读序或读序对如何能跨越两个以上的访问单元。当相应的基因组记录长度大于读序或读序对映射位置离包含该记录的AU的终点的距离时，将发生这种情况。即使读序被编码在不覆盖指定基因组区域的访问单元中，一致的选择性访问算法也应该能够取得与指定基因组区域重叠的所有碱基。
图5示出了所映射的三个读序对的第一读序的位置如何被编码在“pos”描述符中。
图6示出了读序1源自链1并且读序2源自链2的读序对的示例。
图7示出了如果使用链1作为参考，如何对读序2的反向互补进行编码。
图8示出了构成读序对的读序的四种可能组合以及“rcomp”描述符中的相应编码。
图9示出了如何计算三个读序对的配对距离。
图10示出了用于对在参考序列(P类)上完全映射的固定长度读序对进行编码的描述符。
图11示出了用于对在参考序列(P类)上完全映射的可变长度读序进行编码的描述符。
图12示出了用于对相对于参考序列(N类)具有未知碱基的固定长度读序对进行编码的描述符。
图13示出了用于对在具有至少一个替换的参考序列(M类)上映射的固定长度读序对进行编码的描述符。
图14示出了用于对映射到具有至少一个插入、缺失或软剪辑碱基的参考序列(I类)上的固定长度读序对进行编码的描述符。
图15示出了针对最左侧读序具有N重比对并且针对最右侧读序具有M重比对的读序对。
图16例示了需要如何根据给定的约束来映射读序对，以将最左侧读序(读序1)分成两个拼接段并且将最右侧读序(读序2)分成四个拼接段。
图17示出了没有拼接的多重比对。对于每个读序或读序对，N是mmap描述符的第一值，并且通报了第一读序的比对数目。在双端读序的情况下，使用mmap描述符的以下N个值来计算作为第二读序的比对数目的P。
图18示出了如何使用pos、pair和mmap描述符对不具有拼接段的多重比对进行编码。
图19示出了具有拼接段的多重比对的示例。N表示第一读序的拼接段数目，并且被编码为mmap描述符的第一个值。P表示第二读序的拼接段数目，并且使用mmap描述符的接下来的N个值来计算。N1和N2表示第一读序和第二读序的比对数目，并且使用splen描述符的N+P个值来进行计算。
图20示出了如何使用pos、pair、mmap和splen描述符来表示具有拼接段的多重比对。
图21示出了可以用于在根据本发明公开内容的存储介质上写入编码基因组信息的两种编码模式。当采用访问单元毗连(Access Units Contiguous，AUC)模式时，把各个访问单元存储在存储介质的毗连区域上。当采用描述符流毗连(Descriptor StreamContiguous，DSC)模式时，相同类型的所有描述符作为单个块存储在存储介质的毗连区域上。因此，每个AU“分布”在存储介质的不同片段中。
图22示出了如何由编码描述符块构成访问单元。为了能够在网络上传输，块被进一步分解为分组。
图23示出了由报头以及不同流(每个描述符对应一个)的数据块的复用构成的访问单元。几个相同类型的分组被封装在一个块中，并且在一个访问单元中复用多个块。
图24示出了如何能通过对读序进行聚类并组装从各个群集(cluster)取得的片段来构造“内部”参考。每个群集贡献于所构造的参考的一个片段。
图25示出了构造参考的策略在于，一旦将特定排序(例如，词典顺序)应用于读序，就存储最近的读序。
图26示出了可以使用在相应流中存储或携带的六个描述符对属于未映射读序类(U类)的读序进行编码。
图27示出了U类读序的替代编码，其中，使用带符号的pos描述符来对读序在所构造的参考上的映射位置进行编码。
图28示出了实施本发明公开内容的基因组编码器的架构。
图29示出了实施本发明公开内容的基因组解码器的架构。
图30示出了具有读序1的三重比对以及读序2的四重比对的读序对。所有的比对都在同一条染色体(例如，染色体1)上，而读序1的第三比对在染色体X上。使用pair描述符的第一个值来通报在与其它不同的参考序列上存在一个比对。
图31示出了使用mmsc描述符来表示一个不保持主比对的映射毗连性的辅比对的示例。
图32示出了如何使用半映射读序对(HM类)通过用未映射读序组装更长序列(也称为“重叠群(contigs)”)来填充参考序列的未知区域。
发明内容
权利要求的特征通过提供以下项解决了现有技术方案的问题：
一种用于编码基因组序列数据的方法，所述基因组序列数据包括核苷酸序列的读序，所述方法包括以下步骤：
将所述读序与一个或更多个参考序列进行比对，从而创建比对读序，
根据指定的与所述一个或更多个参考序列匹配的规则将所述比对读序分类，从而创建比对读序的类，
将所述分类比对读序编码为语法元素的多个流，
其中，将所述分类比对读序编码为语法元素的多个流的步骤包括以下步骤：根据所述比对读序的类来选择所述语法元素，
向所述语法元素流提供报头信息，从而创建连续的数据块，以将所述基因组数据块熵编码成能够单独访问的数据单元。
在另一方面中，所述编码方法还包括以下步骤：
将不满足所述指定匹配规则的所述读序分类成未映射读序的类，
将所述分类后的未映射读序编码为语法元素的多个流，
向所述语法元素流和所述编码参考序列提供报头信息，从而创建连续的访问单元。
在另一方面中，所述编码方法还包括：所述分类步骤包括辨别在用于比对的参考序列上具有多个比对位置的基因组读序。
在另一方面中，所述编码方法还包括：所述分类步骤包括以下步骤：辨别需要分成称为拼接段的多个片段以满足比对的匹配规则的基因组读序。
在另一方面中，所述编码方法还包括：待编码的基因组序列的读序是成对的。
在另一方面中，所述编码方法还包括以下步骤：
辨认读序的比对数目并且用特定语法元素表示该数目，
对于每个比对辨认对应的映射位置并且用特定语法元素表示各个映射位置，
将比对分数分配给各个比对，以标识主比对和辅比对，
将具有最高分数的比对标识为主比对，
辨别是否在不同于主比对的参考上找到任何比对，并且使用特定描述符来表示这个信息，
辨别是否有任何比对未保留所述主比对在参考序列上的不同毗连性，并且使用特定语法元素来表示这个信息。
在另一方面中，所述编码方法还包括以下步骤：
辨别需要分成两个以上的拼接段以根据规定与所述一个或更多个参考序列的匹配的预定匹配规则比对到参考序列上的读序，
使用全局配置参数通报存在拼接读序，
使用特定语法元素来表示拼接段的数目，
使用特定语法元素来表示各个拼接段的长度。
在另一方面中，所述编码方法还包括以下步骤：
辨别读序对中的各个读序的比对数目并且用特定语法元素表示该数目，
对于读序对中的最左侧读序的每个比对，辨别对应的映射位置，并且用特定语法元素表示各个映射位置，
对于最左侧读序的每个比对，辨别读序对中的最右侧读序的关联的比对，并且用特定语法元素表示该关联，
将比对分数分配给各对比对，以标识主比对和辅比对，
将具有最高分数的一对比对标识为主比对，
辨别是否在不同于主比对的参考上找到任何比对，并且使用特定描述符来表示这个信息，
辨别是否有任何比对在参考序列上表现出不同于主比对的毗连性，并且使用特定语法元素来表示这个信息。
在另一方面中，所述编码方法还包括以下步骤：
辨别需要分成两个以上的拼接段以便根据预定匹配规则比对到参考序列上的读序，
使用全局配置参数通报存在拼接读序，
使用特定语法元素表示读序对中的最左侧读序的拼接段的数目，
用特定语法元素的向量表示与最左侧读序的各个比对相关联的最右侧读序的拼接段的数目，
使用特定语法元素来表示各个拼接段的长度。
在另一方面中，所述编码方法还包括：
所述语法元素流包括具有以下项的基因组数据集报头：
用于唯一地标识各个数据集组的数据集组标识符，
用于唯一地标识各个数据集的基因组数据集标识符，
用于标识数据集所遵循的数据格式规范的品牌标识符，
用于标识数据集所遵循的数据格式规范的次要版本号，
用于通报固定长度读序的以核苷酸为单位的编码基因组读序长度，
通报存在双端读序的标志，
通报存在块报头的标志，
通报访问单元被存储在存储介质上的模式以在解码所述访问单元时便于数据访问的标志，
用于对读序相对于参考序列的失配进行编码的字母表类型，
用于编码数据集的参考序列的数目，
用于唯一地标识各个参考序列的对于每个参考序列的数字标识符，
用于唯一地标识各个参考序列的对于每个参考序列的字符串标识符，
用于对与各个参考序列相关联的访问单元进行计数的对于每个参考序列的编码访问单元数目，
用于区分比对读序、未比对读序、未映射读序和参考序列的编码基因组数据类型，
编码在数据集中的数据类的数目，
在解码处理时使用的数据集中所编码的针对每个数据类使用的描述符的数目，
用于对编码未映射读序进行索引的群集的总数，
用于表示如下整数值的比特数，该整数值用于对编码群集签名的解码所用的群集签名进行编码，
通报是否所有的群集签名就核苷酸数目而言具有相同长度的标志，
群集签名的长度。
在另一方面中，所述编码方法还包括：
所述语法元素流包括主索引表，该主索引表针对每个类包含一个区段，并且包含比对读序的子类，所述区段包括：
各个数据类或数据子类的各个访问单元的最左侧读序的主比对在所述一个或更多个参考序列上的映射位置，
各个数据类或数据子类的各个访问单元的所有主比对中的最右侧映射碱基在所述一个或更多个参考序列上的位置，
构成各个访问单元的各个编码语法元素块的以字节为单位的偏移。
在另一方面中，所述编码方法还包括：所述主索引表还包含各个编码描述符流的大小。
在另一方面中，所述编码方法还包括：所述主索引表还包含各个访问单元的大小。
在另一方面中，所述编码方法还包括：基因组读序具有多重比对，并且所述主索引表包含：
各个数据类或数据子类的各个访问单元的所有读序中的最左侧比对在所述一个或更多个参考序列上的映射位置，
各个数据类或数据子类的各个访问单元的所有比对中的最右侧映射碱基在所述一个或更多个参考序列上的位置，
构成各个访问单元的各个编码语法元素块的以字节为单位的偏移。
在另一方面中，所述编码方法还包括：所述访问单元包含编码读序对。
在另一方面中，所述编码方法还包括：所述主索引表与所述访问单元数据一起编码。
在另一方面中，所述编码方法还包括：所述基因组数据集报头与所述访问单元数据一起编码。
在另一方面中，所述编码方法还包括：所述语法元素流还包括与下列项相关的信息：所使用的参考的类型(已有的或构造的)，以及读序的未匹配到参考序列上的片段。
在另一方面中，所述编码方法还在于，将所述分类后的比对读序编码为多个语法元素流的步骤包括以下步骤：将特定源模型和特定熵编码器与每个描述符流关联起来。
在另一方面中，所述编码方法在于，所述熵编码器是上下文自适应算术编码器、可变长度编码器和golomb编码器中的一种。
一种用于对编码基因组数据进行解码的方法，所述方法包括以下步骤：
解析包含所述编码基因组数据的访问单元，以利用报头信息来提取多个语法元素流，
对所述多个语法元素流进行解码，以根据特定匹配规则来提取比对读序，所述匹配规则定义了所述比对读序相对于一个或更多个参考序列的分类。
在另一方面中，所述解码方法还包括对未映射基因组读序进行解码。
在另一方面中，所述解码方法还包括对包含全局配置参数的基因组数据集报头进行解码。
在另一方面中，所述解码方法还包括对主索引表进行解码，所述主索引表对于每类读序包含一个区段，并且包含相关联的相关映射位置和编码块偏移。
在另一方面中，所述解码方法还包括对与所使用的参考的以下类型有关的信息进行解码：已有的、变换的或构造的。
在另一方面中，所述解码方法还在于，所述基因组读序是成对的。
在另一方面中，所述解码方法还在于，所述基因组数据是被熵解码的。
在另一方面中，所述解码方法还在于，多重比对信息的解码包括以下步骤：
对每个读序的比对数目进行解码，
对每个比对的位置进行解码，
通过解码与每个比对相关联的分数来辨别主比对，
通过解码对应的语法元素来辨别是否有任何辅比对相对于参考序列具有与主比对不同的毗连性。
在另一方面中，所述解码方法还包括以下步骤：
辨别是否有任何编码读序被分成两个以上的拼接段，
对每个拼接段的长度进行解码，
对各个拼接段的映射位置进行解码。
在另一方面中，所述解码方法还在于，编码基因组读序是双端读序，并且所述解码方法还包括以下步骤：
对与最左侧读序的各个比对相关联的最右侧读序的比对的数目进行解码，
对如下的配对信息进行解码：该配对信息将最左侧读序的各个比对与最右侧读序的一个或更多个比对关联起来。
在另一方面中，所述解码方法还在于，编码基因组读序被分成两个以上的拼接段，并且所述解码方法还包括以下步骤：
对各个编码拼接段的长度进行解码，
对各个拼接段的映射位置进行解码。
本发明还提供了一种用于压缩基因组序列数据(289)的基因组编码器(2810)，所述基因组序列数据(289)包括核苷酸序列的读序，所述基因组编码器(2810)包括：
比对单元(281)，其被配置为将所述读序与一个或更多个参考序列进行比对，从而创建比对读序，
构造参考生成单元(282)，其被配置为产生构造参考序列，
数据分类单元(284)，其被配置为根据与一个或更多个已有参考序列或构造参考序列进行匹配的指定匹配规则对所述比对读序进行分类，从而创建比对读序的类(288)，
一个或更多个描述符流编码单元(285-287)，其被配置为通过根据所述比对读序的类选择所述语法元素来将所述分类比对读序编码为语法元素流，
复用器(2816)，其用于复用压缩基因组数据和元数据。
在另一方面中，所述基因组编码器还包括：
数据分类单元(284)包含N类、M类和I类数据的编码器，该编码器配置有生成N类、M类和I类数据的子类的阈值向量。
在另一方面中，所述基因组编码器还包括用于执行前述编码方法的所有方面的特征。
本发明还提供了一种用于对压缩基因组流(291)进行解压的基因组解码器(298)，所述基因组解码器(298)包括：
解复用器(290)，其用于分解压缩基因组数据和元数据，
解析装置(292-294)，其被配置为将所述压缩基因组流解析成语法元素流(295)，
一个或更多个语法元素流解码器(296-297)，其被配置为将描述符流解码为核苷酸序列的分类读序(2911)，
基因组数据类解码器(299)，其被配置为选择性地解码一个或更多个参考序列上的核苷酸序列的所述分类读序，以产生核苷酸序列的未压缩读序。
在另一方面中，所述基因组解码器还包括将所述一个或更多个参考序列存储在压缩基因组流(291)中。
在另一方面中，所述基因组解码器还包括通过带外机制将所述一个或更多个参考序列提供给所述解码器。
在另一方面中，所述基因组解码器还包括在所述解码器处构造所述一个或更多个参考序列。
本发明还提供了一种包括指令的计算机可读介质，该指令在被执行时使得至少一个处理器执行前述编码方法的所有方面。
本发明还提供了一种包括指令的计算机可读介质，该指令在被执行时使得至少一个处理器执行前述解码方法的所有方面。
本发明还提供了存储有根据前述编码方法的所有方面编码的基因组的支持数据。
具体实施方式
本发明所指的基因组或蛋白质组序列例如包括但不限于，核苷酸序列、脱氧核糖核酸(DNA)序列、核糖核酸(RNA)和氨基酸序列。虽然本文的说明相当详细地描述了核苷酸序列形式的基因组信息，但应理解，用于压缩的方法和系统也可实施于其它基因组或蛋白质组序列，虽然如本领域技术人员所能理解的那样将有一些改动。
高通量测序(High Throughput Sequencing，HTS)机产生形式为核苷酸序列(也称为“碱基”)的基因组测序信息，该核苷酸序列由所规定的词汇表中的字母的串表示。最小的词汇表是由五个符号{A，C，G，T，N}表示的，其表示DNA中存在的四种核苷酸，即腺嘌呤、胞嘧啶、鸟嘌呤和胸腺嘧啶。在RNA中，胸腺嘧啶被尿嘧啶(U)取代。N表示测序机不能检测到任何碱基，因此该位置处的核苷酸的真实性质是不确定的。在测序机采用IUPAC模糊码作为词汇表的情况下，符号所用的字母表由以下符号组成：{A，C，G，T，U，W，S，M，K，R，Y，B，D，H，V，N或-}。
在本发明的语境中，“基因组数据集”被定义为任何结构化的基因组数据的集合，该基因组数据包括例如活生物体的基因组数据、通过活生物体的基因组测序或通过对原始测序数据执行的任何其他基因组数据处理步骤生成的一个或更多个序列和元数据。
测序机产生的核苷酸序列称为“读序”。读序可以由几十个至几千个核苷酸构成。一些测序技术产生由“读序对(pairs)”构成的读序，其中一个读序源自一条DNA链，另一个读序源自另一条链。在产生读序对的测序过程中与另一个读序相关联的读序被称为该另一个读序的“配对读序”。
在本公开中，参考序列是与一维整数坐标系相关的核苷酸序列，对于该一维整数坐标系，各个整数坐标与单个核苷酸相关联。坐标值只能等于或大于零。在本发明的语境中，这个坐标系是基于零的(即，第一个核苷酸的坐标为0，并且它被称为处在位置0处)并且从左到右线性增加。
当在参考序列上映射读序时，所述参考序列被用作最左位置被标为位置0的一维坐标系的轴。在映射到参考序列的读序中，构成在由最小坐标编号标识的参考序列位置处映射的读序的核苷酸通常被称为“最左侧”核苷酸，而构成在由最大坐标编号标识的参考序列位置处映射的读序的核苷酸被称为“最右侧”核苷酸。这在图1中示出。在本公开中，核苷酸也称为“碱基”。
当读序被映射到参考序列时，最左侧的映射碱基的坐标被认为表示该读序在参考序列上的映射位置。
参考基因组由一个或更多个参考序列组成，并由科学家组装成物种基因组的代表性示例。例如，GRCh37(基因组参考协会人类基因组(第37版))源自来自纽约布法罗的13名匿名志愿者。然而，参考序列也可以由合成序列组成，该合成序列是被构想出来的，并且只是被构成为考虑读序的进一步处理而改善其压缩性。
在本公开中，基因组记录被定义为以下项的编码表示：
·可选地与比对信息、读序标识符和质量值相关联的单读序或读序对
·参考序列(例如染色体)或其一部分。
在本公开中，基因组记录位置被定义为在基因组记录中编码的读序或读序对的最左侧映射碱基在参考序列上的位置。
读序或读序对记录中的映射碱基包括：
·与参考序列上的对应碱基匹配的比对读序的碱基
·与对应碱基不匹配(也称作单核苷酸多态性)的比对读序的碱基
存在于比对读序中但不存在于参考序列中(也称作插入)的碱基以及比对过程保留的未映射到参考序列上(也称作软剪辑)的碱基没有映射位置。
在本公开中，将构成具有映射在参考序列的最小坐标上的碱基的读序对的读序称为“读序1”，而将其配对读序称为“读序2”。
由使用了当今技术的测序机作为一对而生成的两个读序之间隔开的由核苷酸(或碱基)的数目表示的距离是未知的，通过将组成读序对的两个读序映射到参考序列(即，使适宜的匹配函数最小化)来确定该距离。
在本公开中，基因组记录长度被定义为在记录中编码的最左侧映射碱基与在记录中编码的最右侧映射碱基之间的坐标位置的数目。
在本公开中，配对距离被定义为在记录中编码的最左侧映射碱基与在记录中编码的读序2的最左侧映射碱基之间的坐标位置的数目。图1示出了配对距离的示例。
在本公开中，在单读序的情况下，基因组记录长度(GRL)被计算为通过从最右侧映射碱基的映射位置减去最左侧碱基的映射位置并加上“1”而获得的整数。
GRL＝(最右侧碱基位置)-(最左侧碱基位置)+1
在本公开中，在读序对的情况下，基因组记录长度(GRL)被计算为通过从配对读序(读序2)的最右侧碱基映射位置的坐标中减去映射在参考序列(读序1)的最小位置处的读序的最左侧碱基映射位置并加上“1”而获得的整数。图3例示了基因组记录长度的这种定义。
在本公开中，在编码了参考序列或其一部分的基因组记录的情况下，基因组记录长度被定义为组成编码序列的核苷酸的数目。
在本公开中，基因组范围被定义为由起始坐标S和终止坐标E(S≤E)定义的参考序列上的连续区间。基因组范围的起始和终止位置总是包括在该范围中。
本发明旨在定义一种新方法，使得在通过包括在多个被称为访问单元的数据块中的描述符集来压缩读序时，能够高效地访问映射到任何基因组区域的比对基因组读序。
在本公开中，访问单元(Access Unit，AU)被定义为逻辑数据结构，该逻辑数据结构包含基因组信息或相关元数据的编码表示以便于比特流访问和操作。访问单元是可由实施本公开所描述的发明的解码装置解码的最小数据组织。
根据编码信息的类型，可以独立于任何其它AU或者使用包含在其它AU中的信息来解码AU。
根据编码数据的性质，AU可以是多种类型的。访问单元包含参考序列或其一部分或属于单个数据类的编码读序或编码读序对。例如，访问单元可以包含GRCh37(基因组参考协会人类基因组(第37版))的整个染色体1。另一访问单元可以包含GRCh37的染色体1的位于坐标50'000与150'000之间的核苷酸的编码表示。另一访问单元可以仅包含完美地而没有任何失配地映射到参考序列上的读序或读序对。另一访问单元可以包括仅包含作为相对于参考序列的失配的“N”符号的读序或读序对。另一访问单元可以包括具有任何类型的替换(例如，存在于读序或读序对中的一个碱基与参考序列中的对应映射位置处的碱基不同)的读序或读序对。另一访问单元可以包括具有失配、插入、缺失和软剪辑碱基的读序或读序对。另一访问单元可以仅包含未映射到参考序列上的读序或读序对。另一访问单元可以仅包含这种读序对：一个读序被映射至参考序列上且另一读序未被映射至参考序列上。另一种类型的访问单元可以仅包含由一个或更多个参考序列(例如染色体)组成的参考基因组的编码片段。
访问单元的关键特征在于，它以压缩形式包含重构基因组信息(读序或读序对、参考序列)所需的所有元素、相关联的比对信息以及该访问单元表示的读序或读序对的元数据。换句话说，为了完全重构由访问单元携带的读序、读序对或参考序列以及相关联的信息，仅需要取得访问单元本身，并且在需要时取得包含所参考的参考序列的访问单元。
在每个访问单元中，在下一节中列出的表示编码读序或读序对的描述符被聚集在单独的数据块(一种类型对应一个数据块)中，以便在应用熵编码来压缩它们时利用它们的同种统计特性。
每个访问单元包含表示如下读序或读序对的描述符的压缩子集：该读序或读序对属于映射到参考序列上的基因组区域的同一类。参考序列上的这种基因组区域由起始坐标(或起始位置)和终止坐标(或终止位置)定义。
图22例示了访问单元的示例。访问单元由编码描述符的块(在下一节中描述)组成。为了能够在网络上传输，块被进一步分解为分组。
描述符是表示重建(即，解码)编码参考序列、编码读序及相关映射信息或编码读序对及相关映射信息所需的信息的一部分的语法元素。不同类型的描述符被定义为表示以下项：
·读序在参考序列上的映射位置，
·读序与其配对读序之间的距离，
·读序的长度，
·比对读序中的相对于参考序列的失配的位置，
·相对于关联位置处的参考序列的失配的类型，
·不能通过映射过程映射在参考序列上并且被分类为“软剪辑“碱基的碱基，
·读序长度，
·SAM规范指定的映射标志，
·通过映射过程与单读序或读序对相关联的多个映射位置，
·识别拼接读序(即，这样的读序：与作为映射在参考序列的单个位置上的单个连接读序进行映射的情况相比，当分成几段时，能够获得更高匹配精度的映射位置)的存在，
·所使用的参考序列的具体类型，比如：
ο参考基因组，例如由如基因组参考协会(例如，GRCh37)、加州大学圣克鲁斯分校(例如，hg19)的协会公布的基因组，
ο使用指定的读序集和指定的组装规则集构建的参考序列，
·应用于参考序列的修订的位置和类型(为了减小用于对映射在这样修订的参考序列上的读序的失配进行表示的描述符的熵)，
·不能以指定的匹配精度映射在参考序列的任何位置处的读序的表示，
·整个参考序列或其部分的表示。
下面的章节提供了本公开中提及并且被本发明使用的各个描述符的完整列表和精确定义。
如上所述，本发明使用多种描述符来表示压缩形式的参考序列、读序或读序对，使得可以用它们的的相关信息对其进行完全重构。在读序或读序对也被分类且分至不同的类的情况下，根据在参考序列上的映射结果以及被编码成单独数据块的熵，不同的描述符子集用于表示各类读序或读序对。因此，访问单元仅包含表示参考序列(或其部分)或者属于同一类的读序或读序对所需的那些熵编码描述符。图11示出了针对可变长度的读序的这种情况，图12、图13和图14示出了针对固定长度的读序对的这种情况。
在本公开中，相同类型的熵编码描述符被认为构成描述符流。
在本公开中，访问单元起始位置被定义为访问单元包含的所有基因组记录中的最左侧基因组记录位置。
在本公开中，访问单元终止位置被定义为访问单元包含的所有基因组记录的所有映射碱基中的最右侧映射碱基的位置。
在本公开中，访问单元范围被定义为包含在访问单元中的所有基因组记录中的AU起始位置与最右侧基因组记录位置之间的基因组范围。可以通过从AU终止位置减去AU起始位置并加“1”来计算其位置数量的值。
在本公开中，访问单元覆盖区域被定义为包括在AU起始位置与AU终止位置之间的基因组范围。
在本公开中，访问单元也被认为覆盖起始位置与其终止位置之间的基因组区域。
在AU中编码的一些基因组记录可以在距AU终止位置一定距离(小于基因组记录长度)处具有映射位置。这意味着属于在基因组记录中编码的读序或读序对的一些碱基被映射在由后续AU中的一个覆盖的基因组区域中。这在图4中示出。
根据上文提供的定义，可以辨认构建访问单元的两种有用的方式：
1.所谓的“非重叠模式”，其中相同数据类的访问单元的基因组范围从不重叠，
2.所谓的“重叠模式”，其中相同数据类的访问单元的基因组范围可以重叠。
“非重叠模式”在基因组数据被压缩并作为文件存储在存储器中的情况下以及在存储文件从一个存储装置流传输到另一个存储装置的流传输情况下是优选的，因为它使得能够实现高效的选择性访问过程。“重叠模式”支持这样的流传输情形：在整个基因组序列数据于发送装置处可用之前，可以把基因组数据集的部分编码到访问单元中并进行传输。
所公开的发明的主要创新方面如下：
1通过表示编码读序或读序对如何映射到参考序列以及如何根据其压缩表示进行完全重构的描述符子集来实现将各个读序或读序对编码到基因组记录中。
2使用特定的描述符子集来表示属于各个数据类(其中根据参考序列上的映射结果对读序或读序对进行分割)的读序或读序对。这种表示支持多个映射位置和拼接读序。
3特定的描述符子集用于表示参考序列或其部分。
4在读序对(每个读序具有相同的长度)的情况下，每个读序对被熵编码为与单个描述符子集相关联的基因组记录中的单个实体，以使得编码性能最佳。
5在基因组记录中编码的压缩形式的(即，熵编码)描述符子集包括在具有如下范围的访问单元中：该范围包括读序的最左侧映射碱基或者读序对的读序1的最左侧映射碱基。
6定义了如下编码方法：该方法确定属于一个读序对的读序是否必须编码为一体并且表示该读序对的描述符子集必须以压缩形式包括在单个访问单元中，或者该读序对是否必须被分开并且以压缩形式表示为两个不同访问单元中的两个单独的读序。该方法定义了使用读序对中的两个读序之间的距离与访问单元范围作为输入参数的分割规则。
7在编码器侧应用该新的编码方法并且在压缩比特流中传输编码参数使得在解码器侧能够确定在需要访问映射到(即，属于)特定基因组区域的所有读序和读序对时所要解码的访问单元的最小数量。
8编码方法的定义、编码方法的参数在解码器侧的传输、最小访问单元数的确定以及为取得映射在任何给定基因组区域上的所有读序而对访问单元进行的辨别使得能够对以压缩形式存储或传输(即流传输)的读序和读序对进行高性能的选择性访问。
9由于使用了描述符并在访问单元内将其划分为同种数据块，可以应用高效的熵编码来压缩基因组信息表示。
10被称为主索引表(Master Index Table，MIT)的索引工具的定义包括：
a.对于各类数据，包括了基因组数据集中包含的每个AU的起始位置和终止位置的索引。在每个索引中，AU起始位置以升序进行排序。
b.对于各类数据，指向属于各类数据的编码AU在存储介质上的物理位置的指针的矢量。每个矢量被作为AU起始位置的相应索引进行排序。
MIT通过将参考序列上的基因组区域与包含读序或读序对的压缩表示的访问单元在存储装置上的相应位置关联起来，实现了对映射到指定基因组区域的读序或读序对的高效随机访问。
11在具有可变长度的单读序(即长度可为任意值的读序)的情况下，评估在每个访问单元中编码的读序的长度分布及其映射坐标的新编码方法确定属于访问单元中包含的读序的碱基的映射位置(AU终止位置)的坐标的最大值、将这样的值的表示写入MIT中，以使得解码器可以实现对特定基因组区域的高效的选择性访问操作。通过仅检查MIT以取得覆盖指定的选择性访问区域的各个AU的终止位置的值，解码器可以确定需要解码的访问单元的最小数量(该最小数量用以保证取得覆盖需要选择性访问的基因组区域的所有读序)。
在下文中，将进一步详细地描述上述创新方面中的每一个。
根据匹配规则对读序进行分类
根据相对于一个或更多个“已有”参考序列的比对的匹配结果，所公开的发明将测序机产生的读序分类为六个不同的“类”。
当相对于参考序列比对核苷酸的DNA序列时，可以辨别以下情况：
1.找到参考序列中的区域来毫无误差地匹配读序(即完全映射)。这样的核苷酸序列称为“完美匹配读序”或标为“P类”。
2.找到参考序列中的区域来与这样的读序进行匹配：该读序具有仅由生成该读序的测序机无法检测到任何碱基(核苷酸)的位置的数量决定的失配类型和数目。这种类型的失配用“N”表示，该字母用于表示未定义的核苷酸碱基。在本文中，这种类型的失配被称为“n型”失配。这样的序列属于“N类”读序。一旦读序被分类为属于“N类”，那么将匹配不精确度限制到给定的上限并且在被认为有效匹配与被认为无效匹配之间设定分界是有用的。因此，也通过设定阈值(MAXN)来约束分配至N类的读序，该阈值(MAXN)规定了读序可以包含的未定义碱基(即，被称为“N”的碱基)的最大数目。这种分类隐含地定义了当与对应的参考序列进行参照时属于N类的所有读序共有的所需最小匹配精度(或最大失配度)，这构成了用于对压缩数据应用选择性数据搜索的有用标准。
3.找到参考序列中的区域来与这样的读序进行匹配：该读序具有由生成该读序的测序机无法检测到任何(如果有的话)核苷酸碱基(即“n型”失配)的位置的数量决定的失配类型和数目，加上已经检测到与参考序列中所存在的不同的碱基的失配的数量。标为“替换”的这种类型的失配也称为单核苷酸变异(SNV)或单核苷酸多态性(SNP)。在本文中，这种类型的失配也称为“s型”失配。然后将读序称为“M失配读序”并分配给“M类”。就像“N类”的情况一样，对于属于“M类”的所有读序而言，将匹配不精确度限制到给定的上限并且在被认为有效匹配与被认为无效匹配之间设定分界也是有用的。因此，分至M类的读序也通过定义一组阈值而受到约束，一个阈值用于“n型”失配的数目“n”(MAXN)(如果有的话)，另一阈值用于替换数“s”(MAXS)。第三约束是由数字“n”和“s”二者的任意函数f(n，s)定义的阈值。这种第三约束使得能够根据任何有意义的选择性访问标准生成具有匹配不精确度的上限的类。例如但不作为限制，f(n，s)可以是(n+s)1/2或(n+s)或任何线性或非线性表达式，它们对属于“M类”的读序所允许的最大匹配不精确度水平设定了边界。当出于各种目的而分析读序时，这样的边界构成了用于对压缩数据应用期望的选择性数据搜索的非常有用的标准，因为这使得可以在应用于一种类型或另一种类型的简单阈值之外对“n型”失配和“s型”失配(替换)的数目的任何可能组合设定进一步的边界。
4.第四类由表现出“插入”，“缺失”(也称为indel)和“剪辑”中的任何类型的至少一种失配加上属于N类或M类的任何失配类型(如果有的话)的读序构成。这种序列称为“I失配读序”并被分至“I类”。插入由不存在于参考序列中但存在于读序中的一个或更多个核苷酸的额外序列构成。在本文中，这种类型的失配被称为“i型”失配。在文献中，当插入序列位于序列的边缘时，它也被称为“软剪辑”(即，核苷酸与参考序列不匹配，但被保留在比对读序中，这与被丢弃的“硬剪辑”核苷酸相反)。在本文中，这种类型的失配被称为“c型”失配。缺失是读序中的相对于参考读序的“空洞”(缺少核苷酸)。在本文中，这种类型的失配被称为“d型”失配。与“N类”和“M类”的情况一样，为匹配不精确度规定一个极限是可行而且恰当的。“I类”约束集的定义基于“M类”所使用的相同原则并在表1的最后一个表格行中进行了报告。除了I类数据允许的每种类型的失配的阈值之外，还通过由失配数“n”、“s”、“d”、“i”和“c”的任何函数w(n，s，d，i，c)确定的阈值定义了另一约束。这种附加约束使得可以根据任何有意义的用户定义的选择性访问标准来生成具有匹配不精确度上限的类。例如但不是作为限制，w(n，s，d，i，c)可以是(n+s+d+i+c)1/5或(n+s+d+i+c)或任何线性或非线性表达式，它们对属于“I类”的读序所允许的最大匹配不精确度水平设定了边界。当出于各种目的而分析读序时，这样的分界构成了用对对压缩数据应用期望的选择性数据搜索的非常有用的标准，因为这使得可以在应用于各种类型的可允许失配的简单阈值之外对“I类”读序允许的失配数目的任何可能组合设定进一步的边界。
5.第五类包括在对参考序列进行参照时对于每个数据类都不能找到任何被认为有效的映射的所有读序(即，不满足规定了表1指定的最大匹配不精确度的上限的匹配规则集)。当对参考序列进行参照时，这些序列被认为“未映射”，并且被分类为属于“U类”。
根据匹配规则对读序对进行分类
上一节中说明的分类涉及单读序。在成对地生成读序(其中已知两个读序被长度可变的未知序列分开)的测序技术(即Illumina Inc.)的情况下，将整个读序对分类到单个数据类是合适的。与另一读序相结合的读序被称为其“配对读序”。
如果配对的两个读序属于相同的类，则整个读序对的分类是显然的：对于任何类(即，P、N、M、I、U)，将整个读序对分至同一类。如果两个读序属于不同的类，但它们都不属于“U类”，则将整个读序对分至具有根据下式规定的最高优先级的类：
P<N<M<I
其中，“P类”的优先级最低，“I类”的优先级最高。
如果只有一个读序属于“U类”且其配对读序属于P、N、M、I类中的任何一类，则定义代表“半映射”的第六类“HM类”。
规定这种特别的读序类别是想到它可用于试图确定参考基因组中存在的缺口或未知区域(也称为鲜为人知或未知的区域)。通过使用可以映射在已知区域上的读序对在边缘处映射读序对来重构这样的区域。然后使用未映射的配对读序来构造如图32所示的未知区域的所谓的“重叠群(contigs)”。因此，仅对这种类型的读序对提供选择性访问大大减少了相关的计算负担，从而能够更高效地处理由使用现有技术解决方案将需要完全检查的大量数据集所产生的这种数据。
下表总结了应用于读序的匹配规则，以便定义各个读序所属的数据类。在表的前五列中就失配类型(n型失配、s型失配、d型失配、i型失配和c型失配)的存在或不存在定义了规则。第六列就各个失配类型的最大阈值和可能的失配类型的任何函数f(n，s)和w(n，s，d，i，c)提供了规则。
表1各个读序必须满足以被分类在本发明公开定义的数据类中的失配类型和约束集
基因组读序和参考序列的压缩表示
基因组读序压缩的高效方法的共同要素是利用序列数据相对于参考序列的相关性。即使人群的体况极其多样化，人与人之间不同的核苷酸数的实际部分仅约占构成整个基因组的核苷酸总数的0.1％。因此，表征各个个体的特定基因组信息相对于整个基因组携带的全部信息非常有限。当已有的参考基因组可用(无论是用于先前的测序还是作为公开的“平均”共识参考)时，编码实际信息的最高效方式是仅对相对于参考基因组的差异进行辨别和编码。
为了以FASTQ数据形式的原始读序进行这种操作，执行在可用参考基因组上进行映射的初步预处理步骤。在参考基因组不可用的情况下或者如果不希望由于使用特定参考序列而引入偏差，则通过将可用读序组装成更长的序列来构造新参考序列是可能的替代方案。
当已经相对于已有的或构造的参考序列映射了读序时，每个读序可以由本公开中标为“读序描述符”或简称为“描述符”的多个元素完全表示。
例如，在读序与参考序列的片段完美匹配的情况下，表示该读序所需的描述符子集只是由参考序列上的映射位置的坐标(通常是读序的最左侧碱基的映射位置的坐标)、该读序本身的长度和指示该读序相对于参考序列链是映射在正还是负DNA链上的信息构成。
在不能找到任何读序的所有碱基与参考序列的所有碱基匹配的映射位置的情况下，保留具有最小失配数的一个或多个映射。在这种情况下，需要不同的描述符子集来进一步表示对应于具有最小失配数或接近最小失配数的映射位置而可能发生的替换、插入、缺失和剪辑碱基。利用这样的描述符子集，可以使用由描述符携带的信息和由参考序列携带的信息来重构读序。
映射处理还可以产生其它类型的信息，例如：多个可能的映射位置和相关分数，映射质量，拼接读序的规范，属于一个读序对的读序在两个不同的参考序列(通常是染色体)上的映射，测序处理(例如PCR或光学重复)的特征。所有这些信息需要特定的附加描述符，该附加描述符扩展各个子集，然后通过对各个描述符子集应用恰当的熵编码算法来压缩各个子集。
由于以下原因，基因组测序处理可以生成读序副本(即，同一基因组序列的两个或更多个精确的副本)：
·基因组测序处理的化学(聚合酶链式反应副本)，
·数据采集处理(光学重复)。如果读序对在同一个区块上，并且读序之间的距离小于取决于实验的给定配置参数，则读序被称为光学副本。
因此，各个读序或读序对可以根据映射处理的结果由特定的描述符子集唯一地表示。
诸如SAM和CRAM的常用方法不根据表示其映射信息所需的特定描述符子集来编码读序或读序对。SAM和CRAM不会根据它们所包含的相对于它们所映射到的参考序列的失配的数量和类型将读序分类为数据类。此外，这些格式不将读序分别编码为仅包含压缩形式的属于单个数据类的读序的访问单元中。在成对生成的读序的情况下，现有技术的方法不将它们编码为根据它们相对于参考序列的映射精确度分成类的单个元素。这种现有技术方法的特征在于以下限制和缺点：
1.在不根据对于参考序列的映射结果将读序分类至单独的数据类的情况下对读序或读序对进行编码并且使用唯一的描述符超集是低效的方法，会产生较差的压缩性能。
2.将读序对作为单独的读序进行编码需要复制携带相同信息的若干描述符，因此造成低效且产生较差的压缩性能。
3.取得重构读序对所需的信息是复杂且低效的，因为该处理可能需要在整个数据集中进行强力顺序搜索，而在下一代测序(NGS)技术的情况下，整个数据集可能是非常大的。
4.对映射到特定基因组区域的读序或读序对的选择性访问需要搜索整个数据集，以保证取得所有的读序或读序对。
当通过单个描述符子集对读序对进行编码时，以下技术优点对于本领域技术人员是显然的：
1.不把读序对作为单个元素(例如，读序对标识符、映射距离、映射参考标识符、当前由SAM文件格式中的特定标志编码的各种映射质量信息)进行编码来复制显然是冗余的两个读序共有的信息
2.相互配对信息(即，规定哪个读序是任意当前读序的配对读序的信息)的取得是直接的并且不需要任何进一步的处理。相反，在现有技术的方法中，可能需要解析整个数据。
为了能够高效地选择性地访问测序数据的特定部分并能够在数字数据网络上传输它们，用于表示比对至参考序列的读序的描述符集被构造为在逻辑上分离且独立的数据块(称为访问单元(AU))。各个访问单元仅包含单个数据类的压缩表示，并且可以独立于任何其它访问单元或仅使用携带用于映射的参考序列区域的编码表示的访问单元进行解码。这实现了选择性访问和乱序传输功能。
为了提高压缩效率，本发明使得无需为两个读序都映射在同一参考序列上的各个读序对指定“映射参考标识符”描述符。每个访问单元只能包含映射在同一参考序列上的读序或读序对。通过使用这样的解决方案，表示参考序列标识符的描述符仅需要针对每个访问单元或访问单元组编码一次(并且不像目前在SAM/BAM格式中所做的那样针对每个读序进行重复)。
上述规则的唯一例外是两个读序映射在不同参考序列(例如染色体)上的读序对。在这种情况下，对读序对进行分割并且把两个读序作为两个单独的基因组记录进行编码，且每个编码读序包含其配对读序所映射到的参考序列的标识符。
本领域技术人员知道，对具有同种统计特性的元素组中的信息进行分类相对于使用应用于异种数据集的通用压缩器(例如，LZ型算法)提供了更优的压缩性能。因此，当通过特定的描述符子集成对地对基因组序列进行编码时，由于各个单独的描述符子集所特有的低熵以及重构和提取读序对时的高处理效率，实现了更好的压缩。
读序描述符
本节介绍指定用于表示映射在参考序列上的基因组读序的描述符。用于表示各个读序或读序对的特定描述符子集取决于相对于参考序列的映射结果(即，读序或读序对与参考序列之间存在或不存在失配)。
位置
读序或读序对位置被定义为读序或读序对中的最左侧碱基在参考序列上的映射位置。每个读序或读序对需要一个“位置”类型的描述符。“位置”描述符的值表示：
·读序或读序对中的最左侧碱基在参考序列上的坐标的值
·或者相对于编码在同一访问单元中的先前读序或读序对的坐标的差。
需要“位置”描述符来表示各个编码读序或读序对。
在本发明公开中，这种描述符将被称为pos描述符。
读序配对
在读序对的情况下，表示每个读序如何与其配对读序在读序对中相关联的描述符可以由若干语法元素表示，例如它可以表示：
·读序的碱基与配对读序中的相应碱基(例如，读序中的最左侧映射碱基与配对读序中的最左侧映射碱基)之间的坐标差。在本发明公开中，这种描述符将被称为pair描述符。
·配对读序在参考序列上的绝对坐标与配对序列所映射到的参考序列的标识符。在以下情况下使用此表示选项：
ο读序对的两个读序映射在不同的参考序列(例如，染色体)上，或者
ο当读序对的两个读序映射在同一参考序列上但是被超过给定值(该给定值被指定为最大可能的允许配对距离)的多个碱基分隔开时。
在本发明公开中，这种描述符将被称为abspair描述符。在配对读序映射在不同参考序列上的情况下，标识参考序列的描述符被称为refid。
·在成对读序的情况下，将读序与其配对读序分隔开的编码读序的数量。在本发明公开中，这种描述符将被称为pcount描述符。
读序长度
在具有可变长度的读序的情况下，每个读序的描述符用于表示被表达为构成该读序的核苷酸的数量的长度。显然，在可变读序长度的情况下，每个读序仅需要一个读序长度描述符。
在本公开中，这个描述符也称为rlen描述符。
反向互补
DNA由双螺旋构成，其中每条链是另一条链的互补，因为腺嘌呤('A')仅与胸腺嘧啶('T')结合且胞嘧啶('C')仅与鸟嘌呤('G')结合。因此，仅需要表示一条链便可获知另一条链的核苷酸构成。这就是参考序列总是由单个序列表示并且映射工具能够找到属于两个链的读序的映射位置的原因。如果读序被映射在DNA螺旋的互补链上，则称其为“反向互补”。描述符有必要携带这样的信息并且携带表明原始读序是否是它所映射到的参考序列的反向互补的信息。
每个读序需要反向互补描述符。
在本公开中，这种描述符也称为rcomp描述符。
未知碱基的位置
在测序处理中，可能发生机器无法在读序或构成读序对的读序的给定位置处检测到任何碱基的情况。这种事件由相应读序位置处的特殊符号“N”标识。因此需要辨别“N”符号在读序中的每一次出现的描述符。
这种描述符可以表示：
·表示为参考序列的坐标的、‘N’符号在读序中或在读序对的读序中的绝对位置，或者
·前一个‘N’在同一读序或读序对的读序中的相对位置。
在本公开中，这种描述符也被称为nmis描述符。
失配位置和类型
映射在参考序列上的读序可能相对于所映射到的参考序列片段具有失配。根据以下情况，这些失配可以被分类并标注为替换、缺失或插入：
·相对于参考序列存在不同的核苷酸(碱基)(替换)
·映射读序中缺少核苷酸(缺失)
·读序中存在参考序列中所不存在的核苷酸(插入)
每种失配类型的表示意味着使用三个描述符，一个表示读序或读序对中的失配位置(也称为mmpos)，另一个表示仅存在替换时的失配类型(也称为subtype)，另一个表示存在替换、插入和缺失时的失配类型(也称为mmtype)。
软剪辑
映射到参考序列上的基因组读序可能在它们的边缘部分存在与参考序列在映射位置处具有的核苷酸序列不匹配或者稍微匹配的核苷酸序列。这些序列部分被称为软剪辑，并且可以由如下描述符来表示：该描述符只是由表示序列部分的碱基的符号串构成。
读序只能在读序的起始处和/或终止处允许一个或两个软剪辑。
在本文中，这种描述符还被称为sclips描述符。
映射标志
映射标志用于携带如下的与比对处理相关的特定信息：
·读序或读序对的多个映射位置的存在
·拼接读序的存在
·PCR(聚合酶链反应)或光学重复(optical duplicate)的存在
·补加比对(在比对器针对同一读序或读序对产生了几个可能的映射位置的情况下使用)
·读序未通过质量检查(即，特定于技术供应商的用于测量测序处理的质量的手续)
在本文中，这种描述符也被称为flags描述符。
未映射读序
在读序未被映射到参考序列的任何位置的情况下，该读序被分类为未映射。在这种情况下，根据一些共有特性将所有的未映射读序分组。这种处理被称为“聚类”。具有相同特性的读序组被称为群集。在本发明公开中，在属于同一群集的读序之间共享的特性被称为群集签名或签名。
签名可以由从2到数千的任意数目的核苷酸构成，并且签名可以针对所有群集具有固定长度或可以具有可变长度。可以属于签名的符号的字母表取决于已被测序以产生正被处理的读序的特定基因组样本。作为示例而非限制，可以使用以下字母表：
·对于DNA
ο{A,G,C,T,N}
ο{A,G,C,T,R,Y,S,W,K,M,B,D,H,V,N,.,-}(IUPAC符号)
·对于RNA
ο{A,G,C,U,N}
ο{A,G,C,U,R,Y,S,W,K,M,B,D,H,V,N,.,-}(IUPAC符号)
·对于氨基酸
ο{A,C,D,E,F,G,H,I,K,L,M,N,P,Q,R,S,T,V,W,Y}
用于计算群集签名的字母表的类型由在本公开中描述的称为基因组数据集报头的数据结构所携带的参数Alphabet_ID来标识。
属于同一基因组数据集的群集的签名可以是固定长度或可变长度的。编码在基因组数据集报头中的全局参数用于通报签名长度是固定的还是可变的。如果签名长度是固定的，则第二全局参数表示集群签名的符号长度。在可变签名长度的情况下，此值为0。
每个未映射读序的集群被编码在一个或更多个访问单元中。
对不具有失配的读序或读序对进行编码
在读序或读序对完全地映射到参考序列上(即，没有任何失配)的情况下，需要以下描述符子集来重构该读序和相关联的映射信息：
·每个读序或读序对的一个位置描述符(pos)
·每个读序的一个，或每个读序对的两个反向互补描述符(rcomp)
·每个读序的一个长度描述符(rlen)(仅在可变长度读序的情况下)
·每个读序对的一个配对描述符(pair)
·映射标志(flags)集
在本发明中，这种读序或读序对被分类为属于P类。
位置描述符pos表示读序或读序对的最左侧映射碱基在参考基因组上的位置。它的使用如图5所示。
反向互补描述符rcomp指示读序是映射到参考序列的正链上还是负链上。这个描述符的含义和语法如图7和图8所示。
在可变长度读序的情况下，描述符rlen对读序长度进行编码。
pair描述符携带重构整个读序对所需的信息。图9提供了描述符的语法。
图10提供了对属于P类的读序对进行编码的示例，图11提供了可变长度的单读序的对应示例。
对具有仅由未知碱基表示的失配的读序或读序对进行编码
在读序或读序对映射到参考序列上但包含至少一个未知碱基的情况下，需要以下描述符子集来重构该读序和相关联的映射信息：
·每个读序或读序对的一个位置描述符(pos)
·每个读序的一个，或每个读序对的两个反向互补描述符(rcomp)
·每个未知碱基的位置(nmis)
·每个读序的长度描述符(rlen)(仅在可变长度读序的情况下)
·每个读序对的配对描述符(pair)
·映射标志(flags)集
已经存在于P类子集中的描述符具有相同的语法和语义。附加描述符nmis提供由测序处理检测为“未知”的碱基(符号'N')在读序(读序对)中的位置。
在本发明中，这种读序或读序对被分类为属于N类。
图12提供了对N类读序对进行编码的示例。
对具有未知碱基和替换的读序或读序对进行编码
在读序或读序对映射在参考序列上并且表现出至少一个替换但没有缺失或插入的情况下，需要以下描述符子集来重构该读序和相关联的映射信息：
·每个读序或读序对的位置描述符(pos)
·每个读序的一个，或每个读序对的两个反向互补描述符(rcomp)
·每个替换位置的描述符(mmpos)
·每个替换类型的描述符(subtype)
·每个未知碱基的位置(nmis)
·每个读序的长度描述符(rlen)(仅在可变长度读序的情况下)
·每个读序对的配对描述符(pair)
·映射标志(flags)集
已经存在于P类子集中的描述符具有相同的语法和语义。用于这种读序数据类的附加描述符是表示替换位置的mmpos和表示替换类型的subtype。
图13提供了对这种类型的读序对进行编码的示例。
在本发明公开内容中，这样的读序或读序对被认为属于M类。
对具有至少一个插入、缺失或软剪辑的读序或读序对进行编码
当读序或读序对在具有至少一个插入、缺失或软剪辑的情况下映射在参考序列上时，定义以下描述符子集：
·每个读序或读序对的位置描述符(pos)
·每个读序的一个，或每个读序对的两个反向互补描述符(rcomp)
·失配(插入、缺失、替换)的各个位置的描述符(mmpos)
·失配(插入、缺失、替换)的各种类型的描述符(mmtype)
·每个软剪辑序列的描述符(sclips)
·每个未知碱基的位置(nmis)
·每个读序的长度描述符(rlen)(仅在可变长度读序的情况下)
·每个读序对的配对描述符(pair)
·映射标志(flags)集
已经存在于M类子集中的描述符具有相同的语法和语义。在这种情况下使用的附加描述符是表示替换、插入和缺失的位置的mmpos、表示失配(替换、插入或缺失)的类型的mmtype以及表示软剪辑碱基的sclips。
在本发明公开内容中，这样的读序或读序对被认为属于I类。图14提供了对属于I类的读序对进行编码的示例。
仅一个读序被映射到参考序列的读序对
在读序对由映射读序(属于P类、N类、M类或I类中的任一个)和未映射读序构成的情况下，该读序对被分类为属于被称为HM(半映射)类的单独类。
映射到参考序列上的读序可以是上述任何类(P、N、M和I)，并且将使用已经针对每个类描述的描述符子集来编码。将通过使用恰当的熵编码器压缩表示未映射读序的符号的串来对其进行编码。
对未映射读序或读序对进行编码并构造“内部”参考
属于U类的读序或者属于HM类的读序对的未映射配对读序不能映射到满足指定的匹配精确度约束集的任何“已有”参考序列。本发明公开了一种构造要用于比对和压缩属于这些数据类的读序的一个或更多个“内部”参考序列的方法。
有几种方法可以构造恰当的“内部”参考，例如但不限于：
·将未映射读序分割为包含如下读序的群集：该读序共有至少最小尺寸的共同连续基因组序列(签名)。每个群集可以由其签名唯一地标识，如图24所示。
·按任何有意义的顺序(例如，字典式顺序)对读序进行排序，并使用最后N个读序作为对第N+1个读序进行编码的“内部”参考。这种方法如图25所示。
·对U类读序的子集执行所谓的“重新组装(de-novo assembly)”，以便能够根据指定的匹配精度约束或新的约束集来比对和编码所有或相关的属于所述类的读序的子集。
如果正被编码的读序可被映射到满足指定的匹配精度约束集的“内部”参考上，则使用可为以下类型的语法元素来编码在压缩之后重构读序所需的信息：
1.就内部参考中的读序编号而言，匹配部分在内部参考上的起始位置(pos描述符)。这个位置可以被编码为绝对值或相对于先前编码的读序的差分值。
2.起始位置相对于内部参考的相应读序的起点的偏移(pair描述符)。例如，在固定读序长度的情况下，实际位置是pos*length+pair。
3.被编码为失配位置(snpp描述符)和失配类型(snpt描述符)的可能存在的失配
4.读序的与内部参考不匹配(或者与内部参考匹配，但是有超过规定阈值的数量的失配)的那些部分(通常是被读序对辨别的边缘)被编码在indc描述符中。如图26所示，可以对所使用的内部参考的一部分的边缘执行填充操作，以减小编码在indc描述符中的失配的熵。编码器可以根据正在处理的基因组数据的统计特性来选择最恰当的填充策略。可能的填充策略包括：
a.无填充
b.根据在当前编码数据中的频度选择的固定填充模式
c.与依据最近N个编码读序限定的当前场景的统计特性相应的可变填充模式
填充策略的具体类型将由indc描述符流报头中的特殊值通报。
5.指示是否已使用内部自我生成的参考、外部参考或不使用参考对读序进行了编码的标志(rtype描述符)
6.逐字地被编码的读序(ureads描述符)。
图26提供了这种编码过程的示例。
图27示出了内部参考上的未映射读序的替代编码，其中，pos+pair语法元素由带符号的pos取代。在这种情况下，pos将按照参考序列上的位置表示读序n的最左侧核苷酸位置相对于读序n-1的最左侧核苷酸位置的距离。
在U类读序具有可变长度的情况下，使用附加描述符rlen来存储各个读取长度。
这种编码方法可被扩展，以对于每个读序支持N个起始位置，使得可在两个以上参考位置处分割读序。这对编码由产生非常长的读序(50K+个碱基)的那些测序技术(例如，来自Pacific Bioscience)产生的读序(通常存在由于测序方法中的循环产生的重复模式)特别有用。同样的方法也可用于编码嵌合读序，该嵌合读序被定义为在没有或几乎没有重叠的情况下比对至基因组的两个不同部分的读序。
上述方法显然不仅可以应用于U类，而且可以应用于包含与读序位置相关的语法元素的任意流(pos流)。
具有多个映射位置的读序和读序对
在诸如SAM或CRAM的现有技术方法中，通过复制编码数据记录以及临时追加可选字段来支持读序或读序对被映射到参考序列的多个坐标的情况，这会导致数据表示中的分散，从而由于引入的冗余而造成压缩效率的明显损失。此外，双端读序的一些映射配置(其中一个读序配对到相应配对读序的多个映射位置)仅由用户定义的可选字段支持。对于需要被分割成两个或更多个子序列以便找到对参考序列的恰当映射的读序和读序对的表示也可以如此。这种方法的主要结果是在压缩效率方面导致低效表示的激增，对于数据可重用性也有很大的影响，并且在信息从一个部分表示传递到另一个部分表示时会有损失。
本公开描述的方法支持多个映射位置和拼接读序两者的表示，并且在压缩和压缩数据可访问性方面都比现有解决方案更高效。由于能够将具有同种统计特性的描述符分组，并且有着具有恰当场景的不同熵编码器，提供了更优的压缩。通过定义这样的索引机制：该索引机制使得能够仅解压并提取特定类型的基因组数据而无需解压并访问整个编码信息，提供了更好的文件可访问性。例如，本公开中描述的发明提供了如下的索引机制：该索引机制使得能够提取带有或不带有拼接读序、具有多个映射位置的读序或读序对。这在现有技术的基因组信息表示格式(如SAM和CRAM)中是不可能的。
当对具有多重比对的读序对进行编码时，SAM和CRAM不能支持读序对中的一个读序的与另一读序的单个比对相关联的多重比对的表示。本领域技术人员知道，这在诸如RNA-seq和ChIP-seq的实验中是常见的情况，并且当今SAM和CRAM没有办法支持读序对中的读序的多重比对之间的结合的所有可能组合。本公开中提出的解决方案能够支持目前在高通量测序(HTS)机生成的基因组数据中发现的所有可能的配置。
当将基因组读序映射到参考序列时，以下结果是可能的：
1.读序根据建立的约束映射到参考序列的单个区域。在这种情况下，映射被认为是“唯一的”；
2.读序与参考序列的一个以上的基因组区域完全匹配；
3.比对工具报告多个可能的映射位置。每个映射表现出相对于相应区域的不同编辑距离，但是比对器认为它们都是有效的。比对器通常以完全依赖于实现的方式为每个比对分配分数。根据分配的分数，通常标记一个“主”比对。比对可以表现出与一些其它比对完全相同的编辑距离和相同的失配。
情况1仅要求对比对的唯一性进行编码。
情况2需要确认不存在“主”或“辅”比对，因为所有的映射都是同样可能的。要编码的唯一附加信息是所有映射位置的矢量。
情况3要求把所有映射位置作为编码读序的阵列进行编码。当需要时，须针对每个映射位置复制本发明中公开的所有所需描述符。不同的映射位置可以表现相对于参考的不同程度的误差(替换、索引、剪辑碱基)。
拼接读序
拼接读序(spliced read)被定义为如下的读序：该读序需要被分割成两个或更多个子读序，以便在参考序列上找到合适的映射区域。在这种情况下，子读序之间的距离(称为“接界(junction)”)通常太大而不能被认为是缺失。拼接读序的映射可以涉及到参考序列的正链或负链，因此必须保留和编码被称为“绝对链”的这种信息。
对拼接读序进行编码需要保留每个拼接段的映射位置，该拼接段必须被认为是可变长度读序、编码在仅包含固定长度读序的基因组数据集中。
多重比对和拼接读序的描述符
在以下说明中，如SAM规范中那样使用术语模板，该模板辨别在测序机上测序了一部分或者由原始序列组装而成的核苷酸序列。根据所使用的测序技术，模板的测序可以产生单个核苷酸序列(一个读序)或被认为是“成对的”的两个序列。在本语境中，片段被定义为连续序列或子序列。
当对具有多个比对位置的读序或读序对进行编码时，要由编码数据传送的信息如下：
·每个片段在模板(例如，读序1和读序2)中的所有映射位置
·模板(例如，读序1和读序2)的任何主映射的存在
在本发明公开中，对多个映射位置和拼接读序进行编码是由以下项支持的：
1.在编码数据集的级别上定义的两个全局标志，这两个全局标志对构成数据集的所有访问单元有效：
·ma：当设定时，它指示数据集中存在多重比对
·sr：当设定时，它指示数据集中存在拼接读序
2.以下描述符：
·mmap描述符用于通报在多少个位置处对读序或读序对中的最左侧读序进行了比对；第一数目(N)指的是作为单个片段的读序(如果sr标志＝0)，或者指的是为实现多个可能的比对而将读序拆分成的所有片段。N的值表示针对这个记录中的模板编码了多少pos描述符的值。N之后是一个或更多个数字Mi(如下面针对没有拼接段和具有拼接段的情况所描述的)。
·splen描述符用于拼接比对的情况。它是通报有多少个碱基构成了当前记录的每个拼接段的长度列表。由于每个读序的长度是已知的(在固定读序长度的情况下根据全局读序长度得知，或者在可变长度读序的情况下根据len描述符得知)，所以解码器使用该描述符来计算每个比对使用了多少片段。
·mmscore描述符提供每个比对的分数值和主比对的指示。分数提供由比对工具分配给每个特定映射配置的置信度水平。它通常是浮点值。例如，主比对(即，被认为是最佳候选的比对)可以是得分最高的比对。
·mmsc描述符支持不保持相同的主比对映射毗连性(contiguity)的替代辅比对(最常见的示例是CIGAR字符串包含indel而主比对仅包含M的情况)。在这种情况下，有必要对主比对和其它比对之间的毗连差异进行编码。
不具有拼接段的多重比对
如果数据集中不存在拼接段，则全局sr标志不被置位，并且不使用splen描述符。
在双端测序中，mmap描述符被构造成一个值N，在N之后的是一个或更多个数字(Mi)，其中i从1到整个第一(最左侧)读序比对的数目(N1)。如果针对该读序不存在拼接，则N1＝N。对于每个第一读序比对，无论拼接与否，都使用一个值Mi来通报使用了多少片段来比对第二读序(在这种情况下等于比对数目)，然后是对于第一读序的该比对编码了pair描述符的多少个值。
Mi的值将被用于计算
这种方法在图17中示出。
Mi＝0的特定值表示最左侧读序的第i个比对与最右侧读序的一个比对成对，最右侧读序的这个比对已经与最左侧读序的第k个比对成对，其中k<i(于是检测不到新的比对，这与上式一致)。
例如，在最简单的情况下：
1.如果最左侧读序存在单个对比且最右侧读序存在两个替代比对，则N取值1且M1取值2。
2.如果对于最左侧读序检测到两个替代比对，但对于最右侧读序仅检测到一个比对，则N取值2，M1取值1且M2取值0。
当Mi等于0时，pair描述符的关联值指向已有的第二读序比对。
图17例示了在没有拼接段的多重比对的情况下N、P和Mi的含义，图18示出了如何使用pos、pair和mmap描述符来编码多重比对信息。
具有拼接段的多重比对
在拼接读序的情况下，使用splen描述符。它对于各个记录由N+P个值的阵列构成。前N个值表示单读序或读序对的第一(例如，最左侧)读序的每个比对片段的长度。在双端测序的情况下，以下P值表示第二读序的每个比对片段的长度。像计算
记录的splen描述符的前N个值使得能够计算N1，其表示第一读序的比对数目。如果N1＝N，则第一读序不存在拼接。
记录的splen描述符的接下来的P值使得能够计算N2，N2通报第二读序的比对数目。如果N2＝P，则第二读序不存在拼接。
所定义的mmap和splen描述符使得能够唯一地辨别有多少读序或读序对表示多个映射以及有多少片段构成每个读序或读序对映射。这在图19和图20中示出。
比对得分
mmscore描述符允许通报比对的映射分数。在单端测序中，每个模板具有N1个值；在双端测序中，整个模板的每个比对都有一个值。换句话说，由读序1的一个比对和读序2的一个比对构成的每一对可以具有相关联的分数。在双端读序的情况下，总分数的数目被计算为
N.分数＝MAX(N1,N2)+M0
其中，N1是读序1的总比对数，N2是读序2的总比对数，M0是等于0的Mi值数。
在单端读序的情况下，分数的数目等于N1
不具有拼接描述符的多重比对
下表总结了在不具有拼接读序的多重比对的情况下，在本发明公开内容中定义的mmap和mmscore描述符的使用的语义和作用。
具有拼接描述符的多重比对
下表总结了在具有拼接读序的多重比对的情况下，在本发明公开内容中定义的mmap和mmscore描述符的使用的语义和作用。
不同染色体上的多重比对
可能发生如下情况：比对处理找到了对于与主映射所在的参考序列不同的另一参考序列的替代映射。
在这种情况下，就应用而言，最重要的是在两个或更多个必然不同的访问单元之间维持快速(就随机访问复杂性而言)链接，在该访问单元中对模板的若干替代映射进行了编码。
对于唯一比对的读序对，本发明公开定义了用于表示嵌合比对(其中，读序对中的两个读序被映射在不同的染色体上)的命名为pair的描述符。这种描述符可以用于通报包含同一模板的其它比对的下一记录的参考和位置。这在图30中示出。
在读序对的最左侧读序的一个或更多个比对存在于与当前编码AU相关的参考序列不同的参考序列上的情况下，应使用pair描述符的保留值(与在唯一比对的情况下对于另一参考的比对所使用的值不同)。保留值后面应跟随下一个AU中包含的所有比对中最左侧比对的参考和位置(即，该记录的pos描述符的第一个解码值)。
具有插入、缺失、未映射部分的多重比对
在一些情况下，多重比对可以表现出匹配碱基和失配碱基、插入、缺失和软剪辑的不同配置。例如，虽然主比对仅具有匹配碱基或失配碱基并因此被映射为核苷酸的连续序列，但辅比对可以表现出插入、缺失、软剪辑或拼接，这并不罕见。在本公开中，当映射读序不包含插入、缺失或软剪辑碱基时，认为具有映射毗连性。本发明公开定义了第三描述符，该第三描述符通报辅比对是属于保留了主比对的相同映射毗连性的数据类(就像是P、M或N)还是属于不保留主比对的相同映射毗连性的数据类(I、U、拼接)。这种描述符mmsc(用于多个映射子类)原则上仅是对于每个比对的标志。当设定mmsc时，其之后是由比对器生成的SAM cigar字符串的逐字表示，其扩展有表示未映射核苷酸的附加符号“U”，并且之后是未映射核苷酸的串。图31例示了使用这种描述符的示例。
这种描述符的语法如下：
·如果所有比对共有相同的映射毗连性，则存在单个值并将其设置为0(N＝0)。
·如果至少一个辅比对不保持映射毗连性，则描述符构造如下：
ο第一值N>0表示有多少辅比对不保留映射毗连性。然后，对于不保留映射毗连性的各个辅比对，重复N次以下元素：
ο一个值表示与不保留映射毗连性的辅比对相关的配对值。
ο接下来的值表示该读序对中的哪个读序不保留映射毗连性
ο接下来的值包含表示不保留映射毗连性的读序的映射特性的cigar字符串
ο仅当cigar字符串包含符号“S”时才存在的一个可选字段包含软剪辑核苷酸的逐字字符串。
参考序列描述符
参考序列通常被表示为表示可以在相应生物样本中找到的核苷酸的符号串。在DNA的情况下，核苷酸有四种并且由符号A、C、G和T表示。在RNA的情况下，T由U代替。添加第五符号以表示序列中的测序装置不能根据实验所需的置信度确定核苷酸类型的坐标。在本发明公开中，参考序列可以整体地编码在一个访问单元中，或者分割成两个或更多个子序列。
用于表示待进行熵编码的参考序列或子序列的描述符是依照相应字母表的允许符号的序列或子序列的逐字表示。
源模型、熵编码器和编码模式
对于本发明公开的基因组数据结构的每个数据类、子类和相关描述符流，可以根据每个流携带的数据或元数据的具体特征及其统计特性采用不同的编码算法。“编码算法”必须作为描述符流的特定“源类型”与特定“熵编码器”的关联。可以指定和选择特定的“源模型”以在源熵的最小化方面获得数据的最高效编码。可以考虑编码效率和/或概率分布特征以及相关的实现问题来选择熵编码器。特定“编码算法”(也称为“编码模式”)的每个选择可以应用于与整个数据集的数据类或子类相关联的整个“描述符流”，或者可以对分割为访问单元的描述符的每个部分应用不同的“编码模式”。
与编码模式相关联的每个“源模型”的特征在于：
·由每个源(即，用于对诸如读序位置、读序配对信息、相对于表2中定义的参考序列的失配的数据类进行表示的描述符集)发出的语法元素的定义。
·关联概率模型的定义。
·关联熵编码器的定义。
其它优点
将序列数据分类成所定义的数据类和子类使得能够充分利用较低的信源熵(特征在于由单个单独的数据源(例如，距离、位置等)对语法元素的序列进行建模)来实现高效的编码模式。
本发明的另一优点是可以仅访问所关注的数据类型的子集。例如，基因组学中最重要的应用之一在于找到基因组样本相对于参考(SNV)或群体(SNP)的差异。当今，这种类型的分析需要处理完整的读序，而通过采用本发明公开的数据表示，失配已被隔离到仅一至三个数据类(取决于是否有意考虑“n型”和“i型”失配)。
另一优点是当发布新的参考序列时或者当对已经映射的数据执行重新映射(例如，使用不同的映射算法)以获得新的比对时，可以把参照特定"外部"参考序列压缩的数据和元数据高效地转码为参照另一不同的"外部"参考序列。
图28示出了根据本发明原理的编码装置287。编码装置287接收例如由基因组测序装置280产生的原始序列数据289作为输入。基因组测序装置280在本领域中是已知的，如Illumina HiSeq 2500、Thermo-Fisher Ion Torrent设备或Oxford Nanopore MinION。将原始序列数据289馈送到比对器单元281，该比对器单元281通过将读序与参考序列2820进行比对来准备用于编码的序列。另选地，可以利用专用模块282，通过使用不同策略(如在本文的“构造用于U类和HM类未映射读序的内部参考”章节中所述的)由可用读序生成参考序列。在由参考生成器282处理之后，可以将读序映射到所获得的更长序列上。然后由数据分类模块284对比对后的序列进行分类。然后，数据类288被馈送到描述符编码器285-287。然后将基因组描述符流2811馈送到算术编码器2812-2814，算术编码器2812-2814根据层所携带的数据或元数据的统计特性来对层进行编码。结果得到一个或更多个基因组流2815。
图29示出了根据本发明的原理的解码装置298。解码装置298从网络或存储元件接收复用的基因组比特流2910。将复用的基因组比特流2910馈送到解复用器290，以产生单独的基因组比特流291，然后将该基因组比特流291馈送到熵解码器292-294，以产生基因组描述符流295。提取的基因组描述符流被馈送到描述符解码器296-297，以进一步将描述符解码成读序类。类解码器299进一步处理基因组描述符2911和变换后的参考2914，并合并结果以产生未压缩的读序，然后该读序可以进一步以本领域已知的格式(例如文本文件或zip压缩文件，或FASTQ或SAM/BAM文件)进行存储。
类解码器299能够通过操控由一个或更多个基因组比特流携带的关于原始参考序列的信息来重构该原始基因组序列。在参考序列不通过基因组流传输的情况下，它们必须在解码侧可用并且可由类解码器访问。
本文公开的发明技术可以用硬件、软件、固件或其任意组合来实现。当以软件实现时，这些可以存储在计算机介质上并由硬件处理单元执行。硬件处理单元可以包括一个或更多个处理器、数字信号处理器、通用微处理器、专用集成电路或其它离散逻辑电路。
本公开的技术可以在各种设备或装置中实现，包括移动电话、台式计算机、服务器、平板电脑和类似设备。
访问单元类型
在本发明公开中，被分类在数据类中并被构造在压缩或未压缩层中的基因组数据被组织成如上文所定义的不同的访问单元。
访问单元由以下项目区分：
·类型，其表征访问单元携带的基因组数据和数据集的性质以及访问该访问单元的方式。
·顺序，针对属于相同类型的访问单元提供唯一的顺序。
可以将任何类型的访问单元进一步分类为不同的“类别”。
以下是不同类型的基因组访问单元的定义的非穷尽列表：
1)类型0的访问单元不需要参照任何来自其它访问单元的信息来进行访问或者解码并访问。它们所包含的数据或数据集携带的全部信息可以由解码装置或处理应用独立地读取和处理。作为示例但不作为限制，类型0的访问单元可以用于携带编码参考序列(例如染色体或整个参考基因组或其部分)。
2)类型1的访问单元包含参照由类型0的访问单元携带的数据的数据。对包含在类型1的访问单元中的数据进行读取或解码和处理需要对一个或更多个类型0的访问单元进行访问。类型1的访问单元编码了与“P类”读序相关的基因组数据
3)类型2的访问单元包含参照由类型0的访问单元携带的数据的数据。对包含在类型2的访问单元中的数据进行读取或解码和处理需要对一个或更多个类型0的访问单元进行访问。类型2的访问单元编码了与“N类”读序相关的基因组数据
4)类型3的访问单元包含参照由类型0的访问单元携带的数据的数据。对包含在类型3的访问单元中的数据进行读取或解码和处理需要对一个或更多个类型0的访问单元进行访问。类型3的访问单元编码了与“N类”读序相关的基因组数据
5)类型4的访问单元包含参照由类型0的访问单元携带的数据的数据。对包含在类型4的访问单元中的数据进行读取或解码和处理需要对一个或更多个类型0的访问单元进行访问。类型4的访问单元编码了与“I类”读序相关的基因组数据
6)类型5的访问单元包含无法映射到任何可用参考序列上并使用内部构造的参考序列进行了编码的读序(“U类”)。类型5的访问单元包含参照由类型0的访问单元携带的数据的数据。对包含在类型5的访问单元中的数据进行读取或解码和处理需要对一个或更多个类型0的访问单元进行访问。
7)类型6的访问单元包含读序对，在该读序对中，一个读序可以属于P、N、M、I四个类中的任一者，且另一读序不能映射到任何可用的参考序列上(“HM类”)。类型6的访问单元包含参照由类型0的访问单元携带的数据的数据。对包含在类型6的访问单元中的数据进行读取或解码和处理需要对一个或更多个类型0的访问单元进行访问。
8)类型7的访问单元包含元数据(例如，质量分数)和/或与包含在类型1的访问单元中的数据或数据集相关联的注释数据。可以对类型7的访问单元进行分类并标记在不同层中。
9)类型8的访问单元包含被分类为注释数据的数据或数据集。类型8的访问单元可以按层分类和标记。
10)附加类型的访问单元可以扩展本文描述的结构和机制。作为示例而非限制，基因组变异检测、结构和功能分析的结果可以被编码在新类型的访问单元中。本文描述的访问单元中的数据组织不妨碍把任何类型的数据封装在访问单元中，是对于编码数据的性质而言完全透明的机制。
对类型0的访问单元进行排序(例如，编号)，但是它们不需要以有序的方式存储和/或传输(技术优点：并行处理/并行流传输、复用)。
类型1、2、3、4、5和6的访问单元不需要排序，并且不需要以有序的方式存储和/或传输(技术优点：并行处理/并行流传输)。
图22示出了访问单元是如何由报头以及一个或更多个同种数据流构成的。每个流可以由一个或更多个块构成。每个块包含几个分组，分组是以上介绍的用于表示例如读序位置、配对信息、反向互补信息、失配位置和类型等的描述符的结构化序列。
每个访问单元可以在各个块中具有不同数目的分组，但是在访问单元内，所有块具有相同数目的分组。
每个数据分组可以由3个标识符X、Y、Z的组合来辨别，其中：
·X辨别它所属的访问单元
·Y辨别其所属的块(即，它封装的数据类型)
·Z是表示相对于同一块中的其它分组的分组顺序的标识符
图23示出了访问单元和分组标记的示例，其中AU_T_N是具有标识符N的类型T的访问单元，根据访问单元类型，标识符N可以意味着或不意味着顺序的概念。标识符用于唯一地将一种类型的访问单元与完全解码所携带的基因组数据所需的其它类型的访问单元关联起来。
任何类型的访问单元都可以根据不同的测序处理被进一步分类和标记为不同的“类别”。例如但不作为限制，可以在以下时候进行分类和标记：
1.在不同时间对同一生物体进行测序(访问单元包含具有“时间”含义的基因组信息)
2.对相同生物体的不同性质的有机样本(例如，人体样本的皮肤、血液、头发)进行测序。这些是具有“生物”内涵的范围单元。
数据存储模式
当在存储介质上存储编码描述符时，本发明公开中描述了两种方法：
1.访问单元毗连(Access Unit Contiguous，AUC)模式
2.描述符流毗连(Descriptor Stream Contiguous，DSC)模式
当应用AUC时，属于同一访问单元(但属于不同描述符流)的编码数据块被存储在存储介质的毗连区域中。
AUC模式可以在数据存储方面以两种不同的方式实现：
a)类毗连(Class Contiguous，CC)，其中，邻接地存储在存储介质上的访问单元的顺序可以是基于类的，即，同一类(Class_ID)的所有访问单元被毗连地存储在存储介质中。
b)基因组区域毗连(Genomic Region Contiguous，GRC)，其中，访问单元的顺序是基于参考的，即，毗连地存储具有相同AU_ID(即，映射到同一基因组区域)的所有访问单元。
AUC/CC模式方法在访问单个类的数据时更高效。当访问映射到同一基因组区域的任意类的数据时，AUC/GRC模式更高效。本公开内容中描述的发明和相关语法支持所有模式DSC、AUC/CC和AUC/GRC方法，并且使得编码器能够自由地根据期望的选择性访问性能选择任何模式。表2中列出的基因组数据集报头携带的名为CC_Mode_Flag的标志通报是使用了AUC/CC模式还是AUC/GRC模式。
当应用DSC时，属于同一描述符流的块被存储在存储介质的毗连区域中。基因组数据实际上是按照描述符流(由一个或更多个块构成)来组织的，就熵编码而言表示同种数据。
在表2列出的基因组数据集报头中存储的被称为Block_Header_Flag的标志通报基因组数据集在编码方法之间使用的存储备选。
AUC模式和DSC模式之间的差异在图21中示出，其中，访问单元构造在坐标系中，该坐标系在垂直维度上具有对编码在数据块Bn_m中的描述符的类型进行辨别的Descriptor_ID并且在水平维度上具有访问单元ID。数据块Bn_m包含访问单元m的类型(即，标识符)n的编码描述符。
对压缩基因组数据进行索引以实现高效的选择性访问
为了支持对比对数据的特定区域的选择性访问，在本公开中描述了两种数据结构：携带全局参数的基因组数据集报头以及在编码和解码处理期间使用的称为主索引表(MIT)的索引工具。Genomic Dataset Header的语法在其中提供，Master索引表的语法在其中提供。
基因组数据集报头
基因组数据集报头是承载由编码器和解码器用来操纵编码基因组信息的全局参数的数据结构。
基因组数据集报头携带的信息包括：
·用于唯一地辨别每个数据集组的数据集组标识符
·用于唯一地辨别每个数据集的基因组数据集标识符
·用于辨别数据集遵循的数据格式规范的品牌标识符
·用于辨别数据集遵循的数据格式规范的次要版本号
·用于通报固定长度读序的以核苷酸为单位的编码基因组读序长度，
·通报存在双端读序的标志
·通报存在块报头的标志
·通报在数据集中使用了哪种AU编码模式(AUC/CC或AUC/GRC)的标志
·用于对读序相对于参考序列的失配进行编码的字母表类型，
·用来对数据集进行编码的参考序列的数量
·用于唯一地辨别每个参考序列的针对每个参考序列的数字标识符
·用于唯一地辨别每个参考序列的针对每个参考序列的字符串标识符
·用于对与各个参考序列相关联的访问单元进行计数的针对每个参考序列的编码访问单元数目
·用于区分比对读序、未比对读序、未映射读序和参考序列的编码基因组数据类型
·编码在数据集中的数据类的数目
·在解码处理时使用的数据集中所编码的针对每个数据类使用的描述符的数目
·用于索引编码未映射读序的群集总数
·用于表示整数值(用于对编码群集签名的解码所用的群集签名进行编码)的比特数，
·通报是否所有的群集签名就核苷酸数目而言具有相同长度的标志，
·群集签名的长度。
基因组数据集标题的每个元素的语法和语义在下面列出。
表2–基因组数据集报头语法
主索引表
本发明公开了一种称为主索引表(MIT)的索引工具。
主索引表(MIT)是基于多维阵列的数据结构，包括：
·包含在访问单元(来自不同描述符流的块的集合)中的所有读序或读序对的主比对中的最左侧匹配碱基相对于参考序列的位置(核苷酸的数目)。在表3中，这由Start_AU_Ref_Position[Sequence_ID][Class_ID][AU_ID]表示。
·包含在访问单元(来自不同描述符流的块的集合)中的所有读序或读序对的主比对中的最右侧匹配碱基相对于参考序列的位置(核苷酸的数目)。在表3中，这由End_AU_Ref_Position[Sequence_ID][Class_ID][AU_ID]表示。
·构成相对于各个参考序列编码的各个类的各个AU的各个编码描述符块的第一个字节的字节偏移。相对于Dataset Payload(基于0)的第一个字节计算偏移。如果块是空的并且(1)Block_Header_Flag被设定，则它等于0xFFFFFFFF。如果块是空的并且(2)Block_Header_Flag未被设定，则它等于描述符流中的下一块的Block_Byte_Offset值，或者对于描述符流中的最后块，等于Descriptor_Stream_Size[Class_ID]。在表3中，这由Block_Byte_Offset[Sequence_ID][Class_ID][AU_ID][Descriptor_ID]表示。
·根据所使用的的编码方法(由全局配置参数通报)，存在两个备择信息块：
ο如果各个访问单元作为毗连的数据块存储在存储介质上，则是以字节为单位的每个访问单元的大小，或者
ο如果同一类型的所有描述符被编码并作为毗连的数据块存储在存储介质上，则是编码描述符的每个块的大小。
MIT的最后一部分包含两个备择部分，这些备择部分根据描述符的各个编码块前面的报头的存在而使用。如果存在块报头(设定了Block_Header_Flag)，则MIT包含每个描述符流的以字节为单位的大小。如果不存在块报头(未设定Block_Header_Flag)，则MIT包含各个访问单元的以字节为单位的大小。
两种编码方法之间的选择由表2中的称为Block_Header_Flag的标志通报。
表3–主索引表
主索引表和多重比对
在存在多重比对的情况下，复制上述MIT以提供考虑到了编码成访问单元的读序或读序对的多重比对的索引工具。扩展主索引表包含
·包含在访问单元(来自不同描述符流的块的集合)中的所有读序或读序对的所有比对中的最左侧匹配碱基相对于参考序列的位置(核苷酸的数目)。在表3中，这由Start_AU_Ref_Position[Sequence_ID][Class_ID][AU_ID]表示。
·包含在访问单元(来自不同描述符流的块的集合)中的所有读序或读序对的所有比对中的最右侧匹配碱基相对于参考序列的位置(核苷酸的数目)。在表3中，这由End_AU_Ref_Position[Sequence_ID][Class_ID][AU_ID]表示。
