CN114058683B - 一种核酸扩增结果判定方法 - Google Patents

Abstract

本发明涉及核酸检测技术领域中的一种核酸扩增结果判定方法，包括以下步骤：获取不同循环数的荧光强度数据，并进行曲线拟合，得到拟合函数；根据粒子群优算法计算拟合函数的最优参数，并代入拟合函数，对拟合函数进行一阶导数求解；通过求导计算荧光信号上升或下降的连续次数，并连续次数确定扩增区间；计算二阶导数，得出二阶导数的最大值的解和对应的时间T；判断时间T是否在扩增区间内，若是，则判定时间T即为C_t值，并根据样本定义阈值，判断扩增结果；若否，则根据时间T与扩增区间的位置关系判断扩增结果，具有灵敏度高的优点，突破了实际数据在进行移动平移法时具有滞后性，不能很好反应变化趋势的瓶颈。

Description

一种核酸扩增结果判定方法

技术领域

本发明涉及核酸检测技术领域，具体涉及一种核酸扩增结果判定方法。

背景技术

核酸扩增和检测技术是当今生物学研究的重要手段之一，包括生物技术、基础科学、医学研究等在内的各领域的科学家们用这种方法进行广泛的应用研究。核酸扩增技术包括常规PCR、实时荧光PCR、数字PCR以及目前临床伴随诊断(POCT)技术流行的恒温扩增技术，如LAMP、CPA技术等。这类方法会产生大量的原始数据值，数据处理方法的不同可能会显著的影响最终的结果。

核酸扩增技术常用于用于基因表达的定量分析，或者是常见传染病的定性检测，核酸检测目前是新型冠状病毒检测的“金标准”，具有早期诊断、灵敏度和特异性高等特点。在实时荧光PCR中，模板的定量有两种方法：相对定量和绝对定量。

绝对定量指用已知的标准曲线来推算未知的样本量。分析的结果是给定数量的样品中(给定数量的细胞，每微克总RNA)中核酸的量(拷贝数、微克)。相对定量的分析结果是在相当量的试验组(样品A)和对照组(样品B)中一个靶基因的相对比率(倍数差异)。绝对定量中，未知样本的量(即拷贝数或单位数量)可以通过从已知量的标准品范围中推算得出。为了建立标准曲线，需要已知浓度的模板。模板稀释后，用这些稀释样品作为标准样品，将未知的待测样本和标准样置于同一次实验中进行反应，用稀释的标准样品构建标准曲线，通过推算来确定末知样本中目的基因的量。这与用分子量标记来判定琼脂糖凝胶上未知DNA条带分子量大小的原理相似。

相对定量指在一定样品中靶序列相对于另一参照序列的量的变化。相对定量分为两种,一是以质量单位作为标准进行相对定量，用质量单位(如纽胞数目或核酸的微克数)而不用参照基因作为标准的优点是，实验设计概念简单，数据数学分析处理简单易学。二是以参照基因作为标准进行相对定量，用参照基因而不用质量单位作为标准。这个方法能准确量化出事材料的使用量，当样本用量受限时比较方便。数据分析主要使用双标曲线法和2^-△△Ct法。

无论是定量实验还是定性实验，数据处理都需要进行其中两部核心步骤，数据噪声过滤和Ct值确认。噪声过滤目前主流的方法有三个步奏，数据平滑、背景扣除和幅度归一化。平滑采用移动平均法，简单的三点移动能够体现曲线的趋势并且抹去一些非特异性的峰值。背景扣除是PCR数据处理中的常见步骤，一般是减去运行中观察最小值，再数据平滑处理后进行基线减法能够有效降低背景噪音。幅度归一化统一了不同样本的中心位置，能够降低批次效应带来的影响。其次，Ct值确认主流是选择通过阈值线与扩增曲线交点确认法，阈值线可以人为设定或者默认设置为3-15个循环的荧光信号的标准偏差的10倍。

移动平均法能够体现出扩增曲线的趋势，计算速度快，相对简单，但是对预测实际数据不敏感，具有一定的滞后性，而且并不能总是很好的反应出变化趋势，特别是初始浓度较低的扩增，预测值可能活总停留再过去水平上而无法预计会导致将来更低或者更高的波动。阈值线与扩增曲线交点确认法对阈值线设定区间的噪音敏感。

