CN112685686B - 一种熔解曲线平滑方法 - Google Patents

一种熔解曲线平滑方法 Download PDF

Info

Publication number
CN112685686B
CN112685686B CN202110288235.5A CN202110288235A CN112685686B CN 112685686 B CN112685686 B CN 112685686B CN 202110288235 A CN202110288235 A CN 202110288235A CN 112685686 B CN112685686 B CN 112685686B
Authority
CN
China
Prior art keywords
sliding window
polynomial
fluorescence intensity
derivative
fitting
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.)
Active
Application number
CN202110288235.5A
Other languages
English (en)
Other versions
CN112685686A (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.)
Hangzhou Bori Technology Co Ltd
Original Assignee
Hangzhou Bori 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 Hangzhou Bori Technology Co Ltd filed Critical Hangzhou Bori Technology Co Ltd
Publication of CN112685686A publication Critical patent/CN112685686A/zh
Application granted granted Critical
Publication of CN112685686B publication Critical patent/CN112685686B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Investigating, Analyzing Materials By Fluorescence Or Luminescence (AREA)

Abstract

本发明提供了一种熔解曲线平滑方法，包括获取荧光定量熔解实验中非等间隔荧光强度数据；基于类Savitzky‑Golay拟合方法，确定最优滑动窗口大小以及多项式次数；基于最优滑动窗口大小和多项式次数，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵；基于导数矩阵，计算各滑动窗口内各采样点对应的负导数数值，获得平滑后的熔解曲线。本发明通过采用类Savitzky‑Golay求导方法，获得了非等间隔采样平滑的熔解曲线，简化了非等间隔数据的求导过程，易于理解，容易实现。

Description

