【教程】解决Beast计算分化时间迟迟无法收敛的问题

(2018-03-25 16:57:07)

标签：

杂谈

分类：

教程

通过修改计算参数，加快beast计算分化时间(dating)的速度

在使用beast计算复杂的系统树（starbeast算多基因树），并进行dating时，常常碰到各参数迟迟无法收敛的问题。

例：通过7个基因，计算几十个样本之间的系统树，并通过两个化石点进行dating。进行了2亿次迭代，每隔1000代记录一次，得到约20万棵树。treeAnotator得到一致树后，发现各节点的时间远远小于标定的化石点时间(只有正确时间的十分之一左右)。

尝试用tracer查看日志，结果日志文件过大，无法打开（似乎tracer占用的内存不能超过约1024G）。改用R语言查看日志，绘制所有样本最近共祖时间，如下：

http://s15/bmiddle/0022vasxzy7j2WHBWQmde&690

可见在20万棵树中，所有标本最近共祖的时间tmrca(all)是一路上扬的，远没有收敛。查看posterior，如下：

http://s12/bmiddle/0022vasxzy7j2WTy7LJcb&690

发现后验概率（posterior）似乎已经收敛得很好了。查看likelihood值，结果也是类似的。这说明，其实树形已经稳定了，但是各节点时间一直没能达到稳定值。查看.trees文件中不同位置的树，发现确实如此。理论上，等待足够久，各节点的时间也能稳定。但是，从上图看，2亿次迭代，tmrca(all)只上升了3左右。从文献看，tmrca(all)大概在20附近。要等到它稳定，恐怕要几十亿次迭代。这太夸张了！

把Log文件最后几十次记录找出来，截取一些与节点时间相关的参数，如下：

speciesTree.rootHeight Ac-adr24.treeModel.rootHeight

Bc.treeModel.rootHeight Fc2.treeModel.rootHeight

Ic.treeModel.rootHeight Lc.treeModel.rootHeight

Nc.treeModel.rootHeight tmrca(all)

2.87284582

3.270353743 3.317153064 3.279253223 4.005006013 3.759006593

3.378087205 3.63844149

2.550139809

3.270353743 3.339587371 3.440925336 4.003686067 3.663386896

3.546820897 3.627999962

2.550139809

3.270353743 3.339587371 3.293584023 3.870297429 3.663386896

3.546820897 3.61437376

2.366707754

3.270353743 3.58416859 3.293644167 3.874334339 3.663386896

3.380451897 3.619065406

2.366707754

3.270353743 3.58416859 3.276359533 3.881699061 3.543363065

3.498089838 3.530135811

2.366707754

3.270353743 3.584579341 3.276359533 3.302606932 3.543363065

3.498089838 3.546468572

2.366707754

3.270353743 3.584579341 3.276359533 3.323872233 3.300125918

3.490858461 3.72622604

2.806350094

3.270353743 3.579092863 3.282278839 3.321780638 3.300125918

3.569832108 3.373355094

2.538145158

3.270353743 3.706591285 3.282278839 3.317632789 3.300125918

3.325473337 3.373355094

2.538145158

3.270353743 3.568135219 3.292361814 3.675793114 3.300125918

3.440940347 3.373355094

2.377656423

3.270353743 3.67627336 3.28509724 3.474112887 3.30035021

3.393225199 3.8194239

2.377656423

3.270353743 3.277554017 3.325755914 3.276206485 3.32601057

3.303980335 3.35344485

2.17634813

3.270353743 3.275420573 3.36747482 3.661961262 3.35395715

3.663208007 3.359986572

3.036623483

3.270353743 3.272765522 3.364206264 3.674404662 3.35395715

3.663208007 3.36298017

2.496049934

3.270353743 3.272765522 3.364206264 3.64905042 3.35395715

3.470479974 3.384123441

2.268274248

3.270353743 3.526368583 3.351095394 3.648898065 3.345651559

