CN110246545A - 一种序列的校正方法及其校正装置 - Google Patents

一种序列的校正方法及其校正装置 Download PDF

Info

Publication number
CN110246545A
CN110246545A CN201910493581.XA CN201910493581A CN110246545A CN 110246545 A CN110246545 A CN 110246545A CN 201910493581 A CN201910493581 A CN 201910493581A CN 110246545 A CN110246545 A CN 110246545A
Authority
CN
China
Prior art keywords
sequence
mer
base
site
correction
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
CN201910493581.XA
Other languages
English (en)
Other versions
CN110246545B (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.)
Wuhan Future Group Biological Science And Technology Co Ltd
Original Assignee
Wuhan Future Group Biological Science And Technology Co Ltd
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 Wuhan Future Group Biological Science And Technology Co Ltd filed Critical Wuhan Future Group Biological Science And Technology Co Ltd
Priority to CN201910493581.XA priority Critical patent/CN110246545B/zh
Publication of CN110246545A publication Critical patent/CN110246545A/zh
Application granted granted Critical
Publication of CN110246545B publication Critical patent/CN110246545B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30—Detection of binding sites or motifs
    • 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

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

Abstract

本申请公开了一种序列的校正方法及其校正装置，涉及生物信息技术领域，具体而言，该校正方法通过对待校正序列进行校正，通过k‑mer方式对待校正序列进行划分，从正序统计每个最后位点上每一种碱基对应的k‑mer片段的概率，从反序，将待校正序列上最后一个位点的碱基替换为获得最高概率的k‑mer片段的碱基，根据最后一个位点的碱基对应的最高概率的k‑mer片段，确认倒数第二个位点的碱基，依次类推，得到第一校正序列。该校正方法能过在有效的时间内，对待校正序列进行校正，得到精确度更高的基因序列。

Description

