Time Calibrated Phylogenies

Construction and evaluation of time calibrated trees for each serotype

Serotype A

Check tree model statistics the get_tracer_stats() function calculates similar summary statistics as the Tracer software typically used with BEAST.

Hide code

check_stats <- get_tracer_stats(here("local/beast/a_1/sero_a.log.txt"))

keep_stats <- c("joint", "prior", "likelihood", "treeModel.rootHeight", "age.root.",
                "treeLength", "tmrca.fmd_a_rev.", "clock.rate")

check_stats %>%
  filter(Parameter %in% keep_stats) %>%
  gt() %>%
  tab_header(
    title = md("Serotype A Stats")) %>%
  cols_width(starts_with("Parameter") ~ px(200),
             starts_with("label") ~ px(60),
             everything() ~ px(95)) %>%
  tab_options(table.font.size = "small",
              row_group.font.size = "small",
              stub.font.size = "small",
              column_labels.font.size = "medium",
              heading.title.font.size = "large",
              data_row.padding = px(2),
              heading.title.font.weight = "bold",
              column_labels.font.weight = "bold") %>%
  opt_stylize(style = 6, color = 'gray')

Parameter	Mean	Median	Q_0.025	Q_0.975	ESS
Serotype A Stats
joint	-5660.087	-5658.760	-5681.411	-5645.418	1264
prior	-71.918	-70.707	-91.298	-58.708	1065
likelihood	-5588.169	-5587.791	-5598.043	-5580.312	8107
treeModel.rootHeight	14.227	12.977	9.476	25.213	858
age.root.	1998.715	1999.965	1987.729	2003.467	858
treeLength	61.955	55.221	37.227	120.591	846
clock.rate	0.004	0.004	0.002	0.007	1161

Load MCC tree

Hide code

sero_A.tree <- read.nexus(here("local/beast/a_1/sero_a.mcc.tre"))

plot_time_tree(sero_A.tree, check_stats)

FMDV Effective Population Size
Trees were run using several different clock and prior choices, all showed flat-line Ne. The one here was the simplest, coalescent with contant population size.

Hide code

phylodynamic_process(sero_A.tree, check_stats, x_limits = c("2011-01-01", "2013-01-01"))


 *** inla.core.safe:  rerun to try to solve negative eigenvalue(s) in the Hessian

Serotype Asia1

Hide code

check_stats <- get_tracer_stats(here("local/beast/asia1_1/sero_asia1.log.txt"))

keep_stats <- c("joint", "prior", "likelihood", "treeModel.rootHeight", "age.root.",
                "treeLength", "tmrca.fmd_a_rev.", "clock.rate")

check_stats %>%
  filter(Parameter %in% keep_stats) %>%
  gt() %>%
  tab_header(
    title = md("Serotype Asia1 Stats")) %>%
  cols_width(starts_with("Parameter") ~ px(200),
             starts_with("label") ~ px(60),
             everything() ~ px(95)) %>%
  tab_options(table.font.size = "small",
              row_group.font.size = "small",
              stub.font.size = "small",
              column_labels.font.size = "medium",
              heading.title.font.size = "large",
              data_row.padding = px(2),
              heading.title.font.weight = "bold",
              column_labels.font.weight = "bold") %>%
  opt_stylize(style = 6, color = 'gray')

Parameter	Mean	Median	Q_0.025	Q_0.975	ESS
Serotype Asia1 Stats
joint	-5992.181	-5991.690	-6019.456	-5968.308	1583
prior	-151.880	-151.805	-175.397	-129.265	1063
likelihood	-5840.301	-5840.022	-5856.071	-5826.286	3980
treeModel.rootHeight	10.005	9.950	9.102	11.222	1483
age.root.	2007.157	2007.212	2005.940	2008.060	1483
treeLength	37.323	36.894	30.145	47.059	978
clock.rate	0.006	0.006	0.004	0.007	1163

Load MCC tree

Hide code

sero_Asia1.tree <- read.nexus(here("local/beast/asia1_1/sero_asia1.mcc.tre")) 
plot_time_tree(sero_Asia1.tree, check_stats, legend_pos = c(0.2, 0.5))

FMDV Effective Population Size
Flat Ne, like the other serotypes

Hide code

phylodynamic_process(sero_Asia1.tree, check_stats, x_limits = c("2011-01-01", "2013-01-01"))


 *** inla.core.safe:  rerun to try to solve negative eigenvalue(s) in the Hessian

Serotype O

Hide code

check_stats <- get_tracer_stats(here("local/beast/o_1/sero_o.log.txt"))

keep_stats <- c("joint", "prior", "likelihood", "treeModel.rootHeight", "age.root.",
                "treeLength", "tmrca.fmd_a_rev.", "clock.rate")

check_stats %>%
  filter(Parameter %in% keep_stats) %>%
  gt() %>%
  tab_header(
    title = md("Serotype O Stats")) %>%
  cols_width(starts_with("Parameter") ~ px(200),
             starts_with("label") ~ px(60),
             everything() ~ px(95)) %>%
  tab_options(table.font.size = "small",
              row_group.font.size = "small",
              stub.font.size = "small",
              column_labels.font.size = "medium",
              heading.title.font.size = "large",
              data_row.padding = px(2),
              heading.title.font.weight = "bold",
              column_labels.font.weight = "bold") %>%
  opt_stylize(style = 6, color = 'gray')

Parameter	Mean	Median	Q_0.025	Q_0.975	ESS
Serotype O Stats
joint	-5238.961	-5236.541	-5269.064	-5224.695	279
prior	-80.068	-77.731	-110.322	-67.692	235
likelihood	-5158.893	-5158.525	-5170.130	-5149.610	7552
treeModel.rootHeight	29.087	14.897	11.792	38.048	4372
age.root.	1987.033	2001.224	1978.072	2004.328	4372
treeLength	116.062	52.651	40.024	155.606	4425
clock.rate	0.004	0.004	0.001	0.005	455

Load MCC tree

Hide code

sero_O.tree <- read.nexus(here("local/beast/o_1/sero_o.mcc.tre"))

plot_time_tree(sero_O.tree, check_stats, legend_pos = c(0.2, 0.5))

FMDV Effective Population Size
Flat Ne, like the other serotypes

Hide code

phylodynamic_process(sero_O.tree, check_stats, x_limits = c("2011-01-01", "2013-01-01"))


 *** inla.core.safe:  rerun to try to solve negative eigenvalue(s) in the Hessian