宏基因检测病原已有多篇专家共识陆续发表，管理规定相对成熟。除了宏基因技术，多重PCR扩增结合高通量测序因具有高性价比、分析速度快、灵敏度等优势在病原进行检测中具有重要地位。

设计合适的多重引物是多重PCR难点之一。primer3和SADDLE是两个设计多重引物的工具，两篇文献中均提及如下文献（2004年，IF=12.4）：

article.png

primer3基于该文献计算热力学数据。因此，小编探索该文献后而有此文。术业有专攻，小编也是浅尝辄止。本文主要目的仅仅是介绍相关概念，方便读者自行研究相关文献。本文也只涉及DNA部分，RNA部分未提。

生物概念

计算自由能涉及单个碱基化学结构，因此首先简要回顾生物相关概念，便于大家阅读相关文献时知其所云。

核酸和蛋白组成染色体：

chromosome.png

核苷酸结构如下：

nucletide.png

DNA二级结构包含三种螺旋，主要为BDNA螺旋：

dna_duplex.png

BDNA右手螺旋和Watson–Crick碱基对：

watson-crick.png

影响DNA结构稳定的因素主要有如下5种：

dna_structure.png

以上5种因素直接影响熔链温度（Tm）:

通过核酸的紫外吸收特性，可以计算出不同温度下未变性核酸的含量，从而得出Tm：

tm.png

自由能

专业讲解可参考华中科技大学64 生物能学（1）_哔哩哔哩_bilibili。

个人理解如下：碱基互补配对为放热反应，需要消化自由能。因此，delta G为负数。放热越多，系统越稳定，delta G越小。

NN模型和二级结构

通常通过NN模型计算双链DNA（duplex）的 $\Delta G$ （N表示为1个碱基）。

NN_model.png

<center>DOI: https://doi.org/10.1038/s41467-022-29500-4</center>

上图所示，用法类似于k-mer（k=2）。将每一个 $NN$ 的 $\Delta G$ 进行求和，同时包含“惩罚项”，即为最后DNA双链的 $\Delta G$ 。

注意：因为DNA右手螺旋性质， $AT/TA$ 和 $TA/AT$ 的二级结构分布稍有不同，因此 $\Delta G(AT/TA) \neq \Delta(TA/AT)$ 。详细解释可参考Nucleic Acid Structure 2 – Stability of DNA - YouTube.

DNA具有不同的二级结构，不同的二级结构计算 $\Delta G$ 方式不同：

dna_structure2.png

后文基于上图不同的二级结构分别计算 $\Delta G$ 。

Watson-Crick碱基对

$\Delta G$ 计算

Watson-Crick碱基对计算自由能公式如下：
$\Delta G^0_{37}(total)=\Delta G^0_{37}(initiation)+\Delta G^0_{37}(symmetry)+\sum \Delta G^0_{37}(stack)+\Delta G^0_{37}(AT \ terminal) \\$

delta g symmetry = +0.43 kcal/mol if the duplex is self-complementary

delta G symmetry = 0 if it is non-selfcomplementary.

10种堆积碱基（stack）自由能见下表：

NN_model_watson-crick.png

作者提供案例如下：

watson-crick-deltag.png

TABLE 1是基于 $37 ^{\circ}C$ 计算结果，不同温度下同一个stack $\Delta G$ 也不同，计算公式如下：
$\Delta G^0_T = \Delta H^0-T \Delta S^0$

焓变（ $\Delta H^0$ ）和熵变（ $\Delta S^0$ ）均为热力学相关概念，单位分别为 $cal/mol$ 和 $cal/K/mol$ ；
$T$ 为Kelvin温度。

$T_m$ 计算

$T_m$ 计算公式如下：
$T_m=\frac{\Delta H^0*1000}{\Delta S^0+R*\ln(\frac{C_T}{x})}-273.15$

CT is the total molar strand concentration,

R is the gas constant 1.9872 cal/K-mol

x equals 4 for nonself-complementary duplexes and equals 1 for self-complementary duplexes

此处 $\Delta H^0$ 应该指所有stack之和，即 $\Delta H^0_{NN}$ 。

作者提供案例如下：

For a nonself-complementary duplex with delta H =−43.5 kcal/mol , delta S = −122.5 e.u., and strand concentrations of 0.2 mM for each strand:
$T_m=\frac{-43.5*1000}{-122.5+1.9872*\ln(\frac{0.0004}{4})}-273.15=35.8^{\circ}C$

$T_m$ 计算公式有很多，最常见的怕是：
$T_m=2(A+T)+4(G+C)$
正如不同的算法有不同的适用条件，适用上述公式同样需要留意：

It is important to realize that this relationship holds only for short oligonucleotides (11–27 bases long) and is based on hybridization in 1 M Na+ (6× SSC), in the complete absence of any organic solvents. --Melting Temperature - an overview | ScienceDirect Topics.

更多 $T_m$ 计算公式可参见Melting Temperature - an overview | ScienceDirect Topics。

在多重PCR中， $T_m$ 精确度更为重要：

