babette: Step by Step

Richèl J.C. Bilderbeek

2022-06-21

Introduction

This step-by-step demo shows how to run the babette pipeline in detail.

First, load babette:

library(babette)

In all cases, this is done for a short MCMC chain length of 10K:

inference_model <- create_inference_model()
inference_model$mcmc$chain_length <- 10000
inference_model$mcmc$tracelog$filename <- normalizePath(
  get_beautier_tempfilename(
    pattern = "tracelog_", fileext = ".log"
  ),
  mustWork = FALSE
)
inference_model$mcmc$treelog$filename <- normalizePath(
  get_beautier_tempfilename(
    pattern = "treelog_", fileext = ".trees"  
  ),
  mustWork = FALSE
)

Create a ‘BEAST2’ input file

This step is commonly done using BEAUti. With babette, this can be done as follows:

beast2_input_file <- tempfile(pattern = "beast2_", fileext = ".xml")
create_beast2_input_file_from_model(
  input_filename = get_babette_path("anthus_aco.fas"),
  inference_model = inference_model,
  output_filename = beast2_input_file
)

Display (part of) the ‘BEAST2’ input file

print(head(readLines(beast2_input_file)))
#> [1] "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?><beast beautitemplate='Standard' beautistatus='' namespace=\"beast.core:beast.evolution.alignment:beast.evolution.tree.coalescent:beast.core.util:beast.evolution.nuc:beast.evolution.operators:beast.evolution.sitemodel:beast.evolution.substitutionmodel:beast.evolution.likelihood\" required=\"\" version=\"2.4\">"
#> [2] ""                                                                                                                                                                                                                                                                                                                                                                                   
#> [3] ""                                                                                                                                                                                                                                                                                                                                                                                   
#> [4] "    <data"                                                                                                                                                                                                                                                                                                                                                                          
#> [5] "id=\"anthus_aco\""                                                                                                                                                                                                                                                                                                                                                                  
#> [6] "name=\"alignment\">"

This file can both be loaded by BEAUti and be used by ‘BEAST2’.

The file can be checked if it is indeed a valid input file:

if (is_beast2_installed()) {
  is_beast2_input_file(beast2_input_file)
}
#> [1] TRUE

Run MCMC

This step is commonly done using ‘BEAST2’ from the command-line or using its GUI. With babette, this can be done as follows:

if (is_beast2_installed()) {
  beast2_options <- create_beast2_options(
    input_filename = beast2_input_file
  )
  beastier::check_can_create_file(beast2_options$output_state_filename)
  beastier::check_can_create_treelog_file(beast2_options)
  run_beast2_from_options(
    beast2_options = beast2_options
  )
  testthat::expect_true(file.exists(beast2_options$output_state_filename))
}

Display (part of) the ‘BEAST2’ output files

The .log file contains the model parameters and parameter estimates:

if (is_beast2_installed()) {
  print(head(readLines(inference_model$mcmc$tracelog$filename)))
  print(tail(readLines(inference_model$mcmc$tracelog$filename)))
}
#> [1] "#"                                                                                                                                  
#> [2] "#model:"                                                                                                                            
#> [3] "#"                                                                                                                                  
#> [4] "#<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?><input id=\"posterior\" spec=\"beast.core.util.CompoundDistribution\">"
#> [5] "#    <distribution id=\"prior\" spec=\"beast.core.util.CompoundDistribution\">"                                                     
#> [6] "#        <distribution id=\"YuleModel.t:anthus_aco\" spec=\"beast.evolution.speciation.YuleModel\">"                                
#> [1] "5000\t-1774.9515925560909\t-1862.4118939484665\t87.46030139237557\t-1862.4118939484665\t0.013369563898357617\t87.46030139237557\t240.08126598376356"
#> [2] "6000\t-1761.1391498169985\t-1858.227997698216\t97.08884788121755\t-1858.227997698216\t0.01061639423553826\t97.08884788121755\t291.3859079115836"    
#> [3] "7000\t-1766.4632245951523\t-1861.751041828024\t95.28781723287172\t-1861.751041828024\t0.010073248440069679\t95.28781723287172\t276.74423804585615"  
#> [4] "8000\t-1765.0294949371294\t-1859.3990403191162\t94.36954538198698\t-1859.3990403191162\t0.013601238462251404\t94.36954538198698\t369.81468762036104"
#> [5] "9000\t-1773.0480492221716\t-1867.0580166185284\t94.0099673963568\t-1867.0580166185284\t0.012943639524126681\t94.0099673963568\t203.7827813866175"   
#> [6] "10000\t-1760.5283506732276\t-1856.486521492948\t95.95817081972041\t-1856.486521492948\t0.012709005917065797\t95.95817081972041\t355.1409921542682"

