CN110663022B - 使用基因组描述符紧凑表示生物信息学数据的方法和设备 - Google Patents

使用基因组描述符紧凑表示生物信息学数据的方法和设备 Download PDF

Info

Publication number
CN110663022B
CN110663022B CN201880012026.4A CN201880012026A CN110663022B CN 110663022 B CN110663022 B CN 110663022B CN 201880012026 A CN201880012026 A CN 201880012026A CN 110663022 B CN110663022 B CN 110663022B
Authority
CN
China
Prior art keywords
class
reads
read
information
descriptor
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201880012026.4A
Other languages
English (en)
Other versions
CN110663022A (zh
Inventor
克劳迪奥·阿尔贝蒂
乔吉奥·索亚
丹尼尔·伦齐
穆罕默德·霍索·巴鲁克
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Genomsys SA
Original Assignee
Genomsys SA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from PCT/EP2016/074301 external-priority patent/WO2018068828A1/en
Priority claimed from PCT/EP2016/074297 external-priority patent/WO2018068827A1/en
Priority claimed from PCT/EP2016/074311 external-priority patent/WO2018068830A1/en
Priority claimed from PCT/US2017/017842 external-priority patent/WO2018071055A1/en
Application filed by Genomsys SA filed Critical Genomsys SA
Priority claimed from PCT/US2018/018092 external-priority patent/WO2018152143A1/en
Publication of CN110663022A publication Critical patent/CN110663022A/zh
Application granted granted Critical
Publication of CN110663022B publication Critical patent/CN110663022B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

用于压缩由基因组测序机器所产生的基因组序列数据的方法和设备。通过将序列读段相对于预先存在或构建的参考序列进行比对来编码序列读段,编码处理包括将读段分类为数据类别,然后根据多个描述符块对每个类别进行编码。特定的源模型和熵编码器用于划分数据的每个数据类别,以及每个相关联的描述符块。

Description

使用基因组描述符紧凑表示生物信息学数据的方法和设备
交叉引用相关申请
本申请要求2017年2月14日提交的专利申请PCT/US2017/017842和2017年7月11日提交的专利申请PCT/US2017/041591的优先权和权益。
技术领域
本公开提供了一种表示基因组测序数据的新方法,该方法通过提供已知现有技术表示方法不可用的新功能,减少了所利用的存储空间并提高了访问性能。
背景技术
基因组测序数据的适当表示是实现高效基因组分析应用的基础,该应用诸如是基因组变体调用以及通过处理测序数据和元数据以各种目的进行的所有其他分析。
高通量低成本测序技术的出现使人类基因组测序变得负担得起。这样的机会在从癌症的诊断和治疗到遗传疾病的识别、从用于抗体识别的病原体监测到新疫苗、药物的创建和个性化治疗的定制等多个领域中开辟了新的前景。
医院、基因组数据分析提供商、生物信息学家和大型生物数据存储中心正在寻找负担得起、快速、可靠和互联的基因组信息处理解决方案,这些解决方案能够将基因组医学扩展到世界范围。因为测序处理中的一个瓶颈已经变为数据存储,所以用于以压缩形式来表示基因组测序数据的方法越来越受到研究。
测序数据最常用的基因组信息表示是基于压缩FASTQ和SAM格式。目标是压缩传统使用的文件格式(分别用于不比对和比对数据的FASTQ和SAM)。这样的文件由纯文本字符构成,并且如上所述,通过使用通用方法被压缩,通用方法诸如是LZ(来自发布了第一个版本的作者Lempel和Ziv)方案(众所周知的zip、gzip等)。当使用诸如gzip的通用压缩器时,压缩的结果通常是二进制数据的单个块。特别是当像在高通量测序的情况下,数据量非常大时,这样的单体形式的信息很难存档、传输和阐述。BAM格式的特征是压缩性能差,这是由于集中于压缩低效和冗余的SAM格式,而不是提取由SAM文件所传递的实际基因组信息,并且由于采用诸如gzip的通用文本压缩算法,而不是利用每个数据源(基因组数据本身)的特定性质。
CRAM是一种更复杂的基因组数据压缩方法,使用较少,但比BAM更有效。CRAM为采用相对于参考的差分编码提供了更有效的压缩(CRAM部分利用了数据源冗余),但是CRAM仍然缺乏诸如增量更新、流的支持以及对特定类别压缩数据的选择性访问的特征。
这些方法生成较差的压缩率和数据结构,一旦压缩,这些数据结构就难以导航和操作。由于即使为了进行简单的操作或访问基因组数据集的选定区域,也需要处理大而刚性的数据结构,因此下游分析可能非常慢。CRAM依赖于CRAM记录的概念。每个CRAM记录通过对重建它所必需的所有元素进行编码来表示单个映射或未映射的读段。
CRAM提出了以下缺点和限制,这些缺点和限制通过本文档中所描述的发明得以解决和克服:
1.CRAM不支持数据索引和对共享特定特征的数据子集的随机访问。数据索引超出规范(参见CRAM规范v3.0的第12节)的范围并且作为单独的文件来实施。相反,本文档中所描述的发明的方法采用与编码处理集成的数据索引方法,并且索引被嵌入到编码的(即,压缩的)比特流中。
2.CRAM是由核心数据块构建的,这些数据块可以包含任何类型的映射的读段(完全匹配的读段、仅具有替换的读段、带有插入或删除(也被称为“插入缺失”)的读段)。根据相对于参考序列的映射结果,不存在数据分类和读段分组的概念。这意味着即使只搜索具有特定特征的读段,也需要检查所有数据。本发明通过在编码之前对数据进行分类和划分来解决这样的限制。
3.CRAM是基于将每个读段都封装到“CRAM记录”中的概念。这意味着在搜索由特定生物特征所表征的读段(例如,具有替换但没有“插入缺失”的读段或完全映射的读段)时,需要检查每个完整的“记录”。
相反,在本发明中,存在在分离的信息块中单独编码的数据类的概念,并且不存在封装每个读段的记录的概念。这使得能够更有效地访问具有特定生物学特征的读段组(例如,具有替换但没有“插入缺失”的读段或完全映射的读段),而不需要对每个读段(块)进行解码以检查其特征。
4.在CRAM记录中,每个记录字段与特定的标志相关联,并且每个标志必须始终具有相同的含义,因为不存在情境的概念,这是因为每个CRAM记录可以包含任何不同类型的数据。该编码机制引入了冗余信息,并阻止了基于情境的高效熵编码的使用。
相反,在本发明中,不存在表示数据的标志的概念,因为这是由数据所属的信息“块”内在定义的。这意味着要被使用的符号数量大大减少,随之而来的信息源熵减少,其导致更有效的压缩。这样的改进是可能的,因为不同“块”的使用使得编码器能够根据情境在具有不同含义的每个块上再次使用相同的符号。在CRAM中,每个标志必须始终具有相同的含义,因为不存在情境的概念,并且每个CRAM记录可以包含任何类型的数据。
5.在CRAM替换中,通过使用不同的描述符来表示插入和删除,这种选项增加了信息源字母表的大小,并且产生了更高的源熵。相反,所公开的发明的方法使用单个字母表和编码以用于替换、插入和删除。这使得编码和解码处理更简单,并且产生更低的熵源模型,该模型编码产生由高压缩性能所表征的比特流。
本发明旨在通过分类和划分测序数据来压缩基因组序列,使得要被编码的冗余信息最小化,并且在压缩域中直接启用诸如选择性访问和支持增量更新的特征。
所提出的方法的一个方面是定义在不同块中结构化并且分别编码的数据和元数据的类。这样的方法相对于现有方法的更相关的改进在于:
1.压缩性能的提高,这是由于通过为每类数据或元数据提供有效的源模型而构成的信息源熵的降低;
2.直接在压缩域中为任何进一步的处理目的对压缩数据和元数据的部分进行选择性访问的可能性;
3.用与特定序列读段集相关联的新序列数据和/或元数据和/或新分析结果来递增(即,不需要解码和重新编码)更新压缩数据和元数据的可能性。
附图说明
图1示出了映射读段对的位置如何在“pos”块中被编码为与第一映射读段的绝对位置的差异。
图2示出了配对中的两个读段如何来自两个DNA链。
图3示出了如果链1被用作参考,则读段2的反向补码如何被编码。
图4示出了组成读段对的四种可能的读段组合以及“rcomp”块中的相应编码。
图5示出了在三个读段对的读段长度不变的情况下,如何计算配对距离。
图6示出了在“pair”块中所编码的配对错误如何使解码器能够使用编码的“MPPPD”来重建正确的读段配对。
图7示出了当读段被映射到与其配对不同的参考上时配对距离的编码。在该情况下,额外的描述符被添加到配对距离。一个是表示标志,第二个是参考标识符,然后是配对距离。
图8示出了“nmis”块中“n型”错配的编码。
图9示出了给出相对于参考序列的替换的映射的读段对。
图10示出了如何将替换的位置计算为绝对值或差值。
图11示出了当不使用IUPAC代码时,如何计算编码替换类型的符号。这些符号以圆形替代向量来表示读段中存在的分子和该位置上的参考上存在的分子之间的距离。
图12示出了如何将替换编码到“snpt”块中。
图13示出了当使用IUPAC模糊代码时如何计算替换代码。
图14示出了当使用IUPAC代码时如何对“snpt”块进行编码。
图15示出了对于I类读段,所使用的替换向量与M类如何相同,只是增加了插入符号A、C、G、T、N的特殊代码。
图16示出了在IUPAC模糊代码的情况下,错配和插入缺失的编码的一些示例。在该情况下,替换向量要长得多,因此,可能被计算的符号比五个符号的情况多。
图17示出了错配和插入缺失的不同源模型,其中,每个块包含单一类型的错配或插入的位置。在该情况下,不会为错配或插入缺失类型编码任何符号。
图18示出了错配和插入缺失编码的示例。当读段不存在给定类型的错配或插入缺失时,在相对应的块中编码0。0充当每个块中的读段分隔符和终止符。
图19示出了参考序列中的修改可以如何将M读段转换为P读段。该操作可以降低数据结构的信息熵,尤其是在高覆盖数据的情况下。
图20示出了根据本发明的一个实施例的基因组编码器2010。
图21示出了根据本发明的一个实施例的基因组解码器218。
图22示出了可以如何通过聚类读段并且组装从每个聚类中所获取的片段来构建“内部”参考。
图23示出了构建参考的策略如何在对读段已经应用特定排序(例如,字典顺序)后存储最近的读段。
图24示出了可以如何使用在相对应的块中所存储或携带的六个描述符对属于“未映射”的读段类(U类)的读段进行编码。
图25示出了属于U类的读段的替代编码,其中,有符号的pos描述符被用来编码读段在构造的参考上的映射位置。
图26示出了可以如何应用参考转换来消除读段中的错配。在一些情况下,参考转换可能生成新的错配,或者更改在已经应用转换之前参考该参考时所发现的错配的类型。
图27示出了当移除错配的所有或子集时,参考转换可以如何改变类读段(即,在已经应用参考转换之后,在转换被分配给P类之前属于M类的读段)。
图28示出了半映射的读段对(HM类)可以如何通过用未映射的读段组装较长重叠群来填充参考序列的未知区域。
图29示出了如何用阈值向量来配置N、M和I类数据的编码器并且生成N、M和I数据类的单独子类。
图30示出了所有的数据类可以如何使用相同的被转换的参考进行重新编码,或者不同的转换可以用于每个N、M和I类或者其任意组合。
图31示出了基因组数据集报头的结构。
图32示出了主索引表的一般结构,其中,每一行包含几类数据P、N、M、I、U、HM的基因组间隔以及指向元数据和注释的其他指针。列指代与所编码的基因组数据有关的参考序列上的特定位置。
图33示出了一行MIT的示例,该行包含与P类的读段有关的基因组间隔。与不同参考序列有关的基因组区域由特殊标志(在该示例中为“S”)分隔。
图34示出了本地索引表(LIT)的一般结构以及如何使用它来存储指向所存储或传输的数据中编码的基因组信息的物理位置的指针。
图35示出了被用来访问块有效载荷中的7号和8号访问单元的LIT示例。
图36示出了基因组块报头中所包含的MIT和LIT的多个行之间的功能关系。
图37示出了访问单元是如何由包含属于不同类的数据的不同基因组流所携带的多个基因组数据块组成的。每个块还由被用作数据传输单元的数据包组成。
图38示出了访问单元是如何由报头和属于同质数据的一个或多个块的多路复用块组成的。每个块可以由包含基因组信息的实际描述符的一个或多个包组成。
图39示出了没有拼接的多个比对。最左侧的读段有N个比对。N是要被解码的mmap的第一值,并且表示第一读段的比对的数量。mmap描述符的以下N个值被解码并且被用来计算P,P是第二读段的比对的数量。
图40示出了如何使用pos、pair和mmap描述符来编码没有拼接的多个比对。最左侧的读段有N个比对。
图41示出了具有拼接的多个比对。
图42示出了使用pos、pair、mmap和msar描述符来表示具有拼接的多个比对。
发明内容
以下权利要求的特征通过提供以下内容解决了现有技术解决方案的问题,
一种用于编码基因组序列数据的方法,所述基因组序列数据包括核苷酸序列的读段,所述方法包括以下步骤:
将所述读段与一个或多个参考序列进行比对,从而创建比对的读段,
根据具有所述一个或多个参考序列的指定匹配规则对所述比对的读段进行分类,从而创建比对的读段的类别,
将所述分类的比对的读段编码为多个描述符块,
其中,将所述分类的比对的读段编码为多个描述符块,包括根据所述比对的读段的类别选择所述描述符,
用报头信息来构建所述描述符块,从而创建连续的访问单元。
在另一方面,编码方法还包括:进一步将不满足所述指定匹配规则的所述读段分类为一类未映射的读段,
使用至少一些未映射的读段来构建一组参考序列,
将所述一类未映射的读段与所述一组构建的参考序列进行比对,
将所述分类的比对的读段编码为多个描述符块,
编码所述一组构建的参考序列,
用报头信息来构建所述描述符块和所述编码的参考序列,从而创建连续的访问单元。
在另一方面,编码方法还包括将所述参考序列中没有任何错配的基因组读段识别为第一“P类”。
在另一方面,编码方法还包括当仅在测序机器不能调用任何“碱基”的位置发现错配,并且每个读段中的所述错配的数量不超过给定阈值时,将基因组读段识别为第二“N类”。
在另一方面,编码方法还包括当在所述测序机器不能调用任何“碱基”,被称为“n型”错配,和/或调用与所述参考序列不同的“碱基”,被称为“s型”错配,的位置发现错配,并且所述错配的数量不超过所述“n型”、“s型”错配数量的给定阈值和从基于“n型”和“s型”错配数量计算的给定函数(f(n,s))获得的阈值时,将基因组读段识别为第三“M类”。
在另一方面,编码方法还包括当基因组读段可能具有相同类型的“M类”错配,以及另外以下类型的至少一个错配:“插入”(“i型”)、“删除”(“d型”)、软剪切(“c型”),并且其中,所述每个类型的错配的数量不超过相对应的给定阈值和由给定函数(w(n,s,i,d,c))提供的阈值时,将基因组读段识别为第四“I类”。
在另一方面,编码方法还包括将基因组读段识别为第五“U类”,所述第五“U类”包括在所述P、N、M、I类中没有找到任何分类的所有读段,如上所述。
在另一方面,编码方法还包括要被编码的所述基因组序列的所述读段是配对的。
在另一方面,编码方法还包括,所述分类还包括将基因组读段识别为第六“HM类”,所述第六“HM类”包括所有读段对,其中,一个读段属于P、N、M或I类,另一读段属于“U类”。
在另一方面,编码方法还包括以下步骤:识别所述两个配对读段是否被分类在相同的所述类别(以下中的每一个:P、N、M、I、U)中,然后将所述配对分配给所述相同的识别的类别,
识别所述两个配对读段是否被分类在不同的类别中,并且在它们都不属于“U类”的情况下,将所述读段对分配给具有根据以下表达式定义的最高优先级的所述类别:
P<N<M<I
其中,“P类”具有最低优先级,“I类”具有最高优先级;
识别所述两个配对读段中是否只有一个已经被分类为属于“U类”,并且将所述读段对分类为属于所述“HM类”序列。
在另一方面,编码方法还包括,根据通过所述“n型”错配的数量(292)、所述函数f(n,s)(293)和所述函数w(n,s,i,d,c)(294),分别为每个N、M和I类定义的阈值的向量(292,293,294),读段N、M、I的每个类别被进一步分成两个或多个子类(296,297,298)。
在另一方面,编码方法还包括以下步骤:
识别所述两个配对读段是否被分类在相同的子类中,然后将所述配对分配给所述相同的子类,
识别所述两个配对读段是否被分类为不同类别的子类,然后根据以下表达式将所述配对分配给属于较高优先级类别的所述子类:
N<M<I
其中,N具有最低优先级,I具有最高优先级;
识别所述两个配对读段是否被分类在相同的类别中,并且这样的类别是N或M或I,但是被分类在不同的子类中,然后根据以下表达式将所述配对分配给具有最高优先级的所述子类:
N1<N2<…<Nk
M1<M2<…<Mj
I1<I2<…<Ih
其中,最高索引具有最高优先级。
在另一方面,通过“pos”描述符块来编码关于每个读段的所述映射位置的信息。
在另一方面,通过rcomp描述符块来编码关于每个读段的所述链特异性(即,所述读段序列来自的DNA链)的信息。
在另一方面,通过“pair”描述符块来编码配对端读段的配对信息。
在另一方面,通过“flags”描述符块来编码附加比对信息,所述附加比对信息诸如是所述读段是否被映射在适当的配对中、所述读段未通过平台/供应商质量检查、所述读段是PCR或光学复制、或者所述读段是补充的比对。
在另一方面,通过“nmis”描述符块来编码关于未知碱基的信息。
在另一方面,通过“snpp”描述符块来编码关于替换的所述位置的信息。
在另一方面,通过特定的“snpt”描述符块来编码关于替换的所述类型的信息。
在另一方面,通过“indp”描述符块来编码关于替换、插入或删除类型的错配的所述位置的信息。
在另一方面,通过“indt”描述符块来编码关于诸如替换、插入或删除的错配的信息。
在另一方面,通过“indc”描述符块来编码关于映射的读段的剪切的碱基的信息。
在另一方面,通过“ureads”描述符块来编码关于未映射的读段的信息。
在另一方面,通过“rtype”描述符块来编码关于用于编码的参考序列的所述类型的信息。
在另一方面,通过“mmap”描述符块来编码关于所述映射的读段的多个比对的信息。
在另一方面,通过“msar”描述符块和“mmap”描述符块来编码关于相同读段的拼接的比对和多个比对的信息。
在另一方面,通过“mscore”描述符块来编码关于读段比对分数的信息。
在另一方面,通过“特定rgroup”描述符块来编码关于读段所属的所述组的信息。
在另一方面,编码方法还包括,所述描述符块包括“主索引表”,所述主索引表包含用于比对的读段的每个类别和子类的一个部分,所述部分包括每个类别或子类数据的每个访问单元的第一读段的所述一个或多个参考序列上的所述映射位置;联合编码所述“主索引表”和所述访问单元数据。
在另一方面,编码方法还包括,所述描述符块还包括关于所使用的参考的所述类型(预先存在或构建的)的信息,以及不映射到所述参考序列上的所述读段的所述片段。
在另一方面,编码方法还包括,所述参考序列首先通过应用替换、插入、删除和剪切被转换为不同的参考序列,然后将所述分类的比对的读段编码为参考所述转换的参考序列的多个描述符块。
在另一方面,编码方法还包括,相同的转换被应用于用于所有数据类别的所述参考序列。
在另一方面,编码方法还包括,不同的转换被应用于用于每个数据类别的所述参考序列。
在另一方面,编码方法还包括,所述参考序列转换被编码为描述符块,并且用报头信息来构建,从而创建连续的访问单元。
在另一方面,编码方法还包括,将所述分类的比对的读段和所述有关的参考序列转换编码成多个描述符块,包括将特定源模型和特定熵编码器与每个描述符块相关联的步骤。
在另一方面,编码方法还包括,所述熵编码器是情境自适应算术编码器、可变长度编码器或golomb编码器中的一个。
本发明还提供了一种用于解码编码的基因组数据的方法,包括以下步骤:
解析包含所述编码的基因组数据的访问单元,以通过使用报头信息来提取多个描述符块,
解码所述多个描述符块,以根据相对于一个或多个参考序列定义读段的分类的特定匹配规则来提取读段。
在另一方面,解码方法还包括解码未映射的基因组读段。
在另一方面,解码方法还包括解码分类的基因组读段。
在另一方面,解码方法还包括解码主索引表,所述主索引表包含每个读段类别的一个部分和所述相关联的相关映射位置。
在另一方面,解码方法还包括解码与所使用的参考的所述类型有关的信息:预先存在的、转换的或构建的。
在另一方面,解码方法还包括解码与要被应用于所述预先存在的参考序列的一个或多个转换有关的信息。
在另一方面,解码方法还包括配对的基因组配对。
在另一方面,解码方法还包括所述基因组数据被熵解码。
本发明还提供了一种用于压缩基因组序列数据209的基因组编码器(2010),所述基因组序列数据209包括核苷酸序列的读段,
所述基因组编码器(2010)包括:
比对器单元(201),其被配置为将所述读段与一个或多个参考序列进行比对,从而创建比对的读段,
构建的参考生成器单元(202),其被配置为产生构建的参考序列,
数据分类单元(204),其被配置为根据具有所述一个或多个预先存在的参考序列或构建的参考序列的指定匹配规则对所述比对的读段进行分类,从而创建比对的读段的类别(208);
一个或多个块编码单元(205-207),其被配置为通过根据所述比对的读段的类别选择所述描述符来将所述分类的比对的读段编码为描述符块,
多路复用器(2016),用于多路复用压缩的基因组数据和元数据。
在另一方面,基因组编码器还包括
参考序列转换单元(2019),其被配置为将所述预先存在的参考和数据类别(208)转换成转换的数据类别(2018)。
在另一方面,基因组编码器还包括
数据分类单元(204),其包含数据类别N、M和I的编码器,所述编码器配置有生成数据类别N、M和I的子类的阈值的向量。
在另一方面,基因组编码器还包括以下特征:所述参考转换单元(2019)对所有类别和子类的数据应用相同的参考转换(300)。
在另一方面,基因组编码器还包括以下特征:所述参考转换解码器(2019)对不同类别和子类的数据应用不同的参考转换(301、302、303)。
在另一方面,基因组编码器还包括适于执行前述编码方法的所有方面的特征。
本发明还提供了一种用于解压缩压缩的基因组流(211)的基因组解码器(218),所述基因组解码器(218)包括:
多路分解器(210),用于多路分解压缩的基因组数据和元数据,
解析装置(212-214),其被配置为将所述压缩的基因组流解析成描述符(215)的基因组块,
一个或多个块解码器(216-217),其被配置为将所述基因组块解码成核苷酸序列的分类的读段(2111),
基因组数据类别解码器(219),其被配置为选择性地解码一个或多个参考序列上的所述核苷酸序列的分类的读段,以便产生核苷酸序列的未压缩的读段。
在另一方面,基因组解码器还包括参考转换解码器(2113),其被配置为解码参考转换描述符(2112),并且产生要由基因组数据类别解码器(219)使用的转换的参考(2114)。
在另一方面,基因组解码器还包括,所述一个或多个参考序列被存储在压缩的基因组流(211)中。
在另一方面,基因组解码器还包括,所述一个或多个参考序列经由带外机制被提供给所述解码器。
在另一方面,基因组解码器还包括,在所述解码器处构建所述一个或多个参考序列。
在另一方面,基因组解码器还包括,一个或多个参考序列在所述解码器处由参考转换解码器(2113)转换。
本发明还提供了一种计算机可读介质,包括指令,当指令被执行时,使得至少一个处理器进行前述编码方法的所有方面。
本发明还提供了一种计算机可读介质,包括指令,当指令被执行时,使得至少一个处理器进行前述解码方法的所有方面。
本发明还提供了一种支持数据,存储根据执行前述编码方法的所有方面所编码的基因组。
具体实施方式
本发明中提到的基因组或蛋白质组序列包括例如但不限于核苷酸序列、脱氧核糖核酸(DNA)序列、核糖核酸序列和氨基酸(RNA)序列。尽管本文的描述相当详细地涉及核苷酸序列形式的基因组信息,但是应当理解,尽管有一些变化,但是如本领域技术人员将理解的,用于压缩的方法和系统也可以被实施用于其他基因组或蛋白质组序列。
基因组测序信息是由高通量测序(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或-)。
由测序机器所产生的核苷酸序列被称为“读段”。序列读段可以在几十到几千个核苷酸长之间。一些技术产生“配对”的序列读段,其中,一个读段来自一个DNA链,第二读段来自另一链。在基因组测序中,术语“覆盖”被用来表达序列数据相对于“参考序列”的冗余水平。例如,为了在人类基因组(32亿个碱基长)上达到30×的覆盖,测序机器将产生总共30×32亿个碱基,使得参考中的每个位置平均被“覆盖”30次。
在本公开中,参考序列是由测序机器所产生的核苷酸序列比对/映射的任何序列。序列的一个示例实际上可以是“参考基因组”,科学家作为物种的基因组的代表性示例所组装的序列。例如,GRCh37,基因组参考联盟人类基因组(版本37)来自纽约水牛城的十三名匿名志愿者。然而,参考序列也可以由合成序列组成,考虑到它们的进一步处理,该合成序列的构思和构建仅仅是为了提高读段的可压缩性。这在“U类描述符”以及“U类”和“HM类”的未映射的读段的“内部”参考的构造部分中更详细地描述,并且在图22和23中示出。
测序设备可能在序列读段中引入错误,诸如:
1.由于对调用任何特定碱基缺乏信心而跳过碱基调用的决定。这被称为未知碱基并且被标记为“N”(被表示为“n型”错配);
2.使用错误的符号(即,表示不同的核酸)来代表测序样品中实际存在的核酸;这通常被称为“替换错误”(被表示为“s类型”错配);
3.在一个序列读段中插入不指代任何实际存在的核酸的额外符号;这通常被称为“插入错误”(被表示为“i型”错配);
4.从一个序列读段中删除表示测序样品中实际存在的核酸的符号;这通常被称为“删除错误”(被表示为“d型”错配);
5.一个或多个片段重组合为单个片段,该单个片段不反映起始序列的真实情况;这通常导致比对器决定裁剪碱基(被表示为“c型”错配)。
文献中使用术语“覆盖”来量化可以由可用序列读段所覆盖的参考基因组或其部分的程度。覆盖被称为:
●当参考基因组的一些部分没有被任何可用序列读段映射时,部分(小于1X);
●当参考基因组的所有核苷酸被序列读段中存在的一个且只有一个符号映射时,单个(1X);
●当参考基因组的每个核苷酸被多次映射时,多个(2X、3X、NX)。
本发明旨在定义基因组信息表示格式,其中,相关信息是有效可访问和可运输的,并且冗余信息的权重被降低。
所公开发明的主要创新方面如下。
1.根据相对于参考序列的比对结果,序列读段被分类和划分为数据类。这样的分类和划分使得能够根据与比对结果和匹配精确性有关的标准来选择性地访问所编码的数据。
2.所分类的序列读段和相关联的元数据由描述符的同质块表示,以获得由低信息熵所表征的不同信息源。
3.用被适用于每个类的统计特征的不同源模型对每个分离的信息源建模的可能性以及为每个分离的可访问数据单元(访问单元)在每个读段类和每个描述符块内改变源模型的可能性。根据每个源模型的统计特征,采用适当的情境自适应概率模型以及相关联的熵编码器。
4.描述符块之间对应性和依赖性的定义,以便能够选择性地访问测序数据和相关联的元数据,如果不需要所有信息则不需要解码所有描述符块。
5.相对于“预先存在的”(也被称为“外部”)参考序列,或相对于通过对“预先存在的”参考序列应用适当的转换以减少描述符块信息源的熵而获得的“转换的”参考序列,对每个序列数据类和相关联的元数据块进行编码。所述描述符表示被划分成不同数据类的读段。在参考“预先存在的”参考或“转换的”“预先存在的”参考序列使用相对应的描述符对读段进行任何编码之后,各种错配的出现可以被用来定义对参考序列的适当转换,以便找到具有低熵并且实现更高的压缩效率的最终编码的表示。
6.一个或多个参考序列(也被称为“内部”参考,以区别于“预先存在的”参考序列,这里也被称为“外部”参考序列)的构建被用来对读段类进行编码,该读段类相对于不满足一组约束的预先存在的参考序列呈现一定程度的匹配精确性。设置这样的约束的目的是,以压缩的形式表示相对于“内部”参考序列所比对的读段类的编码成本和表示“内部”参考序列本身的成本,低于逐字编码未比对的读段类,或者使用不具有转换或具有转换的“外部”参考序列。
在下文中,将进一步详细描述上述各个方面。
根据匹配规则对序列读段进行分类
根据相对于一个或多个“预先存在的”参考序列的比对的匹配结果,由测序机器所生成的序列读段被本公开发明分类成六个不同的“类”。
当核苷酸的DNA序列相对于参考序列比对时,可以识别以下情况:
1.发现参考序列中的区域与序列读段匹配,没有任何错误(即,完美的映射)。这样的核苷酸序列被称为“完全匹配的读段”,或者被称为“P类”。
2.发现参考序列中的区域与序列读段匹配,具有错配的类型和数量,错配的类型和数量仅由生成读段的测序机器不能够调用任何碱基(或核苷酸)的位置的数量决定。这样的错配类型用“N”表示,该字母被用来指示未定义的核苷酸碱基。在本文档中,该错配类型被称为“n型”错配。这样的序列属于“N类”读段。一旦读段被分类为属于“N类”,则将匹配不精确性的程度限于给定的上限并且在被认为是有效匹配的和不是有效匹配的之间设置边界是有用的。因此,被分配给N类的读段还受到设置阈值(MAXN)的限制,该阈值定义读段可以包含的未定义的碱基(即,被称为“N”的碱基)的最大数量。这样的分类隐含地定义当参考相对应的参考序列时,属于N类的所有读段共享的所需最小匹配精确性(或最大错配程度),这构成了对压缩的数据应用选择性数据搜索的有用标准。
3.发现参考序列中的区域与序列读段匹配,具有错配的类型和数量,错配的类型和数量由生成读段的测序机器不能够调用任何核苷酸碱基(如果存在的话)的位置的数量(即“n型”错配)加上已经调用不同于参考中存在的碱基的错配的数量决定。这样的被称为“替换”的错配类型也被称为单核苷酸变异(SNV)或单核苷酸多态性(SNP)。在本文档中,该错配类型也被称为“s型”错配。序列读段被称为“M错配的读段”,并且被分配给“M类”。与在“N类”的情况下类似,对于属于“M类”的所有读段,将匹配不精确性的程度限于给定的上限并且在被认为是有效匹配的和不是有效匹配的之间设置边界是有用的。因此,被分配给M类的读段也通过定义一组阈值来限制,如果“n型”错配存在的话,则一个阈值用于“n型”(MAXN)错配的数量“n”,并且另一阈值用于替换“s”(MAXS)的数量。第三约束是由数量“n”和“s”的任何函数所定义的阈值,f(n,s)。这样的第三约束使得能够根据任何有意义的选择性访问标准来生成具有匹配不精确性的上限的类。例如,非限制性地,f(n,s)可以是(n+s)1/2或(n+s)或任何线性或非线性表达式,该表达式将边界设置为属于“M类”的读段所允许的最大匹配不精确性水平。当出于各种目的分析序列读段时,这样的边界构成将期望的选择性数据搜索应用于压缩的数据的非常有用的标准,因为除了被应用于一个类型或另一类型的简单阈值之外,有可能对“n型”错配和“s型”错配(替换)的数量的任何可能组合设置进一步的边界。
4.第四类是通过测序读段构成的,该测序读段在“插入”、“删除”(又名插入缺失)和“剪切”之间呈现至少一个任何类型的错配,如果存在的话,加上属于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公司)的情况下,已知两个读段被未知的可变长度序列分开,适当地考虑将整个对分类为单个数据类。读段和另一读段耦合,被称为其“配对”。
如果两个配对的读段属于相同的类,则分配给整个配对的类是显然的:对于任何类(即,P、N、M、I、U),整个配对被分配给相同的类。如果两个读段属于不同的类,但是都不属于“U类”,则根据以下表达式所定义的最高优先级将整个配对分配给类:
P<N<M<I
其中,“P类”具有最低优先级,并且“I类”具有最高优先级。
如果只有一个读段属于“U类”,并且其配对属于P类、N类、M类、I类中的任何一个,则第六类被定义为“HM类”,其代表“半映射的”。
这种特定类别的读段的定义的动机是,用于试图确定参考基因组中存在的间隙或未知区域(也称为鲜为人知或未知的区域)。这种区域通过使用可以被映射到已知区域上的配对读段在边缘映射配对来重建。然后,未映射的配对用于构建未知区域的所谓“重叠群”,如图28中所示。因此,仅提供对这种类型的读段对的选择性访问大大减少了相关联的计算负担,使得能够非常有效地处理由大量数据集所产生的这种数据,需要使用现有技术的解决方案完全检查这些数据集。
下表总结了被应用于读段的匹配规则,以便定义每个读段所属的数据类。规则在表的前五列中根据错配(n、s、d、i和c类型错配)的类型的存在与否来定义。第六列提供了关于每个错配类型的最大阈值以及可能错配类型的函数f(n,s)和w(n,s,d,i,c)的规则。
/>
表1.每个序列读段必须满足的错配类型和约束集被分类到本发明公开中所定义的数据类中。
序列读段数据类N、M和I的匹配规则划分为不同匹配精确性的子类
先前部分中所定义的类型为N、M和I的数据类可以进一步被分解为任意数量的不同子类,这些子类具有不同程度的匹配精确性。这种选项在提供更精细的粒度并且因此提供对每个数据类的高效得多的选择性访问方面是重要的技术优势。作为示例而非限制,为了将类别N划分成k个子类(子类别N1、…、子类别Nk),需要定义具有相对应的组件的向量MAXN1、MAXN2、…、MAXN(k-1)、MAXN(k),条件是MAXN1<MAXN2<…<MAXN(k-1)<MAXN,并且当为向量的每个元素评估时,将每个读段分配给满足表1中所指定的约束的最低等级的子类。这在图29中示出,其中,数据分类单元291包含P、N、M、I、U、HM类编码器和用于注释和元数据的编码器。N类编码器配置有阈值向量MAXN1至MAXNk292,N类编码器生成N数据的k个子类(296)。
在类型M和类型I的类别的情况下,通过分别为MAXM和MAXTOT定义具有相同属性的向量来应用相同的原理,并且使用每个向量组件作为阈值来检查函数f(n,s)和w(n,s,d,i,c)是否满足约束。与在类型为N的子类的情况下类似,分配给满足约束的最低子类。每个类别类型的子类的数量是独立的,并且任何细分的组合都是可以接受的。这在图29中示出,其中,M类编码器293和I类编码器294分别配置有阈值向量MAXM1至MAXMj和阈值向量MAXTOT1至MAXTOTh。两个编码器分别生成M数据的j子类(297)和I数据的h子类(298)。
当配对中的两个读段被分类在相同的子类中时,则该配对属于相同的子类。
当配对中的两个读段被分类为不同类别的子类时,则根据以下表达式,该配对属于较高优先级类别的子类:
N<M<I
其中,N的优先级最低,并且I的优先级最高。
当两个读段属于类别N或M或I中的一个类别的不同子类时,则根据以下表达式,该配对属于具有最高优先级的子类:
N1<N2<…<Nk
M1<M2<…<Mj
I1<I2<…<Ih
其中,最高索引具有最高优先级。
“外部”参考序列的转换
为被分类在类别N、M和I中的读段所找到的错配可以用于创建“转换的”参考,以用于更有效地压缩读段表示。
被分类为属于类别N、M或I(相对于被表示为RS0的“预先存在”(即“外部”)参考序列)的读段可以根据与“转换的”参考的实际错配的出现,相对于“转换的”参考序列RS1进行编码。例如,如果属于M类的readM in(被表示为M类的第i读段)包含相对于参考序列RSn的错配,则在“转换”之后,可以用A(Refn)=Refn+1获得readM in=readP i(n+1),其中,A是从参考序列RSn到参考序列RSn+1的转换。
图19示出了示例,可以如何通过修改对应于错配位置的碱基,将包含相对于参考序列1(RS1)的错配(属于M类)的读段转换为相对于从RS1所获得的参考序列2(RS2)的完全匹配的读段。这些读段保持被分类,并且在相同的数据类访问单元中与其他读段一起被编码,但是仅使用P类读段所需的描述符和描述符值来进行编码。这种转换可以被表示为:
RS2=A(RS1)
当在被应用于RS1时产生RS2的转换A的表示加上读段相对RS2的表示,对应于比M类的读段相对RS1的表示更低的熵时,传输转换A的表示和读段相对RS2的对应表示是有利的,因为实现了数据表示的更高压缩。
用于在压缩的比特流中传输的转换A的编码需要定义如下表中所定义的两个附件描述符。
图26示出了示例,如何应用参考转换来减少要被编码到映射的读段上的错配的数量。
必须注意,在一些情况下,转换被应用于参考:
●可能在应用转换之前参考参考时不存在的读段的表示中引入错配;
●可能修改错配的类型,读段可能包含A而不是G,而所有其他读段包含C而不是G,但是错配保持在相同的位置;
●不同的数据类和每个数据类的数据子集可能指代相同的“转换的”参考序列,或者指代通过对相同的预先存在的参考序列应用不同的转换而获得的参考序列。
图27还示出了示例,在应用参考转换并且使用“转换的”参考来表示读段之后,读段可以如何通过适当的描述符集(例如,使用P类的描述符来编码来自M类的读段)将编码的类型从数据类改变为另一数据类。例如,当转换改变与读段中实际存在的碱基中的读段的错配相对应的所有碱基时,发生这种情况,从而虚拟地将属于M类的读段(当参考原始的非“转换的”参考序列时)转换成P类的虚拟读段(当参考“转换的”参考时)。以下部分中提供用于每类数据的描述符集的定义。
图30示出了不同类别的数据可以如何使用相同的“转换的”参考R1=A0(R0)(300)来重新编码读段,或者不同的转换AN(301)、AM(302)、AI(303)可以分别被应用于每个类别的数据。
定义表示到描述符块的序列读段所需的信息
一旦用类别的定义完成了读段的分类,则进一步的处理包括定义一组不同的描述符,这些描述符表示当被表示为被映射到给定的参考序列上时能够重建读段序列的剩余信息。这些描述符的数据结构需要存储由解码引擎所使用的全局参数和元数据。这些数据在下表中所描述的基因组数据集报头中被结构化。数据集被定义为重建与单个基因组测序运行和所有后续分析相关的基因组信息所需的编码元素的集合。如果相同的基因组样品在两个不同的运行中测序两次,则将在两个不同的数据集中编码获得的数据。
/>
/>
表1-基因组数据集报头结构。
参考给定参考序列的序列读段(即,DNA片段)可以通过以下方式完全表达:
●参考序列上的起始位置(pos)
●表示读段是否必须被视为与参考相反的补码的标志(rcomp)。
●在配对读段的情况下到配对的距离(pair)。
●在测序技术产生可变长度读段的情况下,读段长度的值(len)。在读段长度恒定的情况下,与每个读段相关联的读段长度显然可以被省略,并且可以被存储在主文件报头中。
●对于每个错配:
○错配位置(N类的nmis,M类的snpp,I类的indp)
○错配类型(N类中不存在,M类中的snpt,I类中的indt)
●指示序列读段的特定特征的标志,例如,
○测序中具有多个片段的模板
○根据比对器被适当比对的每个片段
○未映射的片段
○模板中未映射的下一片段
○第一或最后片段的信号通知
○质量控制失败
○PCR或光学复制
ο二次比对
ο补充比对
●当存在时,软剪切的核苷酸串(I类中的indc)
●如果适用(描述符rtype),指示用于比对和压缩的参考(例如,U类的“内部”参考)的标志
●对于U类,描述符indc识别具有一组指定的匹配精确性约束的与“内部”参考不匹配的那些读段的部分(通常是边缘),。
●描述符ureads用于逐字编码不能被映射到任何可用参考上的读段,该参考是预先存在的(即,“外部的”,例如,实际的参考基因组)或“内部的”参考序列。
该分类创建了可以用来唯一地表示基因组序列读段的一组描述符(描述符)。下表总结了与“外部的”(即,“预先存在的”)或“内部的”(即,“构造的”)参考所比对的每一类读段所需的描述符。
P N M I U HM
pos X X X X X X
pair X X X X
rcomp X X X X X X
flags X X X X X X
rlen X X X X X X
nmis X
snpp X X
snpt X X
indp X X
indt X X
indc X X X
ureads X X
rtype X
rgroup X X X X X X
mmap X X X X X
msar X X X X X
mscore X X X X X
表2-每类数据的定义的描述符块。
表征属于P类的读段,并且在已经通过产生配对、一些标志和读段长度的测序技术所获得的情况下,可以仅通过位置、反向补码信息和配对之间的偏移来完全地重建属于P类的读段。
下一部分进一步详细说明对于P、N、M和I类,如何定义这些描述符,而对于U类,在以下部分中描述。
HM类仅适用于读段对,这是特殊情况,其中一个读段属于P、N、M或I类,另一读段属于U类。
位置描述符
在位置(pos)块中,只有第一编码的读段的映射位置作为绝对值被存储在参考序列上。所有其他位置描述符都采用表示相对于先前位置的差异的值。由读段位置描述符的序列所定义的信息源的这种建模通常以熵减少为特征,特别是对于生成高覆盖结果的测序处理。
例如,图1示出了在将第一比对的起始位置描述为参考序列上的位置“10000”之后,在位置10180开始的第二读段的位置如何被描述为“180”。对于高覆盖(>50x),位置向量的大多数描述符呈现非常高的低值出现率,低值诸如是0和1以及其他小整数。图1示出了如何在pos块中描述三个读段对的位置。
反向补码描述符
由测序技术所产生的读段对中的每个读段可以来自测序的有机样品的基因组链。然而,两个链中只有一个链被用作参考序列。图2示出了在读段对中,一个读段(读段1)可以来自一个链,并且另一读段(读段2)可以来自另一链。
当链1被用作参考序列时,读段2可以被编码为链1上相对应片段的反向补码。这在图3中示出。
在耦合的读段的情况下,四个是直接和反向补码配对的可能组合。这在图4中示出。rcomp块编码四种可能的组合。
相同的编码用于属于N、M、P和I类的读段的反向补码信息。为了能够选择性地访问不同的数据类,如表2中所示,在不同的块中编码属于四个类别的读段的反向补码信息。
配对信息描述符
配对描述符被存储在配对块中。当所采用的测序技术按对产生读段时,这种块存储编码重建原始读段对所需信息的描述符。尽管在本发明公开之日,绝大多数测序数据是通过使用生成配对读段的技术产生的,但并非所有技术都是如此。这就是为什么如果所考虑的基因组数据的测序技术没有生成配对的读段信息,则该块的存在对于重建所有测序数据信息是不必要的。
定义:
●配对:与读段对中的另一读段相关联的读段(例如,在前面的示例中,读段2是读段1的配对)
●配对距离:参考序列上的核苷酸位置的数量,该参考序列将第一读段中的一个位置(配对锚,例如,第一读段的最后核苷酸)与第二读段的一个位置(例如,第二读段的第一核苷酸)分开
●最可能配对距离(MPPD):这是以核苷酸位置的数量所表示的最可能配对距离。
●位置配对距离(PPD):PPD是以将一个读段与其在特定位置描述符块中存在的相应配对分开的读段的数量来表示配对距离的方式。
●最可能位置配对距离(MPPPD):是将一个读段与其在特定位置描述符块中存在的配对分开的最可能读段的数量。
●位置配对错误(PPE):被定义为MPPD或MPPPD与配对的实际位置之间的差异。
●配对锚:配对中第一读段的最后核苷酸的位置,被用作参考,以根据核苷酸位置的数量或读段位置的数量来计算配对的距离。
图5示出了如何计算读段对之间的配对距离。
配对描述符块是配对错误的向量,其被计算为相对于所定义的解码配对距离要被跳过以到达配对的第一读段的配对的读段的数量。
图6示出了如何计算配对错误的示例,作为绝对值和差分向量(以高覆盖的更低熵为特征)。
相同的描述符用于属于N、M、P和I类的读段的配对信息。为了能够选择性地访问不同的数据类,属于四个类别的读段的配对信息在不同的块中被编码,如图8(N类),图10、12和14(M类)以及图15和16(I类)中所示。
在读段被映射到不同参考序列上的情况下的配对信息
在将序列读段映射到参考序列上的处理中,将配对中的第一读段映射到一个参考序列(例如,染色体1)上,并且将第二读段映射到不同的参考序列(例如,染色体4)上,这并不罕见。在该情况下,上述配对信息必须通过与用于映射一个读段的参考序列有关的附加信息来整合。这是通过编码以下内容来实现的:
1.指示配对被映射到两个不同的序列上的保留值(标志)(不同的值指示读段1或读段2被映射到当前未编码的序列上)。
2.指代如表1所述在主报头结构中所编码的参考标识符的唯一参考标识符。
3.第三元素,该第三元素包含在点2所识别的参考上的映射信息,并且被表示为相对于最后编码的位置的偏移。
图7提供了该场景的示例。
在图7中,因为读段4没有被映射到当前编码的参考序列上,所以基因组编码器通过在配对块中构建附加描述符来表示该信息。在下面所示的示例中,配对2的读段4被映射到第4参考上,而当前编码的参考是第1参考。该信息使用3个组件被编码:
1)一个特殊的保留值被编码为配对距离(在该情况下,0xffffff)。
2)第二描述符提供了如主报头中所列出的参考ID(在该情况下,4)。
3)第三元素包含有关参考上的映射信息(170)。
N类读段的错配描述符
N类包括只存在“n型”错配的所有读段,在A、C、G或T碱基的位置,发现N被称为碱基。读段的所有其他碱基与参考序列完全匹配。
图8示出了如何:
将读段1中“N”的位置编码为:
●读段1中的绝对位置或
●作为在相同读段中相对于先前“N”的不同位置
将读段2中“N”的位置编码为
●读段2中的绝对位置+读段1长度或
●相对于先前N的不同位置
在nmis块中,每个读段对的编码由特殊的“分隔符”符号终止。
编码替换(错配或SNP)、插入和删除的描述符
替换被定义为在映射的读段中相对于在相同位置的参考序列中存在的核苷酸碱基不同的核苷酸碱基的存在。
图9示出了映射的读段对中的替换的示例。每个替换被编码为“位置”(snpp块)和“类型”(snpt块)。根据替换、插入或删除的统计发生,可以定义相关联的描述符的不同源模型,并且在相关联的块中编码所生成的符号。
源模型1:作为位置和类型的替换
替换位置描述符
类似于nmis块的值,计算替代位置,即,
在读段1中,替换被编码为
●作为读段1中的绝对位置或
●作为在相同读段中相对于先前替换的不同位置
在读段2中,替换被编码为
●作为读段2中的绝对位置+读段1长度或
●作为相对于先前替换的不同位置
图10示出了如何将替换(其中,在给定的映射位置,读段中的符号不同于参考序列中的符号)编码为
1.错配的位置
■相对于读段的开始或
■相对于先前错配(不同编码)
2.错配的类型,被表示为如图10中所示计算的代码
在snpp块中,每个读段对的编码由特殊的“分隔符”符号终止。
替换类型描述符
对于M类(和下一部分中所描述的I类),错配通过从参考中存在的实际符号到读段{A,C,G,T,N,Z}中存在的相对应替换符号的索引(从右向左移动)被编码。例如,如果比对的读段存在C而不是在参考中相同位置的T,则错配索引将被表示为“4”。解码处理读取所编码的描述符,在参考上给定位置的核苷酸,并且从左向右移动,以检索所解码的符号。例如,对于参考中存在G的位置所接收的“2”将被解码为“N”。图11示出了所有可能的替换和相应的编码符号。显然不同的和情境自适应的概率模型可以根据每个数据类的每个替换类型的统计属性被分配给每个替换索引,以最小化描述符的熵。
在采用IUPAC模糊代码的情况下,替换机制结果是完全相同的,但是替换向量被扩展为:S={A,C,G,T,N,Z,M,R,W,S,Y,K,V,H,D,B}。
图12提供了snpt块中替换类型的编码的示例。
图13中提供了当采用IUPAC模糊代码时替换编码的一些示例。图14中提供了替换索引的进一步示例。
插入和删除的编码
对于I类,错配和删除通过从参考中存在的实际符号到读段:{A,C,G,T,N,Z}中存在的相对应替换符号的索引(从右向左移动)被编码。例如,如果比对的读段存在C而不是在参考中相同位置的T,则错配索引将为“4”。在读段存在删除而参考中存在“A”的情况下,编码符号将为“5”。解码处理读取所编码的描述符,在参考上给定位置的核苷酸,并且从左向右移动,以检索所解码的符号。例如,对于参考中存在G的位置所接收的“3”将被解码为“Z”。
对于所插入的A、C、G、T、N,插入分别被编码为6、7、8、9、10。
图15示出了如何在I类的读段对中编码替换、插入和删除的示例。如用于错配的先前段落中所描述的,为了支持完整的IUPAC模糊代码集,替换向量S={A、C、G、T、N、Z}应该被替换为S={A、C、G、T、N、Z、M、R、W、S、Y、K、V、H、D、B}。在该情况下,插入代码需要具有不同的值,即,在替换向量具有16个元素的情况下,16、17、18、19、20。机制被示出在图16中。
源模型2:每个替换类型和缺失的一个块
对于一些数据统计,可以开发不同于先前部分中所描述的编码模型的编码模型,用于替换和缺失,从而得到熵较低的源。这样的编码模型是上述仅用于错配以及错配和缺失的技术的替代。
在该情况下,为每个可能的替换符号(没有IUPAC代码则为5,具有IUPAC代码则为16)定义一个数据块,加上一个删除块和4个插入块。为了解释简单,但不是对模型应用的限制,下面的描述将集中在不支持IUPAC代码的情况。
图17示出了每个块如何包含单一类型的错配或插入的位置。如果编码的读段对中不存在该类型的错配或插入,则在相对应的块中编码0。为了使得解码器能够开始本部分中所描述的块的解码处理,每个访问单元的报头包含表示要被解码的第一块的标志。在图18的示例中,要被解码的第一元素是C块中的位置2。当读段对中不存在给定类型的错配或缺失时,0被添加到相对应的块中。在解码侧,当每个块的解码指针指向0值时,解码处理移动到下一读段对。
附加表示标志的编码
上面所介绍的每个数据类(P、M、N、I)可能需要对关于编码读段性质的附加信息进行编码。该信息可以例如与测序实验(例如,指示一个读段的重复的概率)有关,或者可以表达读段映射(例如,配对中的第一或第二)的一些特征。在本发明的情境中,对于每个数据类,该信息在单独的块中被编码。这样的方法的主要优点是,只有在需要的情况下,并且只有在所需的参考序列区域中,才有可能选择性地访问该信息。使用这样的标志的其他示例是:
●配对的读段
●适当配对中的映射的读段
●未映射的读段或配对
●来自反向链的读段或配对
●配对中的第一/第二
●不是主要比对
●读段未通过平台/供应商质量检查
●读段是PCR或光学复制
●补充比对
U类的描述符以及“U类”和“HM类”的未映射的读段的“内部”参考的构建
在属于U类的读段或“HM类”的未映射的对的情况下,因为它们不能被映射到满足用于属于P、N、M或I类中的任何类的匹配精确性约束的指定集的任何“外部”参考序列,所以一个或多个“内部”参考序列被“构建”并且用于属于这些数据类的读段的压缩的表示。
多个方法可以构建适当的“内部”参考,例如但不限于:
●将未映射的读段划分成包含共享至少最小尺寸(签名)的公共连续基因组序列的读段的聚类。如图22中所示,每个聚类可以通过其签名被唯一地识别。
●以任何有意义的顺序(例如,字典顺序)对读段进行排序,并且将最后的N个读段用作用于编码N+1的“内部”参考。该方法被示出在图23中。
●对U类读段的子集进行所谓的“重新组装”,以便能够根据所指定的匹配精确性约束或新的约束集,对属于所述类别的读段的所有或相关子集进行比对和编码。
如果被编码的读段可以被映射到满足所指定的匹配精确性约束集的“内部”参考上,则使用可以是以下类型的描述符对压缩后重建读段所需的信息进行编码:
1.根据内部参考中的读段数量,匹配部分在内部参考上的起始位置(pos块)。相对于先前编码的读段,该位置可以被编码为绝对值或不同值。
2.起始位置相对于内部参考中相对应读段的开始的偏移(pair块)。例如,在恒定读段长度的情况下,实际位置是pos*长度+pair。
3.可能存在被编码为错配位置(snpp块)和类型(snpt块)的错配
4.与内部参考不匹配(或这样做,但错配的数量超过所限定的阈值)的那些读段部分(通常是由配对所识别的边缘)在indc块中被编码。如图24中所示,可以对所使用的内部参考的部分的边缘进行填充操作,以降低indc块中所编码的错配的熵。编码器可以根据正在被处理的基因组数据的统计属性来选择最合适的填充策略。可能的填充策略包括:
a.没有填充
b.根据当前所编码的数据中的频率所选择的恒定填充模式。
c.根据当前情境的统计属性的可变填充模式,当前情境根据最新的N编码读段被定义
特定类型的填充策略将由indc块报头中的特殊值来表示。
5.指示已经使用内部自生成、外部或无参考(rtype块)来编码读段的标志
6.逐字编码的读段(ureads)。
图24提供了这样的编码过程的示例。
图25示出了内部参考上未映射的读段的替代编码,其中,pos+pair描述符被有符号的pos替换。在该情况下,pos将根据参考序列上的位置,表示读段n的最左侧核苷酸位置相对于读段n-1的最左侧核苷酸位置的距离。
在U类读段存在可变长度的情况下,使用附加描述符rlen来存储每个读段长度。
该编码方法可以被扩展为支持每个读段的N个起始位置,从而可以将读段分割到两个或多个参考位置。这对于编码由那些产生非常长的读段(50K+碱基)的测序技术(例如,来自Pacific Bioscience)所生成的读段特别有用,这些读段通常存在由测序方法中的循环所生成的重复模式。相同的方法也可以用于编码嵌合序列读段,嵌合序列读段被定义为与基因组的两个不同部分比对、很少或没有重叠的读段。
上述方法可以明显地应用于简单的U类之外,并且可以应用于包含与读段位置有关的描述符的任何块(pos块)。
比对分数描述符
mscore描述符提供每个比对的分数。在本发明的情境中,它用于表示由基因组序列读段比对器所产生的每个读段的映射/比对分数。
分数使用指数和小数部分来表示。用于表示指数和小数部分的位数作为配置参数被传输。作为示例,但不是限制,表2示出了对于11位指数和52位小数部分,这在IEEE RFC754中是如何被指定的。
每个比对的分数可以由以下内容来表示:
●一个符号位(S)
●11位指数(E)
●53位小数(M)
表2.比对分数可以被表示为64位双精度浮点值
用于计算分数的碱基(基数)为10,因此:
分数=-1sx10ExM
读段组
在测序处理期间,可以产生不同类型的所测序的读段。作为示例,但不是限制,类型可以与不同的所测序的样品、不同的实验、测序机器的不同配置有关。根据所公开的发明,在测序和比对之后,通过名为rgroup的专用描述符来保存该信息。rgroup是与每个所编码的读段相关联的标签,并且使得解码装置能够在解码后将所解码的读段分组。
多个比对的描述符
为支持多个比对而指定以下描述符。在存在拼接的读段的情况下,本发明定义了将被设置为1的全局标志spliced_reads_flag。
mmap描述符
mmap描述符用于指示配对中的读段或最左侧的读段已经被比对了多少个位置。包含多个比对的基因组记录与一个多字节mmap描述符相关联。mmap描述符的前两个字节表示无符号的整数N,N将读段称为单个片段(如果编码的数据集中不存在拼接),或者称为对于多个可能的比对已经被拼接的读段的所有片段(如果数据集中存在拼接)。N值表示在该记录中为模板编码了多少个pos描述符的值。如下所述,N后面跟着一个或多个无符号的整数Mi
多个比对链
本发明中所描述的rcomp描述符用于使用本发明中所指定的语法来指定每个读段比对的链。
多个比对的分数
在多个比对的情况下,将如本发明中所指定的一个mscore分配给每个比对。
没有拼接的多个比对
如果访问单元中不存在拼接,则spliced_reads_flag是未设置的。
在配对末端测序中,mmap描述符由16位无符号的整数N以及一个或多个8位无符号的整数Mi组成,i假设值从1到完整的第一(这里是最左侧)读段比对的数量。对于每个第一读段比对,无论是否被拼接,Mi都用来表示多少个片段被用于比对第二读段(在该情况下,没有拼接,这等于比对的数量),然后多少个pair描述符的值被编码用于第一读段的比对。
Mi的值应该用于计算其指示第二读段的比对的数量。
Mi的特殊值(Mi=0)指示最左侧读段的第i比对与最右侧读段的比对配对,该最右侧读段已经与最左侧读段的第k比对配对,k<i(则没有检测到新的比对,这与上面的等式一致)。
例如,在最简单的情况下:
1.如果最左侧的读段有单个比对,最右侧的读段有两个替代的比对,则N将为1,M1将为2。
2.如果检测到最左侧的读段的两个替代的比对,但是仅检测到最右侧的读段的一个比对,则N将为2,M1将为1,M2将为0。
当Mi为0时,配对的相关联值应该链接到现有的第二读段比对;否则将引发语法错误,并且比对被视为中断。
示例:如前所述,如果第一读段有两个映射位置,第二读段仅有一个映射位置,则N为2,M1为1,M2为0。如果接下来是整个模板的另一替代的次要映射,则N将为3,M3将为1。
图39示出了在没有拼接的多个比对的情况下,N、P和Mi的含义,并且示出了如何使用pos、pair和mmap描述符来编码多个比对信息。
相对于图40,以下内容适用:
●最右侧的读段具有个比对
●当最左侧读段的第i比对与最右侧读段的比对配对时,Mi的一些值可以是=0,该最右侧读段的比对已经与最左侧读段的第k比对配对,k<i
●pair描述符的一个保留值可以被呈现给属于其他AU范围的信号比对。
如果存在,则它总是当前记录的第一pair描述符
具有拼接的多个比对
如果数据集用拼接的读段来编码,则msar描述符能够表示拼接长度和链。
在已经解码了mmap和msar描述符后,解码器知道已经编码了多少个读段或读段对以表示多个映射,以及多少个片段组成每个读段或读段对映射。这在图41和图42中被示出。
参考图41,以下内容适用:
●最左侧读段具有带N个拼接的N1个比对(N1≤N)。
●N表示最左侧读段的所有比对中存在的拼接的数量,并且它被编码为mmap描述符的第一值。
●最右侧的读段具有个拼接,其中,Mi是最右侧读段的拼接的数量,该最右侧读段与最左侧读段的第i比对配对关联(1≤i≤N1)。换言之,P表示最右侧读段的拼接的数量,并且使用mmap描述符的第一值之后的N个值来计算。
●N1和N2表示第一和第二读段的比对的数量,并且使用msar描述符的N+P个值来计算。
参考图42,以下内容适用:
●最左侧具有带N个拼接的N1个比对(N1≤N),如果N1=N,并且N2=P,则不存在拼接;
●最右侧读段具有个拼接,tj 1≤j≤P和N2(N2≤P)个比对;
●pair描述符的数量可以被计算为NP=Max(N1,P)+M0,其中,
○M0是值为0的Mi的数量
○在一个特殊的pair描述符指示在其他AU中存在比对的情况下,则NP必须增加1。
比对分数
mscore描述符允许表示比对的映射分数。在单端测序中,它将具有每个模板的N1值;在配对端测序中,它将具有整个模板的每个比对的值(第一读段的不同比对的数量可能+进一步的第二读段比对的数量,即,当Mi-1>0时)。
分数=Max(N1,N2)+M0
其中,M0表示Mi=0的总数。
在本发明中,一个以上的分数值可以与每个比对相关联。比对的数量由配置参数as_depth表示。
没有拼接的多个比对的描述符
表3.在没有拼接的多个比对的情况下,确定在一个基因组记录中表示多个比对所需的描述符的数量
具有拼接的多个比对的描述符
表4示出了在具有拼接的多个比对的情况下,确定在一个基因组记录中表示多个比对所需的描述符的数量。
表4.用于表示多个比对和相关联的分数的描述符。
不同序列上的多个比对
可能发生的情况是,比对处理找到了与主映射所在的序列不同的另一参考序列的替代映射。
对于唯一比对的读段对,当例如与另一染色体上的配对嵌合比对时,应该使用pair描述符来表示绝对读段位置。pair描述符应该用于表示包含相同模板进一步比对的下一记录的参考和位置。最后一个记录(例如,如果替代映射在3个不同的AU中被编码,则为第三个)应该包含第一记录的参考和位置。
在一对中最左侧读段的一个或多个比对存在于与当前编码的AU有关的序列不同的参考序列上的情况下,则保留值用于pair描述符。保留值之后是参考序列标识符和下一AU中所包含的所有序列中最左侧比对的位置(即,该记录的pos描述符的第一解码值)。
具有插入、删除、未映射部分的多个比对
当替代的次映射不保留序列被比对的参考区域的邻接性时,可能无法重建由比对器所生成的精确映射,因为实际序列(以及与诸如替换或插入缺失的错配有关的描述符)仅被编码用于主比对。msar描述符应该用于表示在参考序列包含插入缺失和/或软剪切的情况下,次比对如何映射在参考序列上。如果msar由次比对的特殊符号“*”表示,则解码器将从主比对和次比对映射位置重建次比对。
msar描述符
msar(多个片段比对记录)描述符支持拼接的读段和包含插入缺失或软剪切的替换的次比对。
msar旨在传达以下信息:
●映射的片段长度
●用于次比对和/或拼接的读段的不同映射邻接性(即,插入、删除或剪切的碱基的存在)
msar用于下面描述的扩展CIGAR字符串的语法加上表5中所描述的附加符号。
表5.除了表6中所描述的语法之外,用于msar描述符的特殊符号。
扩展cigar语法
该部分为与序列相关联的字符串以及有关的错配、插入缺失、剪切的碱基和关于多个比对和拼接的读段的信息指定了扩展的CIGAR(E-CIGAR)语法。
表6中列出了本发明中所描述的编辑操作。
/>
表6.MPEG-G E-CIGAR字符串的语法。
源模型、熵编码器和编码模式
对于本发明中所公开的基因组数据结构的每个数据类、子类和相关联的描述符块,可以根据由每个块所携带的数据或元数据的特定特征及其统计特性采用不同的编码算法。“编码算法”必须旨在作为描述符块的特定“源模型”与特定“熵编码器”的关联。可以指定和选择特定的“源模型”,以在源熵最小化方面获得最有效的数据编码。熵编码器的选择可以由编码效率考虑和/或概率分布特征以及相关联的实施问题来驱动。特定“编码算法”,也被称为“编码模式”的每个选择可以被应用于与整个数据集的数据类或子类相关联的整个“描述符块”,或者不同的“编码模式”可以被应用于被划分为访问单元的描述符的每个部分。
与编码模式相关联的每个“源模型”由以下内容来表征:
●由每个源所发出的描述符的定义(即,被用来表示一类数据的描述符集,该类数据诸如是读段位置、读段配对信息、相对于如表2中所定义的参考序列的错配)。
●相关联的概率模型的定义。
●相关联的熵编码器的定义。
进一步的优点
将序列数据分类到定义的数据类和子类,通过由单个单独的数据源(例如,距离、位置等)对描述符序列建模来表征,允许利用较低的信息源熵来实施有效的编码模式。
本发明的另一优点是可以只访问兴趣数据类型的子集。例如,基因组学中最重要的一个应用在于,发现基因组样品相对于参考(SNV)或群体(SNP)的差异。如今,这样类型的分析需要处理完整的序列读段,而通过采用由本发明所公开的数据表示,错配已经被隔离成仅一到三个数据类(取决于对考虑“n型”和“i型”错配的兴趣)。
进一步的优点是,当发布新的参考序列时,或者当对已经映射的数据进行重新映射(例如,使用不同的映射算法)以获得新的比对时,可以进行从参考特定“外部”参考序列到另一不同“外部”参考序列所压缩的数据和元数据的高效代码转换。
图20示出了根据本发明的原理的编码设备207。编码设备207接收例如由基因组测序设备200所产生的原始序列数据209作为输入。基因组测序设备200在本领域中是已知的,就像Illumina HiSeq 2500或者Thermo-Fisher Ion Torrent装置。原始序列数据209被馈送到比对器单元201,比对器单元201通过将读段与参考序列2020比对来准备用于编码的序列。或者,专用模块202可以用于通过使用如本文档中“对于未映射的U类读段构建内部参考”和“HM类”部分中所述的不同的策略从可用读段生成参考序列。在已经被参考生成器202处理之后,读段可以被映射到所获得的较长序列上。然后,数据分类模块204对所比对的序列进行分类。然后,参考转换的进一步的步骤被应用于参考,以便减少由数据分类单元204所生成的数据的熵。这意味着将外部参考2020处理成参考转换单元2019,参考转换单元2019产生转换的数据类2018和参考转换描述符2021。然后,所转换的数据类2018与参考转换描述符2021一起被馈送到块编码器205-207。然后,基因组块2011被馈送到算术编码器2012-2014,算术编码器2012-2014根据由块所携带的数据或元数据的统计属性对块进行编码。结果是基因组流2015。
图21示出了根据本公开的原理的解码设备218。解码设备218从网络或存储元件接收多路复用的基因组比特流2110。多路复用的基因组比特流2110被馈送到多路分解器210,以产生单独的流211,然后,单独的流211被馈送到熵解码器212-214,以产生基因组块215和参考转换描述符2112。所提取的基因组块被馈送到块解码器216-217,以进一步将块解码成数据类,并且参考转换描述符被馈送到参考转换单元2113。类解码器219进一步处理基因组描述符2111和所转换的参考2114,并且合并结果,以产生序列的未压缩读段,然后,该读段可以进一步以本领域已知的格式被存储,该格式例如是文本文件或zip压缩文件、或FASTQ或SAM/BAM文件。
类解码器219能够通过利用由一个或多个基因组流所携带的关于原始参考序列的信息和所编码的比特流中所携带的参考转换描述符2112来重建原始基因组序列。在参考序列不是由基因组流传输的情况下,则它们必须在解码侧可用,并且可由类解码器访问。
此处所公开的发明技术可以用硬件、软件、固件或其任意组合来实现。当以软件实现时,这些可以被存储在计算机介质上并且由硬件处理单元进行。硬件处理单元可以包括一个或多个处理器、数字信号处理器、通用微处理器、专用集成电路或其他分立逻辑电路。
本公开的技术可以在各种装置或设备中实现,包括移动电话、台式计算机、服务器、平板电脑和类似装置。
文件格式:通过使用主索引表来选择性访问基因组数据的区域
为了支持对所比对的数据的特定区域的选择性访问,本文档中所描述的数据结构实现了被称为主索引表(MIT)的索引工具。这是多维数组,包含相关联的参考序列上特定读段映射的轨迹。MIT中所包含的值是每个pos块中第一读段的映射位置,使得支持对每个访问单元的非顺序访问。MIT包含每类数据(P、N、M、I、U和HM)和每个参考序列的一个部分。MIT被包含在所编码的数据的基因组数据集报头中。图31示出了基因组数据集报头的结构,图32示出了MIT的一般视觉表示,图33示出了所编码的读段的P类的MIT的示例。
图33所示的MIT中所包含的值用于直接访问压缩域中的兴趣区域(以及相对应的AU)。
例如,参考图33,如果需要访问参考2上的位置150,000和250,000之间所包括的区域,则解码应用将跳到MIT中的第二参考,并且将寻找两个值k1和k2,使得k1<150,000并且k2>250,000。其中,k1和k2是从MIT所读取的两个索引。在图33的示例中,这将导致MIT的第二向量的第三和第四位置。然后,解码应用将使用这些返回值从如下一部分中所述的pos块本地索引表中获取适当数据的位置。
连同指向包含属于上述四类基因组数据的数据块的指针,MIT可以被用作在其生命周期期间被添加到基因组数据的附加元数据和/或注释的索引。
本地索引表
每个基因组数据块以被称为本地报头的数据结构为前缀。本地报头包含块的唯一标识符、每个参考序列的访问单元计数器的向量、本地索引表(LIT)和可选的一些块特定元数据。LIT是指向块有效载荷中属于每个访问单元的数据的物理位置的指针的向量。图34描述了一般块报头和有效载荷,其中,LIT用于以非顺序方式访问所编码的数据的特定区域。
在先前的示例中,为了访问在第2参考序列上所比对的读段的区域150,000到250,000,解码应用从MIT检索位置3和4。解码处理应该使用这些值来访问LIT相对应部分的第3和第4元素。在图35中所示的示例中,块报头中所包含的总访问单元计数器用于跳过与参考1(示例中为5)有关的AU有关的LIT索引。因此,包含所编码的流中所请求的AU的物理位置的索引被计算为:
属于所请求的AU的数据块的位置=属于要被跳过的参考1的AU的数据块+使用MIT所检索的位置,即,
第一块位置:5+3=8
最后块位置:5+4=9
使用被称为本地索引表的索引机制所检索的数据块是所请求的访问单元的一部分。
图36示出了MIT表中所包含的块如何对应于每个数据类或子类的LIT的块。
图37示出了使用MIT和LIT所检索的数据块如何组成一个或多个如下面部分中所定义的访问单元。
在本发明的实施例中,LIT可以被集成为MIT的子结构。这样的方法的优点是在所压缩的文件的顺序解析的情况下访问所索引的数据的速度。如果LIT被集成在文件报头的MIT中,则解码装置只需要解析一小部分数据,以在选择性访问的情况下检索所请求的压缩信息。对于本领域技术人员来说,在网络上流式传输的情况下,当MIT和LIT中所包含的索引信息将在第一数据块之间被传送时,另一优点是显而易见的,因此使得接收装置能够在整个数据传输完成之前进行诸如分类和选择性访问的操作。
访问单元
在数据类中被分类并且在压缩或未压缩的块中被结构化的基因组数据被组织成不同的访问单元。
基因组访问单元(AU)被定义为基因组数据的部分(以压缩或非压缩的形式),该部分重建核苷酸序列和/或相关元数据、和/或DNA/RNA的序列(例如,虚拟参考)和/或由基因组测序机器和/或基因组处理装置或分析应用所生成的注释数据。图37中提供了访问单元的示例。
访问单元是可以通过仅使用全局可用数据(例如,解码器配置)独立于其他访问单元或者通过使用被包含在其他访问单元中的信息来进行解码的数据块。
访问单元的区别在于:
●类型,表征它们携带的基因组数据和数据集的性质以及可以访问它们的方式,
●顺序,为属于相同类型的访问单元提供唯一顺序。
任何类型的访问单元都可以被进一步分类为不同的“类别”。
以下是不同类型的基因组访问单元的定义的非穷尽列表:
1)类型0的访问单元不需要参考来自其他访问单元的任何信息来进行访问或者解码和访问。由它们包含的数据或数据集所携带的全部信息可以由解码装置或处理应用独立地读取和处理。
2)类型1的访问单元包含参考由类型0的访问单元所携带的数据的数据。读取或解码和处理被包含在类型1的访问单元中的数据需要访问一个或多个类型0的访问单元。类型1的访问单元编码与“P类”序列读段有关的基因组数据。
3)类型2的访问单元包含参考由类型0的访问单元所携带的数据的数据。读取或解码和处理被包含在类型2的访问单元中的数据需要访问一个或多个类型0的访问单元。类型2的访问单元编码与“N类”的序列读段有关的基因组数据。
4)类型3的访问单元包含参考由类型0的访问单元所携带的数据的数据。读取或解码和处理被包含在类型3的访问单元中的数据需要访问一个或多个类型0的访问单元。类型3的访问单元编码与“M类”的序列读段有关的基因组数据。
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的访问单元不需要被排序,并且不需要以有序的方式被存储和/或传输(技术优势:并行处理/并行流)。
图37示出了访问单元如何由报头以及一个或多个同质数据的块组成。每个块可以由一个或多个块组成。每个块包含多个数据包,并且数据包是上面所介绍的描述符的结构化序列,以表示例如读段位置、配对信息、反向补码信息、错配位置和类型等。
每个访问单元在每个块中可以具有不同数量的数据包,但是在访问单元内,所有块具有相同数量的数据包。
每个数据包可以通过3个标识符X Y Z的组合来识别,其中:
●X识别它所属的访问单元
●Y识别它所属的块(即,它封装的数据类型)
●Z是表达数据包相对于相同块中其他数据包的顺序的标识符
图38示出了访问单元和数据包标签的示例,其中,AU_T_N是具有标识符N的类型T的访问单元,该标识符N可以暗示或者不暗示根据访问单元类型的顺序的概念。标识符用于唯一地将一个类型的访问单元与那些需要完全解码所携带的基因组数据的其他类型的访问单元相关联。
根据不同的测序处理,任何类型的访问单元都可以进一步在不同的“类别”中被分类和标记。例如,但不限于,可以在以下情况下进行分类和标记:
1.在不同时间对相同生物体进行测序(访问单元包含具有“时间”内涵的基因组信息),
2.对相同生物体的不同性质的有机样品(例如,人类样品的皮肤、血液、毛发)进行测序。这些是具有“生物”内涵的访问单元。

