WO2018218788A1 - 一种基于全局种子打分优选的三代测序序列比对方法 - Google Patents

一种基于全局种子打分优选的三代测序序列比对方法 Download PDF

Info

Publication number
WO2018218788A1
WO2018218788A1 PCT/CN2017/098122 CN2017098122W WO2018218788A1 WO 2018218788 A1 WO2018218788 A1 WO 2018218788A1 CN 2017098122 W CN2017098122 W CN 2017098122W WO 2018218788 A1 WO2018218788 A1 WO 2018218788A1
Authority
WO
WIPO (PCT)
Prior art keywords
seed
sequence
block
matching
module
Prior art date
Application number
PCT/CN2017/098122
Other languages
English (en)
French (fr)
Inventor
肖传乐
Original Assignee
肖传乐
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 肖传乐 filed Critical 肖传乐
Publication of WO2018218788A1 publication Critical patent/WO2018218788A1/zh

Links

Images

Classifications

    • 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
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression

Definitions

  • the present invention belongs to the field of gene sequencing. Specifically, the present invention relates to a three-generation sequencing (PacBio SMRT and Oxford nanopore sequencing) sequence alignment method, and more particularly to a three-generation sequencing sequence alignment method based on global seed scoring preferred candidate alignment regions.
  • a three-generation sequencing PacBio SMRT and Oxford nanopore sequencing
  • the current three-generation sequencing technology mainly includes PacBio's single molecule (real-time sequencing) technology and the nanopore sequencing technology of the Oxford Nanopore formula.
  • the three-generation sequencing data has the characteristics of long read length (or sequencing sequence) (long read (average 10-15 kb) and no preference for sequencing sequences. These data features can make up for a generation and
  • the second-generation sequencing technology has many flaws, which makes it a widely used market: in genome sequencing, researchers have used three-generation sequencing sequences to complete large-genome assembly, deep analysis of genomic complex regions, 150 gap regions of human genomes, and structural variations.
  • the researchers used the sequencing sequence to contain complete cDNA information to analyze the whole transcriptome alternative splicing and subtype; in DNA modification sequencing, the researchers used the template to modify the base to reduce the polymerase synthesis rate effectively. Detection of unknown modifications of the DNA (eg DNA methylation).
  • the three-generation sequencing technology will become a powerful complement or replacement for the second-generation sequencing technology. In the past two years, it has been widely used in genome assembly, long-segment indel detection and correction, and detection of methylation modification.
  • the high sequencing error rate of the three generations of sequencing data has brought enormous challenges to the processing of three generations of sequencing data.
  • the sequencing data of the three-generation sequencing has the characteristics of high read length (14kbp) and high error rate (up to 15% error rate, which is mainly 10% insertion or 4% deletion, with less 1% substitution), and the second generation.
  • the sequencing and sequencing data has the characteristics of short read length (50-200 bp) and low error rate (error rate is about 1%, mainly due to substitution). Since the data characteristics of the three-generation sequencing and the second-generation sequencing are significantly different, the second-generation sequencing calculation method is obviously not used for the third-generation sequencing data analysis.
  • PacBio has continuously developed the SMART Analysis data analysis platform for the characteristics of three generations of sequencing data.
  • the reference genome alignment (BLASR) and the genome assembly process are very resource-intensive.
  • the 40X human genome three-generation sequencing data, the system software BLASR Completing the reference genome sequence alignment requires 200G memory and tens of thousands of core hours; completing the human genome assembly requires tens of thousands of cores to run for more than three months. That is to say, it is also very challenging for Tianhe No. 2 to complete such assembly calculations. Two or two of them accounted for more than 98% of the total time in the assembly process.
  • the reference genome alignment and the pairwise alignment process require a large amount of computation, which constrains the wide application and development of three generation sequencing.
  • sequence alignment method an efficient two-generation sequencing pairwise alignment method and reference genome alignment method (sequence alignment method) has high practical application value.
  • the speed of the preferred reference genome alignment method based on the global seed score is 5-100 times higher than the current three-generation sequencing reference genome comparison software BLASR and BWA-mem; based on the global seed score on the 54X human genome (preferably pairwise alignment method)
  • the speed is currently three generations of sequencing pairwise software MHAP and Daligner is 20-100 times.
  • the system and method of the present invention can greatly reduce the computation time and resources required for the current three generations of sequencing, and has good commercial value.
  • the present invention provides a three-generation sequencing sequence alignment system based on global seed scoring, which comprises module 1, module 2, module 3, module 4 and module 5, and module 1 is configured to quickly find significant candidate overlaps.
  • module 5 is fitted with a global seed scoring model based on the block data structure, where module 5 contains module 5.1, module 5.2 and module 5.3, and module 5.1 fits the distance difference factor between the two sequences between the two seed pairs.
  • Module 5.2 is fitted with two kinds of sub-voting scoring to obtain the core matching seed position pair rule, and module 5.3 chimeric extended voting scoring to obtain the global seed voting score rule of the core position pair.
  • the above system further comprises a module 6 which is fitted with preferences and usage rules based on global seed scoring.
  • the module 4 fits the block seed matching number and the sensitivity mathematical model, the matching block seed matching number and the sensitivity mathematical model include the reference genome block seed matching number and the sensitivity mathematical model and the pairwise alignment. Block seed matching number and sensitivity mathematical model.
  • Step 2.5 Map all seed sequences into the block data structure of the Z-fold reference sequence
  • the block data structure method for constructing the link reference sequence described in step 2.3 is:
  • the core seed position covered by the sequencing sequence is taken out to the seed position pair in the adjacent block structure, and the core seed pair source can be located through the starting position of each sequencing sequence and the number of the significant block matching on the 2G link reference sequence.
  • the number of the sequencing sequence in 2G is obtained according to the overlapping condition of the two sequencing sequences, the range of matching of the adjacent blocks is obtained, and the core position of the overlapping area is scored by one-way voting, and the core position is obtained for the global seed score.
  • the number of seeds of the 80% seed position of the adjacent block structure to the block structure supporting the core seed position pair will be set to 0.
  • the matching block is a matching sequence after the seed of the sequencing sequence is mapped to the reference genomic block data structure.
  • the block data structure is a matching block of the sequencing sequence, and Said to be a significant match block.
  • the global seed score of each candidate alignment region represents the overlap length of the candidate region, and the global candidate can effectively optimize the candidate region with a long overlapping region, thereby greatly reducing Entering the candidate region of two-two local alignment; based on the global seed scoring model, three generations of sequencing pairwise alignment method and reference genome alignment method are designed, which greatly accelerates the three-generation sequencing sequence alignment process and calculation. Resource consumption.
  • the rules of the system and system thereof of the present invention enable the method of the present invention to greatly reduce the computation time and resources required for the current three generations of sequencing, and have good commercial value.
  • Figure 1 Schematic diagram of block data structure model in module 1.
  • FIG. 1 Schematic diagram of the reference sequence index in Module 2
  • FIG. 1 Schematic diagram of the seed sequence sampling rule in Module 2
  • Example 1 Reference genomic alignment method based on global seed scoring
  • i corresponds to the position of the base in the sequence
  • NC i is the number corresponding to the corresponding position letter, and converts it into decimal data.
  • Step 4 Align all seed sequences into a Z-fold reference genomic block data structure
  • the block structure seed counter When the seed of a sequencing sequence is aligned to the region of the CR block data structure, the block structure seed counter will be incremented by 1, and the seed matching position pair of the structure will record the position of the seed in the sequencing sequence and the relative position in the CR block region. position.
  • the reference genome candidate position (SL) of each sequencing sequence seed is mapped into the reference genome block data structure according to the above rules and formulas, and all matching seed block data structure numbers (CR) are recorded using a look-up table.
  • Step 5 Select the starting seed position pair (core seed position pair) of the local sequence alignment from the block data structure:
  • Step 6 Get the global voting score for the starting seed location pair
  • the number VL and VR of the block data structures that the sequencing sequence can span across the left and right sides is estimated according to Equation 5 and Equation 6. All seed position pairs of the sequencing sequence covering the adjacent block data structure will vote against the starting seed position pair in accordance with the DF formula to obtain a global voting score for the starting position pair.
  • the block structure seed number will be set to 0 and is no longer considered as a candidate block structure.
  • Step 7 According to the global score, select the top 10 start position pairs for local two-two sequence comparison
  • the global voting scores of the starting seed position pairs of each block data structure higher than the seed threshold are obtained in descending order for the above steps 5 and 6. Select the global voting score of up to 20 starting position pairs to complete the local two-two sequence alignment by the modified diff algorithm.
  • the starting position corresponds to two conditions for the sequence alignment result: when a start position pair meeting the overlapping long reading >1000 and the false matching rate ⁇ 0.20 is encountered, the sequencing sequence alignment process is terminated, and the result is used as a sequence of the sequencing sequence. Compare the results output.
  • Step 8 Clean up the sequencing sequence
  • the seed counters of all block structures matched by the sequencing sequence seed are reset to 0, and the look-up table records are emptied. Repeat steps 3 through 8 for reading the next sequencing sequence. The reference genome alignment was completed until all sequencing data was completed.
  • Step 9 Secondary Accurate Search Sequence Alignment Analysis
  • steps 2 to 8 are established by the multithreaded package based on the shared memory variable space pthread, wherein the reference genome index of step 1 will be placed in the multicore shared memory.
  • Example 2 Pairwise alignment method based on global seed score preference
  • the pairwise alignment method based on the global seed score preference is basically similar to the reference genome implementation process of Example 1, with the following differences:
  • Step 1 Data Blocking and 2G Link Sequence Acquisition: Scan the entire three generations of sequencing data files, segment the three generations of sequencing data according to the 2G file size, link the two sequencing sequences with N, and record each sequencing sequence in the 2G link reference sequence. The starting position and termination are performed, and the 2G file is indexed for each sequencing sequence file location, which facilitates subsequent calculation of the position on the linked reference sequence to be translated into the absolute position of each sequencing sequence.
  • Step 2 Same as step 1 of the embodiment 1.
  • Step 3 Similar to step 2 of Embodiment 1, except that the block structure size Z is changed to 2000.
  • Step 4 Similar to step 3 of Example 1, except that the seed (k-mer) step ST is changed to 10.
  • Step 5 Same as step 4 of the embodiment 1.
  • Step 6 Similar to step 5 of Example 1, except that the two seed positions are changed to support conditions to DF ij ⁇ 0.3.
  • the position of the core position on the upper reference genome is converted to the number of the sequencing sequence (read) where the position is located and the sequence is absolutely determined by the start and end positions of each sequencing link sequence on each sequencing sequence position index of the 2G file. position.
  • Step 7 Similar to step 6 of Embodiment 1, it is necessary to modify the two seed positions to support the condition to DF ij ⁇ 0.3.
  • the overlapping region range of two sequencing sequences is obtained by: according to the block number of the significant matching pair and the starting position of the sequencing sequence of the linked reference sequence, the matching block positioning sequencing sequence number and the starting position (S 1 , E 1 ) can be obtained.
  • the position information of the core position pair (reference genomic position P 1 , the position of the sequencing sequence is P 2 ), and the length of the sequencing sequence to be compared is L, it can be concluded that the core position of the matching sequence on the linked reference sequence is the length L l on the left side .
  • the length of the left side of the sequence to be aligned is P+
  • the length of the right side is L-P+
  • the length of the two left sides is the length of the left overlapping area
  • the shorter length on the right side is the length on the right side
  • the two length ranges are the range of the extension block structure.
  • Step 8 The output of the global scoring of the two pairs: in the pairwise comparison, only the core seed position information of the highest 100 global seed scores in 2G is obtained, and no local sequence alignment is needed, and the highest 100 core seed position information is obtained. It is converted to absolute position information that will be converted into two sequencing sequences, and finally the number of the two overlapping sequencing sequences, the absolute position information of the core position pair, and the global voting score are output.
  • Step 9 is similar to step 8 of Example 1, except that after reading the next sequencing sequence, steps 4 through 9 are performed.
  • Step 10 Program Parallelization: Steps 4 through 9 above are used to build a parallelization program based on the shared memory variable space pthread multithreading package, where the reference genome index of step 2 will be placed in multicore shared memory.
  • Step 11 Pairwise alignment of each data block: Data block 1 is to be compared with data block 1-n, data block 2 is to be compared with data block 2-n, and then all sequences are analogized. Pairwise alignment, two pairs of two pairs are matched to match the same two sequences. In the comparison process, the subsequent global seed voting analysis is performed after the sequencing sequence number of the sequencing sequence is larger than that of the significant matching block.
  • the pairwise comparison software compares randomly extracted 500M data.
  • the speed of our software MECAT pairwise comparison is 2-8 times that of MHAP and Daligner software; in the nanopore dataset, the MECAT speed is MHAP and Daligner. 5-10 times.
  • the reference genome software comparison uses the entire data set for comparison.
  • the speed of our software MECAT pairwise alignment is 5-70 times that of BLASR and BWA software; in the nanopore dataset, MECAT speed BLASR and BWA 4-5 Times.
  • the above table time unit is nuclear time

Abstract

一种基于全局种子打分优选的三代测序序列比对方法，该方法通过基于全局种子打分优选的三代测序序列比对系统实现，系统包含模块1、模块2、模块3、模块4和模块5。所述系统及其系统中的规则和应用方法，大幅降低目前三代测序需要的计算时间和资源，具有良好商业价值。

Description

一种基于全局种子打分优选的三代测序序列比对方法 技术领域
本发明属于基因测序领域，具体的，本发明涉及三代测序(PacBio SMRT和Oxford nanopore测序)序列比对方法，特别是涉及一种基于全局种子打分优选候选比对区域的三代测序序列比对方法。
背景技术
目前三代测序技术主要包含PacBio公司的单分子实时测序(single molecule,real-time，SMRT)测序技术和Oxford Nanopore公式的纳米孔(Nanopore)测序技术。与二代测序技术相比，三代测序数据具有读长(或测序序列)很长(long read，平均10-15kb左右)和测序序列无GC偏好性等特点，这些数据特征可以有力弥补了一代和二代测序技术很多缺陷，从而使其具有广泛应用市场：在基因组测序方面，研究者利用三代测序的测序序列完成了大基因组组装、基因组复杂区深度解析、人类基因组150个gap区域和结构变异的解析；在转录组测序方面，研究者利用测序序列已包含完整cDNA信息深入分析全转录组可变剪接和亚型；在DNA修饰测序方面，研究者利用模板修饰碱基降低聚合酶合成速率来有效检测DNA未知的修饰(例如DNA甲基化)。目前，三代测序技术将成为二代测序技术的有力补充或替代，近两年来广泛应用于基因组组装、长片段indel检测和矫正、以及甲基化修饰的检测等研究中。
三代测序数据高测序错误率给三代测序数据处理带来了巨大的挑战。三代测序的测序数据具有测序序列高读长(14kbp)和错误率高(错误率高达15％，其主要是10％插入或4％缺失，有较少的1％替代)等特点，而二代测序测序数据具有短读长(50-200bp)和错误率低(错误率约1％，主要是替代产生)等特点。由于三代测序和二代测序的数据特征有显著不同，因此二代测序计算方法显然不实用于三代测序数据分析。目前PacBio公司针对三代测序数据特征不断开发SMART Analysis数据分析平台，然而该系统参考基因组比对(BLASR)和基因组组装流程两两比对计算十分耗资源：40X人类基因组三代测序数据，该系统软件BLASR完成参考基因组序列比对需要200G内存和几万核小时；完成人类基因组组装，需要几万个核运行三个月以上才能完成，也就是说，天河二号完成这样组装计算量也很有挑战，其中两两比对在拼装流程中占了总时间的98％以上。参考基因组比对和两两比对过程需要大量计算量，约束了三代测序广泛应用和发展。
因此，创建一种高效三代测序的两两比对方法和参考基因组比对方法(序列比对方法)具有很高实际应用价值。基于全局种子打分优选参考基因组比对方法的速度是目前三代测序参考基因组比对软件BLASR和BWA-mem的5-100倍；在54X人基因组上，基于全局种子打分(优选两两比对方法)的速度是目前三代测序两两比对软件MHAP和 Daligner的20-100倍。本发明的系统和方法可以大幅降低目前三代测序需要的计算时间和资源，具有良好商业价值。
发明内容
为解决上述技术问题，本发明提供了基于全局种子打分优选的三代测序序列比对系统，该系统包含模块1、模块2、模块3、模块4和模块5，模块1嵌合快速查找显著候选重叠区域的block数据结构模型，模块2嵌合参考基因组block数据结构的映射规则，模块3嵌合参考基因组索引和测序序列(read)种子序列抽样规则，模块4嵌合匹配块(block)种子匹配数与灵敏度数学模型，模块5嵌合基于块数据结构的全局种子打分模型，其中模块5包含模块5.1、模块5.2和模块5.3，模块5.1嵌合两个种子对之间的两个序列的距离差异因子，模块5.2嵌合两两种子投票打分获取核心匹配种子位置对规则，模块5.3嵌合延伸投票打分获取核心位置对的全局种子投票得分规则。
优选的，上述系统还包含模块6，模块6嵌合基于全局种子打分的优选和使用规则。
上述系统中，模块1嵌合快速查找显著候选重叠区域的block数据结构模型，所述快速查找显著候选重叠区域的block数据结构模型为：设Z为block数据结构的块比例，即块大小，对于参考基因组每Z个碱基建立一个块(block)数据结构，并顺序编号，用于比对过程中快速将种子序列定位到候选比对区域。每个块数据结构中包含种子计数器、p个种子位置对组成：种子位置对记录某一种子在测序序列的位置和该种子在参考基因组的匹配位置；种子计数器用来记录候选块结构比对的种子个数，同时表示块结构热点区域的得分。通常情况，某个特定块结构种子数得分越高，表示测序序列落在此块区间可能性越大(参见附图1)。
上述系统中，模块2嵌合参考基因组block数据结构的映射规则，参考基因组block数据结构的映射规则为：
通过测序序列(read)每个种子(k-mer)编码查询参考基因组索引获得每个种子基因组的精确位置，并用每个种子的精确位置按照公式1的规则映射到上述块结构：
Figure PCTCN2017098122-appb-000001
其中Z表示块结构碱基区域大小，CR表示块结构的序号，CL为种子在参考基因组的在块结构的相对准确位置，SLi表示参考基因组候选位置。
其中，当种子序列比对到第CR号块结构中时，该块结构的种子计数器个数加1，并建立查询表记录所有可能候选区域block结构的位置和该区域种子数。
上述系统中，模块3嵌合参考基因组索引和测序序列(read)种子序列抽样规则，参考基因组索引和测序序列(read)种子序列抽样规则为：
以参考基因组每个位点为起始，取k＝13个碱基长度的片段作为种子序列(k-mer)，建立种子(k-mer)的4进制编码与其对应起始位置的哈希表(参见附图2)。
哈希表中记录每个种子(k-mer)的编码和该种子(k-mer)在参考基因组上所有位置，即参考基因组索引，通过种子序列可以查找基因组中相同序列片段的所有位置；测序序列(read)，每隔特定步长(ST)取k个碱基长度的种子序列，并顺序标号记录，用来寻找序列和参考序列中完全匹配的种子信息(参见附图3)。
上述系统中，模块4嵌合匹配块(block)种子匹配数与灵敏度数学模型，匹配块(block)种子匹配数与灵敏度数学模型包括参考基因组块种子匹配数与灵敏度数学模型和两两比对中块种子匹配数与灵敏度数学模型。
其中，参考基因组块种子匹配数与灵敏度数学模型如下：
假设所有种子(k-mer)比对是独立事件，种子(k-mer)的匹配概率初步符合二项式分布，在参考基因组比对过程中，种子(k-mer)匹配概率用如下公式2计算：
Povl＝(1-e)k  (公式2)
公式2中，当e为0.15，块大小(Z)为1000，种子(k-mer)抽样步长(ST)为20和k为13时，每个块的抽样数为
Figure PCTCN2017098122-appb-000002
两个匹配块平均种子(k-mer)匹配个数为
Figure PCTCN2017098122-appb-000003
现将块匹配阈值设为6，由累计概率公式可知，两个匹配块小于6个种子(k-mer)匹配的概率为26.67％.，假设读长平均重叠长度14kbp(三个重叠块)，则参考基因组比对灵敏度为99.99％。
其中，两两比对中块种子匹配数与灵敏度数学模型如下：
假设所有种子(k-mer)比对是独立事件，种子(k-mer)的匹配概率初步符合二项式分布，在两两比对过程中，种子(k-mer)匹配概率通过如下公式3计算：
Figure PCTCN2017098122-appb-000004
公式3中，当e为0.15，块大小(Z)为2000，种子(k-mer)抽样步长(ST)为5和k为13时，每个块的抽样数为
Figure PCTCN2017098122-appb-000005
两个匹配块平均种子(k-mer)匹配个数为
Figure PCTCN2017098122-appb-000006
现将块匹配阈值设为5，由累计概率公式可知，两个匹配块小于5个种子(k-mer)匹配的概率为0.2，假设两个读长平均重叠长度>6000(三个重叠块)，则两两比对灵敏度为99.2％。
上述系统中，模块5嵌合基于块数据结构的全局种子打分模型，基于块数据结构的全局种子打分模型如下：
对参考基因组(reference)和三代测序序列(read)分别建立种子(k-mer)(k＝13)的哈希表(参见附图2)，同时将基因组和序列分成大小为1000bp的数据块。如果基因 组和测序序列的两个块共享的种子(k-mer)大于阈值(6)时，这两个块就称为一个显著匹配块。
上述系统中，所述模块5中的全局种子打分从显著匹配块开始,其过程包含如下模块5.1，模块5.2和模块5.3，其中模块5.1嵌合两个种子对之间的两个序列的距离差异因子，模块5.2嵌合两两种子投票打分获取核心匹配种子位置对规则，模块5.3嵌合延伸投票打分获取核心位置对的全局种子投票得分规则。
优选的，模块5.1嵌合两个种子对之间的两个序列的距离差异因子，两个种子对之间的两个序列的距离差异因子的计算方法如下：
为有效过滤假阳性的匹配块，引入序列差异因子(DFF)：对位于(s1，t1)位置对的种子(k-mer)匹配和位于(s2，t2)位置对的种子(k-mer)匹配，(s1和s2是参考基因组位置，t1和t2是三代测序序列(read)位置)，两个种子位置对之间序列最短编辑距离是两个种子位置对之间序列长度之差，定义两个匹配种子对之间的长度之差序列差异因子用公式4计算：
Figure PCTCN2017098122-appb-000007
若DFF≤e，则为两个种子位置对相互支持，位置对各加一分。
优选的，模块5.2嵌合两两种子投票打分获取核心匹配种子位置对(起始匹配种子位置对)规则，两两种子投票打分获取核心匹配种子位置对规则为：
在显著匹配块结构中，通过种子位置对的两两投票打分获取核心位置对：当两个种子位置对之间符合DDF＜0.3时，两个种子位置对各自加一分，所有位置对获得块匹配内另外种子的两两投票打分，投票得分最高的位置对即为核心位置对，当一个多个位置对获得相同分数时，选择第一个位置对为核心位置对。
优选的，模块5.3嵌合延伸投票打分获取核心位置对的全局种子投票得分规则，延伸投票打分获取核心位置对的全局种子投票得分规则如下：
测序序列长度(平均长度14kbp)远远大于Z bp，因此测序序列的种子序列通常覆盖多个相邻的块结构区域，核心种子对相邻块数据结构将对核心种子对进行单向种子投票打分。
根据开始种子位置对(SLk，SNk)估算测序序列可以跨越左边和右边相邻的块结构的范围：测序序列在开始种子位置对左右和右边的长度分别是SNk和LL-SNk，可以通过下列公式5和公式6计算测序序列覆盖左和右相邻块数据结构数据的范围(VL和VR)：
Figure PCTCN2017098122-appb-000008
Figure PCTCN2017098122-appb-000009
在公式中，LL是测序序列的长度。测序序列覆盖相邻块结构的所有种子位置对将按照公式4对核心位置对进行种子对投票打分，从而获得核心位置对的全局投票打分(10D)，而获得核心位置对的全局分数，即该候选区域全局得分。(参见附图4)。
当一个相邻块结构中的80％种子位置对符合DFF≤e，即支持开始种子位置对，该块结构种子数将被设置为0，并且该块结构编号将在块结构查阅表中删除，也就是，该块结构不再被考虑为候选块结构。
优选的，上述系统中还包含模块6，所述模块6嵌合基于全局种子打分的优选和使用规则，基于全局种子打分的优选和使用规则如下：
对所有显著匹配块候选区域进行全局种子投票打分获取每个显著匹配块的核心种子对和全局种子得分，根据每个候选区域的核心种子对和全局种子得分判定核心种子位置对进入后续局部序列比对分析，其中全局种子投票打分方法为基于全局种子打分优选的参考基因组比对方法或基于全局种子打分优选的两两比对方法，其中判定方法如下：
(1)应用基于全局种子打分优选的参考基因组比对方法获得的结果，当参考基因组选择最高10个核心种子位置对所在的区域有效候选区域，这些核心种子位置对可以进入后续局部序列比对分析；
(2)应用基于全局种子打分优选的两两比对方法获得结果，选择最高100个核心种子位置对所在的区域有效候选区域，这些核心种子位置对可以进入后续局部序列比对分析。
本发明还提供了一种基于全局种子打分优选的三代测序序列比对方法，所述三代测序序列比对方法为基于全局种子打分优选的参考基因组比对方法和基于全局种子打分优选的两两比对方法中的一种或两种，所述基于全局种子打分优选的参考基因组比对方法和基于全局种子打分优选的两两比对方法执行模块1、模块2、模块3、模块4和模块5中的至少2个以上的模块。
优选的，所述方法还包括执行模块6。
优选的，本发明提供了一种基于全局种子打分优选的三代测序序列比对方法，所述方法为基于全局种子打分优选的参考基因组比对方法，所述参考基因组比对方法包括如下步骤：
步骤1.1：建立参考基因组索引
步骤1.2：构建参考基因组块数据结构
步骤1.3：分割测序序列序列成若干个种子序列
步骤1.4：将所有种子序列映射到Z倍参考基因组块数据结构中
步骤1.5：获取显著块匹配区域的核心种子位置对
步骤1.6：获取核心种子位置对的全局投票打分
步骤1.7：选择最高n个核心位置对进行局部两两序列比对
步骤1.8：二次精准参考基因组序列比对。
上述参考基因组比对方法中，优选的，步骤1.1所述的建立参考基因组索引方法为：
应用模块3，从参考基因组中每个碱基位置提取k(k-mer)长度种子序列，也就是，相邻的种子(k-mer)之间没有间隔。参考基因组所有的碱基将被建立种子(k-mer)索引。
上述参考基因组比对方法中，优选的，步骤1.2所述的构建参考基因组块数据结构方法为：
应用模块1，将参考基因组每Z个碱基区域建立一个块数据结构，每个块数据结构用于记录测序序列种子在该结构代表参考基因组区域的匹配情况。优选的，每个块数据结构由种子匹配计算器、40个种子匹配候选种子位置对组成。
上述参考基因组比对方法中，优选的，步骤1.3所述的分割测序序列序列成若干个种子序列方法为：
应用模块3，在测序序列中按照ST＝20步长提取种子(k-mer)的种子序列，每个种子有k个碱基组装，并按照测序序列顺序进行编码(SN)。
上述参考基因组比对方法中，优选的，步骤1.4所述的将所有种子序列映射到Z倍参考基因组块数据结构中的方法为：
一个测序序列种子序列的所有参考基因组候选位置(SLi,i＝1,2,...n)可以从步骤1.1的参考基因组索引中查找，应用模块2，将每个种子所有候选位置映射到参考基因组块数据结构中存储。并且用查阅表(look-up table)记录所有匹配种子块数据结构编号(CR)，查阅表记录着种子匹配的块区域编号和对应块区域的种子匹配数，每个块区域在查阅表中唯一记录.
上述参考基因组比对方法中，优选的，步骤1.5所述的获取显著块匹配区域的核心种子位置对方法为：
当一个块数据结构的种子数大于7时，该块结构被认为显著块匹配结构，显著块匹配结构的局部比对的核心匹配位置对将通过该块结构中所有种子对两两投票打分确定，按照模块5.2进行两两投票打分获取该显著匹配块的核心种子位置对。
上述参考基因组比对方法中，优选的，步骤1.6所述的获取核心种子位置对的全局投票打分方法为：
应用模块5.3，将测序序列所覆盖的核心种子位置对相邻块结构中位置对取出，对核心位置对进行单向投票打分，获取核心位置对全局种子得分。并将相邻块结构的80％ 种子位置对支持核心种子位置对的块结构的种子数将被设置为0.
上述参考基因组比对方法中，优选的，步骤1.7所述的选择最高n个核心位置对进行局部两两序列比对方法为：
通过步骤1.5和步骤1.6，获得每个高于种子阈值的块数据结构的核心种子位置对和全局投票打分。之后，对所有核心位置对的全局投票得分进行降序排序，选择全局投票得分最高10个核心位置对通过diff方法完成局部两两序列比对，对用nanopore，采用smith-waterman方法进行局部两两比对.如果核心位置对序列比对结果符合两个条件：重叠长读>1000和错误匹配率<0.20，认为该测序序列已找到正确参考基因组匹配位置。按照全局比对得分顺序进行两两序列比对，当遇到符合上述条件的核心位置对时，终止该测序序列序列比对过程，将该结果作为测序序列的序列比对结果输出。
上述参考基因组比对方法中，优选的，步骤1.8所述的二次精准参考基因组序列比对方法为：
针对少数测序序列的块匹配种子量较少，而且布局均一，不能被上述步骤1.4参数搜索到。如果上述过程ST步长分割和Z数据结构没有获得搜索结果输出，将执行步骤1.3的ST变成ST/2步长(10)，之后的块大小为2Z(2000)，其它参数不变，重复上述步骤3到步骤1.7进行更精确的序列比对过程。
优选的，本发明提供了一种基于全局种子打分优选的三代测序序列比对方法，所述方法为基于全局种子打分优选的两两比对方法，所述两两比对方法包括如下步骤：
步骤2.1：三代测序数据分块和测序序列链接成类似参考基因组
步骤2.2：建立参考基因组索引
步骤2.3：构建链接参考序列的块数据结构
步骤2.4：分割测序序列成若干个种子序列
步骤2.5：将所有种子序列映射到Z倍链接参考序列的块数据结构中
步骤2.6：获取显著块匹配区域的的核心种子位置对
步骤2.7：获取核心种子位置对的全局投票打分
步骤2.8：选择最高n个核心位置对的候选区域输出结果。
上述两两比对方法中，优选的，步骤2.1所述的三代测序数据分块和测序序列链接成类似参考基因组方法为：
将三代测序数据集分成2G大小数据块，链接2G数据块内的测序序列(read)成2G的一条参考序列，两条测序序列链接出添加一个N字母，记录每个测序序列在2G参考序列上的起始位置，方便后续寻找两个测序序列重叠的起始位置。
上述两两比对方法中，优选的，步骤2.2所述的建立参考基因组索引方法为：
应用模块3，从链接后的2G参考序列中每个碱基位置提取k(k-mer)长度种子序列，也就是，相邻的种子(k-mer)之间没有间隔。参考基因组所有的碱基将被建立种子(k-mer)索引。
上述两两比对方法中，优选的，步骤2.3所述的构建链接参考序列的块数据结构方法为：
应用模块1，链接参考序列每Z(Z＝2000)个碱基区域建立一个块数据结构，每个块数据结构用于记录测序序列种子在该结构代表建链接参考序列区域的匹配情况。每个块数据结构由种子匹配计算器、40个种子匹配候选种子位置对组成。
上述两两比对方法中，优选的，步骤2.4所述的分割测序序列序列成若干个种子序列方法为：
应用模块3，在测序序列中按照ST＝10步长提取种子(k-mer)的种子序列，每个种子有k个碱基组装，并按照测序序列顺序进行编码(SN)。
上述两两比对方法中，优选的，步骤2.5所述的将所有种子序列映射到Z倍链接参考序列的块数据结构中的方法为：
一个测序序列种子序列的所有参考基因组候选位置(SLi，i＝1,2,…,n)从步骤2.1的参考基因组索引中查找，应用模块2，将每个种子所有候选位置映射到链接参考序列块数据结构中存储。并且用种子映射查阅表(look-up table)记录所有测序种子映射到的块数据结构(block)编号(CR)。种子映射查阅表由两个数据构成：1)测序序列种子映射到的块数据结构(block)编号，每个块数据结构编号在查阅表中唯一记录；2)该种子映射的块数据结构的种子匹配数。
上述两两比对方法中，优选的，步骤2.6所述的获取显著块匹配区域的的核心种子位置对方法为：
当一个块数据结构的种子数大于7时，该块结构被认为显著块匹配结构，显著块匹配结构的局部比对的核心匹配位置对将通过该块结构中所有种子对两两投票打分确定，应用模块5.2进行两两投票打分获取该显著匹配块的核心种子位置对。
上述两两比对方法中，优选的，步骤2.7所述的获取核心种子位置对的全局投票打分方法为：
应用模块5.3，将测序序列所覆盖的核心种子位置对相邻块结构中种子位置对取出，通过2G链接参考序列上每个测序序列的起始位置和显著块匹配的编号可以定位核心种子对来源2G中的测序序列的编号，根据两个测序序列重叠情况，获取相邻块匹配的范围，对重叠区域的核心位置对进行单向投票打分，获取核心位置对全局种子得分。并将相邻块结构的80％种子位置对支持核心种子位置对的块结构的种子数将被设置为0.
上述两两比对方法中，优选的，步骤2.8所述的选择最高n个核心位置对的候选区域输出结果的方法为：
通过步骤2.5和步骤2.6，获得每个高于种子阈值的块数据结构的核心种子位置对和全局投票打分。然后，对所有核心位置对的全局投票得分进行降序排序，选择全局投票得分最高100个核心位置对，将每个核心位置对信息通过测序序列在2G链接参考序列起始位置转化成两条测序序列的绝对位置信息，最后输出两个重叠测序序列的编号、核心位置对的绝对位置信息和全局投票得分。
上述方法中，其中种子(种子序列，seed sequence)是测序序列中的k长度的子序列，种子序列与参考基因组匹配需要长度相同且无空位的完全匹配序列，以种子作为参考，寻找测序序列与参考基因组序列的匹配时分值超过一定阈值的相似性片段。
其中块(block)数据结构(也称块结构)是指对于参考基因组每Z个碱基建立一个块(block)数据结构，并顺序编号，用于比对过程中快速将种子序列定位到候选比对区域。每个块数据结构中包含种子计数器、p个种子位置对组成(参见附图1)。
其中匹配块(matched block)是测序序列种子映射到参考基因组块数据结构后，当一个块数据结构的种子计数器的数值大于阈值(7)时，则该块数据结构为测序序列的匹配块，也称是显著匹配块。
有益效果：
本发明的基于全局种子打分的候选比对区域优选方法中，每个候选比对区域的全局种子得分代表候选区域的重叠长度，通过全局打分可以有效优选重叠区域较长的候选区域，从而大幅降低进入两两局部比对的候选区域量；在全局种子打分模型基础上，设计了三代测序两两比对方法和参考基因组比对方法，这两种方法大大加速了三代测序序列比对过程和计算资源消耗量。
本发明的系统及其系统中的规则，实现了本发明的方法大幅降低目前三代测序需要的计算时间和资源，具有良好商业价值。
附图简要说明
图1：模块1中的block数据结构模型示意图
图2：模块2中的参考序列索引示意图
图3：模块2中的种子序列抽样规则示意图
图4：模块5中的全局种子投票打分模型示意图
具体实施方式
实施例1：基于全局种子打分优选的参考基因组比对方法
通过下面具体实施例来解释基于全局种子打分优选的参考基因组比对方法的相关步骤操作。
步骤1：建立参考基因组索引：
步骤1.1：以参考序列每个位点为起始，取13(种子(k-mer)，k＝13)个碱基长度的片段作为候选的种子序列，建立种子(k-mer)索引，参考基因组是由ATCG四个字母组成的一长串序列，实际长度可达10^9bp以上，为了方便统计采用编码的原理是用数字0替换字符A，数字1替换字符T，数字2为C和3为G。于是参考基因组转换成了由数字0，1，2构成一长串排列。于是ATCG字符组成的序列可以看成是4进制数据，依次从右到左各个字母编号(1,2,3…i)，通过如下公式计算得到：
Figure PCTCN2017098122-appb-000010
公式中i对应于序列中碱基的位置，NCi为相应位置字母对应的数字，将其转化为十进制数据，编码就反应seed序列的特征，如:CTTAACCGGAAAGG对应十进制2*4^13+1*4^12+1*4^11+0*4^10++0*4^9+2*4^8+2*4^7+3*4^6+3*4^5+0*4^4+0*4^3+0*4^2+3*4^1+3*4^0＝4624294.
步骤1.2：建立一个4^13大小整数数组用来记录参考基因组包含该种子(k-mer)数字编码个数SC[3^i]。SC[]数组的下标代表着种子(k-mer)的数字编码。SC每个元素的值代表着参考基因组中包含该元素下标数字编码序列的个数。数组中每个元素的初始值为0。
步骤1.3：逐步扫描参考基因组序列每个位置获得的种子(k-mer)，将其字符按步骤1.1要求转换成十进制编码，记录在步骤1.2中对应数组下标的值中，每记录一次则累加1。统计数组SC中的最大值，记为SC_MAX。
步骤1.4：建立一个指针数组*SI[4^i]指向AL地址，建立一个存储种子(k-mer)位置信息的数组AL[sum sc]，其中SI[i]＝AL+SC[i]，sum sc＝∑SC[n]，再次扫描参考基因组，那么SI[i][SC[i]]＝Location，其中Location代表所有候选位置信息，这是由于指针数组引用了AL中元素地址，Location最终存储在AL数组中，通过查找种子(k-mer)编码对应SI的下标就能在AL中找到参考基因组中seed候选位置信息和个数。记录了参考基因组中所有种子(k-mer)编码，出现次数及对应位置信息(如下表)。
SC Ref.index position
0 12，1001，10003，…
1 101,145,1193,…
2 144,1098,10129,…
3 132,13799,144353,…
步骤2：构建Z倍参考基因组块数据结构：
根据参考基因组长度L，分配一个L/1000+1的结构体数组，每个结构体包含种子匹配数和40个种子匹配位置对。并把每个块结构的种子计数设置为0，并分配一个L/1000+1的查阅表二维数组，里面记录后续测序种子匹配的块结构编号和该块种子匹配量。
步骤3：分割测序序列成若干个种子序列，种子序列提取规则为：
在测序序列中每隔20个碱基提取13个碱基长度的片段作为种子序列，并按照测序序列顺序进行编码。
步骤4：将所有种子序列比对到Z倍参考基因组块数据结构中
从步骤1的参考基因组索引中查到一个测序序列种子序列的所有参考基因组候选位置：(SLi,i＝1,2,...n)。按照步骤3将每个种子所有候选位置映射到参考基因组块数据结构中存储。当一个测序序列的种子比对到CR块数据结构的区域时，该块结构种子计数器将加1，并且该结构的种子匹配位置对将记录该种子在测序序列的位置和在CR块区域的相对位置。每个测序序列种子的参考基因组候选位置(SL)按照上述规则和公式映射到参考基因组块数据结构中，并且用查阅表(look-up table)记录所有匹配种子块数据结构编号(CR)。
步骤5：从块数据结构中选择局部序列比对的开始种子位置对(核心种子位置对)：
块数据结构查阅表将按照每个块结构的种子匹配数降序排序。当一个块数据结构的种子数大于一定阈值，该块结构被认为候选块结构。计算候选块结构中所有种子位置对之间的测序序列和参考基因组序列的长度之差(Dij)。根据公式计算两个种子位置对序列差异因子(DFij)，当序列差异因DFij＜0.2，子，块数据结构中的i和j号种子位置对互相投票支持，其投票得分各加一分。候选块数据结构所有种子位置对都进行上述两两投票打分后，投票得分最高种子对为局部比对开始种子对(核心种子位置对)。
步骤6：获取开始种子位置对的全局投票打分
根据开始种子位置对(SLk，SNk)，根据公式5和公式6估算测序序列可以跨越左边和右边相邻的块数据结构的数目VL和VR。测序序列覆盖相邻块数据结构的所有种子位置对将按照DF公式对开始种子位置对进行投票，从而获得开始位置对的全局投票 打分。当一个相邻块结构的80％种子位置对符合DFij＜0.2，即支持开始种子位置对，该块结构种子数将被设置为0，不再被考虑为候选块结构。
步骤7：根据全局得分，选择最高10个开始位置对进行局部两两序列比对
对上述步骤5和步骤6获得每个高于种子阈值的块数据结构的开始种子位置对的全局投票打分进行降序排序。选择全局投票得分最高20个开始位置对通过修改后的diff算法完成局部两两序列比对。开始位置对序列比对结果符合两个条件：当遇到符合重叠长读>1000和错误匹配率<0.20的开始位置对时，终止该测序序列序列比对过程，将该结果作为测序序列的序列比对结果输出。
步骤8：清理测序序列计算过程留下痕迹
将测序序列种子匹配的所有块结构的种子计数器重新设置为0，，并且查阅表(look-up table)记录清空。读取下一条测序序列重复步骤3到步骤8。直到所有测序数据完成参考基因组比对。
步骤9：二次精准搜索序列比对分析
提取没有匹配测序序列数据，通过二次搜索完成没有匹配数据搜索，将上述过程将ST变成10步长，之后块结构大小改为2000，其它参数不变，重复上述3-8步进行没有匹配序列的更精确的序列比对过程。
步骤10：程序并行化
将上述步骤2至步骤8通过基于共享内存变量空间pthread多线程包建立并行化程序，其中步骤1的参考基因组索引将放在多核共享内存。
实施例2：基于全局种子打分优选的两两比对方法
基于全局种子打分优选的两两比对方法基本上与实施例1的参考基因组实施过程类似，其不同之处如下：
第1步：数据分块和2G链接序列获取：扫描整个三代测序数据文件，按照2G文件大小分割三代测序数据，两个测序序列之间用N链接，并记录每个测序序列在2G链接参考序列上的起始位置和终止，并输出2G文件每个测序序列文件位置索引，方便后续计算链接参考序列上位置转化为每个测序序列的绝对位置。
第2步：与实施例1的步骤1相同。
第3步：与实施例1的步骤2类似，只是将块结构大小Z改为2000.
第4步：与实施例1的步骤3类似，只是将取种子(k-mer)步长ST改为10.
第5步：与实施例1的步骤4相同。
第6步：与实施例1的步骤5类似，只是两个种子位置对支持条件改为DFij<0.3。通过2G文件每个测序序列位置索引上每个测序链接序列的起始和终止位置，将核心位 置对上参考基因组上的位置转化为该位置所在的测序序列(read)的编号和测序序列上绝对位置。
第7步：与实施例1的步骤6类似，需要修改两个种子位置对支持条件改为DFij<0.3。两个测序序列重叠区域范围获取方式:根据显著匹配对的块编号和链接参考序列的测序序列起始位置，可以获取该匹配块定位测序序列编号和起始位置(S1，E1)，通过核心位置对的位置信息(参考基因组位置P1，待比对测序序列位置为P2)，待比对测序序列长度为L，可以得出，链接参考序列上匹配序列核心位置对左边长度Ll＝P1-S1和右边长度Lr＝E1-P1，待比对序列左边长度为P+，右边长度为L-P+，取两个左边长度较短者为左边重叠区域的长度，取右边长度较短者为右边长度，两个长度范围即延伸块结构的范围。
第8步：两两比对全局打分的输出：两两比对中只需获取2G中最高100全局种子得分的核心种子位置信息，不需要做局部序列比对，将最高100个核心种子位置信息转化为将转化成两条测序序列的绝对位置信息，最后输出两个重叠测序序列的编号、核心位置对的绝对位置信息和全局投票得分。
第9步与实施例1的步骤8类似,只是读取下一条测序序列后，执行第4步到第9步。
第10步程序并行化：将上述第4步到第9步通过基于共享内存变量空间pthread多线程包建立并行化程序，其中第2步的参考基因组索引将放在多核共享内存。
第11步每个数据块的两两比对：数据块1要与数据块1-n进行两两比对，数据块2要与数据块2-n进行两联比对，依次类推完成所有序列两两比对，为了匹配相同两个序列进行两个两两比对，比对过程中，待比对测序序列标号要大于显著匹配块的测序序列的编号才进行后续全局种子投票分析。
通过实施例1和实施例2的方法，下载五个真实物种(E.coli,Yeast,A.Thaliana,D.Melanogaster和Human)的PacBio数据集和三个真实物种(E.coli,B.anthracis和Y.pestis)的nanopore数据集进行测试我们MECAT效果。两两比对软件比较随机提取500M数据进行比较，在PacBio数据集中，我们软件MECAT两两比对的速度是MHAP和Daligner软件的2-8倍；在nanopore数据集中，MECAT速度是MHAP和Daligner的5-10倍。参考基因组软件比较使用整个数据集进行比较，在PacBio数据集中，我们软件MECAT两两比对的速度是BLASR和BWA软件的5-70倍；在nanopore数据集中，MECAT速度BLASR和BWA的4-5倍。
Figure PCTCN2017098122-appb-000011
上述表格时间单位是核时

