用于压缩基因组序列数据的方法与流程

1.本发明整体涉及表示由测序机器产生的基因组测序数据的方法，并且更具体地涉及用于压缩这样的基因组测序数据的计算机实现的方法。本公开提供了一种基于参考的压缩方法，其允许快速压缩和解压缩，同时不引起信息丢失，并且具有高压缩比。


背景技术：

2.下一代测序机器现在以可承受的价格产生海量的测序数据。最新的系统在36小时的单次运行中产生超过60亿个150个核苷酸长的序列，足以对20个完整的人类基因组进行测序。这为遗传疾病的诊断和个性化医疗的发展开辟了许多新的视角，旨在适应基于人基因组特异性的治疗。
3.然而，这也带来了新的挑战，特别是与存储海量数据相关的成本。用于原始(未比对)序列数据的最常用的文件格式为fastq格式，其保存序列数据(a、c、t、g核苷酸的字符串，也称为读段)、质量值(测序平台对每个核苷酸造成测序错误的概率)和序列名称。这是普通ascii文本文件，通常用通用文本压缩方案lz(lempel-ziv方案，在gzip软件中实现)来压缩。然而，使用此类压缩方法带来了几个问题：
[0004]-由于数据的冗余未被完全利用而导致的低压缩比
[0005]-压缩和解压缩缓慢
[0006]
还存在专用于fastq编码的压缩方法，分为基于参考的方法或基于非参考的方法。然而，它们中没有一个完全令人满意，原因是：a)基于参考的方法具有良好的压缩比，但速度缓慢，b)基于非参考的方法速度较快，但压缩比较低。这种基于非参考的方法的一个示例由软件spring提供，该软件是用于fastq文件的无参考压缩器(万维网地址：github.com/shubhamchandak94/spring)。然而，由软件spring提供的压缩方法的压缩比低。
[0007]
在基于参考的压缩方法中，已经提出了一些使用序列比对并且旨在具有更快的速度且具有良好压缩比的方法。然而，此类方法遭遇了几个问题，值得注意的一个主要问题是，它们并非完全无损。这种已知的基于参考的压缩方法例如描述于专利文献wo 2018/068829 a1中。在该所描述的方法中，在与一个或多个参考序列比对之后，根据匹配准确度将核苷酸序列分类(从而创建比对读段的类别)，然后使用不同的源模型和熵编码器针对其中数据被分区的每一层将这些核苷酸序列编码为多层语法元素。因此，数据的类别被单独编码，并且被构造在语法元素的不同层中，每个层均包含描述符，这些描述符单义地表示所述层的经分类和比对的读段。该方法旨在获得具有简化的信息熵的不同信息源，从而允许提高压缩性能以及对特定类别压缩数据的选择性访问。然而，这种压缩方法以与在读段比对步骤结束时获得的顺序不同的顺序对读段进行重新排序(即，根据其类别对读段进行重新排序)。于是在压缩过程(尤其是初始序列排序)中丢失了一些信息。因此，一些分析结果的再现性可能受到影响，因为一些下游分析软件可能依赖于读段的顺序。此外，以与读段的初始顺序不同的顺序解压缩数据使得检查未压缩文件是否与初始文件相同更加困难。此外，这种压缩方法速度相对缓慢，尤其是在与最先进的基于非参考的压缩方法相比时。


技术实现要素：