poor thermodynamic predictions may be tolerated for single target PCR primer design, because even if the predicted TM is 10 $^{\circ}C$ inaccurate, one has the luxury of experimentally optimizing the annealing temperature. In more complicated assays, such as multiplex PCR, however, all the amplifications must occur under the same conditions and inaccuracies in TM predictions result in poor primer designs that cause failed amplifications.

--doi: 10.1146/annurev.biophys.32.110601.141800

salt correction

前文已知，溶液中的离子浓度影响DNA结构的稳定性。实际生产中，溶液离子浓度各不同。因此为了提高上述公式的准确度，可对离子浓度（主要 $Na^+$ ）进行矫正：
$\Delta G^0_{37}[Na^+]=\Delta G^0_{37}[1M \ NaCl]-0.114*N/2*\ln[Na^+] \\ \Delta S^0[Na^+]=\Delta S^0[1M \ NaCl]+0.368*N/2*\ln[Na^+] \\$

N is the total number of phosphates in the duplex,

[Na+] is the total concentration of monovalent cations from all sources

$\Delta H$ 不进行矫正。

作者提供案例：

前文中，CGTTGA/GCAACT计算delta g=-5.35 kcal/mol。假设：

0.1 M NaCl, 10 mM sodium phosphate, pH 7 (gives a total of 0.115 M Na+ because at pH 7 there are 1.5 equivalents of sodium for each phosphate) gives
$\begin{aligned} \Delta G^0_{37}[0.115M \ Na^+] &=-5.35 kcal/mol - 0.114*10/2*\ln[0.115] \\ &=-4.12 kcal/mol \end{aligned}$

下文计算delta G时，也可以以此公式进行salt correction。

mismatch

internal mismatch

duplex内部单个碱基错配计算方式同于Watson-Crick，只是基于不同的数据库：

internal_mismatch_db.png

作者提供案例：

internal_mismatch_case.png

terminal mismatch

假设错配在末端1bp或2bp，如下图：

terminal_mismatch.png

这样的mismatch更多认为是terminal mismatch, 不是internal mismatch。teriminal mismatch使得duplex更稳定，即 $\Delta G$ 更小。化学反应往往趋向于更稳定的状态，因此，通常计算上图右侧的 $\Delta G$ 作为最终结果。

Dangling Ends

长链和短链互补配对时，长链多处的片段即Dangling Ends。Dangling Ends第一个碱基对结构稳定影响更大，因此通常只考虑第一个碱基。

Dangling Ends数据库如下：

dangling_ends.png

loop

loop分为internal loop, bulge loop和hairpin loop。

30bp范围内的loop delta g数据库如下：

loop.png

当loop大于30bp时，可通过以下公式推算（Jacobson-Stockmayer entropy extrapolation）：
$\Delta G^0_{37}(\text{loop-n})=\Delta G^0_{37}(\text{loop-x})+2.44*R*310.15*\ln(n/x)$

delta G (loop-n) is the free energy increment of a loop of length n,

delta G (loop-x) is the free-energy increment of the longest loop of length x for which there are experimental data

R is the gas constant

hairpin loops

$\Delta G$ 计算

hairpin loop长度不同， $\Delta G$ 计算公式也不同。

hairpin长度为3

$\Delta G^0_{37}(\text{total})=\Delta G^0_{37}(\text{hairpin of 3})+\Delta G^0_{37}(\text{triloop bonus})+\text{closing AT penalty}$

closing AT penalty = 0.5 kcal/mol

applied only to hairpin sequences that are closed by AT.

delta G (triloop bonus) values are given in the supplementary material. 查看附件补充

备注：2bp默认不形成hairpin loop.

hairpin长度为4

$\Delta G^0_{37}(\text{total})=\Delta G^0_{37}(\text{hairpin of 4})+\Delta G^0_{37}(\text{triloop bonus})+\Delta G^0_{37}(\text{terminal mismatch})$

delta G 37 (terminal mismatch) is the increment for terminal mismatches

作者提供案例如下。假设有一段序列C<u>GCAA</u>G，delta g:
$\Delta G^0_{37}=3.5-1.6-1.23=0.67 \ kcal/mol$

hairpin长度大于4

$\Delta G^0_{37}(\text{total})=\Delta G^0_{37}(\text{hairpin of N})++\Delta G^0_{37}(\text{terminal mismatch})$

$T_m$ 计算

前文可知，watson-crick Tm计算公式为：
$T_m=\frac{\Delta H^0*1000}{\Delta S^0+R*\ln(\frac{C_T}{x})}-273.15$
hairpin Tm计算公式稍有不同：
$T_m=\frac{\Delta H^0*1000}{\Delta S^0}-273.15$

Internal loops

公式为：
$\begin{aligned} \Delta G^0_{37}(\text{loop total}) &=\Delta G^0_{37}(\text{internal loop of N})+\Delta G^0_{37}(\text{asymmetry}) \\ &\quad +\Delta G^0_{37}(\text{left terminal mismatch}) \\ &\quad +\Delta G^0_{37}(\text{right terminal mismatch}) \\ \Delta G^0_{37}(\text{asymmetry}) &= |\text{length A}− \text{length B}| × 0.3 \ \text{kcal/mol} \end{aligned}$