CN107958138B

CN107958138B - 一种从高通量dna测序的原始信号中读取序列信息的方法

Info

Publication number: CN107958138B
Application number: CN201610899880.XA
Authority: CN
Inventors: 周文雄; 陈子天; 康力; 乔朔; 段海峰; 黄岩谊
Original assignee: Saina Biological Technology (beijing) Co Ltd
Current assignee: Saina biological technology (Beijing) Co., Ltd.
Priority date: 2016-10-14
Filing date: 2016-10-14
Publication date: 2019-06-18
Anticipated expiration: 2036-10-14
Also published as: CN107958138A

Abstract

本发明提供了一种测序结果中序列信息错误校正的方法。本发明利用次级超前校正测序结果中超前量。待测的核酸序列进行测序，检测测序产生的对应于核酸序列的信号；测序结果中，通过次级超前校正该信号。本发明同时考虑初级超前、次级超前和滞后现象，将衰减、失相、整体偏移等问题所造成的信号偏差作为一个整体用于校正测序序列信息。

Description

一种从高通量DNA测序的原始信号中读取序列信息的方法

技术领域

本发明涉及一种从高通量DNA测序的原始信号中读取序列信息的方法；特别是从二代测序的原始信号中读取校正的序列信息的方法，属于基因测序领域。

背景技术

在高通量DNA测序中，理想情况下，每一次测序反应所释放出的原始信号强度与被掺入 DNA新生链的碱基个数成正比。而实际情况中，由于若干原因，该正比关系并不总是成立，例如：1.由于流体冲刷、DNA模板水解、碱基错配等原因，原始信号强度总体上呈衰减趋势；2.由于测序反应不完全、副反应、碱基错配等原因，DNA新生链的长度会随着测序反应的进行而逐渐变得不一致(失相现象)，进而导致原始信号强度发生偏差；3.由于核苷酸自发水解、测序芯片本底荧光等原因，原始信号强度会整体偏高。这些因素导致无法根据正比关系从原始信号强度中直接读出待测DNA的序列信息。

现有从原始测序信号中读取序列信息的方法均只考虑了上述部分原因，例如454的专利仅仅考虑了失相现象，并利用矩阵变换的方法来校正失相造成的信号偏差。而实际上，上述原因同时存在，如果仅仅考虑失相，或简单将失相和衰减、整体偏高等因素剥离开，将影响读取DNA序列信息的准确性。而且454的专利仅仅考虑了失相中的初级超前现象，忽略了次级超前现象，这也影响到了最终结果的准确性。此外，454的专利的实际效果还受许多人为设置的参数的影响，实际使用很不方便。

Ion torrent的专利则试图通过改变核苷酸的进样顺序来缓解上述原因所造成的信号偏差。但首先该方法仅能缓解、而不能真正地校正信号偏差，其次改变后的进样顺序减小了每次测序反应的平均测序读长。

本发明同时考虑造成原始信号发生偏差的上述所有因素，并进行综合校正，从而从原始测序信号中读取准确的DNA序列信息。本发明不影响正常的测序反应流程。本发明包括对单色测序信号和多色测序信号的处理，每种信号的处理包括参数估计和信号校正两部分。

发明内容

本发明涉及一种高通量测序中序列数据误差的校正方法；或者说是高通量DNA测序的原始信号中读取序列信息的方法。

本发明公开了一种高通量测序中序列数据误差的校正方法，其包括以下步骤：

A通过已知参考核酸序列在测序中所产生的核酸序列信号，利用参数估计的方法，获得反应的超前量和滞后量信息；

B对待测的核酸序列进行测序，获得对应于核酸序列的信号；

C利用步骤A的参数估计获得的超前量和滞后量信息，以及步骤B产生的核酸序列的信号，获得次级超前累积量；

D利用步骤B产生的核酸序列的信号和次级超前累积获得相位失配量；

E利用相位失配量修正步骤B产生的核酸序列的信号，推算待测的核酸序列；

F步骤C到E循环，并且用上一轮推算产生的核酸序列信号替代步骤C到E中的核酸序列的信号，直到推算的待测核酸序列收敛为止；

其中，所述的参数估计指的是根据参考核酸序列及其测序信号，推断出超前量、滞后量的方法；

其中所述的次级超前指的是测序中，发生了与该待测的核酸序列位置非期望的延伸，在此非期望延伸的基础上，又继续发生了期望的延伸；

相位失配量指的是，由于超前和滞后导致的测序结果的变化。

参考核酸序列通常也被称为参考序列。

在常见的测序过程中，测序试剂中含有核苷酸底物分子。核苷酸底物分子与待测核酸序列发生反应。在一步反应中，如果发生了超前，则检测到的信号比原本应该收集到的信号偏大。

步骤F中，其推算的待测核酸序列收敛指的是能够获得确定的待测序列。其中所述的收敛也是一般数学意义上的收敛。其可以是数列的形式，也可以是其它的形式。

期望的延伸指的是在测序的过程中，比如一般的利用化学反应进行的测序，当测序试剂中含有能与某个或某几个特定碱基发生反应的成分，则在该待测的核酸序列位置上应该匹配延伸该某个或某几个碱基；也就是正常的测序应该进行的延伸，可以被称为期望的延伸。更简单的说，期望的延伸就是测序应该进行的延伸或者正常的延伸。相反的，如果发生了与该正常的测序应该进行的延伸不相同的延伸，则可以被称为不期望的延伸。

测序中，超前和滞后是常见的现象。次级超前是超前现象的一种。在以前的研究中，超前一词被广泛的应用在基因测序领域中。一般的意义上，发生了不希望的延伸即被称为超前。一般的意义上，期望的延伸没有发生即被称为滞后。

最简单的，当测序试剂中含有与碱基A发生反应的成分，当该待测的核酸序列位置上碱基A发生延伸的时候，是期望的延伸；当其他碱基发生延伸的时候，则是非期望的延伸。一般的测序反应中，只延伸A是正常的延伸。在延伸A的基础上，又发生了非期望的延伸，并且在这个非期望延伸的基础上，又延伸了A，则被称为次级超前。如果在A正常延伸的基础上，仅仅发生了非期望的延伸，则被称为初级超前。某些情况下，这种非期望的延伸可以是测序试剂杂质或者其他测序方法等影响因素造成的。