一种序列的校正方法及其校正装置
技术领域
本发明涉及生物信息技术领域，具体而言，涉及一种序列的校正方法及其校正装置。
背景技术
从1977年，第一代DNA测序技术(Sanger)发展至今三十多年时间，测序技术已经取得了相当大的发展，现有通用的测序技术有：一代测序、二代测序和三代测序。
通过测序技术得到的初步结果通常为多个测序片段，其中，二代测序具有准确性高，耗时短的优点，但现有基于二代测序数据的序列校正方法存在耗时长，准确度不高等问题。
发明内容
本发明提供了一种序列的校正方法，该校正方法能够有效避免或减少现有技术在对待校正序列的校正过程中，出现的耗时长，准确度不高等技术问题。
具体地，本申请提供的一种序列的校正方法，其包括：将测序得到的多个测序片段与待校正序列进行比对，将多个测序片段排列至待校正序列上相对应的位置，获得第一排列结果。
需要说明的是，多个测序片段可以为二代序列、三代测序或其测序结果优化后的结果，还可以为后续发展的其他测序技术对基因序列的测序结果。
待校正序列可以包括任何需要校正的序列，例如，通过使用现有测序技术获得的多个测序片段组装得到的基因序列；多个测序片段组装得到的基因序列的部分有效序列；或校正后的校正序列等。现有测序技术优选为二代或三代测序技术。
在本发明的优选实施例中，待校正序列包括：测序片段组装后的基因序列、含有低质量区域的序列、含有低覆盖区域的序列、第一校正序列、第二校正序列、第三校正序列或第四校正序列。
进一步地，以k-mer的方式对上述待校正序列进行划分，得到多个k-mer片段，基于第一排列结果，对测序片段进行对应的k-mer划分。
确认每个k-mer片段的最后位点对应的测序片段，基于第一排列结果，从正序统计每个k-mer片段最后位点上每一种碱基对应的k-mer片段的概率。
具体地，序列上每一个位点，都存在5种可能，分别为A、T、C、G或缺失(使用“-”表示)。测序结果中存在上百万条测序片段(reads)，而测序是通常随机打断的，在这些reads的位置不清楚的情况下，需要根据测序片段重叠的部分来尽量还原测序序列。k-mer，monomeric unit(mer)，相当于nt或者bp，100mer DNA相当于每一条链有100nt，那么整条链就是100bp。一般长短为m的reads可以分成m-k+1个k-mers。在本发明实施例中，以k-mer为3举例说明，以k-mer为3对上述待校正序列进行划分，得到多个3-mer片段，比如位点1、2、3为第一个k-mer，位点2、3、4为第二个k-mer，依次类推。
本申请根据碱基组合的形式来校正序列，从正序确认每个k-mer片段的最后位点，最后位点上可能存在A、T、C和G四种碱基，以及缺失的情况。
从正序统计每个位点上每一种碱基对应的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片段的分值，分值越高，概率越大，计算公式包括：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C；
其中，p为最后位点在待校正序列上的位置，b为碱基A、T、G、C或缺失，score(p，b)为p位置上碱基b的分值，score(p-1，b)为p-1位置上碱基b的分值；countk-mer为该k-mer对应的特定碱基组合出现的次数；C为k-mer区域的有效测序深度。
优选地，将第一校正序列作为待校正序列，进行迭代校正。待校正序列的序列不同，可能会导致多个测序片段在待校正序列上排列分布的不同，从而影响后续的校正结果。提高待校正序列的精确性，能提升排列结果的准确度，从而进一步提高后续的校正结果。
优选地，在一些优选的实施例中，以k-mer方式对待校正序列进行划分为：以k-mer按照预设值为3的方式对待校正序列进行划分。k-mer的预设值不同，可能对后续的校正的结果产生影响，本发明的校正方法中，k-mer为3时，校正结果更加精确。
进一步地，在一些实施例中，本发明提供的校正方法还低质量校正。
具体地，低质量校正包括：在待校正序列上，根据第一排列结果，将所有最后位点上，获得最高概率k-mer片段对应的碱基占该位点上k-mer片段总数的比值小于低质量占比预设值的位点，标记为低质量位点。需要说明的是，在低质量校正中，k-mer片段的概率采用的是由上述正序统计中计算出来的概率。
进一步地，低质量校正包括对低质量位点进行低质量区间校正：将待校正序列上出现两个以上的低质量位点，且两个以上低质量位点之间最大间隔长度小于或等于低质量间隔预设值的区间划分为低质量区间，将低质量区间的序列替换为该区间内出现概率最高的测序片段的序列，得到第二校正序列。如果出现并列最高测序片段的情况，则随机选择。
进一步地，将低质量区间的序列替换为该区间内出现概率最高的测序片段的序列具体包括：将多个测序片段排列至第一校正序列上相对应的位置，获得第二排列结果；
基于第二排列结果，确定低质量区间内重复次数最多的测序片段，将待校正序列上低质量区间内的序列校正为重复次数最多的测序片段对应的序列，得到第二校正序列。
优选地，低质量占比预设值为80％；
优选地，低质量间隔预设值小于等于测序片段的长度；
优选地，低质量间隔预设值为50碱基。
优选地，校正方法还包括将第二校正序列作为待校正序列，进行迭代校正。
进一步地，低质量区间的校正具有快速、高效的优点，能够在短时间内大量将可能发生的错误的位点继续校正。
在一些优选实施例中，上述校正方法还包括：将与前后相邻低质量位点的间隔距离大于低质量间隔预设值的低质量位点的碱基进行低质量位点校正，获得第三校正序列；
优选地，低质量位点校正包括：
将与前后相邻低质量位点的间隔距离大于低质量间隔预设值的低质量位点的碱基串联成低质量长序列；
将多个测序片段排列至第二校正序列上相对应的位置，挑选出与低质量长序列对应的排列，获得第三排列结果；
以k-mer方式对低质量长序列进行划分，得到多个k-mer片段，基于第三排列结果，对测序片段进行对应的k-mer划分；
确认每个k-mer片段的最后位点对应的测序片段，基于第三排列结果，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率。
从反序，将低质量长序列上最末碱基替换为获得最高概率的k-mer片段的碱基，根据最末碱基对应的最高概率的k-mer片段，确认倒数第二位点的碱基；将倒数第二位点的碱基作为下一个k-mer片段的最后位点，根据倒数第二位点碱基对应的最高概率的k-mer片段，确认倒数第三位点的碱基，依次类推，使用低质量长序列校正确认后的碱基替换待校正序列中的对应碱基，得到第三校正序列。
优选地，低质量间隔预设值小于等于测序片段的长度。
优选地，低质量间隔预设值为50碱基。
优选地，统计每个最后位点上每一种碱基对应的k-mer片段的概率包括：按照计算公式，计算每个最后位点上每一种碱基对应的k-mer片段的分值，分值越高，概率越大，计算公式为：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C；
其中，p为最后位点在待校正序列上的位置，b为碱基A、T、G、C或缺失，score(p，b)为p位置上碱基b的分值，score(p-1，b)为p-1位置上碱基b的分值；countk-mer为该k-mer对应的特定碱基组合出现的次数；C为k-mer区域的有效测序深度。
优选地，以k-mer方式对低质量长序列进行划分为：以k-mer为3对方式对低质量长序列进行划分。
优选地，将第三校正序列作为待校正序列，进行迭代校正。
进一步地，在待校正序列或第三校正序列中，还可能存在一些区域，对应匹配的测序片段较少，针对这些采用上述正序统计、反序校正的方法或上述低质量校正方法无法进一步确认的位点，本发明实施例设置了低覆盖校正，获得第四校正序列。
具体地，低覆盖校正包括：
将多个测序片段排列至第二校正序列或第三校正序列上相对应的位置，挑选出低覆盖区域对应的排列，获得第四排列结果；
以k-mer方式对低覆盖区域进行划分，得到多个k-mer片段，基于第四排列结果，对测序片段进行对应的k-mer划分；
基于第四排列结果，确认每个k-mer片段的最后位点对应的测序片段，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率，基于k-mer片段的概率，校正低覆盖区域；
从反序将低覆盖区域上最末碱基替换为获得最高k-mer片段的概率的碱基，根据最末碱基对应的k-mer片段，确认倒数第二位点的碱基；将倒数第二位点的碱基作为下一个低覆盖k-mer片段的最后位点，确认倒数第三位点的碱基，依次类推，获得第四校正序列。
优选地，第四校正序列为将低覆盖校正后的结果与上述中的低质量区间校正后的结果的结合。
优选地，统计每个最后位点上每一种碱基对应的低覆盖k-mer片段的概率包括：按照计算公式，计算每个最后位点上每一种碱基对应的k-mer片段的分值，分值越高，概率越大，计算公式为：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C；
其中，p为最后位点在待校正序列上的位置，b为碱基A、T、G、C或缺失，score(p，b)为p位置上碱基b的分值，score(p-1，b)为p-1位置上碱基b的分值；countk-mer为该k-mer对应的特定碱基组合出现的次数；C为k-mer区域的有效测序深度。
优选地，以k-mer方式对低覆盖序列进行划分为：以k-mer为3对方式对低覆盖序列进行划分；
优选地，将第四校正序列作为待校正序列，进行迭代校正。
进一步地，在序列比对前，校正方法包括将多个测序片段进行预处理；
预处理包括将多个测序片段进行以下一种或多种方式：
1)过滤掉没有跨过校正区域的测序片段；
2)过滤掉比对质量低于比对阈值的测序片段；
3)当测序片段只有部分区域能与待校正匹配时，切除测序片段两端无法匹配的序列，过滤掉当切除的序列长度大于切除预设值的测序片段，切除预设值优选为10％以上。
具体地，测序片段的比对质量的计算公式为：
Score比对质量＝-10×log10Q；
其中，其中Q为该校正reads比对位置错误的概率，Score比对质量即为比对质量分值。
此外，本发明实施例还提供了一种序列的校正装置，校正装置用于执行如上述的校正方法。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将对本发明实施例中的技术方案进行清楚、完整地描述。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。
以下结合实施例对本发明的特征和性能作进一步的详细描述。
实施例1
本实施例提供了一种序列的校正方法，具体如下。
将测序得到的多个测序片段与待校正序列进行比对，将多个测序片段排列至待校正序列上相对应的位置，获得第一排列结果；
以k-mer为3对待校正序列进行划分，根据第一排列结果将测序片段按照对应的3-mer进行划分，得到多个3-mer片段；
确认每个3-mer片段的最后位点，基于第一排列结果，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率；
统计每个最后位点上每一种碱基对应的k-mer片段的概率包括：按照计算公式，计算每个最后位点上每一种碱基对应的k-mer片段的分值，分值越高，概率越大，计算公式为：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C；
其中，p为最后位点在待校正序列上的位置，b为碱基A、T、G、C或缺失，score(p，b)为p位置上碱基b的分值，score(p-1，b)为p-1位置上碱基b的分值；countk-mer为该k-mer对应的特定碱基组合出现的次数；C为k-mer区域的有效测序深度。
从反序，将待校正序列上最后一个位点的最末碱基替换为获得最高概率的k-mer片段的碱基，根据最后一个位点最末的碱基对应的最高概率的k-mer片段，确认倒数第二个位点的碱基；将倒数第二个位点的碱基作为下一个k-mer片段的最后一个位点，根据倒数第二个位点碱基对应的最高概率的k-mer片段，确认倒数第三个位点的碱基，依次类推，得到第一校正序列。
为了获得更准确的校正序列，在本实施例中，采用优选方案，使用高准确度的测序片段进行校正，高准确度的测序片段由二代测序获得。
示例1
为了更清楚的说明校正方法，通过示例进行演算，具体如下。
表1为示例数据，测序片段包括有Read1～Read5，设定位置1为起始位置，位置7为结束位置，使用本实施例1的方法进行校正。其中，k-mer设定为3-mer，各位置有效覆盖度均是5。
表1数据信息
使用本实施例的方法，获得的校正链从位置1至7依次为ACAGACC。
若不进行反向回溯，仅考虑单个碱基，而没有考虑相邻碱基之间的连接关系，即通过最高得分进行选择。位置3，最高得分碱基为A(-3)；位置4，最高得分碱基为G(-6)；位置5，最高得分碱基为T(-9)；位置6，最高得分碱基为G(-12)；位置7，最高得分碱基为C(-15)。最后确定的得分链位置3至位置7为AGTGC，与回溯后的结果不一致，体现了回溯的效果。
示例2
进行k-mer划分时，k值越大，k-mer片段的特异性越好，校正越准确。但是，当k值越大，完全一致k-mer的概率越小。当测序准确率为p时，出现完全一致k-mer的概率是pk。
通过前述score(p,b)的计算公式可知，得分与完全一致k-mer数目相关，若数目过低，即计入校正的序列过少，会严重影响校对效果。以现阶段三代测序平均准确率为85％为例，2-mer完全一致的概率是72.25％，3-mer完全一致的概率是61.41％，4-mer完全一致的概率是52.20％。显而易见的，4-mer及大于4-mer时，概率过低。为了进一步说明k-mer划分的技术效果，通过示例进行演算，比较2-mer与3-mer的校正效果，具体如下。
表2为示例数据，测序片段包括有Read1～Read7，设定位置1为起始位点，位置5为结束位点，进行处理。
表2数据信息
示例2-1
设定k-mer为3，按照本实施例方法进行校正，通过回溯确定的得分链从位置1至5依次为GGAGT。
示例2-2
采用k-mer为2进行划分，参照本实施例方法进行校正，通过回溯确定的校正链从位置1至5依次为GGCGT。
示例2-1、2-2两校正结果相比，位置3的校正结果不同。如前，k值越大，k-mer的特异性越好，校正结果越准确。因此，3-mer是最佳选择。
示例2-3
简单的按单碱基的比对次数进行校正，每个位置选择比对次数最高的作为校正后序列，各位置单碱基的比对次数统计如表3。
表3单碱基比对次数
示例2-3的结果为GGCGT，该结果与示例2-2结果一致，与示例2-1在第3位置存在差异。可见只统计单碱基比对次数无法进行准确的校正。
由上述示例可见，采用本实施例1方法进行校正时，考虑前向序列与当前碱基的组合情况，即考虑最可能存在的连续小片段作为校正结果，有效地避免了只统计单碱基比对次数对于校正结果的偏好性影响。
实施例2
本实施例提供了一种序列的校正方法，与实施例1提供的校正方法大致相同，区别在于：将第一校正序列作为待校正序列，进行迭代处理，迭代2次。
迭代处理具体包括：当进行迭代处理时，将第一校正序列作为待校正序列，再进行一次校正，得到新的第一校正序列。
实施例3
本实施例提供了一种序列的校正方法，与实施例1提供的校正方法大致相同，区别在于还包括有低质量区间校正，区别如下：
将实施例1获得的第一校正序列作为待校正序列，在待校正序列上，根据第一排列结果，将每个最后位点上，获得最高概率的k-mer片段对应的碱基占该位点上k-mer片段总数的比值小于低质量占比预设值80％的位点，标记为低质量位点。
将待校正序列上出现两个以上低质量位点，且两个以上低质量位点之间最大间隔长度小于或等于低质量间隔预设值50碱基的区间划分为低质量区间。
将低质量区间的序列替换为该区间内出现概率最高的测序片段的序列，具体包括：
将多个测序片段排列至第一校正序列上相对应的位置，获得第二排列结果；基于第二排列结果，确定低质量区间内重复次数最多的测序片段，将待校正序列上低质量区间内的序列校正为重复次数最多的测序片段对应的序列，得到第二校正序列。
为了获得更准确的校正序列，在本实施例中，采用优选方案，使用高准确度的测序片段进行校正，高准确度的测序片段由二代测序获得。
实施例4
本实施例提供了一种序列的校正方法，与实施例3提供的校正方法大致相同，区别在于在进行低质量校正时还包括低覆盖校正，区别如下：
将实施例1获得的第一校正序列作为待校正序列，将待校正序列上，对应的测序片段的数量低于覆盖阈值3的区域标记为低覆盖区域，对低覆盖区域进行低覆盖校正，获得第四校正序列。覆盖阈值为位点上对应排列的测序片段的数量。
低覆盖校正包括：
将多个测序片段排列至第二校正序列上相对应的位置，挑选出低覆盖区域对应的排列，获得第四排列结果；
以k-mer方式(k-mer优选为3)对低覆盖区域进行划分，得到多个k-mer片段，基于第四排列结果，对测序片段进行对应的k-mer划分；
基于第四排列结果，确认每个k-mer片段的最后位点对应的测序片段，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率，基于k-mer片段的概率，校正低覆盖区域；正序统计每个最后位点上每一种碱基对应的k-mer片段的概率的计算公式同实施例1中的计算公式。
从反序将低覆盖区域上最末碱基替换为获得最高k-mer片段的概率的碱基，根据最末碱基对应的k-mer片段，确认倒数第二位点的碱基；将倒数第二位点的碱基作为下一个低覆盖k-mer片段的最后位点，确认倒数第三位点的碱基，依次类推，获得第四校正序列。
为了获得更准确的校正序列，在本实施例中，采用优选方案，使用高准确度的测序片段进行校正，高准确度的测序片段由二代测序获得。
本实施例的低覆盖校正可与实施例3的低质量校正并行，实现对实施1获得的第一校正序列的低质量、低覆盖的同时校正。
本实施例的低覆盖校正仅是对待校正序列中低覆盖部分使用较为复杂的算法，而对待校正序列中的低质量部分使用实施例3中计数的简单算法。采用此方法，既能获得准确的校正结果也能提高校正速度。
将低覆盖校正后的结果与实施例3中低质量区间校正后的结果结合，得到第四校正序列。
实施例5
本实施例提供了一种序列的校正方法，与实施例3提供的校正方法大致相同，区别如下：
将第四校正序列作为待校正序列，进行迭代处理，迭代2次。
迭代处理具体包括：当进行迭代处理时，将第四校正序列作为待校正序列，先使用如实施例1方法，然后并行使用实施例3低质量校正、实施例4低覆盖校正，再进行一次校正，共迭代2次，得到新的第四校正序列。
实施例6
本实施例提供了一种序列的校正方法，与实施例3提供的校正方法大致相同，区别在于，在进行低质量校正时还包括低质量位点校正，区别如下：
低质量位点校正包括：
将与前后相邻低质量位点的间隔距离大于低质量间隔预设值50碱基的低质量位点的碱基串联成低质量长序列；
将多个测序片段排列至第二校正序列上相对应的位置，挑选出与低质量长序列对应的排列，获得第三排列结果；
以k-mer方式(k-mer优选为3)对低质量长序列进行划分，得到多个k-mer片段，基于第三排列结果，对测序片段进行对应的k-mer划分；
确认测序片段每个k-mer片段的最后位点，基于第三排列结果，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率(概率的计算同实施例1中提供的计算公式)；
从反序，将低质量长序列上最末碱基替换为获得最高概率的k-mer片段的碱基，根据最末碱基对应的最高概率的k-mer片段，确认倒数第二位点的碱基；将倒数第二位点的碱基作为下一个k-mer片段的最后位点，根据倒数第二位点碱基对应的最高概率的k-mer片段，确认倒数第三位点的碱基，依次类推，使用低质量长序列校正确认后的碱基替换待校正序列中的对应碱基，得到第三校正序列。
为了获得更准确的校正序列，在本实施例中，采用优选方案，使用长读长的测序片段进行校正，长读长的测序片段由三代测序获得。
实施例7
本实施例提供了一种序列的校正方法，与实施例5提供的校正方法大致相同，区别在于：
对实施例5中获得的第四校正序列中低深度区间，使用长读长的测序序列，例如通过三代测序获得的测序序列，采用如实施例4中的低覆盖校正方法进行校正；可选的，进行迭代校正。
验证例
使用公开的测序数据测试本申请技术方案的技术效果。
拟南芥(Arabidopsis thaliana)的测序数据来自于文献：
Michael,T.P.et al.(2018)High contiguity Arabidopsis thaliana genomeassembly with a single nanopore flow cell.Nature communications,9,541；其中包括二代测序数据，及使用PacBio及Nanopre平台获得的三代测序数据。
人(Homo sapiens)的测序数据均来自于文献:
Shi,L.et al.(2016)Long-read sequencing and de novo assembly of aChinese genome.Nature communications,7,12065；其中包括二代测序数据，及使用PacBio平台获得的三代测序数据。
将三代测序数据组装获得的基因序列作为原始的待校正序列。因二代数据准确度高，所以将二代数据作为标准，评价校正结果的准确性。
检测项目：与二代数据相比SNP错误数目、与二代数据相比InDel错误数目、二代数据与校正序列匹配数目、二代数据与校正序列完全匹配数目、校正序列中影响蛋白序列的转录本数目以及CPU时间。其中，CPU时间是进行一轮校正的单次时间。
验证例1
使用拟南芥二代测序数据、OxfordNanopre平台获得的三代测序数据，验证实施例1、3中校正方法的技术效果，具体校正结果见下表4。
使用实施例1提供的校正方法进行一轮校正。校正序列1为利用二代测序数据作为测序片段校正的结果，校正序列2为利用三代测序数据作为测序片段校正的结果。
校正序列3是利用二代测序数据、使用实施例3校正的结果。与校正序列1相比，反映低质量校正的校正效果。
表4校正结果
由表4可知，校正数据的准确度越高，校正结果越好。即便是使用现阶段准确度较差的三代测序数据进行校正，虽然SNP错误数因校正序列本身的准确性问题而增加，但其他指标均得到了改善，特别是匹配Reads数、完全匹配Reads数的增加有利于继续校正时校正数据的利用。
验证例2
验证实施例4、5提供的校正方法的效果，其中实施例4对应迭代次数1，实施例5对应迭代次数2。结果请参照表5和表6。
与现有技术Pilon校正结果进行对比，其参考文献：
Walker,Bruce J.,et al.“Pilon:an integrated tool for comprehensivemicrobial variant detection and genome assembly improvement.”PloS one 9.11(2014):e112963。
表5拟南芥的比对结果
由表5和表6可知，本申请提供的校正方法校正效率高，得到的序列准确度显著优于现有技术。
验证例3
验证实施例5校正方法的校正准确度。
调取拟南芥中的一个数据作为示例，说明技术效果。
待校正序列：
Pilon迭代4次校正结果为：
实施例5校正方法的校正结果为：
通过与二代数据比对，待校正序列存在3处碱基错误(单下划线突出碱基)，Pilon校正了1处，本发明实施例4校正了3处。
对于需要分型的位点(双下划线及加粗突出碱基)，通过调取对应的二代数据，发现仅存在A-C、T-G的组合方式，Pilon校正后是实际不存在的组合方式A-G，本发明的校正方法校正后的组合方式为A-C与实际情况相吻合。说明在校正效果上，本发明的校正方法更为准确。
综上，本申请公开了一种序列的校正方法及其校正装置，涉及生物信息技术领域，具体而言，该校正方法通过对待校正序列进行校正，通过k-mer方式对待校正序列进行划分，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率，从反序，将待校正序列上最后一个位点的碱基替换为获得最高概率的k-mer片段的碱基，根据最后一个位点的碱基对应的最高概率的k-mer片段，确认倒数第二个位点的碱基，依次类推，得到第一校正序列。该校正方法能过在有效的时间内，对待校正序列进行校正，得到精确度更高的基因序列。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
SEQUENCE LISTING
<110> 武汉未来组生物科技有限公司
<120> 一种序列的校正方法及其校正装置
<160> 3
<170> PatentIn version 3.5
<210> 1
<211> 56
<212> DNA
<213> Arabidopsis thaliana
<400> 1
tgagaacgag tagtttggtt gagtattagt gatgatttta aaaacccaaa aatttt 56
<210> 2
<211> 54
<212> DNA
<213> 人工序列
<400> 2
tgagatcgag tagtttggtt gagtattagt gatgatttta aaacccaaaa tttt 54
<210> 3
<211> 55
<212> DNA
<213> 人工序列
<400> 3
tgagatcgag tagtttggtc gagtattact gatgatttta aaaacccaaa atttt 55

