In richelbilderbeek/babettes: Experimental babette features

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

file.copy(babette::get_babette_path("anthus_aco.fas"), "test_output_0.fas")
file.copy(babette::get_babette_path("anthus_aco.fas"), "my_fasta.fas")
file.copy(babette::get_babette_path("anthus_aco.fas"), "my_alignment.fas")
file.copy(babette::get_babette_path("anthus_aco.fas"), "anthus_aco.fas")
file.copy(babette::get_babette_path("anthus_nd2.fas"), "anthus_nd2.fas")

Examples

For all examples, do load babette:

library(babette)

All these examples assume that BEAST2 is installed at the default location at r get_default_beast2_path().

All examples read the alignment from a FASTA file (usually my_fasta.fas). Instead of a full run, the MCMC chain length is shortened to 10K states, with a measurement every 1K states:

mcmc <- create_mcmc(chain_length = 10000, store_every = 1000)

Example #1: all default

Using all default settings, only specify a DNA alignment.

posterior <- bbt_run(
  "test_output_0.fas",
  mcmc = mcmc
)

All other parameters are set to their defaults, as in BEAUti.

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$test_output_0_trees, width = 2)

Example #2: fixed crown age

Using all default settings, only specify a DNA alignment.

[No screenshot, as this cannot be done in BEAUti yet]

posterior <- bbt_run(
  "my_fasta.fas",
  posterior_crown_age = 15,
  mcmc = mcmc
)

fasta_to_phylo creates a random phylogeny from a FASTA file of a certain crown age.

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_fasta_trees, width = 2)

An alternative is to date the node of the most recent common ancestor of all taxa:

posterior <- bbt_run(
  "my_fasta.fas",
  mcmc = mcmc,
  mrca_priors = create_mrca_prior(
    taxa_names = sample(get_taxa_names("my_fasta.fas"), size = 3),
    alignment_id = get_alignment_id("my_fasta.fas"),
    is_monophyletic = TRUE,
    mrca_distr = create_normal_distr(
      mean = 15.0,
      sigma = 0.025
    )
  )
)

Here we use an MRCA prior with fixed (non-estimated) values of the mean and standard deviation for the common ancestor node's time.

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_fasta_trees, width = 2)

Example #3: JC69 site model

posterior <- bbt_run(
  "my_alignment.fas",
  site_models = create_jc69_site_model(),
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_alignment_trees, width = 2)

Example #4: Relaxed clock log normal

posterior <- bbt_run(
  "my_alignment.fas",
  clock_models = create_rln_clock_model(),
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_alignment_trees, width = 2)

Example #5: Birth-Death tree prior

posterior <- bbt_run(
  "my_alignment.fas",
  tree_priors = create_bd_tree_prior(), 
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = BDBirthRate))
plot_densitree(posterior$my_alignment_trees, width = 2)

Example #6: Yule tree prior with a normally distributed birth rate

posterior <- bbt_run(
  "my_alignment.fas",
  tree_priors = create_yule_tree_prior(
    birth_rate_distr = create_normal_distr(
      mean = 1.0,
      sigma = 0.1
    )
  ),
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_alignment_trees, width = 2)

Thanks to Yacine Ben Chehida for this use case

Example #7: HKY site model with a non-zero proportion of invariants

posterior <- bbt_run(
  "my_alignment.fas",
  site_models = create_hky_site_model(
    gamma_site_model = create_gamma_site_model(prop_invariant = 0.5)
  ),
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_alignment_trees, width = 2)

Thanks to Yacine Ben Chehida for this use case

Example #8: Strict clock with a known clock rate

posterior <- bbt_run(
  "my_alignment.fas",
  clock_models = create_strict_clock_model(
    clock_rate_param = 0.5
  ), 
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate))
plot_densitree(posterior$my_alignment_trees, width = 2)

Thanks to Paul van Els and Yacine Ben Chehida for this use case.

Example #9: Two alignments

posterior <- bbt_run(
  c("anthus_aco.fas", "anthus_nd2.fas"),
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate.aco), color = "red") +
  ggplot2::geom_line(ggplot2::aes(y = birthRate.nd2), color = "green")
plot_densitree(posterior$anthus_aco_trees, width = 2)
plot_densitree(posterior$anthus_nd2_trees, width = 2)

Thanks to Paul van Els for this use case and supplying these FASTA files.

Example #10: Two alignments, different site models

posterior <- bbt_run(
  c("anthus_aco.fas", "anthus_nd2.fas"),
  site_models = list(
    create_hky_site_model(), 
    create_tn93_site_model()
  ),
  mcmc = mcmc
)

ggplot2::ggplot(
  data = posterior$estimates,
  ggplot2::aes(x = Sample)
) + ggplot2::geom_line(ggplot2::aes(y = birthRate.aco), color = "red") +
  ggplot2::geom_line(ggplot2::aes(y = birthRate.nd2), color = "green")
plot_densitree(posterior$anthus_aco_trees, width = 2)
plot_densitree(posterior$anthus_nd2_trees, width = 2)

Since v1.12 this it is supported to have two alignments with different site models, clock models and tree priors.

Thanks to Paul van Els for this use case.

file.remove("test_output_0.fas")
file.remove("my_fasta.fas")
file.remove("my_alignment.fas")
file.remove("anthus_aco.fas")
file.remove("anthus_nd2.fas")

richelbilderbeek/babettes documentation built on May 5, 2019, 7:10 a.m.