贝叶斯天际线图（Bayesian Skyline Plot，BSP）的方法是基于溯祖理论（Coalescent theory）来量化种群大小与地理历史的关系。采用分子钟或化石校正，借助BEAST系列软件，基于马尔科夫链的蒙特卡洛算法（Monte Carlo Markov chain，MCMC）求出有效种群大小随时间推移的变化。运算的结果通过Tracer v1.6中的Bayesian Skyline重建功能进行读取。
近些年，随着DNA测序成本的降低以及计算机计算效率的提高，多分子标记联合使用的方法在谱系生物地理学研究中越来越受青睐。而BEAST中的EBSP（Extended Bayesian Skyline Plot）正是一种能使用多个基因片段来探究种群大小与时间变化关系的有力工具，尤其针对不同性质的基因以及个体数较少的分析，该方法能更好地估测有效种群大小。简单来讲，BSP是用单个标记进行运算，EBSP是利用多个标记联合分析。
图1 全球气候变化下鸟类种群的快速扩张
本教程利用BEAST软件通过EBSP方法来探究种群大小与时间变化关系，估测有效种群大小。
数据准备
1、用于分析的三个软件
BEAST v2.3.2：用于计算有效种群大小随时间推移的变化；
Tracer v1.7.2：检验各参数的合理性；
R语言：用于画图，自行安装。
图2 软件安装
2、测试数据
来自某种哺乳动物的三个基因：一个核基因、一个线粒体基因和一个来自X染色体的基因。用自己的数据要先把序列格式转化为nex格式，MEGA/DAMBE这两个软件都可以转成，建议用DAMBE转。
图3 本教程测试数据
每个分子标记的每条序列中的分类群名称必须是唯一的，但是不同分子标记间则可以重复。一个小提示：用MEGA或者DAMBE转成nex格式时，要注意最后这个end一定要加上英文的分号（;）。
公众号后台回复“BEAST”，即可获得软件下载地址以及本次教程中使用的测试数据。接下来就是作者才艺展示时间啦！
种群扩张检测流程
1、数据加载
打开BEAUti，点击File menu，选择Import Alignment，将mystery-mammal文件夹中的3个基因文件（mt.nex.nex、nuclear.nex.nex和X.nex.nex）加载到BEAUti中。因为EBSP是一种多位点检测，因此首先需要将所有的位点/基因加载到BEAUti中。
图4 数据加载
2、选择合适的进化模型
选择Site Model面板，为每个基因选择合适的进化模型，我们这里使用HKY替代模型，并且将“Frequencies”更改为“Empirical”。这里的模型大家可以用mrmodeltest或者jmodeltest，进行计算然后进行选择。
图5 选择合适的进化模型
刚才只更改了X.nex.nex这个文件，其他两个mt和nulear也有要更改，如图所示：
图6 mt和nulear文件修改
3、选择分子钟模型
点击Clock model面板，为每个基因选择分子钟模型并且设置分子钟突变速率，做这一步之前，首先点击mode并取消选择自动设置时钟速率（Automatic set clock rate），如下所示：
图7 取消选择自动设置时钟速率
均使用严格分子钟Strict Clock，替代速率的设置一般根据参考文献中设置（一般的话会在10-7-10-9每年，X），mt和nuclear分别设置为0.005、0.05和0.005，即替代速率为0.005每个百万年，此外，对X和nuclear勾选estimate。
图8 选择严格分子钟Strict Clock并设置替代速率
4、设置先验参数
点击Priors面板，设置先验参数。第一步我们需要为每个基因树设置“Coalescent Extended Bayesian Skyline”。下拉“Tree.t:X”，“Tree.t:mt”和“Tree.t:nuclear”右侧的选项，将Yule Model更改为Coalescent Extended Bayesian Skyline，如下所示：
图9 设置基因树Coalescent Extended Bayesian Skyline
然后，单击Tree.t:X左侧的箭头，然后单击出现在Population Model右侧的铅笔图标。如下图所示，这将弹出一个对话框，其中包含一个文本字段，允许您编辑应用于这个特定基因树的种群大小的比例因子。
这一设置是我们解释每个位点的倍性。由于X基因位于X染色体上，我们对X基因设置为1.5的比例因子，因为种群中的每个个体平均有1.5个该基因拷贝。执行相同的步骤将mtDNA和核基因的比例因子分别设置为0.5和2。本教程中我们用0.5表示mtDNA，因为只有雌性mtDNA对有效种群大小有贡献。
图10 设置基因树的种群大小的比例因子
因为我们在clock model模块给X基因和nuclear基因替代速率选择了“estimate”，所以这里我们要给他们设一个上限，避免估计错误。下拉clockRate. c : X和clockRate. c : nuclear左侧的黑色箭头，并将Upper的infinity改为1.0。
图11 参数设置
接下来，我们需要修改种群规模分布均值的超先验。首先，将1/X更改为Normal，然后将mean设置为1，Sigma设置为0.1。
图12 修改种群规模分布均值的超先验
5、设置权重值
我们最后一步要设置权重值。
点击View，然后单击Show Operators panel，就会跳出Operators面板。
图13 Operators面板
然后将以下这些参数更改一下：
图14 参数设置
6、MCMC参数设置
点击MCMC面板，点击tracelog左侧的箭头，将输出文件名修改为mystery.log，如下图所示。
图15 tracelog参数设置
同样，单击EBSPLogger左侧的箭头，确保EBSP输出文件名为“EBSP.log”。对于本教程，我们链长（Chain Length）设置默认的107步，但是需要注意，对于真正的分析来说，太短了，我们为了加快分析速度所以用了此设置。然后点击file，点击save，保存在BEAST文件夹中，命名为test。你就会在BEAST文件夹中发现一个名为test.xml的文件。
图16 输出test.xml文件
7、BEAST软件运算
进入BEAST文件夹，打开BEAST.exe，然后点击Choose File选择test.xml文件，点击run则可以开始运行。
图17 运行BEAST软件
我们链数设置的少，大概5min就可以跑完，下图为运行中和运算后的截图：
图18 BEAST运算结果
会生成两个我们接下来需要的文件：mystery.log和EBSP.log。
图19 BEAST运算结果文件
8、Tracer软件检验后验参数
利用Tracer软件检验后验参数，该软件需要在Java环境下操作。
点击File，选择Import Trace File，选择mystery.log文件，查看各个后验概率参数，一般来说ESS大于200结果最好。
图20Tracer软件检验结果（ESS大于200结果最好）
9、R语言绘图
打开Rstudio，输入命令getwd()，这个命令用于获取当前工作目录，输入之后点击run，即可以执行，你会发现在左下角会跳出“c:/user/adminPC/Documents”，这是我这台电脑的当前工作路径，大家各自找到自己的，并把BEAST软件运行出来的EBSP.log和我附件里面我给你们的脚本plotEBSP.R文件放在当前工作路径。
图21 R绘图准备工作目录
我的工作路径为“c:/user/adminPC/Documents”，所以把这两个文件（EBSP.log和plotEBSP.R）放进去。
图22 R绘图准备工作目录文件
然后回到Rstudio界面，输入命令 source("plotEBSP.R") 点击run，加载脚本，
然后输入命令plotEBSP("EBSP.log", useHPD=FALSE, plotPopFunctions=TRUE, log="y")，点击run，进行绘图。
右下角会出现绘制出来的图片，调整该框大小，点击Export输出图片。
图23 R绘图命令及结果
输出的图片如下所示，时间轴为横坐标，从右往左指的是时间从以前到现在，纵轴为种群数量，从下往上为种群数量增加。从图中可以看出该物种在0.35MYA以前种群趋于稳定，有略增长的趋势，而在0.35MYA以后发生了种群的扩张。
图24 种群扩张分析结果展示
通过扩展的贝叶斯天际线图（EBSP）分析种群扩张的方法就介绍到这里，更多疑问欢迎大家加入凌波微课交流群进行讨论~
网友评论