一种熔解曲线平滑方法
技术领域
本发明涉及荧光定量检测技术领域，尤其是涉及一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法。
背景技术
目前，诊断传染病标准的核酸定量方法仍主要采用荧光实时定量PCR，其能够量化样品模板的初始值，常被用在基因分析表达、转基因食物检测和癌症检测中。在 PCR 扩增反应完成后，常通过逐渐增加温度、使扩增产物发生降解来获得荧光强度变化值即熔解曲线，从而对扩增产物的特异性进行考察。可见，要获得熔解曲线，需要获得荧光强度的变化情况，即求取荧光强度的导数。
对于分辨率较高、采样点数多且等间隔的荧光强度值，直接差分法与实际相差不大，唯一需要解决的问题是如何消除差分带来的噪声，通常的方法有两种：其一是在差分之前对原始强度数据进行平滑，这种方法最早由 Norris 等人提出，常称为Norris求导法，当然也可采用其他的平滑如Savitzky-Golay平滑方法等；其二是在直接差分之后再进行平滑滤波。
但对于分辨率较低、采样点数较少，特别是不等间隔采样数据（注：为了节约时间，在做熔解实验时，得到的随温度变化的荧光强度数据常常是不等距的），直接差分法误差较大。此时，通常的做法是对非等间隔采样数据进行曲线拟合，然后再进行均匀上采样插值，从而一方面将不等间隔采样数据变成了等间隔，另一方面提高采样的分辨率，再在此基础上继续采用平滑+直接差分法即可。但这样处理会带来以下问题：增加了求导过程的复杂程度。
为此，我们提出一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法。Savitzky-Golay拟合是一种用多项式实现滑动窗内最小二乘拟合的平滑方法，其可通过最小二乘计算得到与平滑系数类似的导数系数。但Savitzky-Golay方法通常针对的是规则等间隔数据，对于非等间隔数据的求导，会出现求取的导数不能正确缩放的问题。为了解决这一问题，我们仿照Savitzky-Golay平滑求导的思路，首先利用最小二乘法对各滑动窗口内荧光强度数据进行多项式拟合，获得多项式系数矩阵，再求取相应的导数矩阵（注：我们称之为类Savitzky-Golay求导方法），然后基于导数矩阵，计算滑动窗口内各采样点对应的负导数数据，获得熔解曲线值。因为在多项式拟合中，滑动窗口大小和拟合阶次对算法的处理效果有很大影响，我们在该过程引入决定系数R2作为评价指标，以便获得最优的拟合窗口大小及拟合阶次。本发明通过采用类Savitzky-Golay求导方法，获得了非等间隔采样平滑的熔解曲线，简化了非等间隔数据的求导过程，易于理解，容易实现。
发明内容
本发明的目的在于提供一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法。
第一方面，本发明实施例提供一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法，包括：
获取荧光定量熔解实验中非等间隔荧光强度数据；
确定最优滑动窗口大小以及多项式次数；
基于所述最优滑动窗口大小和多项式次数，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵；
基于导数矩阵，计算各滑动窗口内各采样点对应的负导数数值，获得平滑后的熔解曲线。
在可选的实施方式中，获取荧光定量熔解实验中非等间隔荧光强度数据的步骤，包括：
获取荧光定量之后扩增子熔解实验中非等间隔的初始荧光强度数据；
针对初始荧光强度数据，去除荧光强度均值，得到非等间隔荧光强度数据。
在可选的实施方式中，确定最优滑动窗口大小以及多项式次数的步骤，包括：
基于类Savitzky-Golay拟合方法，确定荧光强度曲线最优滑动窗口大小及多项式次数。
在可选的实施方式中，基于类Savitzky-Golay拟合方法，确定荧光强度曲线最优滑动窗口大小及多项式次数的步骤，包括：
基于类Savitzky-Golay拟合方法，获得不同滑动窗口大小和多项式次数下的荧光强度拟合曲线，计算多个决定系数；
根据多个决定系数中最大决定系数，确定荧光强度曲线最优滑动窗口大小和多项式次数。
在可选的实施方式中，基于所述最优滑动窗口大小和多项式次数，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵的步骤，包括：
基于最优滑动窗口大小和多项式次数，利用最小二乘法对各滑动窗口内荧光强度数据进行多项式拟合，获得多项式系数矩阵；
对拟合多项式进行求导，获得导数矩阵。
在可选的实施方式中，基于导数矩阵，计算各滑动窗口内各采样点对应的负导数数值，获得平滑后的熔解曲线的步骤，包括：
基于导数矩阵求得的导数数据确定负导数数值；
取各滑动窗口内中心值对应的负导数数值为各滑动窗口下中心采样温度对应的熔解曲线数值。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本申请实施例提供的一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法流程示意图；
图2为本申请实施例获取的原始荧光强度数据图。
图3为本申请实施例平滑前后的熔解曲线图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。
下面结合附图，对本发明的一些实施方式作详细说明。在不冲突的情况下，下述的实施例及实施例中的特征可以相互组合。
图1为本发明实施例提供的一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法流程示意图。如图1所示，该方法应用于计算机设备，该方法可以包括如下步骤：
S110，获取荧光定量熔解实验中非等间隔荧光强度数据；
其中，可以获取荧光定量之后扩增子熔解实验中非等间隔的初始荧光强度数据；针对初始荧光强度数据，去除荧光强度均值，得到非等间隔荧光强度数据。即在原始荧光强度的基础上减去均值。
其中，在获取非等间隔荧光强度数据时，按照实际需要，记录下不等间距温度下的荧光强度数据即可，不需刻意保证等温度间隔。
本实施例采用某阳性实验样本，熔解段反应程序为：95℃ 15s，65℃ 1min，97℃15s，采样温度间隔 0.1~0.2℃不等。如图2所示，选取该阳性实验样本熔解段荧光强度数据，针对该初始荧光强度数据，去除荧光强度均值。
S120，基于类Savitzky-Golay拟合方法，确定最优滑动窗口大小以及多项式次数；
作为一个示例，基于类Savitzky-Golay拟合方法，获得不同滑动窗口大小和多项式次数下的荧光强度拟合曲线，计算多个决定系数；根据多个决定系数中最大决定系数，确定荧光强度曲线最优滑动窗口大小和多项式次数。
按照Savitzky-Golay拟合原理，在不同滑动窗口（注：根据 Savitzky-Golay拟合原理，窗口大小为奇数）、不同多项式阶次下对不等间距下的荧光强度曲线进行拟合，对整体平滑过后的拟合曲线计算决定系数 R2：
Figure 429150DEST_PATH_IMAGE001
其中，
Figure 283973DEST_PATH_IMAGE002
表示实际的荧光强度值，
Figure 798131DEST_PATH_IMAGE003
表示拟合值，
Figure 380291DEST_PATH_IMAGE004
表示实际值的平均值。
决定系数最大值对应的滑动窗口大小、多项式阶次分别为最优的滑动窗口大小及多项式阶次，此时的荧光强度拟合曲线即为最优的 Savitzky-Golay 拟合结果。
本实施例中，我们分别计算了窗口大小5、7、9、11、13及多项式阶次1、2、3次下的决定系数R2，如表1所示。
表1 不同窗口大小及多项式阶次下的决定系数R2
Figure 866767DEST_PATH_IMAGE006
由表1可见，决定系数R2最大值为0.99，与决定系数R2最大值0.99对应的拟合窗口大小为9，多项式阶次为2。此即为最优的滑动窗口大小和多项式阶次。
S130，基于最优滑动窗口大小和多项式次数，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵；
可以基于最优滑动窗口大小和多项式次数，利用最小二乘法对各滑动窗口内荧光强度数据进行多项式拟合，获得多项式系数矩阵；对拟合多项式进行求导，获得导数矩阵。
此处的窗口大小、多项式次数即为根据最优的类Savitzky-Golay 拟合结果获得的窗口大小及多项式阶次，再利用最小二乘法进行拟合、获得多项式系数矩阵。本实施例中基于步骤S120计算得到的最优滑动窗口大小9及多项式阶次2，以中心采样点温度65.8℃对应的滑动窗口内二次多项式拟合函数为例：y =1.7939×107-5.3996×105x +4.0804×103x2，其系数矩阵即为：[1.7939×107, -5.3996×105, 4.0804×103]。
因为是多项式拟合，求取导数矩阵是很容易的。实施例中求取的导数矩阵为：[-5.3996×105 ,8.1608×103]，对应的一次多项式函数为：y =-5.3996×105+8.1608×103x；
依此方法，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵。
S140，基于导数矩阵，计算各滑动窗口内各采样点对应的负导数数值，获得平滑后的熔解曲线。
可以基于导数数据确定负导数数值；取各滑动窗口内中心值对应的负导数数值为各滑动窗口下中心采样温度对应的熔解曲线数值。
例如，将各滑动窗口内各温度代入导数矩阵对应的多项式函数，即可求得各采样点对应的导数数据，再乘“-1”，即得负导数数值；取窗口内中心值对应的负导数数值为各滑动窗口下中心采样温度对应的熔解曲线数值，综合即可得整体熔解曲线。本实施例中，以滑动窗口中心采样点温度65.8℃对应的熔解曲线数值为例，其熔解曲线数值为：y =-1×(-5.3996×105+8.1608×103×65.8)=2979.36。
图3为本实施例平滑前后的熔解曲线（注：平滑前为采用直接差分法计算结果）。由图3可见，平滑效果明显，熔点峰对应Tm值约为86℃。
本申请的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备执行本申请各个实施例移动控制方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器（Read-Only Memory，简称ROM）、随机存取存储器（Random Access Memory，简称RAM）、磁碟或者光盘等各种可以存储程序代码的介质。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释，此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
最后应说明的是：以上实施例，仅为本申请的具体实施方式，用以说明本申请的技术方案，而非对其限制，本申请的保护范围并不局限于此，尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本申请实施例技术方案的范围。都应涵盖在本申请的保护范围之内。

