浅谈PCR引物设计中的吉布斯自由能

写在前面的
最近因为手里一个项目的需求，我开始接触PCR引物设计。因为脑子里的分子生物学基础还在，我误以为引物设计是一个很简单的过程。然而，当我真正开始学习引物设计后，我意识到关于引物设计有太多太多的底层问题被大多数人忽略了。调研了几篇文献后，我越发觉得引物设计是一项非常重要的技术，引物设计不仅与分子生物学、遗传学、热力学等息息相关，适当衍生后，我发现引物设计中涉及的很多思想都与计算机科学、数学等有密不可分的关系。

什么是PCR？什么引物？

        根据yourgenome名词解释PCR (Polymerase Chain Reaction) is a common tool used in medical and biological research labs. It is used in the early stages of processing DNA for sequencing, for detecting the presence or absence of a gene to help identify pathogens during infection, and when generating forensic DNA profiles from tiny samples of DNA.（https://www.yourgenome.org/）
        PCR过程中，待检测的模板DNA会随着温度变化被DNA聚合酶及能与模板互补配对的引物以扩增从而产生大量模板。因此常被用于测序、检测等实验中。因为PCR技术的出现，生命科学的数据产量正式迎来了一次革命性的提升。

什么是吉布斯自由能

        吉布斯自由能(Gibbs free energy，缩写G)又称为吉布斯函数，是在热力学或热化学中用于判断系统或化学反应进行的方向而引入的状态函数。因此，G被定义为

其中

U: 系统内能
T: 温度，以开尔文（K）为单位
S: 系统的熵
p: 系统压强
V: 系统体积
H: 系统的焓

        根据吉布斯函数及热力学定律，我们可以判断，在化学反应中可以根据吉布斯自由能变化量ΔG来判断化学反应进行方向。

其中

ΔG > 0：化学反应向着逆反应方向进行
ΔG < 0：化学反应向着正反应方向进行
ΔG = 0: 反应平衡

PCR中的ΔG

        PCR本质上也是一种化学反应。如果常做PCR的朋友们应该有很深刻的映像：当一对引物被primer程序设计好后，都会给出多个ΔG。例如下图中我随便找的一对引物，在表格中各自有一个ΔG，这个ΔG分别对应它们的最小ΔG，代表的是当前引物与模板上这段引物对应的正确位置的结合能力。（因为ΔG是反应化学反应的方向的，引物与模板的结合与不结合就是两个对立的化学反应方向）因此，ΔG小于0，代表引物在PCR过程中能够与模板结合。但是只有这一个ΔG远远不够。因为在真实的PCR反应系统中，还有其他变量与模板的特征片段争夺这段引物。

        例如，模板上其他的具有和引物局部相似的片段、另一条引物和这条引物的局部相似区（如果存在），引物自身与自身的配对（如上图右下角）。这些非模板的特征片段与对应引物的结合在分子生物学实验中通常被称为非特征/特异匹配。因此，计算这些片段与引物的ΔG可以表征引物的非特异结合能力，这里将模板特征片段与对应引物的吉布斯自由能变定义为ΔGs，其余非特征片段与引物的吉布斯自由能变定义为ΔG1，ΔG2等。一般情况下，只有上下游引物的吉布斯自由能满足如下公式，才能被认为这一对引物在实际的PCR实验中存在很小可能性发生非特异性扩增。

ΔG是如何被计算的

        事实上，当我在学习引物设计的时候经常看见一些推送、帖子写着引物的发卡结构、二聚体对应的ΔG需要大于-10，最好大于-7.但是几乎没有相关推送介绍这些ΔG是如何计算的。写到这里，可能有人会问我，你现在在这里说怎么算的不知道，那你前面的公式写的是啥？的确，前面的公式可以算。但是用那个公式的前提是你必须充分了解你要进行的PCR反应体系中各个状态函数对应的参数、常数都非常清楚。你才能用这个最原始的公式自己计算。然而，现在的PCR反应中大量的仪器、试剂耗材等都被大型公司垄断，这里面的成分你很难知道。甚至你做实验的屋子温度、压强你也不一定知道。那么，还有什么方法可以计算ΔG呢？从理论上说，引物是由一个一个碱基组成，如果我能知道每个碱基或者每单位长度个碱基的ΔG，那么岂不是能算出总的ΔG？事实上，科学家的确是这么干的。该领域的鼻祖级元老告诉了我他们当年搞出了下面的公式：

