CN104919466A

CN104919466A - 数据库驱动的原始测序数据的初步分析

Info

Publication number: CN104919466A
Application number: CN201380065692.1A
Authority: CN
Inventors: L·戈蒂埃; O·伦德
Original assignee: Danmarks Tekniskie Universitet
Current assignee: Danmarks Tekniskie Universitet
Priority date: 2012-10-15
Filing date: 2013-10-11
Publication date: 2015-09-16
Also published as: EP2915084A1; WO2014060305A1; JP2016502162A; US20150294065A1

Abstract

本发明涉及一种识别含有原始测序读段的样品中的生物序列的来源的方法。该方法可以用于识别未知DNA的来源并且可以用于诊断、生物防御、食品安全和质量以及卫生应用。在另一方面，本发明涉及可以用于本发明方法的参考序列数据库。

Description

数据库驱动的原始测序数据的初步分析

发明领域

本发明涉及一种用于识别生物序列的可能来源的方法。在进一步的方面，本发明涉及一种适合用于此目的的数据库。

发明背景

DNA测序是识别碱基(A、T、C或G)顺序的实验过程。截至今天，没有任何技术能够测序超过几千碱基的DNA的完整分子，大多数技术测序100和200个碱基之间。细菌基因组可容易地含有以几百万计的碱基。在过去几年中，测序成本显著降低从而使以诸如人类健康、食品质量控制或微生物群落研究为目的对样品中的DNA大规模测序越来越普遍。可以想象的是，全人类基因组测序将更频繁地用于治疗以尽可能地使治疗个性化，并且将执行常规测序来控制特定活生物体的存在或不存在。无论是作为最终目的本身或作为更复杂的数据分析或在采取更昂贵的分析之前对测序数据的质量控制步骤的基石，快速识别可能的起源DNA都正在迅速成为一种必然。

初级分析包括使从测序获得的相对短的序列(称为短读段)有意义，该相对短的序列或者是通过将它们与参考基因组比对(这需要参考物种的序列是已知的)，或者是通过不使用模型来重构拼图(所谓的测序标签的从头组装-识别未知样品的含量将需要补充步骤)来获得的。与参考比对被认为是比从头组装在计算上更容易的任务。

在可负担非特异性或全基因组测序之前，首先精心地对特定区域测序和组装，识别感兴趣的预测区域。最简单的方法是通过找到由RNA翻译成蛋白质的起始密码子(ATG/AUG)和一个终止翻译的终止密码子(TAG/UAG、TAA/UAA、TGA/UGA)所限定的间隔，寻找开放阅读框架(ORF)。该ORF随后与所有已知基因名单比对。比对方法包括比对算法和程序如Smith和Waterman算法、BLAST算法和程序、SSAHA和BLAT。它们的目的是要在索引序列的数据库中找到优化比对，并通过对所有比对的得分排名找到最佳匹配以及从而找到查询序列的最可能的功能。为了功能注释的目的，通过建立“最佳-匹配基因的组”，或直系同源基因(COG)簇，越来越多的具有不同生物学功能的类似匹配导致这一原理的扩张。随着慢慢地可以利用更多的完整基因组，设计了Mummer算法来比对完整基因组对和可视化如何在遗传相关的物种之间比较整体的基因组结构。

由于目前在数据库中可用的序列数目，针对巨大的已知序列池的新序列的比对可能花费相对长的时间，BLAST在这个意义上是个突破，它在找到几乎最佳结果的同时加速了以前的算法。然而，在基于网络的搜索引擎可以几乎立即返回搜索结果的时代，针对所有已知序列的搜索仍相对缓慢。

Ning等，2001，(Genome:11:1725-1729)，描述了一种算法SSAHA(通过哈希算法的序列搜索和比对)，对含有几千兆碱基的DNA的数据库执行快速比对。SSAHA是一个比对器；因此，其任务是向每个全长查询序列报告它们在何处以及如何匹配参考序列集合中的每个条目。该SSAHA方法是在全长查询序列上寻找尽可能多的匹配。在数据库中的序列通过分解为k个连续碱基的连续k-字被预处理，然后使用哈希表来存储每个k-字每次出现的位置。在数据库中搜索查询序列是通过从哈希表取得对查询序列中每个k-字的“命中”，然后对结果进行分选而完成。该SSAHA算法用于高通量单核苷酸多态性检测和超大规模序列组装。在SSAHA中，每个k-字的存在和位置被存储在同一查找结构中，该结构加载到计算机系统的存储设备中。

已知的映射或比对算法和程序包括诸如Erland、Corona、BFAST、Bowtie、BWA、NovoAlign的方法。它们的目标是在已知的参考中找到读段位置。推而广之，可以将无法找到匹配的读段标记为并非来自该序列。这些程序和算法也存在搜索时间长的缺点，因为它们都评估查询集中的每个序列(即每一测序读段)，并且因为它们试图找到对于它们所有的最佳比对(在用短读段进行时往往被称作比对)。有趣的是，因为上面的程序都使用启发式来以精确性换取速度，因此它们所找到的结果并不相同。

US 2006286566公开了使用k-mer检测突变的方法。该方法涉及通过比较靶核酸序列的一部分与第二序列节段，检测与靶核酸序列部分的匹配来检测靶核酸序列中的明显突变。

US2012000411中公开了能够表征样品内的有机体群体的系统和方法，这是基于对短序列信息串的匹配以从参考基因组数据库中识别基因组。该专利申请没有公开这样的方法，即其中在一个参考序列中的短串集合中搜索短串的存在和在参考序列中的位置的另一集合中搜索位置。

发明概述

本发明提供了一种用于识别原始序列来源，例如从测序仪获得的DNA读段(或短读段)或者从N-或C-端测序或从质谱获得的蛋白质序列的新方法。该方法依赖于预先索引的参考序列的集合和对传入的生物序列(如来自测序仪的读段)的查询集评分的系统，和依赖于提交部分查询集的系统。这可以通过使用基于客户端-服务器的方法，以服务器实体容纳参考的集合并在客户端提交查询序列的子集同时进行评分来完成。

由本发明提供的方法，允许快速确定样品中发现的不同DNA来源，并且不依赖于关于来源序列的给定基因的完整序列以及参考序列的知识。

短读段，尽管并不代表其起源的完整参考，但带有该参考的标签性信号。短读段可以进一步被分解成子序列(称为k-mer或k-字(k-tuple))并且那些k-mer搜索索引k-mer的集合，以确定原始测序数据的来源。

在第一方面，本发明涉及识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自该子集的序列分段成k-mer，

c)针对包含参考序列的k-mer的数据库，查询来自所述子集的k-mer，

d)确定哪个(哪些)参考包含该k-mer，和

e)返回对可能的来源参考的描述。

该方法比传统的比对和映射算法具有若干优点，传统的比对和映射算法关注于比对全部的查询集并因此需要将来自输入装置(如客户端)的全序列传输到可以执行比对的数据库和评分单元(如服务器)。根据本发明，只有序列的子集进行分段并且查询从而使数据传输的需要最小化。传输的子集可以是例如，但不限于，固定大小的随机子集、过滤子集、自适应采样、输入和评分实体之间的迭代同步或异步对话，或者是它们的任意组合。

相比于基于测序读段组装，或基因组构建，随后搜索的方法或者相比于在参照集合中映射所有读段的方法，本方法通过不试图执行全长比对以及通过在数据的子集上工作而对电脑处理能力的需要大大减少，并且因此可以在几秒钟内获得结果。因此，本发明的方法可使用客户端-服务器的方式运行，例如以具有较小的计算机处理能力(例如移动电话)的平板或手持装置作为客户端。因为可以相对快地获得对一个数据子集的结果，用于搜索另外的数据子集所需的时间大大减少。这样，相比于基于整个序列比对的常规方法，可以显著减少的时间段来确定样品中的不同DNA来源的识别。

在其最广泛的方面，本发明涉及只查询在数据库中的存在。然而，在优选的实施方式中，也查询数据库中k-mer在参考序列中的位置，从而允许计算源k-mer的连续性并使评估更加精确。生物体常常彼此遗传相关，本发明也能在参考序列的集合中找到近缘亲本。

在两个单独的数据库或集合中编译数据允许将在参考中搜索k-mer的存在与搜索位置去关联(decoupling)，并考虑优化，例如将尽可能多的对存在的搜索缓存到存储器中，在存储器中的搜索可能比在持久性存储设备中更快。如果发现存在k-mer，那么可进行位置搜索，并且如果在给定参考中存在足够的时间，那么在辅助优化步骤中进行，。因此本发明的一个优选实施方式涉及一种识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列的子集，

b)将来自该子集的序列分段成k-mer，

c)针对包括参考序列的k-mer的第一集合，查询来自所述子集的k-mer，

d)针对包括k-mer在参考序列中的位置的第二集合，查询来自所述子集的k-mer，

e)确定哪个(哪些)参考包含该k-mer，和

f)返回对可能的来源参考的描述，

其中包括参考序列的k-mer的集合与包括k-mer在参考序列中的位置的集合是分开的。

因此本发明的一个优选实施方式涉及一种识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自该子集的序列分段成k-mer，

e)确定哪个(哪些)参考包含该k-mer，和

f)返回对可能的来源参考的描述，

其中包括参考序列的k-mer的集合与包括k-mer在参考序列中的位置的集合是分离的。

本发明的一个显著的特点是，一旦识别出可能的参考，那么将关于该可能的参考的信息返回给用户。返回的信息可以例如是关于可能的物种以及其起源或来源和/或该可能的物种的全长基因组序列的信息。这允许用户使用现有技术的比对或基因组构建算法将来自未知样品的其余的原始读段与参考序列进行比对，以便识别小的变化，如突变和插入。

在进一步的方面，本发明涉及一种包括有参考序列的k-mer的数据库，所述数据库包括：

a)来自参考序列的k-mer的第一集合，以及

b)每一k-mer在该参考序列中的位置的第二集合。

编译两个单独的数据库或集合中的数据允许将搜索k-mer在参考中的存在与搜索位置去关联，并考虑优化，例如将对存在的搜索尽可能多地缓存到存储器中，在存储器中搜索可能比在持久性存储设备中更快。如果发现存在k-mer，那么可进行位置搜索，并且如果在给定参考中存在足够的时间，那么在辅助优化步骤中进行。

在第三个方面，本发明涉及一种用于识别来源序列的可能来源的数据处理系统，该系统优选包括：输入装置，中央处理单元，存储器以及输出装置，其中所述数据处理系统中存储有表示当执行时使得根据本发明的方法得以被执行的系列指令的数据，该存储器还包括根据本发明的数据库。

图3示出本发明的系统的一个实施方式的要点。要点是采样是在“客户端”执行，导致最少量的信息被传输。在图中并未示出对最可能的参考的描述符的应用。