Claims (32)

1.一种用于编码基因组序列数据的计算机实施的方法,其特征在于,所述基因组序列数据包括核苷酸序列的读段,所述方法包括以下步骤:
将所述读段与一个或多个第一参考序列进行比对,从而创建比对的读段,
根据具有所述一个或多个第一参考序列的指定匹配规则将所述比对的读段分类成不同的类别,从而创建比对的读段的类别,所述分类包括:
当相对于用于映射的所述参考序列在映射的读段中不存在错配时,将所述参考序列中没有任何错配的基因组读段识别为第一类别,P类;
当仅在测序机器不能调用任何“碱基”的位置发现错配,并且每个读段中的所述错配的数量不超过给定阈值时,将基因组读段识别为第二类别,N类;
当在所述测序机器不能调用任何“碱基”,被称为“n型”错配,和/或调用与所述参考序列不同的“碱基”,被称为“s型”错配,的位置发现错配,并且所述错配的数量不超过所述“n型”、“s型”错配数量的给定阈值和从给定函数,f(n,s)获得的阈值时,将基因组读段识别为第三类别,M类;
当基因组读段可能具有相同类型的所述第三类别,M类错配,以及另外以下类型的至少一个错配:插入,i型、删除,d型、软剪切或硬剪切,c型,并且其中,所述每个类型的错配的数量不超过相对应的给定阈值和由给定函数,w(n,s,i,d,c)提供的阈值时,将基因组读段识别为第四类别,I类;
将所述分类的比对的读段编码为用于相应类别且在相应类别内同质的多个描述符块,
其中,将所述分类的比对的读段编码为多个描述符块包括,根据所述比对的读段的类别选择所述描述符,
用报头信息来构建所述描述符块,从而创建连续的访问单元;
其中,使用用于关于映射位置的信息的描述符块、用于关于链特异性,即,所述读段序列来自的DNA链的信息的描述符块以及用于关于序列读段的特征的信息的“flags”来构建所述第一类别,P类访问单元;并且其中,在所述P类访问单元中,使用相应描述符块来编码配对端读段的配对信息;
其中,使用所述第一类别,P类访问单元的相同描述符块加上用于关于未知碱基的所述位置的所述信息的描述符块来构建所述第二类别,N类访问单元;
其中,使用所述第一类别,P类访问单元的相同描述符块加上用于关于替换的位置和类型的信息的描述符块来构建所述第三类别,M类访问单元;
其中,使用所述第一类别,P类访问单元的相同描述符块加上用于关于替换、插入、删除和剪切的碱基的位置和类型的信息的描述符块来构建所述第四类别,I类访问单元;
所述方法还包括:
将在所述第一至所述第四类别,P、N、M、I类中没有找到任何分类的基因组读段识别为第五类别,U类,
使用至少一些所述第五类别的所述读段来构建一组第二参考序列,
将所述第五类别的所述读段与所述一组第二参考序列进行比对,
基于相对于所述第二参考序列的指定的匹配精确性约束,将所述第五类别的所述读段编码为相应描述符,
用报头信息来构建所述相应描述符,从而创建第五类别的访问单元。
2.根据权利要求1所述的方法,其特征在于,使用以下中的一个或多个来构建所述第五类别的所述访问单元:
用于关于所述映射位置的信息的描述符块;
用于关于所述链特异性,即,所述读段序列来自的DNA链的信息的描述符块以及用于关于所述序列读段的特征的信息的“flags”;并且其中,使用相应描述符块来编码配对端读段的配对信息;
用于关于替换的所述位置和类型的信息的描述符块;
用于关于与所述第二参考序列不匹配的所述读段的部分的信息的描述符块,
逐字编码无法映射到任何参考序列上的所述读段的描述符块。
3.根据权利要求1或2所述的方法,其特征在于:
要被编码的所述基因组序列的所述读段是配对的;以及
所述分类还包括将基因组读段识别为第六“HM类”,所述第六“HM类”包括所有读段对,其中,一个读段属于P、N、M或I类,另一读段属于“U类”。
4.根据权利要求3所述的方法,其特征在于,还包括以下步骤:
识别两个配对读段是否被分类在相同的所述类别,P、N、M、I或U类中,然后将所述对分配给所述相同的识别的类别,
识别所述两个配对读段是否被分类在不同的类别中,并且在它们都不属于“U类”的情况下,将所述读段对分配给具有根据以下表达式定义的最高优先级的所述类别:
P<N<M<I
其中,“P类”具有最低优先级,“I类”具有最高优先级;以及
识别所述两个配对读段中是否只有一个已经被分类为属于“U类”,并且将所述读段对分类为属于所述“HM类”序列。
5.根据权利要求4所述的方法,其特征在于,根据通过所述“n型”错配的数量(292)、所述函数f(n,s)(293)和所述函数w(n,s,i,d,c)(294),分别为每个N、M和I类定义的阈值的向量(292,293,294),读段N、M、I的每个类别被进一步分成两个或多个子类(296,297,298)。
6.根据权利要求4所述的方法,其特征在于,
还包括以下步骤:
识别所述两个配对读段是否被分类在相同的子类中,然后将所述对分配给所述相同的子类,
识别所述两个配对读段是否被分类为不同类别的子类,然后根据以下表达式将所述对分配给属于较高优先级类别的所述子类:
N<M<I
其中,N具有最低优先级,I具有最高优先级;以及
识别所述两个配对读段是否被分类在相同的类别中,并且这样的类别是N或M或I,但是被分类在不同的子类中,然后根据以下表达式将所述对分配给具有最高优先级的所述子类:
N1<N2<…<Nk
M1<M2<…<Mj
I1<I2<…<Ih
其中,最高索引具有最高优先级。
7.根据权利要求6所述的方法,其特征在于:
通过“pos”描述符块来编码关于每个读段的所述映射位置的信息;
通过rcomp描述符块来编码关于每个读段的所述链特异性,即,所述读段序列来自的DNA链的信息;以及
通过“pair”描述符块来编码配对端读段的配对信息。
8.根据权利要求7所述的方法,其特征在于,通过“flags”描述符块来编码附加比对信息,所述附加比对信息诸如是所述读段是否被映射在适当的对中、所述读段未通过平台/供应商质量检查、所述读段是PCR或光学复制、或者所述读段是补充的比对。
9.根据权利要求8所述的方法,其特征在于,通过“nmis”描述符块来编码关于未知碱基的信息。
10.根据权利要求9所述的方法,其特征在于,
通过“snpp”描述符块来编码关于替换的所述位置的信息;
通过“snpt”描述符块来编码关于替换的所述类型的信息。
11.根据权利要求10所述的方法,其特征在于:
通过“indp”描述符块来编码关于替换、插入或删除类型的错配的所述位置的信息;
通过“indt”描述符块来编码关于诸如替换、插入或删除的错配的所述类型的信息;
通过“indc”描述符块来编码关于映射的读段的剪切的碱基的信息。
12.根据权利要求11所述的方法,其特征在于:
通过“ureads”描述符块来编码关于未映射的读段的信息;
通过“rtype”描述符块来编码关于用于编码的参考序列的所述类型的信息;
通过“mmap”描述符块来编码关于所述映射的读段的多个比对的信息;
通过“msar”描述符块和“mmap”描述符块来编码关于相同读段的拼接的比对和多个比对的信息;
通过“mscore”描述符块来编码关于读段比对分数的信息;以及
通过“rgroup”描述符块来编码关于读段所属的所述组的信息。
13.根据权利要求12所述的方法,其特征在于,使用用于所述映射的读段的I类访问单元的相同描述符块,并且使用用于所述未映射的读段的所述“ureads”描述符块,来构建HM类访问单元。
14.根据权利要求13所述的方法,其特征在于:
使用P类访问单元的相同描述符块加上用于关于替换、插入、删除和剪切的碱基的位置和类型的所述信息的所述“indp”、“indt”和“indc”描述符块来构建I类访问单元,其中,使用所述“mmap”和“msar”描述符块来传送关于多个比对的信息。
15.根据权利要求14所述的方法,其特征在于,
使用扩展cigar字符串来传送关于拼接的比对的信息,包括:
符号=表示匹配碱基,
符号+表示插入,
符号-表示删除,
符号/表示前向链上的拼接,
符号%表示反向链上的拼接,
符号*表示无方向的拼接,
来自用于DNA的IUPAC代码的文本字符表示替换,
符号(n)表示n个软剪切的碱基,其中,n是整数,
符号[n]表示n个硬剪切的碱基,其中,n是整数。
16.根据权利要求15所述的方法,其特征在于:
所述描述符块包括“主索引表”,所述主索引表包含用于比对的读段的每个类别和子类的一个部分,所述部分包括每个类别或子类数据的每个访问单元的第一读段的所述一个或多个参考序列上的所述映射位置;
联合编码所述“主索引表”和所述访问单元数据。
17.根据权利要求16所述的方法,其特征在于:
所述描述符块还包括关于所使用的参考的所述类型的信息,即,预先存在或构建的,以及不映射到所述参考序列上的所述读段的片段。
18.根据权利要求17所述的方法,其特征在于,
所述参考序列首先通过应用替换、插入、删除和剪切被转换为不同的参考序列,然后将所述分类的比对的读段编码为参考所述转换的参考序列的多个描述符块。
19.根据权利要求18所述的方法,其特征在于:
相同的转换被应用于用于所有数据类别的所述参考序列;或者
不同的转换被应用于用于每个数据类别的所述参考序列,其中,所述参考序列转换被编码为描述符块,并且用报头信息来构建,从而创建连续的访问单元。
20.根据权利要求15所述的方法,其特征在于,不同的转换被应用于用于每个数据类别的所述参考序列,其中,所述参考序列转换被编码为描述符块,并且用报头信息来构建,从而创建连续的访问单元,并且其中,将所述分类的比对的读段和有关的参考序列转换编码成多个描述符块包括,将源模型和熵编码器与每个描述符块相关联的步骤。
21.根据权利要求20所述的方法,其特征在于,所述熵编码器是情境自适应算术编码器、可变长度编码器或golomb编码器中的一个。
22.一种用于解码编码的基因组数据的计算机实施的方法,其特征在于,包括以下步骤:
解析包含所述编码的基因组数据的访问单元,以通过使用报头信息来提取多个描述符块,以及
解码所述多个描述符块,以根据相对于一个或多个参考序列定义读段的分类的匹配规则来提取读段,
其中,如果访问单元是第一、第二、第三或第四类别,则所述描述符块根据指定匹配规则描述所述读段相对于第一参考序列的所述匹配,
其中,当访问单元表示相对于用于映射的所述第一参考序列不存在错配的基因组读段时,它是所述第一类别,P类;
其中,当访问单元表示仅在测序机器不能调用任何“碱基”的位置发现错配,并且每个读段中的所述错配的数量不超过给定阈值的基因组读段时,它是所述第二类别,N类;
其中,当访问单元表示在所述测序机器不能调用任何“碱基”,被称为“n型”错配,和/或调用与所述参考序列不同的“碱基”,被称为“s型”错配,的位置发现错配,并且所述错配的数量不超过所述“n型”、“s型”错配数量的给定阈值和从给定函数,
f(n,s)获得的阈值的基因组读段时,它是所述第三类别,M类;
其中,当访问单元表示可能具有相同类型的所述第三类别,M类错配,以及另外以下类型的至少一个错配:插入,i型、删除,d型、软剪切或硬剪切,c型,并且其中,所述每个类型的错配的数量不超过相对应的给定阈值和由给定函数,w(n,s,i,d,c)提供的阈值的基因组读段时,它是所述第四类别,I类;
其中,使用用于关于映射位置的信息的描述符块、用于关于链特异性,即,所述读段序列来自的DNA链的信息的描述符块以及用于关于序列读段的特征的信息的“flags”来构建所述第一类别,P类访问单元;并且其中,在所述P类访问单元中,使用相应描述符块来编码配对端读段的配对信息;
其中,使用所述第一类别,P类访问单元的相同描述符块加上用于关于未知碱基的所述位置的所述信息的描述符块来构建所述第二类别,N类访问单元;
其中,使用所述第一类别,P类访问单元的相同描述符块加上用于关于替换的位置和类型的信息的描述符块来构建所述第三类别,M类访问单元;
其中,使用所述第一类别,P类访问单元的相同描述符块加上用于关于替换、插入、删除和剪切的碱基的位置和类型的信息的描述符块来构建所述第四类别,I类访问单元;
其中,在访问单元是第五类别中,所述描述符块根据指定匹配规则描述所述读段相对于第二参考序列的所述匹配。
23.根据权利要求22所述的方法,其特征在于,使用以下中的一个或多个来构建所述第五类别的所述访问单元:
用于关于所述映射位置的信息的描述符块;
用于关于所述链特异性,即,所述读段序列来自的DNA链的信息的描述符块以及用于关于所述序列读段的特征的信息的“flags”;并且其中,使用相应描述符块来编码配对端读段的配对信息;
用于关于替换的所述位置和类型的信息的描述符块;
用于关于与所述第二参考序列不匹配的所述读段的部分的信息的描述符块,
逐字编码无法映射到任何参考序列上的所述读段的描述符块。
24.根据权利要求23所述的方法,其特征在于,还包括解码主索引表,所述主索引表包含每个读段类别的一个部分和相关联的相关映射位置。
25.根据权利要求24所述的方法,其特征在于,还包括:
解码与所使用的参考的所述类型有关的信息:预先存在的、转换的或构建的。
26.根据权利要求25所述的方法,其特征在于,还包括:
解码与要被应用于所述预先存在的参考序列的一个或多个转换有关的信息;其中,所述描述符块被熵解码。
27.根据权利要求25所述的方法,其特征在于:
通过解码类型为“pos”、“rcomp”、“flags”和“rlen”的描述符块来获得P类读段,
通过解码类型为“pos”、“rcomp”、“flags”、“rlen”和“nmis”的描述符块来获得N类读段,
通过解码类型为“pos”、“rcomp”、“flags”、“rlen”、“snpp”和“snpt”的描述符块来获得M类读段,
通过解码类型为“pos”、“rcomp”、“flags”、“rlen”、“indp”、“indt”和“indc”的描述符块来获得I类读段,
通过解码类型为“pos”、“rcomp”、“flags”、“rlen”、“snpp”、“snpt”、“indc”、“ureads”和“rtype”的描述符块来获得U类读段。
28.根据权利要求27所述的方法,其特征在于,也通过解码类型为“pair”的描述符块来获得P、N、M和I类的配对读段,通过解码类型为“pos”、“rcomp”、“flags”、“rlen”、“indp”、“indt”、“indc”和“ureads”的描述符块来获得HM类的配对读段。
29.一种用于编码基因组序列数据的基因组编码器,其特征在于,被配置为进行根据权利要求1至21中任一项所述的方法。
30.一种用于解码基因组数据的基因组解码器,其特征在于,被配置为进行根据权利要求22至28中任一项所述的方法。
31.一种计算机可读介质,其特征在于,包括当被执行时导致至少一个处理器进行根据权利要求1至21中任一项所述的方法的指令。
32.一种计算机可读介质,其特征在于,包括当被执行时导致至少一个处理器进行根据权利要求22至28中任一项所述的方法的指令。
CN201880012026.4A 2016-10-11 2018-02-14 使用基因组描述符紧凑表示生物信息学数据的方法和设备 Active CN110663022B (zh)

