babette demo

Richèl J.C. Bilderbeek

2022-08-12

Introduction

This vignette briefly demonstrates multiple features of babette, without going into much detail.

First, load the library:

library(babette)

This vignette shows how to:

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

inference_model <- create_test_inference_model()

Also, in all cases, we use the same BEAST2 options:

beast2_options <- create_beast2_options(verbose = TRUE)

Let babette run ‘BEAST2’

For an alignment, we’ll use a babette example alignment.

if (is_beast2_installed()) {
  out <- bbt_run_from_model(
    fasta_filename = get_babette_path("anthus_aco_sub.fas"),
    inference_model = inference_model,
    beast2_options = beast2_options
  )
  bbt_delete_temp_files(
    inference_model = inference_model,
    beast2_options = beast2_options
  )
}
#> |parameter             |value                                                       |
#> |:---------------------|:-----------------------------------------------------------|
#> |input_filename        |/home/richel/.cache/beastier/beast2_3dfaf36441195.xml       |
#> |output_state_filename |/home/richel/.cache/beastier/beast2_3dfaf211b1690.xml.state |
#> |rng_seed              |NA                                                          |
#> |n_threads             |NA                                                          |
#> |use_beagle            |FALSE                                                       |
#> |overwrite             |TRUE                                                        |
#> |beast2_path           |/home/richel/.local/share/beast/lib/launcher.jar            |
#> |verbose               |TRUE                                                        |
#> Running command: '/usr/lib/jvm/java-11-openjdk-amd64/bin/java -cp /home/richel/.local/share/beast/lib/launcher.jar beast.app.beastapp.BeastLauncher -validate /home/richel/.cache/beastier/beast2_3dfaf36441195.xml'
#>                         BEAST v2.6.6, 2002-2021             Bayesian Evolutionary Analysis Sampling Trees                       Designed and developed by Remco Bouckaert, Alexei J. Drummond, Andrew Rambaut & Marc A. Suchard                                                       Centre for Computational Evolution                         University of Auckland                       r.bouckaert@auckland.ac.nz                        alexei@cs.auckland.ac.nz                                                       Institute of Evolutionary Biology                        University of Edinburgh                           a.rambaut@ed.ac.uk                                                        David Geffen School of Medicine                 University of California, Los Angeles                           msuchard@ucla.edu                                                          Downloads, Help & Resources:                           http://beast2.org/                                      Source code distributed under the GNU Lesser General Public License:                   http://github.com/CompEvol/beast2                                                               BEAST developers:   Alex Alekseyenko, Trevor Bedford, Erik Bloomquist, Joseph Heled,  Sebastian Hoehna, Denise Kuehnert, Philippe Lemey, Wai Lok Sibon Li, Gerton Lunter, Sidney Markowitz, Vladimir Minin, Michael Defoin Platel,           Oliver Pybus, Tim Vaughan, Chieh-Hsi Wu, Walter Xie                                                                   Thanks to:          Roald Forsberg, Beth Shapiro and Korbinian StrimmerRandom number seed: 1660292517403File: beast2_3dfaf36441195.xml seed: 1660292517403 threads: 1Loading package NS v1.1.0Loading package BEAST v2.6.6Loading package BEASTLabs v1.9.7Loading package MODEL_SELECTION v1.5.3Loading package BEAST v2.6.661430_aco: 78 4626029_aco: 78 4630116_aco: 78 4630210_aco: 78 4B25702_aco: 78 4Alignment(anthus_aco_sub)  5 taxa  78 sites  9 patternsFailed to load BEAGLE library: no hmsbeagle-jni in java.library.path: [/usr/lib/jvm/java-11-openjdk-amd64/lib/server, /usr/lib/jvm/java-11-openjdk-amd64/lib, /usr/lib/jvm/java-11-openjdk-amd64/../lib, /usr/lib/R/lib, /usr/lib/x86_64-linux-gnu, /usr/lib/jvm/default-java/lib/server, /usr/lib/R/lib, /usr/lib/x86_64-linux-gnu, /usr/lib/jvm/default-java/lib/server, /usr/lib/R/lib, /usr/lib/x86_64-linux-gnu, /usr/lib/jvm/default-java/lib/server, /usr/java/packages/lib, /usr/lib/x86_64-linux-gnu/jni, /lib/x86_64-linux-gnu, /usr/lib/x86_64-linux-gnu, /usr/lib/jni, /lib, /usr/lib]TreeLikelihood(treeLikelihood.anthus_aco_sub0) uses BeerLikelihoodCore4  Alignment(anthus_aco_sub): [taxa, patterns, sites] = [5, 9, 78]===============================================================================Citations for this model:Bouckaert, Remco, Timothy G. Vaughan, Joëlle Barido-Sottani, Sebastián Duchêne, Mathieu Fourment, Alexandra Gavryushkina, Joseph Heled, Graham Jones, Denise Kühnert, Nicola De Maio, Michael Matschiner, Fábio K. Mendes, Nicola F. Müller, Huw A. Ogilvie, Louis du Plessis, Alex Popinga, Andrew Rambaut, David Rasmussen, Igor Siveroni, Marc A. Suchard, Chieh-Hsi Wu, Dong Xie, Chi Zhang, Tanja Stadler, Alexei J. Drummond   BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis.   PLoS computational biology 15, no. 4 (2019): e1006650.===============================================================================Done!
#> cmd: /usr/lib/jvm/java-11-openjdk-amd64/bin/java -cp /home/richel/.local/share/beast/lib/launcher.jar beast.app.beastapp.BeastLauncher -statefile /home/richel/.cache/beastier/beast2_3dfaf211b1690.xml.state -overwrite /home/richel/.cache/beastier/beast2_3dfaf36441195.xml
#> Creating the beautier temprary folder '/home/richel/.cache/beautier'for BEAST2 tracelog, treelog and screenlog files
#> Creating folder '/home/richel/.cache/beastier'for BEAST2 .xml.state output file
#> error_code: 0
#> stdout: 
#> 
#>                         BEAST v2.6.6, 2002-2021
#>              Bayesian Evolutionary Analysis Sampling Trees
#>                        Designed and developed by
#>  Remco Bouckaert, Alexei J. Drummond, Andrew Rambaut & Marc A. Suchard
#>                                     
#>                    Centre for Computational Evolution
#>                          University of Auckland
#>                        r.bouckaert@auckland.ac.nz
#>                         alexei@cs.auckland.ac.nz
#>                                     
#>                    Institute of Evolutionary Biology
#>                         University of Edinburgh
#>                            a.rambaut@ed.ac.uk
#>                                     
#>                     David Geffen School of Medicine
#>                  University of California, Los Angeles
#>                            msuchard@ucla.edu
#>                                     
#>                       Downloads, Help & Resources:
#>                            http://beast2.org/
#>                                     
#>   Source code distributed under the GNU Lesser General Public License:
#>                    http://github.com/CompEvol/beast2
#>                                     
#>                            BEAST developers:
#>    Alex Alekseyenko, Trevor Bedford, Erik Bloomquist, Joseph Heled, 
#>  Sebastian Hoehna, Denise Kuehnert, Philippe Lemey, Wai Lok Sibon Li, 
#> Gerton Lunter, Sidney Markowitz, Vladimir Minin, Michael Defoin Platel, 
#>           Oliver Pybus, Tim Vaughan, Chieh-Hsi Wu, Walter Xie
#>                                     
#>                                Thanks to:
#>           Roald Forsberg, Beth Shapiro and Korbinian Strimmer
#> 
#> Writing state to file /home/richel/.cache/beastier/beast2_3dfaf211b1690.xml.state
#> Random number seed: 1660292518660
#> 
#> 61430_aco: 78 4
#> 626029_aco: 78 4
#> 630116_aco: 78 4
#> 630210_aco: 78 4
#> B25702_aco: 78 4
#> Alignment(anthus_aco_sub)
#>   5 taxa
#>   78 sites
#>   9 patterns
#> 
#> TreeLikelihood(treeLikelihood.anthus_aco_sub0) uses BeerLikelihoodCore4
#>   Alignment(anthus_aco_sub): [taxa, patterns, sites] = [5, 9, 78]
#> ===============================================================================
#> Citations for this model:
#> 
#> Bouckaert, Remco, Timothy G. Vaughan, Joëlle Barido-Sottani, Sebastián Duchêne, Mathieu Fourment, 
#> Alexandra Gavryushkina, Joseph Heled, Graham Jones, Denise Kühnert, Nicola De Maio, Michael Matschiner, 
#> Fábio K. Mendes, Nicola F. Müller, Huw A. Ogilvie, Louis du Plessis, Alex Popinga, Andrew Rambaut, 
#> David Rasmussen, Igor Siveroni, Marc A. Suchard, Chieh-Hsi Wu, Dong Xie, Chi Zhang, Tanja Stadler, 
#> Alexei J. Drummond 
#>   BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. 
#>   PLoS computational biology 15, no. 4 (2019): e1006650.
#> 
#> ===============================================================================
#> Start likelihood: -412.0608952631166 
#> 
#> Operator                                                    Tuning    #accept    #reject      Pr(m)  Pr(acc|m)
#> ScaleOperator(YuleBirthRateScaler.t:anthus_aco_sub)        0.75000        110          9    0.04000    0.92437 Try setting scaleFactor to about 0.562
#> ScaleOperator(YuleModelTreeScaler.t:anthus_aco_sub)        0.50000         75         56    0.04000    0.57252 Try setting scaleFactor to about 0.25
#> ScaleOperator(YuleModelTreeRootScaler.t:anthus_aco_sub)    0.50000         51         65    0.04000    0.43966 Try setting scaleFactor to about 0.272
#> Uniform(YuleModelUniformOperator.t:anthus_aco_sub)               -        643        563    0.40000    0.53317 
#> SubtreeSlide(YuleModelSubtreeSlide.t:anthus_aco_sub)       1.00000          8        594    0.20000    0.01329 Try decreasing size to about 0.5
#> Exchange(YuleModelNarrow.t:anthus_aco_sub)                       -        220        369    0.20000    0.37351 
#> Exchange(YuleModelWide.t:anthus_aco_sub)                         -         10        124    0.04000    0.07463 
#> WilsonBalding(YuleModelWilsonBalding.t:anthus_aco_sub)           -          8         96    0.04000    0.07692 
#> 
#>      Tuning: The value of the operator's tuning parameter, or '-' if the operator can't be optimized.
#>     #accept: The total number of times a proposal by this operator has been accepted.
#>     #reject: The total number of times a proposal by this operator has been rejected.
#>       Pr(m): The probability this operator is chosen in a step of the MCMC (i.e. the normalized weight).
#>   Pr(acc|m): The acceptance probability (#accept as a fraction of the total proposals for this operator).
#> 
#> 
#> Total calculation time: 0.495 seconds
#> stderr: 
#> File: beast2_3dfaf36441195.xml seed: 1660292518660 threads: 1
#> Loading package NS v1.1.0
#> Loading package BEAST v2.6.6
#> Loading package BEASTLabs v1.9.7
#> Loading package MODEL_SELECTION v1.5.3
#> Loading package BEAST v2.6.6
#> Failed to load BEAGLE library: no hmsbeagle-jni in java.library.path: [/usr/lib/jvm/java-11-openjdk-amd64/lib/server, /usr/lib/jvm/java-11-openjdk-amd64/lib, /usr/lib/jvm/java-11-openjdk-amd64/../lib, /usr/lib/R/lib, /usr/lib/x86_64-linux-gnu, /usr/lib/jvm/default-java/lib/server, /usr/lib/R/lib, /usr/lib/x86_64-linux-gnu, /usr/lib/jvm/default-java/lib/server, /usr/lib/R/lib, /usr/lib/x86_64-linux-gnu, /usr/lib/jvm/default-java/lib/server, /usr/java/packages/lib, /usr/lib/x86_64-linux-gnu/jni, /lib/x86_64-linux-gnu, /usr/lib/x86_64-linux-gnu, /usr/lib/jni, /lib, /usr/lib]
#> WARNING: If nothing seems to be happening on screen this is because none of the loggers give feedback to screen.
#> WARNING: This happens when a filename  is specified for the 'screenlog' logger.
#> WARNING: To get feedback to screen, leave the filename for screenlog blank.
#> WARNING: Otherwise, the screenlog is saved into the specified file.
#> Writing file /home/richel/.cache/beautier/tracelog_3dfaf72f775a7.log
#> Writing file /home/richel/.cache/beautier/screenlog_3dfafd2a8f52.csv
#> Writing file /home/richel/.cache/beautier/treelog_3dfaf6d352e33.trees
#> End likelihood: -126.77270757772963
#> [1] TRUE