装置(输入、输出、存储器、CPU)可以是手持式、台式、基于云和/或联机的。

优选地，数据库存储在服务器上，并且输入和输出装置是一个或多个客户端，客户端和服务器经由数据通信连接相连并且服务器的共享允许将参考的集合集中并且如果在单独的处理或者甚至是单独的机器上运行时，在整个客户端的服务器中分配计算能力。在这样的实施方式中，客户端可以包括使客户端采样来源序列的子集，将这些分段成k-mer，并将这些传输到服务器的一系列指令。

客户端可以进一步包括一系列指令，所述指令允许客户端与服务器对话以适应或中断采样程序或基于从服务器传送到客户端的序列，执行将来源序列组装成一个或多个更大的序列。

在一个实施中，系统经由数据连接与测序设备相连接。

在进一步的方面，本发明涉及一种计算机软件产品，含有一系列在执行时使得本发明的方法得以被执行的指令，以及涉及一种集成电路产品，含有系列在执行时使得本发明的方法得以被执行的指令。

附图说明

图1、“存在”和“位置”数据库的构建。

图2、评分一组查询DNA片段，通常是来自测序的原始读段。

图3、本发明的系统的构架的一般描述。

图4、根据改变读段大小(行)和随机置换率(列)，在数据库中用作查询的747个细菌基因组的平均排名(x轴)和排名的标准偏差(y轴)。

图5、索引和评分程序的具体实例的概述，其也被用于实施例1和2。(A)对参考序列的集合索引过程中，非重叠的k-mer被索引为两个不同的键值对存储，一个将k-mer与发现有它们('存在')的参考关联，一个将k-mer与其在发现有该k-mer的参考中的位置('位置')关联。(B)当处理查询集中的测序读段时，重叠的k-mer在“存在”存储中查找。使用重叠的k-mer允许相对迅速地解决读段开始和参考序列开始(虚线)之间的错配。在图中，只有k-mer的子集与索引步骤同相，因此仅有它们能在“存在”中找到。(C)对于给定的读段，将阈值仅仅施加到潜在足够匹配该读段的保留参考。使用例如在参考中的最小区域内最高集中度的k-mer，在查询了“位置”存储的最后步骤中解决了在大量含有不相交的分散k-mer的参考的情况，例如针对哺乳动物基因组的细菌读段。

图6、细菌读段。对于在一组747个基因组中的每一个细菌基因组，模拟了几个读段长度(50个核苷酸(nt)、75nt、100nt、150nt、200nt、250nt)和几个置换错误率(0％、1％、5％、10％)。将100个随机读段用于每个查询并记录列表中的正确参考的排名分布；排名1意味着正确的参考是在列表的最上方。返回的命中的列表被设定为最大长度25并且如果根本不在列表中那么就将该参考算作“未找到”。正确的测试细菌基因组的百分比表示为嵌入每个板的右侧的条。该图表明，正如所料，随着错误率增加，性能下降，同时也表明，长度为50的读段似乎具有相对降低的性能。相比于100个核苷酸的读段，将读段长度增加至超过100个核苷酸仅带来小的改进，并且对错误率的补偿效应有限。

图7、细菌读段(读段数目)。对于在一组747个基因组中的每一个细菌基因组，模拟了几个读段长度(50nt、75nt、100nt、150nt、200nt、250nt)和几个置换错误率(0％、1％、5％、10％)。将100、200或300个随机读段用于每个查询并记录列表中的正确参考的排名分布；排名1意味着正确的参考是在列表的最上方。曲线表示100、200和300个读段。由此可以看出，随机样品中的读段数目从100个读段增加到300个读段使得性能相对较小的增加。错误率或读段长度具有更大的作用。

图8、细菌读段，对于747个测试细菌基因组，执行5次识别程序的一次迭代时，真实参考的性能变异性，平均排名(排名，x轴)和排名的标准差(S_rank，y轴)。平均排名越接近1，越接近完美的性能，并且排名的标准差越小，采样效果的敏感性越小。为了提高清晰度，当测试的大量细菌基因组在散点图上产生等于或接近坐标时，使用六边形分选并相应地着色该区域。各散点图的右侧竖条表示没有前25个匹配内的测试基因组的数目，并且与六边形分选以相同分级着色。尝试了不同的读段大小(行)和错误率(随机置换，列)，产生散点图矩阵。

图9、细菌读段，同一物种。给出正确物种，也就是在我们的集合中属于相同物种的细菌的参考，而不是图7中示出的完全相同的参考的匹配百分比，以及正确物种不在顶部25个匹配的情形下的百分比。对于短读段(50nt)性能相对较差，噪音进一步降低它(第一行的条线图)，但从100nt变得非常好，并保持对噪声的稳健性。

本发明的详细描述

本发明平衡了对来自样品中发现的蛋白质、DNA或RNA的生物序列信息的可能来源执行识别的速度和精度。

在本发明方法中要使用的序列信息可以例如是来自核酸测序仪或者来自蛋白质的C-或N-端测序或来自质谱蛋白质测序的原始读段。因此，在本发明的语境中的措词采样序列是指这种也被称为短读段的原始读段。

在图2中描述的本发明一个具体实施方式可以涉及：

·用参考DNA创建数据库(参见图1)。该数据库是两部分：1)相对于参考来索引的所有参考DNA的k-mer的数据库和2)来自数据库1的k-mer与在该参考序列中的位置之间的关联性的数据库。因而参考k-mer ID和位置被存储在两个不同的数据库中。

图1示出了数据库构建的一个实施方式。用来创建数据库的输入是来自公共或专有数据库的DNA。这些然后都被分成k-mer，其可以优选是不重叠的以节省空间。k-mer还可以是2比特位类型的(2-bit bit packed)，这意味着每个碱基只占用存储器的2个比特。为了加快存储k-mer，优选在插入数据库之前对这些k-mer进行分选。此外在该k-mer所来源的参考序列中的名称和位置可以存储在单独的数据库中。

·针对参考数据库，搜索来自一来源的查询序列的分解成k-mer的读段的选择。

·由来自该查询序列的k-mer数计算主要得分，该k-mer可以在数据库中的给定参考序列中找到。

·将建议的序列返回用户，并且可以用于更大量且传统的分析。

本发明的这一实施的特征是：

·在搜索期间，只有完全匹配的k-mer被登记。

·查询读段分解成多个k-mer，如长度为16。每个k-mer的起点增量为1。

·不是“传统”的从头、比对或映射方法。

图2示出了用于搜索k-mer数据库的一个可能的算法。使用具有一个步长的滑动窗将读段分成k-mer。如果在当前搜索中已经遇到(访问)k-mer，那么选择下一个k-mer。随后在k-mer数据库中查找该k-mer。如果它在数据库中，那么提取在参考序列中的身份和位置。然后计算读段的近似连续性，并且如果最大连续节段超过阈值，那么命中计数增加。这对于读段中的所有k-mer进行重复。对于每个读段，将得分计算为命中数(命中计数)除以查询序列的长度，并且计算命中计数除以匹配参考序列的长度。这对多个读段进行重复，这可以依赖于所获得的得分被先验或动态定义。得分被分选并且将最佳匹配返回给用户。

完全匹配不是在读段水平进行。评分允许沿读段缺失k-mer匹配(因此确保了对生物样品中的测序错误和突变的稳健性)。

该系统的概述是：

·索引所有已知的参考DNA序列为k-mer，存储参考(例如物种)和在参考序列中的位置。这一步优选只在参考DNA序列通过增加新序列或通过添加更多的序列信息来更新时进行。

·客户端，可以通过将DNA的短序列分成k-mer，针对数据库匹配它们并计数对参考序列的命中数，优选以位置信息精化该匹配，来存储它们。

获得的参考随后可用于：

·过滤掉匹配参考的读段，并且找到是否存在来自另一不同参考但低丰度的DNA

·针对该参考执行比对，或利用数据库中的参考反复构建更大的片段，这导致比通过利用预先组装的参考进行从头组装更好的性能；此外，该性能将随着数据库大小增大而提高并且加入了更多的组装好的参考

·识别各种生物体或基因(例如与用于诊断目的相关)的可能存在。

因为只有原始读段的子样品是必要的，这可以减少为执行基本的诊断(如识别传染剂)而被传输的数据量。在更小的序列实验的情况下，这也允许一些分析通过客户端在商品化硬件上进行。

随着低通量桌面测序仪(或一次性测序单元)的发展以及更便宜的GPU或FPGA的兴起，联合了允许对测序数据进行实时或接近与实时的初步分析的技术。

算法

在一个方面，本发明涉及识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自该子集的序列分段成k-mer，

b)针对包含参考序列的k-mer的数据库，查询来自所述子集的k-mer，

d)确定哪个(哪些)参考包含该k-mer，和

e)返回对可能的来源参考的描述。

术语“从一来源的序列”被用于指定从包含生物序列的样品中获得的序列。样品可以是环境样品、来自受试者(如患者)的样品、来自犯罪现场的样品、食品样品、水样品等。样品进行现有技术的DNA/RNA或蛋白质分离和测序方法。结果是一组表征该样品的序列(也称为读段)。该序列通常是在一定间隔内的随机长度。该序列也通常是随机重叠的。来自样品(称为来源序列)的每个序列可进行本发明的方法。

根据本发明的术语“参考”包括存储在数据库中的序列的描述符。参考的典型例子是特定物种或栽培物，或分离株的全长基因组序列。参考也可以由特定物种或特定条件的物种的转录组或蛋白质组组成。物种的转录组和蛋白质组可以响应于年龄和环境条件随时间而改变，而例如物种的基因组序列仍随时间或多或少恒定。数据库可存储关于参考的额外信息。

本发明的方法可以适用于任何生物序列信息，例如氨基酸序列和核苷酸序列如DNA和RNA序列。在优选的实施方式中，序列是DNA序列。

在最广泛的方面，本发明仅依赖于从查询或来源序列中识别k-mer的存在。在这种情况下，从算法输出的是参考列表和在参考中识别的相应命中数目。然而，由于一些基因组，例如人类基因组并且特别是一些植物基因组的大小，许多k-mer可能偶然存在于这些非常大的基因组中。因此，在优选的实施方式中，查询还包括确定k-mer在参考序列中的位置。这允许使用存在和位置来确定查询k-mer在参考序列中的连续性。这使得查询更精确，因为可使用同时基于k-mer在参考中的存在和位置，或近似连续性的得分。因此本发明的一个优选实施方式涉及一种识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自该子集的序列分段成k-mer，

c)针对包括参考序列的k-mer的第一集合，查询来自所述子集的一个或多个k-mer，