Claims (10)

1.一种序列的校正方法，其特征在于，其包括：将多个测序片段与待校正序列进行比对，将多个所述测序片段排列至所述待校正序列上相对应的位置，获得第一排列结果；
以k-mer方式对所述待校正序列进行划分，得到多个k-mer片段；基于所述第一排列结果，对所述测序片段进行对应的k-mer划分；
确认每个所述k-mer片段的最后位点对应的测序片段，基于所述第一排列结果，从正序统计每个所述最后位点上每一种碱基对应的k-mer片段的概率；
从反序，将所述待校正序列上最末碱基替换为获得最高所述概率的k-mer片段的碱基，根据最末碱基对应的最高概率的k-mer片段，确认倒数第二位点的碱基；将所述倒数第二位点的碱基作为下一个k-mer片段的最后位点，根据所述倒数第二位点碱基对应的最高概率的k-mer片段，确认倒数第三位点的碱基，依次类推，得到第一校正序列。
2.根据权利要求1所述的校正方法，其特征在于，将所述第一校正序列作为所述待校正序列，进行迭代校正。
3.根据权利要求1所述的校正方法，其特征在于，所述统计每个所述最后位点上每一种碱基对应的k-mer片段的概率包括：按照计算公式，计算每个所述最后位点上每一种碱基对应的k-mer片段的分值，所述分值越高，所述概率越大，所述计算公式为：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C；
其中，p为所述最后位点在待校正序列上的位置，b为碱基A、T、G、C或缺失，score(p,b)为p位置上碱基b的分值，score(p-1,b)为p-1位置上碱基b的分值；countk-mer为该k-mer对应的碱基组合出现的次数；C为k-mer区域的有效测序深度。
4.根据权利要求3所述的校正方法，其特征在于，所述以k-mer方式对所述待校正序列进行划分为：以k-mer按照预设值为3的方式对所述待校正序列进行划分。
5.根据权利要求1所述的校正方法，其特征在于，所述校正方法还包括低质量校正：
所述低质量校正包括：
在所述待校正序列上，根据第一排列结果，将每个所述最后位点上，获得最高所述概率的k-mer片段对应的碱基占该位点上k-mer片段总数的比值小于低质量占比预设值的位点，标记为低质量位点；
低质量校正包括对所述低质量位点进行区间校正：将所述待校正序列上出现两个以上所述低质量位点，且两个以上所述低质量位点之间最大间隔长度小于或等于低质量间隔预设值的区间划分为低质量区间，将所述低质量区间的序列替换为该区间内出现概率最高的测序片段的序列，得到第二校正序列；
优选地，所述将所述低质量区间的序列替换为该区间内出现概率最高的测序片段的序列包括：
将多个所述测序片段排列至所述第一校正序列上相对应的位置，获得第二排列结果；基于所述第二排列结果，确定所述低质量区间内重复次数最多的测序片段，将所述待校正序列上所述低质量区间内的序列校正为所述重复次数最多的测序片段对应的序列；
优选地，所述低质量占比预设值为80％；
优选地，所述低质量间隔预设值小于等于所述测序片段的长度；
优选地，所述低质量间隔预设值为50碱基。
6.根据权利要求5所述的序列的校正方法，其特征在于，将所述第二校正序列作为所述待校正序列，进行迭代校正。
7.根据权利要求5所述的序列的校正方法，其特征在于，所述低质量校正还包括：
将与前后相邻所述低质量位点的间隔距离大于所述低质量间隔预设值的低质量位点的碱基进行低质量位点校正，获得第三校正序列；
优选地，所述低质量位点校正包括：
将所述与前后相邻所述低质量位点的间隔距离大于所述低质量间隔预设值的低质量位点的碱基串联成低质量长序列；
将多个所述测序片段排列至所述第二校正序列上相对应的位置，挑选出与所述低质量长序列对应的排列，获得第三排列结果；
以k-mer方式对所述低质量长序列进行划分，得到多个k-mer片段，基于所述第三排列结果，对所述测序片段进行对应的k-mer划分；
确认每个所述k-mer片段的最后位点对应的所述测序片段，基于所述第三排列结果，从正序统计每个所述最后位点上每一种碱基对应的k-mer片段的概率；
从反序，将所述低质量长序列上最末碱基替换为获得最高所述概率的k-mer片段的碱基，根据最末碱基对应的最高概率的k-mer片段，确认倒数第二位点的碱基；将所述倒数第二位点的碱基作为下一个k-mer片段的最后位点，根据所述倒数第二位点碱基对应的最高概率的k-mer片段，确认倒数第三位点的碱基，依次类推，使用所述低质量长序列校正确认后的碱基替换待校正序列中的对应碱基，得到第三校正序列；
优选地，所述低质量间隔预设值小于等于所述测序片段的长度：
优选地，所述低质量间隔预设值为50碱基；
优选地，所述统计每个所述最后位点上每一种碱基对应的k-mer片段的概率包括：按照计算公式，计算每个所述最后位点上每一种碱基对应的k-mer片段的分值，所述分值越高，所述概率越大，所述计算公式为：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C
优选地，所述以k-mer方式对所述低质量长序列进行划分为：以k-mer为3对所述低质量长序列进行划分；
优选地，将所述第三校正序列作为所述待校正序列，进行迭代校正。
8.根据权利要求5～7任一项所述的序列的校正方法，其特征在于，所述校正方法还包括：将所述待校正序列上，对应的所述测序片段的数量低于覆盖阈值的区域标记为低覆盖区域，对所述低覆盖区域进行低覆盖校正，获得第四校正序列；
优选的，所述低覆盖校正包括：
将多个所述测序片段排列至所述第二校正序列或第三校正序列上相对应的位置，挑选出所述低覆盖区域对应的排列，获得第四排列结果；
以k-mer方式对所述低覆盖区域进行划分，得到多个k-mer片段，基于所述第四排列结果，对所述测序片段进行对应的k-mer划分；
基于所述第四排列结果，确认每个所述k-mer片段的最后位点对应的所述测序片段，从正序统计每个最后位点上每一种碱基对应的k-mer片段的概率，基于所述k-mer片段的概率，校正所述低覆盖区域；
从反序将所述所述低覆盖区域上最末碱基替换为获得最高所述k-mer片段的概率的碱基，根据最末碱基对应的k-mer片段，确认倒数第二位点的碱基；将倒数第二位点的碱基作为下一个低覆盖k-mer片段的最后位点，确认倒数第三位点的碱基，依次类推，获得第四校正序列；
优选地，所述第四校正序列为将所述低覆盖校正后的结果与如权利要求5中所述的低质量区间校正后的结果的结合；
优选地，所述统计每个所述最后位点上每一种碱基对应的k-mer片段的概率包括：按照计算公式，计算每个所述最后位点上每一种碱基对应的k-mer片段的分值，所述分值越高，所述概率越大，所述计算公式为：
score(p，b)＝max{score(p-1，b∈{A，T，C，G，-})+countk_mer}-C
优选地，所述以k-mer方式对所述低覆盖序列进行划分为：以k-mer为3对所述低覆盖序列进行划分；
优选地，将所述第四校正序列作为所述待校正序列，进行迭代校正。
9.根据权利要求1所述的校正方法，其特征在于，在所述序列比对前，所述校正方法包括将所述多个测序片段进行预处理；
所述预处理包括将所述多个测序片段进行以下一种或多种方式：
1)过滤掉没有跨过校正区域的测序片段；
2)过滤掉比对质量低于比对阈值的测序片段；
3)当所述测序片段只有部分区域能与所述待校正匹配时，切除所述测序片段两端无法匹配的序列，过滤掉当切除的序列长度大于切除预设值的测序片段。
10.一种测序数据的校正装置，其特征在于，所述校正装置用于执行如权利要求1～9任一项所述的校正方法。
CN201910493581.XA 2019-06-06 2019-06-06 一种序列的校正方法及其校正装置 Active CN110246545B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910493581.XA CN110246545B (zh) 2019-06-06 2019-06-06 一种序列的校正方法及其校正装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910493581.XA CN110246545B (zh) 2019-06-06 2019-06-06 一种序列的校正方法及其校正装置