Plot the posterior estimates

if (is_beast2_installed()) {
  library(ggplot2)
  p <- ggplot(
    data = out$estimates,
    aes(x = Sample)
  )
  p + geom_line(aes(y = TreeHeight), color = "green")
  p + geom_line(aes(y = YuleModel), color = "red")
  p + geom_line(aes(y = birthRate), color = "blue")
}

Show the effective sample sizes (ESS)

Effective sample sizes, with 20% burn-in removed:

if (is_beast2_installed()) {
  traces <- remove_burn_ins(
    traces = out$estimates,
    burn_in_fraction = 0.2
  )
  esses <- t(
    calc_esses(
      traces,
      sample_interval = inference_model$mcmc$tracelog$log_every
    )
  )
  colnames(esses) <- "ESS"
  knitr::kable(esses)
}
ESS
posterior 4
likelihood 4
prior 4
treeLikelihood 4
TreeHeight 4
YuleModel 4
birthRate 4

For a reliable inference, use an ESS of at least 200.

Show the summary statistics

if (is_beast2_installed()) {
  sum_stats <- t(
    calc_summary_stats(
      traces$posterior,
      sample_interval = inference_model$mcmc$tracelog$log_every
    )
  )
  colnames(sum_stats) <- "Statistic"
  knitr::kable(sum_stats)
}
Statistic
mean -199.683
stderr_mean 61.32501
stdev 141.624
variance 20057.37
median -130.172
mode n/a
geom_mean n/a
hpd_interval_low -412.0609
hpd_interval_high -126.3272
act 1000
ess 4

Plot the posterior phylogenies

if (is_beast2_installed()) {
  plot_densitree(out$anthus_aco_sub_trees, width = 2)
}