3.521645549 3.444205473

2.380965897

3.270353743 3.549281496 3.350351928 3.443127002 3.345651559

3.521645549 3.499641133

2.526794081

3.270353743 3.545205542 3.321538166 3.443127002 3.348887242

3.521645549 3.499641133

3.14250919

3.270353743 3.280998357 3.540919514 3.285594343 3.332906827

3.521645549 3.499641133

3.206266398

3.270353743 3.441995702 3.417000592 3.285481004 3.471752498

3.525447784 3.435868393

3.206266398

3.270353743 3.440097321 3.391648391 3.627980998 3.471752498

3.294510195 3.293774095

3.206266398

3.270353743 3.440097321 3.389027704 3.684485934 3.292713159

3.429966598 3.293774095

每两行数据之间，隔了1000次迭代。但是，仍然有数据连着几行没有变化。可以理解为，程序搜寻最佳系统树所作的努力，多数花在了调整树形上，少数花在了调整节点时间（枝长/height）上。

相应的，解决的办法就是调整几个平时不被关注的计算用参数，将更多的精力放在调整节点时间上。

在beauti中Operators标签页下，把所有和枝长相关的参数的权重（weight）全部调高。

http://s15/bmiddle/0022vasxzy7j2XefLJAce&690

包括：把每个基因对应的xx.ucld.mean由3调整到30，xx.treeModel.rootHeight由3调整到30，xx.Internal.node.heights由30调整到60，最后把species.nodeReHeight由94调整到180.

再次尝试，结果如下：

http://s10/bmiddle/0022vasxzy7j2Y5I0Wd89&690

第一次记录的树（1000次迭代）后

http://s6/bmiddle/0022vasxzy7j2Y5Sj2t05&690

第2次记录的树（2000）次迭代后

http://s15/bmiddle/0022vasxzy7j2Y5XNue0e&690

第3次记录的树（3000次）迭代后

计算刚刚开始，还远没有出burning

in阶段，节点时间就已经很让人满意了。

附：R语言查看log文件中各参数（与tracer类似）

#install

data.table 安装处理大文件需要的包'data.table'。

install.packages('data.table')

#load data.table载入包

library('data.table')

#set working space设定工作文件夹

setwd(choose.dir())

#choose log file在新出现的窗口出选择beast生成的log文件

choose.files()->logfile

#read log file读取刚刚选择的文件。有进度显示。

fread(input = logfile, sep =

'\t', showProgress = T) -> logdata

#查看log文件一共几行几列。

dim(logdata)

#结果如下：

[1] 199792 360

#即，有199792行，360列

#check all available values查看各参数的名称和序号

names(logdata)

#结果如下：

……………………

……………………

[271] "Bc.treeModel.rootHeight"

"Fc2.treeModel.rootHeight"

"Ic.treeModel.rootHeight"

[274] "Lc.treeModel.rootHeight"

"Nc.treeModel.rootHeight"

"tmrca(all)"

[277] "tmrca(arte)"

"tmrca(outgroup)"

"Ac-adr24.kappa1"

[280] "Ac-adr24.kappa2"

"Ac-adr24.frequencies1"

"Ac-adr24.frequencies2"

……………………

……………………

#绘制某个参数的图。如绘制tmrca(all)：它在上面是第276个参数，则下面第一行中数字为276，其他不变。

logdata[,276]->temp

plot(1:dim(temp)[1],unlist(temp),ylab =

names(temp))

#将特定行的结果保存到工作文件夹下test.csv文件中。例如保存第199593-199792行的数据

write.table(logdata[199593:199792,], file = 'test.csv',

sep = ',')

之后的.csv文件可以用excel查看。

分享：

喜欢

0

赠金笔

阅读┊

收藏

┊

喜欢▼

┊打印┊举报/Report

加载中，请稍候......

前一篇：【教程】win10下安装Biolinux双系统

后一篇：【教程】使用LAMINA测量叶型数据