d)针对包括k-mer在参考序列中的位置的第二集合，查询来自所述子集的一个或多个k-mer，

e)确定哪个(哪些)参考含有该k-mer，和

f)返回对可能的来源参考的描述，

在本发明的一个甚至更优选的实施方式中，针对包括k-mer在参考序列中的位置的第二集合的查询仅仅是当在包括参考序列的k-mer的第一集合中已经发现(即存在)给定k-mer时才进行(参见图2)。

在本发明的一个优选实施方式中，当使用上述步骤a)到f)时，在查询随后的k-mer之前，确定给定k-mer的存在和位置。因此本发明的一个优选实施方式涉及一种识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自该子集的序列分段成k-mer，

d)针对包括k-mer在参考序列中的位置的第二集合，查询来自所述子集的所述k-mer，

e)确定哪个(哪些)参考含有该k-mer，和

f)返回对可能的来源参考的描述，

本发明的一个显著的特点是，只将从测序获得的序列的子集用于查询数据库。这使使得在测序和查询非常大的基因组时可能是限速步骤的数据传输最小化。因此序列的子集可以包括离散序列中的至少1％，例如至少2％，如至少4％，例如至少5％，如至少6％，例如至少7.5％，例如至少10％，如至少15％，例如至少25％，如至少30％，例如至少35％，如至少40％，例如至少50％。

本发明的一个特征在于k-mer查询涉及确定查询k-mer和参考k-mer之间的精确匹配。

当查询来源序列或短读段时，查询优选涉及查询来自至少一个来源序列的所有k-mer。这允许对连续性或近似连续性的最佳计算。优选地，查询来自至少50个来源序列的所有k-mer，例如来自至少100，如至少150个，例如至少200个，如至少250个，例如至少300个，如至少400个，例如至少500个，如至少750，例如至少1000个，例如至少1500个，如至少2000个，例如至少2500个，如至少5000或更多个序列。所查询的来源序列的确切数目特别是通过网络和计算能力、时间限制、统计要求和全长来源序列的大小和与不同参考的源的相关性来确定。

如在实施例中所示，每个来源序列优选是给定的最小长度以得到源生物体、品种、栽培物或分离株的特性指纹。在来源序列是核苷酸序列的情况下，来源序列优选至少50个核苷酸碱基，更优选至少75个核苷酸碱基，例如75至200个核苷酸碱基，如75个核苷酸碱基至100个核苷酸碱基，或100个核苷酸碱基至125个核苷酸碱基，或125个核苷酸碱基至150个核苷酸碱基，或150个核苷酸碱基至175个核苷酸碱基，或175个核苷酸碱基至200个核苷酸碱基，甚至更优选至少100个核苷酸碱基，例如100-300个核苷酸碱基如，如100个核苷酸碱基至150个核苷酸碱基，或150个核苷酸碱基至200个核苷酸碱基，或200个核苷酸碱基至250个核苷酸碱基，或250个核苷酸碱基至300个核苷酸碱基，如至少100个核苷酸碱基，例如100个核苷酸碱，例如200个核苷酸碱基，如至少250个核苷酸碱基，例如300个核苷酸碱基，如400个核苷酸碱基，至少500或更多个核苷酸碱基。

在许多实际实施中，初始查询序列的一个子集。如果这还不足以以足够高的确定性来确定参考，该方法可以进一步包括选择序列的一个或多个另外的子集并将它们进行本发明方法的步骤a)至f)。

原则上，该方法允许使用任何大小的k-mer或k-字。然而，在一个优选实施方式中，k-mer的大小可以被4整除。因此，k-mer可以是大小为4、8、12、16、20、24、28、32、36、40、44、48、52、56、60、64或更长。更优选地，k-mer的长度为16和64之间，更优选16和32之间。更长的k-mer使该方法对测序错误更敏感而更短的k-mer增加随机命中的数量，从而提供噪声。

在一个实施方式中，k-mer是连续的，并且优选地，存储在数据库中的k-mer是连续的，以便覆盖整个参考序列。

优选地，来自来源序列的k-mer是重叠的并且增量为至少一个碱基或氨基酸，例如至少两个，如至少3个，例如至少4个，如至少5个，例如至少6个或更多。这相当于横跨序列滑动宽度为k的窗口。横跨序列可以滑动一、二或更多个碱基/氨基酸的窗口。通过从来源序列进行重叠和增量k-mer，该方法对测序错误或点突变变得不敏感，因为将在查询中识别在例如单个碱基突变/错误的任一侧的k-mer。因此，可以以更高的精度计算连续性。

使用从来源序列的不相交子序列的串联得到的不相交的k-mer也是可能的。

优选地，根据该方法，针对数据库查询来自给定序列的k-mer以确定在一个或多个参考序列中k-mer的存在以及该k-mer在所述一个或多个参考序列中的位置。为了优化数据库使用，优选仅当k-mer存在于数据库中时查询位置。

为了允许定量评价所述查询，该方法涉及为所识别的参考序列计算得分，所述得分与在给定的参考序列中找到的来自一个或多个序列的k-mer数目相关联。这个得分可以例如被来源序列的长度整除。可为识别的参考计算进一步的得分，所述进一步的得分与在参考序列中找到的来自一个或多个序列的k-mer的连续性相关联。例如，该得分可以是来自在数据库中找到的一个来源序列的k-mer和在该数据库的一个参考序列中找到的k-mer的最长序列的百分比。

类似地，对于每个识别的参考序列，可以为识别的参考计算得分，所述得分与在参考序列中的k-mer数目相关联，所述k-mer也存在于来自所述来源的k-mer的子集中。一个实例可以是数据库中来自来源序列中发现的一个参考的k-mer的百分比。在许多实际应用中，查询数百个来源序列并评分，以获得满意的确定性。这个得分还可以包括基于所识别的k-mer的连续性的得分。

优选地，为每个不同的来源序列计算这些得分，例如其中查询来自一个来源序列的所有k-mer并且为所述来源序列计算一或多个得分。优选地，该方法还涉及查询来自第二来源序列，优选来自第三来源序列的所有k-mer等。对于不同的来源序列的得分可例如通过将它们以来源序列的长度加权而合并。

在本发明的一个实施方式中，一旦对为读段所生成的所有k-mer进行处理，在参考中匹配的连续位置的数目被用来隔离最大的匹配簇，即，在所有匹配参考中源自相同读段的最大集中度的匹配k-mer。对于每个这样的簇，计数是通过将在一个簇中的k-mer数量加上给定的参考序列的计数来计算的。当该方法是对给定样品的一个以上的读数迭代时，该计数可通过将在一个簇中的k-mer的数量加上从先前的读段获得的参考序列的计数来更新。即，计数可以通过加上用于该参考的k-mer数量来更新并且已经计数的k-mer列表被更新。然后可以处理下一个序列或读段。获得与发现匹配的k-mer计数相关联的参考列表。对于每对<参考，计数>，计数除以查询集中独一无二的k-mer的数量，提供了在给定参考匹配的被查询子集中DNA量的粗略得分。如果被查询的子集是完全匹配该序列，那么得分将是1，否则会降低；例如，如果所查询的子集是两个参考的等比例混合物，那么对两个参考的得分都将是0.5左右。该计数也可以除以参考的大小(或在参考序列中独一无二的k-mer的数量)，得到由所查询的子集表示的参考部分的粗略得分；该第二得分有助于对匹配参考进行分选，并避免偏向最大参考。最终的得分是这两个得分的加权和，例如其中对每个得分使用相等的加权。

在本发明的一个实施方式中，查询预先选定数目的来源序列并返回结果。然而，在其他实施方式中，一旦以预定义的统计概率识别出参考生物体，那么就可以停止数据库查询。同样，如果在数据库中没有找到预定义分数的k-mer或扩展更多的来源序列，或以松弛参数计算的得分，那么可以停止数据库查询。这在垃圾序列，具有许多测序错误的序列或完全未知序列的情况时可能如此。

来自查询过程的输出可以是根据一个或多个所述得分排名的可能的来源参考列表。数据库输出的其它实例包括关于一个或多个可能参考的下列信息中的一个或多个：可能参考的分类学名称、所述可能参考的近亲、所述参考的来源、遗传连锁信息、关于SNP、基因在序列中的位置和注释的信息。

在一个具体的实施方式中，数据库输出最可能的参考的序列，优选地，其中所述数据库输出最可能的参考物种的全基因组序列。这允许用户使用现有技术的比对算法，针对最可能的物种的全基因组序列比对来源序列，以进一步调查是否有突变或插入或染色体异常(anomaly)，异常(abnormality)或畸变。然而，在本发明的一个实施方式中，本发明的方法不涉及对序列数据使用比对算法，如，例如使用得分矩阵的算法，如，例如Smith-Watermanalgorithm[14],BLAST[1],BLAT[5],Bowtie,BWA,SHRiMP[16]，或熟练技术人员已知的其他比对算法。

在许多情况下，例如，当查询微生物的序列时，该数据库可包括许多密切相关的序列，如来自相同物种的不同分离株的序列。在这样的情况下，来自具有非常相似的序列的参考的结果可在输出中被分组。这还可以允许用户更容易地识别来自以较低量存在的另一物种或不同物种中的小片插入的DNA。

在许多情况下，样品包含物种的混合种群并且全基因组测序这将导致来自几个物种的基因组DNA的混合物。在这种情况下，该方法可以涉及执行该方法的几次迭代，例如在第一次迭代中识别最丰富的参考。在第二迭代中，从来源序列中除去来自最丰富的物种的序列，然后查询数据库或者该方法可以涉及忽略来自该物种的进一步的结果。

可替代地，来自本发明的方法的一个迭代的输出可以包括用于所识别的所有参考的信息和得分。在这种情况下的得分可以包括在不同参考之间的百分比分布。

本实施方式也可以用于识别插入物的参考，例如病毒插入物、转基因或来自另一细菌物种的插入物。

在许多实施方式中，用户将最初知晓在样品中存在来自一个参考的序列或短读段，并且随后的任务是识别在该样品中存在的任何其它序列或短读段的可能的参考。这在诊断学的情况下会是如此，其中样品包含人类DNA和来自可能的病原体的DNA两者。其它实例包括对食物样品中有害细菌的识别，其中，已知样品含有来自食物源(如沙拉、番茄、黄瓜、来自特定物种的肉)的DNA并且任务是识别任何污染的DNA的存在和身份。在这样的方法中，该方法可以包括最初去除与来自预先定义的参考的序列对准的来源序列。或者，该方法可以涉及忽略来自一个或多个预先定义的参考的k-mer。

在一个实施方式中，该方法涉及在从核酸测序仪获得原始读段时采样并查询它们。