The .trees file contains the alignment, taxa and posterior trees:

if (is_beast2_installed()) {
  print(head(readLines(inference_model$mcmc$treelog$filename)))
  print(tail(readLines(inference_model$mcmc$treelog$filename)))
}
#> Warning in readLines(inference_model$mcmc$treelog$filename): incomplete final
#> line found on '/home/richel/.cache/beautier/treelog_7fdd4c97265.trees'
#> [1] "#NEXUS"                ""                      "Begin taxa;"          
#> [4] "\tDimensions ntax=22;" "\t\tTaxlabels"         "\t\t\t61430_aco"
#> Warning in readLines(inference_model$mcmc$treelog$filename): incomplete final
#> line found on '/home/richel/.cache/beautier/treelog_7fdd4c97265.trees'
#> [1] "tree STATE_6000 = ((((21:0.001383345997232411,6:0.001383345997232411):0.0028981830302774974,(12:0.0018871343596790142,18:0.0018871343596790142):0.002394394667830894):0.002194452356756976,(10:0.003104928866312069,(2:0.0019804972342437636,7:0.0019804972342437636):0.0011244316320683052):0.0033710525179548156):0.004140412851271376,((((22:0.00510796025121596,(4:4.5206157997187827E-4,15:4.5206157997187827E-4):0.004655898671244082):6.238358991964617E-5,(16:0.0010563793717392003,(8:7.522098547140065E-4,13:7.522098547140065E-4):3.041695170251938E-4):0.004113964469396407):7.482356146952035E-4,(((((17:1.1912682526888038E-4,1:1.1912682526888038E-4):9.54678330156633E-5,(20:1.8469646955634786E-4,5:1.8469646955634786E-4):2.9898188728195822E-5):3.756560327299596E-4,3:5.902506910145033E-4):1.5099662934881094E-4,11:7.412473203633142E-4):0.0029935405238445756,19:0.00373478784420789):0.00218379161162292):3.560769360924122E-5,(9:0.005428951653624951,14:0.005428951653624951):5.252354958151005E-4):0.004662207086098209):0.0;"     
#> [2] "tree STATE_7000 = ((19:0.00821632996422691,(((10:0.002615288603575628,(2:0.001939737818190603,7:0.001939737818190603):6.755507853850253E-4):0.003785511883347587,(12:0.00269163458730889,18:0.00269163458730889):0.0037091658996143253):3.296704820905617E-4,(21:0.0014869194159949908,(6:4.536321951204295E-4,22:4.536321951204295E-4):0.0010332872208745612):0.005243551553018786):0.0014858589952131329):0.0018569184758427688,(((((1:9.026145841297384E-4,3:9.026145841297384E-4):8.124291381029175E-4,((5:1.1622720810375009E-4,20:1.1622720810375009E-4):5.833613522902529E-4,17:6.995885603940029E-4):0.001015455161838653):0.0012897890261606357,11:0.0030048327483932916):0.0030378557359765305,(14:0.005972704545081083,(((8:8.142673017185204E-4,13:8.142673017185204E-4):5.791347537157843E-4,16:0.0013934020554343047):0.002774940860389075,(4:4.171073977816861E-4,15:4.171073977816861E-4):0.003751235518041694):0.0018043616292577032):6.998393928873883E-5):4.05427476884798E-4,9:0.00644811596125462):0.0036251324788150584):0.0;"          
#> [3] "tree STATE_8000 = (((14:0.006183078764884041,((4:0.001478422787530654,15:0.001478422787530654):0.003632745472428552,((8:8.013459043763553E-4,13:8.013459043763553E-4):5.60513701742108E-4,16:0.0013618596061184633):0.0037493086538407426):0.0010719105049248352):6.136042001583832E-4,((22:0.004032730418039219,9:0.004032730418039219):0.0018520328672933005,(((17:1.1766179438714649E-4,1:1.1766179438714649E-4):7.753498859570909E-4,(3:8.159545522950728E-5,5:8.159545522950728E-5):8.114162251147302E-4):4.3641413991878607E-4,11:0.0013294258202630235):0.0045553374650694955):9.119196797099053E-4):0.006804555497208979,((((19:6.107677109427863E-4,18:6.107677109427863E-4):0.001213654402103842,12:0.0018244221130466283):0.003792687591296369,(10:0.0019659701677494667,(2:0.0011661545442610446,7:0.0011661545442610446):7.998156234884221E-4):0.0036511395365935308):7.372605997755938E-4,((21:6.671870770783519E-4,20:6.671870770783519E-4):4.429637507157722E-4,6:0.0011101508277941241):0.005244219476324467):0.0072468681581328125):0.0;"   
#> [4] "tree STATE_9000 = (((14:0.006179279968460013,((9:0.004811967944344875,((22:2.6787265271728495E-4,(((5:1.3015819804490392E-4,3:1.3015819804490392E-4):7.99793537428323E-5,17:2.1013755178773621E-4):3.846722830005408E-5,1:2.486047800877903E-4):1.9267872629494655E-5):5.994987530313692E-4,11:8.673714057486541E-4):0.003944596538596221):0.0010549878770816437,(4:0.0017790804346977972,15:0.0017790804346977972):0.0040878753867287215):3.12324147033494E-4):0.0024034917779707406,((8:3.543759115429029E-4,13:3.543759115429029E-4):0.001442520292514876,16:0.0017968962040577788):0.006785875542372975):0.004360867777695928,((((18:0.0019663579981994885,12:0.0019663579981994885):1.5117944505992726E-4,21:0.0021175374432594158):0.005601525144721408,6:0.007719062587980824):1.3568330731526627E-5,((10:0.003321437191065138,(2:0.0019184799794055335,7:0.0019184799794055335):0.0014029572116596046):8.593196856970394E-4,(20:6.638852552048523E-4,19:6.638852552048523E-4):0.003516871621557325):0.0035518740419501734):0.00521100860541433):0.0;" 
#> [5] "tree STATE_10000 = (((((11:4.5990924915845547E-4,(((3:7.248705302978052E-5,1:7.248705302978052E-5):5.863404981634565E-5,17:1.3112110284612617E-4):8.582197253301807E-6,5:1.3970330009942797E-4):3.202059490590275E-4):0.004124168574810248,9:0.004584077823968704):6.257852339361398E-4,(22:7.218508995477771E-4,14:7.218508995477771E-4):0.004488012158357067):0.0012367108606534723,((4:1.3578175433307278E-4,15:1.3578175433307278E-4):0.004617164824919563,((8:0.0014041874521398324,13:0.0014041874521398324):0.0010064013909101902,16:0.0024105888430500226):0.002342357736202614):0.0016936273393056795):0.006262431998507481,((((2:0.0017295746413483815,7:0.0017295746413483815):5.733850865214005E-4,20:0.002302959727869782):0.0010952628458942987,(19:7.699211514630165E-4,10:7.699211514630165E-4):0.0026283014223010643):0.003954499882785811,(((18:8.029814271400781E-4,21:8.029814271400781E-4):0.002554793662073097,12:0.003357775089213175):0.0020938597100662528,6:0.005451634799279428):0.0019010876572704638):0.005356283460515906):0.0;"
#> [6] "End;"