在某些情况下，比如化学反应进行测序，每次进入的测序试剂中，测序核苷酸底物分子的种类数为2或者3，则次级超前发生的十分频繁，明显影响测序的信号。

根据本发明优选的技术方案，所述参数估计中，根据需要，还包括获得衰减系数、偏移量信息、单位信号信息中的一种或者多种。

其中单位信号信息的获得方式可以有多种方式。待测序列上连接有已知序列的碱基，测序的时候，可以通过该已知序列的碱基的信号获得单位信号。当进行高通量测序的时候，每个采样点或者说反应室内的单位信号可以是不相同的，这并不影响测序反应本身

根据优选的技术方案，所述步骤A中，所述参数估计获得超前和滞后信息指的是，利用参数估计的方法，获得对应于碱基的超前和滞后的常数。对于简单的情况，在化学测序反应中，根据加入的反应液的不同，超前和滞后的常数是不相同的，因此超前和滞后的常数是可以对应于碱基的。

根据优选的技术方案，所述步骤A中，所述参数估计获得超前和滞后信息指的是，利用参数估计的过程，获得包含每轮超前和滞后。

A通过已知参考核酸序列在测序中所产生的核酸序列信号，进行参数估计；

B对待测的核酸序列进行测序，获得对应于核酸序列的信号；

C.利用参数估计获得的超前和滞后信息，以及步骤B产生的核酸序列的信号，获得次级超前累积量；

其中，所述的参数估计指的是根据参考序列及其测序信号，推断出超前、滞后、衰减系数、偏移量的方法；其中所述的次级超前指的是测序中，发生了与该待测的核酸序列位置非期望的延伸，在此非期望延伸的基础上，又继续发生了期望的延伸。相位失配量指的是，由于超前和滞后导致的测序结果的变化。

本发明公开了一种利用次级超前校正测序结果中超前量的方法，其特征在于：

待测的核酸序列进行测序，检测测序产生的对应于核酸序列的信号；测序结果中，通过次级超前校正该信号；所述的次级超前指的是测序中，发生了与该待测的核酸序列位置非期望的延伸，在此非期望延伸的基础上，又继续发生了期望的延伸。

根据本发明优选的实施方式，测序结果中，还包括初级超前量；其中，所述的初级超前量指的是测序中，与核苷酸测序底物不匹配的延伸。简单的说，正常延伸的基础上，仅仅发生了非期望的延伸，则被称为初级超前。

根据本发明优选的实施方式，除了第一次次级超前以外，以后的超前的影响，包括次级超前和初级超前，将累积到以后的测序反应中。

根据本发明优选的实施方式，测序结果中，如果该核苷酸位置获得的信号与单位信号接近，则通过次级超前校正该信号；所述的获得的信号与单位信号接近指的是，反应获得的信号接近单位信号；优选反应获得信号的强度信息同单位信息有60％以内的偏差，优选反应获得信号的强度信息同单位信息有50％以内的偏差，优选反应获得信号的强度信息同单位信息有40％以内的偏差，进一步优选有30％以内的偏差，进一步优选有20％以内的偏差，进一步优选有10％以内的偏差，进一步优选有5％以内的偏差。

根据本发明优选的实施方式，测序中，当获得第n个测序信号的时候，利用其前面的测序信号，通过反馈迭代从已知参开核酸序列的测序数据产生误差的方法获得校正的测序信号；然后再判断该位置是否存在次级超前。

根据本发明优选的实施方式，所述测序指的是常见的化学测序，向待测的核酸序列中加入核苷酸底物分子、酶等测序试剂的反应液进行反应的过程。

根据本发明优选的实施方式，所述测序中，每一次反应加入的核苷酸底物分子可以是一种或两种或三种。

根据本发明优选的实施方式，所述测序，指的是3端开放的测序过程；测序反应加入的核苷酸种类可以是一种或两种或三种。

根据本发明优选的实施方式，所述测序中，反应加入的核苷酸底物分子可以是A、G、C、 T中的一种或多种，或者是A、G、C、U中的一种或多种。

根据本发明优选的实施方式，所述测序中，检测的信号可以是电信号、生物荧光信号、化学荧光信号或者他们的组合。

根据本发明优选的实施方式，在参数估计的过程中，首先根据参考DNA分子的序列推断理想信号h，根据预设的参数，依次计算失相信号s和预测原始测序信号p；计算p和实际原始测序信号f之间的相关系数c。

根据本发明优选的实施方式，使用最优化方法，找到一组参数，使得相关系数c达到最优值；所找到的参数包括超前量、滞后量；或者还包括衰减系数、偏移量、单位信号中的一种或多种。

根据本发明优选的实施方式，所述的超前量滞后量指的是，在测序中，由于超前和滞后导致的失相的程度。

根据本发明优选的实施方式，所述的测序中，将核苷酸分子分为两组，每次测序加入包含其中一组核苷酸分子的测序反应液。两组测序反应液循环加入，进行测序。

根据本发明优选的实施方式，所述的测序中，将参考核酸序列和待测核酸序列同时放入测序；参考核酸序列通过参数估计获得反应的超前量、滞后量、衰减系数、偏移量、单位信号信息；通过参数估计获得的信息，校正待测核酸序列信号，获得校正的核酸序列。

根据本发明优选的实施方式，所述的测序中，待测序列上连接有已知序列的碱基，测序的时候，可以通过该已知序列的碱基的信号获得单位信号。

根据本发明优选的实施方式，每个采样点的单位信号的是不相同的。

本发明公开了一种基因测序系统，包括计算机，其特征在于，利用前面所述的校正从基因测序产生的序列信息误差的方法获得校正的核酸序列。

本领域中，一般的超前指的是，在某个位置，发生了与预期不相同的，非期望的向前的延伸。

本发明所设计到的所有名词，均为基因测序领域的常用含义。

一轮测序指的是对于待测的核酸序列，进行一次测序。

附图说明

以下各图中，正方形代表组成模板DNA的核苷酸，圆形代表组成DNA新生链的核苷酸；画有斜线的图形代表测序引物区域，白色或灰色填充的图形代表不同种类的核苷酸。图1‐3 均为示意图，其并不表示特定的具体核苷酸序列。

图1高通量DNA测序中的失相现象示意图。

