harrietr: An R package for various phylogenetic and evolutionary biology data manipulations

Why harrietr:

Harriet is believed to be Charles Darwin’s pet
giant turtle. It is thought
that Harriet spent the latter part of her life at in
Brisbane, Australia. Thus, Harriet satisfies the three criteria to name this
package: (1) it is somehow evolutionarily related; (2) it has an Australian
connection; and (3) it avoids the prefix phylo used by many R evolutionarily-
relavant packages. The appended r just helps differentiate to make it easier to
search in Google, and aludes to the fact that it is related to the programming
language R.

How to get it

From CRAN:

1. Add Bioconductor to your list of default repositories:

   setRepositories(ind = 1:2)

2. Install harrietr:

   install.packages("harrietr", dependencies = TRUE)

Latest and gratest version from GitHub:

You must use devtools:

1. If you don’t have devtools installed:

   install.packages('devtools')

2. Add Bioconductor to your list of default repositories:

   setRepositories(ind = 1:2)

3. Install harrietr:

   devtools::install_github("andersgs/harrietr")

How to use it

Three functions are provided at this time:

dist_long — This function takes as input an alignment in DNA.bin format
calculates all the pairwise distances, and returns all the unique pairwise distances
as a data.frame in a long format. For example, in the case of three samples:

You can give dist_long a tree (object of class phylo), and it will add a
fourth column with the pairwise distance obtained from the tree:

melt_dist — This function is used by dist_long, but it takes as input a
distance matrix. This might be useful if you alredy have a distance matrix that
is imported into R.

get_node_support — This function is written to work with trees generated by
IQTREE. In particular, if the tree was generated when
calculating node support by both ultrafast bootstrap and SH approximate likelihood
ratio test, IQTREE writes the support as the node label in the Newick file in the
following format: "SH-aLRT/uBS". In other words, it is a string with two values
separated by a slash. The first value is the SH-aLRT support (as a percentage)
and the second value is the ultrafast bootstrap support (also as a percentage).

The output is a data.frame with each row representing an internal node, with
information that can be used to plot support information layers on a tree.

Some use cases

1. Comparing distances
2. Plotting node support
3. Getting group level stats

library(harrietr)
library(ggplot2)
data("woodmouse")
data("woodmouse_iqtree")
dist_df <- dist_long(aln = woodmouse, order = woodmouse_iqtree$tip.label, dist = "N", tree = woodmouse_iqtree)
ggplot(dist_df, aes(x = dist, y =  evol_dist)) + 
  geom_point() + stat_smooth(method = 'lm') +
  ylab("Evolutionary distance") +
  xlab("SNP distance")

library(ggtree)
library(dplyr)
library(harrietr)
data("woodmouse_iqtree")
p1 <- ggtree(woodmouse_iqtree)
node_support <- get_node_support(woodmouse_iqtree)
p1 + 
  geom_point(data = node_support %>% dplyr::filter(`SH-aLRT` >= 50 & uBS >= 50), aes(x = x, y = y), colour = 'darkgreen', size = 3) +
  geom_point(data = node_support %>% dplyr::filter(`SH-aLRT` >= 50 & uBS >= 50), aes(x = x, y = y), colour = 'darkgreen', size = 5, pch = 21) +
  geom_point(data = node_support %>% dplyr::filter(`SH-aLRT` >= 50 & uBS >= 50), aes(x = x, y = y), colour = 'darkgreen', size = 7, pch = 21)

library(ggplot2)
library(dplyr)
library(harrietr)
data("woodmouse")
data("woodmouse_iqtree")
data("woodmouse_meta")
dist_df <- dist_long(aln = woodmouse, order = woodmouse_iqtree$tip.label, dist = "N", tree = woodmouse_iqtree)
dist_df <- add_metadata(dist_df, woodmouse_meta, isolate = 'SAMPLE_ID', group = 'CLUSTER', remove_ind = TRUE)
dist_df %>%
  dplyr::group_by(CLUSTER) %>%
  dplyr::summarise(q50 = median(dist),
  q25 = quantile(dist, prob = c(0.25)),
  q75 = quantile(dist, prob = c(0.75)),
  min_dist = min(dist),
  max_dist = max(dist)) %>%
  ggplot( aes( x = CLUSTER, y = q50)) +
  geom_errorbar( aes(ymin = q25, ymax = q75),width = 0.25 ) +
  geom_point(size = 3, colour = 'darkred') +
  geom_point( aes( y = min_dist), colour = 'darkgreen', size = 3) +
  geom_point( aes( y = max_dist), colour = 'darkgreen', size = 3) +
  ylab("Pairwise SNP difference") +
  xlab("Groups")