The .xml.state file contains the final state of the MCMC run and the MCMC operator acceptances thus far:

if (is_beast2_installed()) {
  print(head(readLines(beast2_options$output_state_filename)))
  print(tail(readLines(beast2_options$output_state_filename)))
}
#> [1] "<itsabeastystatewerein version='2.0' sample='10000'>"                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
#> [2] "<statenode id='Tree.t:anthus_aco'>(((((10:4.5990924915845547E-4,(((2:7.248705302978052E-5,0:7.248705302978052E-5)32:5.863404981634565E-5,16:1.3112110284612617E-4)40:8.582197253301807E-6,4:1.3970330009942797E-4)39:3.202059490590275E-4)31:0.004124168574810248,8:0.004584077823968704)34:6.257852339361398E-4,(21:7.218508995477771E-4,13:7.218508995477771E-4)28:0.004488012158357067)23:0.0012367108606534723,((3:1.3578175433307278E-4,14:1.3578175433307278E-4)25:0.004617164824919563,((7:0.0014041874521398324,12:0.0014041874521398324)29:0.0010064013909101902,15:0.0024105888430500226)24:0.002342357736202614)38:0.0016936273393056795)30:0.006262431998507481,((((1:0.0017295746413483815,6:0.0017295746413483815)37:5.733850865214005E-4,19:0.002302959727869782)41:0.0010952628458942987,(18:7.699211514630165E-4,9:7.699211514630165E-4)33:0.0026283014223010643)26:0.003954499882785811,(((17:8.029814271400781E-4,20:8.029814271400781E-4)22:0.002554793662073097,11:0.003357775089213175)35:0.0020938597100662528,5:0.005451634799279428)27:0.0019010876572704638)36:0.005356283460515906)42:0.0</statenode>"
#> [3] "<statenode id='birthRate.t:anthus_aco'>birthRate.t:anthus_aco[1 1] (-Infinity,Infinity): 355.1409921542682 </statenode>"                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
#> [4] "</itsabeastystatewerein>"                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        
#> [5] "<!--"                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
#> [6] "{\"operators\":["                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                
#> [1] "{\"id\":\"YuleModelSubtreeSlide.t:anthus_aco\",\"p\":1,\"accept\":12,\"reject\":1908,\"acceptFC\":0,\"rejectFC\":0,\"rejectIv\":885,\"rejectOp\":885},"   
#> [2] "{\"id\":\"YuleModelNarrow.t:anthus_aco\",\"p\":\"NaN\",\"accept\":904,\"reject\":1110,\"acceptFC\":0,\"rejectFC\":0,\"rejectIv\":0,\"rejectOp\":0},"      
#> [3] "{\"id\":\"YuleModelWide.t:anthus_aco\",\"p\":\"NaN\",\"accept\":14,\"reject\":385,\"acceptFC\":0,\"rejectFC\":0,\"rejectIv\":199,\"rejectOp\":199},"      
#> [4] "{\"id\":\"YuleModelWilsonBalding.t:anthus_aco\",\"p\":\"NaN\",\"accept\":40,\"reject\":382,\"acceptFC\":0,\"rejectFC\":0,\"rejectIv\":66,\"rejectOp\":66}"
#> [5] "]}"                                                                                                                                                       
#> [6] "-->"