发明内容

本发明针对现有技术中的缺点，提供了一种核酸扩增结果判定方法，具有灵敏度高的优点，突破了实际数据在进行移动平移法时具有滞后性，不能很好反应变化趋势的瓶颈。

为了解决上述技术问题，本发明通过下述技术方案得以解决：

一种核酸扩增结果判定方法，包括以下步骤：

获取不同循环数的荧光强度数据，并对荧光强度数据进行曲线拟合，得到拟合函数；

根据粒子群优算法计算拟合函数的最优参数，将最优参数代入拟合函数，并对拟合函数进行一阶导数求解；

根据所述一阶导数，通过求导计算荧光信号上升或下降的连续次数，并根据荧光信号上升或下降的连续次数确定扩增区间；

对所述一阶导数计算二阶导数，并得出二阶导数的最大值的解和最大解对应的时间T；

判断时间T是否在扩增区间内，若是，则判定时间T即为C_t值，并根据样本定义阈值，判断扩增结果；若否，则判断时间T与扩增区间的位置关系；

若时间T位于扩增区间的起始位置之前，则直接定义C_t值为一分钟位置，且扩增结果为阳性；若时间T位于扩增区间的结束位置之后，则直接定义C_t值为四十分钟位置，且扩增结果为阴性。

可选的，根据粒子群优算法计算拟合函数的最优参数，包括以下步骤：

初始化荧光信号所形成的粒子群，将粒子群分为若干个局部子集，并对局部子集进行迭代计算，得到结果子集，并将所述结果子集存储于最优参数集内。

可选的，对局部子集进行迭代计算，包括以下步骤：

设定学习率α和参数θ的初始值，并通过模型计算粒子当前位置所对应的目标位置，同时计算损失函数；

根据适应度函数计算粒子个体的适应值；

将所述损失函数沿参数θ的负梯度方向更新，并得到损失函数最小值；

计算粒子个体的局部最优位置向量和粒子群的全局最优向量，并根据损失函数分别更新粒子个体的速度向量和位置向量；

根据所述粒子个体的速度向量和位置向量，更新粒子个体的局部位置向量和全局位置向量，并存入结果子集中；

判断更新后的局部子集的损失函数的大小，得到更优解，并将最优解存入最优参数集；

判断粒子群是否满足参数优化后的终止条件，若满足，则停止迭代计算；若不满足，则重新设定学习率并重新计算粒子个体的适应值。

可选的，所述损失函数的公式为：

其中，J(θ)为损失函数；m代表样本个数；i代表m个样本中的第i个；h_θ(t_i)代表第i个样本t时间点拟合函数对荧光强度的预测值；y_i代表第i个样本的真实荧光强度数据。

可选的，将所述损失函数沿参数θ的负梯度方向更新的更新公式为：

可选的，判断粒子群是否满足参数优化后的终止条件，包括以下步骤：

设定参数θ一次迭代的变化值为△θ，若△θ＜0.0001，则结果收敛，判断结果收敛速度，并根据结果收敛速度判断是否满足终止条件；若△θ≥0.0001，则结果不收敛，不满足终止条件，迭代继续，并减小学习率α的初始值，计算损失函数。

可选的，根据结果收敛速度判断是否满足终止条件，包括以下步骤：

若结果收敛速度缓慢，则不满足终止条件，迭代继续，并增大学习率α的初始值，计算损失函数；

若结果收敛速度正常或较快，则满足终止条件，停止迭代并输出拟合函数的最优参数。

可选的，根据荧光信号上升或下降的连续次数确定扩增区间，包括以下步骤：

设定状态阈值G，并设定荧光信号上升或下降的连续次数大于状态阈值G时，荧光信号上升或下降的连续区间为第一区间；

获取所述第一区间的扩增幅度，并设定扩增幅度阈值D，判断是否第一区间的拟合函数的数值大于0，且第一区间的扩增幅度大于扩增幅度阈值D；

若是，则第一区间为扩增区间，且扩增曲线为潜在阳性曲线；若否，则扩增曲线为阴性曲线。

可选的，拟合函数的计算公式为：

