论文多重PCR引物设计的软件开发.docx
《论文多重PCR引物设计的软件开发.docx》由会员分享，可在线阅读，更多相关《论文多重PCR引物设计的软件开发.docx（31页珍藏版）》请在冰豆网上搜索。
论文多重PCR引物设计的软件开发
多重PCR引物设计的软件开发
第一章引言
1.1背景知识简介
聚合酶链式反应（PolymeraseChainReaction,PCR）是体外扩增DNA的一种技术；其能够在短时间内根据极微量的模板序列扩增出大量特异性DNA片断，扩增过程类似于核裂变。
Mullis博士在1983年发明了该技术，并因此于1993年获得诺贝尔化学奖。
PCR技术是现代分子生物学中最有价值的技术之一[1]；例如：
人类基因组计划，亲子鉴定，对罪犯DNA的鉴别，都依赖于这一技术。
在聚合酶链式反应中，原始的DNA片段称为模板序列；待复制的DNA片段称为目标序列，它是模板序列的一部分；引导目标序列合成的寡核苷酸片段称为引物，其长度一般在16-27之间。
每个反应一般含20至40个循环，每次循环目标序列含量翻一倍；每个循环由三个主要步骤构成：
1.变性（90℃-96℃）：
双链DNA模板在热作用下，氢键断裂，形成单链DNA；2.退火（25℃-65℃）：
系统温度降低，引物与DNA模板结合，形成局部双链；3.延伸（70℃-75℃）：
在Taq酶的作用下，以dNTP为原料，从引物的5′端→3′端延伸，合成与模板互补的DNA链。
扩增一段目标序列需要两条引物，分别和目标序列的两端相匹配。
引物的作用类似于电路里的开关，引物设计不好会导致产物杂乱无章（即扩增失败），因此引物设计是PCR成功的关键之一。
引物设计需要综合考虑引物的退火温度，GC含量，引物长度，特异性等。
现有的PCR引物设计软件已经非常成熟，可以达到较好的效果,常用软件有PRIMER3,PRIMER5,OLIGO6.0等。
1.1.1多重PCR引物设计简介
多重聚合酶链式反应（多重PCR）是指在同一试管中同时进行多个PCR，从而一次性扩增多个目标片段，因此节省了大量的时间和金钱，具有巨大的经济和时间效应。
由于引物之间以及引物与模板序列之间的相互干扰，一方面导致了多重PCR引物设计考虑因素增加，另一方面极大地增加了计算的空间和时间复杂度。
同时多重PCR引物设计问题需要被当作一个整体考虑，不能被分割为多个PCR引物设计问题的简单组合[2]。
现尚无完善的理论解决这一问题。
综上所述，要实现多重PCR引物设计的自动化，在实践中需要解决两个问题：
一、多重PCR引物设计的约束条件与评价标准；二、载体算法的选择。
1.1.2遗传算法简介
遗传算法是一种随机优化算法，用于求解问题的全局最优解，尤其是对非线性问题的全局搜索和最优化。
该算法由美国密歇根大学的JohnHolland教授于80年代提出，其借鉴了达尔文优胜劣汰的思想，以及遗传过程中的染色体交叉与变异的概念。
虽然遗传算法尚未被完全从理论上证明，但遗传算法已得到比较成熟的发展，并广泛地被成功应用于各类实际问题；Schema定理部分证明了该算法的有效性[3]。
遗传算法主要步骤包括初始化，可行解的评估，交叉，变异，和收敛。
每个步骤根据问题的性质与规模不同，均有多种策略可以选择。
在实践中，要将遗传算法和问题结合需要解决如下几个问题：
一、建立问题的数学模型并表示为数据结构；二、初始化策略的选择；三、挑选双亲与交叉策略的选择；四、变异策略的选择；五、收敛条件的选择；六、交叉率，变异率，初始解数目的选择。
参见图1.1遗传算法流程图.
图1.1遗传算法流程图
1.2相关文献综述
由于多重PCR引物设计问题尚未在算法理论上得到完善解决，根据对文献的总结，一方面、国内对引物设计的自动化问题关注较少，鲜有文献记载；国外近几年已提出一些启发性的算法和相应的软件，例如MultiPLX[4][5]，但尚无宣称成熟的算法与软件。
另一方面、为了满足现实需要国外也有比较成熟的大型系统采用常规的算法，使用大量并行计算机以及较长的时间来为大规模多重PCR设计引物，虽然较好地解决了问题，但这需要以大量的时间和金钱为代价，同时也将应用限制在大规模PCR引物设计问题上。
例如日本的PrimerStation,使用了100个CPU（SunFire15K）来为人类全基因扩增问题设计引物，耗时为3个月[6]。
在调查总结文献的基础上，本研究选择了遗传算法作为解决问题的载体，具体思路将在第二章中进行阐述。
台湾国立中央大学与国立中山大学共有两篇论文运用了遗传算法进行了多重PCR引物设计的研究[7][8]；随后由李宗南教授总结后几位作者共同在05年国际遗传算法年会上公开提出用遗传算法来解决多重PCR引物设计这一观点[9]。
以上两项工作基本类似，因此共同存在三个较大的问题：
首先、虽然作者从计算机专家的角度出发强调了遗传算法的应用，但对生物方面的约束条件未进行详细的调查研究而考虑不周全，甚至存在一些生物概念的误解。
例如:
退火温度TM是设计引物时关键的约束条件之一，而计算TM时作者使用了60年代提出的理论估算公式：
TM=（A+T）+4（G+C），该公式忽略了影响退火温度的许多重要因素，已被实践证明无法准确估算退火温度，甚至无法适应单个PCR引物的设计[2]；又例如：
作者将引物自互补和发夹结构两个概念误解为同一个，前者是指完全相同的两条引物互补从而产生二聚物,而后者是指一条引物折叠后和自身产生互补从而形成发夹二级结构；其次、作者对多重PCR生物参数部分未做较好的建模从而导致了生物定义上的错误；例如：
作者因未定义合理的引物评价制度，从而特异性条件无法满足，所以作者自行提出了二次特异性的概念来使条件得到满足，但这是不合理的；又例如：
作者只区分了引物的合不合适，却未分出合适的引物之间的好坏。
最后，作者虽然提供了三个范围的输入给用户，但未给用户提供反馈接口，实验者只能被动接受计算结果，而无法将信息反馈到机器中来设计合适的产物；此外、作者未对遗传算法各参数和各部分未做较好的选择与调整。
以上三个问题，前两个问题导致了程序不能成功地设计可以在实际实验中使用的多重PCR引物；最后一个问题导致了用户缺乏自主权，无法按照实际情况主导程序来选择合适的引物。
综上所述，以上两项工作将遗传算法与引物设计问题结合从而给出了解决问题的新方法，但由于忽略了问题的实用方面从而导致问题未能实际解决。
1.3研究目的与成果
工作目的归纳为：
开发适用于中、小规模生物实验室的多重PCR引物设计软件；并将研究目的推广为：
归纳总结开发多重PCR引物设计软件的一般性思路，提出启发性的方法，并将以上两项付诸实践检验；最后将研究思路概括为：
针对性、实用性和一般性。
为了开发能够在实际生物实验中使用的应用软件，经调查研究后借鉴了应用遗传算法设计多重PCR引物的思想，对多重PCR引物设计软件开发问题的其它部分进行了独立的调查研究与设计。
调查总结了多重PCR引物设计中的生物参数，参考了部分成熟的单重PCR引物设计软件，研究了多目标优化算法；最终建立了解决问题的数学模型，选择了遗传算法作为载体并优化了其各部分的参数和选择计算策略。
开发出在个人电脑（DELLINSPRION630m;CPU:
PentiumIV,1.89GHz;RAM:
512MB）上可运行的多重PCR引物设计软件GAPrimer并撰写了相关用户操作手册；为了使软件在以后能在一个更广泛的领域被应用，以及受到可利用文献资源的限制，开发过程为全英文环境。
GAPrimer使用了23个参数来约束14个生物学条件，通过14个权值来控制上述14个条件，并DOS下操作系统给用户提供了简易的人机接口使得以上37个参数均对用户开放；另外，GAPrimer可以用于对40个引物以下的引物集合评估，还可以用于格式化从网上下载的序列。
1.4本章小结
本章的目的是为下文建立背景框架：
首先回顾了相关的背景知识，随后根据研究目的和工作情况总结了相关文献，最后阐述了研究目的成果；下文中将以GAPrimer的建立过程为基础，对多重PCR引物设计的软件开发问题进行系统的阐述，并对在本章中提出的问题做出详细的解答。
下文将按研究思路发展的顺序分七个章节进行阐述:
问题分析、相关生物学参数简介及其数学模型、引物评价体系、引物设计的遗传算法实现、GAPrimer简介、实验结果及数据分析、总结与展望。
第二章问题分析
2.1多重PCR引物设计问题的数学模型
多重PCR引物设计问题是建立在PCR引物设计的基础上的多目标优化问题。
一方面需要让每对引物满足单重PCR引物的约束条件，另一方面要考虑引物之间的相互影响；与此同时，多重PCR引物设计问题需要被当作一个整体考虑，不能被分割为多个PCR引物设计问题的简单组合[2]。
由此可见，多重PCR引物设计是个全局优化问题，同时也是多目标优化问题。
目前，多重PCR引物设计问题在算法上已被证明为完全非确定多项式（NP-COMPLETE）问题,这意味着解决该问题的高效算法是基本上不存在的[10]。
对于NP-COMPLETE问题，找到全局最优解的唯一方法是穷举所有可能的情况，由于多重PCR引物设计考虑因素较多、数据量较大，因此依靠穷举来实现是多重PCR引物设计基本不可能的。
综上所述，可以认为多重PCR引物设计问题的可行数学模型有两个特点：
一、该模型应适用于多目标优化算法；二、该模型可以较完整地表达NP-COMPLETE问题；同时从实用角度出发，该模型应允许使用者来定义什么是合适的引物。
2.2载体算法的选择
选择载体算法主要需要考虑可行性。
根据对文献的总结，GAPrimer选择了遗传算法作为解决问题的载体，关于遗传算法的简介和其他相关多重PCR引物设计算法的介绍请参见第一章。
选择的原因如下：
首先、遗传算法已经比较成熟，并已被成功运用到许多NP-COMPLETE问题上，这和多重PCR引物设计问题性质相符合。
其次、遗传算法可利用的方法和参考资料较丰富，这为问题解决提供了更多的途径；再次、遗传算法长于解决全局优化问题，而且多目标优化问题是其主要应用对象之一，这使得该算法能够从全局上解决多重PCR引物设计问题。
最后、第一章中提及的两篇相关文献的结果[8][9]证实了遗传算法在多重PCR引物设计问题上的可行性。
2.3本章小结
本章对问题性质进行了分析，并阐述了选择遗传算法的原因，一方面为了启发读者选择并应用新的算法到多重PCR引物设计问题上，从而使问题得到更好的解决。
另一方面、本章中仅对问题进行了分析，为下文提供相关背景，关于建立数学模型的具体内容将在三、四两章中进行阐述。
遗传算法的应用将在第五章中进行讨论。
第三章相关生物学参数简介及其数学模型
设计多重PCR引物，首先要解决引物评估问题，即回答什么是合适的引物，以及这些引物的相互干扰是否在可接受的范围内。
因此多重PCR引物设计问题中要考虑的生物参数分为两部分：
首先是确定单个引物的约束参数，据此可以区分引物的优劣；其次是确定引物间的约束参数，据此来评价引物之间的干扰是否在可接受的范围内。
下文就以上两部分参数分别进行了讨论，并给出了相应的数学模型。
同时，本章也中给出了一些关键生物参数在GAPrimer中的具体评估算法，和对算法的改进意见。
3.1单个引物的约束参数
根据对文献的总结，GAPrimer对单个引物考虑了如下10个约束条件：
引物长度、引物退火温度、引物的GC含量、二级发夹结构、寡核苷酸长度、自互补、3’端自互补、GC夹[11]、特异性、3’端稳定性[12]。
（1）引物长度（PrimerLength）：
引物长度指引物上的碱基数目。
模型中使用引物长度上界，引物长度下界和引物长度理想值三个参数来约束；默认范围为[18，27]，理想值为20。
（2）引物退火温度（MeltingTemperature,TM）：
引物退火温度指引物与模板序列结合的温度。
模型中使用退火温度上界，退火温度下界和退火温度理想值三个参数来约束。
默认范围为[57,63],理想值为60，单位为摄氏度；此外模型中使用了单价阳离子浓度和配对引物浓度两个参数来对退火温度进行盐浓度校正，默认值分别为50mM、50nM[12]。
退火温度的计算使用了NearestNeighbor（NN）模型以及SantaLucia的NN模型参数及其相应盐浓度校正公式[13]。
也有其他NN模型参数可以用来计算退火温度，但根据文献记载，以上NN参数是较新、预测较为准确的模型，因此本程序选择此模型为默认方法。
程序中采用了最基础的扫描算法，但也有较高效的算法可以完成这一任务,例如分式编程[14]。
此外，TM的定义是溶液中50%引物与模板发生结合的温度，虽然计算TM是目前最普遍的观点，并且在目前的单重、多重PCR中应用得较好，但事实上更为准确的约束应该是计算引物的熵值（deltaG）[2]。
因此，为了提高多重PCR引物设计的效果，可以将此项约束转化为引物熵值的计算。
（3）引物GC含量（GCContent）：
引物GC含量指引物中含有鸟嘌呤（G）和胞嘧啶（C）的百分比。
模型中使用GC含量上界，GC含量下界，GC含量理想值三个参数来约束。
默认范围为[40，60]，默认理想值为0，0表示不考虑默认理想值，通过设定一个非零值可以激活次参数。
（4）自互补（SelfComplement）：
自互补指引物自身和自身发生互补，结果是形成非特异性的二聚物。
模型中使用最大互补值这个参数来约束，默认值为7。
GAPrimer借鉴了Primer3的评分系统来预测形成二聚物的可能性，对于每对互补的碱基对评分为+1，对于每对错配的碱基对，分数为-1，总分不可为负，总分为0表示形成二聚物的可能性很小[12]。
图3.1中的总分为2。
在这种评分体制下，7分以下（含7分）的引物基本可以成功引导实际多重PCR实验[15]。
然而更精确的评分方法是用NN模型估测熵值[16]，但此方法会很大程度加重程序的计算复杂度，现有引物设计软件基本不用。
图3.1引物互补评分体系举例
图3.2引物互补扫描方法
GAPrimer中的算法，按照常规扫描了所有的可能组合，图3.1中引物长度为20，有39种可能组合。
引物在程序中按图3.2中的方式移动，取所有可能组合中的最大值作为引物的该项的分数。
（5）3’端自互补（3’EndSelfComplement）：
3’端自互补指引物在其3’端与自身发生互补配对，由于延伸是在3’端开始的，所以应该避免引物间3’端的互补。
模型中使用最大末端互补值这个参数来约束。
默认值为3。
评分体系和扫描算法同第四项（自互补），但3’端自互补仅考虑在3’端开始的不间断的互补。
（6）发夹二级结构（HairpinSecondaryStructure）：
发夹结构是经引物自身互补折叠后形成，其包括两部分：
主干（stem）和环（loop），参见图3.3。
模型中使用最大发夹结构互补值这个参数来约束，默认值为3。
发夹二级结构的评分体系与自互补相同，图3.2中引物的总分为3。
据文献记载，最稳定的发夹结构为环中含有4个或5个碱基的结构，主干中至少含有两对碱基[15]；据此，扫描算法扫描所有可能的环中含4个或5个碱基的结构，例如图3.3中的引物，长度为23，需要比较176次。
最后取所有可能组合中的最大值作为该项分数。
图3.3引物发夹二级结构
（7）寡核苷酸长度（PolyX）：
寡核苷酸长度指引物中相同碱基连续出现的次数,例如寡核苷酸片段TTTTT，其长度为5；在设计引物是应避免此现象出现。
模型中使用了最大寡核苷酸长度这个参数来约束，默认值为5。
（8）GC夹（GCClamp）：
GC夹指在引物3’端左后几个碱基中含有几个G或C，GC夹可以增加3’端稳定性。
模型中使用最少GC夹个数这个参数来约束。
默认值0，0代表默认不考虑这一项。
GAPrimer中考虑了3’端最后的3个碱基，因此最少GC夹个数最高可设为3。
（9）特异性（Specificity）：
特异性指引物在模板上有唯一的结合点，而不会在错误的结合点上和模板发生结合。
模型中使用最大误配值这个参数来约束，默认值为9。
评分体系和自互补中的体系相同；扫描算法考虑所有可能的组合情况：
引物相对于模板的移动与图3.3引物自互补扫描算法中的类似；取所有可能组合的最大值作为该项的评分。
（10）3’端稳定性（3’EndStability）：
3’端稳定性指引物3’端最后5个碱基的熵值，熵值越小代表越稳定。
模型中使用最大三端稳定性这个参数来约束。
3’端稳定的组合为‘GCGCG’,使用NN模型以及SantaLucia的NN模型参数计算，可得其熵值为-6.86kcal/mol[12]。
因此取-6.9kcal/mol为最大三端稳定性默认值，该值不允许用户改变。
3.2引物间的约束参数
根据对文献总结，GAPrimer对多个引物之间，以及引物对的考虑了如下6个约束条件：
引物长度差、退火温度差、引物间互补、引物间3’端互补、产物长度、产物间长度差。
（1）引物长度差（PrimerLengthDifference）：
指引物之间的长度差，引物设计需要是所有引物的长度尽量接近。
模型中使用引物最大长度差这个参数来约束，默认值为3。
算法中需评估所有可能的引物对，最后取所有长度差之和的平均值作为整个引物集的长度差。
n重PCR需要2n个引物，因此需评估2n2-n种组合。
（2）退火温度差（TMDifference）：
指引物间退火温度之差，由于所有的引物在同一试管中发生反应，因此退火温度要尽可能相似。
模型中使用最大退火温度差这个参数来约束，默认值2，单位为摄氏度。
算法中需评估所有可能的引物对，最后取所有退火温度差之和的平均值作为整个引物集的退火温度差。
n重PCR需要2n个引物，因此需评估2n2-n种组合。
（3）引物间互补（Complement）：
引物间互补指不同引物之间的互补，结果是形成二聚物。
模型中使用最大互补值这个参数来约束，默认值为7。
评估所用评分体系和扫描算法均与评估自互补中所用的相同。
算法需要评估所有可能的引物对，最后取所有互补分数之和加上自互补分数的平均值作为整个引物集的互补分数。
n重PCR需要2n个引物，因此需评估2n2-n种组合[17][18]，此外需加上对2n个引物自互补的组合的评估。
（4）引物间3’端互补（EndComplement）：
引物间3’端互补指不同引物之间的3’端互补，结果是形成二聚物。
模型中使用最大末端互补值这个参数来约束，默认值为3。
评估所用评分体系和扫描算法均与3’端自互补中所用的相同。
算法需要评估所有可能的引物对，最后取所有3’端互补分数之和加上3’端自互补分数的平均值作为整个引物集的3’端自互补分数。
n重PCR需要2n个引物，因此需评估2n2-n种组合，此外需加上对2n个引物3’端自互补的组合的评估。
（5）产物长度（ProductLengthDifference）：
产物长度指最终产生的产物的长度，一方面，产物序列中应该包含目标序列；另一方面，产物越长成本越高，因此需限制产物长度。
程序中使用最大产物与目标序列长度差这个参数来约束，默认值为150，代表产物的长度最大超过目标序列长度150个碱基。
（6）产物间长度差（ProductLength）：
多重PCR的产物在同一试管中，为了在凝胶电泳中区分不同的产物，需要使产物间有一定的长度差。
模型中使用产物长度差上界和产物长度差下界这两个参数来约束，默认范围是[50,100]。
算法需要评估所有可能的产物对，最终取长度差之和的平均值作为整个引物集的产物间长度差。
n重PCR需要需评估（n2-n）/2种组合
3.3本章小结
本章中介绍了14个独立的生物学约束条件，及其相应的23个约束参数；同时讨论了相关的评估模型与算法；并提出了一些改进意见。
这14个约束条件为引物集综合评价体系的建立提供了基础，对引物集的综合评价体系的建立将在第四章中进行讨论。
第四章引物集评价体系
4.1评价体系的数学模型
首先，由于多重PCR引物设计是多目标优化问题，不能只简单地评价某个引物集是否合适,也基本不可能完全满足每一个条件；一个比较合理的可行方法是评估引物集合所有的属性，然后提供一个综合评价体系框架，然后把决定权交给使用者。
其次，遗传算法需要一个宽松的限制条件，这有助于算法搜索到全局的最优解[19]；因此该评价体系应该提供一个约束较小的环境，从而为遗传算法的应用提供基础。
GAPrimer中所采用的数学模型首先评估上述14个独立生物学约束条件，然后通过14个权值来平衡这14个约束条件，并依此建立起引物集评价体系。
模型首先对每一个约束条件归一化，给出位于区间[0,1]的评分Si，参见公式4.1；分数越低表示该项约束条件满足得越好，0代表该项条件完全被满足。
为了提供宽松的约束环境，程序并不抛弃某项分数超过1的引物，除非这些引物其它所有约束项综合起来有助于提高整个引物集合的综合质量，否则这些分数超过1的引物会在程序中自动被淘汰。
此外，是否保留这样的引物集合最终是由用户通过权值设定来决定的。
在各项分数都给出后，针对每一项模型中给出位于区间[0,100]的权值Wi，Wi为浮点数，所有14个权值的和为100，参见公式4.1。
以上模型代表对于多重PCR引物设计问题的14个约束条件用户给以的考虑程度在整个问题中所占的百分比。
例如，Wi=20，代表用户在整个问题用户对i条件给以20%的考虑。
从另一个角度理解，实际上用户通过对模型中权值的设定定义了一个平衡点；因此这种模型允许用户来决定什么才是选择的标准。
另外，程序允许用户设定的所有权值的和不等于100，程序会将这些权值归一化后放大100倍，从而保证在任何权值系统下所有权值的和均位于区间[0,100]；无论用户采用何种权值设定，用户永远是在设定整个问题的一个平衡点。
最后，模型综合了权值和各项分数，给出了位于区间[0,100]的综合分数Score，参见公式4.2。
由于对14个约束条件施加了宽松的评价体系，因此某个引物集的综合分数可以超过100，这代表这个引物集不在用户设定得范围内，但这样的引物集一定会在遗传算法中被淘汰，除非它是最优解。
综合分数为0的引物集是理想引物集，这样的引物集不存在任何干扰，当然这是基本上不可达到的；综合分数小于100的引物集都是满足条件的引物集，综合分数越低引物集的综合质量越好；例如，引物集的综合分数为40，代表了这个引物集已经达到了理想引物集60%的效果。
根据GAPrimer运行的结果统计，一般分数在50以下的引物集都具有较好的质量，从在理论上而言完全能够成功地引导多重PCR；当然分数低于100的引物集都是满足条件的，但程序的任务是根据用户设定的平衡点找出解决问题的全局最优解。
（4.1）
（4.2）
综上所述，GAPrimer所使用的数学模型中满足条件的引物集的分数小于100，分数越低引物集质量越好；模型也帮助用户定义了一个解决多目标优化问题的框架，并通过14个权值使得用户可以设定整个问题的平衡点。
在下一届中将对14个权值进行详细的描述，指出其含义及其约束的生物学条件，并给出默认值。
4.2权值及其约束对象
（1）约束引物长度的权值（WeightforPrimerLength）：
默认值为10，提高此权值可以使引物长度更接近用户设定的引物长度理想值。
（2）约束引物长度差的权值（WeightforPrimerDeltaLength）：
默认值为10，提高此权值可以使引物集长度差缩小。
（3）约束退火温度的权值（WeightforTm）：
默认值为10，提高此权值可以使退火温度更接近用户设定的退火温度理想值。
（4）约束退火温度差的权值（WeightforDeltaTm）：
默认值为10，提高此权值可以使引物集退火温度差缩小。
（5）约束GC含量的权值（WeightforGCContent）：
默认值为5。
若用户已设定GC含量理想值，则提高此权值可以使引物GC含量更接近理想值；若用户未设定GC含量理想值，则提高此权值意味着对引物GC含量给以更多的考虑，从而保证引物集中
升级会员 