测定PCR的直线-S形生长曲线的循环阈值的方法

摘要

本发明公开一种测定PCR的直线‑S形生长曲线的循环阈值的方法，包括以下步骤：S1：获取代表PCR生长曲线的数据点的集合；S2：基于由线性函数和S形函数组成的分段函数对数据集进行拟合，搜索确定最优分段点，以将PCR生长曲线分为直线段和S形曲线段；S3：分别计算得到的直线段的拟合函数，和S形曲线段的拟合函数。本发明采用由线性函数和S形函数组成的分段函数对其进行全局搜索最优分段点拟合，可提高拟合函数的拟合程度，进而再经二次求导法可求得更为准确的循环阈值，进而有利于临床中依据样本的循环阈值大小来正确做出阴阳性的判断。

权利要求

1.一种测定PCR的直线‑S形生长曲线的循环阈值的方法，其特征在于，包括以下步骤：S1：获取代表PCR生长曲线的数据点的集合 ，其中， 代表循环数， 代表荧光强度；

S2：基于由线性函数和S形函数组成的分段函数对数据集进行拟合，搜索确定最优分段点，以将PCR生长曲线分为直线段和S形曲线段；

S3：分别计算得到的直线段的拟合函数，和S形曲线段的拟合函数；

S4：基于二次求导法，计算该分段函数的二阶导数；

S5：比较S4中计算的二阶导数的大小，其最大值点对应的循环数即为循环阈值；

所述S2中采用的S形函数为Logistic函数，且搜索确定最优分段点具体包括以下步骤：S21：分别写出线性函数 ，和Logistic函数 的表达式，其中，其中， 是线性函数 的第一个系数；

是线性函数 的第二个系数；

是S形函数 的第一个系数；

是S形函数 的第二个系数；

是S形函数 的第三个系数；

是S形函数 的第四个系数；

S22：设分段点为 ，并以该分段点为界限将数据集P分成 和 两组，其中，S23：基于 对 进行拟合，其拟合直线的误差平方和为S24：对 再扩充如下所示的 个元素：得到扩充后的数据集

并基于 对 进行拟合，其拟合曲线在分段点 之后部分的误差平方和为S25：构建模型拟合程度的决定系数 ，其中， 是均值；

S26：求出 的最大值点 ，此时，最优分段点即为 。

2.根据权利要求1所述的测定PCR的直线‑S形生长曲线的循环阈值的方法，其特征在于，所述S3具体包括以下步骤：S31：联立线性方程组计算 的参数 和 ；

S32：基于莱文贝格－马夸特算法计算 的参数 和 。

3.根据权利要求2所述的测定PCR的直线‑S形生长曲线的循环阈值的方法，其特征在于，所述S31具体包括以下步骤：S311：将数据集 代入 ，可得到：S312：令

则 ，即可求得 ，其中，是由数据点 组成的矩阵；

是由数据点 组成的矩阵；

是由参数 和参数 组成的矩阵；

是矩阵 的转置。

说明书

测定PCR的直线‑S形生长曲线的循环阈值的方法

技术领域

[0001] 本发明涉及PCR技术领域，尤其涉及一种测定PCR的直线‑S形生长曲线的循环阈值的方法。

背景技术

[0002] 循环阈值（Ct值），一般落在PCR生长曲线的指数阶段，在阈值交点法中表示为PCR反应使得探测到的荧光信号强度达到设定阈值时所对应的循环数，而在二次求导法中表示
为PCR反应的扩增速度的变化率最大值点。待测样品中模板的起始拷贝数越大，Ct值越小，
所以Ct值可以用来关联模板起始拷贝数。在临床中依据样本的Ct值大小来判断阴阳性。

[0003] PCR生长曲线不是标准的指数曲线，这是因为随着PCR进行的循环数增加，DNA聚合酶活性变低、反应物枯竭、反应副产物焦磷酸阻碍合成反应等等因素，使得PCR不是指数扩
增，典型的PCR生长曲线呈“S形”，包括基线期、指数期、线性期、平台期。基线期通常是在基
准区域，指数期开始于基线期的末端，且比基线期有明显更高的荧光数值。在典型PCR曲线
中，确定循环阈值(Ct值)表征PCR反应的数量特征。