当为了诊断目的识别DNA数据的查询集(如来自测序仪的短读段或原始读段)时，我们认为包括映射或比对针对于全面的参考数据库的所有读段的蛮力方法有两个主要的缺点：首先是几百兆或千兆字节的数据的大部分或者从测序设施转移到计算中心，其次是执行该任务必须的计算资源是至关重要的。假设一个参考集合包含10,000大肠杆菌规模的细菌并且需要30秒优化比对器，如BWA和bowtie2来处理250Mb的原始测序数据(如果基因组是4Mb大小，平均覆盖为约60x)，尽管可以在多个CPU上并行进行处理，但是这样的CPU处理会花费3天半。可以进行精化使得基因组串联但代价是需要日益增加的存储器容量，向初始参考基因组分配映射位置的后处理运算，以及参考接近的基因组时不可避免的多个匹配，这是短读段比对器往往不适的方面。使用FM-索引在大小为u的参考中定位长度p的串的n个出现次数具有上限O(p+n logεu)，这意味着，尽管随着参考大小(以logε定义的术语)的增加，复杂性增长缓慢，但它与高度相似的基因组的数量呈线性增长。我们的方法包含了巨大的参考数据库的观点，并且不试图在一台计算机的所有RAM中容纳它们。

数据库

在一个方面，本发明涉及包含参考序列的k-mer的数据库，所述数据库包括：

a.来自参考序列的k-mer的第一集合，和

b.每个k-mer在该参考序列中的位置的第二个集合。

数据库架构允许非常快速地查询来自来源序列的k-mer，如示于所附实施例中，其证明了结果可在几秒钟内被返回。

该数据库可以进一步包括关于与给定参考相关联的全长序列，和/或所述参考的来源，和/或所述参考的一种或多种分类学描述符的信息。可以被存储的附加信息是关于DNA序列中的基因注释的信息。

当建立该数据库时，k-mer可以进行哈希函数，向每个独一无二的k-mer分配唯一的密钥。其他的可能性包括搜索树或哈希函数和搜索树的组合。该唯一的密钥可以与关于k-mer所在的那些参考的信息相关联。

在第二集合中，第二集合中的每个独一无二的k-mer也可以用作密钥，并通过哈希表、搜索树或它们的组合与关于在每个参考中k-mer(k-mer在该参考中存在)的位置的信息相关联。此集合可以包括关于k-mer存在的位置的进一步的信息，例如与序列如编码序列、调节序列等的任何注解的关联性。

关于存在给定k-mer的参考序列的一个或多个另外的信息，例如与序列的任何注释的关联性、编码序列、调控序列、可能参考的分类学名称、所述可能参考的近亲、所述参考的来源、一组另外相关的参考、参考是从哪里获得的(例如土壤、海洋、内脏、下水道等)、参考序列是何时获得的、分类学分类、近缘物种、关于参考序列是从哪个数据库(例如NCBI、EBI/Sanger)下载的信息，或者其他信息都可存储在单独的数据库中，例如SQL数据库，其可以另外用于提取关于根据本发明的参考序列的信息。

术语“一组另外相关的序列”是指来自取自诸如土壤、海洋、内脏、下水道等的相似环境的样品的序列。

因此，在本发明的一个实施方式中，包括参考序列的k-mer的数据库包括：

a)来自参考序列的k-mer的第一集合，和

b)每一k-mer在该参考序列中的位置的第二集合。

c)第三集合或数据库，其具有选自下组的参考标识符和一或多个信息：描述行、数据来源、可能参考的分类学名称、所述可能参考的近亲、所述参考的来源、一组另外相关的参考的信息、参考是从哪里获得的(例如土壤、海洋、内脏、下水道等)、参考序列是何时获得的、分类学分类、近缘物种、关于参考序列是从哪个数据库(例如NCBI、EBI/Sanger或其他数据库)下载的信息。

在一个优选的实施方式中，k-mer的第一集合是关联到每个k-mer(数据库中的密钥)的键值对存储或NoSQL数据库，例如KyotoCabinet，对应于具有该k-mer的参考的标识符列表如图1中所示。k-mer在参考序列中的位置的第二集合也可以被存储在键值对存储或NoSQL数据库，例如KyotoCabinet中(参见图1)。参考标识符和信息条(例如描述行和数据来源)之间的关联性被存储在单独的SQL数据库中。

虽然假定有足够的查找表，但k-mer在数据库中的长度优选匹配k-mer在来源序列中的长度。然而，k-mer在数据库中优选不重叠。使用重叠的k-mer将增加数据处理时间。

根据本发明，在数据库中参考序列的索引的k-mer可以是重叠或非重叠的。在优选实施方式中，索引的参考序列的k-mer是非重叠的。本领域技术人员可以理解，相似的评分原则可以用于参考序列中的非重叠或重叠的k-mer的索引数据库。

如果将树或哈希用于k个索引和查找，在以k-mer索引的大小为u的参考中定位长度p的串的n个出现的时间复杂度具有O(p+n log u)或O(p+n)的复杂度。

这并不排除这样的实施方式，其中k-mer是重叠的并且增量为至少一个碱基或氨基酸，例如至少两个，如至少3个，例如至少4个，如至少5，例如至少6个或更多。

在优选的实施方式中，给定参考的完整基因组序列被分段成k-mer并上传到数据库中。还可以设想仅基于给定参考的转录组或给定参考的蛋白质组建立数据库。

如果目的只是为了识别来源序列的可能的参考，该数据库不必是完整的。提供对来自特定参考的基因组DNA的随机选择可能就足够。该选择也可以是非随机的，例如排除重复DNA以及所谓的垃圾DNA的拉伸。

对于生物序列、蛋白质、RNA、DNA每种类型，可以构建一个包含所有可用信息的数据库。在其他实施方式中，可以构建用于专门目的的专门的数据库，例如目的仅仅是为了从来源序列识别给定参考序列的存在或不存在。例如，该数据库可以包括来自人类、动物、哺乳动物、鸟类、鱼类、真菌、昆虫、植物、细菌、古细菌、病毒和/或质粒的序列信息。如果没有找到具有足够高的得分的匹配参考，也可以通过一个服务器向一个或几个其他服务器转发的关于读段的请求来建立数据库网络。

为了不损害速度而最佳利用硬件资源，该数据库可以被划分成存储在几个不同的服务器中的子数据库。

在其他实施方式中，根据选自门、纲、目、科、属和种的一种或多种分类学描述符，或一种或多种环境描述符，如来源、分布、起源和通常的搜索频率，使所述数据库组织成子数据库。

该数据库可以按照在图1中的描述来构建并使用被称为键值对存储的数据库引擎(例如BSDDB、KyotoCabinet、LevelDB、MongoDB和其他)存储。因此，在本发明的一个实施方式中，使用选自BSDDB、KyotoCabinet、LevelDB、MongoDB的组中的键值对存储来存储数据库。

算法的应用

本发明的方法和系统可以在需要识别在样品中发现的DNA的可能来源的多种应用中使用。

诊断

在药物治疗中，需要快速识别感染的可能来源。这可以使用根据本发明的方法来完成。由此可以选择将以最有效的方式和最少的副作用治疗感染的合适的治疗。

进一步的诊断应用涉及识别癌细胞中的病毒插入物。在本申请中，有利的是从获自原始读段的序列过滤全长人类序列或者简单地忽略在数据库中识别的所有人类命中。这将允许识别人基因组中的相对小的病毒插入物。

生物防御

在生物防御应用中，有必要进行快速可靠地识别所遇到的传染性物种或致病剂。本发明提供了对该来源没有先验知识时快速识别该来源的可能性。本发明的方法允许没有对病原体物种的先验知识时区别物种。

生物防御的其它应用包括识别转基因病原体，其中已经插入了例如有毒的转基因。该数据库有利也包含来自现有技术的质粒的序列信息。这将允许容易地识别插入物的侧翼区。如果转基因是来自在数据库中发现的生物体，它也将可能识别转基因的来源。在这种情况下，数据库可以返回病原体的名称、该转基因来源的生物体的名称、该转基因编码的基因以及用于插入该转基因的质粒。

食品安全和质量

目前用于鉴定食品中潜在的有害传染的方法缓慢(基于传染性生物体的分离和生长)或需要传染源的以前的知识(基于PCR的方法)。本方法不需要如此，并允许当局和制造商简单地分离基因组DNA，测序DNA并上传原始读段到能够操作本发明的方法的系统。

当在食物样品中寻找细菌、真菌或病毒时，有利的可能是查询只包含来自细菌、真菌或病毒的序列的数据库的一小部分。以这种方式，来自食物(蔬菜、水果、肉)的任何基因组序列都将被识别为在数据库中不存在，从而提高了该方法的性能。

其他应用包括质量控制。一种可能的应用是识别肉的种类，如肉末、肉酱、即食餐、方便食品。有许多企图诈骗的实例，其中昂贵的肉，如牛或羊肉被用更便宜的肉类，如猪肉替换或“稀释”。

其他可能的质量控制应用包括确定植物的品种，如葡萄、苹果、土豆等。

再其他的可能性包括水质的控制。

卫生和预防

本发明提供了通过能够快速鉴定取自与清洁程序相连的样品中的DNA来源来进行卫生控制的可能性。进一步的应用包括识别污染的可能来源从而能够应用最适于消除特定传染原的卫生技术。

项目

现在，本发明描述为任意编号的项目1到56，这将被视为本发明的实施方式。进一步参照所附的权利要求定义本发明。

1.一种识别生物序列的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自该子集的序列分段成k-mer，

c)针对包括参考序列的k-mer的数据库，查询来自所述子集的k-mer，

d)确定哪个(哪些)参考含有该k-mer，和

e)返回对可能的来源参考的描述。

2.项目1的方法，其中生物序列或短读段是氨基酸序列。

3.项目1的方法，其中所述生物序列或短读段是DNA或RNA序列。

4.任一前述项目的方法，其中k-mer查询涉及确定查询k-mer和参考k-mer之间的完全匹配。

5.任一前述项目的方法，其中所述查询进一步包括确定k-mer在参考序列中的位置。

6.任一前述项目的方法，其中存在和位置被用来确定在参考序列中的查询k-mer的连续性。

7.任一前述项目的方法，其中查询涉及从至少一个来源序列或短读段查询所有的k-mer，优选从至少50，例如从至少100，如从至少150，例如从至少200，如从至少250，例如从至少300，如从至少400，例如从至少500，如从至少750，例如从至少1000，例如从至少1500，如从至少2000，例如从至少2500，如从至少5000或更多个序列。

8.任一前述项目的方法，其中所述来源序列是至少50个碱基的核苷酸序列，优选至少100个碱基，如至少150个碱基，例如至少200个碱基，如至少250个碱基，例如至少300个碱基，如至少400，至少500或更多个碱基。

