Biogeography

Extend Bayesian Skyline Plot (EBSP)

Tags: reading; summary

EBSP^[Heled J, Drummond AJ (2008) Bayesian inference of population size history from multiple loci. BMC Evolutionary Biology, 8, 289.] 是一种常用的种群历史动态推算方法，常用于处理 RAD-seq 所得的次世代测序 SNP 数据。RAD-seq 可得到大量的基因位点，常常从中随机选择若干（一般 50 个位点以上）、包含较多 SNP 多态性的位点（SNP 数量需大于 3 ），使用 BEAST 2 ^[Bouckaert R, Heled J, Kuhnert D, Vaughan T, Wu C-H, Xie D, Suchard MA, Rambaut A, and Drummond AJ. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Computation Biology, 10(4), e1003537, Apr 2014. doi: 10.1371/journal.pcbi.1003537. URL http://dx.doi.org/10.1371/journal.pcbi.1003537.] 软件进行 EBSP 的推算。较详细的教程详见 Trucchi & Cristofari 编写的 Extended Bayesian Skyline Plot using highly variable RADseq loci ^[Trucchi E et al. (2014) King penguin demography since the last glaciation inferred from genome- wide data. Proceedings of the Royal Society B: Biological Sciences, 281, 20140528.]。本文就自己运行该计算时遇到的细节问题做一记录。使用的处理平台是 Mac 系统。

RAD-seq 结果的处理

目前有两个常用的软件，ipyrad 和 STACK ，二者在已发表的文献中使用偏好不分伯仲，各有优点。本人的 RAD-seq 结果数据是使用 ipyrad 的前身 pyrad 处理的，pyrad 当时还是在 Python 2 平台上编写的，在下游数据处理上不如最新的 ipyrad 方便，但是实质上二者接近相同。STACK 软件处理结果的文件格式可能略有差异，在前述英文教程中，作者即使用的 STACK 处理 RAD-seq 数据，可能是因为在进行 EBSP 计算前，就处理数量庞大的位点、统计各个位点的 SNP 数来说，STACK 也许更方便（本人未加验证）。因此建议计算 EBSP 的用户首选 STACK 处理 RAD-seq 的结果。

本文中使用的是 pyrad 处理的 RAD-seq 数据，在 EBSP 计算前，处理 RAD-seq 结果的过程相对比较复杂，目前也没有找到方便使用的脚本，有待后续开发吧。因此，手动处理这些庞大的基因位点数据，着实花费不少功夫。

1) RAD-seq 数据质量控制

这一步骤其实在运行 pyrad 或 STACK 时的参数设定中就必须考虑到，具体方法可参考[处理SNP数据：质量控制]，更详细可参见已发表的各类文献。

2) 各位点统计 SNP 数

该步骤的目的，是将 pyrad 输出的 .loci 文件中各个位点的 SNP 数量统计出来。

  • 使用TextWrangler 软件打开 .loci 文件，在菜单 Search 中查找标记每个位点的符号 //，点击 Find all。

  • 拷贝查找的结果到 Excel ，包含位点编号，以及代表该位点 SNP 的符号 -(2个核苷酸基差异) 和 *(2个以上核苷酸差异）。在 Excel 中使用如下图的公式，分别统计每个位点的 -和 * 的数量，相加起来就是每个位点的 SNP 数，分别对应其位点编号（在 Excel 中使用分列功能，以编号左右的 | 符号分列即得）。

  • 根据 SNP 数量将位点进行分类，如 4 个 SNPs 的位点归为 4-SNPs 一类。在 EBSP 计算中，一般使用大于 3 个 SNPs 的位点。

以后应该开发一个针对 .loci 文件自动统计 SNP 数量的脚本。

3) 分割位点并去除 SNPs 过少位点

该步骤的目标，是将 pyrad 输出的 .loci 文件，根据各个位点分割成单独的 .fasta 文件。

使用 python 脚本 loci2fasta.py 将 .loci 文件转换为每个位点上、多个样本个体的 .fasta 文件（有多少个基因座，就生成多少个 .fasta 文件）。在终端中运行：