Claims (4)

1.一种非等间隔采样类Savitzky-Golay求导的熔解曲线平滑方法，其特征在于，包括：
获取荧光定量熔解实验中非等间隔荧光强度数据；
基于类Savitzky-Golay拟合方法，获得不同滑动窗口大小和多项式次数下的荧光强度拟合曲线，计算多个决定系数，根据所述多个决定系数中最大决定系数，确定最优滑动窗口大小以及多项式次数；
基于所述最优滑动窗口大小和多项式次数，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵；
基于所述导数矩阵，计算各滑动窗口内各采样点对应的负导数数值，获得平滑后的熔解曲线。
2.根据权利要求1所述的方法，其特征在于，获取荧光定量熔解实验中非等间隔荧光强度数据的步骤，包括：
获取荧光定量之后扩增子熔解实验中非等间隔的初始荧光强度数据；
针对所述初始荧光强度数据，去除荧光强度均值，得到所述非等间隔荧光强度数据。
3.根据权利要求1所述的方法，其特征在于，基于所述最优滑动窗口大小和多项式次数，对各滑动窗口内荧光强度数据进行多项式拟合，确定导数矩阵的步骤，包括：
基于所述最优滑动窗口大小和多项式次数，利用最小二乘法对各滑动窗口内荧光强度数据进行多项式拟合，获得多项式系数矩阵；
对拟合多项式进行求导，获得导数矩阵。
4.根据权利要求1所述的方法，其特征在于，基于导数矩阵，计算各滑动窗口内各采样点对应的负导数数值，获得平滑后的熔解曲线的步骤，包括：
基于所述导数矩阵求得的导数数据确定负导数数值；
取各滑动窗口内中心值对应的负导数数值为各所述滑动窗口下中心采样温度对应的熔解曲线数值。
CN202110288235.5A 2020-08-25 2021-03-18 一种熔解曲线平滑方法 Active CN112685686B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202010867905.4A CN111931118A (zh) 2020-08-25 2020-08-25 荧光定量熔解实验中熔解曲线测量方法
CN2020108679054 2020-08-25

Publications (2)

Publication Number Publication Date
CN112685686A CN112685686A (zh) 2021-04-20
CN112685686B true CN112685686B (zh) 2021-08-06

Family

ID=73304478

