Comparative Phylogenetics in R

• Introduction
• Loading and Plotting Phylogenies
• Character Data
• Calculating Phylogenetic Signal

• Introduction
  
• Loading and Plotting Phylogenies
• Character Data
• Calculating Phylogenetic Signal

• Loading tree files
• Parts of tree files
• Displaying trees

The libraries we're using

library(ape)
  #The Swiss army knife of comparative phylogenetics R packages
library(geiger)
  #A package for "macroevolutionary simulation and estimating parameters related to diversification"
library(phytools)
  #Very useful for visualization particularly, great blog support

• Can read in a single tree or multiple trees.
• Takes “.tre” OR “.nex” nexus files
• You can also load in newick files with read.tree().

istioTree <- read.nexus("Data/Fish_12Tax_time_calibrated.tre")

istioTree

Phylogenetic tree with 12 tips and 11 internal nodes.

Tip labels:
    Istiompax_indica, Istiophorus_albicans, Istiophorus_platypterus, Kajikia_albida, Kajikia_audax, Makaira_nigricans, ...

Rooted; includes branch lengths.

str(istioTree)

List of 4
 $ edge       : int [1:22, 1:2] 13 14 15 16 17 17 16 15 18 19 ...
 $ edge.length: num [1:22] 19.98 39.23 6.64 5.47 8.79 ...
 $ Nnode      : int 11
 $ tip.label  : chr [1:12] "Istiompax_indica" "Istiophorus_albicans" "Istiophorus_platypterus" "Kajikia_albida" ...
 - attr(*, "class")= chr "phylo"
 - attr(*, "order")= chr "cladewise"

What does the file look like?

head(istioTree$edge) #Two columns, start node and end node/tip; tips labeled first

     [,1] [,2]
[1,]   13   14
[2,]   14   15
[3,]   15   16
[4,]   16   17
[5,]   17    3
[6,]   17    2

head(istioTree$edge.length) #Gives lengths of branches by edge

[1] 19.976147 39.227748  6.635844  5.473333  8.786535  8.786535

What does the file look like?

istioTree$tip.label #Gives a list of tip labels, order corresponds to edge number

 [1] "Istiompax_indica"           "Istiophorus_albicans"      
 [3] "Istiophorus_platypterus"    "Kajikia_albida"            
 [5] "Kajikia_audax"              "Makaira_nigricans"         
 [7] "Tetrapturus_angustirostris" "Tetrapturus_belone"        
 [9] "Tetrapturus_georgii"        "Tetrapturus_pfluegeri"     
[11] "Trachurus_trachurus"        "Xiphias_gladius"           

is.ultrametric(istioTree) #Checks to see if the tree is ultrametric

[1] TRUE

The basics

plot(istioTree, cex = 4) #Basic view (phylogram)
nodelabels(cex = 4)

Finding tips that match

• not specific to comparative phylogenetics
  

colors <- c("black", "red")
matches <- grep("Kajikia", istioTree$tip.label, value = T)
matches

[1] "Kajikia_albida" "Kajikia_audax" 

Highlighting matching taxa

plot(istioTree, cex = 4, 
     tip.col = ifelse(istioTree$tip.label %in% matches,
                      'red','black'))

Finding the node with a most-recent common ancestor.

#Remember, these are our matches
matches

[1] "Kajikia_albida" "Kajikia_audax" 

#What is the node number?
findMRCA(tree = istioTree, tips = matches, type = "node")

[1] 20

We have given you a phylogeny of clupeiform fishes from Bloom and Lovejoy, 2014. Load it into R and display the tree. The file is called “Clupeiform phylogeny.nex”.

Bonus challenge: Which node number represents the most recent common ancestor for the Anchoa clade?

The node of the Anchoa MRCA is 178.

anchovyTree <- read.nexus("Data/Clupeiform phylogeny.nex")
matches <- grep("Anchoa", anchovyTree$tip.label, value = T)
plot(anchovyTree, cex = 0.5, 
     tip.col = ifelse(anchovyTree$tip.label %in% matches,'red','black'))
findMRCA(tree = anchovyTree, tips = matches, type = "node")

[1] 178

nodelabels(node = 178, text = "MRCA", bg = F, col = "red")

Getting fancy…

plot(istioTree, type = "cladogram", use.edge.length = F, 
     edge.width = 2, font = 4, cex = 2.5)

Getting fancy…

plot(istioTree, type = "fan", rotate.tree = 30, 
     main = "Phylogeny of Billfishes", cex.main = 3.5, 
     label.offset = 2, cex = 2.5, 
     tip.color = c("blue", "red", "red", "yellow2",
                   "yellow2", "orange1", "chartreuse4", 
                   "chartreuse4", "chartreuse4", 
                   "chartreuse4", "black", "purple"))

Try to imitate this tree. It is a subset of the anchovy tree trimmed using extract.clade() like this: subsetAnchovyTree <- extract.clade(anchovyTree, node = 160). Use ?ape::plot.phylo() for help.