# 用法：loci2fasta.py [-h] -L file.LOCI

cd /path/to/file_converters-master
loci2fasta.py -L input.loci

# 在目的路径中输出名为 loci 的文件夹，包含所有位点的 .fasta 文件
# 注意，包含各基因座的单独 .fasta 文件的文件夹中， .fasta 文件的数量不能超过8000个，如果基因座数过多，
# 可将其分为若干部分，从每一部分中随机抽取相应比例个数的文件，否则，下面的随机选取过程运行不成功。

按照 SNPs 数量，去掉 SNPs 数量小于 3 的位点。

4) 随机选取若干位点

该步骤的目标，是从所有 SNPs 数量大于 3 的位点中，随机选取若干个文件（本例中为300个）。

cd
cd /path/to/folder/of/seperated/.fasta"
gshuf -zn300 -e *.fasta | xargs -0 -I{} cp -v {} /path/to/folder/for/random/selected/files
# 注意文件夹路径必须完整

这时，随机选取的 300 个 .fasta 文件就复制到了目标新文件夹中，如 loci_300 。

再次对选取的 300 个位点，按照 SNPs 数量进行分类（可根据 SNPs 数重命名 .fasta 文件）。

5) 选择特定样本个体组成数据集

该步骤的目标是，从多个样本中，选择某些你想要分析的样本（如在同一种群或集团的个体），使每个基因位点都包含相同的目标样本个体。

实际中，往往是这种情况，因为你在做 RAD-seq 时，几十甚至上百的样本个体往往不可能属于同一个种群或集团，需要将个体分开处理，特别是计算 Bayesian Skyline Plot 推断种群历史的时候，必须使用同一物种的同一种群来分析。

在 .fasta 文件中个体名称与序列是相邻的行，序列在个体名的下一行，因此要想从众多的、均包含有某个个体的 .fasta 文件中，提取出该个体及其序列，就要提取其个体名和相应的下一行：

# 如要从前面随机选取的 300 个文件中，
# 提取出所有位点中的 t53 个体及其序列：

find path/to/folder/loci_300 -type f | xargs grep -n -A1 "t53" > t53_allloci.txt
# 这时结果中包含来源文件的路径和文件名

find path/to/folder/loci_300 -type f | xargs grep -h -A1 "t53" > t53_allloci.txt
# 这时仅包含想要的个体和序列，不包含来源文件的路径和文件名

# 将这两者一结合，在 Excel 中处理，既能将来源的位点名称
# 与相应的目标个体及序列对应，以供后续处理。

提取出所有目标个体的序列之后，再在 Excel 中，按照位点名称筛选出所有目标个体的序列，这样每个位点就包含了所选择的目标个体的序列，而不是最初的所有个体的序列。

最终的结果，是得到 300 个 .fasta 文件，每个文件代表每个位点，每个位点包含着相同的所选目标个体。此外每个位点的 SNPs 数也要显示在文件名上，以便后续处理。

这些位点的 .fasta 文件，就组成了 BEAST 2 运行 EBSP 的数据集。上述过程相当繁琐耗时，但目前尚未找到简便的脚本，以后可待开发。

BEAST 运行 EBSP

