Skip to content

Instantly share code, notes, and snippets.

@darmitage darmitage/NaiiveD
Created Jul 25, 2012

Embed
What would you like to do?
Naiive Diversity Profiles + Rarefaction
require(picante)
require(gtools)
require(reshape)
###########################################################
### For calculation of Naiive (Z = I) diversity profiles ###
###########################################################
NaiiveD = function(comm, phy){
#Begin
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)
p=matrix(0,taxa,samples)
for (k in 1:samples){
p[,k]<-samp1[,k]/sum(samp1[,k])
}
Z <- diag(taxa)
lenq = 50
qq <- seq(length=lenq, from=0, by=.11)
Zp=matrix(0,taxa,samples)
for (k in 1:samples) {
for (i in 1:taxa){
for (j in 1:taxa){
Zp[i,k]<-Zp[i,k]+Z[i,j]*p[j,k]
}
}}
Dqz = matrix(0, lenq ,samples)
for (k in 1:samples) {
for (iq in 1:lenq) {
q<-qq[iq];
for (zpi in 1:length(Zp[,k])){
if (Zp[zpi,k]>0)(
Dqz[iq,k]<-Dqz[iq,k]+ p[zpi,k]*(Zp[zpi,k])^(q-1))
}
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)
NaiiveRarefy=function (samp,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<-NaiiveD(rare_sample)
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
NaiiveD(samp)
NaiiveRarefy(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.