CN116130000B - 引物探针序列组合设计方法和装置 - Google Patents

引物探针序列组合设计方法和装置 Download PDF

Info

Publication number
CN116130000B
CN116130000B CN202211533134.0A CN202211533134A CN116130000B CN 116130000 B CN116130000 B CN 116130000B CN 202211533134 A CN202211533134 A CN 202211533134A CN 116130000 B CN116130000 B CN 116130000B
Authority
CN
China
Prior art keywords
primer
primer probe
probe sequence
score
sequence combination
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
CN202211533134.0A
Other languages
English (en)
Other versions
CN116130000A (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.)
Wuhan Aimisen Life Technology Co ltd
Original Assignee
Wuhan Aimisen Life 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 Wuhan Aimisen Life Technology Co ltd filed Critical Wuhan Aimisen Life Technology Co ltd
Priority to CN202211533134.0A priority Critical patent/CN116130000B/zh
Publication of CN116130000A publication Critical patent/CN116130000A/zh
Application granted granted Critical
Publication of CN116130000B publication Critical patent/CN116130000B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • G16B25/20—Polymerase chain reaction [PCR]; Primer or probe design; Probe optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/20—Sequence assembly
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Molecular Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • Epidemiology (AREA)
  • Public Health (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Genetics & Genomics (AREA)
  • Databases & Information Systems (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioethics (AREA)
  • Analytical Chemistry (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本申请涉及一种引物探针序列组合设计方法，包括以下步骤：获取基于目标模板序列生成的引物探针序列组合集，引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合，每个引物探针序列组合用相应的参数信息表示，参数信息用于确定相应的引物探针序列组合的序列特征；构建评分模型，评分模型用于根据引物探针序列组合的序列特征计算引物探针序列组合的评分；采用优化算法在引物探针序列组合集中搜索评分模型的最优解处理评分模型，根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。本申请能够有效筛选扩增性能高、特异性高的引物对和探针检测甲基化的DNA，节约筛选的人工成本和耗材成本。

Description

引物探针序列组合设计方法和装置
技术领域
本申请涉及生物信息学技术领域，具体而言，涉及一种引物探针序列组合设计方法和装置。
背景技术
DNA甲基化指在DNA甲基转移酶的作用下，以S-腺苷甲硫氨酸作为甲基供体，将一个甲基基团共价转移到核苷酸上，最常见的是将甲基基团转移到胞嘧啶碱基第5位的碳原子上，生成5-甲基胞嘧啶。在哺乳动物的个体发育过程中，DNA的甲基化在调控基因表达、形成基因组印记及X染色体失活中起着重要作用。DNA甲基化不仅能够维持正常的生物功能，其异常可导致多种疾病的发生，如生长发育、行为异常以及癌症的发生和发展。因此，对基因中的特异位点进行甲基化检测有利于深入了解和阐明甲基化的生物学功能，也有利于筛查、诊断和评估疾病及其进展情况。
甲基化特异性PCR是一种经典的检测基因特定位点甲基化水平的方法，其引物设计的主要方法为：DNA首先用亚硫酸氢盐处理，并针对DNA序列中富含C的区域设计引物，同时CpG二核苷酸位点尽量位于引物的3’端，以区分甲基化和未甲基化的模板。甲基化特异性PCR法需要电泳检测或桑格尔测序来分析结果，步骤较为繁琐。
利用甲基化特异性引物和TaqMan MGB探针来覆盖甲基化位点，通过荧光定量PCR的方法对甲基化模板/非甲基化模板进行定量或定性分析成为最常用的甲基化检测的替代方法。具体地，DNA模板经历亚硫酸氢盐转化后会形成新的DNA序列，原序列中未甲基化的胞嘧啶C会转变为尿嘧啶U且在随后的PCR过程中转变为胸腺嘧啶T，而甲基化的胞嘧啶C仍为胞嘧啶。如此以来，原来互补配对的双链DNA不再配对，成为两条单链，且DNA的复杂度降低，GC含量降低，AT含量增加，极大地增加了引物和探针的设计难度。
因此，通过荧光定量PCR进行甲基化检测需要对模板目标区域同时设计几对甚至几十对甲基化/非甲基化引物和探针，再分别进行荧光定量PCR实验测试每一对引物和探针的性能，从而筛选用于甲基化检测的引物和探针。然而，上述引物和探针筛选方法十分繁琐，且工作量巨大。
如何有效筛选荧光定量PCR的检测引物和探针是对甲基化模板/非甲基化模板进行定量或定性分析的难点。
发明内容
为了解决上述问题，有效筛选用于甲基化检测的荧光定量PCR引物和探针，本申请的第一目的在于提供一种引物探针序列组合的设计方法，包括以下步骤：
获取基于目标模板序列生成的引物探针序列组合集，引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合，每个引物探针序列组合用相应的参数信息表示，参数信息用于确定相应的引物探针序列组合的序列特征；
构建评分模型，评分模型用于根据引物探针序列组合的序列特征计算引物探针序列组合的评分；
采用优化算法处理评分模型，以获得引物探针序列组合集中评分满足算法终止条件的目标引物探针序列组合。
本申请的引物探针序列组合设计方法通过设置引物探针序列组合的参数信息，构建引物探针序列组合的评分模型，结合优化算法的处理，获得评分满足算法终止条件的目标引物探针序列组合，可以有效筛选扩增性能良好的引物探针序列组合，扩增效率可达90％以上，用于甲基化检测，实现对甲基化模板/非甲基化模板进行定量或定性分析，减少复杂的人工劳动，同时还可以节约试剂、耗材，节约成本。
在其中一个实施例中，引物探针序列组合的序列特征包括引物对序列特征、探针序列特征、引物探针相互作用的序列特征，根据引物探针序列组合的序列特征计算引物探针序列组合的评分具体包括：
根据引物对的序列特征确定引物探针序列组合的引物评分；
根据探针的序列特征确定引物探针序列组合的探针评分；
根据引物探针相互作用的序列特征确定引物探针序列组合的引物探针相互作用评分；
根据引物探针序列组合的引物评分、探针评分、引物探针相互作用评分计算引物探针序列组合对应的评分。
在其中一个实施例中，评分模型为：
其中，TS表示每个引物探针序列组合对应的评分，S1表示根据引物对序列特征确定的引物评分，S2表示根据探针序列特征确定的探针评分，S3表示根据引物探针相互作用的序列特征确定的引物探针相互作用评分，Wx表示Sx对应的权重，且
在其中一个实施例中，引物对序列特征包括扩增子长度、上下游引物Tm值差值的绝对值、上下游引物各自形成发夹结构的长度之和、上下游引物之间形成二聚体的长度和引物对的CG位点数量中的至少一种；及/或
探针序列特征包括探针的5’端是否以G碱基开头、探针长度和探针内CG位点数量中的至少一种；及/或
引物探针相互作用的序列特征包括探针分别与上下游引物形成二聚体的长度之和以及探针Tm值与上下游引物的Tm均值的差值中的至少一种。
在其中一个实施例中，S1的计算公式为：S1＝100*exp(-(Saβ1+Sbβ2+Scβ3+Sdβ4+Seβ5；
其中，
Sa为扩增子长度，β1为Sa的系数；
Sb为上下游引物Tm值的差值的绝对值，β2为Sb的系数；
Sc为上下游引物各自形成发夹结构长度之和，β3为Sc的系数；
Sd为上下游引物之间形成二聚体的长度，β4为Sd的系数；
Se为上下游引物的CG位点数量，β5为Se的系数；
及/或
S2的计算过程包括：
判断探针的5’端是否以G碱基开头；
若探针的5’端以G碱基开头，则S2的计算公式为：S2＝0；
若探针的5’端不以G碱基开头，则S2的计算公式为：S2＝100*exp(-(Sfβ6+Sgβ7))；
其中，
Sf为探针的长度，β6为Sf的系数；
Sg为探针内CG位点数量，β7为Sg的系数；
及/或
S3的计算公式为：S3＝100*exp(-k1*(d+1)/ΔT)；
其中，k1为系数，d表示探针分别与上下游引物形成二聚体的长度之和，ΔT为探针的Tm值与上下游引物的Tm均值之间的差值。
在其中一个实施例中，评分模型中的系数使用多元回归模型确定。
在其中一个实施例中，评分模型中，β1＝0.001，β2＝0.03，β3＝0.05，β4＝0.08，β5＝-0.1，β6＝0.8，β7＝-4，k1＝1.57，W1＝0.5，W2＝0.3，W3＝0.2。
在其中一个实施例中，参数信息包括上游引物位置、上游引物长度、下游引物位置、下游引物长度、探针位置、探针长度以及探针方向。
在其中一个实施例中，优化算法包括模拟退火算法、遗传算法、爬山算法、粒子群算法、蚁群算法。
在其中一个实施例中，优化算法为模拟退火算法。
在其中一个实施例中，采用优化算法在引物探针序列组合集中搜索评分模型的最优解，最优解包括满足优化算法终止条件的评分及其对应的参数信息，根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合具体包括：
获取评分模型的初始解，评分模型的初始解包括引物探针序列组合集中任意一个引物探针序列组合的评分及其对应的参数信息；
利用模拟退火算法中的Metropolis准则，重复获取评分模型的新解并基于评分模型的新解对评分模型的初始解进行更新，以获得评分模型满足算法终止条件的最优解；
根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
在其中一个实施例中，利用模拟退火算法中的Metropolis准则，重复获取评分模型的新解并基于评分模型的新解对评分模型的初始解进行更新，以获得评分模型满足算法终止条件的最优解包括：
S1：获取初始退火温度T0和终止退火温度为Tn，将评分模型的初始解作为评分模型的当前最优解；
S2：基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤，对评分模型的当前最优解进行更新；
S3：判断评分模型的当前最优解对应的引物探针序列组合的评分是否满足算法终止条件；
S4：若评分模型的当前最优解对应的引物探针序列组合的评分不满足算法终止条件，则根据公式Ti＝αTi-1降低当前退火温度，i为正整数；
S5：判断当前退火温度是否小于Tn，若当前退火温度不小于Tn，返回执行步骤S2；
S6：若评分模型的当前最优解对应的引物探针序列组合的评分满足算法终止条件，终止搜索评分模型的最优解。
在其中一个实施例中，基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤中，每次获取新解和接受新解的步骤具体包括：
获取评分模型的新解，比较评分模型的当前最优解对应的引物探针序列组合的评分TSk和新解对应的引物探针序列组合的评分TSk+1；
若TSk＜TSk+1，则接受新解作为评分模型的当前最优解，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解。
在其中一个实施例中，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解具体包括：
若TSk＞TSk+1，获取当前退火温度并根据公式计算接受概率P，其中，ΔTS为TSk+1与TSk的差值，Ti为当前退火温度；
当接受概率P小于概率阈值时，接受新解作为评分模型的当前最优解。
在其中一个实施例中，判断评分模型的当前最优解对应的引物探针序列组合的评分是否满足算法终止条件具体包括：
判断评分模型的当前最优解对应的引物探针序列组合的评分是否小于评分阈值；
若评分模型的当前最优解对应的引物探针序列组合的评分不小于评分阈值，则确定评分模型的当前最优解对应的引物探针序列组合的评分满足优化算法终止条件；
若评分模型的当前最优解对应的引物探针序列组合的评分小于评分阈值，则确定评分模型的当前最优解对应的引物探针序列组合的评分不满足优化算法终止条件。
本申请的第二目的在于提供一种根据上述方法设计的引物探针组合产品，引物探针组合产品包括对目标模板执行甲基化检测的引物和探针。
本申请的第三目的在于提供一种引物探针序列组合设计装置，包括：
引物探针序列组合集获取模块：用于获取基于目标模板序列生成的引物探针序列组合集，引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合，每个引物探针序列组合用相应的参数信息表示，参数信息用于确定相应的引物探针序列组合的序列特征；
评分模型构建模块：用于构建评分模型，评分模型用于根据引物探针序列组合的序列特征计算引物探针序列组合的评分；
目标引物探针序列组合输出模块：用于采用优化算法在引物探针序列组合集中搜索评分模型的最优解处理评分模型，最优解包括满足优化算法终止条件的评分及该评分对应的参数信息，根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
本申请的第四目的在于提供一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现上述方法的步骤。
本申请的第五目的在于提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法的步骤。
本申请的第六目的在于提供一种计算机程序产品，包括计算机程序，其特征在于，该计算机程序被处理器执行时实现上述方法的步骤。
附图说明
为了更清楚地说明本申请具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本申请的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本申请实施例提供的引物探针序列组合的设计方法的流程图；
图2为本申请实施例提供的模拟退火算法原理示意图；
图3为本申请实施例提供的引物探针序列组合设计装置的结构框图；
图4为本申请实施例1中采用模拟退火算法搜索评分模型最优解的评分动态变化过程示意图；
图5为本申请实施例1中核苷酸组合1中引物对扩增目的片段的熔解曲线；
图6为本申请实施例1中核苷酸组合2中引物对扩增目的片段的熔解曲线；
图7为本申请实施例1中核苷酸组合3中引物对扩增目的片段的熔解曲线；
图8为本申请实施例1中核苷酸组合4中引物对扩增目的片段的熔解曲线；
图9为本申请实施例1中核苷酸组合5中引物对扩增目的片段的熔解曲线；
图10为本申请实施例1中核苷酸组合6中引物对扩增目的片段的熔解曲线；
图11为本申请实施例1中核苷酸组合7中引物对扩增目的片段的熔解曲线；
图12为本申请实施例1中核苷酸组合8中引物对扩增目的片段的熔解曲线；
图13为本申请实施例1中核苷酸组合9中引物对扩增目的片段的熔解曲线；
图14为本申请实施例1中核苷酸组合10中引物对扩增目的片段的熔解曲线。
具体实施方式
现将详细地提供本申请实施方式的参考，其一个或多个实例描述于下文。提供每一实例作为解释而非限制本申请。实际上，对本领域技术人员而言，显而易见的是，可以对本申请进行多种修改和变化而不背离本申请的范围或精神。例如，作为一个实施方式的部分而说明或描述的特征可以用于另一实施方式中，来产生更进一步的实施方式。
因此，旨在本申请覆盖落入所附权利要求的范围及其等同范围中的此类修改和变化。本申请的其它对象、特征和方面公开于以下详细描述中或从中是显而易见的。本领域普通技术人员应理解本讨论仅是示例性实施方式的描述，而非意在限制本申请更广阔的方面。
为了解决上述技术问题，如图1所示，本申请的第一方面提供了一种引物探针序列组合设计方法，包括以下步骤：
S10：获取基于目标模板序列生成的引物探针序列组合集，引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合，每个引物探针序列组合用相应的参数信息表示，参数信息用于确定相应的引物探针序列组合的序列特征；
具体地，目标模板是指用于甲基化检测的DNA模板。目标模板序列可以是人工预设的，也可以根据使用的目标探针对应的基因组位置信息提取基因组上的碱基序列得到，其中，目标探针可以是已有的DNA甲基化芯片上固定的探针。
可以理解的是，为了检测甲基化位点，需将目标模板转换成甲基化模板和非甲基化模板，并设计扩增甲基化模板的引物探针序列组合以及扩增非甲基化模板的引物探针序列组合。甲基化模板是将目标模板的DNA序列中所有的非CG位点中的C碱基替换为T碱基的DNA模板。非甲基化模板是将目标模板的DNA序列中所有的C碱基替换为T碱基的DNA模板，包括CG位点中的C碱基。此外，目标模板反向互补序列上所有的非CG位点中的C碱基替换为T碱基，得到反向互补的甲基化模板序列；与目标模板反向互补序列上所有的C碱基替换为T碱基，包括CG位点中的C碱基，得到反向互补的非甲基化模板序列。
确定模板以后，即可进行引物和探针的设计。具体地，在进行引物和探针的设计时，输入计划设计引物的DNA序列的上下游坐标，可以分别选择原始DNA序列、甲基化的DNA序列、非甲基化的DNA序列、原始DNA序列的反向互补序列、甲基化的反向互补序列、非甲基化的反向互补序列中的任意一种作为模板进行引物和探针的设计。在确定上述某一种DNA序列作为模板以后，将其称作上端序列，与上端序列反向互补的DNA序列称为下端序列。
一般而言，默认的引物长度为20个碱基，探针的长度为15个碱基，可以自行根据需求进行调节。在设计上游引物时，可将鼠标点击在上端序列的任意位置，获取鼠标点击位置的坐标，定义为F0，然后以F0为参考，取上端序列中F-19到F0之间20bp长度的序列作为上游引物，并显示在上端序列对应位置的下方。相应的，在设计下游引物时，可将鼠标点击在下端序列的任意位置，获取鼠标点击位置的坐标，定义为R0，然后以R0为参考，取下端序列中R0到R19之间20bp长度的序列作为下游引物，并显示在下端序列对应位置的上方。设计探针时可以有两种情况，一种是以上端序列为模板，点击上端序列与中间横线之间空白位置，获取鼠标点击位置的坐标，定义为P0，然后以P0位参考，取上端序列P-14到P0之间15bp长度的序列作为上游探针，并显示在上端序列对应位置的下方。另一种是以下端序列为模板，点击下端序列与中间横线之间空白位置，获取鼠标点击位置的坐标，定义为P1，然后以P1位参考，取下端序列P1到P15之间15bp长度的序列作为下游探针，并显示在下端序列对应位置的上方。
在一些实施方案中，完成引物设计以后，可根据实际需求在引物中引入突变碱基，引入突变碱基可以改变引物的退火温度(Tm值)，也可以消除引物自身形成的发卡结构以及引物二聚体。本申请引入碱基突变的原理如下所示：假设上下游引物长度分别为FL和RL，在上游引物中引入突变碱基的位置为{MF1，MF2，...，MFn}，MFn的取值范围为[1，FL]之间整数；在下游引物中引入突变的碱基位置为{MR1，MR2，...，MRm}，MRm的取值范围为[1，RL]之间整数。对于每一个位置的碱基，引入突变后都会有三种突变情况，这是因为DNA主要有T、C、G、A四种碱基组成。如MF1位置原有碱基为T，则引入突变的情况有C、G和A三种，加上原有的T，共有4种类型的上游引物。即每个位置在引入突变后，将对应4种类型的引物。当在上游引物中引入一个碱基突变时(F0位置不能引入突变)，将会产生4n-1种不同类型的上游引物，当在下游引物中引入一个碱基突变时(R0位置不能引入突变)，将会产生4m-1种不同类型的下游引物。当上下游引物中分别引入一个碱基突变时，将会有4(n+m-2)种不同类型的引物对。给出所有类型引物对的序列，并同时计算每个引物的Tm值和二聚体形成的情况，后续筛选步骤可以参考这些引物探针序列组合的参数信息挑选符合要求的引物。
一些实施方案中，每个引物探针序列组合的参数信息包括上游引物位置、上游引物长度、下游引物位置、下游引物长度、探针位置、探针长度和探针方向，上游引物位置用FP表示，上游引物长度用FL表示，下游引物位置用RP表示，下游引物长度用RL表示，探针位置用Mp表示、探针长度用ML表示以及探针方向用MO表示，探针方向用于表示探针是上游探针或下游探针。
当以上7个参数确定后，即可确定一组适用于荧光定量PCR的甲基化/非甲基化引物对和对应的检测探针，进而可根据引物对和对应的检测探针的序列确定每个引物探针序列组合的序列特征。为便于描述，这些参数统称为P，相应的参数空间称为当目标模板序列确定时，可利用计算机程序随机生成一组参数集合SP，设集合大小为n，则集合中每一个参数Pi对应一组引物对和探针，引物对和探针的序列可根据具体的参数信息从模板序列上获取。
一些实施方案中，引物探针序列组合的序列特征包括引物对序列特征、探针序列特征、引物探针相互作用的序列特征。具体地，引物对序列特征包括扩增子长度、上下游引物Tm值差值的绝对值、上下游引物各自形成发夹结构的长度之和、上下游引物之间形成二聚体的长度和引物对的CG位点数量中的至少一种，用于后续计算引物探针序列组合的评分；探针序列特征包括探针的5’端是否以G碱基开头、探针长度和探针内CG位点数量中的至少一种，用于后续计算引物探针序列组合的评分；引物探针相互作用的序列特征包括探针分别与上下游引物形成二聚体的长度之和以及探针Tm值与上下游引物的Tm均值的差值中的至少一种，用于后续计算引物探针序列组合的评分。
S20：构建评分模型，评分模型用于根据引物探针序列组合的序列特征计算引物探针序列组合的评分；
具体地，评分模型用于根据引物对序列特征、探针序列特征、引物探针相互作用的序列特征计算引物探针序列组合的评分。具体地，根据引物对的序列特征确定引物探针序列组合的引物评分，根据探针的序列特征确定引物探针序列组合的探针评分，根据引物探针相互作用的序列特征确定引物探针序列组合的引物探针相互作用评分，进一步，根据引物探针序列组合的引物评分、探针评分、引物探针相互作用评分计算引物探针序列组合的评分。
一些具体实施方案中，通过对引物探针序列组合的引物评分、探针评分、引物探针相互作用评分设置相应的权重，根据引物探针序列组合的引物评分、探针评分、引物探针相互作用评分的加权确定引物探针序列组合的评分。相应的，用于计算引物探针序列组合评分的评分模型可表示为：
其中，TS表示引物探针序列组合的评分，S1表示根据引物对序列特征确定的引物评分，S2表示根据探针序列特征确定的探针评分，S3表示根据引物探针相互作用的序列特征确定的引物探针相互作用评分，Wx表示Sx对应的权重，且
一些实施例中，S1的计算公式为：S1＝100*exp(-(Saβ1+Sbβ2+Scβ3+Sdβ4+Seβ5))；
其中，
Sa为扩增子长度，β1为Sa的系数；
Sb为上下游引物Tm值的差值的绝对值，β2为Sb的系数；
Sc为上下游引物各自形成发夹结构长度之和，β3为Sc的系数；
Sd为上下游引物之间形成二聚体的长度，β4为Sd的系数；
Se为上下游引物的CG位点的数量，β5为Se的系数。
可以理解的是，通常扩增子长度越小则引物对越优，但最小为60bp，对于游离的cfDNA，其片段长度通常在小于200bp处富集，因此扩增子长度与S1成反比。PCR扩增时，上下游引物的Tm值越接近越好，对应的|ΔTm|值越小，因此|ΔTm|与S1成反比。引物自身形成的发卡结构会影响引物与模板的结合稳定性，发夹结构越长，影响越大，因此发夹结构长度与S1成反比。引物之间的二聚体长度对引物与模板的结合带来不利影响，二聚体越长，这种影响也越大，因此二聚体长度与S1也成反比。引物中CG位点提供了区分甲基化模板与非甲基化模板的信息，CG位点数量越多，甲基化模板与非甲基化模板也越容易被区分，表现为甲基化引物的高特异性扩增，因此CG位点数量与S1成正比。
一些实施例中，S2的计算过程包括：
判断探针的5’端是否以G碱基开头；
若探针的5’端以G碱基开头，则S2的计算公式为：S2＝0；
若探针的5’端不以G碱基开头，则S2的计算公式为：S2＝100*exp(-(Sfβ6+Sgβ7))；
其中，
Sf为探针的长度，β6为Sf的系数；
Sg为探针内CG位点数量，β7为Sg的系数。
一些实施例中，S3的计算公式为：S3＝100*exp(-k1*(d+1)/ΔT)；
其中，k1为系数，d表示探针分别与上下游引物形成二聚体的长度之和，ΔT为探针的Tm值与上下游引物的Tm均值之间的差值。
一些实施方案中，使用多元回归模型确定评分模型中的系数值和权重值。具体地，随机选择一个基因CpG岛内一段长度为200bp的DNA序列作为模板，使用计算机随机数生产的方法生成100对引物和探针的位置。针对每一对引物和探针，都可以根据其参数信息确定其序列特征，如S1中Sa，Sb，Sc，Sd，Se的值、S2中Sf，Sg的值、S3中(d+1)/ΔT的值。人工合成甲基化的目标序列，并将其克隆至载体上作为模板，再分别人工合成100对引物探针组合，使用荧光定量PCR法分别检测100对引物探针组合的扩增效率。最后，以扩增效率为因变量(Y)，以S1、S2和S3中确定的序列特征的值为自变量(X)，并加入3个权重因子W1、W2和W3，建立回归模型，回归模型的公式为：
其中ε为构建模型而引入的随机因子，防止模型过拟合。根据建立的模型及100对X和Y的值，求出各系数值和权重值具体为：β1＝0.001，β2＝0.03，β3＝0.05，β4＝0.08，β5＝-0.1，β6＝0.8，β7＝-4，k1＝1.57，W1＝0.5，W2＝0.3，W3＝0.2。
S30：采用优化算法在引物探针序列组合集中搜索评分模型的最优解，最优解包括满足优化算法终止条件的评分及该评分对应的参数信息，根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
一些实施方案中，优化算法包括模拟退火算法、遗传算法、爬山算法、粒子群算法、蚁群算法。
一些具体实施方案中，优化算法为模拟退火算法。
具体地，模拟退火算法(simulated annealing，SA)来源于固体退火原理，是一种基于概率的算法。固体退火原理是将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温度升高变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。
类似的，模拟退火算法的原理是：先从一个较高的初始温度出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，局部最优解能概率性地跳出并最终趋于全局最优。即在每个温度下进行一轮搜索，每轮搜索时对旧解添加随机扰动生成新解，并按一定规则接受新解。
如图2所示，模拟退火算法本质是双层循环，外层循环控制温度由高向低变化，温度计算公式中α为取值在[0，1]上的温度衰减系数，如0.95；内层循环中，温度固定，对旧解添加随机扰动得到新解并按一定规则接受新解。内层循环中某一固定退火温度下的总迭代次数称为马尔科夫链长度，如图2中的马尔科夫链的长度为1000。
本申请根据模拟退火算法的原理，外层循环控制温度由高向低变化，温度计算公式中α为取值在[0，1]上的温度衰减系数，内层循环中，温度固定，在引物探针序列组合集中搜索评分模型的新解，并按一定规则接受新解，以得到用于对目标模板执行甲基化检测的引物探针序列组合。
相应的，一些实施方案中，采用优化算法在引物探针序列组合集中搜索评分模型的最优解，最优解包括满足优化算法终止条件的评分及其对应的参数信息，根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合具体包括：
S301：获取评分模型的初始解，评分模型的初始解包括引物探针序列组合集中任意一个引物探针序列组合的评分及其对应的参数信息；
S302：利用模拟退火算法中的Metropolis准则，重复获取评分模型的新解并基于评分模型的新解对评分模型的初始解进行更新，以获得评分模型满足算法终止条件的最优解；
S303：根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
其中，Metropolis准则是一种用概率来接受新解的方式，而不是使用完全确定的规则，从而降低搜索评分模型最优解的计算量。
具体地，利用模拟退火算法中的Metropolis准则，重复获取评分模型的新解并基于评分模型的新解对评分模型的初始解进行更新，以获得评分模型满足算法终止条件的最优解包括：
S1：获取初始退火温度T0和终止退火温度为Tn，将评分模型的初始解作为评分模型的当前最优解；
S2：基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤，对评分模型的当前最优解进行更新；
S3：判断评分模型的当前最优解对应的引物探针序列组合的评分是否满足算法终止条件；
S4：若评分模型的当前最优解对应的引物探针序列组合的评分不满足算法终止条件，则根据公式Ti＝αTi-1降低当前退火温度，i为正整数；
S5：判断当前退火温度是否小于Tn，若当前退火温度不小于Tn，返回执行步骤S2；
S6：若评分模型的当前最优解对应的引物探针序列组合的评分满足算法终止条件，终止搜索评分模型的最优解。
具体地，引物探针序列组合集中任意一个引物探针组合的评分可以根据评分模型表示的函数关系以及引物探针序列组合对应的参数信息确定的引物探针组合序列特征计算得到。模拟退火算法中，初始温度T0的值可以自行设置，一些具体实施方案中，T0设置为100℃。马尔科夫链长度L其没有固定值，可以自行设置，但是普遍接受的范围是500-2000之间，一些具体实施方案中，L的值取1000。
当评分模型的当前最优解完成L次迭代后，若未得到满足算法终止条件的最优解，进入外层循环降低当前退火温度，再重复评分模型的当前最优解的L次迭代过程，直至得到满足算法终止条件的最优解。具体地，根据公式Ti＝αTi-1降低当前退火温度，其中，α为温度衰减系数，i为正整数。
另一些实施方案中，将当前退火温度小于Tn作为模拟退火算法的外层循环的终止条件，若当前退火温度小于Tn，则终止搜索评分的最优解。
进一步，基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤中，每次获取新解和接受新解的步骤具体包括：
获取评分模型的新解，比较评分模型的当前最优解对应的引物探针序列组合的评分TSk和新解对应的引物探针序列组合的评分TSk+1，若TSk＜TSk+1，则接受新解作为评分模型的当前最优解，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解。
其中，评分模型的新解是指引物探针组合集中任意一个不同于评分模型的已知解的解，评分模型的新解包括引物探针组合集中新的引物探针序列组合及其对应的参数信息，TSk+1表示评分模型的新解对应的引物探针序列组合的评分，TSk表示评分模型的当前最优解对应的引物探针序列组合的评分。
一些实施方案中，为了实现对评分模型当前最优解的L次迭代，基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤，对评分模型的当前最优解进行更新具体包括：
S21：令m＝0；
S22：获取评分模型的新解，比较评分模型的当前最优解对应的引物探针序列组合的评分TSk和新解对应的引物探针序列组合的评分TSk+1；
S23：若TSk＜TSk+1，则接受新解作为评分模型的当前最优解，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解；
S24：令m＝m+1，判断m是否大于L；
S25：若m不大于L，返回执行步骤S22；若m大于L，执行步骤S3。
其中，m表示评分模型的当前最优解的迭代次数，当迭代次数未达到L时，则继续执行获取评分模型的新解和接受新解的步骤，当迭代次数达到L时，则执行判断评分模型的当前最优解是否符合算法终止条件的步骤。
一些具体实施方案中，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解具体包括：
若TSk＞TSk+1，获取当前退火温度并根据公式计算接受概率P；其中，ΔTS为TSk+1与TSk的差值，Ti为当前退火温度；
当接受概率P小于概率阈值时，接受新解作为评分模型的当前最优解。
一些具体实施方案中，概率阈值取0.05。
一些实施方案中，判断当前最优解对应的引物探针序列组合的评分是否满足算法终止条件具体包括：
S31：判断当前最优解对应的引物探针序列组合的评分是否小于评分阈值；
S32：若当前最优解对应的引物探针序列组合的评分不小于评分阈值，则确定当前最优解对应的引物探针序列组合的评分满足优化算法终止条件；
S33：若当前最优解对应的引物探针序列组合的评分小于评分阈值，则确定当前最优解对应的引物探针序列组合的评分不满足优化算法终止条件。
具体地，通过模拟退火算法搜索评分模型的当前最优解，当评分模型的当前最优解达到评分阈值时，即终止模拟退火算法的搜索过程，并根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
一些具体实施方案中，评分阈值为85。
本申请的引物探针序列组合设计方法通过设置引物探针序列组合的参数信息，构建引物探针序列组合的评分模型，结合模拟退火算法，输出评分满足算法终止条件的引物探针序列组合，可以有效筛选扩增性能良好的引物探针序列组合，用于甲基化检测，对甲基化模板/非甲基化模板进行定量或定性分析，扩增效率可达90％以上，可以减少复杂的人工劳动，同时还可以节约试剂、耗材，节约成本。
本申请的第二方面还提供了一种根据上述方法设计的引物探针组合产品，引物探针组合产品包括对目标模板执行甲基化检测的引物和探针。可以理解的是，在得到本申请的引物探针序列组合后可以通过常规的合成方法合成引物和探针，得到本申请用于对目标模板执行甲基化检测的引物探针组合产品。
相应的，如图3所示，本申请的第三方面提供了一种引物探针序列组合设计装置，包括：
引物探针序列组合集获取模块100：用于获取基于目标模板序列生成的引物探针序列组合集，引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合，每个引物探针序列组合用相应的参数信息表示，参数信息用于确定相应的引物探针序列组合的序列特征；
评分模型构建模块200：用于构建评分模型，评分模型用于根据引物探针序列组合的序列特征计算引物探针序列组合的评分；
目标引物探针序列组合输出模块300：用于采用优化算法在引物探针序列组合集中搜索评分模型的最优解处理评分模型，最优解包括满足优化算法终止条件的评分及该评分对应的参数信息，根据最优解对应的参数信息输出对应的引物探针序列组合作为目标引物探针序列组合。
一些实施方案中，引物探针序列组合的序列特征包括引物对序列特征、探针序列特征、引物探针相互作用的序列特征，根据引物探针序列组合的序列特征计算引物探针序列组合的评分具体包括：
根据引物对的序列特征确定引物探针序列组合的引物评分；
根据探针的序列特征确定引物探针序列组合的探针评分；
根据引物探针相互作用的序列特征确定引物探针序列组合的引物探针相互作用评分；
根据引物探针序列组合的引物评分、探针评分、引物探针相互作用评分计算引物探针序列组合对应的评分。
一些实施方案中，评分模型为：
其中，TS表示每个引物探针序列组合对应的评分，S1表示根据引物对序列特征确定的引物评分，S2表示根据探针序列特征确定的探针评分，S3表示根据引物探针相互作用的序列特征确定的引物探针相互作用评分，Wx表示Sx对应的权重，且
一些实施方案中，引物对序列特征包括扩增子长度、上下游引物Tm值差值的绝对值、上下游引物各自形成发夹结构的长度之和、上下游引物之间形成二聚体的长度和引物对的CG位点数量中的至少一种；及/或
探针序列特征包括探针的5’端是否以G碱基开头、探针长度和探针内CG位点数量中的至少一种；及/或
引物探针相互作用的序列特征包括探针分别与上下游引物形成二聚体的长度之和以及探针Tm值与上下游引物的Tm均值的差值中的至少一种。
一些实施方案中，S1的计算公式为：S1＝100*exp(-(Saβ1+Sbβ2+Scβ3+Sdβ4+Seβ5))；
其中，
Sa为扩增子长度，β1为Sa的系数；
Sb为上下游引物Tm值的差值的绝对值，β2为Sb的系数；
Sc为上下游引物各自形成发夹结构长度之和，β3为Sc的系数；
Sd为上下游引物之间形成二聚体的长度，β4为Sd的系数；
Se为上下游引物的CG位点数量，β5为Se的系数；
及/或
S2的计算过程包括：
判断探针的5’端是否以G碱基开头；
若探针的5’端以G碱基开头，则S2的计算公式为：S2＝0；
若探针的5’端不以G碱基开头，则S2的计算公式为：S2＝100*exp(-(Sfβ6+Sgβ7))；
其中，
Sf为探针的长度，β6为Sf的系数；
Sg为探针内CG位点数量，β7为Sg的系数；
及/或
S3的计算公式为：S3＝100*exp(-k1*(d+1)/ΔT)；
其中，k1为系数，d表示探针分别与上下游引物形成二聚体的长度之和，ΔT为探针的Tm值与上下游引物的Tm均值之间的差值。
一些实施方案中，评分模型中的系数使用多元回归模型确定。
一些实施方案中，评分模型中，β1＝0.001，β2＝0.03，β3＝0.05，β4＝0.08，β5＝-0.1，β6＝0.8，β7＝-4，k1＝1.57，W1＝0.5，W2＝0.3，W3＝0.2。
一些实施方案中，参数信息包括上游引物位置、上游引物长度、下游引物位置、下游引物长度、探针位置、探针长度以及探针方向。
一些实施方案中，优化算法包括模拟退火算法、遗传算法、爬山算法、粒子群算法、蚁群算法。
一些实施方案中，目标引物探针序列组合输出模块300具体包括：
求解单元：用于获取评分模型的初始解，评分模型的初始解包括引物探针序列组合集中任意一个引物探针序列组合的评分及其对应的参数信息；
更新单元：用于利用模拟退火算法中的Metropolis准则，重复获取评分模型的新解并基于评分模型的新解对评分模型的初始解进行更新，以获得评分模型满足算法终止条件的最优解；
输出单元：用于根据最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
在其中一个实施例中，更新单元具体包括：
获取子单元：用于获取初始退火温度T0，将评分模型的初始解作为评分模型的当前最优解；
更新子单元：用于基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤，对评分模型的当前最优解进行更新；
第一判断单元：用于判断评分模型的当前最优解对应的引物探针序列组合的评分是否满足算法终止条件；
温度控制单元：用于在评分模型的当前最优解对应的引物探针序列组合的评分不满足算法终止条件时，根据公式Ti＝αTi-1降低当前退火温度，i为正整数；
第二判断单元：用于判断当前退火温度是否小于Tn，若当前退火温度不小于Tn，则触发更新子单元；
终止单元：用于当评分模型的当前最优解对应的引物探针序列组合的评分满足算法终止条件时，终止搜索评分模型的最优解。
一些实施方案中，基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤中，每次获取新解和接受新解的步骤具体包括：
获取评分模型的新解，比较评分模型的当前最优解对应的引物探针序列组合的评分TSk和新解对应的引物探针序列组合的评分TSk+1；
若TSk＜TSk+1，则接受新解作为评分模型的当前最优解，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解。
一些实施方案中，若TSk＞TSk+1，则以概率P接受新解作为评分模型的当前最优解具体包括：
若TSk＞TSk+1，获取当前退火温度并根据公式计算接受概率P，其中，ΔTS为TSk+1与TSk的差值，Ti为当前退火温度；
当接受概率P小于概率阈值时，接受新解作为评分模型的当前最优解。
一些实施方案中，第一判断单元具体包括：判断子单元：用于判断评分模型的当前最优解对应的引物探针序列组合的评分是否小于评分阈值；
第一确定子单元：用于在判断结果为评分模型的当前最优解对应的引物探针序列组合的评分不小于评分阈值时，确定评分模型的当前最优解对应的引物探针序列组合的评分满足优化算法终止条件；
第二确定子单元：用于在在判断结果为评分模型的当前最优解对应的引物探针序列组合的评分小于评分阈值时，确定评分模型的当前最优解对应的引物探针序列组合的评分不满足优化算法终止条件。
本申请还涉及一种计算机可读存储介质，计算机存储介质用于存储计算机指令、程序、代码集或指令集，当其在计算机上运行时，使得计算机执行以实现上述引物探针序列组合的设计方法。
计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括——但不限于——电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括——但不限于——无线、电线、光缆、RF等等，或者上述的任意合适的组合。
可以以一种或多种程序设计语言或其组合来编写用于执行本公开操作的计算机程序代码，程序设计语言包括面向对象的程序设计语言-诸如Java、Smalltalk、C++，还包括常规的过程式程序设计语言-诸如”C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)-连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
本申请还涉及一种计算机设备，包括：
一个或多个处理器；以及
存储装置，存储一个或多个程序，
当一个或多个程序被一个或多个处理器执行，使得一个或多个处理器实现上述引物探针序列组合的设计方法。
可选的，计算机设备还可以包括收发器。处理器和收发器相连，如通过总线相连。需要说明的是，实际应用中收发器不限于一个，该计算机设备的结构并不构成对本申请实施例的限定。
其中，处理器可以是CPU，通用处理器，DSP，ASIC，FPGA或者其他可编程逻辑器件、晶体管逻辑器件、硬件部件或者其任意组合。其可以实现或执行结合本申请公开内容所描述的各种示例性的逻辑方框，模块和电路。处理器也可以是实现计算功能的组合，例如包含一个或多个微处理器组合，DSP和微处理器的组合等。
总线可包括一通路，在上述组件之间传送信息。总线可以是PCI总线或EISA总线等。总线可以分为地址总线、数据总线、控制总线等。存储器可以是ROM或可存储静态信息和指令的其他类型的静态存储设备，RAM或者可存储信息和指令的其他类型的动态存储设备，也可以是EEPROM、CD-ROM或其他光盘存储、光碟存储(包括压缩光碟、激光碟、光碟、数字通用光碟、蓝光光碟等)、磁盘存储介质或者其他磁存储设备、或者能够用于携带或存储具有指令或数据结构形式的期望的程序代码并能够由计算机存取的任何其他介质，但不限于此。
本申请还涉及一种计算机程序产品，包括计算机程序，该计算机程序被处理器执行时实现上述引物探针序列组合设计方法。
需要说明的是，本申请计算机程序产品所涉及的用户信息(包括但不限于用户设备信息、用户个人信息等)和数据(包括但不限于用于分析的数据、存储的数据、展示的数据等)，均为经用户授权或者经过各方充分授权的信息和数据，且相关数据的收集、使用和处理需要遵守相关国家和地区的相关法律法规和标准。
下面将结合实施例对本申请的实施方案进行详细描述。
实施例1
以Illumina Human Methylation 450k探针cg04261408为例，该探针识别SDC2基因正链DNA Chr8：96494447位置的胞嘧啶，该位置上下游各600bp碱基序列的DNA位置为：Chr8：96493847-96495047。取探针cg04261408识别的胞嘧啶位点的上下游各600bp碱基序列(总长度1201bp)作为模板，采用本申请构建的评分模型对所设计的甲基化引物、探针进行评分，并筛选最优的甲基化引物探针的组合。经过约80000次循环，此时的引物、探针得分超过85分，最高可达89.18分。图4表示使用模拟退火算法搜索最优的引物、探针时评分的动态变化过程示意图。
根据评分结果随机选取10种引物探针序列组合，如表1所示，其中，得分为89.19分的引物探针序列组合共6种，得分大于85分而小于89.19分的引物探针序列组合共2种，得分低于85分的引物探针共2种。
表1
将Chr8：96493847-96495047正链DNA经亚硫酸氢盐转化后的序列人工合成，同时合成完全甲基化且经亚硫酸氢盐转化后的序列和完全未甲基化且经亚硫酸氢盐转化后的序列，并分别克隆至载体，形成甲基化人工合成质粒1和非甲基化人工合成质粒2。将103拷贝/微升的人工合成质粒1和103拷贝/微升的人工合成质粒2等体积混合，作为模板备用，阴性对照管的模板为超纯水。人工合成表1中所列出的10种引物探针序列组合备用。表1中列举的探针5’端的荧光报告基团均为FAM，3’端的荧光淬灭基团均为MGB。
首先，验证引物的特异性：在每个PCR反应管中加入天根SuperReal PreMix Plus(SYBR Green)(FP205-02)试剂，再加入制备的模板、一对SDC2基因的甲基化检测引物对和探针组合、超纯水等，进行荧光定量PCR反应。每对核苷酸组合共检测3个复孔。根据熔解曲线评估扩增产物中是否存在非特异性扩增，若熔解曲线导数图谱中存在多个杂乱的峰，表明出现非特异性扩增，若熔解曲线为单一的尖锐峰，表明引物的特异性较高，未出现非特异性扩增。核苷酸组合1～10中引物对扩增甲基化的目标片段的熔解曲线如图5-图14所示。由图5～图14可以看出，评分较低的核苷酸组合7～10中的引物对均出现非特异性扩增的峰，而核苷酸组合1～6中的引物对均展示出特异性扩增，且其仅扩增甲基化的质粒模板，而不扩增未甲基化的质粒模板。
然后，使用TaqMan扩增体系验证核苷酸组合的扩增效果，所用的热启动DNA聚合酶和PCR缓冲液购自Invitrogen(Cat：14966005)。将质粒模板1依次稀释至107、106、105、104、103、102、101、0拷贝/微升，分别以上述质粒作为模板，再分别使用表1的引物对和探针进行荧光定量PCR扩增，扩增结束后统计Ct值，进而计算扩增效率。使用不同的核苷酸组合扩增的扩增效率及线性相关系数如表2所示。
表2
核苷酸组合 扩增效率(％) 线性相关系数R2值
1 98.01 0.996
2 96.33 0.985
3 102.29 0.994
4 105.87 0.984
5 101.54 0.999
6 102.62 0.988
7 84.22 0.993
8 120.47 0.989
9 81.55 0.982
10 86.28 0.991
由表2可以看出，线性相关系数R2值均在0.98以上，表明本次扩增实验数据可靠性高。另外，当使用核苷酸组合7～10进行PCR扩增时，其扩增效率或低于90％，或大于110％，不满足PCR扩增时扩增效率应在90％～110％之间的要求，因此，软件评分低而不合格的引物探针在实际实验过程中同样展示出很差的扩增效果。相反，软件评分较高的引物对和探针组合，如核苷酸组合1～6，其扩增效率满足90％～110％的要求，具有良好的扩增性能。
由以上结果可以看出，利用本申请的引物探针序列组合评分模型搜索最佳的(即评分最高的)的引物对和探针作为甲基化检测引物和探针的方法是可行的，利用该方法设计引物和探针可以减少复杂的人工劳动，同时还可以节约试剂、耗材，进而节约人工成本和耗材成本。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (18)

1.一种引物探针序列组合设计方法，其特征在于，包括以下步骤：
获取基于目标模板序列生成的引物探针序列组合集，所述引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合；
基于引物探针序列组合的序列特征构建评分模型，所述评分模型用于根据引物探针序列组合集中每个引物探针序列组合的序列特征计算每个引物探针序列组合的评分；
采用优化算法在所述引物探针序列组合集中搜索所述评分模型的最优解，所述最优解包括满足所述优化算法终止条件的评分及其对应的引物探针序列组合，根据所述最优解输出相应的引物探针序列组合作为目标引物探针序列组合；
其中，所述引物探针序列组合的序列特征包括引物对序列特征、探针序列特征、引物探针相互作用的序列特征，所述根据引物探针序列组合集中每个引物探针序列组合的序列特征计算每个引物探针序列组合的评分具体包括：
根据所述引物对的序列特征确定引物探针序列组合的引物评分；
根据所述探针的序列特征确定引物探针序列组合的探针评分；
根据所述引物探针相互作用的序列特征确定引物探针序列组合的引物探针相互作用评分；
根据引物探针序列组合的引物评分、探针评分、引物探针相互作用评分计算引物探针序列组合的评分。
2.根据权利要求1所述的方法，其特征在于，所述评分模型为：
；
其中，表示每个引物探针序列组合对应的评分，/>表示根据引物对序列特征确定的引物评分，/>表示根据探针序列特征确定的探针评分，/>表示根据引物探针相互作用的序列特征确定的引物探针相互作用评分，/>表示/>对应的权重，且/>。
3.根据权利要求2所述的方法，其特征在于，所述引物对序列特征包括扩增子长度、上下游引物Tm值差值的绝对值、上下游引物各自形成发夹结构的长度之和、上下游引物之间形成二聚体的长度和引物对的CG位点数量中的至少一种；及/或
所述探针序列特征包括探针的5’端是否以G碱基开头、探针长度和探针内CG位点数量中的至少一种；及/或
所述引物探针相互作用的序列特征包括探针分别与上下游引物形成二聚体的长度之和以及探针Tm值与上下游引物的Tm均值的差值中的至少一种。
4.根据权利要求3所述的方法，其特征在于，S 1的计算公式为：；
其中，
为扩增子长度，/>为/>的系数；
为上下游引物Tm值的差值的绝对值，/>为/>的系数；
为上下游引物各自形成发夹结构长度之和，/>为/>的系数；
为上下游引物之间形成二聚体的长度，/>为的/>系数；
为上下游引物的CG位点数量，/>为/>的系数；
及/或
S2的计算过程包括：
判断探针的5’端是否以G碱基开头；
若探针的5’端以G碱基开头，则的计算公式为：/>；
若探针的5’端不以G碱基开头，则的计算公式为：；
其中，
为探针的长度，/>为/>的系数；
为探针内CG位点数量，/>为/>的系数；
及/或
所述的计算公式为：/>；
其中，k1为系数，d表示探针分别与上下游引物形成二聚体的长度之和，为探针的Tm值与上下游引物的Tm均值之间的差值。
5.根据权利要求4所述的方法，其特征在于，所述评分模型中的系数使用多元回归模型确定。
6.根据权利要求5所述的方法，其特征在于，所述评分模型中， 。
7.根据权利要求1~6任一项所述的方法，其特征在于，每个引物探针序列组合用相应的参数信息表示，所述参数信息包括上游引物位置、上游引物长度、下游引物位置、下游引物长度、探针位置、探针长度以及探针方向。
8.根据权利要求7所述的方法，其特征在于，所述优化算法包括模拟退火算法、遗传算法、爬山算法、粒子群算法、蚁群算法中的至少一种。
9.根据权利要求7所述的方法，其特征在于，所述优化算法为模拟退火算法。
10.根据权利要求9所述的方法，其特征在于，所述采用优化算法在所述引物探针序列组合集中搜索所述评分模型的最优解，所述最优解包括满足所述优化算法终止条件的评分及其对应的引物探针组合，根据所述最优解输出相应的引物探针序列组合作为目标引物探针序列组合具体包括：
获取所述评分模型的初始解，所述评分模型的初始解包括所述引物探针序列组合集中任意一个引物探针序列组合的评分及其对应的参数信息；
利用模拟退火算法中的Metropolis准则，重复获取所述评分模型的新解并基于所述评分模型的新解对所述评分模型的初始解进行更新，以获得所述评分模型满足算法终止条件的最优解；
根据所述最优解对应的参数信息输出相应的引物探针序列组合作为目标引物探针序列组合。
11.根据权利要求10所述的方法，其特征在于，利用模拟退火算法中的Metropolis准则，重复获取所述评分模型的新解并基于所述评分模型的新解对所述评分模型的初始解进行更新，以获得所述评分模型满足算法终止条件的最优解包括：
S1：获取初始退火温度T0，将所述评分模型的初始解作为所述评分模型的当前最优解；
S2：基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤，对所述评分模型的当前最优解进行更新；
S3：判断所述评分模型的当前最优解对应的引物探针序列组合的评分是否满足算法终止条件；
S4：若所述评分模型的当前最优解对应的引物探针序列组合的评分不满足算法终止条件，则根据公式降低当前退火温度，i为正整数；
S5：判断当前退火温度是否小于终止退火温度Tn，若当前退火温度不小于Tn，返回执行步骤S2；
S6：若所述评分模型的当前最优解对应的引物探针序列组合的评分满足算法终止条件，终止搜索所述评分模型的最优解。
12.根据权利要求11所述的方法，其特征在于，所述基于模拟退火算法中的Metropolis准则，在当前退火温度下重复L次获取新解和接受新解的步骤中，每次获取新解和接受新解的步骤具体包括：
获取所述评分模型的新解，比较所述评分模型的当前最优解对应的引物探针序列组合的评分和新解对应的引物探针序列组合的评分/>；
若，则接受新解作为所述评分模型的当前最优解，若/>，则以概率P接受新解作为所述评分模型的当前最优解。
13.根据权利要求12所述的方法，其特征在于，所述若，则以概率P接受新解作为所述评分模型的当前最优解具体包括：
若，获取当前退火温度并根据公式/>计算接受概率P，其中，为/>与/>的差值，/>为当前退火温度；
当接受概率P小于概率阈值时，接受新解作为所述评分模型的当前最优解。
14.根据权利要求10~13任一项所述的方法，其特征在于，判断评分模型的当前最优解对应的引物探针序列组合的评分是否满足算法终止条件具体包括：
判断评分模型的当前最优解对应的引物探针序列组合的评分是否小于评分阈值；
若评分模型的当前最优解对应的引物探针序列组合的评分不小于评分阈值，则确定评分模型的当前最优解对应的引物探针序列组合的评分满足优化算法终止条件；
若评分模型的当前最优解对应的引物探针序列组合的评分小于评分阈值，则确定评分模型的当前最优解对应的引物探针序列组合的评分不满足优化算法终止条件。
15.根据权利要求1~14任一项所述的方法设计的引物探针组合产品，所述引物探针组合产品包括用于对目标模板执行甲基化检测的引物和探针。
16.一种引物探针序列组合设计装置，其特征在于，包括：
引物探针序列组合集获取模块：用于获取基于目标模板序列生成的引物探针序列组合集，所述引物探针序列组合集包括至少一个用于对目标模板执行甲基化检测的引物探针序列组合；
评分模型构建模块：用于基于引物探针序列组合的序列特征构建评分模型，所述评分模型用于根据引物探针序列组合集中每个引物探针序列组合的序列特征计算每个引物探针序列组合的评分；
目标引物探针序列组合输出模块：用于采用优化算法在所述引物探针序列组合集中搜索所述评分模型的最优解，所述最优解包括满足所述优化算法终止条件的评分及其对应的引物探针组合，根据所述最优解输出相应的引物探针序列组合作为目标引物探针序列组合；
其中，所述引物探针序列组合的序列特征包括引物对序列特征、探针序列特征、引物探针相互作用的序列特征，所述根据引物探针序列组合集中每个引物探针序列组合的序列特征计算每个引物探针序列组合的评分具体包括：
根据所述引物对的序列特征确定引物探针序列组合的引物评分；
根据所述探针的序列特征确定引物探针序列组合的探针评分；
根据所述引物探针相互作用的序列特征确定引物探针序列组合的引物探针相互作用评分；
根据引物探针序列组合的引物评分、探针评分、引物探针相互作用评分计算引物探针序列组合的评分。
17.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1~14中任一项所述方法的步骤。
18.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1~14中任一项所述方法的步骤。
CN202211533134.0A 2022-12-02 2022-12-02 引物探针序列组合设计方法和装置 Active CN116130000B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211533134.0A CN116130000B (zh) 2022-12-02 2022-12-02 引物探针序列组合设计方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211533134.0A CN116130000B (zh) 2022-12-02 2022-12-02 引物探针序列组合设计方法和装置

Publications (2)

Publication Number Publication Date
CN116130000A CN116130000A (zh) 2023-05-16
CN116130000B true CN116130000B (zh) 2023-12-12

Family

ID=86303502

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211533134.0A Active CN116130000B (zh) 2022-12-02 2022-12-02 引物探针序列组合设计方法和装置

Country Status (1)

Country Link
CN (1) CN116130000B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107447004A (zh) * 2017-08-11 2017-12-08 北京呈诺医学科技有限公司 Dna甲基化pcr检测引物或探针的特异性检测方法
CN110951828A (zh) * 2016-04-29 2020-04-03 广州市康立明生物科技有限责任公司 一种引物和探针的设计方法
CN112779320A (zh) * 2020-12-04 2021-05-11 深圳市易基因科技有限公司 多区域dna甲基化检测探针设计及其检测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3828289A4 (en) * 2018-07-26 2022-06-08 Excellen Medical Technology Co., Ltd. METHOD AND KIT FOR IDENTIFYING THE CONDITION OF GASTRIC CANCER

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110951828A (zh) * 2016-04-29 2020-04-03 广州市康立明生物科技有限责任公司 一种引物和探针的设计方法
CN107447004A (zh) * 2017-08-11 2017-12-08 北京呈诺医学科技有限公司 Dna甲基化pcr检测引物或探针的特异性检测方法
CN112779320A (zh) * 2020-12-04 2021-05-11 深圳市易基因科技有限公司 多区域dna甲基化检测探针设计及其检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
ACustomized Simulated Annealing Suitable for Primer Design in Polymerase Chain Reaction Processes;Luciana Montera，et al.;《Simulated Annealing， Theory with Applications》;全文 *
SARS病毒再测序DNA微阵列的探针设计;谢建明, 方辉, 陆祖宏, 孙啸;生物技术(第01期);全文 *
仿生计算在生物信息学中的应用;宋杰;《中国优秀博硕士学位论文全文数据库（博士） 信息科技辑》;第I138-68页 *

Also Published As

Publication number Publication date
CN116130000A (zh) 2023-05-16

Similar Documents

Publication Publication Date Title
LaFleur et al. Automated model-predictive design of synthetic promoters to control transcriptional profiles in bacteria
JP7385686B2 (ja) 無細胞核酸の多重解像度分析のための方法
CN112639982A (zh) 使用神经网络调用倍性状态的方法和系统
KR102028375B1 (ko) 희귀 돌연변이 및 카피수 변이를 검출하기 위한 시스템 및 방법
Uricchio et al. Selection and explosive growth alter genetic architecture and hamper the detection of causal rare variants
JP2023022220A (ja) 遺伝子の変動の非侵襲的評価のための方法および処理
JP7237003B2 (ja) 遺伝子片の評価のための方法およびプロセス
TW202039860A (zh) 游離dna末端特徵
US20230002817A1 (en) Method and kit for detecting genome editing and application thereof
WO2016001760A2 (en) Determining mutation burden in circulating cell-free nucleic acid and associated risk of disease
CN110246543B (zh) 基于二代测序技术利用单样本检测拷贝数变异的方法和计算机系统
Masharing et al. ddRAD sequencing based genotyping of six indigenous dairy cattle breeds of India to infer existing genetic diversity and population structure
Sun et al. HS-MMGKG: a fast multi-objective harmony search algorithm for two-locus model detection in GWAS
Wei et al. CONY: A Bayesian procedure for detecting copy number variations from sequencing read depths
CN116130000B (zh) 引物探针序列组合设计方法和装置
He et al. Characterization and machine learning prediction of allele-specific DNA methylation
Choi et al. Quantification and sequencing of crossover recombinant molecules from Arabidopsis pollen DNA
JP2023506084A (ja) ゲノム瘢痕アッセイ及び関連する方法
Maeda et al. Development of a DNA barcode tagging method for monitoring dynamic changes in gene expression by using an ultra high-throughput sequencer
JP7470787B2 (ja) 単一試料からの腫瘍純度の推定
Sun et al. CRISPR-M: Predicting sgRNA off-target effect using a multi-view deep learning network
WO2022020346A1 (en) Cancer detection, monitoring, and reporting from sequencing cell-free dna
CN110462063B (zh) 一种基于测序数据的变异检测方法、装置和存储介质
CN114517223A (zh) 一种用于筛选snp位点的方法及其应用
Zhu et al. Genomic microsatellite characteristics analysis of Dysommaanguillare (Anguilliformes, Dysommidae), based on high-throughput sequencing technology

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