Publications (2)

Publication Number Publication Date
CN110246545A true CN110246545A (zh) 2019-09-17
CN110246545B CN110246545B (zh) 2021-04-13

Family

ID=67886387

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910493581.XA Active CN110246545B (zh) 2019-06-06 2019-06-06 一种序列的校正方法及其校正装置

Country Status (1)

Country Link
CN (1) CN110246545B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111477276A (zh) * 2020-04-02 2020-07-31 上海之江生物科技股份有限公司 微生物的种特异共有序列的获得方法、装置及应用
CN114520024A (zh) * 2022-01-17 2022-05-20 浙江天科高新技术发展有限公司 一种基于k-mer的序列联配方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105574361A (zh) * 2015-11-05 2016-05-11 上海序康医疗科技有限公司 一种检测基因组拷贝数变异的方法
CN106874709A (zh) * 2015-12-12 2017-06-20 北京大学 测序结果中序列数据错误的检测和校正方法
CN107229842A (zh) * 2017-06-02 2017-10-03 肖传乐 一种基于局部图的三代测序序列校正方法
CN108573127A (zh) * 2017-03-14 2018-09-25 深圳华大基因科技服务有限公司 一种核酸第三代测序原始数据的处理方法及其应用
CN108595915A (zh) * 2018-04-16 2018-09-28 北京化工大学 一种基于dna变异检测的三代数据校正方法
WO2019086900A1 (en) * 2017-11-03 2019-05-09 Oxford University Innovation Limited Computer-implemented method and system for determining a disease status of a subject from immune-receptor sequencing data

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105574361A (zh) * 2015-11-05 2016-05-11 上海序康医疗科技有限公司 一种检测基因组拷贝数变异的方法
CN106874709A (zh) * 2015-12-12 2017-06-20 北京大学 测序结果中序列数据错误的检测和校正方法
CN108573127A (zh) * 2017-03-14 2018-09-25 深圳华大基因科技服务有限公司 一种核酸第三代测序原始数据的处理方法及其应用
CN107229842A (zh) * 2017-06-02 2017-10-03 肖传乐 一种基于局部图的三代测序序列校正方法
WO2019086900A1 (en) * 2017-11-03 2019-05-09 Oxford University Innovation Limited Computer-implemented method and system for determining a disease status of a subject from immune-receptor sequencing data
CN108595915A (zh) * 2018-04-16 2018-09-28 北京化工大学 一种基于dna变异检测的三代数据校正方法

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111477276A (zh) * 2020-04-02 2020-07-31 上海之江生物科技股份有限公司 微生物的种特异共有序列的获得方法、装置及应用
CN114520024A (zh) * 2022-01-17 2022-05-20 浙江天科高新技术发展有限公司 一种基于k-mer的序列联配方法
CN114520024B (zh) * 2022-01-17 2024-03-22 浙江天科高新技术发展有限公司 一种基于k-mer的序列联配方法