Applications Claiming Priority (8)

Application Number Priority Date Filing Date Title
PCT/EP2016/074301 WO2018068828A1 (en) 2016-10-11 2016-10-11 Method and system for storing and accessing bioinformatics data
PCT/EP2016/074297 WO2018068827A1 (en) 2016-10-11 2016-10-11 Efficient data structures for bioinformatics information representation
PCT/EP2016/074311 WO2018068830A1 (en) 2016-10-11 2016-10-11 Method and system for the transmission of bioinformatics data
USPCT/US2017/017842 2017-02-14
PCT/US2017/017842 WO2018071055A1 (en) 2016-10-11 2017-02-14 Method and apparatus for the compact representation of bioinformatics data
PCT/US2017/041591 WO2018071080A2 (en) 2016-10-11 2017-07-11 Method and systems for the representation and processing of bioinformatics data using reference sequences
USPCT/US2017/041591 2017-07-11
PCT/US2018/018092 WO2018152143A1 (en) 2017-02-14 2018-02-14 Method and apparatus for the compact representation of bioinformatics data using multiple genomic descriptors

Publications (2)

Publication Number Publication Date
CN110663022A CN110663022A (zh) 2020-01-07
CN110663022B true CN110663022B (zh) 2024-03-15

Family

ID=90195449

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201880012026.4A Active CN110663022B (zh) 2016-10-11 2018-02-14 使用基因组描述符紧凑表示生物信息学数据的方法和设备