9.任一前述项目的方法，其中序列的子集包括至少1％的离散序列，例如至少2％，如至少4％，例如至少5％，如至少6％，例如至少7、5％，例如至少10％，如至少15％，例如至少25％，如至少30％，例如至少35％，如至少40％，例如至少50％。

10.任一前述项目的方法，进一步包括选择序列的一个或多个另外的子集并将它们进行项目1的步骤a)至步骤e))。

11.任一前述项目的方法，其中所述子集是随机的或过滤的。

12.任一前述项目的方法，其中k-mer的大小为4、8、12、16、20、24、28、32、36、40、44、48、52、56、60、64或更长。

13.任一前述项目的方法，其中k-mer是连续的。

14.任一前述项目的方法，其中k-mer是重叠的并且增量为至少一个碱基或氨基酸，例如至少两个，如至少3个，例如至少4个，如至少5个，例如至少6个或更多。

15.任一前述项目的方法，其中k-mer是不相交的子序列的串联。

16.任一前述项目的方法，其中针对数据库查询来自给定序列的k-mer以确定k-mer在一个或多个参考序列中的存在和k-mer在所述一个或多个参考序列中的位置。

17.项目16的方法，其中只在k-mer存在时才查询位置。

18.任一前述项目的方法，其中为返回的参考计算得分。

19.任一前述项目的方法，其中为识别的参考序列计算得分，所述得分与在给定的参考序列中发现的来自一个或多个序列的k-mer数目相关联。

20.任一前述项目的方法，其中为识别的参考计算得分，所述得分通过在参考序列中发现的来自一个或多个序列的k-mer的平均局部集中度与连续性或近似连续性相关联。

21.任一前述项目的方法，其中为识别的参考计算得分，所述得分与在参考序列中的k-mer数目相关联，所述k-mer也存在于来自所述来源的k-mer的子集中。

22.任一项目18至21的方法，其中根据所述得分对可能的来源参考排名。

23.任一前述项目的方法，其中查询来自一个来源序列或短读段的所有k-mer并且为所述来源序列或短读段计算一个或多个得分。

24.项目23的方法，还包括：查询来自第二来源序列或短读段，优选来自第三来源序列或短读段的所有k-mer等。

25.任一前述项目的方法，其中一旦以预定义的统计概率识别出一参考生物体，那么可以停止数据库查询。

26.任一前述项目的方法，其中如果在数据库中没找到预定分段的k-mer，那么可以停止数据库查询。

27.任一前述项目的方法，其中数据库输出关于一个或多个可能参考的下列信息中的一个或多个：可能参考的分类学名称、所述可能参考的近亲、所述参考的来源、一组另外相关的参考。

28.任一前述项目的方法，其中所述数据库输出最可能的参考的序列，优选其中所述数据库输出最可能的参考物种的全基因组序列。

29.任一前述项目的方法，其中来自具有非常相似的序列的参考的结果或者来自另外相关的参考的结果在输出中分组。

30.任一前述项目的方法，其中执行该方法的几次迭代，如在第一次迭代中，识别最丰富的参考并从所述来源序列或短读数中去除来自所述最丰富的参考的序列。

31.项目30的方法，还包括在第二次迭代中，识别第二最丰富的参考，去除来自所述第二最丰富的参考的序列等。

32.项目30的方法，还包括在第二次迭代中，识别插入物的可能参考。

33.任一前述项目所述的方法，该方法还包括最初去除与来自预定义的参考的序列对准的来源序列。

34.任一前述项目的方法，其中所述方法包括如果数据库中不存在来自一个来源序列或短读段的预定义数目的k-mer，那么从所述来源序列或短读段忽略k-mer。

35.任一前述项目的方法，其中查询涉及忽略来自一个或多个预定义的参考的k-mer。

36.任一前述项目的方法，其中当从核酸测序仪获得原始序列时查询所述原始序列。

37.一种包括有参考序列的k-mer的数据库，所述数据库包括：

a)来自参考序列的k-mer的第一集合，和

b)每个k-mer在参考序列中的位置的第二集合。

38.项目37的数据库，其中所述数据库还包括关于与给定参考相关联的全长序列的信息，和/或所述参考的来源，和/或所述参考的一个或多个分类学描述符。

39.项目37-38的数据库，其中在数据库中的k-mer进行哈希函数运算，向每个独一无二的k-mer分配唯一密钥。

40.任一项目37-39的数据库，其中在第一集合中的每个独一无二的k-mer由矢量与关于有k-mer存在的那些参考的信息相关联。

41.任一项目37-40的数据库，其中在第二集合中的每个独一无二的k-mer由矢量与关于该k-mer存在时，其在每个参考中的位置信息相关联。

42.项目37-41的数据库，其中k-mer的长度为4、8、12、16、20、24、28、32、36、40、44、48、52、56、60、64或更长。

43.任一项目37-42的数据库，其中k-mer是非重叠的。

44.任一项目37-43的数据库，其中k-mer是重叠的并且增量为至少一个碱基或氨基酸，例如至少两个，如至少3个，例如至少4个，如至少5个，例如至少6个或更多。

45.任一项目37-44的数据库，其中数据库包括来自每个参考的完整序列的k-mer。

46.任一项目37-46的数据库，其中所述数据库包括来自人类、动物、哺乳动物、鸟类、鱼类、真菌、昆虫、植物、细菌、古细菌、病毒和/或质粒的序列信息。

47.任一项目37-46的数据库，其中所述数据库被分成存储在几个不同服务器中的子数据库。

48.任一项目37-47的数据库，其中根据选自门、纲、目、科、属和种的一种或多种分类学描述符，或一种或多种环境描述符，如来源、分布、起源和过去查询的频率，使所述数据库组织成子数据库。

49.一种用于识别来源序列的可能来源的数据处理系统，该系统包括输入装置、中央处理单元、存储器和输出装置，其中，所述数据处理系统中存储有表示在执行时使得项目1-36的方法得以被执行的系列指令的数据，所述存储器还包括根据任一项目37-49的数据库。

50.项目49的系统，其中所述数据库存储在服务器中，并且输入和输出装置是客户端，所述客户端和服务器经由数据通信连接来连接。

51.任一项目49-50的系统，其中所述客户端选自个人计算机、台式PC、便携式PC、手持式计算装置，如智能电话。

52.任一项目49-51的系统，其中所述客户端包括使客户端采样来源序列的子集，将这些来源序列分段成k-mer，并将这些来源序列传送到服务器的一系列指令。

53.项目49-52的系统，所述客户端还包括允许其基于从服务器传送到客户端的序列，执行将来源序列组装成一个或多个更大的序列的一系列指令。

54.任一项目49-53的系统，通过数据连接与测序设备相连接。

55.一种计算机软件产品，含有在执行时使得项目1-36的方法得以被执行的系列指令。

56.一种集成电路产品，含有在执行时使得项目1-36的方法得以被执行的系列指令。

实施例

具有k-mer的序列的快速识别

这里，我们提出了新方法，Tapir，即能够迅速地指向DNA或RNA的可能来源并能够直接在从DNA测序仪获得的原始读段上工作。我们的系统包括参照已知DNA的服务器，以及具有待量化的DNA数据的客户端。为了说明如何使用，我们已参考了数千细菌基因组、噬菌体基因组、噬菌体和质粒，以及人类基因组、小鼠基因组、拟南芥和来自真菌、古细菌的各种序列。我们也已经在Web浏览器上实施了客户端运行，并能够处理来自便携式计算装置的数据中的数亿个碱基的数据。该方法依赖于索引k-mer，以及向服务器传输有限量的数据。它能够在几秒钟内从Android智能手机完成其任务，消耗适度数量的带宽与服务器通信，并且尽我们的知识提供了一种不同于任何现有工具的使用简便性。它在我们的核心设施上使用，以测序运行来进行常规的即时质量检查，并且可获自http://tapir.cbs.dtu.dk。

介绍

DNA测序在过去十年已经越来越实惠[13]，扼要的讲对DNA测序进行叙述再次成为绝对平常化。今天的高端测序仪具有处理相当于几个人类基因组或几百个细菌的能力，并且下一代测序仪已经开始变得可用，因此需要低得多的初始投资，并提供灵活的测序量。完整细菌分离株的测序是一天的事，并且很快就是几小时的事。最近发布的纳米孔测序[12]呈现了一个USB供电的装置，能够直接测序DNA，以及因为测序装置将是一次性的，故而资本投资水平前所未有的低。Oxford Nanopore，在该未来产品背后的公司已经在2012年发布新款[8]。提取DNA是相对简单的过程，并且可以预见，DNA测序将很快在分子生物学中成为常规和廉价的过程。患者将被常规测序，通过他们的DNA追踪传染原的暴发，水和食品的质量也通过DNA测序监测。

在分析方面，利用如Smith-Waterman算法[14]的开拓性的工具进行序列的局部比对已经成为生物信息学的基石。一旦将它应用到查询集合和参考集合之间，将允许对比对排名，使研究者对新测序的DNA或RNA从它与现有序列的相似性来推断其起源和功能。虽然有批评该方法有时是不准确的[2，11]，但其受欢迎程度仍然是不争的，公共数据库中有大量的功能注释提到“通过序列同源性”。然而，在数据库中实现将新获得的DNA与存档的现有数据库比对仍是相对苛刻的计算任务。BLAST[1]和随后的BLAT[5]提高了速度，但在网络搜索引擎几乎立即返回结果的时代，用目前可用的序列数量针对已知序列的池搜索新序列可能需要相当长的时间。设计用于短读段测序的新工具自从被开发出来也仅仅命名了两个，如Bowtie[6]和BWA[7]，但这些工具是被设计成针对给定参考来比对所有测序读段。为了实现速度，这样的工具将参考的索引加载到存储器中，并以此限定可以被处理的参考DNA的量。

我们注意到在找到查询序列和参考的集合之间的绝对最佳比对的计算需求任务与从一组查询序列快速识别与它们最匹配的参考之间的距离。据我们所知，没有任何一种工具采用短DNA或RNA序列集，如来自DNA测序仪的读段，并返回该集所代表的参考的列表，无论是全基因组或个体基因。要做到这一点，我们提出了与在BLAT和SSAHA[9，10]两者中使用来自比对种子以及在MUSCLE[3]中的k-mer计数不同的方式使用k-mer，以几秒或者更短时间相当准确地识别DNA序列的来源。

材料和方法

将获自EBI和NCBI的公众可得的基因组、重叠群、质粒和个体基因下载作为参考DNA。每个参考序列被分成重叠k-mer(on-overlapping k-mer)并且对于遍布所有参考的所有k-mer，创建键值对存储或NoSQL数据库(使用KyotoCabinet[4])，关联到每个k-mer体(数据库中的密钥)，对应于具有该k-mer的参考的标识符列表(图1)。称它为存在数据库。同样地，在参考中发现该k-mer的位置存储在被称为位置数据库的地方(图1)。参考标识符和信息，如描述行和数据源之间的关联性，储存在单独的SQL数据库中。