Parse output

This step is commonly done using Tracer. With babette, this can be done as follows.

Parsing .log file to obtain the parameter estimates:

if (is_beast2_installed()) {
  knitr::kable(head(parse_beast_tracelog_file(inference_model$mcmc$tracelog$filename)))
}
Sample posterior likelihood prior treeLikelihood TreeHeight YuleModel birthRate
0 -7886.241 -7876.842 -9.398615 -7876.842 2.2921346 -9.398615 1.00000
1000 -1957.654 -2022.335 64.680739 -2022.335 0.0692071 64.680739 43.56191
2000 -1765.519 -1856.357 90.837724 -1856.357 0.0140691 90.837724 166.98712
3000 -1765.092 -1856.161 91.069551 -1856.161 0.0133865 91.069551 226.68867
4000 -1771.409 -1860.723 89.313659 -1860.723 0.0123373 89.313659 368.85286
5000 -1774.952 -1862.412 87.460301 -1862.412 0.0133696 87.460301 240.08127

Parsing .trees file to obtain the posterior phylogenies:

if (is_beast2_installed()) {
  plot_densitree(parse_beast_trees(inference_model$mcmc$treelog$filename))
}

Parsing .xml.state file to obtain the MCMC operator acceptances:

if (is_beast2_installed()) {
  knitr::kable(head(parse_beast_state_operators(beast2_options$output_state_filename)))
}
operator p accept reject acceptFC rejectFC rejectIv rejectOp
YuleBirthRateScaler.t 0.75 302 103 0 0 0 0
YuleModelTreeScaler.t 0.50 81 300 0 0 0 0
YuleModelTreeRootScaler.t 0.50 118 295 0 0 59 59
YuleModelUniformOperator.t NaN 2344 1703 0 0 0 0
YuleModelSubtreeSlide.t 1.00 12 1908 0 0 885 885
YuleModelNarrow.t NaN 904 1110 0 0 0 0