Sesame St Phylogeny

library(knitr)
knitr::opts_chunk$set(comment="  ", cache=T, warning=F, message=F)
sst.colors <- c(`0`="#00af42",`1`="#fcd000") #from the sign

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(googlesheets4) # read in the data from google sheets
muppets <- read_sheet("1M6hy0z76JMO6aLbA1nZWfkczF7D2QOeuUPnLIQwDedg") %>% 
  select(-Photo) %>% column_to_rownames("Name")

## ! Using an auto-discovered, cached token.

##   To suppress this message, modify your code or options to clearly consent to
##   the use of a cached token.

##   See gargle's "Non-interactive auth" vignette for more details:

##   <https://gargle.r-lib.org/articles/non-interactive-auth.html>

## ℹ The googlesheets4 package is using a cached token for 'icyjava@gmail.com'.

## ✓ Reading from "Sesame Phylo".

## ✓ Range 'Sheet1'.

Heatmap of character states

library(vegan) # make a heatmap of the character states
library(pheatmap)
pheatmap(muppets, col=sst.colors, legend_breaks=c(0,1), legend_labels=c("No", "Yes"))

Infer the tree with parsimony

library(phangorn) #infer the tree with parsimony
mup <- as.phyDat(as.matrix(muppets), type="USER", c(0,1)) #convert character matrix to phyDat format
tree <- pratchet(mup) #find most parsimonious tree with parsimony ratchet (Nixon, 1999)

   [1] "Best pscore so far: 58"
   [1] "Best pscore so far: 58"
   [1] "Best pscore so far: 58"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"
   [1] "Best pscore so far: 56"

tree <- acctran(tree, mup) # assign branch lengths via accelerated transformation algorithm - 
# changes are assigned along branches of a phylogenetic tree as close to the root as possible
(tree.rooted <- midpoint(tree))# root the tree at its midpoint - could use root() to pick outgroup

   
   Phylogenetic tree with 45 tips and 44 internal nodes.
   
   Tip labels:
     Big Bird, Bip Bippadotta, Barkley, Count Von Count, Don Music, Baby Bear, ...
   Node labels:
     1, 0.25, NA, 0.166666666666667, NA, NA, ...
   
   Rooted; includes branch lengths.

Map traits on the tree

library(phytools) #map a character onto the tree

for(traitname in colnames(muppets)) {
  trt <- setNames(muppets[[traitname]],rownames(muppets)) #choose the "Monster" trait
  mtree <- make.simmap(tree.rooted,trt, nsim=50, message=F) #simulate 100 evolutions of the trait on our tree
  dmtree <- describe.simmap(mtree) #summarize the probability of the state at each node
  plot(dmtree, colors=sst.colors) # plot the phylogeny with the trait probabilities at nodes shown as piecharts
  legend("topright", title=paste0(traitname, "?"), legend=c("No","Yes"), fill=sst.colors)
}