Also Published As

Publication number Publication date
CN110246545B (zh) 2021-04-13

Similar Documents

Publication Publication Date Title
Bobay et al. Recombination events are concentrated in the spike protein region of Betacoronaviruses
Yang et al. The complete chloroplast genome sequence of date palm (Phoenix dactylifera L.)
Marcet-Houben et al. Beyond the whole-genome duplication: phylogenetic evidence for an ancient interspecies hybridization in the baker's yeast lineage
Barrett et al. Plastid genomes and deep relationships among the commelinid monocot angiosperms
Hein et al. Functional single-cell genomics of human cytomegalovirus infection
Liu et al. Validation of reference genes for gene expression studies in virus-infected Nicotiana benthamiana using quantitative real-time PCR
Karol et al. Complete plastome sequences of Equisetum arvense and Isoetes flaccida: implications for phylogeny and plastid genome evolution of early land plant lineages
Tang et al. Unraveling ancient hexaploidy through multiply-aligned angiosperm gene maps
Liao et al. Unraveling a genetic roadmap for improved taste in the domesticated apple
Eichmeier et al. Comprehensive virus detection using next generation sequencing in grapevine vascular tissues of plants obtained from the wine regions of Bohemia and Moravia (Czech Republic)
CN110246545A (zh) 一种序列的校正方法及其校正装置
KR101930253B1 (ko) 공통서열을 포함한 참조표준 게놈지도 구축 장치 및 방법
Hasiów‐Jaroszewska et al. Ratio of mutated versus wild‐type coat protein sequences in P epino mosaic virus determines the nature and severity of yellowing symptoms on tomato plants
CN102899335A (zh) 一种高通量Small RNA测序获得番木瓜环斑病毒基因组序列的方法
Zhang et al. Comprehensive analysis of differentially expressed genes reveals the molecular response to elevated CO2 levels in two sea buckthorn cultivars
James et al. Universal and taxon-specific trends in protein sequences as a function of age
Zhang et al. Two ancient rounds of polyploidy in rice genome
Nip et al. Reference-free assembly of long-read transcriptome sequencing data with RNA-Bloom2
Casola et al. Pinaceae show elevated rates of gene turnover that are robust to incomplete gene annotation
Guo et al. Miniature inverted-repeat transposable elements drive rapid microRNA diversification in angiosperms
Vogl et al. Probabilistic analysis indicates discordant gene trees in chloroplast evolution
CA3069749A1 (en) Systems and methods for targeted genome editing
Cao et al. Pan-genome analyses of peach and its wild relatives provide insights into the genetics of disease resistance and species adaptation
Csűrös et al. Fast mapping and precise alignment of AB SOLiD color reads to reference DNA
Casey Using Micro-Synteny for Phylogenetic Inference and Analysis

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
CB02 Change of applicant information
CB02 Change of applicant information

