inst/doc/rates-selection_MrBayes.R

In EvoPhylo: Pre- And Postprocessing of Morphological Data from Relaxed Clock Bayesian Phylogenetics

## ---- setup, echo=FALSE-------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, collapse = TRUE)

## ---- eval = FALSE------------------------------------------------------------
#  library(EvoPhylo)

## ---- include=FALSE-----------------------------------------------------------
devtools::load_all(".")

## ---- eval = FALSE------------------------------------------------------------
#  ## Import summary tree with three clock partitions produced by
#  ## Mr. Bayes (.t or .tre files) from your local directory
#  
#  tree3p <- treeio::read.mrbayes("Tree3p.t")

## -----------------------------------------------------------------------------
data(tree3p)

## -----------------------------------------------------------------------------
## Get table of clock rates with summary stats for each node in 
## the tree for each relaxed clock partition (3 partitions in this tree file)
RateTable_Means_3p <- get_clockrate_table_MrBayes(tree3p, summary = "mean")

## ---- eval = FALSE------------------------------------------------------------
#  ## Export the rate tables (using Mr. Bayes example with 3 partitions)
#  write.csv(RateTable_Means_3p, file = "RateTable_Means_3p.csv")

## ---- fig.width=8, fig.height=8, fig.align = "center", out.width = "70%", message=FALSE,warning=FALSE----
## Plot tree node labels
library(ggtree)
tree_nodes <- ggtree(tree3p, branch.length = "none", size = 0.05) +
  geom_tiplab(size = 2, linesize = 0.01, color = "black", offset = 0.5) +
  geom_label(aes(label = node), size = 2, color="purple")
tree_nodes

## ---- eval = FALSE------------------------------------------------------------
#  ## Save your plot to your working directory as a PDF
#  ggplot2::ggsave("Tree_nodes.pdf", width = 10, height = 10)

## ----eval=FALSE---------------------------------------------------------------
#  ## Import rate table with clade membership (new "clade" column added)
#  ## from your local directory
#  RateTable_Means_3p_Clades <- read.csv("RateTable_Means_3p_Clades.csv", header = TRUE)

## -----------------------------------------------------------------------------
data(RateTable_Means_3p_Clades)

head(RateTable_Means_3p_Clades)

## ---- eval = FALSE------------------------------------------------------------
#  ## Get summary statistics table for each clade by clock
#  clockrate_summary(RateTable_Means_3p_Clades,
#                    file = "Sum_RateTable_Means_3p.csv")

## ---- echo = FALSE------------------------------------------------------------
t1 <- clockrate_summary(RateTable_Means_3p_Clades, digits = 2) 
kableExtra::kbl(t1, caption = "Rate table summary statistics") |> 
  kableExtra::kable_styling(font_size = 15, full_width = FALSE,
                            bootstrap_options = "striped", "condensed")

## ---- fig.width=10, fig.height=5, fig.align = "center", out.width = "100%"----
## Overlapping plots
clockrate_dens_plot(RateTable_Means_3p_Clades, stack = FALSE,
                    nrow = 1, scales = "fixed")

## ---- fig.width=10, fig.height=5, fig.align = "center", out.width = "100%"----
## Stacked plots
clockrate_dens_plot(RateTable_Means_3p_Clades, stack = TRUE,
                    nrow = 1, scales = "fixed")

## ---- fig.width=10, fig.height=5, fig.align = "center", out.width = "100%"----
## Stacked plots with viridis color scale
clockrate_dens_plot(RateTable_Means_3p_Clades, stack = TRUE,
                    nrow = 1, scales = "fixed") +
  ggplot2::scale_color_viridis_d() +
  ggplot2::scale_fill_viridis_d()

## ---- fig.width=8, fig.height=6, fig.align = "center", out.width = "70%"------
## Plot regressions of rates from two clocks
p12 <- clockrate_reg_plot(RateTable_Means_3p_Clades, clock_x = 1, clock_y = 2)
p13 <- clockrate_reg_plot(RateTable_Means_3p_Clades, clock_x = 1, clock_y = 3)
p23 <- clockrate_reg_plot(RateTable_Means_3p_Clades, clock_x = 2, clock_y = 3)

library(patchwork) #for combining plots
p12 + p13 + p23 + plot_layout(ncol = 2)

## ---- eval = FALSE------------------------------------------------------------
#  ## Save your plot to your working directory as a PDF
#  ggplot2::ggsave("Plot_regs.pdf", width = 8, height = 8)

## ----eval = FALSE-------------------------------------------------------------
#  ## Import summary tree with a single clock partitions produced by
#  ## Mr. Bayes (.t or .tre files) from examples directory
#  tree1p <- treeio::read.mrbayes("Tree1p.t")

## -----------------------------------------------------------------------------
data(tree1p)

## -----------------------------------------------------------------------------
RateTable_Means_1p <- get_clockrate_table_MrBayes(tree1p, summary = "mean")

## ---- eval = FALSE------------------------------------------------------------
#  ## Export the rate tables
#  write.csv(RateTable_Means_1p, file = "RateTable_Means1.csv")
#  
#  ## Import rate table after adding clade membership (new "clade" column added)
#  RateTable_Means_1p_Clades <- read.csv("RateTable_Means1_Clades.csv", header = TRUE)

## -----------------------------------------------------------------------------
#Below, we use the rate table with clade membership `RateTable_Means_1p_Clades` that accompanies `EvoPhylo`.
data(RateTable_Means_1p_Clades)

