Skip to content

Instantly share code, notes, and snippets.

@darmitage darmitage/ChaoPD
Created May 15, 2012

Embed
What would you like to do?
Chao's phylogenetic diversity
require(picante)
require(gtools)
require(reshape)
###########################################################
### Utility functions ###
###########################################################
internal2tips.self = function (phy, int.node){
#require(picante); require(ape)
Ntaxa = length(phy$tip.label)
Nnode = phy$Nnode
if ((Ntaxa + Nnode - 1) != nrow(phy$edge)) {
print("tree structure error")
break
}
if (mode(int.node) == "character"){
nodes = which(phy$node.label == int.node) + Ntaxa
}else nodes = int.node
tips = c()
repeat {
nodes = phy$edge[which(phy$edge[, 1] %in% nodes), 2]
if (length(nodes) == 0)
break
tips = c(tips, nodes)
}
tips = tips[tips <= Ntaxa]
if( int.node <= Ntaxa & length(tips) == 0 ){
tips = int.node
}
tips = phy$tip.label[tips]
return(tips)
}
nodes2root <-
function(phy,int.node) {
Ntaxa = length(phy$tip.label) + phy$Nnode
Nnode = phy$Nnode
tips = c()
nodes = int.node
repeat {
nodes = phy$edge[which(phy$edge[,2]%in%nodes),1]
if (length(nodes)==0) break
tips = c(tips,nodes)
}
tips = tips[tips<=Ntaxa]
return(c(rev(tips),int.node))
}
###########################################################
### Calculates phylogenetic diversity profiles ###
###########################################################
ChaoD = function(comm, phy){
#Begin
phy1 <- prune.sample(comm, phy)
samp1 <- t(comm)
samp1 <- samp1[mixedsort(rownames(samp1)),]
taxa = nrow(samp1)
rownames(samp1) <- rep(1:taxa)
sums <- subset(rowSums(samp1), rowSums(samp1) > 0)
samp1 <- samp1[names(sums),]
samples = ncol(samp1)
taxa = nrow(samp1)
phy1$tip.label <- rep(1:length(phy1$tip.label))
phy1 = prune.sample(t(samp1),phy1)
Bnode = taxa + 1
#Measures overall root-node length for each taxon
Li = t(t(dist.nodes(phy1)[1:taxa, Bnode]))
#Pull individual branch lengths
Lb <- cbind(phy1$edge[,2], phy1$edge.length)
Lb <- rbind(Lb, c(Bnode,0))
Lb <- Lb[order(Lb[,1]),]
#Build taxa abundance matrix
p=matrix(0,taxa,samples)
for (k in 1:samples){
p[,k]<-samp1[,k]/sum(samp1[,k])
}
#Compute T, the total average branch length for a community
That=matrix(0,1,samples)
for (k in 1:samples){
That[,k] <-sum(Li * p[,k])
}
#Compute Pb, the relative abundances of branches and Lbk, the matrix of branches/samples
Pb <- matrix(NA, length(Lb[,2]), ncol = samples)
Lbk <- matrix(NA, length(Lb[,2]), ncol = samples)
for (k in 1:samples){
for(i in 1:length(Lb[,1])){
branch <- internal2tips.self(phy1,Lb[i,1])
j = Lb[i,1]
Pb[j,k] <- sum(p[branch,k])
}}
#Compute Lb, matrix of lengths of branches present in a community
for (k in 1:samples){
Sk <- as.numeric(which(samp1[,k] != 0))
for(i in 1:length(Lbk[,1])) {
Lbk[i,k] <- ifelse(Lb[i,1] %in% nodes2root(phy1, Sk), Lb[i,2], 0)
}}
#Initialize q values
lenq = 50
qq <- seq(length=lenq, from=0, by=.11)
# Initialise the Diversity matrix to zero
Dqz = matrix(0, lenq ,samples)
# Loop to calculate the diversity Dqz for each value of q (iq) and each sample (k)
for (k in 1:samples) {
for (iq in 1:lenq) {
q<-qq[iq]
for (zpi in 1:length(Pb[,k])){
if (Pb[zpi,k]>0)(
Dqz[iq,k]<-Dqz[iq,k] + (Lbk[zpi,k]/That[,k]) * Pb[zpi,k]^q)
}
Dqz[iq,k] <- Dqz[iq,k]^(1/(1-q));
}}
results = Dqz
return(results)
}
###########################################################
### ChaoRarefy- diversity profiles rarefied ###
###########################################################
#This function takes your sample and phylogeny and rarefies to account for differences in sample sizes
# Y is usually the smallest sample in the community = min(rowSums(samp))
# X is the number of sampled communities to average the results over (should be >100)
ChaoRarefy=function (samp,phy,Y,K){
Res<-matrix(NA,ncol=3,nrow=dim(samp)[1])
rowsums<-rowSums(samp)
df.list<-vector("list",K)
for(a in 1:K){
print(a)
samplematrix1<-samp
rare_sample<-rrarefy(samplematrix1,Y)
rare_sample = rare_sample[,colSums(rare_sample)>0]
dqz<-ChaoD(rare_sample,phy)
df.list[[a]]<-dqz
} # close for loop
Res<-Reduce("+",df.list)/length(df.list)
return(Res)
}
###########################################################
### An example ###
###########################################################
data(phylocom)
samp = phylocom$sample
phy = phylocom$phylo
ChaoRarefy(samp = samp, phy = phy, Y = min(rowSums(samp)), K = 100)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.