科学网

 找回密码
  注册

tag 标签: 时间信号

相关帖子

版块 作者 回复/查看 最后发表

没有相关内容

相关日志

应用BETS检测数据集的时间信号
raindyok 2020-8-5 19:25
序言: 准确推断快速进化病原(如 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, M het )和等时模型(Isochronous model, M iso ),采用广义垫脚石采样法(Generalized stepping-stone sampling, GSS)获得两个模型的边际似然估值(marginal likelihood estimation, MLE),最后应用贝叶斯因子法(Bayes Factor, BF)确定哪个模型更适合分析的数据集。其中,异时模型( M het )是应用实际采样时间校准分子钟(Tip date)计算替换速率,而等时模型( M iso )不启用采样时间校准功能(即:without tip date)而假定以一个随意的替换速率估值。 是否有 时间信号的判断依据 : M het 优于 M iso 时,说明数据集含有时间信号。 当logBF=Log(P(Y|Mhet))-Log(P(Y|Miso)) 1时,说明 M het 略优于 M iso ; 当logBF=Log(P(Y|Mhet))-Log(P(Y|Miso)) 1时3时, M het 明显优于 M iso ; 当logBF=Log(P(Y|Mhet))-Log(P(Y|Miso)) 1时5时, M het 绝对优于 M iso 。 准备工作 1. nexus格式的比对序列 为方便提取序列的时间信息,建议每条序列名称命名为类似 Accession_date格式,如:AY99893_2001.67211 2.BEAST 1.x系列 需要提前计算好序列数据的最适替代模型和树先验,本实例为简化直接采用GTR+F+I和Constant size的树先验进行配置。 分析流程: 1. M het 模型的设置 Step 1 :打开BEAUti配置XML文件,在Paritions标签下载入Nexus格式的序列,如下图所示 Step 2 :Tips标签下,勾选“Use tip dates”,点击“Parse Dates”并设置规则提取每个序列中年份信息(示例数据年份前有 _符号间隔) Step 3: 在Site标签下,设置数据集对应的替换模型,示例数据为GTR+F+I,这里+F为Empirical 碱基频率,+I为InvariantSites Step 5: 为演示方便,分子钟头模型用Strict模型,Tree Prior用Coalescent的ConstantSize模型,这里不作截图,直接切换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,如下图所示: 设置完成,可以通过界面右下角的“GenerateBEAST File”生成XML,保存为xxxx_tipdate.xml后用 BEAST打开并运行,直至完成。 2. M iso 模型的设置 同上,将标签切换至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法比较 M het 模型与 M iso 模型 打开运行结果文件xxxxx.mle.result.log,找到尾部复制log marginal likelihood后的数值进行比较。 M het 的MLE值为 -4133.385, M iso 的MLE值为 -4411.406,logBF= 278.021,远大于5 。 4.结果解读 根据计算获得的BF值可知, M het 绝对 优于 M iso 表明数据集具有足够的时间信号,可用于后续的Bayesiandating分析。 参考: 1. BETS简明教程: http://beast.community/bets_tutorial 2.BETS软件文章 (Preprint版): https://www.biorxiv.org/content/10.1101/810697v1 推荐阅读 BET 引用文章 (1)Livia et al., 2020, Nature Microbiology, https://doi.org/10.1038/s41564-020-0706-0 (2)Düx et al., 2020, Science, https://doi:10.1126/science.aba9411 示例数据 请访问我的网盘 http://raindy.ys168.com/ 示例数据目录内
个人分类: 软件教程|6579 次阅读|0 个评论
日期随机化检验(DRTs)之 TipDatingBEAST篇
raindyok 2018-9-26 21:30
絮语 : 病毒(或其他病原)的进化速率(rate)和时间尺度(timescale)可以通过对具有时间结构(Time-structured)的序列进行系统发育分析获得,应用病毒序列的采样时间进行分子钟校准前,通常需要进行日期随机化检验(date-randomization test)判断分析的数据集中是否有足够的时间信号(Temporal signal)可用于校准。 1. 两个检验标准 (CR1和CR2):通过原始(正确的)采样时间进行校准估算获得的进化速率平均值(Mean)与应用日期随机化后的数据获得的进化速率95%置信区间(CI)进行比较(如下图中的小红点与黑色线段的上下限)。如果两者之间没有重合,则通过相对宽松的标准1(即:CR1);通过原始(正确的)采样时间估算获得进化速率的95%CI与随机化后的数据获得进化速率的任何区间(如下图的红色虚线段与黑色线段的任意一个值)均没有重合,则能通过更为严格的标准2(即:CR2); 2. 两个不同方法 :标准的 Bayesian dating permutation tests 和相同采样时间归组的 Clustered permutation tests 。 本文以标准的方法为示例,应用R包TipDatingBeast进行DRTs简明图解分析。 所需工具: 1. R包: 官网: https://cran.r-project.org/ ; 下载安装对应操作系统的版本 ; 2.库文件 TipDatingBeast 安装R后,在界面窗口输入以下命令(红色部分)安装 install.packages(TipDatingBeast) 图解操作: Step 1. 配置用于BEAST运行的XML文件 应用 BEAUti 1.8.x配置时,在Tip标签下,启用Tip-date功能,将每条序列的采样时间根据不同规则提取出来用于校准(如下图所示),其他根据规范逐一配置完成,最好导出配置好的xml文件(本例为PMMoV_DRTs.xml),注意该xml文件中的时期是正确采样时间(非随机化的); Step 2. 随机化 XML 文件的采样日期 打开R-console (Windows下推荐R-gui.exe)前,先将工作目录改为PMMoV_DRTs.xml 所在的目录,可以通过“Misc”菜单下的“Change Working Directory”修改工作目录,也可以通过命令 setwd() 来修改当前工作目录(如: setwd(/Users/Raindy/Documents/R/DRTs) )。 随后输入下列命令分别载入TipDatingBeast文件和生成20个日期随机化的xml文件: library(TipDatingBeast) #载入 TipDatingBeast库 RandomDates(name=PMMoV_DRTs,reps=20,writeTrees=F) #将目录内名称为 PMMoV_DRTs 的xml文件进行日期随机化,并生成20个不同的重复 #注意只需要文件名(不需要扩展名xml) 当出现 Replicates done: 20时,20个日期随机化后的xml文件已经在工作目录生成。 Step 3.应用BEAST分别运行前两步得到的*.xml,获得运行日志文件*.log 在大批量运行这些xml文件前,建议先随机抽个随机化的xml文件,看看程序是否会存在如下类似报错信息。 SEVERE: Parsing error - poorly formed XML (possibly not an XML file): The string -- is not permitted within comments. java.lang.RuntimeException: Terminate 解决办法: 用文本编辑器搜索字符串 !-- write tree log to file ,将它后一行的 !-- 的字符串删除后删除,保存xml 文件 Step 4. 应用 TipDatingBeast绘制DTRs结果 将 Step 3运行获得的21个 *.log复制到当前工作目录,输入以下命令,对DRTs结果进行绘制: library(TipDatingBeast) #载入 TipDatingBeast库 PlotDRT(name= PMMoV_DRTs,reps=20,burnin=0.1) # 对 PMMoV_DRTs.log和 PMMoV_DRTs.rep1.log, PMMoV_DRTs.rep2.log,.... PMMoV_DRTs.rep20.log的参数进行绘制,可以输入VAR查看变量 #示例数据最佳分子钟模型为relax clock with lognormal distribution,所以对应的参数VAR为meanRate 绘制完成,工作目录会生成三个文件, meanRate .stats.csv和两个PDF文件,均为DRTs的结果,一个是原数绘制的,另一个是数值取Log10后绘制的结果。 Step 5.结果美化及解读 可以使用R对 meanRate .stats.csv进行绘制,效果图如下图Figure S2B 所示,DRTs结果表明该数据集通过DRTs检验,具有时间信号,其采样时间可以用于分子钟校准。 Raindy注: 示例图详见参考文献5: https://www.sciencedirect.com/science/article/pii/S0168170218303563 参考文献 : 1. Duchêne, S., Duchêne, D., Holmes, E.C., Ho, S.Y.W., 2015. The performance of the date-randomization test in phylogenetic analyses of time-structured virus data. Molecular biology and evolution 32, 1895-1906. 2. Murray, G.G.R., Wang, F., Harrison, E.M., Paterson, G.K., Mather, A.E., Harris, S.R., Holmes, M.A., Rambaut, A., Welch, J.J., 2016. The effect of genetic structure on molecular dating and tests for temporal signal. Methods in Ecology and Evolution 7, 80-89. 3. Rieux, A., Khatchikian, C.E., 2017. TipDatingBeast: an R package to assist the implementation of phylogenetic tip-dating tests using BEAST. Molecular ecology resources 17, 608-613. 4. Tong, K.J., Duchene, D.A., Duchene, S., Geohegan, J.L., Ho, S.Y.W., 2017. A comparison of mnethods for estimating substitution rates from ancient DNA sequence data. bioRxiv. 5. Guan, X., Yang, C., Fu, J., Du, Z., Ho, S.Y.W., Gao, F., 2018. Rapid evolutionary dynamics of pepper mild mottle virus. Virus research 256, 96-99.
个人分类: 软件教程|9276 次阅读|0 个评论

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-30 23:46

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部