CN101872386B - 利用双s形levenberg-marquardt和稳健线性回归的温度阶跃校正 - Google Patents

Abstract

本发明涉及利用双S形LEVENBERG-MARQUARDT和稳健线性回归的温度阶跃校正。本发明提供用于通过针对在PCR过程期间可能出现的温度变化对PCR数据进行校正来改进PCR扩增曲线中的Ct确定的系统和方法。具有通过Levenberg-Marquardt(LM)回归方法所确定的参数的双S形函数被用于找到在温度变化之后的区域中曲线的一部分的近似，其中出现温度变化的周期称为“CAC”。为在温度变化之前的区域中曲线的一部分确定稳健线性近似。利用线性近似和LM方法来确定周期CAC或CAC+1的荧光强度的值，并且从表示在出现温度变化之前曲线的一部分的数据集的一部分中减去这些值的差值，以产生已校正变化的数据集。该已校正变化的数据集可以被显示或另外被用于进一步的处理。

Description

利用双S形LEVENBERG-MARQUARDT和稳健线性回归的温度阶跃校正

本申请是申请日为2006年12月19日、申请号为200610169383.0以及发明名称为“利用双S形LEVENBERG-MARQUARDT和稳健线性回归的温度阶跃校正”的发明专利申请的分案申请。

技术领域

本发明广泛地涉及用于处理表示S形曲线或生长曲线的数据的系统和方法，并且更具体地涉及用于校正温度变化和用于确定PCR扩增曲线中的特征周期阈值(Ct)或肘值(elbowvalue)的系统和方法。

背景技术

聚合酶链反应(PCR)是一种用于酶合成或扩增限定的核酸序列的离体方法。该反应典型地利用两种寡核苷酸引物，这两种寡核苷酸引物杂交成相对股并且位于模板或要被扩增的目标DNA序列的侧面。通过热稳定DNA聚合酶来催化这些引物的延伸。包括由聚合酶引起的模板变性、引物退火和退火引物的延伸的一系列重复周期导致特定DNA片断的指数累积。荧光探针或标记典型地被用于促进扩增过程的检测和量化的方法中。

在图1中示出了典型的实时PCR曲线，其中针对典型的PCR过程画出了荧光强度值对周期数。在这种情况下，在PCR过程的每一周期中监控PCR产物的形成。通常在温度循环器中测量扩增，该温度循环器包括用于在扩增反应过程中测量荧光信号的部件和装置。这种温度循环器的例子是RocheDiagnosticsLightCycler(Cat.No.20110468)。借助荧光标记的杂交探针来例如检测扩增产物，该荧光标记的杂交探针仅仅在它们被结合到目标核酸上时才发射荧光信号，或者在某些情况下也借助结合到双股DNA上的荧光染料来例如检测扩增产物。

对于典型的PCR曲线来说，识别在基线区域末端处的通常被称为肘值或周期阈值(Ct)的过渡点极其有助于理解PCR扩增方法的特征。该Ct值可被用作PCR过程的效率的量度。例如，针对要被分析的所有反应确定所规定的信号阈值，并且针对目标核酸以及针对例如标准或看家(housekeeping)基因的参考核酸确定用于达到该阈值所需的周期数(Ct)。基于针对目标核酸和参考核酸所获得的Ct值，可以确定目标分子的绝对或相对拷贝数(Gibson等人的GenomeResearch6：995-1001；Bieche等人的CancerResearch59：2759-2765，1999年；WO97/46707；WO97/46712；WO97/46714)。图1中在基线区域15的末端处的区域20中的肘值将在周期数30的区域中。

在一些PCR试验、例如HIV试验中，在PC反应过程中典型地存在退火温度变化。该温度变化引起在出现温度变化的周期数时的荧光信号的随后变化。因此，为了计算正确的Ct值，有必要校正该信号变化。出现温度变化的周期是已知的，并且如果基线是完全平坦的并且不具有峰值，则校正该温度变化将是简单的事情。不幸地，基线常常是倾斜的并且也可能在任何位置处包含信号峰值(异常数据)。如果在温度变化位置处出现峰值，则更难以校正基线曲线。

因此，期望提供用于确定曲线、例如S形曲线或生长曲线、以及尤其是PCR曲线中的肘值的系统和方法，该系统和方法克服上述的和其它的缺点。具体地，该系统和方法应以对于例如异常数据的人为因素(artifact)来说可靠且稳健的方式实施温度阶跃校正。

发明内容

本发明提供用于通过针对在PCR过程期间可能出现的温度变化对PCR数据进行校正来改进PCR扩增曲线中的Ct确定的系统和方法。

根据一个方面，具有通过Levenberg-Marquardt(LM)回归方法所确定的参数的双S形函数被用于找到在温度变化之后的区域中曲线的一部分的近似，其中出现温度变化的周期被称为“CAC”。为在温度变化之前的区域中曲线的一部分确定稳健线性近似。利用线性近似和LM方法来确定周期CAC或CAC+1的荧光强度的值，并且从表示在出现温度变化之前曲线的一部分的数据集的一部分中减去这些值的差值，以产生已校正变化的数据集。然后，该已校正变化的数据集被返回，并且可以被显示或另外被用于进一步的处理。

在本发明的第一方面中，提供一种计算机实现的、校正具有基线部分和生长部分的聚合酶链反应(PCR)生长曲线的数据集中的温度阶跃变化的方法，该方法包括：

-接收聚合酶链反应生长曲线的数据集，其中所述数据集包括动态聚合酶链反应(PCR)过程的多个数据点，每一个数据点具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积；

-计算曲线的第一部分的线性近似，所述第一部分包括数据集中的数据点，这些数据点包括小于或等于在聚合酶链反应过程中出现退火温度变化的周期数(CAC)的周期数；

-利用曲线的第一部分的线性近似针对第一x值估计第一y值；

-通过将Levenberg-Marquardt(LM)回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，数据集的所述第二部分包括具有大于CAC的周期数的数据点；

-利用针对曲线的第二部分所计算的近似针对第一x值估计第二y值；

-确定第一和第二y值之间的差值；以及

-从对应于曲线的第一部分的数据点的每一个y值中减去该差值，以生成修改后的数据集。

在本发明的第二方面中，提供一种计算机可读介质，该计算机可读介质包括用于控制处理器以校正具有基线部分和生长部分的聚合酶链反应(PCR)生长曲线的数据集中的温度阶跃变化的代码，该代码包括指令，该指令用于：

-接收聚合酶链反应生长曲线的数据集，其中所述数据集包括动态聚合酶链反应(PCR)过程的多个数据点，每一个数据点具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积；

-计算曲线的第一部分的线性近似，所述第一部分包括数据集中的数据点，这些数据点包括小于或等于在聚合酶链反应过程中出现退火温度变化的周期数(CAC)的周期数；

-利用曲线的第一部分的线性近似针对第一x值估计第一y值；

-通过将Levenberg-Marquardt(LM)回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，数据集的所述第二部分包括具有大于CAC的周期数的数据点；

-利用针对曲线的第二部分所计算的近似针对第一x值估计第二y值；

-确定第一和第二y值之间的差值；以及

-从对应于曲线的第一部分的数据点的每一个y值中减去该差值，以生成修改后的数据集。

在本发明的又一方面中，提供一种动态聚合酶链反应(PCR)系统，包括：

-动态聚合酶链反应分析模块，该动态聚合酶链反应分析模块产生表示具有基线部分和生长部分的动态聚合酶链反应扩增曲线的聚合酶链反应数据集，其中所述数据集包括多个数据点，每一个数据点都具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积；以及

-智能模块，该智能模块适于通过以下方式来处理聚合酶链反应数据集以校正数据集中的温度阶跃变化：

-计算曲线的第一部分的线性近似，所述第一部分包括数据集中的数据点，这些数据点包括小于或等于在聚合酶链反应过程中出现退火温度变化的周期数(CAC)的周期数；

-利用曲线的第一部分的线性近似针对第一x值估计第一y值；

-通过将Levenberg-Marquardt(LM)回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，数据集的所述第二部分包括具有大于CAC的周期数的数据点；

-利用针对曲线的第二部分所计算的近似针对第一x值估计第二y值；

-确定第一和第二y值之间的差值；以及

-从对应于曲线的第一部分的数据点的每一个y值中减去该差值，以生成修改后的数据集。

附图说明

参照说明书的包括附图和权利要求的剩余部分，将认识到本发明的其它特征和优点。下面根据附图详细地描述本发明的进一步的特征和优点以及本发明的各种实施例的结构和操作。在附图中，相似的参考数字表示相同或功能相似的元件。

图1示出在PCR过程的上下文中的扩增曲线的例子。

图2A和2B(扩展的刻度)示出具有在周期15处的温度变化与在周期14处的峰值的PCR曲线的例子；图2C和2D(扩展的刻度)示出由用于处理整个PCR数据集的LM回归方法所产生的拟合；图2E和2F示出当从周期16至最后一个周期的数据点被用于LM回归时以满刻度和扩展的刻度得到的曲线拟合。

图3示出用于校正温度变化和确定单个S形曲线中的过渡值、例如动态PCR扩增曲线的肘值或Ct值的方法的一个实施例。

图4A针对具有温度变化的PCR数据集示出该数据集从周期1至周期CAC的部分的线性近似的例子并且示出的是稳健的线性估计；图4B针对图4A的PCR数据集示出该数据集从周期CAC+1至末端周期的部分的近似的例子并且示出的是稳健的双S形估计，其中根据本发明利用用于确定双S形函数的参数的LM回归方法来确定近似；图4C示出已校正温度变化的数据，即示出图4A和图4B中所示的数据集的已校正温度变化的数据的例子。

图5示出根据本发明的多个方面用于平分法的函数的典型形状。

图6示出根据本发明的一个实施例的峰值识别和替换方法的详细方法流程。

图7示出包括参数a-g的双S形方程的分解，其中参数a-g定义双S形曲线的形状和位置。

图8示出参数(d)对曲线和(e)的位置、拐点的x值的影响，除了参数d之外，图8中的所有曲线具有相同的参数值。

图9示出针对不同参数集的三个曲线形状的例子。

图10示出根据一个方面的用于确定双S形方程参数(e)和(g)的值的方法。

图11示出用于初始参数集的Levenberg-Marquardt回归方法的方法流程。

图12示出将稳健线性回归应用于从周期1至周期CAC的数据点的结果。

图13A和13B(扩展的刻度)示出在温度变化校正之后的完整的PCR数据集。

图14示出在峰值点去除之后的最终曲线。

图15示出描述软件和硬件资源之间的关系的一般框图。

具体实施方式

本发明提供用于针对在PCR过程期间可能出现的温度变化而校正PCR扩增曲线的系统和方法。温度变化的一个例子是在试验过程中在某一周期退火温度的受控变化。典型地，在通过基线区域所表示的过程的一部分期间出现该温度变化。该温度变化引起在出现变化的周期数时的荧光信号的随后变化。出现温度变化的周期于此将被称作CAC(代表退火变化的周期)。在某些方面中，具有通过Levenberg-Marquardt(LM)回归方法所确定的参数的双S形函数被用于找到在温度变化之后区域中的曲线的一部分的近似。为在温度变化之前区域中的曲线的该部分确定稳健线性近似。利用线性近似和LM方法来确定周期CAC+1的荧光强度的值，然后，从表示在出现温度变化之前曲线的该部分的数据集的一部分中减去这些值的差，以产生已校正变化的数据集。然后，该已校正变化的数据集被返回并且可以被显示或另外被用于进一步的处理。例如，已校正变化的数据集可以被用于确定PCR试验的Ct值。而且，已校正的数据集可以被处理以去除可能存在于数据集的为确定稳健线性近似而被处理的部分中的任何峰值。

在图1中示出了在PCR过程的上下文中的扩增曲线10的一个例子。如所示的那样，曲线10包括停滞期区域15和指数生长期区域25。停滞期区域15一般被称作基线或基线区域。这种曲线10包括连接停滞期和指数生长期区域的感兴趣的过渡区域20。区域20一般被称作肘或肘区域。肘区域20典型地限定基线的终点和基本过程的生长或扩增率的过渡。识别区域20中的特定过渡点可以有助于分析基本过程的特性。在典型的PCR曲线中，识别被称为肘值或周期阈值(Ct)的过渡点有助于理解PCR过程的效率特性。

可以提供类似的S形曲线或生长曲线的其它方法包括细菌法、酶法和结合法。在细菌生长曲线中，例如，感兴趣的过渡点已被称作停滞期中的时间θ。产生可根据本发明被分析的数据曲线的其它特定方法包括链置换扩增(SDA)方法、基于核酸序列的扩增(NASBA)方法和转录介导的扩增(TMA)方法。可以分别在Wang，Sha-Sha等人的“HomogeneousReal-TimeDetectionofSingle-NucleotidePolymorphismsbyStrandDisplacementAmplificationontheBDProbeTecETSystem”(ClinChem200349(10)：1599)和Weusten，JosJ.A.M.等人的“PrinciplesofQuantitationofViralLoadsUsingNucleicAcidSequence-BasedAmplificationinCombinationWithHomogeneousDetectionUsingMolecularBeacons”(NucleicAcidsResearch，200230(6)：26)中找到SDA和NASBA方法和数据曲线的例子。因此，尽管本文的剩余部分将在本发明对PCR曲线的适用性方面来讨论本发明的实施例和方面，但应理解的是，本发明可以被应用于涉及其它方法的数据曲线。

如图1中所示，例如可在二维坐标系统中表示典型的PCR生长曲线的数据，其中PCR周期数定义x轴并且累积的多核苷酸生长的指示器定义y轴。典型地，如图1中所示，累积生长的指示器是荧光强度值，因为荧光标记的使用可能是最广泛使用的标记方案。然而，应理解的是，可以根据所使用的具体标记和/或检测方案使用其它指示器。充裕或累积信号生长的其它有用的指示器的例子包括发光强度、化学发光强度、生物发光强度、磷光强度、电荷转移、电压、电流、功率、能量、温度、粘性、光散射、放射性强度、反射率、透射率和吸收率。周期的定义也可以包括时间、方法周期、单元操作周期和再生周期。

图2A和2B(扩展的刻度)示出具有在周期15处的温度变化与在周期14处的峰值的PCR曲线的例子。如果回归方法、例如LM回归方法被用于处理整个PCR数据集，将如图2C和2D(扩展的刻度)中所示的那样看到双S形对数据的全部拟合。图2C和2D表示双S形总体上拟合数据点，但是当然不能以令人满意的方式处理周期CAC处的不连续，其中在该周期CAC处出现退火温度变化。在该例子中，该CAC是周期15。如果替代地，从周期16至最后一个周期的数据点被用于LM回归，分别在图2E和2F中示出了在满刻度和扩展刻度曲线上所得到的曲线拟合。现在，对于从CAC+1至结束的数据集来说，图2F中的曲线拟合比图2D的曲线拟合更好，其中使用整个数据集。

一般方法概述

根据本发明，可以参照图3简要地描述用于校正温度变化和确定单个S形曲线中的过渡值、例如动态PCR扩增曲线的肘值或Ct值的方法100的一个实施例。在步骤110中，接收或另外获取表示曲线的实验数据集。也识别出现温度变化的周期。典型地，该周期值是先验地已知的，例如由提供数据的装置或设备记录。在图1中示出了所绘制的PCR数据集的例子，其中对于PCR曲线来说y轴和x轴分别表示荧光强度和周期数。在某些方面中，数据集应包括连续的并且沿着轴等间隔的数据。

在本发明的示范性实施例中，可以通过利用传统的个人计算机系统来实施该方法，该个人计算机系统包括、但不限于：用于输入数据集的例如键盘、鼠标等的输入装置；用于表示曲线区域中的感兴趣的特定点的显示装置、例如监控器；对于执行该方法中的每一步骤来说必要的处理装置、例如CPU；例如调制解调器的网络接口；用于存储数据集、在处理器上运行的计算机代码等的数据存储装置。此外，也可以在PCR装置中实施该方法。

在图15中显示了根据本发明的系统。该图示出说明软件和硬件资源之间的关系的一般框图。该系统包括可以位于温度循环器装置中的动态PCR分析模块和是计算机系统的一部分的智能模块。经由网络连接或直接连接将数据集(PCR数据集)从分析模块传输至智能模块，或反之亦然。由在处理器上运行的并存储于智能模块的存储装置上的计算机代码根据如图3中所示的方法来处理数据集，并在处理之后，将数据集传输回分析模块的存储装置，其中可以在显示装置上显示所修改的数据。在具体实施例中，也可以在PCR数据获取装置中实现智能模块。

在方法100在驻留于例如温度循环器的PCR数据获取装置中的智能模块(例如执行指令的处理器)中实施的情况下，在数据正被收集时数据集可以被实时地提供给智能模块，或者数据集可以被存储于存储单元或缓冲器中，并在已经完成实验之后被提供给智能模块。类似地，数据集可以经由通向获取装置的网络连接(例如LAN、VPN、内部网、因特网等)或直接连接(例如USB或其它直接有线或无线连接)被提供给单独的系统、例如桌上计算机系统，或在例如CD、DVD、软盘等的便携式介质上被提供。在某些方面中，数据集包括具有一对坐标值(或二维矢量)的数据点。对于PCR数据来说，该对坐标值典型地表示周期数和荧光强度值。在已经在步骤110中接收或获取数据集之后，可以进一步分析数据集，例如以便校正在试验过程中的温度变化以及确定基线区域的末端。

在步骤115中，根据从周期1至周期CAC的数据点计算稳健线性近似。图4A针对具有温度变化的PCR数据集示出该数据集从周期1至周期CAC的部分的线性近似的例子。该稳健线性近似被用于例如通过外推法估计CAC+1处的荧光水平。在步骤120中，计算曲线在CAC值之后的部分、例如曲线的扩增部分的近似。在该步骤过程中，在一个实施例中，具有通过Levenberg-Marquardt(LM)回归方法或其它回归方法所确定的参数的双S形函数被用于找到表示在数据集中从周期CAC+1至最后一个周期的数据集的部分的曲线的近似。该近似据说是“稳健”的，因为异常数据或峰值点对曲线拟合的质量具有最小的影响。图4B针对图4A的PCR数据集示出该数据集从周期CAC+1至最后周期的部分的近似的例子，其中利用Levenberg-Marquardt回归方法来确定近似，以确定双S形函数的参数。双S形估计步骤返回周期CAC+1处的所估计的荧光水平。它也针对存在于周期CAC+1和最后一个周期之间的曲线上的峰值对数据进行校正。

在步骤125中，通过计算利用稳健线性回归所计算的周期CAC+1处的荧光水平的差值来确定温度变化，并且确定稳健双S形估计。为了在步骤128中的温度变化校正，从起始于周期1至周期CAC的数据点中减去该值。图4C针对图4A和4B中所示出的数据集示出已校正温度变化的数据的例子。

在可选的步骤130中，利用已校正变化的数据来执行峰值去除方法(例如Levenberg-Marquardt峰值去除方法)以便检测并校正任何剩余峰值(具体地，经阶跃校正的曲线的基线区域中的峰值)。如将在下面更详细讨论的，图6示出用于识别并替换表示PCR或其它生长曲线的数据集中的峰值点的方法流程。

在步骤135中，返回修改后的、已校正变化的数据集，例如用于如图4C中所示的显示或进一步的处理。例如，在步骤140中，然后处理标准化曲线，以确定Ct值，并将结果(已校正变化的数据和/或Ct值)例如返回至执行分析的系统或请求分析的单独的系统。可以利用与执行图3的分析的系统相耦合的显示装置、例如监控器屏幕或打印机来再现图形显示，或者可将数据提供给单独的系统以便在显示装置上再现。可以利用已校正变化的数据集根据各种方法来确定Ct值。例如，在一个方面中，US专利申请序列号11/316315和11/349550的教导可以被用于确定Ct值。

稳健线性回归

温度变化出现在PCR曲线的基线区域中。典型地，可以利用形式为f(x)＝slope·x+intercept的线性函数来近似PCR曲线的基线区域，其中x表示周期数。为了找到线性近似的变量slope(斜率)和intercept(截距)的最佳值，稳健线性回归被用于一个方面中。可以使用其它回归方法或拟合函数，然而稳健线性拟合优选于其它方法、例如最小平方拟合，以便最小化例如峰值的人为因素对回归质量的影响。

根据一个方面，稳健线性回归使量

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|y_i-\mathrm{intercept}-\mathrm{slope}\·x_i|---\left(1\right)$

最小化，并且对于斜率的固定值来说，使方程(1)最小化的截距的值是

intercept＝median{y_i-slope·x₁}(2)

于是，slope的值可以作为方程的根被找到：

$0=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i\·\mathrm{sgn}(y_i-\mathrm{intercept}-\mathrm{slope}\·x_i)---\left(3\right)$

其中sgn()是如下被定义的符号函数

x＞0，sgn(x)＝1

x＝0，sgn(x)＝0

x＜0，sgn(x)＝-1。

在将方程(2)代入方程(3)之后，存在一个要求解的未知方程。在某些方面中，平分法被用于对所得到的方程求解。

平分法是一种求根算法，其通过重复地将区间分成两半并且选择存在根的子区间来工作。该平分法被用于用以求解下面的方程的本发明的稳健线性回归：

$0=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i\·\mathrm{sgn}(y_i-\mathrm{median}\{y_i-\mathrm{slope}\·x_i\}-\mathrm{slope}\·x_i)---\left(4\right)$

其中x_i是周期数，并且y_i是周期x_i处的荧光水平，并且N＝CAC。该平分法被用于找到变量slope的值，对于该值来说方程(4)被最小化：该变量表示PCR曲线中的基线的斜率，并且在某些方面中取[-0.2；0.2]区间中的值。在图5中示出了方程(4)中的函数的典型形状。

平分法线性地收敛，这可以是缓慢的，但如果所考虑的区间包含根，则被确保收敛。在一个方面中，初始收敛区间被定义为[-0.2；0.2]，并执行平分算法的100次迭代。应理解的是，可以使用更小或更大的区间，并且可以执行更多或更少次迭代。

下面提供根据本发明的方面的、用于执行平分法的Mathematica^TM代码的例子。

abValue[data_]：＝Module[{}，

RSum[b_]：＝Sum[(data[[i，1]]*Sign[(data[[i，2]]-a[b]-bdata[[i，1]])])，

{i，1，Length[data]}]；

a[b_]：＝Median[Table[(data[[i，2]]-bdata[[i，1]])，{i，1，Length[data]}]]；

Lx＝-0.2；

Rx＝0.2；

RsumL＝RSum[Lx]；

RsumR＝RSum[Rx]；

Nx＝(Lx+Rx)/2；

RsumNx＝RSum[Nx]；

For[i＝l，i≤100，i++，

{

If[RsumL＜0&&RsumNx＞0，Rx＝Nx]；

If[RsumL＜0&&RsumNx＜0，Lx＝Nx]；

If[RsumL＞0&&RsumNx＜0，Rx＝Nx]；

If[RsumL＞0&&RsumNx＞0，Lx＝Nx]；

RsumL＝RSum[Lx]；

RsumR＝RSum[Rx]；

Nx＝(Lx+Rx)/2；

PT[[i，1]]＝i；

PT[[i，2]]＝Lx；

PT[[i，3]]＝Rx；

RsumNx＝RSum[Nx]；

}]

]

LM回归方法

如上面提到的，图6示出用于识别和替换表示PCR或其它生长曲线的数据集中的峰值点的方法流程。图6的步骤502至524也示出用于近似数据集的曲线并确定拟合函数的参数的方法流程。在一个实施例中，Levenberg-Marquardt(LM)方法被用于根据从周期CAC+1至最后一个周期的曲线点来计算数据集的稳健曲线近似。LM回归被用于计算周期CAC+1处的荧光水平的估计，以便使例如峰值的人为因素对估计的质量可能具有的影响最小化。该LM方法是非线性回归方法；它是一种使非线性函数和数据集之间的距离最小化的迭代技术。该方法的作用就像最陡下降方法和Gauss-Newton方法的组合一样：在当前近似没有很好地拟合时，它的作用就像最陡下降方法一样(更慢，但更可靠地收敛)，但在当前近似变得更精确时，它的作用于是将像Gauss-Newtom方法一样(更快，但不太可靠地收敛)。

一般地，LM回归方法包括需要各种输入并提供输出的算法。在一个方面中，输入包括要处理的数据集、被用于拟合数据的函数、以及函数的参数或变量的初始猜测。该输出包括用于使函数和数据集之间的距离最小化的函数的参数集。

根据一个实施例，拟合函数是以下形式的双S形：

$f\left(x\right)=a+\mathrm{bx}+\frac{c}{(1+{\exp }^{-d(x-e)})(1+{\exp }^{-f(x-g)})}.---\left(5\right)$

作为拟合函数的该方程的选择基于其拟合典型的PCR曲线或其它生长曲线可能采用的不同曲线形状的灵活性和能力。本领域的技术人员应理解，可以随意地使用上面的拟合函数或其它拟合函数的变型。

双S形方程(5)具有7个参数：a、b、c、d、e、f和g。该方程可被分解为常数、斜率和双S形的和。双S形自身是两个S形的相乘。图7示出双S形方程(5)的分解。参数d、e、f和g决定两个S形的形状。为了示出它们对最终曲线的影响，考虑单个S形：

$\frac{1}{1+{\exp }^{-d(x-e)}},---\left(6\right)$

其中参数d决定曲线的“锐度”，并且参数e决定拐点的x值。图8示出参数d对曲线的影响以及参数e对拐点的x值的位置的影响。在下面，表1描述参数对双S形曲线的影响。

表1：双S形参数描述

参数	对曲线的影响
		a	在x＝0时y的值
b	基线和曲线的平稳段的斜率
		c	曲线的AFI
d	第一S形的“锐度”(见图9)
		e	第一S形的拐点的位置(见图9)
f	第二S形的“锐度”
		g	第二S形的拐点的位置

在一个方面中，为了防止曲线采用不切实际的形状，应限制双S形方程的“锐度”参数d和f。因此，在一个方面中，d＜-1或d＞1.1或f＜-1或f＞1.1的任何迭代都被认为是不成功的。在其它方面中，可以使用对参数d和f的不同限制。

因为Levenberg-Marquardt算法是迭代算法，所以典型地需要用于拟合的函数的参数的初始猜测。初始猜测越好，近似将越好，并且算法将向局部最小值收敛就越不太可能。由于双S形函数的复杂性和PCR曲线或其它生长曲线的各种形状，所以对每一参数的一次初始猜测可能不足以防止算法有时向局部最小值收敛。因此，在一个方面中，输入多个(例如3个或更多个)初始参数集并保存最佳结果。在一个方面中，大部分参数被保持为在所使用的多个参数集上恒定；仅仅参数c、d和f可以对多个参数集中的每一个来说是不同的。图9示出针对不同参数集的三个曲线形状的例子。这三个参数集的选择指示表示PCR数据的曲线的三种可能的不同形状。应理解的是，可以处理多于三个的参数集并保存最佳结果。

如图6中所示，在步骤510中确定LM方法的初始输入参数。这些参数可以通过操作者输入或可以被计算。根据一个方面，根据如下面所论述的步骤502、504和506确定或设置参数。

初始参数(a)的计算：

参数(a)是基线的高度；它的值对于所有初始参数集来说是相同的。在一个方面中，在步骤504中，参数(a)被赋予数据集中的第三个最低y轴值、例如荧光值。这提供稳健的计算。当然，在其它方面中，参数(a)可被随意地赋予任何其它荧光值、例如最低y轴值、第二最低值等等。

初始参数(b)的计算：

参数(b)是基线和曲线的平稳段的斜率。它的值对于所有初始参数集来说是相同的。在一个方面中，在步骤502中，将0.01的静态值赋予(b)，因为理想地不应存在任何斜率。在其它方面中，参数(b)可被赋予不同值，例如范围从0至大约0.5的值。在一个方面中，值(b)表示从CAC+1至基线末端的基线斜率。

初始参数(c)的计算：

参数(c)表示曲线的平稳段的高度减去基线的高度，这被表示为绝对荧光增加或AFI。在一个方面中，对于第一个参数集来说，c＝AFI+2，而对于最后两个参数来说，c＝AFI。这在图9中被示出，其中对于最后两个参数集来说，c＝AFI。对于第一个参数集来说，c＝AFI+2。该改变是由于通过第一个参数集所建模的曲线的形状，该形状不具有曲线的平稳段。

参数(d)和(f)的计算：

参数(d)和(f)限定两个S形的锐度。因为不存在针对这些参数基于曲线来给出近似的方法，所以在一个方面中在步骤502中使用三个静态代表值。应理解的是，其它静态或非静态值可被用于参数(d)和/或(f)。这些对对所遇到的PCR曲线的最普遍的形状建模。在下面，表2示出针对如图9中所示的不同参数集的(d)和(f)的值。

表2：参数d和f的值

参数集编号	d的值	f的值
			1	0.1	0.7
2	1.0	0.4
			3	0.35	0.25

参数(e)和(g)的计算：

在步骤506中，确定参数(e)和(g)。参数(e)和(g)限定两个S形的拐点。在一个方面中，它们在所有初始参数集中都采用相同的值。参数(e)和(g)可以具有相同或不同的值。为了找到近似，在一个方面中，使用超过强度、例如荧光的平均值的第一点(其不是峰值)的x值。根据该方面的用于确定(e)和(g)的值的方法在图10中被示出并在下面进行论述。

参照图10，最初，确定曲线(例如荧光强度)的平均值。接着确定超过平均值的第一数据点。然后确定是否：

a.该点不位于曲线的起点附近、例如在最初的5个周期内；

b.该点不位于曲线的终点附近、例如在最后的5个周期内；以及

c.该点周围的导数(例如在它周围的2个点的半径中)没有显示出符号的任何改变。如果它们显示出符号的任何改变，则该点可能是峰值并且因此应被丢弃。

在下面，表3示出根据一个方面的如图9中所使用的初始参数值的例子。

表3：初始参数值：

返回图6，一旦在步骤510中设置了所有参数，就利用输入数据集、函数和参数执行LM方法520。传统上，Levenberg-Marquardt方法被用于解决非线性最小平方问题。传统的LM方法计算被定义为曲线近似和数据集之间的误差的平方和的距离量度。然而，当最小化该平方和时，它给予异常数据一个重要的权重，因为它们的距离大于非峰值数据点的距离，这常常导致不合宜的曲线或不太理想的曲线。因此，根据本发明的一个方面，通过最小化绝对误差的和来计算近似和数据集之间的距离，因为这不给予异常数据同样大的权重。在该方面中，通过下式来给出近似和数据之间的距离：

距离＝∑|y_数据-y_近似|。(7)

如上所述，在一个方面中，如在步骤522和524中所示，输入并处理多个(例如3个)初始参数集中的每一个并保存最佳结果，其中最佳参数是在方程(7)中提供最小或最小化距离的参数集。在一个方面中，大部分参数在多个参数集中被保持恒定；仅仅c、d和f可以对于每一个参数集来说是不同的。应理解的是，可以使用任何数量的初始参数集。

图11示出根据本发明的用于参数集的LM方法520的方法流程。如上面所解释的，Levenberg-Marquardt方法的作用可以就像最陡下降方法或像Gauss-Newton方法一样。它的作用取决于衰减因数λ。λ越大，Levenberg-Marquardt算法的作用就将越像最陡下降方法一样。另一方面，λ越小，Levenberg-Marquardt算法的作用就将越像Gauss-Newton方法一样。在一个方面中，λ以0.001开始。应理解的是，λ可以以任何其它值开始，例如从大约0.000001至大约1.0。

如前所述，Levenberg-Marquardt方法是迭代技术。根据一个方面，如图11中所示，在每一次迭代过程中执行下面的操作：

1.计算在先的近似的Hessian矩阵(H)。

2.计算在先的近似的转置Jacobian矩阵(J^T)。

3.计算在先的近似的距离向量(d)。

4.使Hessian矩阵对角线增大当前衰减因数λ：

H_aug＝Hλ(8)

5.求解增大的方程：

H_augx＝J^Td(9)

6.将增大的方程的解x添加到函数的参数中。

7.计算新的近似和曲线之间的距离。

8.如果具有新的参数集的距离小于具有先前的参数集的距离：

·认为该迭代是成功的。

·保存或存储该新的参数集。

·将衰减因数λ例如减小因数10。

如果具有新的参数集的距离大于具有先前的参数集的距离：

·认为该迭代是不成功的。

·丢弃该新的参数集。

·将衰减因数λ例如增大因数10。

在一个方面中，图11的LM方法迭代，直至达到了下列标准之一：

1.已经运行了规定次数的(N次)迭代。该第一标准防止算法无限地迭代下去。例如，在如图10中所示的一个方面中，缺省迭代值N是100。如果算法能够收敛，则100次迭代应足以使算法收敛。一般地，N的范围可以从小于10至100或更大。

2.两次成功迭代之间的距离的差值小于阈值、例如0.0001。当该差值变得非常小时，已经达到所期望的精度并且继续迭代是无意义的，因为解不会变得显著更好。

3.衰减因数λ超过规定值，例如大于10²⁰。当λ变得非常大时，该算法不会比当前解更好地收敛，因此继续迭代是无意义的。一般地，该规定值可以大大小于或大于10²⁰。

例子

对生长曲线的从周期1至周期CAC的数据点应用稳健线性回归导致图12中所示出的图。可以看出，曲线拟合的质量是极好的，因为温度变化和峰值都不影响最终结果。然后，稳健线性回归和双S形被用于估计周期CAC+1处的荧光信号的值，在这种情况下CAC+1是周期16。然后，这些值之间的差值被用于阶跃校正。应注意的是，估计周期CAC处的稳健线性近似并把双S形外推至CAC值以确定差值是同样有效的。在该数据集中，周期16处的稳健线性方程的估计值是：

左＝2.91669

周期16处的稳健双S形的估计值是：

右＝2.11582

并且差值是：

dif＝左-右＝0.800866

然后，从数据集的对应于周期1至CAC的部分的荧光信号中减去该差值。在图13A和13B(扩展的刻度)中示出了在调节之后的完整的数据集与双S形Levenberg-Marquafdt曲线拟合。在一个方面中，LM方法然后被用于去除周期15处的异常数据。在图14中示出了具有在周期15处被插入的替换点的最终曲线。现在校正了由于退火温度变化而导致的峰值和荧光强度变化。

结论

根据本发明的一个方面，提供一种计算机实现的用于校正具有基线部分和生长部分的聚合酶链反应(PCR)生长曲线的数据集中的温度阶跃变化的方法。该方法典型地包括接收PCR生长曲线的数据集，其中该数据集包括动态聚合酶链反应(PCR)过程的多个数据点，每一个数据点具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积。该方法也典型地包括计算曲线的第一部分的线性近似，该第一部分包括数据集中的数据点，该数据点包括小于或等于在PCR过程中出现退火温度变化的周期数(CAC)的周期数；并通过将Levenberg-Marquardt(LM)回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，数据集的第二部分包括具有大于CAC的周期数的数据点。该方法也典型地包括利用曲线的第一部分的线性近似来针对第一x值估计第一y值，利用针对曲线的第二部分所计算的近似来针对第一x值估计第二y值，确定第一和第二y值之间的差值，并从对应于曲线的第一部分的数据点的每一个y值中减去该差值以生成修改后的数据集。在某些方面中，第一x值是CAC。在其它方面中，第一x值是CAC+1。在另一实施例中，CAC数据点被包括在数据的第一部分中。

在某一实施例中，计算线性近似的步骤包括将稳健线性回归应用于曲线的第一部分。在具体实施例中，应用稳健线性回归包括应用求根方法。在又一具体实施例中，求根方法是平分法。

在某些实施例中，该方法进一步包括以下步骤：通过将LM回归方法应用于修改后的数据集来确定对应于数据曲线的第一部分的一个或多个数据点是否是异常峰值；以及从修改后的数据集中去除或替换所识别出的峰值的数据值。在又一确定的实施例中，该方法进一步包括利用修改后的数据集确定周期阈值(Ct)。

在另一实施例中，通过荧光强度值、发光强度值、化学发光强度值、磷光强度值、电荷转移值、生物发光强度值或吸收率值之一来表示扩增多核苷酸的累积。

在某些方面中，双S形函数具有下面的形式：

并且计算包括迭代地确定该函数的参数a、b、c、d、e、f和g中的一个或多个。

根据本发明的另一方面，提供一种计算机可读介质，该计算机可读介质包括用于控制处理器以校正具有基线部分和生长部分的聚合酶链反应(PCR)生长曲线的数据集中的温度阶跃变化的代码。该代码典型地包括指令，该指令用于接收PCR生长曲线的数据集，其中该数据集包括动态聚合酶链反应(PCR)过程的多个数据点，每一个数据点都具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积。该代码也典型地包括用于计算曲线的第一部分的线性近似并通过将Levenberg-Marquardt(LM)回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似的指令，其中该第一部分包括数据集中的数据点，这些数据点包括小于或等于在PCR过程中出现退火温度变化的周期数(CAC)的周期数，数据集的第二部分包括具有大于CAC的周期数的数据点。该代码也典型地包括用于利用曲线的第一部分的线性近似针对第一x值估计第一y值、利用针对曲线的第二部分所计算的线性针对第一x值估计第二y值、确定第一和第二y值之间的差值、并从对应于曲线的第一部分的数据点的每一个y值中减去该差值以生成修改后的数据集的指令。在某些方面中，第一x值是CAC。在其它方面中，第一x值是CAC+1。在另一方面中，CAC数据点被包括在数据的第一部分中。

在某些实施例中，该代码进一步包括用于通过将LM回归方法应用于修改后的数据集来确定对应于数据曲线的第一部分的一个或多个数据点是否是异常峰值、以及用于从修改后的数据集中去除或替换所识别出的峰值的数据值的指令。

在另一实施例中，用于计算线性近似的指令包括用于对曲线的第一部分应用稳健线性回归的指令。在具体实施例中，用于应用稳健线性回归的指令包括用于应用求根方法的指令。在又一具体实施例中，求根方法是平分法。

在某一实施例中，该代码进一步包括用于利用修改后的数据集来确定周期阈值(Ct)的指令。

在另一实施例中，通过荧光强度值、发光强度值、化学发光强度值、磷光强度值、电荷转移值、生物发光强度值或吸收率值之一来表示扩增多核苷酸的累积。

在又一实施例中，双S形函数具有下面的形式：

并且用于计算的指令包括用于迭代地确定该函数的参数a、b、c、d、e、f和g中的一个或多个的指令。

根据本发明的又一方面，提供一种动态聚合酶链反应(PCR)系统，该系统典型地包括动态聚合酶链反应分析模块，该动态聚合酶链反应分析模块产生表示具有基线部分和生长部分的动态PCR扩增曲线的PCR数据集，其中该数据集包括多个数据点，每一个数据点都具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积，以及该系统典型地包括智能模块，该智能模块适于通过以下方式来处理PCR数据集以校正数据集中的温度阶跃变化，即计算曲线的第一部分的线性近似并通过将Levenberg-Marquardt(LM)回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，其中该第一部分包括数据集中的数据点，这些数据点包括小于或等于在PCR过程中出现退火温度变化的周期数(CAC)的周期数，数据集的第二部分包括具有大于CAC的周期数的数据点。该智能模块也典型地适于利用曲线的第一部分的线性近似针对第一x值估计第一y值、利用针对曲线的第二部分所计算的线性针对第一x值估计第二y值、确定第一和第二y值之间的差值、并从对应于曲线的第一部分的数据点的每一个y值中减去该差值以生成修改后的数据集。在某些方面中，第一x值是CAC。在其它方面中，第一x值是CAC+1。在另一方面中，CAC数据点被包括在数据的第一部分中。

在动态PCR系统的某些实施例中，动态PCR分析模块驻留在动态温度循环器装置中，并且智能模块包括可通信耦合至分析模块的处理器。在其它某些实施例中，智能模块包括驻留在通过网络连接或直接连接之一耦合至分析模块的计算机系统中的处理器。

在另一实施例中，该智能模块适于通过将LM回归方法应用于修改后的数据集来确定对应于数据曲线的第一部分的一个或多个数据点是否是异常峰值以及从修改后的数据集中去除或替换所识别出的峰值的数据值。

在某一实施例中，计算线性近似包括对曲线的第一部分应用稳健线性回归。在具体实施例中，应用稳健线性回归包括应用求根方法。在又一具体实施例中，求根方法是平分法。

在另一实施例中，该智能模块适于利用修改后的数据集来确定周期阈值(Ct)。

在又一实施例中，通过荧光强度值、发光强度值、化学发光强度值、磷光强度值、电荷转移值、生物发光强度值或吸收率值之一来表示扩增多核苷酸的累积。

在某些实施例中，双S形函数具有下面的形式：

并且计算包括迭代地确定该函数的参数a、b、c、d、e、f和g中的一个或多个。

应理解的是，包括曲线近似和回归方法的Ct确定方法可以以在计算机系统的处理器上运行的计算机代码来实施。该代码包括用于控制处理器以实施Ct确定方法的各个方面和步骤的指令。该代码典型地被存储于硬盘、RAM或例如CD、DVD等的便携式介质上。类似地，可以在例如温度循环器的PCR装置中实施该方法，该温度循环器包括执行存储于耦合至处理器的存储单元中的指令的处理器。可以通过通向代码源的网络连接或直接连接或利用众所周知的便携式介质将包括这种指令的代码下载至PCR装置存储单元。

本领域技术人员应理解的是，可以利用例如C、C++、C#、Fortran、VisualBasic等的多种编程语言以及例如Mathematica的、提供有助于数据可视化和分析的预先打包的例行程序、函数和程序的应用来对本发明的各种方法进行编码。应用的另一例子是

虽然已经通过例子并根据具体实施例描述了本发明，但应理解的是本发明并不局限于所公开的实施例。相反，旨在覆盖对本领域技术人员来说将是显而易见的各种修改和相似的配置。例如在进行多个温度阶跃校正的情况下，可以对数据集的不同部分执行多个回归。例如，在实现了PCR数据的基线部分的两个温度变化的情况下，线性回归可被应用于包括周期1至第一变化的周期的数据集部分，并且第二线性回归可被应用于第一变化的周期至第二变化的周期的数据集部分，并且LM方法可被应用于第二变化的周期至最后周期的数据集。因此，所附的权利要求的范围应符合最宽泛的解释，以便包括所有的这种修改和相似的配置。

Claims (17)

1.一种用于测量聚合酶链反应PCR过程的效率的方法，该方法包括以下步骤：

－针对要被分析的所有反应确定所规定的信号阈值并且针对目标核酸以及针对参考核酸确定达到该阈值所需的周期数，其中通过以下方式来校正具有基线部分和生长部分的聚合酶链反应PCR曲线的数据集中的温度阶跃变化以及确定周期阈值Ct：

－接收聚合酶链反应生长曲线的数据集，其中所述数据集包括动态聚合酶链反应PCR过程的多个数据点，每一个数据点具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积；

－计算曲线的第一部分的线性近似，所述第一部分包括数据集中的数据点，这些数据点包括小于或等于在聚合酶链反应过程中出现退火温度变化的周期数CAC的周期数；

－利用曲线的第一部分的线性近似针对第一x值估计第一y值；

－通过将Levenberg－MarquardtLM回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，数据集的所述第二部分包括具有大于CAC的周期数的数据点；

－利用针对曲线的第二部分所计算的近似针对第一x值估计第二y值；

－确定第一和第二y值之间的差值；以及

－从对应于曲线的第一部分的数据点的每一个y值中减去该差值，以生成修改后的数据集。

2.如权利要求1所述的方法，其中第一x值是x＝CAC+1。

3.如权利要求1所述的方法，其中第一x值是x＝CAC。

4.如权利要求1所述的方法，进一步包括：

－通过将LM回归方法应用于修改后的数据集来确定对应于数据曲线的第一部分的一个或多个数据点是否是异常峰值；以及

－从修改后的数据集中去除或替换所识别出的峰值的数据值。

5.如权利要求1所述的方法，其中CAC数据点被包括在数据的第一部分中。

6.如权利要求1所述的方法，其中计算线性近似的步骤包括将稳健线性回归应用于曲线的第一部分。

7.如权利要求6所述的方法，其中应用稳健线性回归包括应用求根方法。

8.如权利要求7所述的方法，其中求根方法是平分法。

9.如权利要求1所述的方法，进一步包括利用修改后的数据集来确定周期阈值Ct。

10.如权利要求1所述的方法，其中双S形函数具有下面的形式：

$a+bx+\frac{c}{(1+{\exp }^{-d(x-e)})(1+{\exp }^{-f(x-g)})},$ 并且其中计算包括迭代地确定该函数的参数a、b、c、d、e、f和g中的一个或多个，其中参数a表示在x＝0时y的值，参数b表示基线和曲线的平稳段的斜率，参数c表示曲线的AFI，参数d表示第一S形的锐度，参数e表示第一S形的拐点的位置，参数f表示第二S形的锐度，以及参数g表示第二S形的拐点的位置。

11.一种动态聚合酶链反应PCR系统，包括：

－动态聚合酶链反应分析模块，该动态聚合酶链反应分析模块产生表示具有基线部分和生长部分的动态聚合酶链反应扩增曲线的聚合酶链反应数据集，其中所述数据集包括多个数据点，每一个数据点具有一对坐标值(x，y)，其中x表示周期数，并且y表示扩增多核苷酸的累积；以及

－智能模块，该智能模块适于通过以下方式来处理聚合酶链反应数据集以校正数据集中的温度阶跃变化以及确定周期阈值Ct：

－计算曲线的第一部分的线性近似，所述第一部分包括数据集中的数据点，这些数据点包括小于或等于在聚合酶链反应过程中出现退火温度变化的周期数CAC的周期数；

－利用曲线的第一部分的线性近似针对第一x值估计第一y值；

－通过将Levenberg－MarquardtLM回归方法应用于数据集的第二部分和双S形函数以确定该函数的参数来计算曲线的第二部分的近似，数据集的所述第二部分包括具有大于CAC的周期数的数据点；

－利用针对曲线的第二部分所计算的近似针对第一x值估计第二y值；

－确定第一和第二y值之间的差值；以及

－从对应于曲线的第一部分的数据点的每一个y值中减去该差值以生成修改后的数据集。

12.如权利要求11所述的系统，其中第一x值是x＝CAC+1。

13.如权利要求11所述的系统，其中第一x值是x＝CAC。

14.如权利要求11所述的系统，其中所述动态聚合酶链反应分析模块驻留于动态温度循环器系统中，并且其中所述智能模块包括可通信耦合至分析模块的处理器。

15.如权利要求11所述的系统，其中所述智能模块包括驻留在通过网络连接或直接连接之一耦合至分析模块的计算机系统中的处理器。

16.如权利要求11所述的系统，其中所述智能模块适于

－通过将LM回归方法应用于修改后的数据集来确定对应于数据曲线的第一部分的一个或多个数据点是否是异常峰值；以及

－从修改后的数据集中去除或替换所识别出的峰值的数据值。

17.如权利要求11所述的系统，其中双S形函数具有下面的形式：

$a+bx+\frac{c}{(1+{\exp }^{-d(x-e)})(1+{\exp }^{-f(x-g)})},$ 并且其中计算包括迭代地确定该函数的参数a、b、c、d、e、f和g中的一个或多个，其中参数a表示在x＝0时y的值，参数b表示基线和曲线的平稳段的斜率，参数c表示曲线的AFI，参数d表示第一S形的锐度，参数e表示第一S形的拐点的位置，参数f表示第二S形的锐度，以及参数g表示第二S形的拐点的位置。