本发明涉及生物信息分析技术领域，提供了一种基于混合并行方式的蛋白质热力学分析高效随机模拟方法，包括：步骤A：确定蛋白质能量模型和能量区间；步骤B：确定所述蛋白质能量区间的分段方式；步骤C：模拟及计算蛋白质系统态密度。采用本发明提供的方法，可以高效地分析和研究蛋白质折叠的整个热力学过程，进而对蛋白质折叠过程进行探索和研究。
1.一种基于混合并行方式的蛋白质热力学分析高效随机模拟方法，其特征在于，包括：步骤A：确定蛋白质能量模型和能量区间；
采用ECEPP蛋白质能量模型，ECEPP能量力场的表达形式为：EECEPP＝EC+ELJ+EHB+ETor其中， 是两电荷之间的库伦作用力，rij表示原子i和j之间的距离； 是两原子之间的兰纳-琼斯作用力；
是氢键作用力；ETor＝∑lUl(1±cos(nlξl))是两面角旋转作用力，ξl是第l个两面角；所述步骤A进一步包括：对所使用的蛋白质能量区间进行离散化处理，若取k个能量bin区间值，则对[Emin,Emax]平均划分k个bin区间，用每个bin区间中间的一个能量值代表能量区间值；所述步骤B进一步包括：步骤B1：对能量区间平均分为M段，设相邻子能量区间之间的重合度等于Δ个bin区间，则每一段含有 个bin区间；
步骤B2：依照当前计算得到的蛋白质系统态密度函数的对数S(E)分布特点，自适应地对能量区间分段，若某个子能量区间为[Ebegin,Eend]，则所述步骤C进一步包括：
通过MPI的主从进程模式和OpenMP的多线程并行模式，模拟及计算蛋白质系统态密度；
在所述主从进程模式的N个分进程中，分进程1为主进程，其余分进程均为从进程；
步骤S11：初始化蛋白质系统态密度函数的对数S(E)＝lng(E)＝0，直方图H(E)＝0，Emin≤E≤Emax，修正因子df＝In(f)，定义因子f，则修正因子df定义为f的自然对数，并初始化df＝1；
步骤S13：依照所确定的蛋白质能量区间的分段方式将能量区间Emin≤E≤Emax分成M段，并分配到M个分线程中，t＝1；
步骤S14：在每个分线程中，对原来的构型限制在相应的子能量区间里进行随机变动，产生新的构型，计算能量Enew，根据Metropolis准则确定新构型被接受的概率，t＝t+1；
步骤S15：所有线程间相互通信，综合得到整个区间的S(E)和H(E)，s＝s+1；
步骤S16：所有进程间相互通信，主进程收集所有从进程的Stmp(E)和Htmp(E)并累加计算出全局的S(E)和H(E)，即全局的S(E)＝S(E)+所有从进程的Stmp(E)，全局的H(E)＝H(E)+所有从进程的Htmp(E)，将全局的S(E)和H(E)广播给所有从进程，判断直方图平缓条件：若不满足则返回执行步骤S12继续迭代；若满足则执行步骤S17；
步骤S17：改变修正因子df，再返回执行步骤S12继续迭代，直到满足进程终止条件其中 求得S(E)，得到蛋白质系统相对的态密度g(E)＝eS(E)。
2.如权利要求1所述的方法，其特征在于，在所述步骤S14中，根据Metropolis准则确定新构型被接受的概率进一步包括：若接受新构型，则：
S(Enew)＝S(Enew)+df,H(Enew)＝H(Enew)+1；
S(Eold)＝S(Eold)+df,H(Eold)＝H(Eold)+1。
3.如权利要求1所述的方法，其特征在于，所述从进程包括如下步骤：步骤S21：初始化蛋白质系统态密度函数的对数S(E)＝lng(E)＝0，Stmp(E)＝lngtmp(E)＝0，直方图H(E)＝0，Htmp(E)＝0，Emin≤E≤Emax，修正因子df＝In(f)，定义因子f，则修正因子df定义为f的自然对数，并初始化df＝1；
步骤S23：依照所确定的蛋白质能量区间的分段方式将能量区间Emin≤E≤Emax分成M段，并分配到M个分线程中，t＝1；
步骤S24：在每个分线程中，对原来的构型限制在相应的子能量区间里进行随机变动，产生新的构型，计算能量Enew，根据Metropolis准则确定新的构型被接受的概率，t＝t+1；
步骤S25：所有线程间相互通信，综合得到整个区间的Stmp(E)和Htmp(E)，s＝s+1；
步骤S26：所有进程间相互通信，从进程将Stmp(E)和Htmp(E)发送给主进程，然后接收经主进程计算得出的全局的S(E)和H(E)更新原来的S(E)和H(E)，将Stmp(E)和Htmp(E)初始化为
若不满足则返回执行步骤S22继续迭代；若满足则执行步骤S27；
步骤S27：改变修正因子df，再返回执行步骤S22继续迭代，直到满足进程终止条件其中
4.如权利要求3所述的方法，其特征在于，在所述步骤S24中，根据Metropolis准则确定新的构型被接受的概率进一步包括：若接受新构型，则：
S(Enew)＝S(Enew)+df,H(Enew)＝H(Enew)+1,Stmp(Enew)＝Stmp(Enew)+df,Htmp(Enew)＝Htmp(Enew)+1；
S(Eold)＝S(Eold)+df,H(Eold)＝H(Eold)+1,Stmp(Eold)＝Stmp(Eold)+df,Htmp(Eold)＝Htmp(Eold)+1。
5.如权利要求1所述的方法，其特征在于，在所述步骤S17中，改变修正因子df的方式为：先连续进行N次迭代的，f＝fα，0＜α＜1，再进行1次迭代的 反复重复上述方式。
6.如权利要求3所述的方法，其特征在于，在所述步骤S27中，改变修正因子df的方式为：先连续进行N次迭代的，f＝fα，0＜α＜1，再进行1次迭代的 反复重复上述方式。
基于混合并行方式的蛋白质热力学分析高效随机模拟方法
【技术领域】
[0001] 本发明涉及生物信息分析技术领域，特别是涉及一种基于混合并行方式的蛋白质热力学分析高效随机模拟方法。【背景技术】
[0002] 蛋白质折叠主要研究蛋白质如何在短时间内从一维多肽链折叠为天然三维结构，形成具有生命功能的大分子。生物体的遗传信息（DNA）通过RNA转录和翻译过程传递给蛋白质（即中心法则），因此蛋白质折叠也被称为第二遗传密码，它的研究可以帮助揭示生命遗传信息的表达和功能传递的奥秘。在从一维多肽链到天然三维结构的折叠过程中，蛋白质可发生误折叠或聚集，其结构和功能因此受到破坏，从而引起‘折叠病’，比如老年痴呆症等。因此蛋白质折叠研究对探索多种‘折叠病’机理意义重大。
[0003] 目前，研究蛋白质折叠的算法大多数都在分子动力学模拟和随机模拟中实现。一般而言，分子动力学模拟常用于研究蛋白质系统的动力学过程；而随机模拟则可以研究蛋白质系统的整个热动力学过程。针对使用高精确度的全原子蛋白质模型的模拟，需要计算成千上万个原子之间的多种相互作用力，对于分子动力学模拟只能模拟纳秒级的蛋白质折叠过程，故其在微秒到毫秒时间内的蛋白质折叠研究中具有很大的局限性；此外，分子动力学模拟也受一个初始实验构型的影响。而随机模拟不但能用于微秒到毫秒时间内的蛋白质折叠研究，而且不依赖于一个具体的初始构型，可以更广泛地搜索构型空间。
[0004] 经典的WangLandau算法就是随机模拟领域最吸引人最有发展情景的新算法，它能解决生物信息学、统计物理学等多个领域的很多复杂问题。比如在蛋白质折叠研究中，该算法有两个最显著的优点：第一，蛋白质模拟不会局限在局部最小能量状态，因而能较好地在整个能量区间进行自由行走；第二，通过该算法可模拟和计算出蛋白质系统态密度，因而就能进一步求解得到宽广温度范围内的很多热动力学量如比热等，这样就能高效地分析和研究蛋白质折叠的整个热力学过程。但WangLandau算法在计算精度和速度上还有待于进一步提升。
[0005] 鉴于此，克服该现有技术所存在的缺陷是本技术领域亟待解决的问题。【发明内容】
[0006] 本发明要解决的技术问题是提供一种基于混合并行方式的蛋白质热力学分析高效随机模拟方法。
[0007] 本发明采用如下技术方案：
[0008] 一种基于混合并行方式的蛋白质热力学分析高效随机模拟方法，包括：
[0009] 步骤A：确定蛋白质能量模型和能量区间；
[0010] 步骤B：确定所述蛋白质能量区间的分段方式；
[0011] 步骤C：模拟及计算蛋白质系统态密度。
[0012] 进一步地，所述步骤A进一步包括：
[0013] 采用ECEPP蛋白质能量模型，ECEPP能量力场的表达形式为：
[0014] EECEPP=EC+ELJ+EHB+ETor
[0015] 其中， 是两电荷之间的库伦作用力，rij表示原子i和j之间的距离； 是两原子之间的兰纳-琼斯作用力；
是氢键作用力；ETor=∑lUl(1±cos(nlξl))是两面角旋转作用力，ξl是第l个两面角。
[0016] 进一步地，所述步骤A进一步包括：
[0017] 对所使用的蛋白质能量区间进行离散化处理，若取k个能量bin区间值，则对[Emin,Emax]平均划分k个bin区间，用每个bin区间中间的一个能量值代表能量区间值。
[0018] 进一步地，所述步骤B进一步包括：
[0019] 步骤B1：对能量区间平均分为M段，设相邻子能量区间之间的重合度等于Δ个bin区间，则每一段含有 个bin区间；
[0020] 步骤B2：依照当前计算得到的蛋白质系统态密度函数的对数S(E)分布特点，自适应地对能量区间分段，若某个子能量区间为[Ebegin,Eend]，则
[0021] 进一步地，所述步骤C进一步包括：
[0022] 通过MPI的主从进程模式和OpenMP的多线程并行模式，模拟及计算蛋白质系统态密度。
[0023] 进一步地，在所述主从进程模式的N个分进程中，分进程1为主进程，其余分进程均为子进程。
[0024] 进一步地，所述主进程包括如下步骤：
[0025] 步骤S11：初始化蛋白质系统态密度函数的对数S(E)=lng(E)=0，直方图H(E)=0（Emin≤E≤Emax），修正因子df=1(=lnf=lne)；
[0026] 步骤S12：s=1；
[0027] 步骤S13：依照所确定的蛋白质能量区间的分段方式将能量区间（Emin≤E≤Emax）分成M段，并分配到M个分线程中，t=1；
[0028] 步骤S14：在每个分线程中，对原来的构型限制在相应的子能量区间里进行随机变动，产生新的构型，计算能量Enew，根据Metropolis准则确定新构型被接受的概率，t=t+1；
[0029] 所述步骤S14循环tmax次；
[0030] 步骤S15：所有线程间相互通信，综合得到整个区间的S(E)和H(E)，s=s+1；
[0031] 所述步骤S14和S15循环smax次；
[0032] 步骤S16：所有进程间相互通信，主进程收集所有从进程的Stmp(E)和Htmp(E)并累加计算出全局的S(E)和H(E)，即全局的S(E)=S(E)+所有从进程的Stmp(E)，全局的H(E)=H(E)+所有从进程的Htmp(E)，将全局的S(E)和H(E)的广播给所有从进程，判断直方图平缓条件：
[0033]
[0034] 若不满足则返回执行步骤S12继续迭代；若满足则执行步骤S17；
[0035] 步骤S17：改变修正因子df，再返回执行步骤S12继续迭代，直到满足进程终止条件其中 求得S(E)，得到蛋白质系统相对的态密度g(E)=eS(E)。
[0036] 进一步地，在所述步骤S14中，根据Metropolis准则确定新构型被接受的概率进一步包括：
[0037]
[0038] 若接受新构型，则：
[0039] S(Enew)=S(Enew)+df,H(Enew)=H(Enew)+1；
[0040] 否则：
[0041] S(Eold)=S(Eold)+df,H(Eold)=H(Eold)+1。
[0042] 进一步地，所述从进程包括如下步骤：
[0043] 步骤S21：初始化蛋白质系统态密度函数的对数S(E)=lng(E)=0，Stmp(E)=lngtmp(E)=0，直方图H(E)=0，Htmp(E)=0（Emin≤E≤Emax），修正因子df=1(=lnf=lne)；
[0044] 步骤S22：s=1；
[0045] 步骤S23：依照所确定的蛋白质能量区间的分段方式将能量区间（Emin≤E≤Emax）分成M段，并分配到M个分线程中，t=1；
[0046] 步骤S24：在每个分线程中，对原来的构型限制在相应的子能量区间里进行随机变动，产生新的构型，计算能量Enew，根据Metropolis准则确定新构型被接受的概率，t=t+1；
[0047] 所述步骤S24循环tmax次；
[0048] 步骤S25：所有线程间相互通信，综合得到整个区间的Stmp(E)和Htmp(E)，s=s+1；
[0049] 所述步骤S24和S25循环smax次；
[0050] 步骤S26：所有进程间相互通信，从进程将Stmp(E)和Htmp(E)发送给主进程，然后接收经主进程计算得出的全局的S(E)和H(E)更新原来的S(E)和H(E)，将Stmp(E)和Htmp(E)初始化为0，判断直方图平缓条件：
[0051]
[0052] 若不满足则返回执行步骤S22继续迭代；若满足则执行步骤S27；
[0053] 步骤S27：改变修正因子df，再返回执行步骤S22继续迭代，直到满足进程终止条件其中
[0054] 进一步地，在所述步骤S24中，根据Metropolis准则确定新构型被接受的概率进一步包括：
[0055]
[0056] 若接受新构型，则：
[0057] S(Enew)=S(Enew)+df,H(Enew)=H(Enew)+1,
[0058] Stmp(Enew)=Stmp(Enew)+df,Htmp(Enew)=Htmp(Enew)+1；
[0059] 否则：
[0060] S(Eold)=S(Eold)+df,H(Eold)=H(Eold)+1,
[0061] Stmp(Eold)=Stmp(Eold)+df,Htmp(Eold)=Htmp(Eold)+1。
[0062] 进一步地，在所述步骤S17和S27中，改变修正因子f的方式为：
[0063] 先连续进行N次迭代的f=fα(0<α<1)，再进行1次迭代的
[0064] 反复重复上述方式。
[0065] 与现有技术相比，本发明的有益效果在于：
[0066] 与经典的WangLandau算法相比，本发明使用基于退火机制的更新修正因子方式可提高计算精度和速度，利用一种灵活的能量区间分段方式可使分线程间负载均衡，采用MPI+OpenMP混合编程模型的混合并行方式可大大加快模拟和计算速度。采用本发明提供的方法，可以高效地分析和研究蛋白质折叠的整个热力学过程，进而对蛋白质折叠过程进行探索和研究。【附图说明】
[0067] 图1是本发明实施例基于混合并行方式的蛋白质热力学分析高效随机模拟方法流程图；
[0068] 图2是本发明实施例中蛋白质能量区间的分段方式示意图；
[0069] 图3是模拟及计算蛋白质系统态密度的混合并行方法流程图。【具体实施方式】
[0070] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0071] 此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0072] 本发明提供了一种基于混合并行方式的蛋白质热力学分析高效随机模拟方法，如图1所示，该方法包括：
[0073] 步骤A：确定蛋白质能量模型和能量区间；
[0074] 步骤B：确定蛋白质能量区间的分段方式；
[0075] 步骤C：模拟及计算蛋白质系统态密度。
[0076] 在步骤A中，可采用ECEPP蛋白质能量模型，ECEPP能量力场的表达形式为：
[0077] EECEPP=EC+ELJ+EHB+ETor
[0078] 其中， 是两电荷之间的库伦作用力，rij表示原子i和j之间的距离； 是两原子之间的兰纳-琼斯作用力；
是氢键作用力；ETor=∑lUl(1±cos(nlξl))是两面角旋转作用力，ξl是第l个两面角。由于ECEPP能量力场采用的是角度坐标系，所以其计算效率高于基于笛卡尔坐标系的其他蛋白质模型。
[0079] 为了便于计算机模拟仿真，还可对所使用的蛋白质能量区间进行离散化处理，若取k个能量bin区间值，则对[Emin,Emax]平均划分k个bin区间，用每个bin区间中间的一个能量值代表能量区间值。
[0080] 在步骤B中，为了使分线程间负载均衡，可采用如下具有自适应特点的一种灵活的能量区间分段方式：
[0081] 步骤B1：对能量区间平均分为M段，设相邻子能量区间之间的重合度等于Δ个bin区间，则每一段含有 个bin区间；
[0082] 其中，设相邻子能量区间之间的重合度等于Δ个bin区间，是为了能综合得到整个区间的S(E)和H(E)（即综合所有子能量区间的S(E)和H(E)），需要相邻子能量区间有合适的重合度。
[0083] 步骤B2：依照当前计算得到的蛋白质系统态密度函数的对数S(E)分布特点，自适应地对能量区间分段，若某个子能量区间为[Ebegin,Eend]，则
[0084] 其中，一般地，蛋白质系统态密度函数的对数S(E)=lng(E)是单调递增的上凹函数，自适应地对能量区间分段需使得每一段的 均衡，也即若某个子能量区间为[Ebegin,Eend]，则 如图2所示。
[0085] 在步骤C中，通过MPI的主从进程模式和OpenMP的多线程并行模式，模拟及计算蛋白质系统态密度；在主从进程模式的N个分进程中，分进程1为主进程，其余分进程均为子进程。
[0086] 如图3所示，主进程包括如下步骤：
[0087] 步骤S11：初始化蛋白质系统态密度函数的对数S(E)=lng(E)=0，直方图H(E)=0（Emin≤E≤Emax），修正因子df=1(=lnf=lne)；
[0088] 步骤S12：s=1；
[0089] 步骤S13：依照所确定的蛋白质能量区间的分段方式将能量区间（Emin≤E≤Emax）分成M段，并分配到M个分线程中，t=1；
[0090] 步骤S14：在每个分线程中，对原来的构型限制在相应的子能量区间里进行随机变动，产生新的构型，计算能量Enew，根据Metropolis准则确定新构型被接受的概率：（简称为MCS步）
[0091]
[0092] 若接受新构型，则：
[0093] S(Enew)=S(Enew)+df,H(Enew)=H(Enew)+1；
[0094] 否则：
[0095] S(Eold)=S(Eold)+df,H(Eold)=H(Eold)+1。
[0096] t=t+1，所述步骤S14循环tmax次（也即经过tmax次（如10次）MCS步）；
[0097] 步骤S15：所有线程间相互通信，综合得到整个区间的S(E)和H(E)，s=s+1；
[0098] 所述步骤S14和S15循环smax次（如100次）；
[0099] 步骤S16：所有进程间相互通信，主进程收集所有从进程的Stmp(E)和Htmp(E)并累加计算出全局的S(E)和H(E)，即全局的S(E)=S(E)+所有从进程的Stmp(E)，全局的H(E)=H(E)+所有从进程的Htmp(E)，将全局的S(E)和H(E)的广播给所有从进程，判断直方图平缓条件：
[0100]
[0101] 若不满足则返回执行步骤S12继续迭代；若满足则执行步骤S17；
[0102] 步骤S17：改变修正因子df，再返回执行步骤S12继续迭代，直到满足进程终止条件（也即 如 可取0.0001）；求得S(E)，得到蛋白质系统相对的态密度g(E)=eS(E)。改变修正因子f的方式为：
[0103] 先连续进行N次迭代的f=fα(0<α<1)，再进行1次迭代的
[0104] 反复重复上述方式。
[0105] 从进程包括如下步骤：
[0106] 步骤S21：初始化蛋白质系统态密度函数的对数S(E)=lng(E)=0，Stmp(E)=lngtmp(E)=0，直方图H(E)=0，Htmp(E)=0（Emin≤E≤Emax），修正因子df=1(=lnf=lne)；
[0107] 步骤S22：s=1；
[0108] 步骤S23：依照所确定的蛋白质能量区间的分段方式将能量区间（Emin≤E≤Emax）分成M段，并分配到M个分线程中，t=1；
[0109] 步骤S24：在每个分线程中，对原来的构型限制在相应的子能量区间里进行随机变动，产生新的构型，计算能量Enew，根据Metropolis准则确定新构型被接受的概率：（简称为MCS步）
[0110]
[0111] 若接受新构型，则：
[0112] S(Enew)=S(Enew)+df,H(Enew)=H(Enew)+1,
[0113] Stmp(Enew)=Stmp(Enew)+df,Htmp(Enew)=Htmp(Enew)+1；
[0114] 否则：
[0115] S(Eold)=S(Eold)+df,H(Eold)=H(Eold)+1,
[0116] Stmp(Eold)=Stmp(Eold)+df,Htmp(Eold)=Htmp(Eold)+1。
[0117] t=t+1，所述步骤S24循环tmax次（即经过tmax次（如10次）MCS步）；
[0118] 步骤S25：所有线程间相互通信，综合得到整个区间的Stmp(E)和Htmp(E)，s=s+1；
[0119] 所述步骤S24和S25循环smax次（如100次）；
[0120] 步骤S26：所有进程间相互通信，从进程将Stmp(E)和Htmp(E)发送给主进程，然后接收经主进程计算得出的全局的S(E)和H(E)更新原来的S(E)和H(E)，将Stmp(E)和Htmp(E)初始化为0，判断直方图平缓条件：
[0121]
[0122] 若不满足则返回执行步骤S22继续迭代；若满足则执行步骤S27；
[0123] 步骤S27：改变修正因子df，再返回执行步骤S22继续迭代，直到满足进程终止条件（也即 如 可取0.0001）。改变修正因子f的方式为：
[0124] 先连续进行N次迭代的f=fα(0<α<1)，再进行1次迭代的
[0125] 反复重复上述方式。
[0126] 本发明和经典的WangLandau算法相比，使用基于退火机制的更新修正因子方式可提高计算精度和速度，利用一种灵活的能量区间分段方式可使分线程间负载均衡，采用MPI+OpenMP混合编程模型的混合并行方式可大大加快模拟和计算速度。MPI+OpenMP混合编程模型可以充分利用这两种编程模式的优点，即MPI可以解决多处理器进程间的粗粒度通信，而OpenMP提供轻量级线程可以很好地解决每个多处理器计算机内部各处理器间的交互。
[0127] 本发明能有效模拟和计算蛋白质系统态密度，能进一步求解得到宽广温度范围内的很多热动力学量如比热等，因此能研究蛋白质折叠的整个热力学过程，进而对蛋白质折叠过程进行探索和研究。
[0128] 本领域普通技术人员可以理解实施例的各种方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器（ROM，Read Only Memory）、随机存取存储器（RAM，Random Access Memory）、磁盘或光盘等。
[0129] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