Country Status (1)

Country Link
CN (1) CN110663022B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1443449A3 (en) * 2003-02-03 2006-02-22 Samsung Electronics Co., Ltd. Apparatus, method and computer readable medium for encoding a DNA sequence
CN103336916A (zh) * 2013-07-05 2013-10-02 中国科学院数学与系统科学研究院 一种测序序列映射方法及系统
WO2016202918A1 (en) * 2015-06-16 2016-12-22 Gottfried Wilhelm Leibniz Universität Hannover Method for compressing genomic data

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1443449A3 (en) * 2003-02-03 2006-02-22 Samsung Electronics Co., Ltd. Apparatus, method and computer readable medium for encoding a DNA sequence
CN103336916A (zh) * 2013-07-05 2013-10-02 中国科学院数学与系统科学研究院 一种测序序列映射方法及系统
WO2016202918A1 (en) * 2015-06-16 2016-12-22 Gottfried Wilhelm Leibniz Universität Hannover Method for compressing genomic data

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"CRAM format specification (version 3.0)";无名;《https://samtools.github.io/hts-specs/CRAMv3.pdf》;20160908;第1-39页 *
"基于码书索引变换的高通量DNA序列数据压缩算法";谭丽 等;《电子学报》;20150531;第1-7页 *

Also Published As

