博文

应用BETS检测数据集的时间信号

已有 8389 次阅读 2020-8-5 19:25 |个人分类:软件教程|系统分类:科研笔记| 时间信号, Temporal signal, 快速进化病毒, Tip date

序言：
准确推断快速进化病原（如 RNA病毒等）的进化速率或最近共祖时间，可以通过对含有时间信号（Temporal signal）或时间结构（Time-structured）的序列进行系统发育分析获得。目前判断数据集中是否有足够的时间信号的方法主要有Root-to-tip （RTT）线性回归分析和日期随机化检验（DRT，详见日志https://user.qzone.qq.com/58001704/blog/1537536180）。然而，这两种方法都有各自的优缺点，比如RTT主要基于严格分子钟假设，而DRT需要比较原始数据和随机化数据，运算量非常大。

本文介绍一个完全基于贝叶斯模型比较的一种检测方法BETS，即 Beysian evaluation of Temporal Signal。该方法主要引入一组成对的模型：异时模型（Heterochronous model, Mhet) 和等时模型（Isochronous model, Miso），采用广义垫脚石采样法（Generalized stepping-stone sampling, GSS）获得两个模型的边际似然估值（marginal likelihood estimation, MLE），最后应用贝叶斯因子法（Bayes Factor, BF）确定哪个模型更适合分析的数据集。其中，异时模型（Mhet) 是应用实际采样时间校准分子钟（Tip date）计算替换速率，而等时模型（Miso）不启用采样时间校准功能（即：without tip date）而假定以一个随意的替换速率估值。

准备工作
1. nexus 格式的比对序列
        为方便提取序列的时间信息，建议每条序列名称命名为类似 "Accession_date" 格式，如：AY99893_2001.67211
2. BEAST 1.x 系列
      需要提前计算好序列数据的最适替代模型和树先验，本实例为简化直接采用GTR+F+I和Constant size的树先验进行配置。

分析流程：
1. Mhet 模型的设置
        Step 1：打开BEAUti配置XML文件，在Paritions 标签下载入Nexus格式的序列，如下图所示

        Step 2：Tips 标签下，勾选“Use tip dates”，点击“Parse Dates”并设置规则提取每个序列中年份信息（示例数据年份前有 _ 符号间隔）

        Step 3: 在Site标签下，设置数据集对应的替换模型，示例数据为GTR+F+I，这里+F 为Empirical 碱基频率，+I 为Invariant Sites

        Step 5: 为演示方便，分子钟头模型用 Strict 模型，Tree Prior 用Coalescent 的Constant Size模型，这里不作截图，直接切换MCMC标签，在MLE下拉菜单选择GSS算法后，点击“Settings”进行设置：

        GSS的step 数、链长还有采样频率可以参考PS/SS一文（https://user.qzone.qq.com/58001704/blog/1529390689）设置，注意选择Tree working prior，与Tree Prior的对应，示例为Constant size，故这里选择Matching coalescent model，如下图所示：

        设置完成，可以通过界面右下角的“Generate BEAST File”生成XML，保存为xxxx_tipdate.xml后用 BEAST打开并运行，直至完成。

2. Miso 模型的设置
        同上，将标签切换至Tips，取消勾选“Use tip dates”，如下所示：

        由于没启用tip功能，需要通过Prior标签设置一下固定的速率值，选择Clock.rate参数，Prior Distribution 选择CMTC Rate Reference输入一个固定的速率值，如2.291E-3 （理论上可以任意值，不过建议参考 Root-to-tip等的结果进行设置）。

        设置完毕，将生成并保存的XML通过BEAST运行，直至结果，在BEAST软件目录会生成对应的文件 xxxxx.mle.result.log。

3. BF法比较 Mhet 模型与 Miso 模型