其中

ΔGi(total):一对引物对对应片段的总ΔG；
ΔGj:每个二临近碱基的ΔG，元老们测算出来的；
n:二临近碱基的组合数；
ΔG(init):初始吉布斯自由能常数，元老们测算出来的；
ΔG(sym):系统矫正ΔG常数，元老们测算出来的，不同系统的对应常数不同。

其中每一个二临近碱基组合的ΔG如下表所示

例如两段序列互补配对,那么吉布斯自由能计算如下

                    5' C-G-T-T-G-A 3'
                       | | | | | |
                    3’ G-C-A-A-C-T 5'

        但是，文献里说了Predictions are particularly poor for oligonucleotides shorter than 8 bp。用中文解释就是上面的模型对于核苷酸数量短于8bp的序列预测效果差。那么问题来了，如果是引物二聚体如下面所示，那么ΔG是怎么算出来的呢？如果尊敬的读者您知道该如何计算，欢迎向我的公众号私信计算方法，非常感谢。（因为我的公众号没有留言功能）

               5' TGTTTGAAGAGCAGAATCAATGGAG 3'
                       ||||| |  ||    |||       
               3'      CTTCTAGGTTTTAAAACCCTCTACC 5'

我的感悟

        在这一个星期的学习中，我非常清晰地认识到基础分子生物学实验背后涉及的交叉学科知识太多太深了。而且在大量的文献阅读中我发现相关的研究和知识都是在1980年到1996年之间被美日两国科学家研究推进的。这段时间该领域的元老们不断提出计算方程、不断做实验改进各个参数、常数。于是，到目前为止，总共出现了至少5种版本的参数表。在我学习的这段时间内发现国内甚至是国外的帖子都相当少地关注过这些物理化学量是如何被计算的，只是记住了一些结论引物的发卡结构、二聚体对应的ΔG需要大于-10，最好大于-7。但是就我目前测试的不同参数表的结果给出，不同的引物设计软件（如primer5、Oligo7等）参考的参数表是不同的。不同的参数表对同一个引物及模板计算的ΔG差异甚至能达到-10左右。因此，一些结论性的话如引物的发卡结构、二聚体对应的ΔG需要大于-10，最好大于-7这些不一定是绝对的，你一定要了解清楚这些数据是基于哪套参数表下的。否则你设计的引物不一定是真正最优的。但是当前国内甚至是国外生命科学界很多人都秉持实用主义原则！认为已经被证明、发现的东西没有必要再花精力去理解是怎么被证明发现的。自己只需要记住结论即可。诚然，在当前成果导向的科研环境下，用最少的时间去做出最好的结果对于以科研为职业的人来说很重要。但是我个人认为基础学科、底层科学对于应用学科的推动有极大的作用。如果个人的底层科学知识过少，甚至会限制个人在科学研究中的上限。

写在后面的
生命科学中的底层代码是我考虑作为偏向科学哲学科普的专题。如果从数学、计算机科学的视角看生命科学的PCR，那么可以说PCR的发明之于生命科学的意义犹如牛顿和莱布尼茨发明了微积分，犹如香农提出了信息论，甚至是犹如贝尔提出了贝尔不等式对量子力学的意义。生命科学的数据产出因为PCR的发明而上升了一个数量级，这个时候生命科学与物理学、化学和计算机科学的联系更加紧密。从某种角度说，引物设计已经不单单利用算法挖掘最优字符串了，甚至引物设计中的部分算法在某种程度上推动了语义模型的发展。

参考文献

[1] SantaLucia J Jr. A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics. Proc Natl Acad Sci U S A. 1998 Feb 17;95(4):1460-5. doi: 10.1073/pnas.95.4.1460. PMID: 9465037; PMCID: PMC19045.
[2] SantaLucia J Jr, Allawi HT, Seneviratne PA. Improved nearest-neighbor parameters for predicting DNA duplex stability. Biochemistry. 1996 Mar 19;35(11):3555-62. doi: 10.1021/bi951907q. PMID: 8639506.