图2初级超前和次级超前现象

图3不再发生三级超前

图4参数估计基本过程

图5信号校正基本过程

图6单色2+2测序原始信号

图7单色2+2测序原始信号的参数估计过程中各参数的变化趋势

图8单色2+2测序的原始信号及失相信号

图9单色2+2测序信号的信号校正中的迭代步骤

图10一次双色2+2测序的原始信号

图11双色2+2测序的参数估计中各参数的变化趋势

图12一次双色2+2测序的原始信号及失相信号

图13双色2+2测序的信号校正中的迭代步骤

图14多次单色2+2测序的信号校正的统计结果实施例1(变换矩阵的构造)

具体实施方式

为了进一步阐明本发明，现列出如下具体实施方式。其中所涉及的具体的参数、步骤等，为本领域的常规知识。具体实施方式和实施例并不限制本发明的保护范围。除特殊说明外，本发明涉及到的所有名词均为本领域的常规含义。除特殊说明外，本发明涉及到的所有的基因序列，均为市场上人工合成的序列，例如PCR方法。常见的序列合成的公司有很多，例如 invitrogen。本发明涉及到的所有基因序列，都是invitrogen公司的合成序列。本发明涉及到的所有基因序列，所起到的作用仅仅是对于本发明所述方法的简单说明，并没有特殊含义或限制，简单的替换不影响本发明的实施效果。本发明所涉及到序列、参数、具体步骤可以认为是本领域的常规示例。

本发明典型的测序方法：利用5’多磷酸末端或中间磷酸修饰有具有荧光切换性质荧光团的核苷酸底物分子进行测序；所述的荧光切换性质指的是测序后荧光信号强度相比测序反应前有明显上升；每轮测序使用一套反应液组，每套反应液组至少包括两个反应液，每个反应液包含A、G、C、T核苷酸分子中的至少一种，或者每个反应液包含A、G、C、U核苷酸分子中的至少一种；首先，将待测的核苷酸序列片段固定在反应室中，通入一套反应液组中的一个反应液；检测、记录荧光信息；每次通入一个反应液，将同一反应液组中的其它反应液依次通入，并且每次检测、记录荧光信息；其中，所述反应液组中，至少有一个反应液包含两种或者三种不同的核苷酸分子。

高通量测序通过实施一系列酶促反应、并检测反应所释放出的信号来获得待测DNA的序列信息。若某一DNA新生链已延伸至第n个碱基，当前酶促反应中所加入的核苷酸恰好和待测DNA模板上的第n+1至第n+m个碱基互补配对，则理想情况下，在该次酶促反应中该DNA 新生链将延伸至第n+m个碱基。如果该DNA新生链在该次酶促反应中实际延伸超过了第n+m 个碱基，则称该DNA新生链在该次反应中发生超前现象；如果该DNA新生链在该次酶促反应中实际延伸不足第n+m个碱基，则称该DNA新生链在该次酶促反应中发生滞后现象。超前现象和滞后现象合称失相现象。需要说明的是，该DNA新生链在延伸至第n个碱基时，可能已经发生了多次超前和滞后现象。

如图1示意图所示，在测序反应前，所有的DNA新生链具有相同的长度1。在测序反应后，DNA分子1、3和5被正常延伸，长度为2；DNA分子2由于副反应而发生超前现象，长度为3；DNA分子4因反应不完全而未被延伸，发生滞后现象，长度为1。在测序反应后，各 DNA新生链的长度发生差异。图1中所画5个DNA分子只是简略表示，不代表实际测序中也只有5个DNA分子(实际测序中有多个DNA分子)。

如图2示意图所示，在某个共聚物a被正常延伸后，副反应造成紧随共聚物a之后第一个共聚物b被延伸，该现象称为初级超前现象。如果共聚物b仅有1个核苷酸组成，则紧随b之后的共聚物c也会被进一步延伸，该现象称为次级超前现象。如果共聚物b由不止一个核苷酸组成，则次级超前现象不会发生。

i.测序方法

本发明采用如下方法对DNA进行测序。本发明所涉及的基因测序方法可以参考CN2015108223619。将待测DNA固定在固相表面，杂交上测序引物，不断实施测序反应并检测反应所释放的信号。每一次反应包括如下步骤：向反应器(芯片)加入含有核苷酸、酶等反应所必需试剂的反应液，发生特定生化反应，检测反应所释放的信号，清洗反应器。所加入的核苷酸可以是天然的脱氧核苷酸，也可以是带有化学修饰基团的核苷酸，但其3’端为羟基。每一次反应所加入的核苷酸种类可以为1种、2种或3种，但不能是4种(4种指的是 ACGT或ACGU)。相邻两次反应所加入核苷酸种类的并集包括全部4种核苷酸。

某一次反应若加入2种核苷酸，则反应中这2种核苷酸可以释放出相同类型的信号，也可以释放出不同类型的信号；某一次反应若加入3种核苷酸，则这3种核苷酸可以释放出相同类型的信号，也可以释放出互不相同类型的信号，也可以其中2种释放出相同类型的信号、另1种释放出与之不同类型的信号。这里信号的类型指的是信号的形式(如电信号、生物荧光信号、化学荧光信号等)，或光学信号的颜色(如绿色荧光信号、红色荧光信号等)，或以上的混合。这里为了简便起见，凡是某一次反应中所有核苷酸所释放信号类型全部相同的，称为单色信号；凡是一次反应中所有核苷酸所释放的类型不止一种的，称为多色信号。这里的“色”只是为了简便起见，并不把信号的类型限定在光学信号中。

本发明涉及三种含义不同的信号，分别是：

1.理想信号h，指根据待测DNA的序列及加入核苷酸的顺序，在理想情况下直接推断出的测序信号，直接反映了DNA的序列信息；

2.失相信号s，指理想信号h受到失相现象而产生偏差后形成的信号；

3.预测原始测序信号p，指根据预设的参数，失相信号s在考虑了被延伸碱基数目和测序信号强度的倍比关系(单位信号)、信号衰减、整体偏移等因素后形成的信号，是根据预设的参数、对实际原始测序信号的预测；

4.实际原始测序信号f，指高通量DNA测序中仪器直接测量得到的信号。

ii.参数估计