为了对一组短查询序列或读段评分，对它们的随机样品进行迭代(图2)。对于每个序列，对通过在整个序列上滑动宽度k的窗口所获得的连续k-mer进行迭代。对于每个k-mer，如果它之前并未被计数并且在存在数据库中发现它，那么将查询参考的位置。一旦对读段的所有k-mer进行处理，查找在参考中匹配的连续位置的数目并且只考虑最大匹配簇，即，在所有匹配参考中源自相同读段的匹配k-mer的最大集中度。对于每个这样的簇，将k-mer数加上可能之前用于该参考的加数并且更新已经计数的k-mer列表。然后处理下一个序列或读段。获得与发现匹配的k-mer计数相关联的参考列表。对于每对<参考，计数>，计数除以查询集中独一无二的k-mer的数量，提供了在给定参考匹配的被查询子集中DNA量的粗略得分。如果查询集完全匹配该序列，那么得分将是1，否则会降低；例如，如果查询集是两个参考的等比例混合物，那么对两个参考的得分都将是0.5左右。该计数也可以除以参考的大小(参考序列中的独一无二的k-mer数目)，得到由该查询表示的参考部分的粗略得分；该第二得分有助于对匹配参考进行分选，并避免偏向最大参考。最终的得分被计算为这两个得分的加权和，默认是相等的加权。如果查询集很大，例如，如果考虑来自一轮DNA测序的所有读段，那么仅使用该集的随机样品。

为了方便该服务使用，实施作为在Web浏览器中的网页运行的HTML5/JavaScript客户端。在编写时，Firefox 15.0是实现所有需要的功能的唯一浏览器，并测试在Linux、Mac OS X、微软的Windows和Android 4.0上的工作。

为了基准化该原本设计为识别测序数据中的细菌的系统，反复取得在2012年可从EBI获得的细菌的所有序列，也就是747个细菌基因组。对于每一个基因组，从基因组序列产生随机可能重叠的子序列，以模拟从DNA测序仪获得的读段；使用长度为50、100、150、200和250个碱基的子序列。还以0％(没有错误)、1％、5％和10％的比率引入了碱基的均匀的随机置换，以模拟一类测序错误和在实际样品中的准时突变的存在。对于每个基因组，长度和置换率，取100个子序列或读段的随机样品并且该采样重复十次。

结果

对于每一种细菌基因组，采取了100个随机模拟读段，并使用我们的方法，在参考中，针对于包括那些细菌基因组的数据库对它们评分，记录该查询基因组在25个最佳得分列表中的排名。平均排名和排名的标准差如图4所示。

平均排名越接近1，评分越好，而排名的标准差越小，对采样影响越不敏感。写入每个单个面板的丢失的排名数，对应于没有在25个最高得分的基因组数量。

读段在长度为50个碱基时性能低于最佳，但在100个碱基的读段时已经有显著的改善，所查询基因组在97％和99％的时间具有低置换率，在前5位具有较低置换率而在前15位时具有更高的置换率。读段长度增加至250个碱基有助于补偿更高置换率对平均排名的负面影响。

使用的长度范围和置换率与从新一代测序平台，例如lllumina(100个碱基，具有约0.1-1％的错误率)，Life Technologies的SOLiD 5500(75nt读段，具有0.01％的错误率)，Ion Torrent PGM(200-300个碱基，具有1％的错误率)，或Pacific Bioscience(3,000个碱基，具有15％的错误率)获得的范围是可比的。我们的方法在这些范围内执行良好，并且我们预期因增加对双末端测序(一种用于提供替代较长读段的技术)的支持而实现了性能进一步提高。我们的方法似乎对于测序错误(例如碱基置换)相对不敏感，并且对于我们的测试查询的预期较低的排名随置换率增加而使影响最小。

得益于NoSQL数据库的使用，随着基因组数据越来越丰富，预计规模会扩大，并且不断能够在比较实惠的电脑系统上索引和查询越来越大的参考集合。

为了便于使用我们的方法，开发了一种基于浏览器的客户端。我们用高达2GB的原始FASTQ文件测试，并监控到其在RAM中只用了稍微超过200MB并在20秒内返回结果。

结论

基于TAPIR的概念相当简单。已经宣布了DNA数据库规模的增长，并至少观察了十多年，但最近在DNA测序技术的发展使得快速且经济实惠的数据生成成为现实。我们主张对所有已知的DNA的匹配实验获得的DNA序列是生物信息学中最重要的挑战之一。我们在这里表明，这可以用与互联网网络搜索巨头已经让公众使用相匹配的速度和易用性完成。当考虑用桌面DNA测序仪进行诸如实时监控、患者感染、生物防御或食品安全的任务时，我们的方法提供了一种立即早期步骤，在此期间搜索空间可以缩小并且此后可以执行更先进的分析方法。

实施例2

在本实施例中，参考了来自细菌、病毒、噬菌体、质粒、以及人类、小鼠、植物、真菌和古细菌的数以万计的基因组和基因组区域。还实施了在Web浏览器上运行的客户端，并且证明采用该客户端在几秒钟内处理并识别来自商品化便携式计算装置的千兆字节的原始测序数据，同时消耗适度量的带宽与服务器通信。因此，在本实施例中，示出了来自原始读段的DNA的识别可以与查询搜索引擎一样简单。

查询DNA序列针对参考的全面集合的匹配集

查找比对程序的主观方式是将它们分成两个主要类别：一类试图针对已知参考的集合映射一个查询序列(如BLAST)，和一类试图尽可能快地映射针对一种指定参考的大量短序列(例如，Bowtie或BWA)。我们提出了一个可对大量短序列识别很好的参考的折衷办法；我们匹配了针对于参考序列的集合的若干序列并且选出了哪些参考在查询集中得到最好的体现。

在本实施例中提出的方法不涉及在索引k-mer期间的任何选择步骤，该功能从序列的集合构建时极大地简化了复杂性。这以空间为代价而得到，有潜在的较少信息的k-mer被索引，但是这通过以下好处得到弥补：该过程对于参考集合的总大小是线性的，并且可以平行进行。这使得索引所有已知的DNA最终似是可信的(类似于在互联网上的所有文件的网络搜索引擎的索引。)

在这个例子中，我们的算法并不仅仅是计数k-mer，并且它不执行完整的映射或比对。该算法考虑了在每个读段内的匹配k-mer，以及彼此接近的匹配k-mer的簇。

在本实施例中，使用非重叠的k-mer进行索引，同时在查询中使用重叠的k-mer，如图5所示，但我们认为这是一个实施细节并且可以很容易地使用重叠的k-mer进行索引和在查询中的不重叠的k-mer同时对给定分数保持相同的指导原则以匹配参考。

如果将树或哈希用于k索引和查找，那么在使用k-mer索引的大小为u的参考中确定长度为p的串n次出现的位置的时间复杂度为O(p+n log u)或O(p+n)复杂度。

当为了诊断目的要识别DNA数据的查询集，如来自测序仪的原始读段时，我们认为包括针对于全面的参考数据库映射所有读段的蛮力方法有两个主要的缺点：成百兆或千兆字节的数据的大部分或者从测序设施传送到计算中心，并且执行该任务必须的计算资源是至关重要的。假设一个参考集合包含10,000大肠杆菌规模的细菌并且需要30秒优化比对器，如BWA和bowtie2来处理250Mb的原始测序数据(如果基因组是4Mb大小，平均覆盖为约60x)，尽管这可以在多个CPU上并行进行处理，但是这样的CPU处理会花费3天半。

除了时间复杂度，数据传输将是250Mbases的DNA，将测序数据传送到存放有参考的数据中心。我们基于k-mer的方法减少了对小的参考集的详细的调查，如映射读段或SNP判读，甚至基于模板的从头组装。当评估性能时，随意选择，如果正确的答案是在5个建议的匹配集中，那么最初只是考虑搜索成功。针对那些参考映射所有读段的以准确识别哪个是最佳匹配这项任务可以在12分钟内在相同的CPU上完成，或者如果在规定的3天半的时间内获得了强大的多核体系结构，那么任务可以在更短时间内完成。传输全部的基因组将代表约20 Mbases的DNA，它可以很容易地通过3G移动互联网连接来完成。我们的方法使得移动测序设施，如Ion总线[15]能够在现场的偏远位置执行关键的诊断和科学任务。如果有未映射的读段，因为存在较小的区域，如质粒、毒力基因、病毒或细菌的混合物，这些读段可以被类似地处理，并且全部内容通过几次迭代来识别。

建立基准

为了基准化该原本设计为识别测序数据中的细菌的系统，反复取得在2012年可从EBI获得的细菌的所有序列，也就是747个细菌基因组。除了这些以外，所含的所有的参考数据库为：来自NCBI的细菌参考、噬菌体和病毒、质粒以及人类基因组(见下文表1)。表1示出在2012年开始的基因组参考(参考的来源和数目)的快照。参考是完整基因组或质粒的混合物，和诸如重叠群或基因的基因组片段的混合物。

表1.基因组参考

数据库	参考数	大小(DNA碱基)
			HIV	4053	36471153
噬菌体基因组(Sanger)	1078	59538128
			病毒基因组(Sanger)	3464	64859892
细菌基因组(Sanger)	747	2418028337
			细菌基因(NCBI)	5218077	4963568551
细菌基因组(NCBI)	4693	2418028337
			病毒基因组(NCBI)	1750	60637755
真菌	202270	298736207
			人类微生物组	1653700	1490442185
质粒	159705	132800479
			病毒(Virii)	78630	65110952
智人(Hg19)	3134	2844000504
			小鼠(Mus musculus)	305	2745142291
植物(RefSeq)	558267	8622349159
			无脊椎动物(Genbank)	1123813	18429666992
原生动物(Genbank)	47275	1997449553
			真菌(Genbank)	200	242402709

对于每种基因组，从基因组序列产生随机可能的重叠子序列，以模拟从DNA测序仪获得的读段；使用长度为50、100、150、200和250个碱基的子序列。还以0％(没有错误)、1％、5％和10％的比率引入了均匀的碱基随机置换，以便模拟一类测序错误和实际样品中的准时突变的存在。对于每种基因组、长度和置换率，执行100个子序列或读段的随机样品，并且采样重复5次。

目的是评估当考虑到诸如测序错误或突变的不确定性时，在样品或足够接近的基因组中是否可以找到有哪种已知的DNA。

预测性能