Claims (31)

  1. 一种基于全局种子打分优选的三代测序序列比对系统，该系统包含模块1、模块2、模块3、模块4和模块5，模块1嵌合快速查找显著候选重叠区域的block数据结构模型，模块2嵌合参考基因组block数据结构的映射规则，模块3嵌合参考基因组索引和测序序列种子序列抽样规则，模块4嵌合匹配块种子匹配数与灵敏度数学模型，模块5嵌合基于块数据结构的全局种子打分模型，其中模块5包含模块5.1、模块5.2和模块5.3，模块5.1嵌合两个种子对之间的两个序列的距离差异因子，模块5.2嵌合两两种子投票打分获取核心匹配种子位置对规则，模块5.3嵌合延伸投票打分获取核心位置对的全局种子投票得分规则。
  2. 根据权利要求1所述的系统，其特征在于所述系统还包含模块6，模块6嵌合基于全局种子打分的优选和使用规则。
  3. 根据权利要求1或2所述的系统，其特征在于系统中的模块1嵌合快速查找显著候选重叠区域的block数据结构模型，所述快速查找显著候选重叠区域的block数据结构模型为：设Z为block数据结构的块比例，对于参考基因组每Z个碱基建立一个block数据结构，并顺序编号，用于比对过程中快速将种子序列定位到候选比对区域。
  4. 根据权利要求1或2所述的系统，其特征在于所述系统中，模块2嵌合参考基因组block数据结构的映射规则，参考基因组block数据结构的映射规则为：
    通过测序序列每个种子编码查询参考基因组索引获得每个种子基因组的精确位置，并用每个种子的精确位置按照公式1的规则映射到上述块结构：
    Figure PCTCN2017098122-appb-100001
    其中Z表示块结构碱基区域大小，CR表示块结构的序号，CL为种子在参考基因组的在块结构的相对准确位置，SLi表示参考基因组候选位置。
  5. 根据权利要求1或2所述的系统，其特征在于所述模块3嵌合参考基因组索引和测序序列种子序列抽样规则，参考基因组索引和测序序列种子序列抽样规则为：
    以参考基因组每个位点为起始，取k＝13个碱基长度的片段作为种子序列，建立种子的4进制编码与其对应起始位置的哈希表。
  6. 根据权利要求1或2所述的系统，其特征在于所述模块4嵌合匹配块种子匹配数与灵敏度数学模型，其中匹配块种子匹配数与灵敏度数学模型包括参考基因组块种子匹配数与灵敏度数学模型和两两比对中块种子匹配数与灵敏度数学模型。
  7. 根据权利要求6所述的系统，其特征在于所述系统中的参考基因组块种子匹配数与灵敏度数学模型如下：
    假设所有种子比对是独立事件，种子的匹配概率初步符合二项式分布，在参考基因 组比对过程中，种子匹配概率用如下公式2计算：
    Povl＝(1-e)k   (公式2)
    公式2中，当e为0.15，块大小为1000，种子抽样步长为20和k为13时，每个块的抽样数为
    Figure PCTCN2017098122-appb-100002
    两个匹配块平均种子匹配个数为
    Figure PCTCN2017098122-appb-100003
  8. 根据权利要求6所述的系统，其特征在于所述系统中两两比对中块种子匹配数与灵敏度数学模型如下：
    假设所有种子比对是独立事件，种子的匹配概率初步符合二项式分布，在两两比对过程中，种子匹配概率通过如下公式3计算：
    Figure PCTCN2017098122-appb-100004
    公式3中，当e为0.15，块大小为2000，种子抽样步长为5和k为13时，每个块的抽样数为
    Figure PCTCN2017098122-appb-100005
    两个匹配块平均种子匹配个数为
    Figure PCTCN2017098122-appb-100006
  9. 根据权利要求1或2所述的系统，其特征在于所述模块5嵌合基于块数据结构的全局种子打分模型，基于块数据结构的全局种子打分模型如下：
    对参考基因组和三代测序序列分别建立种子k＝13的哈希表，同时将基因组和序列分成大小为1000bp的数据块。
  10. 根据权利要求9所述的系统，其特征在于所述模块5中的全局种子打分从显著匹配块开始,其过程包含如下模块5.1，模块5.2和模块5.3，其中模块5.1嵌合两个种子对之间的两个序列的距离差异因子，模块5.2嵌合两两种子投票打分获取核心匹配种子位置对规则，模块5.3嵌合延伸投票打分获取核心位置对的全局种子投票得分规则。
  11. 根据权利要求2所述的系统，其特征在于所述系统中模块6嵌合基于全局种子打分的优选和使用规则，基于全局种子打分的优选和使用规则如下：
    对所有显著匹配块候选区域进行全局种子投票打分获取每个显著匹配块的核心种子对和全局种子得分，根据每个候选区域的核心种子对和全局种子得分判定核心种子位置对进入后续局部序列比对分析，其中全局种子投票打分方法为基于全局种子打分优选的参考基因组比对方法或基于全局种子打分优选的两两比对方法，其中判定方法如下：
    (1)应用基于全局种子打分优选的参考基因组比对方法获得的结果，当参考基因组选择最高10个核心种子位置对所在的区域有效候选区域，这些核心种子位置对可以进入后续局部序列比对分析；
    (2)应用基于全局种子打分优选的两两比对方法获得结果，选择最高100个核心种子位置对所在的区域有效候选区域，这些核心种子位置对可以进入后续局部序列比对分析。
  12. 一种基于全局种子打分优选的三代测序序列比对方法，所述三代测序序列比对 方法为基于全局种子打分优选的参考基因组比对方法和基于全局种子打分优选的两两比对方法中的一种或两种，所述基于全局种子打分优选的参考基因组比对方法和基于全局种子打分优选的两两比对方法执行权利要求1所述模块1、模块2、模块3、模块4和模块5中的至少2个以上的模块。
  13. 一种基于全局种子打分优选的三代测序序列比对方法，所述方法为基于全局种子打分优选的参考基因组比对方法，所述参考基因组比对方法包括如下步骤：
    步骤1.1：建立参考基因组索引
    步骤1.2：构建参考基因组块数据结构
    步骤1.3：分割测序序列成若干个种子序列
    步骤1.4：将所有种子序列映射到Z倍参考基因组块数据结构中
    步骤1.5：获取显著块匹配区域的核心种子位置对
    步骤1.6：获取核心种子位置对的全局投票打分
    步骤1.7：选择最高n个核心位置对进行局部两两序列比对
    步骤1.8：二次精准参考基因组序列比对。
  14. 根据权利要求13所述的方法，其特征在于步骤1.1所述的建立参考基因组索引方法为：
    应用权利要求5所述的模块3，从参考基因组中每个碱基位置提取k长度种子序列，也就是，相邻的种子之间没有间隔，参考基因组所有的碱基被建立种子索引。
  15. 根据权利要求13所述的方法，其特征在于步骤1.2所述的构建参考基因组块数据结构方法为：
    应用权利要求3所述的模块1，将参考基因组每Z个碱基区域建立一个块数据结构，每个块数据结构用于记录测序序列种子在该结构代表参考基因组区域的匹配情况。
  16. 根据权利要求15所述的方法，其特征在于每个块数据结构由种子匹配计算器、40个种子匹配候选种子位置对组成。
  17. 根据权利要求13所述的方法，其特征在于所述步骤1.3所述的分割测序序列序列成若干个种子序列方法为：
    应用权利要求5所述的模块3，在测序序列中按照ST＝20步长提取种子的种子序列，每个种子有k个碱基组装，并按照测序序列顺序进行编码。
  18. 根据权利要求13所述的方法，其特征在于所述步骤1.4所述的将所有种子序列映射到Z倍参考基因组块数据结构中的方法为：
    一个测序序列种子序列的所有参考基因组候选位置从步骤1.1的参考基因组索引中查找，应用权利要求4所述的模块2，将每个种子所有候选位置映射到参考基因组块数 据结构中存储。并且用查阅表记录所有匹配种子块数据结构编号，查阅表记录着种子匹配的块区域编号和对应块区域的种子匹配数，每个块区域在查阅表中唯一记录。
  19. 根据权利要求13所述的方法，其特征在于所述步骤1.5所述的获取显著块匹配区域的核心种子位置对方法为：
    当一个块数据结构的种子数大于7时，该块结构被认为显著块匹配结构，显著块匹配结构的局部比对的核心匹配位置对将通过该块结构中所有种子对两两投票打分确定，按照权利要求10所述的模块5.2进行两两投票打分获取该显著匹配块的核心种子位置对。
  20. 根据权利要求13所述的方法，其特征在于所述步骤1.6所述的获取核心种子位置对的全局投票打分方法为：
    应用权利要求10所述的模块5.3，将测序序列所覆盖的核心种子位置对相邻块结构中位置对取出，对核心位置对进行单向投票打分，获取核心位置对全局种子得分。并将相邻块结构的80％种子位置对支持核心种子位置对的块结构的种子数将被设置为0.
  21. 根据权利要求13所述的方法，其特征在于所述步骤1.7所述的选择最高n个核心位置对进行局部两两序列比对方法为：
    通过步骤1.5和步骤1.6，获得每个高于种子阈值的块数据结构的核心种子位置对和全局投票打分，之后，对所有核心位置对的全局投票得分进行降序排序，选择全局投票得分最高10个核心位置对通过字基于最长公共字串的字符差异比对方法(diff)完成局部两两序列比对，对用nanopore，采用经典的局部匹配的方法(smith-waterman)进行局部两两比对，如果核心位置对序列比对结果符合两个条件：重叠长读>1000和错误匹配率<0.20，认为该测序序列已找到正确参考基因组匹配位置，按照全局比对得分顺序进行两两序列比对，当遇到符合上述条件的核心位置对时，终止该测序序列序列比对过程，将该结果作为测序序列的序列比对结果输出。
  22. 根据权利要求13所述的方法，其特征在于所述步骤1.8所述的二次精准参考基因组序列比对方法为：
    针对少数测序序列的块匹配种子量较少，而且布局均一，不能被步骤1.4参数搜索到，如果上述过程ST步长分割和Z数据结构没有获得搜索结果输出，将执行步骤1.3的ST变成ST/2步长，之后的块大小为2Z，其它参数不变，重复上述步骤3到步骤1.7进行更精确的序列比对过程。
  23. 一种基于全局种子打分优选的三代测序序列比对方法，所述方法为基于全局种子打分优选的两两比对方法，所述两两比对方法包括如下步骤：
    步骤2.1：三代测序数据分块和测序序列链接成类似参考基因组
    步骤2.2：建立参考基因组索引
    步骤2.3：构建链接参考序列的块数据结构
    步骤2.4：分割测序序列序列成若干个种子序列
    步骤2.5：将所有种子序列映射到Z倍链接参考序列的块数据结构中
    步骤2.6：获取显著块匹配区域的的核心种子位置对
    步骤2.7：获取核心种子位置对的全局投票打分
    步骤2.8：选择最高n个核心位置对的候选区域输出结果。
  24. 根据权利要求23所述的方法，其特征在于所述步骤2.1所述的三代测序数据分块和测序序列链接成类似参考基因组方法为：
    将三代测序数据集分成2G大小数据块，链接2G数据块内的测序序列成2G的一条参考序列，两条测序序列链接出添加一个N字母，记录每个测序序列在2G参考序列上的起始位置，为后续寻找两个测序序列重叠的起始位置。
  25. 根据权利要求23所述的方法，其特征在于所述步骤2.2所述的建立参考基因组索引方法为：
    应用权利要求5所述的模块3，从链接后的2G参考序列中每个碱基位置提取k长度种子序列，也就是，相邻的种子之间没有间隔。参考基因组所有的碱基将被建立种子索引。
  26. 根据权利要求23所述的方法，其特征在于所述步骤2.3所述的构建链接参考序列的块数据结构方法为：
    应用权利要求3所述的模块1，链接参考序列每Z个碱基区域建立一个块数据结构，每个块数据结构用于记录测序序列种子在该结构代表建链接参考序列区域的匹配情况，每个块数据结构由种子匹配计算器、40个种子匹配候选种子位置对组成。
  27. 根据权利要求23所述的方法，其特征在于所述步骤2.4所述的分割测序序列序列成若干个种子序列方法为：
    应用权利要求5所述的模块3，在测序序列中按照ST＝10步长提取种子的种子序列，每个种子有k个碱基组装，并按照测序序列顺序进行编码。
  28. 根据权利要求23所述的方法，其特征在于所述步骤2.5所述的将所有种子序列映射到Z倍链接参考序列的块数据结构中的方法为：
    一个测序序列种子序列的所有参考基因组候选位置从步骤2.1的参考基因组索引中查找，应用权利要求4所述的模块2，将每个种子所有候选位置映射到链接参考序列块数据结构中存储。并且用查阅表记录所有匹配种子块数据结构编号，查阅表记录着种子匹配的块区域编号和对应块区域的种子匹配数，每个块区域在查阅表中唯一记录。
  29. 根据权利要求23所述的方法，其特征在于所述步骤2.6所述的获取显著块匹配区域的的核心种子位置对方法为：
    当一个块数据结构的种子数大于7时，该块结构被认为显著块匹配结构，显著块匹配结构的局部比对的核心匹配位置对将通过该块结构中所有种子对两两投票打分确定，应用权利要求10所述的模块5.2进行两两投票打分获取该显著匹配块的核心种子位置对。
  30. 根据权利要求23所述的方法，其特征在于所述步骤2.7所述的获取核心种子位置对的全局投票打分方法为：
    应用权利要求10所述的模块5.3，将测序序列所覆盖的核心种子位置对相邻块结构中种子位置对取出，通过2G链接参考序列上每个测序序列的起始位置和显著块匹配的编号可以定位核心种子对来源2G中的测序序列的编号，根据两个测序序列重叠情况，获取相邻块匹配的范围，对重叠区域的核心位置对进行单向投票打分，获取核心位置对全局种子得分。并将相邻块结构的80％种子位置对支持核心种子位置对的块结构的种子数将被设置为0。
  31. 根据权利要求23所述的方法，其特征在于所述步骤2.8所述的选择最高n个核心位置对的候选区域输出结果的方法为：
    通过步骤2.5和步骤2.6，获得每个高于种子阈值的块数据结构的核心种子位置对和全局投票打分，然后，对所有核心位置对的全局投票得分进行降序排序，选择全局投票得分最高100个核心位置对，将每个核心位置对信息通过测序序列在2G链接参考序列起始位置转化成两条测序序列的绝对位置信息，最后输出两个重叠测序序列的编号、核心位置对的绝对位置信息和全局投票得分。