Publication number Publication date
CN110663022A (zh) 2020-01-07

Similar Documents

Publication Publication Date Title
CN110678929B (zh) 用于高效压缩基因组序列读段的方法和系统
CA3039688C (en) Efficient data structures for bioinformatics information representation
CN110168652B (zh) 用于存储和访问生物信息学数据的方法和系统
KR20190113969A (ko) 게놈 서열 리드의 효율적인 압축 방법 및 시스템
JP7362481B2 (ja) ゲノムシーケンスデータをコード化する方法、コード化されたゲノムデータをデコード化する方法、ゲノムシーケンスデータをコード化するためのゲノムエンコーダ、ゲノムデータをデコードするためのゲノムデコーダ、及びコンピュータ読み取り可能な記録媒体
AU2018221458B2 (en) Method and apparatus for the compact representation of bioinformatics data using multiple genomic descriptors
CN110178183B (zh) 用于传输生物信息学数据的方法和系统
CN110663022B (zh) 使用基因组描述符紧凑表示生物信息学数据的方法和设备
JP2020510907A (ja) ゲノムシーケンスリードの効率的圧縮のための方法及びシステム
NZ757185B2 (en) Method and apparatus for the compact representation of bioinformatics data using multiple genomic descriptors
EA043338B1 (ru) Способ и устройство для компактного представления биоинформационных данных с помощью нескольких геномных дескрипторов
NZ753247B2 (en) Efficient data structures for bioinformatics information representation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant