#######################################################################################################
#	Program rarefaction samples communities in order to compare phylogenetic diversity across microbial 
#		communities with different sequencing depths
#	samplematrix: community matrix - rows are samples, columns taxa. see Picante manual for more details
#	Y number of taxa in smallest dataset
#	K is the number times the communities are resampled
#	requires packages VEGAN and PICANTE
#
#  if you have loaded VEGAN and you still get this error message: Error: could not find function "rrarefy" you need
#		to update your VEGAN package.
#  PD is unrooted, so you may not have a community composed of 1 taxon.
#######################################################################################################

PDrarefy=function (samplematrix,tree,Y,K){

pd_res<-matrix(NA,ncol=K,nrow=dim(samplematrix)[1])
rich_res<-matrix(NA,ncol=K,nrow=dim(samplematrix)[1])
Res<-matrix(NA,ncol=3,nrow=dim(samplematrix)[1])

rowsums<-rowSums(samplematrix)
##reads_nodups<-reads_nodups-rowsums

df.list<-vector("list",K)

for(a in 1:K){
print(a)
##samplematrix1<-cbind(samplematrix,reads_nodups)
samplematrix1<-samplematrix
rare_sample<-rrarefy(samplematrix1,Y)


#colnames(rare_sample)[dim(rare_sample)[2]]<-"Yx2Apop2"

#rare_sample[,dim(rare_sample)[2]]<-1
pd<-pd(rare_sample,tree,include.root=FALSE)
pd_res[,a]<-as.numeric(pd[,1])
rich_res[,a]<-as.numeric(pd[,2])

df.list[[a]]<-pd

} # close for loop

#print(pd_res)
#print(rich_res)


Res<-Reduce("+",df.list)/length(df.list)

#colnames(Res)<-c("sample","phylogenetic diversity","taxonomic richness")

return(Res)
####################################################################
# pull out the number of species in each site the sample list      #
####################################################################
}#close program