Hello again,
Today I (José Hidasi-Neto) will show you functions that calculate four measures of phylogenetic diversity that were published in Cadotte et al. (2010).
This is the paper:
http://www.ncbi.nlm.nih.gov/pubmed/19903196
http://www.utsc.utoronto.ca/~mcadotte/pubs/Cadotte%20et%20al.%202010%20Ecol.%20Letts.pdf
In my last post, I showed you a measure named Phylogenetic-Abundance Evenness (PAE), which is one of the measures defined by Cadotte et al.. Some people asked me to create functions for the other methods found in the paper, so today I'll show you them. However, instead of explaining and exemplifying (like I did for PAE), I will only present my functions for the Imbalance of Abundances at Higher Clades (IAC), the Entropy of Evolutionary Distinctiveness (H'ed), the Abundance-weighted Evolutionary Distinctiveness (AED), and the Entropy of Abundance-weighted Evolutionary Distinctiveness (H'aed). Then, you can use the phylogeny and abundance files from the example I showed you in my last post to test them. Therefore, this is another reason to read the paper and understand what the measures represent. Ok?
##################
Here we go to R again
##################
# First of all, install and load these packages:
install.packages("ape")
install.packages("picante")
library(ape)
library(picante)
###############################################
# Now let's begin with the arguments (same for all functions):
# First argument: 'file with species abundances' - a data frame with species in columns and their abundances in the first line. You have an example file (and how to open it) at the end of this post.
# Second argument: 'rooted phylogeny' # Please, try to use a well-resolved phylogeny.
###############################################
###############################################
###############################################
Let's begin with Imbalance of Abundances at Higher Clades (IAC) measure:
where ni is each species abundance, n^ is the mean species abundance in the assemblage, and v is the number of nodes.
__________________________
As stated by Cadotte et al. (2010):
"This metric is the relative per-node imbalance in the distribution of individuals, and is always ≥ 0. If IAC = 0 then the distribution of individuals is balanced among higher clades."
__________________________
Here is the IAC function:
IAC<-function(abund,tree){
spnames<-names(abund)
abd<-numeric(length(abund))
for(i in 1:length(abund)){
sp<-spnames[i]
abd[i]<-abs((abund[sp])-mean(unlist(abund)))
}
return(sum(unlist(abd))/tree$Nnode)
}
# IAC(abund,tree)
###############################################
###############################################
##############################################################################################
The following is the Entropy of Evolutionary Distinctiveness (H'ed) :
with its respective equitability (Eed):
__________________________
As stated by Cadotte et al. (2010):
"entropy measure, H'ed, can be calculated, which represents the evolutionary information content of a randomly sampled species from the community"
__________________________
Here is the function to calculate H'ed and its respective Eed:
### Ok, I just used the 'abund' file here to give continuity to the understanding of the other functions =]
Hed<-function(abund,tree){
PD<-sum(evol.distinct(tree,type="fair.proportion")[,2])
EDPD<-evol.distinct(tree,type="fair.proportion")[,2]/PD
Hed<-sum(-(EDPD*log(EDPD)))
Eed<-Hed/log(length(abund))
return(cbind(Hed,Eed))
}
# Hed(abund,tree)
###############################################
###############################################
__________________________
As stated by Cadotte et al. (2010):
"Our final measure assigns an index of ED to species or individuals based on abundances and phylogenetic distances. We derive this metric from the ED method of Isaac et al."
__________________________
Here is the AED function:
AED<-function(abund,tree){
p.tree<-multi2di(tree, random = TRUE)
uninodes<-p.tree$edge[,2][-which(p.tree$edge[,2] %in% 1:length(p.tree$tip.label))]
infonode<-cbind(p.tree$edge,p.tree$edge.length)
nodevalues<-numeric(length(uninodes))
for(i in 1:length(uninodes)){
temptree<-extract.clade(p.tree,uninodes[i])
abundvalues<-abund[1,which(names(abund) %in% temptree$tip.label)]
nodevalues[i]<-(infonode[which(infonode[,2]==uninodes[i]),3])/(sum(abundvalues))
}
nodeabval<-cbind(uninodes,nodevalues)
tipnums<-data.frame(cbind(p.tree$tip.label,1:length(p.tree$tip.label)))
AED<-numeric(length(abund))
for(j in 1:length(abund)){
sp<-tipnums[,1][j]
nodes <- .get.nodes(p.tree,as.vector(sp))
splength<-infonode[which(infonode[,2]==as.numeric(as.vector(tipnums[which(tipnums[,1]==sp),2]))),3]
splength<-splength/(abund[which(colnames(abund)==sp)])
AED[j]<-(sum(nodeabval[which(nodeabval[,1] %in% nodes),2])) + splength
}
names(AED)<-tipnums[,1]
AED<-unlist(AED)
return(AED)
}
###############################################
###############################################
###############################################
Finally, similarly to the ED, we also have the Entropy of the Abundance-weighted Evolutionary Distinctiveness (H'aed):
with its respective equitability (Eaed):
Here is our last function:
Haed<-function(abund,tree){
AEDres<-AED(abund,tree)
PD<-sum(evol.distinct(tree,type="fair.proportion")[,2])
Haed<- -(sum(((AEDres*abund)/PD)*log(((AEDres*abund)/PD))))
Eaed<-Haed/log(sum(abund))
return(cbind(Haed,Eaed))
}
###############################################
###############################################
###############################################
Github link: https://github.com/hidasi/rfunctions
Thank you very much for reading!
See you next time.
0 Comments