[0004] 而在真实应用场景中，存在反应液蒸发、试剂探针降解、冻干球复溶、临床样本的物质干扰等，而导致扩增曲线异常，因而曲线不典型，即不呈“S形”，常表现为基线期不是水
平波动而是呈倾斜向上或向下直线，需要分别识别区分基线期和指数期。

[0005] 而二次求导法是以扩增速度的变化率最大点计算Ct值，因此，其准确性依赖于拟合函数的拟合程度，此时，对于非典型的S形生长曲线，如果直接使用S形函数拟合，因拟合
程度较低，进而用二次求导法进行求解的Ct值准确性较低，同时，阈值交点法的计算是基于
基线处的信号，容易受到检测系统影响，如基线漂移，且选择阈值和确定Ct值往往由操作人
员人工干预，导致计算的阈值带有主观因素判断，因此，需要对其进行改进。

发明内容

[0006] 本发明的目的是提供一种测定PCR的直线‑S形生长曲线的循环阈值的方法，该方法可在当PCR生长曲线是呈“直线‑S形”，即“直线向上‑S形”或者“直线向下‑S形”时，采用由
线性函数和S形函数组成的分段函数对其进行全局搜索最优分段点拟合，可提高拟合函数
的拟合程度，进而再经二次求导法可求得更为准确的循环阈值，进而有利于临床中依据样
本的循环阈值大小来正确做出阴阳性的判断。

[0007] 为实现上述目的，采用以下技术方案：

[0008] 一种测定PCR的直线‑S形生长曲线的循环阈值的方法，包括以下步骤：

[0009] S1：获取代表PCR生长曲线的数据点的集合 ，其中， 代表循环数， 代表荧光强度；

[0010] S2：基于由线性函数和S形函数组成的分段函数对数据集进行拟合，搜索确定最优分段点，以将PCR生长曲线分为直线段和S形曲线段；

[0011] S3：分别计算得到的直线段的拟合函数，和S形曲线段的拟合函数；

[0012] S4：基于二次求导法，计算该分段函数的二阶导数；

[0013] S5：比较S4中计算的二阶导数的大小，其最大值点对应的循环数即为循环阈值。

[0014] 进一步地，所述S2中采用以下任意一种S形函数进行拟合：Logistic函数、Sigmoid函数、Gompertz函数、Chapman函数。

[0015] 进一步地，所述S2中采用的S形函数为Logistic函数，且搜索确定最优分段点具体包括以下步骤：

[0016] S21：分别写出线性函数 ，和Logistic函数 的表达式，其中，

[0017]

[0018] 其中， 是线性函数 的第一个系数；

[0019] 是线性函数 的第二个系数；

[0020] 是S形函数 的第一个系数；

[0021] 是S形函数 的第二个系数；

[0022] 是S形函数 的第三个系数；

[0023] 是S形函数 的第四个系数；

[0024] S22：设分段点为 ，并以该分段点为界限将数据集P分成 和 两组，其中，

[0025]

[0026] S2 3：基于 对 进行 拟合 ，其拟合直 线的误差 平方 和为

[0027] S24：对 再扩充如下所示的 个元素：

[0028]

[0029] 得到扩充后的数据集

[0030]

[0031] 并基于 对 进行拟合，其拟合曲线在分段点 之后部分的误差平方和

[0032] S25：构建模型拟合程度的决定系数 ，

[0033] 其中， 是 均值；

[0034] S26：求出 的最大值点 ，此时，最优分段点即为 。

[0035] 进一步地，所述S3具体包括以下步骤：

[0036] S31：联立线性方程组计算 的参数 和 ；

[0037] S32：基于莱文贝格－马夸特算法计算 的参数 和 。

[0038] 进一步地，所述S31具体包括以下步骤：

[0039] S311：将数据集 带入 ，可得到：

[0040]

[0041] S312：令

[0042]

[0043] 则 ，即可求得 ，其中，

[0044] 是由数据点 组成的矩阵；

[0045] 是由数据点 组成的矩阵；

[0046] 是由参数 和参数 组成的矩阵；