[0008]
以下独立权利要求的特征通过提供一种用于压缩基因组序列数据的方法解决了现有技术解决方案的问题。在一个方面，公开了用于压缩由测序机器产生的基因组序列数据的计算机实现的方法，所述基因组序列数据包括已与参考序列比对的核苷酸或碱基的序列的读段，从而产生比对读段，所述比对读段作为读段列表存储在初始文件中，所述方法包括以下步骤：
[0009]-针对每个比对读段，确定所述读段与所述参考序列是完全映射的还是不完全映射的，或者所述读段是否与所述参考序列是未映射的，
[0010]-根据所述确定对所述读段进行编码，其中根据第一编码过程对被确定为完全映射的所述读段进行编码，并且根据第二编码过程对被确定为未映射的所述读段进行编码，
[0011]-其中所述确定步骤包括针对每个不完全映射的读段，将所述读段与所述参考序列之间的失配的数目与阈值进行比较，
[0012]-其中，在所述编码步骤中，根据所述第二编码过程或第三编码过程对被确定为不完全映射的所述读段进行编码，当所述失配数目大于所述阈值时，根据所述第二编码过程对所述不完全映射的读段进行编码，当所述失配数目低于所述阈值时，根据所述第三编码过程对所述不完全映射的读段进行编码，
[0013]-其中，在所述第二编码过程中，所述读段的每个核苷酸或碱基被单独编码，
[0014]-其中所述第一编码过程和所述第三编码过程包括不同的描述符集合，每个描述符集合单义地表示与对应的编码过程相关联的读段，所述第一编码过程和所述第三编码过程中的每一者是简化的信息源熵编码过程。
[0015]
本发明克服了现有压缩方法的缺点，方式为：允许快速压缩和解压缩，同时不引起信息丢失，并且提供高压缩比。更具体地，本发明的重点是以最紧凑的方式对最频繁出现的情况进行编码，即使这意味着对于罕见的最不频繁出现的情况采用降级的编码模式。这导致压缩性能的巨大提高。此外，由于本发明中使用的基因组信息表示格式，由根据本发明的方法执行的压缩更快。最后但并非最不重要的是，根据本发明的方法照此保持读段的初始顺序，并且不根据读段的类别对读段进行重新排序。因此，在该过程期间没有信息丢失，这使得能够更容易地进行下游分析，以及在解压缩步骤之后进行有效的一致性检查。
[0016]
根据附图和以下具体实施方式，本发明的这些和其他特征和优点将变得更加明显。此外，尽管阈值在本文中可以被称为被超过或未被超过，但是应当理解，可以概念性地采用此类阈值，使得确定此类阈值是否被满足、符合或以其他方式被检测到，而不管用于实现那些阈值评估的数目或值是使用正值还是负值来描述。
[0017]
根据本公开的一个创新方面，公开了一种用于压缩基因组序列数据的方法。在一个方面，所述方法可以包括经由通过一个或多个计算机执行软件指令来进行一个或多个操作，其中所述操作包括：由所述一个或多个计算机获得读段记录；由所述一个或多个计算机确定所述读段记录是对应于被完全映射到参考序列的读段还是被不完全映射到所述参考序列的读段；基于由所述一个或多个计算机确定所述读段记录对应于被不完全映射到所述参考序列的读段，由所述一个或多个计算机确定所述不完全映射的读段的失配数目是否满足预先确定的失配阈值数目，以及基于确定所述失配数目满足所述预先确定的失配阈值数目，由所述一个或多个计算机将所述不完全映射的读段的每个失配编码为具有1个字节大
小的记录。
[0018]
其他方面包括执行如本文所公开的方法的动作的对应系统、装置和计算机程序，如由在计算机可读存储设备上编码的指令所定义的。
[0019]
这些和其他版本可任选地包括以下特征中的一个或多个特征。例如，在一些具体实施中，由所述一个或多个计算机确定所述不完全映射的读段的失配数目是否满足预先确定的失配阈值数目可以包括由所述一个或多个计算机确定所述不完全映射的读段的所述失配数目是否大于所述预先确定的失配阈值数目。
[0020]
在一些具体实施中，每个读段记录可以包括：指示比对读段相对于所述参考序列的绝对起始位置的数据；指示所述读段的长度的数据；指示所述读段是完全映射的还是不完全映射的数据；指示在所述读段中识别的失配的数目的数据；以及指示所述读段中的所述可能失配中的每一者的相对位置的数据。
[0021]
在一些具体实施中，将所述不完全映射的读段的每个失配编码为具有1个字节大小的记录包括：针对每个特定失配，由所述一个或多个计算机编码所述字节的前两位以包括表示存在于所述读段中的替代核苷酸或碱基而不是所述参考序列中的相应参考核苷酸或碱基的数据；以及由所述一个或多个计算机编码所述字节的剩余六位以包括表示所述参考序列中的所述失配的位置的数据，所述位置被计算为相对于所述读段的前一个失配的偏移。
[0022]
在一些具体实施中，所述方法可以还包括：由一个或多个计算机确定所述偏移是否大于最大可编码值；以及基于确定所述偏移大于所述最大编码值，由一个或多个计算机在所述特定失配与所述前一个失配之间插入至少一个假失配。
[0023]
在一些具体实施中，所述方法可以还包括：基于确定所述失配数目不满足所述预先确定的失配阈值数目，由一个或多个计算机使用简化的信息熵编码过程将对应于所述失配中的每一者的位置的所述参考序列的位置列表编码到所述参考序列中。
[0024]
在一些具体实施中，所述方法可以还包括：基于确定所述读段记录对应于被完全映射到所述参考序列的读段，由一个或多个计算机使用简化的信息熵编码来编码所述读段记录的至少一部分。
[0025]
在一些具体实施中，所述一个或多个计算机可以包括一个或多个硬件处理器。
[0026]
在一些实施方案中，所述一个或多个硬件处理器可以包括一个或多个现场可编程门阵列(fpga)。
[0027]
在一些具体实施中，所述用于压缩基因组序列数据的方法可以由一个或多个硬件处理器执行。在此类具体实施中，所述硬件处理器可以包括被配置为执行一个或多个操作的硬件处理电路系统。在一个方面，所述操作可以包括：由所述硬件处理电路系统获得读段记录；由所述硬件处理电路系统确定所述读段记录是对应于被完全映射到参考序列的读段还是被不完全映射到所述参考序列的读段；基于由所述硬件处理电路系统确定所述读段记录对应于被不完全映射到所述参考序列的读段，由所述一个或多个计算机确定所述不完全映射的读段的失配数目是否满足预先确定的失配阈值数目，以及基于确定所述失配数目满足所述预先确定的失配阈值数目，由所述硬件处理电路系统将所述不完全映射的读段的每个失配编码为具有1个字节大小的记录。
[0028]
根据附图和以下具体实施方式，本发明的这些和其他特征和优点将变得更加明
显。
[0029]
在一些具体实施中，每个读段记录可以包括：指示所述比对读段相对于所述参考序列的绝对起始位置的数据；指示所述读段的长度的数据；指示所述读段是完全映射的还是不完全映射的数据；指示在所述读段中识别的失配的数目的数据；以及指示所述读段中的所述可能失配的相对位置的数据。
[0030]
在一些具体实施中，由所述硬件处理电路系统确定所述不完全映射的读段的失配数目是否满足预先确定的失配阈值数目可以包括由所述硬件处理电路系统确定所述不完全映射的读段的所述失配数目是否大于所述预先确定的失配阈值数目。
[0031]
在一些具体实施中，将所述不完全映射的读段的每个失配编码为具有1个字节大小的记录可以包括：针对每个特定失配，由所述硬件处理电路系统编码所述字节的前两位以包括表示存在于所述读段中的替代核苷酸或碱基而不是所述参考序列中的相应参考核苷酸或碱基的数据；以及由所述硬件处理电路系统编码所述字节的剩余六位以包括表示所述参考序列中的所述失配的位置的数据，所述位置被计算为相对于所述读段的前一个失配的偏移。
[0032]
在一些具体实施中，所述硬件处理器电路系统被进一步配置为执行以下操作，所述操作包括：由所述硬件处理电路系统确定所述偏移是否大于最大可编码值；以及基于确定所述偏移大于所述最大编码值，由所述硬件处理电路系统在所述特定失配与所述前一个失配之间插入至少一个假失配。
[0033]
在一些具体实施中，所述硬件处理器电路系统被进一步配置为执行以下操作，所述操作包括：基于确定所述失配数目不满足所述预先确定的失配阈值数目，由所述硬件处理电路系统使用简化的信息熵编码过程将对应于所述失配中的每一者的位置的所述参考序列的位置列表编码到所述参考序列中。
[0034]
在一些具体实施中，所述硬件处理器电路系统被进一步配置为执行以下操作，所述操作包括：基于确定所述读段记录对应于被完全映射到所述参考序列的读段，由所述硬件处理电路系统使用简化的信息熵编码来编码所述读段记录的至少一部分。
[0035]
在一些具体实施中，该硬件处理电路系统包括一个或多个现场可编程门阵列(fpga)。
[0036]
根据本公开的另一个创新方面，公开了一种用于压缩由测序机器产生的基因组序列数据的计算机实现的方法，所述基因组序列数据包括已与参考序列比对的核苷酸或碱基的序列的读段，从而产生比对读段，所述比对读段作为读段列表存储在初始文件中。在一个方面，所述方法可以包括针对每个比对读段的以下动作：确定所述读段与所述参考序列是完全映射的还是不完全映射的，或者所述读段是否与所述参考序列是未映射的；根据所述确定对所述读段进行编码，其中根据第一编码过程对被确定为完全映射的所述读段进行编码，并且根据第二编码过程对被确定为未映射的所述读段进行编码，其中所述确定步骤包括针对每个不完全映射的读段，将所述读段与所述参考序列之间的失配的数目与阈值进行比较，其中，在所述编码步骤中，根据所述第二编码过程或第三编码过程对被确定为不完全映射的所述读段进行编码，当所述失配数目大于所述阈值时，根据所述第二编码过程对所述不完全映射的读段进行编码，并且当所述失配数目小于所述阈值时，根据所述第三编码过程对所述不完全映射的读段进行编码，其中，在所述第二编码过程中，所述读段的每个核
苷酸或碱基被单独编码，其中所述第一编码过程和所述第三编码过程包括不同的描述符集合，每个描述符集合单义地表示与对应的编码过程相关联的所述读段，所述第一编码过程和所述第三编码过程中的每一者是简化的信息源熵编码过程。
[0037]
其他方面包括执行如本文所公开的方法的动作的对应系统、装置和计算机程序，如由在计算机可读存储设备上编码的指令所定义的。
[0038]
这些和其他版本可任选地包括以下特征中的一个或多个特征。例如，在一些具体实施中，所述确定步骤可以包括当确定读段与所述参考序列不完全映射并且具有低于所述阈值的失配数目时的进一步确定，所述进一步确定有关于所述读段与所述参考序列是全局映射的还是局部映射的，并且其中所述第三编码过程包括第一编码子过程和第二编码子过程，根据所述第一编码子过程对被确定为全局映射的所述读段进行编码，根据所述第二编码子过程对被确定为局部映射的所述读段进行编码，所述第一编码子过程和所述第二编码子过程包括不同的描述符集合，每个描述符集合单义地表示与对应的编码子过程相关联的所述读段。
[0039]
在一些具体实施中，所述第一编码子过程的所述描述符可以包括所述参考序列中的比对起始位置、读段长度和由符号替换表示的失配列表，并且其中所述第二编码子过程的所述描述符包括所述参考序列中的局部比对起始位置、读段长度、由符号替换表示的失配列表，以及不是所述比对的一部分的所述读段的剪切部分的长度。
[0040]
在一些具体实施中，在所述编码步骤中，将要根据所述第二编码子过程编码的读段的所述剪切部分串联，所述剪切部分的每个核苷酸或碱基被单独编码。
[0041]
在一些具体实施中，在所述编码步骤中，在1个字节上编码不完全映射的读段的每个失配。
[0042]
在一些具体实施中，在所述编码步骤中，不完全映射的读段的每个失配如下编码：所述字节的前两位用于编码所述读段中存在的替代核苷酸或碱基而不是所述参考序列中的相应参考核苷酸或碱基，并且所述字节的后六位用于编码所述参考序列中的所述失配的位置，所述位置被计算为相对于所述读段的前一个失配的偏移。
[0043]
在一些具体实施中，在所述编码步骤中，如果在给定失配与所述前一个失配之间计算出的所述偏移大于最大可编码值，则在所述两个失配之间插入至少一个假失配，直到所述失配中的每一者与所述至少一个假失配之间的每个偏移都低于所述最大可编码值，假失配被定义为这样的失配：对于所述失配，所述字节的位用于编码所述失配，或者用于编码与所述参考序列中的相应参考核苷酸或碱基相等的核苷酸或碱基。
[0044]
在一些具体实施中，初始步骤是将所述读段列表划分为读段块，其中每个块开始于包含对所述块进行解码所需的信息的标头，其中所述压缩方法逐块进行。
[0045]
在一些具体实施中，读段块具有相同的块大小。
[0046]
在一些实施方案中，最终步骤是提供包含编码读段列表的压缩文件，所述编码读段以与存储在初始文件中的读段的顺序相同的顺序存储在所述压缩文件中。
[0047]
在一些具体实施中，所述阈值等于31。
[0048]
在一些具体实施中，针对每个比对读段，提供了确定所述读段是否包含至少一个对应于测序机器不能检出任何碱基或核苷酸的情况的失配的步骤。
[0049]
在一些具体实施中，针对包含每个至少一个对应于测序机器不能检出任何碱基或
核苷酸的情况的失配的读段，提供了确定此类失配的数目的步骤，以及将所述数目与参考阈值进行比较的步骤。
[0050]
在一些具体实施中，在所述编码步骤中，如果此类失配的数目大于所述参考阈值，则将要根据所述第二编码过程编码的读段的每个核苷酸或碱基以4位单独编码，并且如果此类失配的数目小于所述参考阈值，则将要根据所述第二编码过程编码的读段的每个核苷酸或碱基以2位单独编码，并且所述编码步骤还包括编码沿着所述参考序列的位置的列表，所述位置对应于所述参考序列中的此类失配的位置。
附图说明
[0051]
图1是示出根据本发明的压缩方法的步骤的流程图。
[0052]
图2是示出用于实现根据本发明的压缩方法的步骤的设备的示意图。
[0053]
图3示出了与参考序列全局映射的读段的第一示例。
[0054]
图4示出了在必须插入假失配的情况下与参考序列全局映射的读段的第二示例。
具体实施方式
[0055]
本发明中提及的基因组序列包括例如但不限于核苷酸序列、脱氧核糖核酸(dna)序列、核糖核酸(rna)序列和氨基酸序列。尽管本文的具体实施方式相当详细地描述了关于核苷酸序列形式的基因组信息，但是应当理解，如本领域技术人员所理解的，根据本发明的压缩方法也可以被实施用于其他基因组序列，尽管有一些变化。
[0056]
基因组测序信息通过测序机器以由来自限定词汇表的字母串表示的核苷酸(或更一般地，碱基)序列的形式产生。最小词汇表由五个符号表示：{a、c、g、t、n}，表示存在于dna中的4种类型的核苷酸，即腺嘌呤、胞嘧啶、鸟嘌呤和胸腺嘧啶。在rna中，胸腺嘧啶被尿嘧啶(u)替代。n指示测序机器不能检出任何碱基，因此该位置的真实性质是不确定的。
[0057]
由测序机器产生的核苷酸序列被称为“读段”。序列读段的长度可以介于几十个核苷酸至几千个核苷酸之间。一些技术成对产生序列读段，其中一个读段来自一条dna链并且第二读段来自另一条链。贯穿本公开，“参考序列”是由测序机器产生的核苷酸或碱基序列在其上比对/映射的任何序列。这种参考序列的一个示例实际上可以是参考基因组，即由科学家组装的作为物种的基因集合的代表性示例的序列。然而，参考序列还可以由合成序列组成，考虑到它们要被进一步处理，该合成序列被设想为仅仅改善读段的可压缩性。测序机器可能在序列读段中引入错误，并且值得注意的是，可能使用错误符号(即，表示不同的核酸)来表示实际存在于经测序的样品中的核酸或碱基；这通常称为替换错误或“失配”。
[0058]
本发明是基于参考的压缩方法，其接收核苷酸或碱基序列的读段作为输入，此类读段先前已与参考序列进行比对，从而产生比对读段。然后将比对读段存储为初始文件中的读段列表。比对读段以及一旦在初始文件中比对就将这些读段存储的方式对于本发明不是至关重要的，并且不是本公开的目的。然后将每个读段编码为参考序列上的位置以及与所述参考序列的差异的列表。然后可以通过根据本发明配置的适当解压缩软件，从比对编码信息和参考序列重建每个读段。
[0059]
优选地，对读段进行处理并在将它们作为输入提供给压缩软件和设备之前将它们与参考序列比对的比对软件不考虑序列读段中引入的某些类型的错误，例如插入错误或缺
失错误。插入错误在于，在一个序列读段中插入一个或多个不涉及任何实际存在的核酸的额外符号。缺失错误在于，从一个序列读段中缺失一个或多个代表实际存在于测序样品中的核酸的符号。更准确地讲，在给定序列读段中存在插入错误或缺失错误的情况下，比对软件随后将所得到的错误核酸视为替换错误，也称为“失配”。比对软件配置的这种优先选择允许更快的后续编码，特别是在速度与压缩比之间提供更好的折衷。
[0060]
对于每个读段，比对软件向压缩软件和设备提供对应的读段记录。每个读段记录至少包含以下信息：比对读段相对于参考序列的绝对起始位置、读段的长度、读段的比对类型、读段中识别的失配的数目以及读段中所述可能失配的相对位置(在适当的情况下)。
[0061]
现在将参考图1描述根据本发明的压缩方法。该方法例如由图2中所示的设备20执行。该设备包括至少一个处理器22和可操作地耦接到处理器22的一个存储器24以形成计算装置。存储器24可以存储计算机程序代码或软件26，所述计算机程序代码或软件包含计算机可执行指令，这些计算机可执行指令在由处理器22执行时，使得处理器22执行包括根据本发明的压缩方法的步骤的操作。
[0062]
其中存储了比对读段作为读段列表的初始文件例如存储在设备20的存储器中。返回图1，该方法优选地包括初始步骤2，其中将比对读段的初始列表划分为读段块。通常，将比对读段列表划分为50,000个读段块，该特定值不被解释为限制可用其它值以相同方式应用的本发明的范围。优选地，读段块具有相同的块大小。每个读段块开始于包含对该块进行解码所需的信息的标头，诸如该块的内容的以字节为单位的大小，以及/或者该块或其内容的标识符，以及/或者包含在该块中的读段的数目。这允许支持压缩文件的串联，以及流式传输能力(每个读段块均包含对该块的读段进行解码所需的所有信息)。此外，由于压缩方法然后可以逐块进行，所以这也允许对读段块进行多线程处理，从而允许并行化和处理时间的一些所得增益。如果给定块的所有读段具有相同的长度，则读段长度也被存储在标头中，否则在压缩方法期间明确地存储每个读段长度的列表。
[0063]
每个读段记录包含关于读段比对的类型的信息。通常，可识别两种主要的比对类型：完全比对和不完全比对，加上对应于“未映射”读段的额外类型。“不完全比对”意味着读段含有至少一个非n失配，而读段的至少一部分与参考序列的一部分匹配(根据该定义，不完全映射的读段可含有一个或多个n，条件是其还含有一个或多个其它失配)。在示例性实施方案中，每个读段记录以以下位标记开始，每个位标记具有在两个可能值之间的一个值：
[0064]-第一位标记，其指示相对于参考序列的正向或反向取向，
[0065]-第二位标记，其指示完全比对或未完全比对，
[0066]-第三位标记，其指示读段是否包含至少一个n，
[0067]-第四位标记，其指示位置信息是以16位还是32位编码。
[0068]
以下步骤4-12以逐读段块进行，并且在块内以逐读段进行。
[0069]
该方法包括下一步骤4，其针对每个比对读段，确定所述读段与参考序列是完全映射的还是不完全映射的，或者所述读段是否与参考序列是未映射的。该确定步骤4包括针对每个不完全映射的读段，将所述读段与参考序列之间的失配的数目与阈值进行比较4a。在优选的实施方案中，尽管不应被解释为限制本发明的范围，但所述阈值等于31。特意选择该特定值，以便提供用于以足够紧凑的方式存储失配数目的最佳可能折衷，如稍后关于步骤12将更好理解的。实际上，已经在统计学上观察到，在绝大多数情况下，不完全映射的读段
具有少于31个失配。该选择背后的原理在于以最紧凑的方式对最频繁出现的情况进行编码，留下一些非常少的降级情况。如果读段被确定为以低于阈值的失配数目发生不完全映射，则确定步骤4包括关于读段与参考序列是全局映射的还是局部映射的进一步确定。“全局映射的读段”是不完全映射的读段，其整个序列(包括该读段的起点和终点)与参考序列是不完全映射的。“局部映射的读段”是含有与参考序列是不完全映射的核苷酸或碱基的区段的不完全映射的读段。因此，核苷酸或碱基的所述区段对应于初始读段的一部分。
[0070]
优选地，该方法还包括步骤6，其针对每个比对读段，确定所述读段是否包含至少一个n，即，所述读段是否包含至少一个对应于测序机器不能检出任何碱基或核苷酸的情况的失配。对于包含至少一个n的每个读段，该方法然后包括步骤8，其确定此类n失配的数目；和步骤10，其将所述n失配的数目与参考阈值进行比较。在优选的实施方案中，尽管不应被解释为限制本发明的范围，但所述参考阈值等于31。
[0071]
无论确定步骤4的结果如何，该方法均包括下一个步骤12，其至少根据所述确定对读段进行编码。更准确地讲，根据第一编码过程对被确定为与参考序列完全映射的读段进行编码，无论所述读段是不包含n，还是具有数目低于参考阈值的n。根据第二编码过程对被确定为未映射的读段或者被确定为完全映射但具有数目大于参考阈值的n的读段进行编码，其中每个核苷酸或碱基被单独编码，而不管所述核苷酸或碱基是比对的还是未比对的。根据第二编码过程或第三编码过程对被确定为不完全映射的读段进行编码。更准确地讲，根据第二编码过程对被确定为以大于阈值的失配数目发生不完全映射的读段进行编码。如果读段被确定为以低于阈值的失配数目发生不完全映射，则如果所述读段不包含n或具有低于参考阈值的n数目，则根据第三编码过程对所述读段进行编码。如果不是这样，即，如果读段具有大于参考阈值的n数目，则根据第二编码过程对所述读段进行编码。
[0072]
不管给定读段是已被确定为完全映射的、不完全映射的还是未映射的，如果所述读段包含至少一个n但具有低于参考阈值的n数目，则编码步骤12包括编码沿着参考序列的位置的列表，所述位置对应于参考序列中的n的位置。然后将该位置列表存储在计算装置的存储器中，所述装置实施所述压缩方法。如果读段包含至少一个n但具有低于参考阈值的n数目，则将根据第二编码过程进行编码，并且该读段的每个核苷酸或碱基以2位单独编码。
[0073]
如果读段包含至少一个n但具有大于参考阈值的n数目，则在任何情况下都将根据第二编码过程对所述读段进行编码，并且该读段的每个核苷酸或碱基以4位单独编码。在这种情况下，编码步骤12不包括对参考序列中的n位置的列表进行编码和存储。实际上，然后根据第二编码过程以与读段的其他核苷酸或碱基极其相同的方式对每个n失配直接进行编码。
[0074]
第一编码过程和第三编码过程包括不同的描述符集合。每个描述符集合单义地表示与对应的编码过程相关联的读段，第一编码过程和第三编码过程中的每一者是简化的信息熵编码过程。更准确地讲，第三编码过程包括第一编码子过程和第二编码子过程。根据第一编码子过程对在步骤4期间被确定为全局映射的不完全映射的读段进行编码。根据第二编码子过程对在步骤4期间被确定为局部映射的不完全映射的读段进行编码。第一编码子过程和第二编码子过程包括不同的描述符集合，每个描述符集合单义地表示与对应的编码子过程相关联的读段。
[0075]
于是，针对每个读段编码并且使得能够在数据解压缩期间重建整个读段序列的比
对信息取决于用于所述读段的对应的编码过程或子过程。例如，用于第一编码过程的描述符可以是：
[0076]
ο完全映射的读段相对于参考序列的绝对起始位置(以16位或32位编码)，以及
[0077]
ο读段的长度(采用差分编码法相对于前一个读段的长度进行编码，其中可变长度码在2位至34位的范围内)。
[0078]
用于第一编码子过程的描述符可以是：
[0079]
ο不完全映射的读段相对于参考序列的绝对起始位置(以16位或32位编码)，
[0080]
ο读段的长度(采用差分编码法相对于前一个读段的长度进行编码，其中可变长度码在2位至34位的范围内)，以及
[0081]
ο读段的失配列表。
[0082]
用于第二编码子过程的描述符可以是：
[0083]
ο读段的不完全映射的部分相对于参考序列的绝对起始位置
–
也称为局部比对起始位置(以16位或32位编码)，
[0084]
ο读段的长度(采用差分编码法相对于前一个读段的长度进行编码，其中可变长度码在2位至34位的范围内)，
[0085]
ο读段的失配列表，以及
[0086]
ο读段的不是比对的一部分的剪切部分的长度(对于每个剪切部分以8位编码)。
[0087]
优选地，在第一子过程和第二子过程中编码的失配列表包括标头(在1个字节上编码的位标记)。该字节的前五位用于编码包含在读段中的失配数目(在优选的实施方案中，其中阈值等于31，所述数目在范围[0-31]内)。然后一位用于编码不完全映射的读段是全局映射还是局部映射。另一位用于编码2位模式是否被激活用于第二编码过程。最后一位用于编码4位模式是否被激活用于第二编码过程。优选地，针对在编码步骤12期间根据第二编码子过程编码的每个读段，将所述读段的剪切部分(即，不是局部比对的一部分的那些部分)串联，并且所述剪切部分的每个核苷酸或碱基被单独编码。在优选的实施方案中，所述读段的这种剪切部分的每个核苷酸或碱基以2位单独编码。
[0088]
在优选的实施方案中，在不完全映射的读段的失配列表中编码(即，根据第一编码子过程或第二编码子过程编码)的每个失配在1个字节上编码。更准确地讲，要根据第一编码子过程或第二编码子过程编码的不完全映射的读段的每个失配可以如下编码：
[0089]
ο字节的前两位用于编码读段中存在的替代核苷酸或碱基，而不是参考序列中的相应参考核苷酸或碱基，
[0090]
ο最后六位用于编码参考序列中的失配的位置，所述位置被计算为相对于读段的前一个失配的偏移(失配的相对位置，除了读段的第一失配被编码的绝对位置之外)。因此，该偏移的范围(以6位编码)为[0-63]。
[0091]
图3提供了根据第一编码子过程对读段的失配进行编码的一个示例。该读段是不完全映射的读段，其与参考序列是全局映射的。该读段具有两个失配：
[0092]
ο第一失配，其位于该读段中的第12位，其在于在该读段中用t核苷酸替换参考序列中的a核苷酸，以及
[0093]
ο第二失配，其位于该读段中的第21位，其在于在该读段中用g核苷酸替换参考序列中的c核苷酸。
[0094]
该读段的失配列表然后被编码为：
[0095]
ο《12,t》，值“12”对应于该读段中的第一失配的绝对位置，以及
[0096]
ο《9,g》，值“9”对应于该读段中的第二失配的相对位置，即第二失配与第一失配之间的偏移。
[0097]
例如，《12,t》可以被转换为值“51”(在1个字节上编码)，并且《9,g》可以被转换为值“38”(在1个字节上编码)。这种字节编码用以下方式获得：
[0098]
偏移位置x4 核苷酸值(其中a＝0，c＝1，g＝2，t＝3)
[0099]
优选地，对于要根据第一编码子过程或第二编码子过程进行编码的每个不完全映射的读段，如果在该读段的给定失配与前一个失配之间计算出的偏移大于最大可编码值，则在所述两个失配之间插入至少一个“假”失配，直到所述失配中的每一者与所述至少一个“假”失配之间的每个偏移都低于所述最大可编码值。“假”失配被定义为这样的失配：对于该失配，字节的位用于编码该失配，或者用于编码与参考序列中的相应参考核苷酸或碱基相等的核苷酸或碱基。在优选的实施方案中，尽管不应被解释为限制本发明的范围，但最大可编码值等于63，对应于以6位编码的最大值。
[0100]
图4提供了在必须插入“假”失配的情况下，根据第一编码子过程对读段的失配进行编码的一个示例。该读段是不完全映射的读段，其与参考序列是全局映射的。该读段具有两个失配：
[0101]
ο第一失配，其位于该读段中的第22位，其在于在该读段中用t核苷酸替换参考序列中的a核苷酸，以及
[0102]
ο第二失配，其位于该读段中的第134位，其在于在该读段中用g核苷酸替换参考序列中的c核苷酸。
[0103]
第二失配与第一失配之间的位置偏移为112，其大于最大可编码值63。因此，“假”失配必须插入所述两个失配之间，使得所述失配中的每一者与该“假”失配之间的每个偏移都低于所述最大可编码值。具有t核苷酸的“假”失配(对应于参考序列中的“真”t核苷酸)例如插入该读段中的第85位。该“假”失配与第一失配之间计算出的位置偏移为63，其对应于最大可编码值。第二失配与该“假”失配之间计算出的位置偏移为49，其低于63。
[0104]
该读段的失配列表然后被编码为：
[0105]
ο《22,t》，值“22”对应于该读段中的第一失配的绝对位置，
[0106]
ο《63,t》，值“63”对应于该读段中的“假”失配的相对位置，即“假”失配与第一失配之间的偏移，以及
[0107]
ο《49,g》，值“49”对应于该读段中的第二失配的相对位置，即第二失配与“假”失配之间的偏移。
[0108]
例如，《22,t》可以被转换为值“91”(在1个字节上编码)，《63,t》可以被转换为值“255”(在1个字节上编码)，并且《49,g》可以被转换为值“198”(在1个字节上编码)。这种字节编码用以下方式获得：
[0109]
偏移位置x4 核苷酸值(其中a＝0，c＝1，g＝2，t＝3)
[0110]
所述方法包括提供包含编码读段列表的压缩文件的最终步骤14。所述编码读段以与存储在初始未压缩文件中的读段的顺序相同的顺序存储在所述压缩文件中。然后可以通过根据本发明配置的适当解压缩软件和/或方法，从比对编码信息和参考序列重建每个读
段。
[0111]
尽管参考计算装置20的示例性架构(出于说明的目的在图2中示出)进行了描述，但是此处公开的本发明的技术可以在硬件、软件、固件或它们的任何组合中实现。当在软件中实现时，计算机程序代码可以存储在计算机介质上，并且由包括一个或多个处理器的硬件处理单元执行，这与使用图2的装置20时的情况相同。应当理解，如本文所用的术语“处理器”旨在包括一个或多个处理装置，包括信号处理器、微处理器、微控制器、专用集成电路(asic)、现场可编程门阵列(fpga)或其他类型的处理电路系统，以及此类电路系统元件的部分或组合。另外，如本文所用的术语“存储器”旨在包括与处理器相关联的电子存储器，诸如以任何组合使用的随机存取存储器(ram)、只读存储器(rom)或其他类型的存储器。
[0112]
因此，用于执行本文所述的方法和协议的软件指令或代码可以存储在相关联的存储器装置(例如，rom、固定存储器或可移动存储器)中的一者或多者中，并且当准备使用时，被加载到ram中并由处理器执行。
[0113]
本公开的技术可以在各种各样的装置或设备中实施，包括例如移动电话、计算机、服务器、平板电脑和类似装置。
[0114]
尽管本文已参考附图描述了本发明的说明性实施方案，但是应当理解，本发明并不限于那些确切的实施方案，并且在不脱离本发明的范围或实质的情况下，本领域的技术人员可以作出各种其他的改变和修改。
[0115]
根据本发明的压缩方法的统计和数值示例
[0116]
以下比较实施例是在未压缩数据文件上进行的，该文件含有4800万个核苷酸读段或序列：
[0117]
ο未压缩数据文件的大小：35,770mb(兆字节)
[0118]
ο已用gzip软件压缩的文件的大小：6,649mb
[0119]
ο已用基于非参考的spring软件压缩的文件的大小：1,402mb
[0120]
ο已用根据本发明的基于参考的压缩方法压缩的文件的大小：1,179mb
[0121]
ο使用基于非参考的spring软件的压缩时间：1,722s
[0122]
ο使用根据本发明的基于参考的压缩方法的压缩时间：181s
[0123]
ο未压缩数据文件的以位/核苷酸计的平均大小(ascii编码)：8位/核苷酸
[0124]
ο已用适于4个可能字符a、t、c、g的编码压缩的文件的以位/核苷酸计的平均大小：2位/核苷酸
[0125]
ο已用根据本发明的基于参考的压缩方法压缩的文件的以位/核苷酸计的平均大小：0.33位/核苷酸
[0126]
上文指出的数值示例说明本发明允许快速压缩和解压缩，同时提供高压缩比。