其中，θ＝(θ₁，θ₂，θ₃，θ₄，θ₅，θ₆，θ₇，θ₈)，θ₁、θ₂、θ₃和θ₄部分函数为对正常的扩增S型曲线进行拟合，对荧光扩增曲线数据的非线性特征参数进行调节；θ₅对荧光扩增曲线数据的不同平移进行拟合；θ₆、θ₇和θ₈对荧光扩增曲线数据的线性特征的斜率进行调节；y为观察值。

可选的，所述非线性特征包括曲线沿X轴的拉伸幅度和沿Y轴的拉伸幅度。

采用本发明提供的技术方案，与现有技术相比，具有如下有益效果：

通过划定扩增区间进行二阶导数最大值计算判定Ct值方法，该方法对真实数据值敏感，能够更好的反映出扩增的变化情况，灵敏度高，其次扩增区间二阶导数法确定Ct值鲁棒性较强，无需进行基线设定，且常规的的PCR检测和恒温扩增技术都可以通用。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提出的一种核酸扩增结果判定方法的流程图；

图2为本发明实施例提出的一种核酸扩增结果判定方法的粒子群优算法流程图；

图3为本发明实施例提出的一种核酸扩增结果判定方法的实验得到的荧光强度数据对应的一阶导数和二阶导数分布曲线图以及通过拟合函数得到的荧光强度数据分布。

具体实施方式

下面结合实施例对本发明做进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例。

如图1所示，一种核酸扩增结果判定方法，包括以下步骤：获取不同循环数的荧光强度数据，其中不同循环数的荧光强度数据需通过实验得到，并对实验得到的荧光强度数据进行曲线拟合，得到拟合函数，其中，拟合函数的计算公式为：

其中θ＝(θ₁，θ₂，θ₃，θ₄，θ₅，θ₆，θ₇，θ₈)，θ₁、θ₂、θ₃和θ₄部分函数为对正常的扩增S型曲线进行拟合，对荧光扩增曲线数据的非线性特征参数进行调节；θ₅对荧光扩增曲线数据的不同平移进行拟合；θ₆、θ₇和θ₈对荧光扩增曲线数据的线性特征的斜率进行调节；y为观察值，非线性特征包括曲线沿X轴的拉伸幅度和沿Y轴的拉伸幅度。

如图2所示，根据粒子群优算法计算拟合函数的最优参数，将最优参数代入拟合函数，并对拟合函数进行一阶导数求解，得到如图3所示的一阶导数曲线图，具体的，根据粒子群优算法计算拟合函数的最优参数，包括以下步骤：初始化荧光信号所形成的粒子群，将粒子群分为若干个局部子集，并对局部子集进行迭代计算，得到结果子集，并将结果子集存储于最优参数集内，

最优参数集用于存储每次迭代计算得到的结果子集，且每次迭代计算仅对局部子集进行更新，而局部子集是指在算法过程中挑选出整个粒子群中的一部分进行寻优和计算，剔除劣解，采用随机初始化每个粒子的位置向量和速度，从而根据位置向量和速度进行迭代寻优。

对局部子集进行迭代计算，包括以下步骤：设定学习率α和参数θ的初始值，并通过模型计算粒子当前位置所对应的目标位置，同时计算损失函数；其中，损失函数的公式为：

其中，J(θ)为损失函数；m代表样本个数；i代表m个样本中的第i个；h_θ(t_i)代表第i个样本t时间点拟合函数对荧光强度的预测值；y_i代表第i个样本的真实荧光强度数据。

进一步的，根据适应度函数计算粒子个体的适应值；将损失函数沿参数θ的负梯度方向更新，并得到损失函数最小值，使用梯度下降法来估计参数θ的取值，由于梯度下降算法每次仅随机使用一组数据进行训练，所以m＝1，将损失函数沿参数θ的负梯度方向更新的更新公式如下：

计算粒子个体的局部最优位置向量和粒子群的全局最优向量，并根据损失函数分别更新粒子个体的速度向量和位置向量；根据粒子个体的速度向量和位置向量，更新粒子个体的局部位置向量和全局位置向量，并存入结果子集中；判断更新后的局部子集的损失函数的大小，得到更优解，并将最优解存入最优参数集；

具体的，判断跟新后的局部自己的损失函数的大小包括以下步骤，计算更新后粒子个体的适应值，其中适应值即为损失函数值，并比较更新后粒子个体的适应值与更新前的粒子个体的适应值的大小，若更新后的粒子个体的适应值大于更新前的粒子个体的适应值，则更新后的粒子个体的适应值为更优解。