[0047] 是矩阵 的转置。

[0048] 采用上述方案，本发明的有益效果是：

[0049] 该方法可在当PCR生长曲线是呈“直线‑S形”，即“直线向上‑S形”或者“直线向下‑S形”时，采用由线性函数和S形函数组成的分段函数对其进行全局搜索最优分段点拟合，可
提高拟合函数的拟合程度，进而再经二次求导法求得更为准确的循环阈值，进而有利于临
床中依据样本的循环阈值大小来正确做出阴阳性的判断。

附图说明

[0050] 图1为本发明的流程性框图；

[0051] 图2为本发明其中一实施例中，对“直线向下‑S形”生长曲线进行拟合的结果示意图；

[0052] 图3为对图2的拟合函数求二阶导数的结果示意图；

[0053] 图4为本发明其中一实施例中，对“直线向上‑S形”生长曲线进行拟合的结果示意图；

[0054] 图5为对图4的拟合函数求二阶导数的结果示意图。

具体实施方式

[0055] 以下结合附图和具体实施例，对本发明进行详细说明。

[0056] 参照图1至5所示，本发明提供一种测定PCR的直线‑S形生长曲线的循环阈值的方法，包括以下步骤：

[0057] S1：获取代表PCR生长曲线的数据点的集合 ，其中， 代表循环数， 代表荧光强度；

[0058] S2：基于由线性函数和S形函数组成的分段函数对数据集进行拟合，搜索确定最优分段点，以将PCR生长曲线分为直线段和S形曲线段；

[0059] S3：分别计算得到的直线段的拟合函数，和S形曲线段的拟合函数；

[0060] S4：基于二次求导法，计算该分段函数的二阶导数；

[0061] S5：比较S4中计算的二阶导数的大小，其最大值点对应的循环数即为循环阈值。

[0062] 进一步地，所述S2中采用以下任意一种S形函数进行拟合：Logistic函数、Sigmoid函数、Gompertz函数、Chapman函数。

[0063] 其中，所述S2中采用的S形函数为Logistic函数，且搜索确定最优分段点具体包括以下步骤：

[0064] S21：分别写出线性函数 ，和Logistic函数 的表达式，其中，

[0065]

[0066] 其中， 是线性函数 的第一个系数；

[0067] 是线性函数 的第二个系数；

[0068] 是S形函数 的第一个系数；

[0069] 是S形函数 的第二个系数；

[0070] 是S形函数 的第三个系数；

[0071] 是S形函数 的第四个系数；

[0072] S22：设分段点为 ，并以该分段点为界限将数据集P分成 和 两组，其中，

[0073]

[0074] S2 3 ：基于 对 进行 拟合 ，其拟 合直 线的误 差平 方和 为

[0075] S24：对 再扩充如下所示的 个元素：

[0076]

[0077] 得到扩充后的数据集

[0078]

[0079] 并基于 对 进行拟合，其拟合曲线在分段点 之后部分的误差平方和为

[0080]

[0081] S25：构建模型拟合程度的决定系数 ，其中，

[0082] 是 均值；

[0083] S26：求出 的最大值点 ，此时，最优分段点即为 。

[0084] 所述S3具体包括以下步骤：

[0085] S31：联立线性方程组计算 的参数 和 ；

[0086] S32：基于莱文贝格－马夸特算法计算 的参数 和 。

[0087] 所述S31具体包括以下步骤：

[0088] S311：将数据集 带入 ，可得到：

[0089]

[0090] S312：令

[0091]

[0092] 则 ，即可求得 ，其中，

[0093] 是由数据点 组成的矩阵；

[0094] 是由数据点 组成的矩阵；

[0095] 是由参数 和参数 组成的矩阵；

[0096] 是矩阵 的转置。

[0097] 本发明工作原理：

[0098] 继续参照图1至5所示，首先获取代表PCR生长曲线的数据点的集合，其中， 代表循环数， 代表荧光强度；随后，基于由线性函数和
S型函数组成的分段函数对数据集进行拟合，搜索确定最优分段点，以将PCR生长曲线分为
直线段和S型曲线段；再分别计算得到的直线段的拟合函数，和S型曲线段的拟合函数；随
后，基于二次求导法，计算其二阶导数；最后，比较计算的二阶导数的大小，其最大值点对应
的循环数即为循环阈值，具体地：