根据已知序列的参考DNA分子及其实际原始测序信号，推断出本次测序中相关参数的过程，称为参数估计。参数估计的基本过程如图4所示。参数估计涉及一组可以描述本次测序相关性质的参数，例如失相系数、单位信号强度、衰减系数、整体偏移系数等。

首先根据参考DNA分子的序列推断理想信号h，根据预设的参数，依次计算失相信号s 和预测原始测序信号p。计算p和实际原始测序信号f之间的相关系数c。使用最优化方法，找到一组参数，使得相关系数c达到最优值。这里的相关系数c包括但不限于皮尔逊相关系数、斯皮尔曼相关系数、平均互信息、欧几里得距离、汉明距离、车比雪夫距离、马哈兰诺比斯距离、曼哈顿距离、明科斯基距离、对应信号差值的绝对值的最大值或最小值等。这里的最优化方法包括但不限于网格搜索法、穷举法、梯度下降法、牛顿法、Hessian矩阵法、启发式搜索等，其中启发式搜索包括但不限于遗传算法、模拟退火算法、蚁群算法、谐和算法、火花算法、粒子群算法、免疫算法等。这里提到的相关系数和最优化方法均为数学中的常规知识。

通过简单的，超前、滞后、偏移量等参数，对于测序信号的影响，就可以得到理想信号和实际测序信号之间的变换。而通过理想信号和实际测序信号之间的拟合的过程，也可以得到超前、滞后、偏移量等参数，这也就是参数估计的过程。对于该拟合过程中方法可以如上一段所述。具体拟合过程的表现形式，可以是矩阵的形式，也可以是函数的形式。

如果测序中采集到的是单色信号，则直接按上述方法计算。如果测序中采集到的是多色信号，则将每种类型的信号单独拆分出来，每种类型的信号单独按照上述方法计算。

利用h计算s的实施方式是根据h的特征及有关参数，构造变换矩阵T，并利用T将h变换为s。利用s计算p的实施方式是根据有关参数，构建变换函数φ，并利用d将s变换为 p。二者的具体实施方式将在下面详述。

iii.信号校正

根据参数估计所得到的参数，以及未知序列的待测DNA的实际原始测序信号，推断出待测DNA序列信息的过程，称为信号校正。信号校正的基本过程如图5所示，大体上可以看作参数估计的逆过程。

首先根据参数估计得到的参数，利用变换函数φ的反函数将实际原始测序信号f变换为失相信号s。将s视为零阶失相信号s₀，根据s₀和有关参数构建变换矩阵T₁，并利用T₁的广义逆矩阵将s₀变换为一阶失相信号s₁；再根据s₁和有关参数构建变换矩阵T₂，并利用T₂的广义逆矩阵将s₁变换为二阶失相信号s₂；以此类推，计算出一系列失相信号s₀,s₁,s₂,…。若计算中发现两个相邻失相信号s_i和s_i+1相等，则停止计算，并返回s_i作为信号校正的结果。

上述的广义逆矩阵也可以用吉洪诺夫正则化(Tikhonov regularization)的方法代替。

上述利用变换函数φ的反函数将f变换为s的过程，及利用T的广义逆矩阵将s_i变换为 s_i+1的过程将在下面详述。

iv.变换矩阵T的构造方法

变换矩阵T的构造依赖于一条测序有关的信号x及与失相参数。在参数估计中，信号x 是理想信号h；在信号校正中，信号x是各阶失相信号s_i。为了提高校正的准确性，可以通过在信号x后添加若干个1来延长信号x；根据本发明优选的方法，通常添加1‐100个1。失相参数包括超前系数ε和滞后系数λ。

变换矩阵T的构造中还需要一个辅助矩阵D。设信号x由m个数值组成，测序反应实际进行了n次，则变换矩阵T和辅助矩阵D均是n行m列的矩阵。辅助矩阵D的第一行中，只有第一列的元素为1，其他元素均为0。

利用辅助矩阵D的第k行来计算变换矩阵T的第k行。对变换矩阵T第k行的第1个元素：

1.若k为奇数，则考虑到滞后现象，令该元素为(1‐λ)D_1i；

2.若k为偶数，则令该元素为0。

对变换矩阵T第k行的第i个元素(第1个元素除外)：

1.若k和i的奇偶性相同，则考虑到滞后现象，令该元素为(1‐λ)D_ki；

2.若k和i的奇偶性不同，则考虑到初级超前现象，令该元素为ε(1‐λ)D_k,i‐1；

3.若信号x的第i‐1个元素小于2，则考虑到次级超前现象，该元素在1和2的计算结果的基础上，还要再加上变换矩阵T同一行的第i‐1个元素T_k,i‐1。

利用变换矩阵T的第k行来计算辅助矩阵的第k+1行。辅助矩阵D的第1行中，只有第1列的元素为1，其他元素均为0。对辅助矩阵的第k行(第1行除外)：

1.第1个元素为辅助矩阵上一行、同一列的元素D_k‐1,i和变换矩阵T中对应元素的上一行、同一列元素T_k‐1,i的差值；

2.第i个元素(第1个元素除外)在辅助矩阵上一行、同一列的元素D_k‐1,i和变换矩阵T中对应元素的上一行、同一列元素T_k‐1,i的差值的基础上，再加上变换矩阵T 中对应元素的上一行、上一列元素T_k‐1,i‐1。

因此，本发明先规定辅助矩阵D的第1行的值，然后根据辅助矩阵D的第1行去计算变换矩阵T的第1行；利用变换矩阵T的第1行去计算辅助矩阵的第2行；利用辅助矩阵D的第2行去计算变换矩阵T的第2行；以此类推，逐步得到辅助矩阵和变换矩阵T的所有元素的值。

辅助矩阵D只是为了计算上的简便而引入的，可以通过常规的数学变形将其消去，从而直接计算变换矩阵T。

在上述计算中，失相参数与核苷酸类型有关，也和被计算的元素所处的行号k和列号i 有关。在实际计算中，既可以为简便起见而使失相系数ε和λ分别保持恒定，也可以为精确起见而使失相系数ε和λ随核苷酸的种类、行号k和列号i而有所变化。

在参数估计中，根据预设的失相系数和理想信号h，按照上述计算方式得到变换矩阵T，则失相信号s为变换矩阵T和理想信号h的积。若理想信号h表示为一个列向量，则s为T乘以h；若理想信号表示为一个行向量，则s为h乘以T的转置矩阵。

在信号校正中，根据预设的失相系数和第i阶失相信号s_i，按照上述计算方式得到变换矩阵T，则第i+1阶失相信号为变换矩阵T的广义逆矩阵T⁺和第i阶失相信号的积。若s_i表示为一个列向量，则s_i+1为T⁺乘以s_i；若s_i表示为一个行向量，则s_i+1为s_i乘以T⁺的转置矩阵。第 i+1阶失相信号si+1在按上述方法计算完毕后，可以再接着取整数值，取整的方式包括但不限于：

1.四舍五入：取为最接近的整数值；

2.向上取整：取为大于自身的最小整数；

3.向下取整：取为小于自身的最大整数；

4.向0取整：若自身大于0，则向下取整；若自身小于0，则向上取整；

5.正取整：按上述任何一种方式取整，然后将所有的非正数改为1。

v.变换函数的构造方法

变换函数φ依赖于若干参数，包括单位信号a(被延伸碱基数目和测序信号强度的倍比关系)、衰减系数b、整体偏移c等。这里的参数a、b、c可以分别是单一系数，也可以是一组系数。例如单位信号a和核苷酸类型、测序反应发生的次数有关。在计算中既可以为简便起见而使这些参数使用单一值，也可以为精确起见而使这些参数随相关因素而发生变化，还可以某些参数使用单一值、某些参数随相关因素发生变化。

变换函数φ(s)的形式包括但不限于：

1.φ(s)＝φ_aφ_bφ_s+φ_c

2.φ(s)＝φ_aφ_b(φ_s+φ_c)

3.φ(s)＝φ_b(φ_aφ_s+φ_c)

4.φ(s)＝φ_a(φ_bφ_s+φ_c)

其中φ_a、φ_b、φ_c和φ_s均分别为和a、b、c有关的数学函数，包括但不限于常函数、幂函数、指数函数、对数函数、三角函数、反三角函数、取整函数、特殊函数，以及上述函数相互运算、复合、迭代、分段所产生的函数等。其中特殊函数包括但不限于椭圆函数、伽马函数、贝塞尔函数、贝塔函数等。

变换函数φ(x)将失相信号s变换为预测原始测序信号p，即p＝φ(s)。变换函数φ(x) 的反函数φ^-1(x)将实际原始测序信号f变换为失相信号s，即s＝φ^-1(f)。这里的反函数取数学中的常规含义。

相比现有方法(主要是454的专利US8364417B2)，本发明主要做了几点改进。

第一，同时考虑初级超前、次级超前和滞后现象来构建变换矩阵，并利用该变换矩阵来校正因失相造成的测序错误。

第二，将衰减、失相、整体偏移等问题所造成的信号偏差作为一个整体来解决，既不是只校正一个问题所造成的信号偏差，也不是简单地一个一个地解决。

第三，改进了信号校正的方法，避免引入需要人为主观因素判断的参数设置，提高了方法的稳健性和可重复性。

第四，既有单色信号的校正，也有双色信号的校正。

本发明的方法具有以下效果，相比于背景技术提到的方法，具有以下优点：

1.在2+2式测序中，次级超前现象非常显著，所造成的信号偏差是没有考虑次级超前现象的454专利所无法校正的。本发明考虑了次级超前现象，可以很好地校正该现象所造成的信号偏差。

2.在实际运用中，如果只用简单的线性拟合方法来从原始测序信号中读出序列信息，将至多准确读到大约100bp左右。而如果对相同的数据采用本发明所描述的方法，将能准确读取到350bp左右，极大地提高测序读长和测序准确率。

3.本发明既能校正单色信号，也能校正双色信号。

4.本发明不影响正常的测序进样顺序。这是ion torrent测序方法不能达到的。

实施例1

采用2+2式测序方法，组合为M/K，即凡奇数轮加入A或C，凡偶数轮加入G或T。当被测DNA序列为CCTGTATGACCGTATTCCGGGTCCTGTCGGTA时，所获得的理想信号为h＝(2,3,1, 2,3,2,1,2,2,4,2,3,1,3,1)。

为计算简便起见，可以在计算中认为M和K的超前系数相同，滞后系数也相同。例如，当超前系数为0.02、滞后系数为0.01、共进行10次测序反应时，根据前述方法构造出的变换矩阵为：

为提高计算准确性起见，可以在计算中认为M和K的超前系数不同，滞后系数也不同。例如，当M的超前系数和滞后系数分别为0.02和0.01，K的超前系数和滞后系数分别为0.01 和0.02，共进行10次测序反应时，根据前述方法构造出的变换矩阵为：

若采用2+2双色测序方法，则变换矩阵的计算方法不变，不同的只是在参数估计和信号校正中的使用方式。

实施例2(单色2+2的参数估计)

一次单色2+2测序实验，核苷酸组合为M/K，即凡奇数轮加入A或C，凡偶数轮加入G或T。被测序列为：

AAGAGCTGGACAGCGATACCTGGCAGGCGGAGCTGCATATCGAAGTTTTCCTGCCTGCTCAGGTGCCGGATTCAGAGCTGGATGCGTGGATGGAGTCCCGGATTTATCCGGTGATGAGCGATATCCCGGCACTGTCAGATT TGATCACCAGTATGGTGGCCAGCGGCTATGACTACCGGCGCGACGATGATGCGGGCTTGTGGAGTTCAGCCG ATCTGACTTATGTCATTACCTATGAAATGTGAGGACGCTATGCCTGTACCAAATCCTACAATGCCGGTGAAAGGT GCCGGGACCACCCTGTGGGTTTATAAGGGGAGCGGTGACCCTTACGCGAATCCGCTTTCAGACGTTGACTGG TCGCGTCTGGCAAAAGTTAAAGACCTGACGCCCGGCGAACTGACCGCTGAGTCCTATGACGACAG

共进行200次测序反应，得到实际原始测序信号如图6所示。

可以看到：原始测序信号的取值范围大约在100～1500之间，整体呈下降趋势，从大约第 80次测序反应开始，信号呈交替波动状，无法从中直接读取序列信息。