## ----eval = FALSE-------------------------------------------------------------
#  ## Get summary statistics table for each clade by clock
#  clockrate_summary(RateTable_Means_1p_Clades,
#                    file = "Sum_RateTable_Medians1.csv")

## ---- echo = FALSE------------------------------------------------------------
t1 <- clockrate_summary(RateTable_Means_1p_Clades, digits = 2) 
kableExtra::kbl(t1, caption = "Rate table summary statistics") |> 
  kableExtra::kable_styling(font_size = 15, full_width = FALSE,
                            bootstrap_options = "striped", "condensed")

## ---- fig.width=6, fig.height=4, fig.align = "center", out.width = "100%"-----
## Stacked plots with viridis color scale
clockrate_dens_plot(RateTable_Means_1p_Clades, stack = TRUE,
                    nrow = 1, scales = "fixed") +
  ggplot2::scale_color_viridis_d() +
  ggplot2::scale_fill_viridis_d()

## ----eval=FALSE---------------------------------------------------------------
#  ## Import all log (.p) files from all runs and combine them, with burn-in = 25%
#  ## and downsampling to 2.5k trees in each log file
#  posterior3p <- combine_log("LogFiles3p", burnin = 0.25, downsample = 1000)

## ---- results='hide'----------------------------------------------------------
data(posterior3p)

## Show first 10 lines of combined log file
head(posterior3p, 10)

## ---- fig.width=8, fig.height=3, fig.align = "center", out.width = "90%"------
library(ggplot2)
B1 <- plot_back_rates (type = "MrBayes", posterior3p, clock = 1,
                      trans = "log10", size = 10, quantile = 1)
B1
B2 <- plot_back_rates (type = "MrBayes", posterior3p, clock = 2,
                      trans = "log10", size = 10, quantile = 1)
B2
B3 <- plot_back_rates (type = "MrBayes", posterior3p, clock = 3,
                      trans = "log10", size = 10, quantile = 1)
B3

## ---- fig.width=8, fig.height=8, fig.align = "center", out.width = "70%"------
## Plot tree using various thresholds for clock partition 1
A1 <- plot_treerates_sgn(
  type = "MrBayes", trans = "none",                            #Indicates software name output and type of transformation
  tree3p, posterior3p,                                         #Summary tree and posterior files
  clock = 1,                                                   #Show rates for clock partition 1
  summary = "mean",                                            #sets summary stats to get from summary tree nodes
  branch_size = 1.5, tip_size = 3,                             #sets size for tree elements
  xlim = c(-450, -260), nbreaks = 8, geo_size = list(3, 3),    #sets limits and breaks for geoscale
  threshold = c("1 SD", "3 SD"))                               #sets threshold for selection mode
A1

## ---- fig.width=20, fig.height=8, fig.align = "default", out.width = "100%"----
## Plot tree using various thresholds for other clock partition and combine them
A2 <- plot_treerates_sgn(
  type = "MrBayes", trans = "none",                            #Indicates software name output and type of transformation
  tree3p, posterior3p,                                         #Summary tree and posterior files
  clock = 2,                                                   #Show rates for clock partition 1
  summary = "mean",                                            #sets summary stats to get from summary tree nodes
  branch_size = 1.5, tip_size = 3,                             #sets size for tree elements
  xlim = c(-450, -260), nbreaks = 8, geo_size = list(3, 3),    #sets limits and breaks for geoscale
  threshold = c("1 SD", "3 SD"))                               #sets threshold for selection mode

A3 <- plot_treerates_sgn(
  type = "MrBayes", trans = "none",                            #Indicates software name output and type of transformation
  tree3p, posterior3p,                                         #Summary tree and posterior files
  clock = 3,                                                   #Show rates for clock partition 1
  summary = "mean",                                            #sets summary stats to get from summary tree nodes
  branch_size = 1.5, tip_size = 3,                             #sets size for tree elements
  xlim = c(-450, -260), nbreaks = 8, geo_size = list(3, 3),    #sets limits and breaks for geoscale
  threshold = c("1 SD", "3 SD"))                               #sets threshold for selection mode

library(patchwork)
A1 + A2 + A3 + plot_layout(nrow = 1)

## ----eval=FALSE---------------------------------------------------------------
#  ## Import rate table with clade membership
#  RateTable_Means_3p_Clades <- read.csv("RateTable_Means_Clades.csv", header = TRUE)

## -----------------------------------------------------------------------------
data(RateTable_Means_3p_Clades)

## ---- eval = FALSE------------------------------------------------------------
#  ## Get table of pairwise t-tests for difference between the posterior
#  ## mean and the rate for each tree node
#  RateSign_Tests <- get_pwt_rates_MrBayes(RateTable_Means_3p_Clades, posterior3p)
#  
#  ## Show first 10 lines of table
#  head(RateSign_Tests, 10)

## ---- echo = FALSE------------------------------------------------------------
RateSign_Tests <- get_pwt_rates_MrBayes(RateTable_Means_3p_Clades, posterior3p)
t3 <- head(RateSign_Tests, 10)
kableExtra::kbl(t3, caption = "Combined log file") |> 
  kableExtra::kable_styling(font_size = 15, full_width = FALSE,
                            bootstrap_options = "striped", "condensed")

## ---- eval=FALSE--------------------------------------------------------------
#  ## Export the table
#  write.csv(RateSign_Tests, file = "RateSign_Tests.csv")