[0099] 由于PCR的标准生长曲线形如S形，而在数学上有一类S形函数，故可以用S形函数拟合标准生长曲线，S形函数包括如下表格中所示的四种：

[0100]

[0101] 表1 四种S型函数的参数及表达式

[0102] 该实施例中，采用Logistic函数。首先要确定最优分段点，对于直线‑S形的生长曲线，需要分直线和S形函数全局搜索最优分段点拟合；对于代表生长曲线的数据集
，设分段拟合函数 由线性函数 和S形函
数 组成。其中，

[0103] 是线性函数 的第一个系数；

[0104] 是线性函数 的第二个系数；

[0105] 是S形函数 的第一个系数；

[0106] 是S形函数 的第二个系数；

[0107] 是S形函数 的第三个系数；

[0108] 是S形函数 的第四个系数。

[0109] 设分段点为 ，将数据集P分成两组，其中，

[0110]

[0111] 对 扩充 个元素： ，得到，

[0112] ，共有 个元素；

[0113] 基 于 对 进 行 拟 合 ，其 拟 合 直 线 的 误 差 平 方 和 为，

[0114] 基于 对 进行拟合，其拟合曲线在分段点 之后部分的误差平方和为。

[0115] 决定系数 可以度量模型拟合的程度，其数值介于0到1。 的值越大，说明模型拟合得越好。 的计算公式为：

[0116] ，其中， 表示 的平均值。

[0117] 求出 的最大值点 ，那么全局最优的分段点为 。

[0118] 此外，上述方法为分一段直线和一段S形函数拟合，可根据具体的实际需求拓展为多段直线和一段S形函数拟合。

[0119] 求得最优分段点后，就需要求 和 的相关参数，其中，对于 ，联立方程组求解即可：

[0120] 可首先将数据集 带入 ，可得到：

[0121]

[0122] 令

[0123]

[0124] 则 ，即可求得 ，其中，

[0125] 是由数据点 组成的矩阵；

[0126] 是由数据点 组成的矩阵；

[0127] 是由参数 和参数 组成的矩阵；

[0128] 是矩阵 的转置。

[0129] 而对于 ，则基于莱文贝格－马夸特算法求其参数，莱文贝格－马夸特算法对于非线性最小化问题能给出数值解。该算法在求解问题迭代过程中可以修改参数，达到结
合高斯‑牛顿算法和梯度下降法的优点，改善两者不足。

[0130] 假设 是一个非线性映射， ，对于 ，有 ,

[0131] 目的是任意给定 和适合的初始值 ，找到一个 ,使误差平方和 尽可能小（局部极小），其中 。

[0132] 求解是一个迭代过程，给出近似的数值解 。根据泰勒公式展开，得到函数的一阶微分近似。该近似的好处有两个：一是线性，二是只需一阶微分。

[0133] ，

[0134] 其中 是 的雅可比矩阵。每次迭代的目标是：假若 是迭代点，要找到 ，使得 最小。

[0135] 由投影公式知有最小误差时满足下面的公式：

[0136]

[0137] 把公式修改为：

[0138]

[0139] 这种计算方法便是莱文贝格‑马夸特算法。当 较大时，该算法与最速下降法接近；当 较小时，该算法与高斯‑牛顿法接近。初始可以先用一个较小的 ，若 长度增大
则增加 ，从而确保每次迭代可使 的长度减少。

[0140] 迭代终止的条件可以设置为：

[0141] 1） 长度变化小于给定值；

[0142] 2） 长度变化小于给定值；

[0143] 3）迭代次数达到设定上限。

[0144] 其中，其中一实施例中，求得：

[0145]

[0146] 对其求二阶导数，其有唯一的最大值点26.7，即循环阈值为26.7。

[0147] 在分别求得 和 的相关参数后，随后，基于二次求导法，分别计算这两个函数的二阶导数，二阶导数最大值点即为循环阈值，实际计算应用结果可见图2至5所示。

[0148] 以上仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。