马尔科夫链蒙特卡洛算法并行化及其应用

马尔科夫链蒙特卡洛算法并行化及其应用
屈志勇;陈亭;王铁强;孙辰军;周纯葆
【摘要】为使高性能计算助力群体遗传学和系统地理学研究，提出一种基于 MPI (message passing interface)的群体遗传学分析软件，利用集群中多个CPU核心的计算能力加速群体遗传学分析。

进行正确性验证，对并行加速比和并行效率进行评估，在保证计算结果正确性前提下，利用256个 CPU 核心时可以得到最好的并行加速比(185.16)，在利用128个CPU核心时可以得到最好的并行效率(93.68%)。

实验结果表明，利用高性能计算能够进行快速有效的群体遗传学分析。

%To accelerate the research of population genetics and phylogeography based on high performance computing (HPC),a parallel implementation of the IM program using the message passing interface (MPI)on multiple CPU cores was proposed.The validity,parallel speed and parallel efficiency were evaluated.Best parallel speedup (185.16)when using 256 CPU cores and best parallel efficiency (93.68%)when using 128 CPU cores were obtained.The results demonstrate that the divergence popula-tion genetics analysis can be done effectively and rapidly based on HPC.
【期刊名称】《计算机工程与设计》
【年(卷),期】2016(037)007
【总页数】7页(P1811-1816,1826)
【关键词】群体遗传学;系统地理学;IM模型;高性能计算;消息传递接口
【作者】屈志勇;陈亭;王铁强;孙辰军;周纯葆
【作者单位】山西省气象信息中心，山西太原 030006;北京计算机技术及应用研
究所，北京 100854;国网河北省电力公司，河北石家庄 050021;国网河北省电力
公司，河北石家庄 050021;中国科学院计算机网络信息中心超级计算中心，北京100190
【正文语种】中文
【中图分类】TP39
IM(isolation with migration)模型是进行群体遗传学研究的一种重要模型。

在之
前的工作中，首先笔者针对MCMC(Markov chain Monte Carlo)算法所用数据集进行了并行化计算[1]；之后针对MC3(Metropolis-coupled Markov chain Monte Carlo)算法所用马尔科夫链进行了并行化计算[2]。

由于在实际计算中数据集的规模和马尔科夫链的数量均有所限制，因此并行扩展性并不理想。

本文基于以上工作基础，提出一种基于MPI(message passing interface)的具有良好并行扩
展性的MCMC并行算法，能够使IM软件在集群中多个CPU核心环境下进行准确、快速的计算。

本文对正确性、并行加速比和并行效率进行了评估。

从结果可以看出，本文提出的软件可以利用多个CPU核心进行快速、有效的群体遗传学分析计算。

在群体水平研究物种的形成或分化群体遗传学研究是进化研究中一个重要领域[3]。

代表性的做法是首先对待研究的群体或者物种中的个体进行基因取样；之后通过分析基因的遗传变异模式来了解群体分化过程。

通过对多个遗传标记进行连续不断的取样，最终可以得到对群体参数充分好的估计。

规模足够大的，多基因位点的数据集可以为研究群体遗传过程提供更多有用的信息。

IM模型是群体分化研究中的一种重要模型，利用IM模型中的多种群体参数，可
以评估不同因素在群体分化过程中的不同作用。

IM模型是一个能够根据已对位排
列的DNA序列，同时对两个群体的分化时间和群体之间个体迁移概率进行评估的模型[4]。

IM软件是利用IM模型进行群体遗传分析的工具。

IM软件已经被广泛的应用在群体遗传学研究、物种形成研究和物种杂交研究中[5-8]。

IM软件中的群体参数估计是通过MCMC算法来进行计算的。

MCMC算法是一种随机游走算法，算法中允
许从后验概率分布中进行取样。

MCMC算法会建立一条马尔科夫链，该链的平稳概率分布就是MCMC算法最终得到的群体参数的后验概率分布，MCMC算法具
有广泛的应用[9-12]。

MCMC算法有一个主要的缺点就是容易陷入局部最优，因
此在实际中普遍使用的是MC3算法。

MC3算法同时对多条马尔科夫链进行计算，这些马尔科夫链具有不同的概率分布，且只有一条链为目标链，其余的链是为了与目标链进行信息交换，使得目标链有可能跳出局部最优。