对每个细菌基因组，采取100个随机模拟的读段，并使用我们的方法，针对包括那些细菌基因组的数据库，从其他细菌、噬菌体、植物、真菌、病毒和哺乳动物序列和基因组的较大集合中，对它们评分，记录查询基因组中25个最佳匹配的参考列表中的排名。为了评估对每种测试细菌基因组的结果的变异性，这对每个基因组重复5次并且平均排名和对该排名的标准偏差如图9所示。

读段在长度为50个核苷酸时性能相对较差，但读段长度增加时观察到了显著提高，在测序碱基中的长度为100的读段已经接近最大性能。最好的结果是示出了正确的基因组在97％的时间是在低错误率的结果列表中，在前5位具有较低置换率，并且在前15位具有较高置换率。增加读段长度达250个碱基帮助补偿错误率增加带来的负面影响。增加被发送用于识别的随机样品中的读段数并没有产生太大的影响，参见图7中的随机样品中：100个读段是少量的数据，但它在大量情形下似乎足以识别DNA。

如前面所详述，我们的方法旨在返回在提出的匹配集中的正确参考并且通过如此做，简化了粗苯的方法需要利用计算指令程序来探索的搜索空间。将我们限制为在前五个结果中找到查询序列很可能比需要的还要严格，因为运行所有25个分析相比于穷举搜索仍将是至关重要的，但是指出该方法已经能够返回在非常小的候选答案集中的正确答案。

在迭代搜索和识别的情况下，可以考虑指出正确的细菌物种，即使不是正确的精确株或基因组参考，已经是比较成功的答案。图6示出了以超过个核苷酸的读段进行的识别过程执行得很好。

使用的长度范围和置换率与从新一代测序平台，例如lllumina(150个碱基，具有约0.1-1％的错误率)，Life Technologies的SOLiD 5500(最大75nt读段，具有0.01％的错误率)，Ion Torrent PGM(最大200-300个碱基，具有1％的错误率)，或Pacific Bioscience(3,000个碱基，具有15％的错误率)获得的范围是可比的。我们的方法在这些范围内执行良好，并且我们预期因增加对双末端测序(一种用于提供替代较长读段的技术)的支持而实现了性能进一步提高。我们的方法似乎对于测序错误(例如碱基置换)相对不敏感，并且对于我们的测试查询的预期较低的排名随置换率增加而使影响最小。

我们也尝试基于来自Ion Torrent PGM的、从病毒和细菌分离株到宏基因组学混合物的范围的测序数据的方法。索引的参考的集合中的非常相似的基因组，如几株同一物种，通过增加密切相关的基因组比正确参考的基因组有更低的排名的概率可以有助于性能的劣化。这是通过考虑物种增加的性能而不是精确的参考所确认的，并且这是中度不便的，即可以在第二迭代期间消除歧义。最后，因为我们已经考虑读段范围内的k-mer，而不是孤立的实体，我们从来自不同哺乳动物的样品测序获得了非常令人鼓舞的结果，并预测能够在不久的将来可靠地识别它们。

计算性能

服务器：

在服务器上的存储器使用量可以保持最低限度使用基于磁盘的键值存储和调优的性能可以通过这些缓存到存储器中运行可在计算机上实现。由于使用的NoSQL数据库，我们也期望能规模达的基因组数据得到越来越丰富，并不断能够引用索引和查询越来越大集合比较实惠的电脑系统上。与当前的实现既索引系统和服务器在Python中实现，在使用8芯(英特尔至强，2.93GHz的)几个小时来执行的参照DNA 44Gbases的索引，和一个输入样品的处理以几秒钟。甲显著加速可与优化的努力来实现，例如移动到C的瓶颈，但它也可以由专更多内核增加的更多的请求的处理全球演出，在需要变得显而易见。

客户：

为了方便使用我们的方法，开发了使用JavaScript和HTML5特征的基于浏览器的客户端，其可以在http://tapir.cbs.dtu.dk访问。该客户端目前正在最新的Firefox版本(版本15或更高)上运行。

随着具有在主频为2.53GHz的Intel Core i5CPU的相对普通的笔记本电脑上Firefox，可在30秒内处理大小高达2Gb的FASTQ文件的原始读段，文件越小越快，使用在RAM中稍小于300Mb，并且与服务器通信几秒。

还实施了基于控制台的命令行工具来执行我们的算法和随后的比对。该实施可获自流行的软件库https://bitbucket.org/lgautier/dnasnout-client。该实施使用我们的算法以抓取参考基因组，并用bowtie2进行所有读段的索引和映射。当考虑前10个读段时，完整的迭代花费不到一分钟并且在98％的情况下一次迭代就足够了。随着浏览器的快速发展，预计很快就能仅使用网络浏览器开展与流行病学实验室用台式测序所作相似的工作流程。

讨论

我们主张针对所有已知的DNA匹配实验获得的DNA序列是生物信息学中最重要的挑战之一。我们在这里表明，这可以用与互联网网络搜索巨头已经使公众习惯相匹配的速度和易用性来完成。当考虑诸如实时监控，患者感染、生物防御或食品安全的任务时，今天的桌面DNA测序如Ion Torrent PGM或lllumina MiSeq已经不能胜任此任务并且我们的方法提供了一种立即早期步骤，在此期间搜索空间可以缩小而且可以事后在本地执行更先进的分析方法，而不需要在执行DNA测序的实验室和计算设施之间传输大量的原始数据。

方法

基因组参考的来源：

下载公众可获得的基因组、重叠群、质粒和可从EBI和NCBI获得的个体基因作为参考DNA。参考的确切组成将随着时间扩大，但在表1中列出了用于本实施例的快照。

参考索引：

每个参考序列被分成不重叠的k-mer并且对于所有参考的所有k-mer，创建键值对存储或NoSQL的数据库(使用KyotoCabinet[4])，关联每个k-mer(在数据库中的密钥)，对应于具有该k-mer的参考的标识符列表。称这为存在数据库。同样地，在该参考中发现该k-mer的位置被存储在称作位置数据库的地方。k被选择为等于16，因为它给出了满意的结果，并且作为4的倍数很好地适用于位包装。参考标识符和信息，如描述行和数据源之间的关联，储存在单独的SQL数据库。

评分：

为了对短查询序列或读段集评分，对它们的随机样品迭代。样品大小越大，可靠的准确性越大。对于每个序列，对在通过在序列上滑动宽度k的窗口所获得的连续k-mer迭代。对于每个k-mer，如果它之前并未被计数并且在存在数据库中发现它，那么将查询参考的位置。一旦对读段的所有k-mer进行处理，查找在参考中匹配的连续位置的数目并且只考虑最大匹配簇，即，在所有匹配参考中源自相同读段的匹配k-mer的最大集中度。对于每个这样的簇，将k-mer数加上可能之前用于该参考的加数并且更新已经计数的k-mer列表。然后处理下一个序列或读段。当处理了所有读段后，获得与发现匹配的k-mer计数相关联的参考列表。对于每对<参考，计数>，计数除以查询集中独一无二的k-mer的数量，提供了在给定参考匹配的被查询子集中DNA量的粗略得分。利用示出的评分原则，如果查询集完全匹配该序列，那么得分将是1，否则会降低；例如，如果查询集是两个参考的等比例混合物，那么对两个参考的得分都将是0.5左右。该计数也可以除以参考的大小，得到由该查询表示的参考部分的粗略得分；该第二得分有助于对匹配参考进行分选，并避免偏向最大参考。最终的得分被计算为这两个得分的加权和，其中使用相等的加权。如果查询集很大，例如，如果考虑来自一轮DNA测序的所有读段，那么仅使用该集的随机样品。

客户端的实施：

为了方便使用该服务，实施在Web浏览器的页面运行的HTML5/JavaScript客户端。对于目前的研究，使用Firefox 15版本，并且测试它在Linux、Mac OS X、微软Windows(各种笔记本和台式机)以及在Android 4.0(平板ASUS TF101-预计它也将在高端智能手机上工作)的运行。然而，本领域技术人员将理解，其他合适的浏览器也可以是有用的。该客户端也被实施为Python库和命令行工具以便于在现有的工作流程和管线中评估和整合。

其它技术规格：

除了结合到诸如KyotoCabinet库以外的所有实施是在服务器端利用Python版本2.7.3进行。网络应用是使用微框架Flask并由lighttp服务。为Python版本3.3开发客户端库和命令行工具。

本领域技术人员将理解，算法或部分算法的实施可以由其它合适的和一般公知的编程语言来实现，例如C编程语言，其可以通过降低用于查询的时间，提高该方法的性能。

参考文献

[1]Stephen F.Altschul,Warren Gish,Webb Miller,Eugene W.Myers,andDavid J.Lipman.Basic local alignment search tool.Journal of Molecular Biology,215(3):403-410,October 1990.

[2]Damien Devos and Alfonso Valencia.Practical limits of functionprediction.Proteins:Structure,Function,and Genetics,41(1):98-107,October2000.

[3]R.C.Edgar.MUSCLE:multiple sequence alignment with high accuracyand high throughput.Nucleic Acids Research,32(5):1792-1797,March 2004.

[4]Mikio Hirabayashi.Kyoto cabinet:a straightforward implementation ofDBM.

http://fallabs.com/kyotocabinet/

[5]W.J.Kent.BLAT—The BLAST-Like alignment tool.Genome Research,12(4):656-664,March 2002.

[6]Ben Langmead,Cole Trapnell,Mihai Pop,and Steven L Salzberg.Ultrafast and memory-efficient alignment of short DNA sequences to the humangenome.Genome Biology,10(3):R25,2009.

[7]H.Li and R.Durbin.Fast and accurate short read alignment withburrows-wheeler transform.Bioinformatics,25(14):1754-1760,May 2009.

[8]Christopher E Mason and Olivier Elemento.Faster sequencers,largerdatasets,new challenges.Genome Biology,13(3):314,2012.

[9]Z.Ning.SSAHA:a fast search method for large DNA databases.GenomeResearch,11(10):1725-1729,October 2001.

[10]Zemin Ning,W.Spooner,A.Spargo,S.Leonard,M.Rae,and A.Cox.The SSAHA trace server,pages 519-520.IEEE.

[11]Burkhard Rost.Enzyme function less conserved than anticipated.Journal of Molecular Biology,318(2):595-608,April 2002.

[12]Nicole Rusk.Cheap third-generation sequencing.Nature Methods,6(4):244-244,April 2009.

[13]Jay Shendure and Hanlee Ji.Next-generation DNA sequencing.NatureBiotechnology,26(10):1135-1145,October 2008.

[14]T.F.Smith and M.S.Waterman.Identification of common molecularsubsequences.Journal of Molecular Biology,147(1):195-197,March 1981.

