Schefferville - transitions; Time profile phylogenetic beta diversity
Tammy L. Elliott
Date: March 26, 2015
R version 3.1.0
# set global chunk options:
library(knitr)
opts_chunk$set(cache=FALSE, fig.align='center')Prepare tree by truncating it at specific times
The following analyzes are based on the three cluster grouping based on environmental distances. I chose to compare only three clusters since the calinski criteria indicates that three are optimal for environmental distances.
All analyzes use Hill Number 1 - the exponent of Shannon’s diversity.
I truncated the tree at 400 million years and at 50 million year intervals after that.
When using the treeSlice function, I chose the tree in the list with the highest number of tip.labels as the clade for that chosen time.
# Truncate phylogenetic tree at different intervals
#Use branching.times(trans.one.tree) to see different branching times.
trans.20<-treeSlice(trans.one.tree, 20, trivial=FALSE)
trans.20<-trans.20[[1]] # this gives spermatophyta
#write.tree(trans.20, file="trans.20.newick")
trans.70<-treeSlice(trans.one.tree, 70, trivial=FALSE)
trans.70<-trans.70[[1]] # this also gives spermatophyta
#write.tree(trans.70, file="trans.70.newick")
trans.120<-treeSlice(trans.one.tree, 120, trivial=FALSE)
trans.120<-trans.120[[2]] # this gives angiosperms
#write.tree(trans.120, file="trans.120.newick")
trans.170<-treeSlice(trans.one.tree, 170, trivial=FALSE)
trans.170<-trans.170[[2]] # this gives angiosperms
#write.tree(trans.170, file="trans.170.newick")
trans.220<-treeSlice(trans.one.tree, 220, trivial=FALSE)
trans.220<-trans.220[[2]] # this gives dicots
#write.tree(trans.220, file="trans.220.newick")
trans.270<-treeSlice(trans.one.tree, 270, trivial=FALSE)
trans.270<-trans.270[[2]] # this gives dicots
#write.tree(trans.270, file="trans.270.newick")
trans.320<-treeSlice(trans.one.tree, 320, trivial=FALSE)
asterids<-trans.320[[2]] # this gives asterids
rosids<-trans.320[[4]]
#write.tree(asterids, file="asterids.newick")
#write.tree(rosids, file="rosids.newick")Remove taxa for each of the time periods so that the data match those in the corresponding tree
The following is sample code showing how I truncated the abundance data frame for each time frame.
abd.sp.allsp.t<-t(abd.sp.allsp)
x<-colnames(abd.sp.allsp.t)
y<-trans.20$tip.label
setdiff(x,y)## [1] "CysMon" "DipCom" "DryExp" "EquArv" "EquSci" "EquSyl" "EquVar"
## [8] "HupApr" "LycAno" "SelSel"sperm.abd.20<-abd.sp.allsp.t[,-match(c("CysMon", "DipCom", "DryExp", "EquArv", "EquSci", "EquSyl", "EquVar", "HupApr","LycAno", "SelSel"), colnames(abd.sp.allsp.t))]Data inputed into following analyzes and set-up of profile analyzes
# The following are the input trees for the analyzes:
#tree #all taxa
#trans.20 this gives spermatophyta
#trans.70 this also gives spermatophyta
#trans.120 this gives angiosperms
#trans.170 this gives angiosperms
#trans.220 this gives dicots
#trans.270 this gives dicots
#asterids this gives asterids
#rosids this gives rosids
# The following are the corresponding data sets
#abd.sp.allsp.t
#sperm.abd.20
#sperm.abd.70
#angio.abd.120
#angio.abd.170
#eudicot.abd.220
#eudicot.abd.270
#asterids.320
#rosids.320
# Set-up time profile analyzes, assume Hill Number 1
part.pr<-phyloDATA2014(tree, abd.sp.allsp)
part<-PhD2014(part.pr$pAbun, part.pr$pLength, 1)
part$pCqN## [1] 0.9090127trans.20<-readLines("trans.20.txt")
part.pr.20<-phyloDATA2014(trans.20, t(sperm.abd.20))
part.20<-PhD2014(part.pr.20$pAbun, part.pr.20$pLength, 1)
part.20$pCqN## [1] 0.9173953Compare among phylogenetic beta diversity values amongst major clades
among<-cbind(part$pCqN, part.20$pCqN, part.70$pCqN, part.120$pCqN, part.170$pCqN, part.220$pCqN, part.270$pCqN,
asterids.320.ph$pCqN)
among## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.9090127 0.9173953 0.9173953 0.9147207 0.8991289 0.8939779 0.8626639
## [,8]
## [1,] 0.853654#dev.new(width=8, height=4)
plot(among[1,], axes=FALSE, type="l",col="black", pch=16, cex=1.25, ylab="Beta diversity", xlab="Time before present (my)",las=1, cex.axis=0.9,
main="Time profile of among group phylogenetic beta diversity", ylim=c(0.45,1),cex.main=1, bty="c")
axis(1,1:8,labels=c("420", "400", "350", "300", "250", "200", "150", "100-asterids"), cex.lab=0.9)
axis(2)
box(bty="l", lwd=3)
Create clusters for each different major clade
The following shows an example of how the data for the clusters were divided.
# Three cluster groupings
abd.sp.groups<- cbind(t(abd.sp.allsp), vasc.env.clust.3$cluster)
colnames(abd.sp.groups)[which(colnames(abd.sp.groups) == "")] <- "Group"
Three.clust.1<-subset(abd.sp.groups,abd.sp.groups[,126]=="1")
Three.clust.2<-subset(abd.sp.groups,abd.sp.groups[,126]=="2")
Three.clust.3<-subset(abd.sp.groups,abd.sp.groups[,126]=="3")
sperm.abd.20.groups<- cbind(sperm.abd.20, vasc.env.clust.3$cluster)
colnames(sperm.abd.20.groups)[which(colnames(sperm.abd.20.groups) == "")] <- "Group"
sperm.20.Three.clust.1<-subset(sperm.abd.20.groups,sperm.abd.20.groups[,116]=="1")
sperm.20.Three.clust.2<-subset(sperm.abd.20.groups,sperm.abd.20.groups[,116]=="2")
sperm.20.Three.clust.3<-subset(sperm.abd.20.groups,sperm.abd.20.groups[,116]=="3")Calculate within clade phylogenetic beta diversity
The following is sample code showing how I calculated the within-clade phylogenetic beta diversity.
# Cluster 1 profile
Three.clust.1.ph.pr<-phyloDATA2014(tree, t(Three.clust.1))
Three.clust.1.ph<-PhD2014(Three.clust.1.ph.pr$pAbun, Three.clust.1.ph.pr$pLength, 1)
Three.clust.1.ph$pCqN## [1] 0.9266243sperm.Three.clust.1.ph.pr<-phyloDATA2014(trans.20, t(sperm.20.Three.clust.1))
sperm.Three.clust.1.ph<-PhD2014(sperm.Three.clust.1.ph.pr$pAbun, sperm.Three.clust.1.ph.pr$pLength, 1)
sperm.Three.clust.1.ph$pCqN## [1] 0.9273743Show phylogenetic beta partitioning within three clusters based on environmental distances
cluster.three<-cbind(Three.clust.3.ph$pCqN, sperm.Three.clust.3.ph$pCqN, angio.Three.clust.3.ph$pCqN, eudicot.Three.clust.3.ph$pCqN, asterids.Three.clust.3.ph$pCqN, rosids.Three.clust.3.ph$pCqN)
among.groups<-cbind(part$pCqN, part.20$pCqN, part.120$pCqN, part.220$pCqN, part.270$pCqN,
asterids.320.ph$pCqN)
#dev.new(width=11.8, height=4)
par(mfrow=c(1,3))
plot(cluster.one[1,], axes=FALSE, col="black", pch=16, cex=1.25, ylab="Phylogenetic beta diversity", xlab="Clade",las=1, cex.axis=0.9,
main="Phylogenetic beta diversity partitioned within \nand among clades for Cluster 1", ylim=c(0.45,1),cex.main=1, bty="c")
axis(1,1:6,labels=c("Vasc.", "Sperm.", "Angio.", "Eudicots", "Asterids", "Rosids"), cex.lab=0.9)
axis(2)
points(among.groups[1,], pch=8, col="grey40", cex=1.25)
box(bty="l", lwd=3)
legend("bottomright", c("Phylogenetic beta diversity (among)", "Phylogenetic beta diversity (within)"), col = c("grey50", "black"), cex=1.25,bty="n", pch=c(8,16))
plot(cluster.two[1,], axes=FALSE, col="black", pch=16, cex=1.25, ylab="Phylogenetic beta diversity", xlab="Clade",las=1, cex.axis=0.9,
main="Phylogenetic beta diversity partitioned within \nand among clades for Cluster 2", ylim=c(0.45,1),cex.main=1, bty="c")
axis(1,1:6,labels=c("Vasc.", "Sperm.", "Angio.", "Eudicots", "Asterids", "Rosids"), cex.lab=0.9)
axis(2)
points(among.groups[1,], pch=8, col="grey40", cex=1.25)
box(bty="l", lwd=3)
legend("bottomright", c("Phylogenetic beta diversity (among)", "Phylogenetic beta diversity (within)"), col = c("grey50", "black"), cex=1.25,bty="n", pch=c(8,16))
plot(cluster.three[1,], axes=FALSE, col="black", pch=16, cex=1.25, ylab="Phylogenetic beta diversity", xlab="Clade",las=1, cex.axis=0.9,
main="Phylotenetic beta diversity partitioned within \nand among clades for Cluster 3", ylim=c(0.45,1),cex.main=1, bty="c")
axis(1,1:6,labels=c("Vasc.", "Sperm.", "Angio.", "Eudicots", "Asterids", "Rosids"), cex.lab=0.9)
axis(2)
points(among.groups[1,], pch=8, col="grey40", cex=1.25)
box(bty="l", lwd=3)
legend("bottomright", c("Phylogenetic beta diversity (among)", "Phylogenetic beta diversity (within)"), col = c("grey50", "black"), cex=1.25,bty="n", pch=c(8,16))