MCMC算法需要得到马尔科夫链的平稳概率分布作为最后的后验概率分布。

由于需要长时间的迭代过程来获得马尔科夫链的平稳概率分布，因此MCMC算法的计算非常耗时。

并且MCMC算法的主要计算是马尔科夫链的计算，而MC3算法中同时对多条马尔科夫链进行计算，计算量会随着所用的马尔科夫链的数量成倍增长。

1.1 IM模型
基本IM模型包含的群体参数包括有效群体规模(N1，N2和NA)，群体间的基因
迁移概率(m1和m2)和群体分化时间(t)，如图1(A)所示。

在基本IM模型中分化
后的群体规模随着时间始终保持不变。

这个假设在扩展的IM模型中进行了放松，扩展IM模型中分化后的子群体规模可以随着时间进行变化。

扩展IM模型比基本IM模型多一个群体参数s，其代表祖先群体形成分化后子群体的比例，如图1(B)
所示。

在扩展IM模型中，分化后子群体的初始群体规模分别为sNA和(1-s)NA，其中0<s<1。

基本IM模型和扩展IM模型中的群体参数在IM软件中都是与中性突变概率u相除来表示的。

1.2 MCMC算法
IM软件中后验概率分布是通过式(1)来计算的，其中D代表DNA序列，T代表系统发育树，τ代表系统发育树中的枝长，θ代表IM模型中的群体参数，并且
θ={θ1,θ2,θA,m1,m2,t}
Metropolis-Hastings算法被应用于MCMC算法中，其计算流程如下，其中x代表马尔科夫链的当前状态，x′代表提出的新状态，g(x→x′)代表从给定的状态x得到新状态x′的条件概率。

(1)根据g(x→x′)随机产生新状态x′。

(2)计算新状态的接受概率，记为R
(3)获得(0,1)之间的均匀随机数U，如果U<R，则接收新状态为下一次迭代的初始状态，否则继续以当前状态为下一次迭代的初始状态。

(4)如果到达最大迭代次数退出，否则跳转至(1)。

MCMC算法的计算过程需要充分多的迭代次数才能够得到马尔科夫链的平稳概率分布。

1.3 MC3算法
MC3算法同时对多条马尔科夫链进行计算，不同的马尔科夫链具有不同的概率分布。

每条马尔科夫链具有一个加热参数β，其中0<β≤1。

只有一条链的加热参数为1，称之为冷链；其余多条链的加热参数为0<β<1，称之为热链。

MC3算法的计算流程如下，其中n为马尔科夫链的总数：
(1)对马尔科夫链i进行如下计算i∈(1,2,…,n)。

1)根据)随机产生一个新状态；
2)计算新状态的接受概率，记为R1
3)获得(0,1)之间均匀分布的随机数U1，如果U1<R1，则接受新状态为下一次迭代的初始状态；否则以当前状态为下一次迭代的初始状态。

(2)经过一定迭代次数之后，随机选择两条马尔科夫链j和k进行状态信息的交换，并计算信息交换的接受概率，记为R2
(3)获得(0,1)之间均匀分布的随机数U2，如果U2<R2，则交换马尔科夫链j和k
的状态信息；否则不发生状态信息的交换。

(4)如果到达最大迭代次数退出，否则跳转至(1)。

MC3算法通过冷链与热链之间的状态信息交换能够帮助冷链跳出局部最优。

但代
价是额外增加了多条马尔科夫链的计算量。

针对MCMC算法的并行化加速，在本文之前的工作中利用多个CPU核心同时对
似然估计进行计算来达到加速计算的效果[1]。

然而由于MCMC算法本身的缺点，容易陷入局部最优，因此MC3算法的使用更加广泛。

因此本文之前的工作中针对MC3算法又进行了并行化加速，利用多个CPU核心同时对多条马尔科夫链进行
计算来达到加速计算的效果[2]。

但是由于MC3算法的并行扩展性并不好，只能利用与马尔科夫链的数量相同的CPU核心来进行计算。

鉴于以上两种并行算法的缺点，本文将MCMC算法与MC3算法的并行化相结合，于此同时对进程的虚拟拓
扑结构、进程间通讯，进程间同步和进程间负载平衡等方面进行了优化设计。