判断粒子群是否满足参数优化后的终止条件，若满足，则停止迭代计算；若不满足，则重新设定学习率并重新计算粒子个体的适应值，其中，判断粒子群是否满足参数优化后的终止条件，包括以下步骤：设定参数θ一次迭代的变化值为△θ，若△θ＜0.0001，则结果收敛，判断结果收敛速度，并根据结果收敛速度判断是否满足终止条件；若△θ≥0.0001，则结果不收敛，不满足终止条件，迭代继续，并减小学习率α的初始值，计算损失函数，其中0.0001为求解最优解与设定目标之间的精度差值。

具体的，判断结果收敛速度的方法为，观察单位时间内，△θ的数值大小，若△θ较小离目标值比较远则代表收敛速度缓慢，反之则收敛正常或较快。

根据结果收敛速度判断是否满足终止条件，包括以下步骤：若结果收敛速度缓慢，则不满足终止条件，迭代继续，并增大学习率α的初始值，计算损失函数；若结果收敛速度正常或较快，则满足终止条件，停止迭代并输出拟合函数的最优参数。

一般在算法开始的时候保持最快的惯性速度，随着迭代次数增加，距离目标越近，减小惯性速度，避免速度过快飞过目标，以达到最好的收敛效果。

根据一阶导数，通过求导计算荧光信号上升或下降的连续次数，并根据荧光信号上升或下降的连续次数确定扩增区间，其中求导计算荧光信号上升或下降的连续次数的求导公式为：

从而通过求导公式得出连续次数，另一方面，根据荧光信号上升或下降的连续次数确定扩增区间，包括以下步骤：设定状态阈值G，并设定荧光信号上升或下降的连续次数大于状态阈值G时，荧光信号上升或下降的连续区间为第一区间；获取第一区间的扩增幅度，并设定扩增幅度阈值D，判断是否第一区间的拟合函数的数值大于0，且第一区间的扩增幅度大于扩增幅度阈值D；若是，则第一区间为扩增区间，且扩增曲线为潜在阳性曲线；若否，则扩增曲线为阴性曲线，其中扩增幅度为根据第一区间的荧光最大值和最小值之间的差值得到。

如图3所示，对一阶导数计算二阶导数，得到二阶导数曲线图，并得出二阶导数的最大值的解和最大解对应的时间T；判断时间T是否在扩增区间内，若是，则判定时间T即为C_t值，并根据样本定义阈值，判断扩增结果；若否，则判断时间T与扩增区间的位置关系；若时间T位于扩增区间的起始位置之前，则直接定义C_t值为一分钟位置，且扩增结果为阳性；若时间T位于扩增区间的结束位置之后，则直接定义C_t值为四十分钟位置，且扩增结果为阴性。

其中，根据样本定义阈值可以人工设定或者按照常规样本的扩增阈值设定，而不同样本的阳性扩增阈值存在差异，此处的扩增阈值指扩增区间的扩增幅度；而上述中的一分钟位置与四十分钟位置为真实的下机数据，一般为四十分钟的数据，按照二阶导数最大值判定C_t值位置可能在0之前或者40之后，即为越界。

从而通过划定扩增区间进行二阶导数最大值计算判定Ct值方法，该方法对真实数据值敏感，能够更好的反映出扩增的变化情况，灵敏度高，其次扩增区间二阶导数法确定Ct值鲁棒性较强，无需进行基线设定，且常规的的PCR检测和恒温扩增技术都可以通用。

此外，需要说明的是，本说明书中所描述的具体实施例，其零、部件的形状、所取名称等可以不同。凡依本发明专利构思所述的构造、特征及原理所做的等效或简单变化，均包括于本发明专利的保护范围内。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离本发明的结构或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。

Claims (6)

1.一种核酸扩增结果判定方法，其特征在于，包括以下步骤：

获取不同循环数的荧光强度数据，并对荧光强度数据进行曲线拟合，得到拟合函数；

根据粒子群优化算法计算拟合函数的最优参数，将最优参数代入拟合函数，并对拟合函数进行一阶导数求解，其中，根据粒子群优化算法计算拟合函数的最优参数，包括以下步骤：初始化荧光信号所形成的粒子群，将粒子群分为若干个局部子集，并对局部子集进行迭代计算，得到结果子集，并将所述结果子集存储于最优参数集内，