PCT/CN2017/098122 2017-06-02 2017-08-18 一种基于全局种子打分优选的三代测序序列比对方法 WO2018218788A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201710412287.2A CN107256335A (zh) 2017-06-02 2017-06-02 一种基于全局种子打分优选的三代测序序列比对方法
CN201710412287.2 2017-06-02

Publications (1)

Publication Number Publication Date
WO2018218788A1 true WO2018218788A1 (zh) 2018-12-06

Family

ID=60023899

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2017/098122 WO2018218788A1 (zh) 2017-06-02 2017-08-18 一种基于全局种子打分优选的三代测序序列比对方法

Country Status (2)

Country Link
CN (1) CN107256335A (zh)
WO (1) WO2018218788A1 (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108460246B (zh) * 2018-03-08 2022-02-22 北京希望组生物科技有限公司 一种基于三代测序平台的hla基因分型方法
CN108776749B (zh) * 2018-06-05 2022-05-03 北京诺禾致源科技股份有限公司 测序数据的处理方法及装置
CN108920902A (zh) * 2018-06-29 2018-11-30 郑州云海信息技术有限公司 一种基因序列处理方法及其相关设备
CN108985008B (zh) * 2018-06-29 2022-03-08 郑州云海信息技术有限公司 一种快速比对基因数据的方法和比对系统
CN109326325B (zh) * 2018-07-25 2022-02-18 郑州云海信息技术有限公司 一种基因序列比对的方法、系统及相关组件
CN110517727B (zh) * 2019-08-23 2022-03-08 苏州浪潮智能科技有限公司 序列比对方法及系统
CN111190915B (zh) * 2020-01-02 2023-05-16 腾讯科技（深圳）有限公司 一种道具标识或角色标识的确定方法、服务器及存储介质
CN111627496B (zh) * 2020-05-09 2022-05-17 苏州浪潮智能科技有限公司 一种哈希表的压缩方法、系统及相关装置
CN114520024B (zh) * 2022-01-17 2024-03-22 浙江天科高新技术发展有限公司 一种基于k-mer的序列联配方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002026934A2 (en) * 2000-09-28 2002-04-04 New York University System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map
CN104951672A (zh) * 2015-06-19 2015-09-30 中国科学院计算技术研究所 一种第二代、三代基因组测序数据联用的拼接方法及系统
CN105389481A (zh) * 2015-12-22 2016-03-09 武汉菲沙基因信息有限公司 一种三代全长转录组中可变剪切体的检测方法
CN106022002A (zh) * 2016-05-17 2016-10-12 杭州和壹基因科技有限公司 一种基于三代PacBio测序数据的补洞方法
CN106021997A (zh) * 2016-05-17 2016-10-12 杭州和壹基因科技有限公司 一种三代PacBio测序数据的比对方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105989249B (zh) * 2014-09-26 2019-03-15 南京无尽生物科技有限公司 用于组装基因组序列的方法、系统及装置
WO2016148650A1 (en) * 2015-03-17 2016-09-22 Agency For Science, Technology And Research Bioinformatics data processing systems

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2002026934A2 (en) * 2000-09-28 2002-04-04 New York University System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map
CN104951672A (zh) * 2015-06-19 2015-09-30 中国科学院计算技术研究所 一种第二代、三代基因组测序数据联用的拼接方法及系统
CN105389481A (zh) * 2015-12-22 2016-03-09 武汉菲沙基因信息有限公司 一种三代全长转录组中可变剪切体的检测方法
CN106022002A (zh) * 2016-05-17 2016-10-12 杭州和壹基因科技有限公司 一种基于三代PacBio测序数据的补洞方法
CN106021997A (zh) * 2016-05-17 2016-10-12 杭州和壹基因科技有限公司 一种三代PacBio测序数据的比对方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
XIAO, C.L. ET AL.: "MECAT: Fast mapping, error correction, and the novo assembly for single-molecule sequencing reads (incl. suppl. data)", NATURE METHODS, vol. 14, no. 11, pages 1072 - 1074, XP055632797, DOI: 10.1038/nmeth.4432 *
XIAO, C.L., 18 December 2017 (2017-12-18), Retrieved from the Internet <URL:https://github.com/xiaochuanle/MECAT> *
XIAO, C.L.: "MECAT 1.0 installing software, read me document and screenshot display of Mecat 1.0 installing software documents (non-official translation)", 18 January 2017 (2017-01-18), Retrieved from the Internet <URL:https://github.com/xiaochuanle/MECAT/bolb/master/Mecat1.0.zip> *

Also Published As

Publication number Publication date
CN107256335A (zh) 2017-10-17

Similar Documents

Publication Publication Date Title
WO2018218788A1 (zh) 一种基于全局种子打分优选的三代测序序列比对方法
US11810648B2 (en) Systems and methods for adaptive local alignment for graph genomes
US10600217B2 (en) Methods for the graphical representation of genomic sequence data
CA2424031C (en) System and process for validating, aligning and reordering genetic sequence maps using ordered restriction map
US20160259880A1 (en) Systems and methods for genomic pattern analysis
CN107403075B (zh) 比对方法、装置及系统
US20110196872A1 (en) Computational Method for Comparing, Classifying, Indexing, and Cataloging of Electronically Stored Linear Information
US10192028B2 (en) Data analysis device and method therefor
CN107133493B (zh) 基因组序列的组装方法、结构变异探测方法和相应的系统
JP2008547080A (ja) ダイタグ配列の処理および／またはゲノムマッピングの方法
WO2018218787A1 (zh) 一种基于局部图的三代测序序列校正方法
CN115631789B (zh) 一种基于泛基因组的群体联合变异检测方法
Gärtner et al. Coordinate systems for supergenomes
EP3938932B1 (en) Method and system for mapping read sequences using a pangenome reference
CN103294932A (zh) 用于碱基序列分析的参考序列处理系统及方法
US11482304B2 (en) Alignment methods, devices and systems
JP2023014025A (ja) 方法、コンピュータプログラム、及びコンピュータシステム（文字列類似度決定）
Chen et al. CGAP-align: a high performance DNA short read alignment tool
CN107688727B (zh) 生物序列聚类和全长转录组中转录本亚型识别方法和装置
Warnke-Sommer et al. Parallel NGS assembly using distributed assembly graphs enriched with biological knowledge
CN115861649A (zh) 基于注意力机制和加权概念格的古建筑图像语义完备方法
Zhang MST Based Ab Initio Assembler of Expressed Sequence Tags
CN115775592A (zh) circRNA检测方法、计算机程序产品及系统
Bao Algorithms for Reference Assisted Genome and Transcriptome Assemblies

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17911456

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 18/06/2020)

122 Ep: pct application non-entry in european phase

Ref document number: 17911456

Country of ref document: EP

Kind code of ref document: A1