Here's how I made that tree. How'd you do?

subsetAnchovyTree <- extract.clade(anchovyTree, node = 160)
plot(subsetAnchovyTree, type = "fan", rotate.tree = 45, open.angle = 180, 
     main = "Phylogeny of Anchovies", cex.main = 5, col.main = "green",
     label.offset = 2, cex = 2.5, tip.color = c("blue", "red"))

• Introduction
• Loading and Plotting Phylogenies
• Character Data
• Calculating Phylogenetic Signal

• Loading data
• Uniting with tree
• Plotting on tree

Multiple characters in a .csv file

characterTable <- read.csv("Data/CodingTableThresh95.csv", row.names = 1)
characterTable[,1:4] #Showing the data

                           bathmin bathmax SSSmin SSS8max
Istiompax_indica             -7000    -485   3416    3651
Istiophorus_albicans         -6275      -1   3393    3741
Istiophorus_platypterus      -7272      -1   3359    3639
Kajikia_albida               -5750   -1267   3561    3681
Kajikia_audax                -7000      -1   3451    3638
Makaira_nigricans            -5750   -1107   3437    3694
Tetrapturus_angustirostris   -7000    -841   3425    3649
Tetrapturus_belone           -1950     -95   3677    3800
Tetrapturus_georgii          -6414   -3392   3458    3741
Tetrapturus_pfluegeri        -6693   -2112   3437    3741
Xiphias_gladius              -6738    -748   3337    3748

Combining the tree and the character data

treeWData <- treedata(istioTree, characterTable, sort = T)

Warning in treedata(istioTree, characterTable, sort = T): The following tips were not found in 'data' and were dropped from 'phy':
    Trachurus_trachurus

  #Trims out taxa that are missing from the tree or the character table

Original tree versus trimmed tree

Combining the tree and the character data

knitr::kable(treeWData$data[,1:4]) #Here's the data

You've already read in the anchovy phylogeny. Now read in the character table, Bloom_etal_2018_Reduced_Dataset.csv, from Bloom et al., 2018.

How many taxa are represented in the phylogeny and the character dataset?

• Discrete traits
• Continuous traits

Plotting discrete traits on tips.

library(RColorBrewer) # Accessory package with better color

#Entering in discrete tip values
discChar <- c(1, 1, 1, 3, 3, 1, 2, 2, 2, 2, 1)

#Plotting
plot(treeWData$phy, main = "Preferred prey", cex.main = 4, 
     label.offset = 2.5, no.margin = F, cex = 3)

#Setting up colors.
reconCol <- brewer.pal(3, "RdYlBu")

#Plot tips
tiplabels(text = rep(" ", 11), tip = seq(1,11,1), cex = 3, 
          frame = "circle", bg = reconCol[discChar])

#Plot legend
legend(x = 0, y = 5, fill = reconCol, bty = "n", cex = 3, 
       legend = c("Sardines", "Pompano", "Hamburgers"))

You have a tree united with characters. It's plotting time!

Your challenge: Color-code the tip labels according to whether a species is diadromous or not.

# Set up colors
labelCols <- c("red", "blue");
diadLabs <- labelCols[as.factor(anchovyTreeWData$data[,3])]

#Plotting
plot(anchovyTreeWData$phy, main = "Diadromy", cex.main = 4, 
     no.margin = F, cex = 3, 
     label.offset = 2.5)

tiplabels(text = rep(" ", length(anchovyTreeWData$phy$tip.label)), 
          cex = 1, frame = "circle", bg = diadLabs)

