babette demo

Richèl J.C. Bilderbeek

2022-06-21

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_7fdd2b819204.xml       |
#> |output_state_filename |/home/richel/.cache/beastier/beast2_7fdd6653ba50.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_7fdd2b819204.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: 1655799456878File: beast2_7fdd2b819204.xml seed: 1655799456878 threads: 1Loading package BEAST v2.6.6Loading 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_7fdd6653ba50.xml.state -overwrite /home/richel/.cache/beastier/beast2_7fdd2b819204.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_7fdd6653ba50.xml.state
#> Random number seed: 1655799457827
#> 
#> 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: -288.63130518219134 
#> 
#> Operator                                                    Tuning    #accept    #reject      Pr(m)  Pr(acc|m)
#> ScaleOperator(YuleBirthRateScaler.t:anthus_aco_sub)        0.75000        107         12    0.04000    0.89916 Try setting scaleFactor to about 0.562
#> ScaleOperator(YuleModelTreeScaler.t:anthus_aco_sub)        0.50000         66         65    0.04000    0.50382 Try setting scaleFactor to about 0.25
#> ScaleOperator(YuleModelTreeRootScaler.t:anthus_aco_sub)    0.50000         49         69    0.04000    0.41525 Try setting scaleFactor to about 0.292
#> Uniform(YuleModelUniformOperator.t:anthus_aco_sub)               -        679        557    0.40000    0.54935 
#> SubtreeSlide(YuleModelSubtreeSlide.t:anthus_aco_sub)       1.00000          8        600    0.20000    0.01316 Try decreasing size to about 0.5
#> Exchange(YuleModelNarrow.t:anthus_aco_sub)                       -        149        403    0.20000    0.26993 
#> Exchange(YuleModelWide.t:anthus_aco_sub)                         -          4        114    0.04000    0.03390 
#> WilsonBalding(YuleModelWilsonBalding.t:anthus_aco_sub)           -          5        114    0.04000    0.04202 
#> 
#>      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.367 seconds
#> stderr: 
#> File: beast2_7fdd2b819204.xml seed: 1655799457827 threads: 1
#> Loading package BEAST v2.6.6
#> 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_7fdd165a3f12.log
#> Writing file /home/richel/.cache/beautier/screenlog_7fdd35e08068.csv
#> Writing file /home/richel/.cache/beautier/treelog_7fdd4a17b1c6.trees
#> End likelihood: -128.9033429365592
#> [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 -169.4615
stderr_mean 34.41272
stdev 79.47277
variance 6315.921
median -130.7262
mode n/a
geom_mean n/a
hpd_interval_low -288.6313
hpd_interval_high -127.7622
act 1000
ess 4

Plot the posterior phylogenies

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