具体地，对局部子集进行迭代计算，包括以下步骤：设定学习率α和参数θ的初始值，并通过模型计算粒子当前位置所对应的目标位置，同时计算损失函数；根据适应度函数计算粒子个体的适应值；将所述损失函数沿参数θ的负梯度方向更新，并得到损失函数最小值；计算粒子个体的局部最优位置向量和粒子群的全局最优位置向量，并根据损失函数分别更新粒子个体的速度向量和位置向量；根据所述粒子个体的速度向量和位置向量，更新粒子个体的局部位置向量和全局位置向量，并存入结果子集中；判断更新后的局部子集的损失函数的大小，得到更优解，并将最优解存入最优参数集；判断粒子群是否满足参数优化后的终止条件，若满足，则停止迭代计算；若不满足，则重新设定学习率并重新计算粒子个体的适应值，

其中，所述拟合函数的计算公式为：

且θ＝(θ₁，θ₂，θ₃，θ₄，θ₅，θ₆，θ₇，θ₈)，θ₁、θ₂、θ₃和θ₄部分函数为对正常的扩增S型曲线进行拟合，对荧光扩增曲线数据的非线性特征参数进行调节；θ₅是对荧光扩增曲线数据的不同平移进行拟合；θ₆、θ₇和θ₈是对荧光扩增曲线数据的线性特征的斜率进行调节；h_θ(t)代表在t时间点的拟合函数对荧光强度的预测值；

根据所述一阶导数，通过求导计算荧光信号上升或下降的连续次数；

设定状态阈值G，并设定荧光信号上升或下降的连续次数大于状态阈值G时，荧光信号上升或下降的连续区间为第一区间；

获取所述第一区间的扩增幅度，并设定扩增幅度阈值D，判断是否第一区间的拟合函数的数值大于0，且第一区间的扩增幅度大于扩增幅度阈值D；

若是，则第一区间为扩增区间，且扩增曲线为潜在阳性曲线；若否，则扩增曲线为阴性曲线；

根据所述一阶导数计算二阶导数，并得出二阶导数的最大值的解和最大值的解对应的时间T；

判断时间T是否在扩增区间内，若是，则判定时间T即为C_t值，并根据样本定义阈值，判断扩增结果；若否，则判断时间T与扩增区间的位置关系；

若时间T位于扩增区间的起始位置之前，则直接定义C_t值为一分钟位置，且扩增结果为阳性；若时间T位于扩增区间的结束位置之后，则直接定义C_t值为四十分钟位置，且扩增结果为阴性。

2.根据权利要求1所述的一种核酸扩增结果判定方法，其特征在于，所述损失函数的公式为：

其中，J(θ)为损失函数；m代表样本个数；i代表m个样本中的第i个样本；h_θ(t_i)代表第i个样本在t时间点的拟合函数对荧光强度的预测值；y_i代表第i个样本的真实荧光强度数据。

3.根据权利要求2所述的一种核酸扩增结果判定方法，其特征在于，将所述损失函数沿参数θ的负梯度方向更新的更新公式为：

4.根据权利要求1所述的一种核酸扩增结果判定方法，其特征在于，判断粒子群是否满足参数优化后的终止条件，包括以下步骤：

设定参数θ一次迭代的变化值为Δθ，若Δθ＜0.0001，则结果收敛，判断结果收敛速度，并根据结果收敛速度判断是否满足终止条件；若Δθ≥0.0001，则结果不收敛，不满足终止条件，迭代继续，并减小学习率α的初始值，计算损失函数。

5.根据权利要求4所述的一种核酸扩增结果判定方法，其特征在于，根据结果收敛速度判断是否满足终止条件，包括以下步骤：

若结果收敛速度缓慢，则不满足终止条件，迭代继续，并增大学习率α的初始值，计算损失函数；

若结果收敛速度正常或较快，则满足终止条件，停止迭代并输出拟合函数的最优参数。

6.根据权利要求1所述的一种核酸扩增结果判定方法，其特征在于，所述非线性特征包括曲线沿X轴的拉伸幅度和沿Y轴的拉伸幅度。

Images