[16].Rumble SM,Lacroute P,Dalca AV,Fiume M,Sidow A,et al.(2009)SHRiMP:accurate mapping of short color-space reads.PLoS ComputationalBiology 5:e1000386.

[17].Li H,Homer N(2010)A survey of sequence alignment algorithms fornext-generation sequencing.Briefings in Bioinformatics 11:473-483.

[18].Babraham bioinformatics-FastQ screen.

http://www.bioinformatics.babraham.ac.uk/projects/fastq screen/.URL

http://www.bioinformatics.babraham.ac.uk/projects/fastq_screen/.

Claims

1.一种识别生物序列，如短读段，的可能来源的方法，该方法包括：

a)从一来源采样序列或短读段的子集，

b)将来自所述子集的序列分段成k-mer，

e)确定哪个/哪些参考含有所述一个或多个k-mer，和

f)返回对可能的来源参考的描述，

其中所述包括参考序列的k-mer的第一集合与包括k-mer在参考序列中的位置的第二集合是分离的。

2.如权利要求1的方法，其中所述方法不涉及对序列数据使用比对算法，例如采用评分矩阵的比对算法。

3.如任一前述权利要求的方法，其中所述查询进一步包括确定k-mer在所述参考序列中的位置。

4.如任一前述权利要求的方法，其中存在和位置被用来确定查询k-mer在参考序列中的的连续性。

5.如任一前述权利要求的方法，其中所述生物序列是氨基酸序列。

6.如权利要求1-4的方法，其中所述生物序列是DNA或RNA序列。

7.如任一前述权利要求的方法，其中k-mer查询涉及确定查询k-mer和参考k-mer之间的完全匹配。

8.如任一前述权利要求的方法，其中查询涉及从至少一个来源序列或短读段查询所有k-mer，优选从至少50，例如从至少100，如从至少150，例如从至少200，如从至少250，例如从至少300，如从至少400，例如从至少500，如从至少750，例如从至少1000，例如从至少1500，如从至少2000，例如从至少2500，如从至少5000或更多个序列。

9.如任一前述权利要求的方法，其中所述源序列是至少50个碱基的核苷酸序列，优选至少100个碱基，如至少150个碱基，例如至少200个碱基，如至少250个碱基，例如至少300个碱基，如至少400，至少500或更多个碱基。

10.如任一前述权利要求的方法，其中序列的子集包括至少1％的离散序列，例如至少2％，如至少4％，例如至少5％，如至少6％，例如至少7、5％，例如至少10％，如至少15％，例如至少25％，如至少30％，例如至少35％，如至少40％，例如至少50％。

11.如任一前述权利要求的方法，进一步包括选择序列的一个或多个另外的子集并将它们执行权利要求1的步骤a)至步骤f)。

12.如任一前述权利要求的方法，其中所述子集是随机的或过滤的。

13.如任一前述权利要求的方法，其中k-mer的大小为4、8、12、16、20、24、28、32、36、40、44、48、52、56、60、64或更长。

14.如任一前述权利要求的方法，其中k-mer是连续的。

15.如任一前述权利要求的方法，其中k-mer是重叠的并且增量为至少一个碱基或氨基酸，例如至少两个，如至少3个，例如至少4个，如至少5个，例如至少6个或更多。

16.如任一前述权利要求的方法，其中k-mer是不相交的子序列的串联。

17.如任一前述权利要求的方法，其中针对数据库查询来自给定序列的k-mer以确定k-mer在一个或多个参考序列中的存在和k-mer在所述一个或多个参考序列中的位置。

18.如权利要求17的方法，其中只在k-mer存在时才查询位置。

19.如任一前述权利要求的方法，其中为返回的参考计算得分。

20.如任一前述权利要求的方法，其中为识别的参考序列计算得分，所述得分与在给定的参考序列中发现的来自一个或多个序列的k-mer数目相关联。

21.如任一前述权利要求的方法，其中为识别的参考计算得分，所述得分通过在参考序列中发现的来自一个或多个序列的k-mer的平均局部集中度与连续性或近似连续性相关联。

22.如任一前述权利要求的方法，其中为识别的参考计算得分，所述得分与在参考序列中的k-mer数目相关联，所述k-mer也存在于来自所述来源的k-mer的子集中。

23.如任一权利要求19至22的方法，其中根据所述得分对可能来源参考排名。

24.如任一前述权利要求的方法，其中查询来自一个来源序列或短读段的所有k-mer并且为所述来源序列或短读段计算一个或多个得分。

25.如任一前述权利要求的方法，其中为参考序列获得匹配k-mer的计数。

26.如任一前述权利要求的方法，其中通过将参考序列的匹配k-mer的计数除以被查询子集中的独一无二的k-mer的数目来获得得分。

27.如权利要求24-26的方法，其中通过将参考序列的匹配k-mer的计数除以参考序列的大小来获得得分。

28.如权利要求24-27的方法，其中参考序列的得分被计算为权利要求26和27的得分的加权和。

29.如任一前述权利要求的方法，还包括查询来自第二来源序列，优选来自第三来源序列的所有k-mer。

30.如任一前述权利要求的方法，其中一旦以预定的统计概率识别出参考生物体，就可以停止数据库查询。

31.如任一前述权利要求的方法，其中如果在数据库中没找到预定分段的k-mer，就可以停止数据库查询。

32.如任一前述权利要求的方法，其中所述数据库输出关于一个或多个可能参考的下列信息中的一个或多个：序列的任何注释、编码序列、调控序列、可能参考的分类学名称、所述可能参考的近亲、所述参考的来源、一组另外相关的参考、参考是从哪里获得的(例如土壤、海洋、内脏或下水道)、参考序列是何时获得的、分类学分类、近缘物种、关于参考序列是从哪个数据库(例如NCBI或EBI/Sanger数据库)下载的信息。

33.如任一前述权利要求的方法，其中所述数据库输出最可能的参考的序列，优选其中所述数据库输出最可能的参考物种的全基因组序列。

34.如任一前述权利要求的方法，其中来自具有非常相似的序列的参考的结果或者来自另外相关的参考的结果在输出中分组。

35.如任一前述权利要求的方法，其中执行该方法的几次迭代，如在第一次迭代中，识别最丰富的参考并从所述来源序列或短读数中去除来自所述最丰富的参考的序列。

36.如权利要求35的方法，还包括在第二次迭代中，识别第二最丰富的参考，去除来自所述第二最丰富的参考的序列等。

37.如权利要求36的方法，还包括在第二次迭代中，识别插入物的可能参考。

38.如任一前述权利要求的方法，该方法还包括最初去除与来自预定的参考的序列对准的源序列。

39.如任一前述权利要求的方法，其中所述方法包括如果数据库中不存在来自一个来源序列的预定义数目的k-mer，就从所述来源序列忽略k-mer。

40.如任一前述权利要求的方法，其中查询涉及忽略来自一个或多个预定的参考的k-mer。

41.如任一前述权利要求的方法，其中当从核酸测序仪获得原始序列时查询这些原始序列。

42.如任一前述权利要求的方法，其中使用适应采样。

43.一种在权利要求1-42所限定的方法中使用的、包括有参考序列的k-mer的数据库，所述数据库包括：

a)来自参考序列的k-mer的第一集合，和

b)每个k-mer在参考序列中的位置的第二集合。

44.如权利要求43的数据库，其中所述数据库还包括关于与给定参考相关联的全长序列的信息，和/或所述参考的来源，和/或所述参考的一个或多个分类学描述符。

45.如任一权利要求43-44的数据库，其中在数据库中的k-mer进行哈希函数运算，向每个独一无二的k-mer分配唯一密钥。

46.如任一权利要求43-45的数据库，其中在第一集合中的每个唯一k-mer由矢量与关于有k-mer存在的那些参考的信息相关联。

47.如任一权利要求43-46的数据库，其中在第二集合中的每个独一无二的k-mer由矢量与关于该k-mer存在时，其在每个参考中的位置的信息相关联。

48.如任一权利要求43-47的数据库，还包括第三集合或数据库，其具有选自下组的参考标识符和信息类型：描述行、数据来源、序列的任何注释、编码序列、调控序列、可能参考的分类学名称、所述可能参考的近亲、所述参考的来源、一组另外相关的参考、参考是从哪里获得的(例如土壤、海洋、内脏或下水道)、参考序列是何时获得的、分类学分类、近缘物种、关于参考序列是从哪个数据库(例如NCBI或EBI/Sanger数据库)下载的信息。

49.如任一权利要求43-48的数据库，其中k-mer的长度为4、8、12、16、20、24、28、32、36、40、44、48、52、56、60、64或更长。

50.如任一权利要求43-49的数据库，其中k-mer是非重叠的。

51.如任一权利要求43-50的数据库，其中k-mer是重叠的并且增量为至少一个碱基或氨基酸，例如至少两个，如至少3个，例如至少4个，如至少5个，例如至少6个或更多。

52.如任一权利要求43-51的数据库，其中所述数据库包括来自每个参考的完整序列的k-mer。

53.如任一权利要求43-52的数据库，其中所述数据库包括来自人类、动物、哺乳动物、鸟类、鱼类、真菌、昆虫、植物、细菌、古细菌、病毒和/或质粒的序列信息。

54.如任一权利要求43-53的数据库，其中所述数据库被分成被存储在几个不同服务器中的子数据库。

55.如任一权利要求43-54的数据库，其中根据选自门、纲、目、科、属和种的一种或多种分类学描述符，或一种或多种环境描述符，如来源、分布、起源和过去查询的频率，使所述数据库组织成子数据库。

56.一种用于识别来源序列的可能来源的数据处理系统，该系统包括输入装置、中央处理单元、存储器和输出装置，其中，所述数据处理系统存储有表示在执行时使得权利要求1-42的方法得以被执行的系列指令的数据，所述存储器还包括根据任一权利要求43-55的数据库。

57.如权利要求56的系统，其中所述数据库存储在服务器中，并且输入装置和输出装置是客户端，所述客户端和服务器经由数据通信连接来连接。

58.如任一权利要求56-57的系统，其中所述客户端选自个人计算机、台式PC、便携式PC、手持式计算装置如智能电话。

59.如任一权利要求56-58的系统，其中所述客户端包括使客户端能够采样来源序列的子集，将这些来源序列分段成k-mer，并将这些来源序列传送到服务器的一系列指令。

60.如权利要求56-59的系统，所述客户端还包括允许其基于从服务器传送到客户端的序列将来源序列组装成一个或多个更大的序列的一系列指令。

61.如任一权利要求56-60的系统，通过数据连接与测序设备相连接。

62.一种计算机软件产品，含有在执行时使得权利要求1-42的方法得以被执行的系列指令。

63.一种集成电路产品，含有在执行时使得权利要求1-42的方法得以被执行的系列指令。