利用前述参数估计方法，根据被测DNA分子的序列及测序方式，可推知理想信号为h＝(2, 1,1,1,1,3,3,1,1,1,1,1,3,3,2,2,1,2,1,1,1,2,2,1,1,1,1,1,2,5,2,2,2,2,1,1,2,4,2,2,1,2, 2,1,1,1,1,3,1,2,1,4,1,3,1,2,3,2,1,3,1,1,2,4,1,2,1,1,1,1,1,1,1,1,3,2,3,3,2,1,1,4,1, 1,5,2,1,6,3,1,1,2,1,1,1,2,2,1,3,2,1,1,1,1,2,1,1,2,1,2,1,3,1,6,1,3,2,1,2,1,1,1,1,2, 2,2,1,3,2,2,3,1,1,2,3,4,1,2,2,1,1,1,1,2,2,3,6,1,2,1,4,2,2,4,3,4,2,3,7,9,1,1,2,4,1, 1,1,4,4,2,2,1,1,1,2,1,2,1,1,3,2,1,2,4,2,4,1,1,1,2,1,3,5,3,3,1,3,2,2,1,3,2,1,1,3,2, 3,1,1,2,1,2,2,1,1,2,2,1,3,1)。

按照前述参数估计的方法来估计此次测序中的相关参数。构造变换矩阵时，为计算准确度起见，认为M和K的超前和滞后系数均不相同。

设t为测序反应的次数。构造变换函数φ(s)＝φ_aφ_bφ_s+φ_c，其中：

1.其中a称为单位信号；

2.其中b称为衰减系数；

3.其中d和e分别称为M和K的整体偏移；

4.其中s为失相信号。

参数估计中，所使用的相关系数为皮尔逊相关系数，所使用的最优化方法为梯度下降。在经过48轮迭代计算后，梯度下降达到收敛条件，得到M的超前系数为0.0117，M的滞后系数为0.0067，K的超前系数为0.0128，K的滞后系数为0.0067，单位信号为519.7，衰减系数为0.9849，M的整体偏移为122.7，K的整体偏移为150.1，相关系数为0.999961。所有参数在迭代计算过程中的变化趋势如图7所示。

实施例3(单色2+2的信号校正)

一次单色2+2测序实验，被测序列未知。其实际原始测序信号f，以及经实施例1中的变换函数φ(s)的反函数和有关参数变换得到的失相信号如图8所示(倒三角符号表示该位置上的信号强度与理想信号不符)。

可见在经过变换函数φ(s)的反函数变换得到的失相信号中依然有许多位置上的信号值与理想信号不符。经过前述信号校正的步骤，共进行4次迭代，分别得到一阶失相信号s₁、二阶失相信号s₂、三阶失相信号s₃和四阶失相信号s₄。在经过四舍五入取整后，s₃和s₄的所有信号值均相等，因此停止迭代，输出s₄作为校正结果。这四阶失相信号如图9所示，其中倒三角表示该位置上的信号强度与理想信号不符。可以看到，随着迭代的进行，倒三角符号逐渐变少，表明校正的精度越来越高，最终输出的校正结果中，前173次测序反应的信号均被校正至完全正确，从第174次反应起才出现校正错误。

实施例4(双色2+2的参数估计)

一次双色2+2测序实验，核苷酸组合为M/K，其中A和G标记相同颜色的荧光基团，C和T标记相同颜色的荧光基团，被测序列为：

共进行200次测序反应，得到实际原始测序信号如图10所示。

可以看到：原始测序信号的取值范围大约在100～1200之间，整体呈下降趋势，从大约第 80次测序反应开始，信号呈交替波动状，无法从中直接读取序列信息。

由于采用了双色测序方法，因此理想信号、失相信号、原始测序信号等均分别有2条，分别对应AG标记的荧光基团和CT标记的荧光基团。

利用前述参数估计方法，根据被测DNA分子的序列及测序方式，可推知AG标记的荧光基团所对应的理想信号为：

h₁＝(2,1,1,1,0,2,2,1,0,1,1,0,1,2,1,2,0,2,1,1,0,1,1,0,1,0,0,1,2,1,0,1,0,1,0,0,1, 3,0,2,1,0,1,1,1,1,0,2,1,1,0,3,1,2,1,1,0,2,1,0,1,0,0,3,1,1,1,1,0,1,1,0,1,0,0,2,1,1, 1,1,1,1,1,0,2,1,1,4,1,1,0,2,0,0,1,1,1,0,1,2,0,1,0,1,1,1,1,1,1,1,0,3,0,3,1,1,1,1,0, 1,1,0,0,1,1,0,1,1,1,0,1,0,1,1,3,2,1,2,1,1,0,0,1,1,0,1,4,0,0,0,3,1,0,3,3,3,0,3,2,4, 1,0,2,4,1,1,0,3,1,0,1,1,0,1,2,0,0,1,0,0,1,1,1,2,1,2,0,1,0,1,0,2,4,1,3,1,1,1,1,1)。

CT标记的荧光基团所对应的理想信号为：

h₂＝(0,0,0,0,1,1,1,0,1,0,0,1,2,1,1,0,1,0,0,0,1,1,1,1,0,1,1,0,0,4,2,1,2,1,1,1,1, 1,2,0,0,2,1,0,0,0,1,1,0,1,1,1,0,1,0,1,3,0,0,3,0,1,2,1,0,1,0,0,1,0,0,1,0,1,3,0,2,2, 1,0,0,3,0,1,3,1,0,2,2,0,1,0,1,1,0,1,1,1,2,0,1,0,1,0,1,0,0,1,0,1,1,0,1,3,0,2,1,0,2, 0,0,1,1,1,1,2,0,2,1,2,2,1,0,1,0,2,0,0,1,0,1,1,0,1,2,2,2,1,2,1,1,1,2,1,0,1,2,0,5,5, 0,1,0,0,0,0,1,1,3,2,1,0,1,0,0,1,2,0,1,3,1,0,1,2,1,2,1,0,1,1,1,1,1,2,0,0,2,1,1,0)。

按照前述参数估计的方法来估计此次测序中的相关参数。构造变换矩阵时，为计算准确度起见，认为M和K的超前和滞后系数均不相同。对于某个根据给定的失相系数构造出的变换矩阵T，认为AG标记的荧光基团所对应的失相信号为s₁＝Th₁，CT标记的荧光基团所对应的失相信号为s₂＝Th₂。