本文将进程建立虚拟的拓扑结构，组成N×M的二维网格形式，如图2所示。

其
中DNA序列数据被分割为M个片段，每个进程都对一个数据片段进行相关计算。

马尔科夫链被分为N组，每个进程都对一组马尔科夫链进行相关计算。

如图2中
所示，同一行中的进程具有相同的马尔科夫链，同一列中的进程具有相同的数据片段。

因此同一行中的进程进行的是MCMC算法的并行计算，而同一列中的进程进行的是MC3算法的并行计算。

如果数据被分为M段，马尔科夫链被分为N组，
则总共需要的进程数量是N×M个。

这里有两种特例，当M=1时，是MC3算法的并行计算；当N=1时，是MCMC算法的并行计算。

为了保证计算的正确性，第一，必须保证同一行中的进程除了需要处理的数据片段
之外其它都是相同的；第二，必须保证同一列中的进程除了马尔科夫链之外其它都是相同的；本文利用两个随机数产生数RNG1和RNG2，其中RNG1用于MCMC算法和MC3算法中产生随机数来决定是否接受新状态，不同行之间可以具有不同的RNG1，但是同一行中的RNG1必须相同；RNG2用于产生随机数来确定需要信息交换的马尔科夫链的标号，所有进程必须具有相同的RNG2。

算法的计算流程如下：
(1)初始化阶段：
初始化并行环境
获取命令行中的相关参数
获取进程数量和进程标号
建立进程虚拟拓扑结构、行和列的子通讯器
读入数据片段
设置随机数产生器RNG1和RNG2
(2)MCMC初始化阶段
压缩DNA序列数据
for 进程中包含的马尔科夫链
设置加热参数的初始值
设置系统发育树中的枝长
选择IM模型中的群体参数
建立初始系统发育树
利用包含的数据片段计算局部似然
end for
同一行进程通过局部似然计算全局似然
(3)MCMC计算阶段
for (t=0; t<迭代次数；t++)
for进程中包含的马尔科夫链
产生新状态x′
利用包含的数据片段计算局部似然
通过局部似然计算全局似然
通过式(3)计算新状态接受概率R1
利用RNG1产生随机数U1
if U1≤R1
x(t+1)=x′
else
x(t+1)=x
end if
end for
利用RNG2产生交换马尔科夫链的标号j和k if 进程包含了马尔科夫链j或者k
if 马尔科夫链j和k在同一个进程中
通过式(4)计算交换接受概率R2
利用RNG2产生随机数U2
if U2≤R2
交换马尔科夫链j和k的状态信息
end if
else
根据马尔科夫链状态计算权重
对同一列中的进程交换加热参数和权重
通过式(4)，利用交换得到的加热参数和
权重计算交换接受概率R3
利用RNG2产生随机数U3
if U3≤R3
交换马尔科夫链j和k的加热参数
end if
end if
else
利用RNG2产生随机数U4
end if
end for
收集进程中冷链状态信息并生成结果文件
3.1 正确性验证
本文利用文献[5]中青蟹(Scylla serrata)的线粒体DNA序列数据对IM软件和本文提出的软件的计算结果进行了比较，来验证计算结果的正确性。

数据集1包含的是Red Sea-South China Sea 和West Indian Ocean两个群体的线粒体DNA序列数据，DNA序列长度为535，群体规模分别为10和89。

数据集2包含的是Red Sea-South China Sea 和West Pacific两个群体的线粒体DNA序列数据，DNA序列长度为524，群体规模分别为11和197。

数据集3包含的是3个地区合并组成的群体(western Indian Ocean，Red Sea-South China Sea和West Pacific)与Northwest Australia 群体的线粒体DNA序列数据，DNA序列长度为524，群体规模分别为291和117。

本文针对上述3个数据集分别利用IM软件和本文提出的软件进行了计算，每个数据集分别利用3个不同的随机数种子进行了计算。

利用IM软件(A，B，C)和本文
提出的软件(D，E，F)进行计算得到的群体规模参数s如图3所示。

3个地理亚种
群分别是Red Sea-South China Sea和West Indian Ocean(A,D)，Red Sea-South China Sea和West Pacific(B,E)和combination of three geographic regions(western Indian Ocean，Red Sea-South China Sea和West Pacific)和Northwest Australia(C,F)。