Family Applications (2)

Application Number Title Priority Date Filing Date
CN202010867905.4A Pending CN111931118A (zh) 2020-08-25 2020-08-25 荧光定量熔解实验中熔解曲线测量方法
CN202110288235.5A Active CN112685686B (zh) 2020-08-25 2021-03-18 一种熔解曲线平滑方法

Family Applications Before (1)

Application Number Title Priority Date Filing Date
CN202010867905.4A Pending CN111931118A (zh) 2020-08-25 2020-08-25 荧光定量熔解实验中熔解曲线测量方法

Country Status (1)

Country Link
CN (2) CN111931118A (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112652363A (zh) * 2020-12-30 2021-04-13 杭州博日科技股份有限公司 熔解曲线异常值处理方法、装置以及电子设备
CN113641667B (zh) * 2021-08-12 2022-05-20 深圳市润迅通投资有限公司 一种分布式大数据采集平台的数据异常监控系统及方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102187215A (zh) * 2008-09-09 2011-09-14 生物辐射实验室股份有限公司 基于回归的多级pcr分析系统
CN108287137A (zh) * 2017-12-22 2018-07-17 必欧瀚生物技术（合肥）有限公司 一种基于分段多项式拟合的基线校正方法
CN109145873A (zh) * 2018-09-27 2019-01-04 广东工业大学 基于遗传算法的光谱高斯峰特征提取算法
CN112331266A (zh) * 2020-11-20 2021-02-05 安图实验仪器(郑州)有限公司 消除pcr荧光基线期波动的方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8374795B2 (en) * 2008-05-13 2013-02-12 Roche Molecular Systems, Inc. Systems and methods for step discontinuity removal in real-time PCR fluorescence data
CN111089856B (zh) * 2019-12-26 2021-05-14 厦门大学 一种拉曼光谱弱信号提取的后处理方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102187215A (zh) * 2008-09-09 2011-09-14 生物辐射实验室股份有限公司 基于回归的多级pcr分析系统
CN108287137A (zh) * 2017-12-22 2018-07-17 必欧瀚生物技术（合肥）有限公司 一种基于分段多项式拟合的基线校正方法
CN109145873A (zh) * 2018-09-27 2019-01-04 广东工业大学 基于遗传算法的光谱高斯峰特征提取算法
CN112331266A (zh) * 2020-11-20 2021-02-05 安图实验仪器(郑州)有限公司 消除pcr荧光基线期波动的方法

Also Published As

Publication number Publication date
CN112685686A (zh) 2021-04-20
CN111931118A (zh) 2020-11-13

Similar Documents

Publication Publication Date Title
CN112685686B (zh) 一种熔解曲线平滑方法
US7720611B2 (en) Baselining amplification data
CN112342282B (zh) 荧光定量的指标确定方法
JP6190419B2 (ja) 情報システムを用いて反応を解析する方法およびシステム
EP1804172B1 (en) PCR elbow determination using curvature analysis of a double sigmoid
EP2107470B1 (en) PCR elbow determination using quadratic test for curvature analysis of a double sigmoid
JP5166608B2 (ja) 多段階回帰ベースｐｃｒ解析システム
US7668663B2 (en) Levenberg-Marquardt outlier spike removal method
JP4885934B2 (ja) リアルタイム核酸増幅データから初期核酸濃度を定量化する方法
US7991562B2 (en) PCR elbow determination using curvature analysis of a double sigmoid
CA2571399C (en) Temperature step correction with double sigmoid levenberg-marquardt and robust linear regression
JP2012519002A (ja) 融解曲線クラスタ化によるｓｎｐ検出
EP2120154A2 (en) Systems and methods for step discountinuity removal in real-time PCR fluorescence data
JP2006068011A (ja) リアルタイム核酸増幅データから初期核酸濃度を定量化する方法
CN110890127B (zh) 酿酒酵母dna复制起始区域识别方法
CN115545123B (zh) 熔解曲线优化方法、装置、电子设备及存储介质
EP1880334A2 (en) Identifying statistically linear data
JP5313663B2 (ja) 増幅データのベースライン化
JP5020728B2 (ja) 品質管理装置、方法及びプログラム
EP1880203A2 (en) Baselining amplification data
JP7345811B2 (ja) 自動閾値化を用いたデータ処理方法及びシステム
US10032001B2 (en) Methods and systems for identifying the quantitation cycle for a PCR amplification reaction
CN115859007B (zh) 一种石化仪表采样数据滑窗约束容错滤波降噪方法和装置
JP2008521104A5 (zh)
CN112652363A (zh) 熔解曲线异常值处理方法、装置以及电子设备

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
GR01 Patent grant
GR01 Patent grant