设t为测序反应的次数。针对AG标记的荧光基团和CT标记的荧光基团，分别构造变换函数φ1(s)＝φ_a1φ_bφ_s+φ_c1和φ2(s)＝φ_a2φ_bφ_s+φ_c2，其中

1.其中a₁和a₂分别是AG和CT标记的荧光基团所释放信号的单位信号；

2.其中b称为衰减系数；

3.其中d₁、e₁、d₂、e₂分别是A、G、C、T的整体偏移；

4.其中s为失相信号。

参数估计中，所使用的相关系数为皮尔逊相关系数，所使用的最优化方法为梯度下降。在经过17轮迭代计算后，梯度下降达到收敛条件，得到M的超前系数为0.0125，M的滞后系数为0.0067，K的超前系数为0.0126，K的滞后系数为0.0068，AG和CT标记的荧光基团所释放信号的单位信号分别为519.8和480.7，衰减系数为0.9860，A的整体偏移为164.5，G的整体偏移为133.2，C的整体偏移为140.7，T的整体偏移为175.7，相关系数为0.999964。所有参数在迭代计算过程中的变化趋势如图11所示。

实施例5(双色2+2的信号校正)

一次双色2+2测序实验，凡奇数轮加入G和T，凡偶数轮加入A和C，其中A和G标记相同颜色的荧光基团，C和T标记相同的另一种颜色的荧光基团。被测序列未知。本次测序中得到的原始测序信号f，以及经实施例4中变换函数φ1(s)和φ2(s)的反函数和有关参数变换得到的失相信号s如图12所示。由于采用了双色测序方法，因此理想信号、失相信号、原始测序信号等均分别有2条，分别对应AG标记的荧光基团和CT标记的荧光基团。图12中大量的倒三角符号表明在失相信号s中依然有许多位置上的信号值与理想信号不符。

经过前述信号校正的步骤，共进行4次迭代，分别得到一阶失相信号s₁、二阶失相信号 s₂、三阶失相信号s₃和四阶失相信号s₄。在经过四舍五入取整后，s₃和s₄的所有信号值均相等，因此停止迭代，输出s₄作为校正结果。这四阶失相信号如图13所示，其中倒三角表示该位置上的信号强度与理想信号不符。可以看到，随着迭代的进行，倒三角符号逐渐变少，表明校正的精度越来越高，最终输出的校正结果中，前166次测序反应的信号均被校正至完全正确，从第167次反应起才出现校正错误。

实施例6(大量序列得出的综合性能)

为综合评估本发明从原始测序信号中读取序列信息的准确度，分别进行了五次单色2+2 测序实验，每次测序均进行500次测序反应。每次测序实验中，一部分被测DNA被作为参照，其序列和原始测序信号被用于参数估计；另一部分被测DNA被作为测序样品，将分别使用两种方式进行信号校正：一种根据本发明所描述的方法，利用参照DNA所估计出的参数对其进行信号校正，另一种简单地假设原始测序信号和理想信号间存在简单的正比关系，以此推断 DNA序列信息。

这五次测序实验中，利用参照DNA的原始测序信号所估计出的失相系数分别为0.001、 0.003、0.005、0.010和0.011(参数估计时设置超前系数和滞后系数相等)。对于信号校正，分别统计两种方法校正得到的信号中信号强度和理想信号强度不符的第一次测序反应的编号 (即完全正确的校正信号的长度)，并绘制成柱状图(如图14所示，误差线为标准差)。可以看到，当失相系数为0.001时，依简单正比关系推算所得到的校正信号在不到100次测序反应时即出现校正错误，而本发明所描述的方法得到了完全正确的校正结果。随着失相系数的增大，两种方法的校正结果的准确率均有所下降，但本发明的校正结果中完全正确的校正信号的长度依然是通过简单正比关系推算的3‐5倍，体现了本发明在提高从原始测序信号中读取DNA序列的准确性和有效读长上的明显优越性。

实施例7

2+2测序，单色：配置3套反应液，每套两瓶，每瓶有两种标记有荧光基团的碱基，荧光基团均为X。一套中的两瓶反应液，恰好包含完整的4种碱基。6瓶溶液互不重复。

	第一瓶	第二瓶
			第一套	AX+CX	GX+TX
第二套	AX+GX	CX+TX
			第三套	AX+TX	CX+GX

完整的测序过程包括三轮，三轮依次进行。每轮的测序过程分别使用上述三套试剂。除此之外完全相同(使用相同的测序引物，反应条件完全相同)。

每轮测序包含：

1.将测序引物杂交在已经制备好的DNA阵列上

2.开始测序过程。重复2.1‐2.4过程有限次。

2.1进第一瓶试剂。反应并采集荧光信号。

2.2清洗flowcell中的全部残留反应液和产生的荧光分子

2.3进第二瓶试剂。反应并采集荧光信号。

2.4 清洗flowcell中的全部残留反应液和产生的荧光分子

3.将延伸过的测序引物解旋。

至此，便可进行下一轮实验。

准备反应液：配制测序反应液洗液，简称洗液，含有:

20mM Tris‐HCl pH 8.8

10mM(NH4)2SO4

50mM KCl

2mM MgSO4

0.1％20

配制测序反应液母液(简称母液)，含有:

20mM Tris‐HCl pH 8.8

10mM(NH4)2SO4

50mM KCl

2mM MgSO4

0.1％20

8000unit/mL Bst polymerase

100unit/mL CIP

配制三组测序反应液，共六瓶。分别为：

1A、母液+20uM dA4P‐TG+20uM dC4P‐TG

1B、母液+20uM dG4P‐TG+20uM dG4P‐TG

2A、母液+20uM dA4P‐TG+20uM dG4P‐TG

2B、母液+20uM dC4P‐TG+20uM dG4P‐TG

3A、母液+20uM dA4P‐TG+20uM dT4P‐TG

3B、母液+20uM dC4P‐TG+20uM dG4P‐TG

配制好的反应液和母液，置于4c冰箱或冰上待用。

杂交测序引物：

将测序芯片内注入测序引物溶液(10uM溶解于1X SSC buffer)，升温至90度，在以5 摄氏度/min的速度降温至40度。用洗液冲洗掉测序引物溶液。