该过程可以详细按照英文教程 Extended Bayesian Skyline Plot using highly variable RADseq loci。这里就不做赘述。要注意的地方有以下几点：

  • 根据 SNPs 数来连锁 site model 和 clock model。tree 全部不连锁。

  • BEAUTi 程序编写 .xml 文件时，在 Mode 菜单中勾选 “Auto set clock rate” 和 “Allow parameter linking”，取消勾选"Automatic set fix mean substitution rate"。

  • 在未知 substitute rate 的情况下，设定 site model 的最后要勾选面板上的 “Fixed mean substitute rate”。

  • 设定 Priors 时，要根据样本的倍型，更改 “Population model”，如二倍体位点时为 2 ，线粒体基因位点为 1.5 ，单倍体位点为 1 。可在最后的 .xml 文件中统一更改。

  • 设定 HKY 模型下的 kappa prior ，将该参数调整到 log-normal 分布下范围为 60-70 。

  • 调出 Operators 面板，更改 Bit Flip， Sample Off Values， Scale， Up Down 四者的值。这四个参数的默认值分别是 30, 15, 15, 5 ，这样的默认值下，运行 EBSP 会非常耗时！必须根据你所选择的位点数，改成相应比例的值，如本例选择了 300 个位点，那么就将这四个值分别改成各自默认值的 600 倍，即 18000, 9000, 9000, 3000 。

  • 最后，删除生成的 .xml 文件中计算各位点的树的命令行，即类似下面的这些行，这样能节省大量磁盘空间：

    <logger id="treelog.t:king_red_snp4_loc40039" fileName="$(tree).trees"
logEvery="1000" mode="tree">
        <log id="TreeWithMetaDataLogger.t:king_red_snp4_loc40039"
spec="beast.evolution.tree.TreeWithMetaDataLogger"
tree="@Tree.t:king_red_snp4_loc40039"/>
  • 运行 BEAST 要达到结果 trace log 文件的各重要参数的 ESS 值大于 100 ，最好大于 200 。否则可适当增加 MCMC 的运行长度，或者重复运行多次。

绘制 EBSP Plot

EBSPAnalyser 处理结果

BEAST 运行所得的另一个 .log 文件即是 EBSP log 文件（默认名称为 EBSP.log），即可用于绘制 EBSP plot 。绘制之前，需用 BEAST 中的 EBSPAnalyser 进行处理，需要使用 beast.jar 调取：

cd BEAST/lib                # 必须 cd 到 beast.jar 所在的文件夹
java -cp beast.jar beast.app.tools.EBSPAnalyser

这时 EBSPAnalyser 软件即被调取，输入文件选择 EBSP.log，再定义输出文件名、重建方式、burnin 比例，即可输出 EBSP plot 的文本文件。

R 绘制 EBSP Plot

如果输出的 EBSPAnalyser.out 文件的标题行后面有空余 tab 空格，需要去掉。

data <- read.table('EBSPAnalyser.out', header = T, sep="\t")
plot(data$time, data$mean, log="y", type="l", col=1, xlim=c(0,0.015),
ylim=c(0.00001,0.05))
polygon(c(data$time, rev(data$time)), c(data$X95HPD.upper,
rev(data$X95HPD.lower)), col = "#98FB98ac", border = NA)

或者使用 ggplot：

require(ggplot2)
read.table('EBSPAnalyser.out', header=T)->data
ggplot(data, aes(x=time)) +
theme_bw() + scale_x_log10() + scale_y_log10() + geom_line(data=ebsp.r, aes(y=median), linetype='dashed') + geom_ribbon(data=ebsp.r, aes(ymin= X95HPD.lower, ymax= X95HPD.upper,
 linetype='dashed', size=.2, alpha=.25)

其他

有一个未解决的问题是，暂时没有找到方法校准 EBSP 中涉及到的分子钟速率，也就是说 RAD-seq 序列的进化速率无法确定，以至于 EBSP plot 的横轴代表的具体历史时间好像无法确定。以往的研究是使用线粒体基因得出的分子钟速率粗略地估算横轴的时间，但是要知道，线粒体基因的数据相比成千个位点的 RAD-seq 数据，信息量太少了！就像用一个精确度不高的尺子来校准精确度高的尺子。因此，如何确定 RAD-seq 序列的分子钟速率，是一个有待解决的问题。

有文献^[Helmstetter, A. J. et al. 2020. The Demographic History of Micro-endemics: Have Rare Species Always Been Rare? bioRxiv: 2020.2003.2010.985853.] 使用 fastsimcoal2 对 EBSP 得出的种群历史进行验证，也就是用 fastsimcoal2 模拟 EBSP 的种群历史场景，验证两者结果的契合度。后续可以参考。

comments powered by Disqus