CN107798216A - 采用分治法进行高相似性序列的比对方法 - Google Patents

采用分治法进行高相似性序列的比对方法 Download PDF

Info

Publication number
CN107798216A
CN107798216A CN201710791282.5A CN201710791282A CN107798216A CN 107798216 A CN107798216 A CN 107798216A CN 201710791282 A CN201710791282 A CN 201710791282A CN 107798216 A CN107798216 A CN 107798216A
Authority
CN
China
Prior art keywords
sequence
comparison
kart
sequencing
comparison area
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.)
Granted
Application number
CN201710791282.5A
Other languages
English (en)
Other versions
CN107798216B (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.)
Academia Sinica
Original Assignee
Academia Sinica
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
Application filed by Academia Sinica filed Critical Academia Sinica
Publication of CN107798216A publication Critical patent/CN107798216A/zh
Application granted granted Critical
Publication of CN107798216B publication Critical patent/CN107798216B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2455—Query execution
    • G06F16/24564—Applying rules; Deductive queries
    • G06F16/24566—Recursive queries
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/28—Databases characterised by their database models, e.g. relational or object models
    • G06F16/284—Relational databases
    • G06F16/285—Clustering or classification
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Medical Informatics (AREA)
  • Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Chemical & Material Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computational Linguistics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提出一种采用分治法进行高相似性序列的比对方法。本发明提出的方法(称为Kart法)是采用分治法将序列切割为数个较小的片段，每个小片段均可个别处理，并且最终序列的全长比对由这些片段的比对组成，因此Kart法可视为能平行进行比对的方法。高通量测序的技术让生物学家得以用精密到核苷酸的分辨率来探讨基因体之间的差异，由于高通量测序可产生巨量的数据，因此高通量测序序列的分析需仰赖快速的比对方法。本发明提出的Kart法可快速地处理短序列与长序列，此外Kart法也可容许更高的测序错误率，根据实验结果，Kart比多数的比对方法要快上许多，即使是错误率高达15％的序列，Kart依然能够产生准确的比对。

Description