进行第一次测序：

将测序芯片置于测序仪上。使用第一组反应液进行测序。遵循如下流程。

1，通入洗液10mL，冲洗芯片

2，将芯片降温至4摄氏度

3，通入100uL反应液1A

4，将芯片升温至65摄氏度

5，等待1min

6，用473nm激光激发，拍摄荧光图像。

7，通入洗液10mL，冲洗芯片

8，将芯片降温至4摄氏度

9，通入100uL反应液1B

10，将芯片升温至65摄氏度

11，等待1min

12，用473nm激光激发，拍摄荧光图像。

重复1‐12的步骤50次，得到100个荧光信号。

Claims

1.一种校正从基因测序产生的序列信息误差的方法，其包括：

A通过已知参考核酸序列在测序中所产生的核酸序列信号，利用参数估计的方法，获得反应的超前和滞后信息；

B对待测的核酸序列进行测序，获得对应于核酸序列的信号；

C利用步骤A的参数估计获得的超前和滞后信息，以及步骤B产生的核酸序列的信号，获得次级超前累积量；

E利用相位失配量修正步骤B产生的核酸序列的信号，推算待测的核酸序列信号；

F步骤C到E循环，并且用上一轮推算产生的核酸序列信号替代步骤C到E中的核酸序列的信号，直到推算的待测核酸序列信号收敛为止，

其中，所述的参数估计指的是根据参考核酸序列及其测序信号，推断出超前、滞后的方法；

2.根据权利要求1所述的方法，其特征在于，

所述步骤A的参数估计中，还包括获得衰减系数信息。

3.根据权利要求1所述的方法，其特征在于，

所述步骤A的参数估计中，包括获得偏移量信息。

4.根据权利要求1所述的方法，其特征在于，

所述步骤A的参数估计中，包括获得单位信号信息。

5.根据权利要求1所述的方法，其特征在于，

所述步骤A中，所述参数估计获得超前和滞后信息指的是，利用参数估计的方法，获得对应于碱基的超前和滞后的常数。

6.根据权利要求1所述的方法，其特征在于，

步骤A中，所述参数估计获得超前和滞后信息指的是，利用参数估计的方法，获得包含每轮的超前和滞后信息。

7.根据权利要求1所述的方法，其特征在于，

在参数估计的过程中，首先根据参考DNA分子的序列推断理想信号h，根据预设的参数，依次计算失相信号s和预测原始测序信号p；计算p和实际原始测序信号f之间的相关系数c。

8.一种校正从基因测序产生的序列信息误差的方法，其包括：

B对待测的核酸序列进行测序，获得对应于核酸序列的信号；

C利用参数估计获得的超前和滞后信息，以及步骤B产生的核酸序列的信号，获得次级超前累积量；

其中，所述的参数估计指的是根据参考序列及其测序信号，推断出超前、滞后、衰减系数、偏移量的方法；

9.根据权利要求8所述的方法，其特征在于，

10.一种利用次级超前校正测序结果中超前量的方法，其特征在于，

待测的核酸序列进行测序，检测测序产生的对应于核酸序列的信号；

通过次级超前校正该信号；

所述的次级超前指的是测序中，发生了与该待测的核酸序列位置非期望的延伸，在此非期望延伸的基础上，又继续发生了期望的延伸。

11.根据权利要求10所述的方法，其特征在于，

测序结果中，还包括初级超前；

其中，所述的初级超前指的是测序中，与核苷酸测序底物不匹配的延伸。

12.根据权利要求10所述的方法，其特征在于,

测序结果中，如果该核酸序列位置获得的信号与单位信号接近，则通过次级超前校正该信号；

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有60％以内的偏差。

13.根据权利要求12所述的方法，其特征在于，

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有50％以内的偏差。

14.根据权利要求12所述的方法，其特征在于，

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有40％以内的偏差。

15.根据权利要求12所述的方法，其特征在于，

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有30％以内的偏差。

16.根据权利要求12所述的方法，其特征在于，

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有20％以内的偏差。

17.根据权利要求12所述的方法，其特征在于，

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有10％以内的偏差。

18.根据权利要求12所述的方法，其特征在于，

所述的获得的信号与单位信号接近指的是，反应获得信号的强度信息同单位信息有5％以内的偏差。

19.根据权利要求10所述的方法，其特征在于，

测序中，当获得第n个测序信号的时候，利用其前面的测序信号，通过反馈迭代从已知参考核酸序列的测序数据产生误差的方法获得校正的测序信号；然后再判断该位置是否存在次级超前。

20.根据权利要求10-18任一项所述的方法，其特征在于，

所述测序中，每一次测序反应试剂中加入的核苷酸底物分子可以是一种或两种或三种。

21.根据权利要求10-18任一项所述的方法，其特征在于，

所述测序，指的是3端开放的测序过程；测序反应加入的核苷酸种类可以是一种或两种或三种。

22.根据权利要求10所述的方法，其特征在于，

所述测序中，反应加入的核苷酸底物分子可以是A、G、C、T中的一种或多种，或者是A、G、C、U中的一种或多种。

23.根据权利要求10所述的方法，其特征在于，

所述测序中，检测的信号可以是电信号、生物荧光信号、化学荧光信号或者它们的组合。

24.根据权利要求10所述的方法，其特征在于，

所述的测序中，将核苷酸底物分子分为互不相同的两组，每次测序加入包含其中一组核苷酸底物分子的测序反应液；两组测序反应液循环加入，进行测序。

25.根据权利要求10所述的方法，其特征在于，

所述的测序中，将参考核酸序列和待测核酸序列同时进行测序；

参考核酸序列通过参数估计获得反应的超前、滞后、衰减系数、偏移量、单位信号信息；

通过参数估计获得的信息，校正待测核酸序列信号，获得校正的核酸序列。

26.根据权利要求10所述的方法，其特征在于，

所述的测序中，待测序列上连接有已知序列和数量的碱基，测序的时候，可以通过该已知序列的碱基的信号获得单位信号。

27.根据权利要求26所述的方法，其特征在于，

每个采样点的单位信号的是不相同的。

28.一种基因测序系统，包括计算机，其特征在于，

利用权利要求1-27任一项所述的方法获得校正的核酸序列。