其中IM软件得到的计算结果如图3(a)、3(b)、3(c)
所示；本文提出的软件得到的计算结果如图3(d)、3(e)、3(f)所示。

从结果可以看出，我们得到的计算结果与IM软件几乎相同，计算结果的正确性得到了验证。

3.2 性能评估
本文利用中国科学院超级计算机中心的深腾7000集群对本文提出软件的性能进行了评估。

深腾7000拥有1140个刀片计算结点，每个计算结点包含两个4核心Intel Xeon E5450的CPU，CPU频率为3.00 GHz，每个计算结点拥有32 G内存，计算结点之间利用InfiniBand相连接。

性能评估所用到的数据是从人类线粒体数据库mtDB下载的人类线粒体DNA序列组成的，DNA序列长度为16567，两个群体分别为African和Asian，群体规模
分别为102和92。

IM模型中对DNA序列数据进行分析所利用的DNA突变模型为HKY(Hasegawa-Kishino-Yano)模型。

IM模型中利用的马尔科夫链的数量为
16条，IM软件的初始计算阶段为10000次迭代，正式计算阶段为10000次迭代。

在每一次迭代计算过程中都需要进行马尔科夫链的状态信息交换。

一次迭代过程中系统发育树的更新次数为1次。

性能评估中，IM软件和本文提出的软件分别对数据集利用5个不同的随机数种子进行了计算，利用5次计算的平均时间来对性能进行评估。

IM软件的平均计算时间为167 127 s。

本文从4个CPU核心到256个CPU核心分别进行了计算。

同
样的，本文对并行效率也进行了评估。

并行加速比的定义如下
式中：Sp为并行加速比，Ts为串行软件的计算时间，Tp为并行软件的计算时间。

并行效率的定义如下
式中：Ep为并行软件的并行效率，P为并行软件计算所用到的CPU核心数量。

并行加速比如图4所示，本文提出的软件得到了近似线性的加速比，并且当
N=16，M=16(256个CPU核心)时，得到了最好的加速比(185.16)。

并行效率如图5所示，本文提出的软件的并行效率都在60%以上，并且当N=16，M=8(128个CPU核心)时，得到了最好的并行效率(93.68%)。

影响并行加速比和并行效率
主要有3个因素，包括进程间通讯，进程间同步和进程间的负载平衡。

进程间通讯对集群中多CPU核心环境下的计算来说是无法避免的，包括马尔科夫链之间的状态信息交换和全局似然的计算。

经过特定迭代次数之后，当两条待交换的马尔科夫链在不同的进程中时，图2中这些进程所在行的相同列之间需要进行
马尔科夫链状态信息的交换。

为了减少进程间的通讯量，本文利用交换加热参数来替代交换马尔科夫链的状态信息。

主要依据是马尔科夫链之间唯一的不同之处就是加热参数的不同。

每次迭代计算过程中，每个进程分别利用各自的数据片段进行计算，之后图2中同一行中的所有进程需要利用局部似然来计算全局似然。

这里本
文利用全局规约操作来进行计算。

进程间同步发生在图2中同一行中的进程和同一列中的进程。

对于同一行中的进程，由于包含了相同的马尔科夫链和不同的数据片段，而用于判断是否接受新状态和状态信息是否交换的RNG1是相同的，因此同一行中的所有进程对于每次MCMC计算都会产生相同的计算结果，因此同一行中的所有进程都是同步的。

而同一行中的所有进程是平均分配的数据片段，最多只有一个DNA位点的差距，因此同一行中的所有进程基本都是同步的。

同一列中的所有进程包含了相同的数据片段，不同的马尔科夫链数量。

当马尔科夫链的数量分配不平均时，相应的计算量成倍增加。

尤其对于需要进行状态信息交换的进程，同步需要等待的时间就会加长。

由于每次迭代过程中需要两条马尔科夫链之间进行状态信息的交换，而其余未参与交换的进程就可以进入下一次迭代。

这个过程是随机的，因此各个进程之间的同步也都是动态的，无法人为控制。

为了能够尽量做好进程间同步，本文将马尔科夫链进行了平均分配，并将数据片段进行平均分配。

由于每个进程分得的数据片段都是平均的，最多只相差一个DNA位点，因此同一行中的所有进程负载基本是平衡的。