采用分治法进行高相似性序列的比对方法
技术领域
本发明涉及一种比对方法，特别涉及一种采用分治法进行高相似性序列的比对方法。
背景技术
高通量测序(Next-generation sequencing,NGS)的技术让生物学家得以用精密 到核苷酸的分辨率来探讨基因体之间的差异，造就许多重大研究的发现。NGS现今已成为 DNA测序以及探讨族群基因体差异的主要方法之一。由于新的测序技术可在一天中产生数 百万，乃至于数十亿以上的核苷酸测序数据，许多NGS应用都需要快速的比对方法来进行大 量序列的分析。传统的序列比对方法，如BLAST[1]或BLAT[2]无法有效率地处理这么庞大的 短序列数据，因此近年来有许多针对NGS短序列的比对方法被发展发出。依据序列建立索引 的方式，这些方法可大致分为两类:哈希表(hash tables)与后缀数组(suffix array)或块 排序压缩(BWT)为主的方法。以哈希表为基础的比对法撷取序列中所有可能的固定长度片 段(k-mer)来取得该片段在数据库中出现的位置信息，而以后缀数组或块排序压缩为基础 的比对法是找寻查询序列与参考序列中的最长一致片段(maximal exact matches,MEM)。 这两类的序列排笔法各有其优缺点，然而以后缀数组或块排序压缩为基础的比对法由于有 较佳的内存配置而较为普及。
以哈希表为基础的比对法包含了CloudBurst[3]、Eland(proprietary)、MAQ[4]、RMAP[5]、SeqMap[6]、SHRiMP[7]、ZOOM[8]、BFAST[9]、NovoAlign(商用软件)、SSAHA[10]以及SOAPv1[11]等等。大多数这类方法都遵循着”种子与延伸”法(seed-and-extendstrategy)[12]，最典型的例子即是BLAST。BLAST记录着数据库序列中所有的固定长度片段的位置信息，使用查询序列的固定长度片段来进行搜寻，从哈希表中找到一致的固定长度片段的纪录。这笔纪录就会被当作种子，并采用Smith-Waterman算法来延伸该种子片段，找找出查询序列与数据库参考序列的相似片段。
采用后缀数组或块排序压缩(BWT)[13]为基础的比对法包含了Bowtie[14,15]、BWA[16]、BWA-SW[17]、BWA-MEM(Heng Li)、SOAPv2[18]、CUSHAW[19]、Subread[20]、HISAT/HISAT2[21]、HPG-aligner[22]以及segemehl[23]。大多数这类的比对法依赖后缀数组来找寻最长的一致片段(称为MEM)，并依此来产生序列比对，而产生的方式也类似于”种子与延伸”法。比较特别的是Subread比对法，它采用的是”种子与投票”(seed-and-vote)法来决定最佳的比对区域。采用后缀数组或块排序压缩的主要优点是重复片段会汇集在一起，因此进行比对时只需要计算一次[12]。
虽然目前已有许多比对法能够处理NGS技术所产生的巨量短序列数据，然而有些方法速度不够快，而有些方法比对不够准确。此外，第三世代的测序技术使得比对更具有挑战性，其测序带来更长的序列以及更高的错误率。例如PacBio RS II系统平均可产生5,500bp至8,500bp长度的测序序列，但是单一测序序列的准确率则平均只有87％，大多数的短序列比对法无法处理这样的测序数据。
发明内容
基于前述种种的困难性，本发明发展了一个新的比对方法，称为Kart法，它同时采用了BWT数组与哈希表来做为索引系统。并且Kart采用分治法来进行序列比对。它将查询序列分切为简易比对区与一般比对区，而每个比对区是各自独立的并且可组合成一个完整的比对。从我们的实验结果中发现，不管原始序列长度为何，需要插入间隔的比对区其平均长度约莫20左右。并且实验结果也显示，Kart产生比对的时间远少于其他方法，而产生的比对准确率也优于或等同于其他方法，即使是处理错误率高达15％的资料仍是如此，这显示Kart能够处理低质量的测序数据。
本发明的主要目的即在提出一种采用分治法进行高相似性序列的比对方法，其是采用分治法将序列切割为数个较小的片段，每个小片段均可个别处理，并且最终序列的全长比对由这些片段的比对组成，因此Kart的算法可视为能平行进行比对的方法。
为达上述目的，本发明的采用分治法进行高相似性序列的比对方法包括以下步骤：提供包含至少一条参考序列的一数据库；以查询序列Q在该数据库中找寻所有区域性的最长一致片段作为简易比对区，并将该些简易比对区根据其序列区块与基因体区块的位置差进行分群，以建构出全局比对的基础架构；以及去除该些简易比对区中的重叠区块，并插入一般比对区以填补相邻简易比对区的间隙，所述一般比对区可以分别并平行地进行比对，所有比对区接合后产生完整的比对。
在本发明的一实施例中，所述简易比对区是通过块排序压缩(BWT)或哈希表(hash tables)搜寻而来，而该块排序压缩与该哈希表是该数据库序列所建立的索引。
在本发明的一实施例中，所述数据库的该参考序列为基因体序列或染色体序列或基因体重组序列，而所述查询序列为基因体序列或染色体序列或基因体重组序列或基因体测序机器所产生的短序列。
在本发明的一实施例中，所述相邻的简易比对区的间隙是两序列片段间的差异点所造成的，该些差异点为两序列所发生的替代错误、插入错误或删除错误所造成的。当相邻的简易比对区的间隙区块为NP-gap free，这部分只需要线性的无间隔比对。如此，大为限缩了耗时的间隔比对区的范围，降低了整体比对的时间。
附图说明
图1为本发明使用找寻长度≥k的所有LMEM的算法，LMEM_Search函数可在BWT数组中切割R[start,stop]，以找出当中可能的最长一致片段，它会回传一个LMEM以及该LMEM在参考序列中出现的位置。
图2为本发明比对方法中的简易比对区A与简易比对区B发生序列重迭的示意图，本发明将较小的简易比对区进行缩减，以消除重叠区域。
图3为本发明中的简易比对区与一般比对区的示意图，本发明利用BWT搜寻算法可将一个测序片段切割为数个子区域，其中简易比对区为完全一致的片段，而一般比对区则含有相异序列的片段。
【符号说明】
A简易比对区
B简易比对区
具体实施方式
下面结合具体实施例来进一步描述本发明，本发明的优点和特点将会随着描述而更为清楚。但这些实施例仅是范例性的，并不对本发明的范围构成任何限制。本领域技术人员应该理解的是，在不偏离本发明的精神和范围下可以对本发明技术方案的细节和形式进行修改或替换，但这些修改和替换均落入本发明的保护范围内。
算法概述
大多数以后缀数组或块排序压缩为基础的比对方法都遵循着”种子与延伸”法，亦即以最长相同片段(MEM)为种子并向左右延伸发展出最终的序列比对，而延伸的方式以动态规划算法来实作。因此种子的搜寻方式与相异处的处理策略对于比对方法的效能有重大的影响，然而这些比对方法本质上都是线性的，本发明所采用的分治法可大大减少间隔比对的长度，进而降低动态规划算法的计算量，这样的算法非常适用于相似性高的序列，如高通量测序的每一条短序列都是基因体中某个特定片段所产生的副本，只是当中可能包含某些测序错误及变异。
简易比对区与一般比对区
由于产生无间隔比对(un-gapped alignment，亦即无插入缺失indels)的速度远快于有间隔比对，因此我们将需要比对的高通量测序短序列与参考序列相对应的区域区分为两类:简易比对区(simple pair)与一般比对区(normal pair)，其中简易比对区具有相同的序列片段，而一般比对区则含有相异点，需要无间隔比对或有间隔比对。找出简易比对区与一般比对区后，这些区域就可以分别并平行地处理，而最终的比对则是这些比对区接合后的结果。
给定短序列R以及参考序列G，我们以G和它的反转序列G’来建构块排序压缩数组(BWT array)，在不失一般性下，我们假设G是G和G’的接合序列。令R[i1]是R的第i1个核苷酸，而R[i1,i2]是介于R[i1]和R[i2]的子序列。同样地，令G[j1]是G的第j1个核苷酸，而G[j1,j2]是介G[j1]和G[j2]的子序列。我们定义一个区域性的最长一致片段(locally maximalexact matches,LMEMs)是R[i1,i2]与G[j1,j2]完全一致的片段，其长度为l，其中R[i1,i2]我们称为查询序列区块，而G[j1,j2]我们称为基因体区块，并且i2–i1＝j2-j1＝l–1，我们以一个四元组(i1,i2,j1,j2)来代表该LMEM，并且以ΔPos＝(j1-i1)代表短序列区块与基因体区块的位置差。
从测序短片段中找寻所有的LMEM
既然一个LMEM代表R与G之间的一段相同序列片段，在本发明中我们将LMEM视为简易比对区，Kart经由搜寻BWT数组来找出所有的LMEM，每一次LMEM搜寻从R[i1]开始并中止于R[i2]，亦即BWT数组的搜寻于R[i2+1]遇到相异点而中止LMEM的延伸。在这种情况下，因R[i2+1]很可能是一个测序错误点或是序列变异点，所以下一次的LMEM搜寻将跳过R[i2+1]，由R[i2+2]开始。只有长度不小于默认值k且出现在参考序列次数低于50的LMEM会被Kart视为合格的简易比对区，k的值域一般介于10和16之间，取决于参考基因体的大小。一般来说较短的LMEM(<k bp)很可能含有测序错误的核苷酸，对应到正确区域的机率值相对来说较低，对于一个大型的基因体来说，需要较大的k值来平衡LMEM的专一性与灵敏性。
图1说明了LMEM的搜寻算法，BWT_search函式是很普通的BWT搜寻法，其输入为查询序列并输出符合资格LMEM以及在参考基因体序列所对应到的位置信息。若查询序列无任何的测序错误点或变异点，那么在BWT数组搜寻后应只输出一个涵盖整个查询序列的LMEM(亦即LMEM.len＝|R|)。然而实际上高通量测序数据中包含许多的错误点与变异点，这些都会导致查询序列被切割为数个不定长度的LMEM，Kart将这些合格的LMEM视为简易比对区，并根据这些简易比对区的分布情况找出一般比对区来产生一个或多个的比对。
判别简易比对区与一般比对区
如果查询序列取自于基因体的重复片段时，该序列就会对应到多个基因体区域，那么简易比对区就有可能分布在多个不同区域，我们定义一个候选比对就是查询序列与参考序列中某个特定区域所对应的一个比对。为了找出所有可能的候选比对，我们根据简易比对区的位置差来进行分群。首先我们排序所有的简易比对区，若相邻的简易比对区的ΔPos差异<g，则我们将这些简易比对区分在同一群中(g的默认值为5)，这些简易比对区就形成一个候选比对。要注意的是，同一个候选比对中，两个简易比对区可能因为串联重复(tandem repeats)、序列变异或重叠的LMEM的关系会有序列重叠的状况。在这种情况下，我们会从较小的简易比对区中删除重叠的片段，以确保在同一个候选比对中的所有简易比对区都没有重叠。图2示范两个重叠的简易比对区，由于序列变异的关系，其中简易比对区A与B发生重叠，因为A较小，所以Kart缩减A的重叠区块，如此一来在同一个候选比对中，任两个简易比对区都不会共享到同一个核苷酸。
接下来我们插入一般比对区来填补两个相邻的简易比对区中的间隙，以产生一个完整的比对。其作法如下。假设两个相邻的简易比对区分别是(i2q-1,i2q,j2q-1,j2q)和(i2q+1,i2q+2,j2q+1,j2q+2)，那么若满足下列条件:i2q+1-i2q>1或j2q+1-j2q>1，那么Kart就会塞入一个一般比对区来填补这中间的间隙，在这种情况下，我们塞入的一般比对区为(ir,ir+1,jr,jr+1)，满足以下条件:ir–i2q＝i2q+1–ir+1＝1if i2q+1-i2q>1，否则令ir＝ir+1＝-1(亦即空核苷酸)。同样地，jr–j2q＝j2q+1–jr+1＝1if j2q+1-j2q>1，否则令jr＝jr+1＝-1。另一方面，若候选比对中的第一个(或最后一个)简易比对区没有涵盖到查询序列的第一个(或最后一个)核苷酸，那么我们也会建立相对应的一般比对区来填补查询序列头尾的空隙。图3示范了简易比对区与一般比对区的概念，在这个范例中，这个查询序列含有三个替代错误点和长度为2的插入点。在LMEM搜寻后，我们可找到四个简易比对区(分别是A,B,C和D)。然而这四个简易比对区并没有涵盖整个查询序列，因此我们检查任两个相邻的简易比对区，依据其间隙产生一般比对区来填补，如此一来，我们在简易比对区A与B之间插入一般比对区(11,11,321,321)，同样地，在简易比对区B与C之间插入一般比对区(23,24,-1,-1)，在C与D之间插入(49,51,357,359)。因此这些简易比对区与一般比对区就可形成一个完整的候选比对。
一般比对区的四种类型
简易比对区来自于完美配对的LMEM，而它们之间的间隙则交错着一般比对区，这些比对区可个别处理并且最终的比对就是这些比对区相接的结果。若进一步来看一般比对区，我们发现一般比对区并不一定需要有间隔比对。当一般比对区的查询序列区块与基因体区块皆超过30个核苷酸，我们就再一次进行序列分割以减少需要有间隔比对的区块大小。这一次我们只在这样的一般比对区中找寻长度超过8个核苷酸的LMEM。由于这样的一般比对区其长度远小于整个基因体，我们采用8-mer索引(在该区域中建置相对应的哈希表)找出查询序列区块与基因体区块共有的8-mer种子，这些种子可接合成较长的LMEM，通过这些LMEM我们可以进一步地将较长的一般比对区切割成数个较小的比对区。
当处理PacBio的测序序列时，若序列两端的一般比对区长度超过5000个核苷酸，则相对应的查询序列区块就直接删除，在这种情况下我们只进行区域性的比对，以避免在低质量测序区域产生不准确的比对。当处理Illumina测序序列时，我们发现若一般比对区的查询序列区块与基因体区块有同样大小时，它们很可能只含有替代错误，那么无间隔比对就是其最佳比对。然而若一个一般比对区含有插入错误，那么无间隔比对就会导致低相似度的比对，因此只要经过线性扫描一遍，我们就可以判定一般比对区是否需要进行有间隔比对。此外，Illumina测序序列也可能含有转接序列(adaptor sequence)或接合序列(chimera)，我们会检视头尾两端的一般比对区，我们将删除序列相似度低于50％的一般比对区。在这种情况下Kart只会产生区域比对。总结上述分析，我们可将一般比对区分成下列四种类型:
1.NP-clip:若(1)PacBio测序序列头尾的比对区长度超过5000个核苷酸或(2)Illumina序测序列头尾的比对区相似度低于50％。
2.NP-gap free:若查询序列区块与基因体区块有相同大小，并且线性扫描后的相异点个数低于区块大小的20％。这部分只需要无间隔比对。
3.NP-indel:若一般比对区的其中一个区块(查询序列区块或基因体区块)为空字符串，另一个区块则含有一个以上的核苷酸。
4.NP-NW:剩下的一般比对区，此类比对区需要由Needleman-Wunsch算法产生有间隔比对。
除了NP-clip类型的一般比对区，可从表格3看出本发明有效地区分出需要进行有间隔比对的NP-NW比对区，在不同测试数据中，这类的一般比对区平均长度只有约20核苷酸，如此我们可大幅地缩短比对时间。
双端测序序列的比对
双端测序序列是两条距离在某个特定范围内的两条测序序列，它们来自于同一个测序片段，可用来帮助判别比对的准确性和可靠性。Kart也支持双端测序序列的比对。要产生双端测序序列的比对，Kart先针对每一条测序序列找出所有可能的候选比对，然后比较这两群候选比对，并找出能够满足双端测序序列距离条件的候选比对。如果没有这样的配对，这表示双端测序序列可能含有较多的测序错误，导致测序序列无法正确地配对到相对应的基因体区域中。在这种情况下，Kart将启动救援程序，根据另一条测序序列所找到的候选比对来产生该条测序序列可能的候选比对，以满足双端测序序列的距离条件。救援程序的具体方式描述如下。
假设G1和G2是双端测序序列R1和R2的两组候选比对，令G1＝{m1,m2,…,mp}且G2＝{n1,n2,…,nq}，其中m1,m2,…,及mp代表R1的候选比对，而n1,n2,…,及nq代表R2的候选比对。对于G1的每一个候选比对m以及其对应的参考基因体坐标c，Kart试着将R2对应到c的下游区域，此时我们建立目标区域的哈希表，并以10-mer来找出可能的LMEM以提升种子搜寻的敏感度。如此一来Kart可能在下游区域搜寻到R2较佳的候选比对n’。Kart重复同样的方式，在G2中的每一个候选比对nj上游中找出R1较佳的候选比对m’。此时mi与n’以及m’与nj就能满足双端测序序列的距离条件，我们从中挑选具有较高序列相似度的比对当作最终的双端测序序列的比对结果。
算法摘要
给定一条测序序列R，Kart先通过BWT搜寻找到所有不定长度的LMEM，每一个LMEM会转换成一个或多个简易比对区。接着Kart根据简易比对区的ΔPos将所有的简易比对区进行分群，去除了相邻简易比对区的重叠区块后，Kart塞入一般比对区以填补简易比对区之间的间隙，产生完整的候选比对。若一般比对区的查询序列区块与基因体区块都有超过30个以上的核苷酸时，Kart就执行第二阶段的序列分割，进一步地将大区块的一般比对区分割成数个较小的比对区，以减少有间隔比对的长度。最终的比对就是这些简易比对区与剩余的一般比对区的接合结果，最后Kart输出具有最高比对分数的候选比对或符合双端测序序列的距离条件的比对。
为证明本发明之技术特征，以下即藉由具体实验及实验结果来验证本发明之功效。
算法实作与实验设计
Kart是在Linux 64-bit环境下以标准C/C++语言发展出来的，它支持了多执行序，Kart读取以BWT为基础的索引档案和一组测序序列数据(单端测序或双端测序)，这些数据以FASTA或FASTQ格式纪录，而Kart以SAM(Sequence Alignment/Map)格式[24]输出比对结果。
由于真实数据缺乏实际答案，难以据此评估序列比对的准确性，因此我们产生仿真数据来评估比对方法的效能。我们利用wgsim(https://github.com/lh3/wgsim)来产生所有的仿真数据，这个仿真程序先以0.1％的突变机率修改原有的基因体序列，其中15％的突变发生INDEL变异，而85％发生SNP变异，接着wgsim再以2％的机率产生测序错误，以仿真Illumina的测序数据。我们同时也用wgsim来仿真PacBio的测序数据，这里的突变机率为13％，并且所有皆为INDEL变异，再加上额外的2％机率模拟测序错误。并且我们可预期未来的测序技术将使测序序列的长度越来越长[17]，例如最新的Illumina MiSeq系统可产生300bp长的测序序列，所以我们各产生一千万条长度为100bp、150bp和300bp的Illumina双端测序序列，以及一百万条长度为7000bp的PacBio单端测序序列。为了便于描述，我们将这些仿真数据命名为Hg19_L100_E02(代表从人体基因组Hg19版本衍生出的仿真数据，长度为100bp，错误率为2％)、Hg19_L150_E02、Hg19_L300_E02以及Hg19_L7000_E15(13％为INDEL错误率，2％为替代错误率)。另外我们也从NCBI SRA数据库与PacBio网站随机下载了四组真实测序资料，分别是SRR622458、SRR826460、SRR826471以及M130929。前三组为Illumina测序数据，最后一组为PacBio测序数据。这些测试资料皆来自于人类基因体样本。
我们采用准确率(precision)与召回率(recall)和运行时间来评估不同方法处理测试数据的效能，一条测序序列若能对应到正确坐标(误差值为30bp)则定义为一个真阳性(TP)，给定一组含有N个序列的测序数据，其中n个测序序列有对应到至少一个基因体区域并产生比对，则准确率与召回率的计算方式如下:
precision＝#of TPs/n×100％；
ecall＝#of TPs/N×100％.
因此，若每一条测序数据皆对应到一个比对，则准确率就会等于召回率。为了避免对应到多个基因体区域而导致效能评估产生偏差，我们只以每一条测序序列的第一条比对来评估效能。对于真实数据，我们则评估其敏感度(sensitivity，亦即对应率)以及运行时间，其中敏感度是有产生比对的测序序列的百分比，定义为
sensitivity＝n/N×100％.
由于缺乏真实数据的实际对应坐标，我们计算每一条比对中的相同核苷酸数来评估比对的准确性，因为最佳比对理论上具有最多相同核苷酸数。
所有的测序数据都在Linux 64-bit环境，配备4颗Intel Xeon E7-48302.13GHzCPU以及2TB内存，Kart与现今普遍使用的比对方法，如BWA-MEM、Bowtie2、Cushaw3、HISAT2、HPG-aligner、Subread、LAST[25]、Minimap[26]和BLASR[27]相互比较，其中BLASR只适用于PacBio数据，而其他方法大多只适用于Illumina数据。其他未考虑的比对方法则因为不支持多线程或不接受测试数据的格式而不在此效能比较中，如Gassst、Ssaha2和NovoAlign。另外有些方法则无法正常运作或花费太多时间也不在此效能比较中，如GEM、hobbes和razers3。
被挑选进效能比较的比对方法皆普遍地使用于NGS数据分析，除非效能不佳，否则我们采用其默认值来分析所有的测试数据，我们也强制所有方法只输出单一最佳比对结果，所有的方法皆使用16个线程来加速比对程序，并且采用其最新的版本来进行比较。
Illumina仿真数据的实验结果
表1显示各比对方法于Illumina仿真资料的评估结果。从表1中我们可看出，大多数的比对方法在各仿真数据可产生相似的比对结果，其中准确率和召回率都介于97-99％。事实上，发生错误比对的因素主要是来自于重复序列区域所导致的模糊性。结果也显示，准确率一般都随着测序序列长度的增加而有所增长，举例来说，Kart在Hg19_L100_E02、Hg19_L150_E02以及Hg19_L300_E02的比对准确度分别为97.8％、98.5％与99.1％，并且BWA-MEM在同测试数据的准确度分别为98.6％、98.9％与99.2％。而Bowtie和HISAT2在较长的测序序列中有较差的敏感度，特别是HISAT2在Hg19_L300_E02的敏感度只有53.6％。
从运行时间来看可以发现，Kart是所有比较方法中执行速度最快的算法，前三个仿真数据的分析结果，Kart的运行时间分别是53、66和113秒。因此我们的分治法提供了NGS测序序列快速比对的解决方案，特别是在处理较长序列速度的优势更加明显。
仿真数据的实验结果
表1同时也揭露了一百万条PacBio仿真数据的评估结果，因为PacBio数据含有较多的插入与删除错误，因此在评估时若比对后的坐标值在实际值的100bp范围内，我们就认为该比对是正确的。此次我们只评估了Kart、BWA-MEM、LAST、Mimimap以及BLASR等方法，而其他的方法并不适用于处理PacBio的长序列上。从表1我们可以看出，Kart、BWA-MEM、LAST和BLASR都产生了相似的召回率，这表示这些方法都可以用来处理高错误率的长序列，然而Kart是当中最快速的方法，其运行时间只有733秒，而BWA-MEM、LAST以及BLASR分别花费4614、78432和9185秒。虽然Minimap只花费了288秒，但它的比对准确率却只有83.4％，并且值得注意的是，Minimap并非是一般的序列比对方法，因为它并不会产生实际的比对，相反地，Minimap只能用来快速地一段找到较长的相似片段。所以我们无法评估Minimap在真实PacBio数据的效能。
表1.Illumina和PacBio的仿真数据测试结果。一千万笔长度分别为100bp、150bp和300bp的双端测序序列以及一百万笔长度为7000bp的单端测序序列，以人类基因体为样本的仿真数据。
真实资料的实验结果
除了仿真数据外，我们从NCBI SRA和PacBio网站分别下载了四组真实数据，分别是
SRR622458(4千万条长度101bp的双端测序序列).
SRR826460(4千万条长度150bp的双端测序序列).
SRR826471(3千4百万条长度250bp的双端测序序列).
M130929(1.2百万长度7118bp的单端测序序列).
表2纪录了各方法在这些真实数据的评估结果。在这个评估中，我们采用敏感度与平均相同核苷酸数来评估比对的质量性。从表2我们可以发现，Kart仍是这些方法中执行速度最快的方法，Kart至少是其他方法数倍快，并且Kart也是这些方法中产生最多相同核苷酸数(敏感度×平均相同核苷酸数)。以SRR622458为例，Kart的敏感度为98.6％，每一个比对平均产生99个相同核苷酸，BWA-MEM、Bowtie2以及Cushaw3可产生与Kart相似的相同核苷酸数，但是这些方法花费却需要花费更多的时间。值得注意的是，由于不明原因，HPG-aligner无法顺利完成SRR622458数据中所有的测序序列的比对。有些方法则遗留较多的序列未完成比对，例如HISAT2在前三个Illumina数据分别只完成86.0％、91.9％和43.9％的测序序列比对。
对于PacBio测序数据M130929来说，Kart和BLASR产生了相近的比对结果，然而BLASR却比Kart花费了更多的时间。虽然BWA-MEM执行速度较BLASR快，但是它的敏感度与平均相同核苷酸数并没有其他方法好。LAST的速度是最慢的，但是它可产生与Kart和BLASR相匹敌的结果。
本发明进一步比较各方法消耗内存的数据。虽然有些比对方法可让使用者设定最大的内存使用量，在测试中我们并没有限制这些方法在内存的使用上限，因此各方法可尽量地使用所有的内存来进行序列比对。在表格2中我们可以看到各方法在不同测试数据中都消耗了类似的内存用量。BWA-MEM、Bowtie2、Cushaw3和HISAT2消耗了较少的内存量(<10GB)，Kart与Subread各需要12GB和18GB，而HPG-aligner和BLASR则消耗了约30GB。
表2.各比对方法于不同长度的测试数据的效能比较
分治法的效率分析
从仿真数据的效能分析中我们可以证明，Kart是一个很有效率的NGS序列比对算法，我们采用了分治法将查询序列切割为简易比对区与一般比对区，并且分别处理各区块的比对。简易比对区本身已是最佳比对，而一般比对区需要花费较多的时间来找到其最佳比对。因此，若一般比对区的比例较低及区块较小的话，那么花费的时间也就会越少。
为了示范在不同长度的序列比对中，Kart所采用的分治法皆能展现独特效率，我们分析了四组仿真数据的简易比对区和一般比对区的平均区块大小。表格3显示了序列分割后的区块平均大小，这些区块分别为LMEM-seed(亦即简易比对区)、8-LMEM-seed(第二阶段简易比对区)、NP-gap free、NP-indels和NP-NW。值得注意的是，前四组区块都不需要进行有间隔的比对，只有最后一组需要。以SRR622458为例，LMEM-seed的平均区块大小为73bp，并且96.5％的核苷酸属于这一群组。当Kart针对较长的一般比对区块进行第二阶段的分割后，我们可找到平均长度为11.4bp的8-LMEM-seed比对区块。最花费计算时间的NP-NW比对区块平均长度为17.5bp并且只有1.9％的核苷酸落在此群组。SRR826460与SRR622458有类似的结果。SRR826471的NP-NW群组有较高的比例，这表示Illumina在长序列(250bp)测序上产生了较多的错误。对于真实的PacBio数据来说，LMEM-seed的平均长度为21.3bp，并且只有13.7％的核苷酸落在此群组，但是第二阶段分割所产生的8-LMEM-seed则有平均12.4bp的长度，并且39.7％的核苷酸落在此群组。值得注意的是，最花费计算时间的NP-NW群组，虽然有44.3％的核苷酸落在此群组，但平均长度只有21.3bp，因此可大幅缩短比对的时间。
表3.四组真实数据的平均简易比对区与一般比对区的长度分析(NP-clip不归在此分析中)。
综上所述，在本发明中提出了Kart，提供NGS序列比对一个高度敏锐、快速且准确的方法。我们利用BWT数组来搜寻出简易比对区，并据此产生出一般比对区，进而重组出最终的比对。每一个简易比对区都代表着查询序列与参考序列之间完全一致的片段，而每一个一般比对区都代表着区块间的差异点。通过实验分析，我们证明了Kart所采用的分治法可大幅减少动态规划的计算量，进而节省了大量的计算时间，特别是在处理较长的测序序列时，效果更加显著。在我们仿真数据与真实数据的效能评估中，Kart可产生最佳或可匹敌的比对，并且是所有方法中花费最少的计算时间。
PacBio测序数据由于有极长的序列和极高的错误率，通常很难找出比对，然而从仿真数据与真实数据的分析中我们发现，Kart不仅可以产生出准确的比对，也比其他方法花费更少的时间。随着测序技术的进步，新的测序机器很可能会产生更长的测序序列，并且包含着更多的错误点。我们的实验结果证明Kart对于各种不同长度和不同质量的测序序列，都可以产生出有效率且准确的序列比对。
在本发明中，我们只示范了分治法在DNA测序序列上的应用，事实上，这个算法对于RNA-seq与基因体序列之间的比对有一样显著的功效，只需要对intron的两端做一些鉴定即可。同时，在更大规模的比对，譬如两个人类基因体序列的比对上，可以将其中一个序列视为查询序列，另一视为参考序列。我们的方法利用平行运算，初步实验，比目前最快方法还要快三百倍。我们方法的特点在于能够快速地找出序列之间相同的片段，并通过这些相同片段定义出其他的无间隔比对片段，藉以分别处理剩余的一般比对区块来达到快速比对的目的。
参考文献：
1.Altschul,S.F.,et al.,Basic local alignment search tool.J Mol Biol,1990.215(3):p.403-10.
2.Kent,W.J.,BLAT--the BLAST-like alignment tool.Genome Res,2002.12(4):p.656-64.
3.Schatz,M.C.,CloudBurst:highly sensitive read mapping withMapReduce.Bioinformatics,2009.25(11):p.1363-9.
4.Li,H.,J.Ruan,and R.Durbin,Mapping short DNA sequencing reads andcalling variants using mapping quality scores.Genome Research,2008.18(11):p.1851-1858.
5.Smith,A.D.,Z.Y.Xuan,and M.Q.Zhang,Using quality scores and longerreads improves accuracy of Solexa read mapping.Bmc Bioinformatics,2008.9.
6.Jiang,H.and W.H.Wong,SeqMap:mapping massive amount ofoligonucleotides to the genome.Bioinformatics,2008.24(20):p.2395-2396.
7.Rumble,S.M.,et al.,SHRiMP:Accurate Mapping of Short Color-spaceReads.Plos Computational Biology,2009.5(5).
8.Lin,H.,et al.,ZOOM！Zillions of oligos mapped.Bioinformatics,2008.24(21):p.2431-2437.
9.Homer,N.,B.Merriman,and S.F.Nelson,BFAST:An Alignment Tool forLarge Scale Genome Resequencing.Plos One,2009.4(11):p.A95-A106.
10.Ning,Z.,A.J.Cox,and J.C.Mullikin,SSAHA:a fast search method forlarge DNA databases.Genome Res,2001.11(10):p.1725-9.
11.Li,R.Q.,et al.,SOAP:short oligonucleotide alignmentprogram.Bioinformatics,2008.24(5):p.713-714.
12.Li,H.and N.Homer,A survey of sequence alignment algorithms fornext-generation sequencing.Brief Bioinform,2010.11(5):p.473-83.
13.Wheeler,M.B.a.D.J.W.a.M.B.a.D.J.,A block-sorting lossless datacompression algorithm.SRC Research Report,1994(124).
14.Langmead,B.,et al.,Ultrafast and memory-efficient alignment ofshort DNA sequences to the human genome.Genome Biol,2009.10(3):p.R25.
15.Langmead,B.and S.L.Salzberg,Fast gapped-read alignment with Bowtie2.Nat Methods,2012.9(4):p.357-9.
16.Li,H.and R.Durbin,Fast and accurate short read alignment withBurrows-Wheeler transform.Bioinformatics,2009.25(14):p.1754-1760.
17.Li,H.and R.Durbin,Fast and accurate long-read alignment withBurrows-Wheeler transform.Bioinformatics,2010.26(5):p.589-95.
18.Li,R.Q.,et al.,SOAP2:an improved ultrafast tool for short readalignment.Bioinformatics,2009.25(15):p.1966-1967.
19.Liu,Y.,B.Schmidt,and D.L.Maskell,CUSHAW:a CUDA compatible shortread aligner to large genomes based on the Burrows-Wheelertransform.Bioinformatics,2012.28(14):p.1830-7.
20.Liao,Y.,G.K.Smyth,and W.Shi,The Subread aligner:fast,accurate andscalable read mapping by seed-and-vote.Nucleic Acids Research,2013.41(10).
21.Kim,D.,B.Landmead,and S.L.Salzberg,HISAT:a fast spliced alignerwith low memory requirements.Nature Methods,2015.12(4):p.357-U121.
22.Tarraga,J.,et al.,Acceleration of short and long DNA read mappingwithout loss of accuracy using suffix array.Bioinformatics,2014.30(23):p.3396-3398.
23.Hoffmann,S.,et al.,Fast Mapping of Short Sequences withMismatches,Insertions and Deletions Using Index Structures.Plos ComputationalBiology,2009.5(9).
24.Li,H.,et al.,The Sequence Alignment/Map format andSAMtools.Bioinformatics,2009.25(16):p.2078-9.
25.Frith,M.C.,R.Wan,and P.Horton,Incorporating sequence quality datainto alignment improves DNA read mapping.Nucleic Acids Research,2010.38(7).
26.Li,H.,Minimap and miniasm:fast mapping and de novo assembly fornoisy long sequences.Bioinformatics,2016.32(14):p.2103-2110.
27.Chaisson,M.J.and G.Tesler,Mapping single molecule sequencing readsusing basic local alignment with successive refinement(BLASR):application andtheory.Bmc Bioinformatics,2012.13.

Claims (5)

1.一种采用分治法进行高相似性序列的比对方法，其特征在于，包括以下步骤：
提供包含至少一条参考序列的一数据库；
以查询序列Q在该数据库中找寻所有区域性的最长一致片段作为简易比对区，并将该些简易比对区根据其序列区块与基因体区块的位置差进行分群，以建构出全长比对的基础架构；以及
去除该些简易比对区中的重叠区块，并插入一般比对区以填补相邻简易比对区的间隙，所述一般比对区可以分别并平行地进行比对，所有比对区接合后产生完整的比对。
2.如权利要求1所述的比对方法，其特征在于，所述简易比对区是通过块排序压缩(BWT)或哈希表(hash tables)搜寻而来，而该块排序压缩与该哈希表是该数据库序列所建立的索引。
3.如权利要求1所述的比对方法，其特征在于，所述数据库的该参考序列为基因体序列或染色体序列或基因体重组序列，而所述查询序列为基因体序列或染色体序列或基因体重组序列或基因体测序机器所产生的短序列。
4.如权利要求1所述的比对方法，其特征在于，所述相邻的简易比对区的间隙是两序列片段间的差异点所造成的，该些差异点为两序列所发生的替代错误、插入错误或删除错误所造成的。
5.如权利要求4所述的比对方法，其特征在于，所述相邻的简易比对区的间隙区块为NP-gap free时，这部分只需要线性的无间隔比对。
CN201710791282.5A 2016-09-07 2017-09-05 采用分治法进行高相似性序列的比对方法 Active CN107798216B (zh)

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
US201662384342P 2016-09-07 2016-09-07
US62/384,342 2016-09-07
US15/694,365 2017-09-01
US15/694,365 US20180067992A1 (en) 2016-09-07 2017-09-01 Divide-and-conquer global alignment algorithm for finding highly similar candidates of a sequence in database

Publications (2)

Publication Number Publication Date
CN107798216A true CN107798216A (zh) 2018-03-13
CN107798216B CN107798216B (zh) 2021-06-04

Family

ID=61281333

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710791282.5A Active CN107798216B (zh) 2016-09-07 2017-09-05 采用分治法进行高相似性序列的比对方法

Country Status (2)

Country Link
US (1) US20180067992A1 (zh)
CN (1) CN107798216B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108776749A (zh) * 2018-06-05 2018-11-09 南京诺禾致源生物科技有限公司 测序数据的处理方法及装置
CN108920902A (zh) * 2018-06-29 2018-11-30 郑州云海信息技术有限公司 一种基因序列处理方法及其相关设备
CN110517728A (zh) * 2019-08-29 2019-11-29 苏州浪潮智能科技有限公司 一种基因序列比对方法及装置
CN110517727A (zh) * 2019-08-23 2019-11-29 苏州浪潮智能科技有限公司 序列比对方法及系统
CN111445952A (zh) * 2020-03-25 2020-07-24 山东大学 超长基因序列的相似性快速比对方法及系统

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114564306B (zh) * 2022-02-28 2024-05-03 桂林电子科技大学 一种基于GPU并行计算的第三代测序RNA-seq比对方法
CN114550820B (zh) * 2022-02-28 2024-05-03 桂林电子科技大学 一种基于WFA算法的第三代测序RNA-seq比对方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102521529A (zh) * 2011-12-09 2012-06-27 北京市计算中心 基于blast的分布式基因序列比对方法
US20120330566A1 (en) * 2010-02-24 2012-12-27 Pacific Biosciences Of California, Inc. Sequence assembly and consensus sequence determination
CN103793628A (zh) * 2012-10-29 2014-05-14 三星Sds株式会社 考虑整个短片段的碱基序列比对系统及方法
CN103793625A (zh) * 2012-10-29 2014-05-14 三星Sds株式会社 碱基序列比对系统及方法
US20140214334A1 (en) * 2013-01-28 2014-07-31 Hasso-Plattner-Institut Fuer Softwaresystemtechnik Gmbh Efficient genomic read alignment in an in-memory database
CN104239749A (zh) * 2013-06-20 2014-12-24 三星Sds株式会社 碱基序列对准系统及方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9679104B2 (en) * 2013-01-17 2017-06-13 Edico Genome, Corp. Bioinformatics systems, apparatuses, and methods executed on an integrated circuit processing platform
US20170270245A1 (en) * 2016-01-11 2017-09-21 Edico Genome, Corp. Bioinformatics systems, apparatuses, and methods for performing secondary and/or tertiary processing

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120330566A1 (en) * 2010-02-24 2012-12-27 Pacific Biosciences Of California, Inc. Sequence assembly and consensus sequence determination
CN102521529A (zh) * 2011-12-09 2012-06-27 北京市计算中心 基于blast的分布式基因序列比对方法
CN103793628A (zh) * 2012-10-29 2014-05-14 三星Sds株式会社 考虑整个短片段的碱基序列比对系统及方法
CN103793625A (zh) * 2012-10-29 2014-05-14 三星Sds株式会社 碱基序列比对系统及方法
US20140214334A1 (en) * 2013-01-28 2014-07-31 Hasso-Plattner-Institut Fuer Softwaresystemtechnik Gmbh Efficient genomic read alignment in an in-memory database
CN104239749A (zh) * 2013-06-20 2014-12-24 三星Sds株式会社 碱基序列对准系统及方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
龚贺华: "LSS-DCA：一个快速的分治多序列对齐算法", 《中国优秀博硕士学位论文全文数据库 (硕士) 信息科技辑》 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108776749A (zh) * 2018-06-05 2018-11-09 南京诺禾致源生物科技有限公司 测序数据的处理方法及装置
CN108776749B (zh) * 2018-06-05 2022-05-03 北京诺禾致源科技股份有限公司 测序数据的处理方法及装置
CN108920902A (zh) * 2018-06-29 2018-11-30 郑州云海信息技术有限公司 一种基因序列处理方法及其相关设备
CN110517727A (zh) * 2019-08-23 2019-11-29 苏州浪潮智能科技有限公司 序列比对方法及系统
CN110517727B (zh) * 2019-08-23 2022-03-08 苏州浪潮智能科技有限公司 序列比对方法及系统
CN110517728A (zh) * 2019-08-29 2019-11-29 苏州浪潮智能科技有限公司 一种基因序列比对方法及装置
CN110517728B (zh) * 2019-08-29 2022-04-29 苏州浪潮智能科技有限公司 一种基因序列比对方法及装置
CN111445952A (zh) * 2020-03-25 2020-07-24 山东大学 超长基因序列的相似性快速比对方法及系统
CN111445952B (zh) * 2020-03-25 2024-01-26 山东大学 超长基因序列的相似性快速比对方法及系统

Also Published As

Publication number Publication date
US20180067992A1 (en) 2018-03-08
CN107798216B (zh) 2021-06-04

Similar Documents

Publication Publication Date Title
CN107798216A (zh) 采用分治法进行高相似性序列的比对方法
Zielezinski et al. Alignment-free sequence comparison: benefits, applications, and tools
US10600217B2 (en) Methods for the graphical representation of genomic sequence data
Chaisson et al. Short read fragment assembly of bacterial genomes
Sundquist et al. Whole-genome sequencing and assembly with high-throughput, short-read technologies
Schbath et al. Mapping reads on a genomic sequence: an algorithmic overview and a practical comparative analysis
Lin et al. Kart: a divide-and-conquer algorithm for NGS read alignment
US20170242958A1 (en) Systems and methods for genotyping with graph reference
Fu et al. MSOAR: a high-throughput ortholog assignment system based on genome rearrangement
Lin et al. AGORA: assembly guided by optical restriction alignment
KR20160073406A (ko) 방향성 비순환 구조에서 쌍형성된-말단 데이터를 사용하기 위한 시스템 및 방법
US9372959B2 (en) Assembly of metagenomic sequences
Alser et al. From molecules to genomic variations: Accelerating genome analysis via intelligent algorithms and architectures
Choi et al. Libra: scalable k-mer–based tool for massive all-vs-all metagenome comparisons
Prezza et al. SNPs detection by eBWT positional clustering
Alser et al. Going from molecules to genomic variations to scientific discovery: Intelligent algorithms and architectures for intelligent genome analysis
Prezza et al. Detecting mutations by ebwt
Saeed et al. A high performance multiple sequence alignment system for pyrosequencing reads from multiple reference genomes
Zerbino Genome assembly and comparison using de Bruijn graphs
Fu et al. A parsimony approach to genome-wide ortholog assignment
Vezzi Next generation sequencing revolution challenges: Search, assemble, and validate genomes
Dai et al. Cloud based short read mapping service
Karimi et al. Binos4dna: Bitmap indexes and nosql for identifying species with dna signatures through metagenomics samples
Sinha et al. A model for optimal assignment of non-uniquely mapped NGS reads in DNA regions of duplications or deletions
Runge et al. RnaBench: A Comprehensive Library for In Silico RNA Modelling

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