数据库驱动的原始测序数据的初步分析的制作方法

数据库驱动的原始测序数据的初步分析的制作方法
【专利说明】数据库驱动的原始测序数据的初步分析 发明领域
[0001] 本发明涉及一种用于识别生物序列的可能来源的方法。在进一步的方面，本发明 涉及一种适合用于此目的的数据库。
[0002] 发明背景
[0003] DNA测序是识别碱基（A、T、C或G)顺序的实验过程。截至今天，没有任何技术能够 测序超过几千碱基的DNA的完整分子，大多数技术测序100和200个碱基之间。细菌基因 组可容易地含有以几百万计的碱基。在过去几年中，测序成本显著降低从而使以诸如人类 健康、食品质量控制或微生物群落研宄为目的对样品中的DNA大规模测序越来越普遍。可 以想象的是，全人类基因组测序将更频繁地用于治疗以尽可能地使治疗个性化，并且将执 行常规测序来控制特定活生物体的存在或不存在。无论是作为最终目的本身或作为更复杂 的数据分析或在采取更昂贵的分析之前对测序数据的质量控制步骤的基石，快速识别可能 的起源DNA都正在迅速成为一种必然。
[0004] 初级分析包括使从测序获得的相对短的序列（称为短读段）有意义，该相对短的 序列或者是通过将它们与参考基因组比对（这需要参考物种的序列是已知的），或者是通 过不使用模型来重构拼图（所谓的测序标签的从头组装-识别未知样品的含量将需要补充 步骤）来获得的。与参考比对被认为是比从头组装在计算上更容易的任务。
[0005] 在可负担非特异性或全基因组测序之前，首先精心地对特定区域测序和组装，识 别感兴趣的预测区域。最简单的方法是通过找到由RNA翻译成蛋白质的起始密码子（ATG/ AUG)和一个终止翻译的终止密码子（TAG/UAG、TAA/UAA、TGA/UGA)所限定的间隔，寻找开放 阅读框架（ORF)。该ORF随后与所有已知基因名单比对。比对方法包括比对算法和程序如 Smith和Waterman算法、BLAST算法和程序、SSAHA和BLAT。它们的目的是要在索引序列的 数据库中找到优化比对，并通过对所有比对的得分排名找到最佳匹配以及从而找到查询序 列的最可能的功能。为了功能注释的目的，通过建立"最佳-匹配基因的组"，或直系同源基 因（COG)簇，越来越多的具有不同生物学功能的类似匹配导致这一原理的扩张。随着慢慢 地可以利用更多的完整基因组，设计了 Mummer算法来比对完整基因组对和可视化如何在 遗传相关的物种之间比较整体的基因组结构。
[0006] 由于目前在数据库中可用的序列数目，针对巨大的已知序列池的新序列的比对可 能花费相对长的时间，BLAST在这个意义上是个突破，它在找到几乎最佳结果的同时加速了 以前的算法。然而，在基于网络的搜索引擎可以几乎立即返回搜索结果的时代，针对所有已 知序列的搜索仍相对缓慢。
[0007] Ning 等，2001，（Genome: 11:1725-1729)，描述了一种算法 SSAHA(通过哈希算法的 序列搜索和比对），对含有几千兆碱基的DNA的数据库执行快速比对。SSAHA是一个比对器； 因此，其任务是向每个全长查询序列报告它们在何处以及如何匹配参考序列集合中的每个 条目。该SSAHA方法是在全长查询序列上寻找尽可能多的匹配。在数据库中的序列通过分 解为k个连续碱基的连续k-字被预处理，然后使用哈希表来存储每个k-字每次出现的位 置。在数据库中搜索查询序列是通过从哈希表取得对查询序列中每个k-字的"命中"，然后 对结果进行分选而完成。该SSAHA算法用于高通量单核苷酸多态性检测和超大规模序列组 装。在SSAHA中，每个k-字的存在和位置被存储在同一查找结构中，该结构加载到计算机 系统的存储设备中。
[0008] 已知的映射或比对算法和程序包括诸如Erland、Corona、BFAST、Bowtie、BWA、 NovoAlign的方法。它们的目标是在已知的参考中找到读段位置。推而广之，可以将无法找 到匹配的读段标记为并非来自该序列。这些程序和算法也存在搜索时间长的缺点，因为它 们都评估查询集中的每个序列（即每一测序读段），并且因为它们试图找到对于它们所有 的最佳比对（在用短读段进行时往往被称作比对）。有趣的是，因为上面的程序都使用启发 式来以精确性换取速度，因此它们所找到的结果并不相同。
[0009] US 2006286566公开了使用k-mer检测突变的方法。该方法涉及通过比较靶核酸 序列的一部分与第二序列节段，检测与靶核酸序列部分的匹配来检测靶核酸序列中的明显 突变。
[0010] US2012000411中公开了能够表征样品内的有机体群体的系统和方法，这是基于对 短序列信息串的匹配以从参考基因组数据库中识别基因组。该专利申请没有公开这样的方 法，即其中在一个参考序列中的短串集合中搜索短串的存在和在参考序列中的位置的另一 集合中搜索位置。
[0011] 发明概述
[0012] 本发明提供了一种用于识别原始序列来源，例如从测序仪获得的DNA读段（或短 读段）或者从N-或C-端测序或从质谱获得的蛋白质序列的新方法。该方法依赖于预先索 引的参考序列的集合和对传入的生物序列（如来自测序仪的读段）的查询集评分的系统， 和依赖于提交部分查询集的系统。这可以通过使用基于客户端-服务器的方法，以服务器 实体容纳参考的集合并在客户端提交查询序列的子集同时进行评分来完成。
[0013] 由本发明提供的方法，允许快速确定样品中发现的不同DNA来源，并且不依赖于 关于来源序列的给定基因的完整序列以及参考序列的知识。
[0014] 短读段，尽管并不代表其起源的完整参考，但带有该参考的标签性信号。短读段可 以进一步被分解成子序列（称为k-mer或k-字（k-tuple))并且那些k-mer搜索索引k-mer 的集合，以确定原始测序数据的来源。
[0015] 在第一方面，本发明涉及识别生物序列的可能来源的方法，该方法包括：
[0016] a)从一来源采样序列或短读段的子集，
[0017] b)将来自该子集的序列分段成k-mer，
[0018] c)针对包含参考序列的k-mer的数据库，查询来自所述子集的k-mer，
[0019] d)确定哪个（哪些）参考包含该k-mer，和
[0020] e)返回对可能的来源参考的描述。
[0021] 该方法比传统的比对和映射算法具有若干优点，传统的比对和映射算法关注于比 对全部的查询集并因此需要将来自输入装置（如客户端）的全序列传输到可以执行比对的 数据库和评分单元（如服务器）。根据本发明，只有序列的子集进行分段并且查询从而使数 据传输的需要最小化。传输的子集可以是例如，但不限于，固定大小的随机子集、过滤子集、 自适应采样、输入和评分实体之间的迭代同步或异步对话，或者是它们的任意组合。
[0022] 相比于基于测序读段组装，或基因组构建，随后搜索的方法或者相比于在参照集 合中映射所有读段的方法，本方法通过不试图执行全长比对以及通过在数据的子集上工作 而对电脑处理能力的需要大大减少，并且因此可以在几秒钟内获得结果。因此，本发明的方 法可使用客户端-服务器的方式运行，例如以具有较小的计算机处理能力（例如移动电话） 的平板或手持装置作为客户端。因为可以相对快地获得对一个数据子集的结果，用于搜索 另外的数据子集所需的时间大大减少。这样，相比于基于整个序列比对的常规方法，可以显 著减少的时间段来确定样品中的不同DNA来源的识别。
[0023] 在其最广泛的方面，本发明涉及只查询在数据库中的存在。然而，在优选的实施方 式中，也查询数据库中k-mer在参考序列中的位置，从而允许计算源k-mer的连续性并使评 估更加精确。生物体常常彼此遗传相关，本发明也能在参考序列的集合中找到近缘亲本。
[0024] 在两个单独的数据库或集合中编译数据允许将在参考中搜索k-mer的存在与搜 索位置去关联（decoupling)，并考虑优化，例如将尽可能多的对存在的搜索缓存到存储器 中，在存储器中的搜索可能比在持久性存储设备中更快。如果发现存在k-mer，那么可进行 位置搜索，并且如果在给定参考中存在足够的时间，那么在辅助优化步骤中进行，。因此本 发明的一个优选实施方式涉及一种识别生物序列的可能来源的方法，该方法包括：
[0025] a)从一来源米样序列的子集，
[0026] b)将来自该子集的序列分段成k-mer，
[0027] c)针对包括参考序列的k-mer的第一集合，查询来自所述子集的k-mer，
[0028] d)针对包括k-mer在参考序列中的位置的第二集合，查询来自所述子集的k-mer，
[0029] e)确定哪个（哪些）参考包含该k-mer，和
[0030] f)返回对可能的来源参考的描述，
[0031] 其中包括参考序列的k-mer的集合与包括k-mer在参考序列中的位置的集合是分 开的。
[0032] 因此本发明的一个优选实施方式涉及一种识别生物序列的可能来源的方法，该方 法包括：
[0033] a)从一来源采样序列或短读段的子集，
[0034] b)将来自该子集的序列分段成k-mer，
[0035] c)针对包括参考序列的k-mer的第一集合，查询来自所述子集的k-mer，
[0036] d)针对包括k-mer在参考序列中的位置的第二集合，查询来自所述子集的k-mer，
[0037] e)确定哪个（哪些）参考包含该k-mer，和
[0038] f)返回对可能的来源参考的描述，
[0039] 其中包括参考序列的k-mer的集合与包括k-mer在参考序列中的位置的集合是分 离的。
[0040] 本发明的一个显著的特点是，一旦识别出可能的参考，那么将关于该可能的参考 的信息返回给用户。返回的信息可以例如是关于可能的物种以及其起源或来源和/或该可 能的物种的全长基因组序列的信息。这允许用户使用现有技术的比对或基因组构建算法将 来自未知样品的其余的原始读段与参考序列进行比对，以便识别小的变化，如突变和插入。
[0041] 在进一步的方面，本发明涉及一种包括有参考序列的k-mer的数据库，所述数据 库包括：
[0042] a)来自参考序列的k-mer的第一集合，以及
[0043] b)每一 k-mer在该参考序列中的位置的第二集合。
[0044] 编译两个单独的数据库或集合中的数据允许将搜索k-mer在参考中的存在与搜 索位置去关联，并考虑优化，例如将对存在的搜索尽可能多地缓存到存储器中，在存储器中 搜索可能比在持久性存储设备中更快。如果发现存在k-mer，那么可进行位置搜索，并且如 果在给定参考中存在足够的时间，那么在辅助优化步骤中进行。
[0045] 在第三个方面，本发明涉及一种用于识别来源序列的可能来源的数据处理系统， 该系统优选包括：输入装置，中央处理单元，存储器以及输出装置，其中所述数据处理系统 中存储有表示当执行时使得根据本发明的方法得以被执行的系列指令的数据，该存储器还 包括根据本发明的数据库。
[0046] 图3示出本发明的系统的一个实施方式的要点。要点是采样是在"客户端"执行， 导致最少量的信息被传输。在图中并未示出对最可能的参考的描述符的应用。
[0047] 装置（输入、输出、存储器、CPU)可以是手持式、台式、基于云和/或联机的。
[0048] 优选地，数据库存储在服务器上，并且输入和输出装置是一个或多个客户端，客户 端和服务器经由数据通信连接相连并且服务器的共享允许将参考的集合集中并且如果在 单独的处理或者甚至是单独的机器上运行时，在整个客户端的服务器中分配计算能力。在 这样的实施方式中，客户端可以包括使客户端采样来源序列的子集，将这些分段成k-mer， 并将这些传输到服务器的一系列指令。
[0049] 客户端可以进一步包括一系列指令，所述指令允许客户端与服务器对话以适应或 中断采样程序或基于从服务器传送到客户端的序列，执行将来源序列组装成一个或多个更 大的序列。
[0050] 在一个实施中，系统经由数据连接与测序设备相连接。
[0051] 在进一步的方面，本发明涉及一种计算机软件产品，含有一系列在执行时使得本 发明的方法得以被执行的指令，以及涉及一种集成电路产品，含有系列在执行时使得本发 明的方法得以被执行的指令。
【附图说明】
[0052] 图1、"存在"和"位置"数据库的构建。
[0053] 图2、评分一组查询DNA片段，通常是来自测序的原始读段。
[0054] 图3、本发明的系统的构架的一般描述。
[0055] 图4、根据改变读段大小（行）和随机置换率（列），在数据库中用作查询的747 个细菌基因组的平均排名（X轴）和排名的标准偏差（y轴）。
[0056] 图5、索引和评分程序的具体实例的概述，其也被用于实施例1和2。（A)对参考序 列的集合索引过程中，非重叠的k-mer被索引为两个不同的键值对存储，一个将k-mer与发 现有它们（'存在'）的参考关联，一个将k-mer与其在发现有该k-mer的参考中的位置（' 位置'）关联。⑶当处理查询集中的测序读段时，重叠的k-mer在"存在"存储中查找。使 用重叠的k-mer允许相对迅速地解决读段开始和参考序列开始（虚线）之间的错配。在图 中，只有k-mer的子集与索引步骤同相，因此仅有它们能在"存在"中找到。（C)对于给定的 读段，将阈值仅仅施加到潜在足够匹配该读段的保留参考。使用例如在参考中的最小区域 内最高集中度的k-mer，在查询了"位置"存储的最后步骤中解决了在大量含有不相交的分 散k-mer的参考的情况，例如针对哺乳动物基因组的细菌读段。
[0057] 图6、细菌读段。对于在一组747个基因组中的每一个细菌基因组，模拟了几个读 段长度（50个核苷酸（nt)、75nt、100nt、150nt、200nt、250nt)和几个置换错误率（0%、1%、 5%、10% )。将100个随机读段用于每个查询并记录列表中的正确参考的排名分布；排名1 意味着正确的参考是在列表的最上方。返回的命中的列表被设定为最大长度25并且如果 根本不在列表中那么就将该参考算作"未找到"。正确的测试细菌基因组的百分比表示为嵌 入每个板的右侧的条。该图表明，正如所料，随着错误率增加，性能下降，同时也表明，长度 为50的读段似乎具有相对降低的性能。相比于100个核苷酸的读段，将读段长度增加至超 过100个核苷酸仅带来小的改进，并且对错误率的补偿效应有限。
[0058] 图7、细菌读段（读段数目）。对于在一组747个基因组中的每一个细菌基因组，模 拟了几个读段长度（50nt、75nt、100nt、150nt、200nt、250nt)和几个置换错误率（0%、1%、 5%、10% )。将100、200或300个随机读段用于每个查询并记录列表中的正确参考的排名 分布；排名1意味着正确的参考是在列表的最上方。曲线表不100、200和300个读段。由 此可以看出，随机样品中的读段数目从100个读段增加到300个读段使得性能相对较小的 增加。错误率或读段长度具有更大的作用。
[0059] 图8、细菌读段，对于747个测试细菌基因组，执行5次识别程序的一次迭代时，真 实参考的性能变异性，平均排名（排名，X轴）和排名的标准差（S rank, 7轴）。平均排名越 接近1，越接近完美的性能，并且排名的标准差越小，采样效果的敏感性越小。为了提高清晰 度，当测试的大量细菌基因组在散点图上产生等于或接近坐标时，使用六边形分选并相应 地着色该区域。各散点图的右侧竖条表示没有前25个匹配内的测试基因组的数目，并且与 六边形分选以相同分级着色。尝试了不同的读段大小（行）和错误率（随机置换，列），产 生散点图矩阵。
[0060] 图9、细菌读段，同一物种。给出正确物种，也就是在我们的集合中属于相同物种的 细菌的参考，而不