#Plot legend
legend(x = 0, y = 50, bty = "n", cex = 2, 
       fill = labelCols, 
       legend = c("Diadromous", "Non-diadromous")

Plotting continuous characters using a color gradient

library(plotrix); #Contains a useful way of producing a color legend

#Set up colors and assigning values
colPal <- brewer.pal(9, "Reds")
normedTips <- round((treeWData$data[,8]-min(treeWData$data[,8]))/
                           (max(treeWData$data[,8])-min(treeWData$data[,8]))*8,0)
reconCol <- colPal[normedTips+1]
reconCol

 [1] "#FB6A4A" "#FC9272" "#FC9272" "#FEE0D2" "#FC9272" "#FEE0D2" "#FB6A4A"
 [8] "#FFF5F0" "#FEE0D2" "#FCBBA1" "#67000D"

Plotting continuous characters using a color gradient

# Do the plot
plot(treeWData$phy, main = "Niche breadth", label.offset = 2.5, 
     no.margin = F)

#Label Tips
tiplabels(text = rep(" ", 11), seq(1,11,1), 
          cex = 0.65, frame = "circle", bg = reconCol[1:11])

#Making a key for the plot
col.labels<-c("Narrow", "Broad")
color.legend(11,0,13,4, col.labels, 
             colPal,gradient="y") #Numbers are plot coordinates

You guessed it. Time for a continuous character plot.

Your challenge: Color-code the tip labels according to log body size.

Bonus challenge: Add another tip label showing trophic position. (Hint: check out ?tiplabels())

# Set up colors
logBody <- as.numeric(anchovyTreeWData$data[,1])
colPalBody <- rev(brewer.pal(9, "Reds"))
normedTipsB <- round((logBody-min(logBody))/
                           (max(logBody)-min(logBody))*8,0)
tipColBody <- colPalBody[normedTipsB+1]

trophPos <- as.numeric(anchovyTreeWData$data[,2])
colPalTroph <- rev(brewer.pal(9, "Blues"))
normedTipsT <- round((trophPos-min(trophPos))/
                           (max(trophPos)-min(trophPos))*8,0)
tipColTroph <- colPalTroph[normedTipsT+1]

#Plotting
plot(anchovyTreeWData$phy, main = "Body Size and Trophic level", 
     cex.main = 4, no.margin = F, cex = 3, label.offset = 5)

tiplabels(cex = 3, adj = 1, pch = 21, bg = tipColBody)

tiplabels(cex = 3, adj = 3.5, pch = 21, bg = tipColTroph)

#Making keys for the plot
col.labels<- round(c(min(
  as.numeric(anchovyTreeWData$data[,1])), 
  max(as.numeric(anchovyTreeWData$data[,1]))),2)
text(5,46, "Body Size", cex = 3)
color.legend(0,40,4,44,col.labels,rev(colPalBody),
             gradient="y", cex = 3, align = "rb")

col.labels<-round(c(min(
  as.numeric(anchovyTreeWData$data[,2])), 
  max(as.numeric(anchovyTreeWData$data[,2]))),2)
text(5,36, "Trophic Level", cex = 3)
color.legend(0,30,4,34,col.labels,rev(colPalTroph),
             gradient="y", cex = 3, align = "rb")

• Introduction
  
• Loading and Plotting Phylogenies
• Character Data
• Modeling Character Evolution
• Calculating Phylogenetic Signal

• The basic function
• Setting up for the loop
• The loop
• The results

• Can evaluate Pagel's Lambda or Blomberg's K
• Lambda: scaling for correlations between species relative to correlation expected under Brownian evolution
• K: scaled ratio of among-species variance to contrasts variance

phytools::phylosig(treeWData$phy, x=treeWData$data[,1],method = "K")

[1] 0.1342713

sigTable <- matrix(data = NA, ncol = 4, nrow = ncol(treeWData$data))
rownames(sigTable) <- colnames(treeWData$data)
colnames(sigTable) <- c("Lambda", "Lambda P-value", "K", "K P-value")

Pagel's Lambda, Blomberg's K, and their respective P-values for each character

count <- 1
while(count <= ncol(treeWData$data)){
  temp <- phytools::phylosig(treeWData$phy, treeWData$data[,count], 
                             method = "lambda", test = T); #Calculate Pagel's lambda
  sigTable[count,1] <- temp$lambda
  sigTable[count,2] <- temp$P
  temp <- phytools::phylosig(treeWData$phy, treeWData$data[,count], 
                             method = "K", test = T); #Calculate Blomberg's k
  sigTable[count,3] <- temp$K
  sigTable[count,4] <- temp$P
  count <- count + 1
}

Et voila!

knitr::kable(sigTable)

What is the phylogenetic signal for log body size and trophic position in anchovies?

Step 1: Set up the data and table for results.

# Data
chovyCharTab <- cbind(as.numeric(anchovyTreeWData$data[,1]),as.numeric(anchovyTreeWData$data[,2]))
rownames(chovyCharTab) <- rownames(anchovyTreeWData$data)
colnames(chovyCharTab) <- colnames(anchovyTreeWData$data)[1:2]

#Set up a table for results
sigTable <- matrix(data = NA, ncol = 4, nrow = ncol(chovyCharTab))
rownames(sigTable) <- colnames(chovyCharTab)
colnames(sigTable) <- c("Lambda", "Lambda P-value", "K", "K P-value")

Step 2: Loop through each column.

#Loop to calculate Blomberg's K and P-value for each character
count <- 1
while(count <= ncol(chovyCharTab)){
  temp <- phytools::phylosig(anchovyTreeWData$phy, chovyCharTab[,count], 
                             method = "lambda", test = T); #Calculate Pagel's lambda
  sigTable[count,1] <- temp$lambda
  sigTable[count,2] <- temp$P
  temp <- phytools::phylosig(anchovyTreeWData$phy, chovyCharTab[,count], 
                             method = "K", test = T) #Calculate Blomberg's k
  sigTable[count,3] <- temp$K
  sigTable[count,4] <- temp$P
  count <- count + 1
}

Step 3: How'd we do?

knitr::kable(sigTable)

Analysis of Phylogenetics and Evolution with R

• By Emmanuel Paradis
  • Lead author of ape package
• http://ape-package.ird.fr/APER.html

Phytools

• Another CompPhy R package popularly used
• http://blog.phytools.org/