同一列中的进程之间马尔科夫链的数量严重影响着进程间的负载平衡，本文通过马尔科夫链的平均分配来人为达到负载平衡。

当每个进程只对一条马尔科夫链进行计算时，马尔科夫链的计算能够达到最好的效率。

当数据片段中的序列长度足够长时，计算时间可以掩盖进程间的通讯时间。

本文中的测试数据DNA序列长度为16567，当N=8时，每个进程分得的数据片段长度约为2000，这种情况下可以得到较理想的并行加速比和并行效率。

本文针对目前比较流行的IM模型进行了分析，并对IM软件进行了基于MPI的并行化加速扩展，形成了可在集群中多CPU核心环境下进行群体遗传学分析的软件。

并且针对进程间通讯，进程间同步和进程间的负载平衡进行了优化，以期达到较好的并行加速比和并行效率。

利用真实数据进行了正确性验证，针对相同的数据，本文提出的软件可以得到与IM软件近似相同的结果。

利用人为生成的测试数据进行了并行加速比和并行效率的评估。

本文提出的软件在使用256个CPU核心时，得到了最好的加速比(185.16)、在使用128个CPU核心时，得到了最好的加速效率(93.68%)。

【相关文献】
[1]ZHOU Chunbao,LANG Xianyu,WANG Yangang,et al.A parallel implementation of the isolation with migration model[J].e-Science Technology & Application,2012,3(1):24-29(in Chinese).[周纯葆,郎显宇,王彦棡,等.隔离迁移(isolation with migration)模型数值计算的并行实现
[J].科研信息化技术与应用,2012,3(1):24-29.]
[2]Zhou C,Lang X,Wang Y,et al.Parallel metropolis coupled Markov chain Monte Carlo for isolation with migration model[J].Applied Mathematics & Information
Sciences,2013,7(1L):219-224.
[3]Hey J,Pinho C.Population genetics and objectivity in species
diagnosis[J].Evolution,2012,66:1313-1429.
[4]Wang Y,Hey J.Estimating divergence parameters with small samples from a large number of loci[J].Genetics,2010,184(2):363-379.
[5]He L,Zhang Aibing,Zhu Chaodong,et al.Phylogeography of the mud crab (Scylla serrata) in the Indo-West Pacific reappraised from mitochondrial molecular and oceanographic clues:Transoceanic dispersal and coastal sequential colonization[J].Marine
Ecology,2011,32(1):52-64.
[6]Sergio Marchant,Amy L Moran,Peter B Marko.Out-of-the tropics or trans-tropical dispersal?The origins of the disjunct distribution of the gooseneck barnacle Pollicipes elegans[J].Frontiers in Zoology,2015,12(1):1-16.
[7]Federica Mattucci,Rita Oliveira,Leslie A Lyons,et al.European wildcat populations are subdivided into five main biogeographic groups: Consequences of Pleistocene climate changes or recent anthropogenic fragmentation[J].Ecology & Evolution,2015,6(3):392-395.
[8]Matthew A Campbell,Naoki Takebayashi,J Andrés López.Beringian sub-refugia revealed in blackfish (Dallia):Implications for understanding the effects of Pleistocene glaciations on Beringian taxa and other Arctic aquatic fauna[J].Bmc Evolutionary Biology,2015(15):1-16.
[9]Mathew B,Bauer AM,Koistinen P,et al.Bayesian adaptive Markov chain Monte Carlo estimation of genetic parameters[J].Heredity,2012,109(4):235-245.
[10]Mathieu Rousseau,James Fraser,Maria A Ferraiuolo,et al.Three-dimensional modeling of chromatin structure from interaction frequency data using Markov chain Monte Carlo sampling[J].BMC Bioinformatics,2011,12:414.
[11]Lisa M Sweeney,Ann Parker,Lynne T Haber,et al.Application of Markov chain Monte Carlo analysis to biomathe-matical modeling of respirable dust in US and UK coal miners[J].Regul Toxicol Pharmacol,2013,66(1):47-58.
[12]Paul K Newton,Jeremy Mason,Kelly Bethel,et al.Spreaders and sponges define metastasis in lung cancer:A Markov chain Monte Carlo mathematical model[J].Cancer Res,2013,73(9):2760-2769.。