PhyloFeast - measuring diet diversity
Coral Levi
Before you get started:
Format your food journal so it looks like the Phylofeast_ExampleData.csv on D2L
Double check that all the genus names EXACTLY match those found in the Phylofeast_TreeSummary.csv file.
Save your formatted data as a .csv file (File –> Save As).
Read through the methods section in the lab handout to calculate richness, and to learn how diversity is calculated.
Enter your manual counts here:
###YOUR TEXT HERE - DESCRIBE HOW YOU CALCULATED RICHNESS
day1 <- 13
day2 <- 6
day3 <- 5
day4 <- 8
day5 <- 13
day6 <- 14
day7 <- 14
###FILL THIS IN SO YOU HAVE RICHNESS FOR ALL DAYS OF YOUR JOURNAL
#make and display a table of your results
richnessDat <- rbind(day1, day2, day3, day4, day5, day6, day7)
#CONTINUE THE LIST, WITH THE REMAINING DAYS)
colnames(richnessDat) <- "Richness"
richnessDat## Richness
## day1 13
## day2 6
## day3 5
## day4 8
## day5 13
## day6 14
## day7 14- Once you have calculated richness by hand, and learned how to calculate diversity, follow the steps below.
Diversity calculations
#install.packages("vegan") #delete this line after successful download
#Load the vegan package into R
library(vegan)## Warning: package 'vegan' was built under R version 3.3.3## Loading required package: permute## Warning: package 'permute' was built under R version 3.3.3## Loading required package: lattice## This is vegan 2.4-3#Load your formatted data file. Remember it must follow the specific formatting, and use correctly spelled and capitalized genus names.
#setwd("C:/Users/levic/Desktop/Evolution")
getwd()## [1] "C:/Users/levic/Desktop/Evolution"my.food <- read.csv("foodjournal.csv", header=T, row.names=1)
#Convert it to a matrix for diversity calculations
my.food <- as.matrix(my.food)
#Calculate the total times each plant genus was consumed
total <- apply(my.food, MARGIN=2, FUN=sum)
###WHAT DOES MARGIN = 2 DO? CONSULT A HELP FILE IN R
#Margin tells the function to look at the columns of the data and make the calculations using that data.
#Add the totals to the bottom row of the data set
my.food <- rbind(my.food, total)
#See what it looks like
my.food## Paracoffea Triticum Brassica Avena Arachis Olea Solanum Cucumis
## day1 1 2 2 1 1 1 2 0
## day2 0 3 1 0 0 0 1 0
## day3 1 3 0 0 0 0 0 0
## day4 1 4 0 1 1 1 0 0
## day5 1 1 1 1 1 3 2 1
## day6 0 5 1 0 0 1 0 0
## day7 1 2 2 1 0 0 2 2
## total 5 20 7 4 3 6 7 3
## Spinacia Salvia Mangifera Ananas Cocos Papaver Citrus Acer Glycine
## day1 1 0 0 1 2 0 0 0 2
## day2 0 0 0 0 1 0 0 0 3
## day3 1 0 0 0 1 0 0 0 1
## day4 1 0 0 0 1 0 0 0 3
## day5 2 0 0 0 1 0 0 0 4
## day6 1 1 1 1 1 1 2 1 1
## day7 1 1 0 0 1 1 1 1 3
## total 7 2 1 2 8 2 3 2 17
## Cicer Vanilla Solanum.1
## day1 0 1 2
## day2 0 0 1
## day3 0 0 0
## day4 0 0 0
## day5 1 0 1
## day6 1 1 1
## day7 0 1 0
## total 2 3 5#Calculate Simpson's diversity metric per day, and for the week (total).
simpsonDat <- diversity(my.food, index="simpson")
simpsonDat## day1 day2 day3 day4 day5 day6 day7
## 0.9141274 0.7800000 0.7346939 0.8165680 0.8950000 0.8950000 0.9150000
## total
## 0.9108661Phylogenetic diversity
#install.packages('picante') #delete this line after successful download
#Load the picante package into R
library(picante)## Warning: package 'picante' was built under R version 3.3.3## Loading required package: ape## Warning: package 'ape' was built under R version 3.3.3## Loading required package: nlme#Load the ape package into R (you don't need to install it again)
library(ape)
#Import the tree file, which has all the plant genera included
tree <- read.tree("PhylofeastJune2016.tre")
#View the generic tree and add tip labels. You may need to wait a moment for the tip labels to show up.
plot.phylo(tree, cex=.2)
head(tree$tip.label)## [1] "Ericameria" "Rigiopappus" "Euthamia" "Amphiachyris"
## [5] "Gundlachia" "Xylothamia"#Don't worry if you can't read the tip labels yet - we'll prune the treeAs you will notice, this tree is enormous. As humans, we have the potential to consume an amazing diversity of plants! Now we will investigate how much of the diversity you incorporated into your diet this past week.
#Prune the full genus level tree to only include the plant genera in your food journal
tips2drop <- tree$tip.label[which(!tree$tip.label %in% colnames(my.food))]
tree.pruned <- drop.tip(tree, tips2drop)
#Look at this new tree, which includes only the plants in your food journal
plot(tree.pruned)
Calculate Faith’s phylogenetic diversity index (PD), per day and for the week (total). The ‘pd’ function also automatically calculates richness (SR). Check the numbers against the calculations you did earlier. Faith’s PD represents the minimum total length of all the phylogenetic branches required to span a given set of taxa on the phylogenetic tree. A larger value indicates a more phylogenetically diverse plant diet.
pdDat <- pd(my.food, tree.pruned)
pdDat## PD SR
## day1 1289.6061 13
## day2 741.3094 6
## day3 740.4719 5
## day4 874.2314 8
## day5 1214.4149 13
## day6 1556.8600 15
## day7 1558.6748 14
## total 1871.5326 20#Put your different diversity measures into a single tidy table
diversityDat <- cbind(pdDat, simpsonDat)
colnames(diversityDat) <- c("PD", "Richness", "Simpson")
#add your name (or some unique string) as a column to your data
diversityDat$student <- rep("Coral", nrow(diversityDat))
#make a new column for day (which includes total)
diversityDat$day <- rownames(diversityDat)
#write your diversity data to a csv file, copy and paste into google doc
write.csv(diversityDat, file = "my-diversity-data.csv",
quote = F, row.names = F)Diversity comparisons
To analyze the class data, you will want to be familiar with tapply() and subset().
Use code from previous labs to analyze the class data set in an interesting way (download the shared google spreadsheet as a CSV, once everyone has entered their data.)
###INSERT YOUR CODE HERE
###How does my average phylogenetic diversity (as a Vegan) compare to the class?
classdata <- read.csv("PhyloFeast_studentData.csv")
tapply(classdata$PD,
list(classdata$day),
FUN=mean)## day 1 day 2 day 3 day 4 day 5 day 6 day 7
## 807.8024 829.9972 743.4309 681.0823 773.2123 754.0295 768.8237
## day1 day2 day3 day4 day5 day6 day7
## 1033.1990 1035.5006 937.9933 946.5396 978.4583 1047.1253 1117.1633
## total
## 1286.0124avg_pd<-subset(classdata$PD,classdata$day=="total")
student_name<-subset(classdata$student,classdata$day=="total")
class(student_name)## [1] "factor"studentmean<-mean(avg_pd)
par(mar=c(7,5,1,1))
barplot(avg_pd,names.arg=student_name,ylab="Average Phylogentetic Diversity", main="Average Phylogenetic Diversity of Each Individual in the Class' Diet", las=2, cex.names=0.5, col=ifelse(student_name=="Coral","purple","turquoise"))
title(xlab="Student Name", line=5.3)
legend(x=0.2, y=4000, c("Coral", "Other Students"), col =c("purple", "turquoise"), pch=15,pt.cex=1)
abline(h=studentmean,cex=2)
In this figure we see the average phylogenetic diversity for each student’s diet in relation to the plant genera they eat in a week. My value is very similar to the trend overall in the class and is onyl a little higher than the mean for the class. We can see that the average phylogenetic divesity for the class is 1286.012 and mine is about 2000.
caption <- "In this figure we see the average phylogenetic diversity for each student's diet in relation to the plant genera they eat in a week. My value is very similar to the trend overall in the class and is onyl a little higher than the mean for the class. We can see that the average phylogenetic divesity for the class is 1286.012 and mine is about 2000."