Address after: 430000 No.01, 19th floor, building C11, Wuhan software new town phase 2, No.8 Huacheng Avenue, Donghu New Technology Development Zone, Wuhan City, Hubei Province

Applicant after: Wuhan hope group Biotechnology Co.,Ltd.

Address before: 430000 No.01, 19th floor, building C11, Wuhan software new town phase 2, No.8 Huacheng Avenue, Donghu New Technology Development Zone, Wuhan City, Hubei Province

Applicant before: WUHAN NEXTOMICS BIOSCIENCES Co.,Ltd.

GR01 Patent grant
GR01 Patent grant
PE01 Entry into force of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A sequence correction method and its correction device

Effective date of registration: 20210728

Granted publication date: 20210413

Pledgee: Guanggu Branch of Wuhan Rural Commercial Bank Co.,Ltd.

Pledgor: Wuhan hope group Biotechnology Co.,Ltd.

Registration number: Y2021420000070

PC01 Cancellation of the registration of the contract for pledge of patent right
PC01 Cancellation of the registration of the contract for pledge of patent right

Date of cancellation: 20220804

Granted publication date: 20210413

Pledgee: Guanggu Branch of Wuhan Rural Commercial Bank Co.,Ltd.

Pledgor: Wuhan hope group Biotechnology Co.,Ltd.

Registration number: Y2021420000070

PE01 Entry into force of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A Sequence Correction Method and Its Correction Device

Effective date of registration: 20221125

Granted publication date: 20210413

Pledgee: Pudong Shanghai Development Bank Limited by Share Ltd. Wuhan branch

Pledgor: Wuhan hope group Biotechnology Co.,Ltd.

Registration number: Y2022420000372