PDcalc: an implementation of the Phylogenetic Diversity (PD) calculus in R

I have started putting my various PD functions together in an R package (PDcalc). You can find a development version of the package here:

https://github.com/davidnipperess/PDcalc

I will keep all current versions of the functions available for download from this site for as long as I can but they will not be updated. All future development will be for the R package.

Thanks to everyone who has been using (and providing feedback on) my functions. I hope you find PDcalc as least as useful.

David

phylocurve: an R function for generating a rarefaction curve of Phylogenetic Diversity

Here's a function I have written for the R statistical environment that calculates a rarefaction curve giving expected phylogenetic diversity (mean and variance) for multiple values of sampling effort. Sampling effort can be defined in terms of the number of individuals, sites or species. Expected phylogenetic diversity is calculated using an exact analytical formulation (Nipperess & Matsen 2013) that is both more accurate and more computationally efficient than randomisation methods. I am providing it for free and without warranty under the GNU General Public License. You need to be familiar with R to use this function. The function also requires that the ape package be installed. To load the function, place the file in your working folder and type ‘source(“phylocurve.R”)’.

UPDATE (8th Jan 2016): I have finally implemented the exact solution for the variance of PD under rarefaction! A special acknowledgement to Florent Mazel for sharing his code with me. The function could probably be more efficient but it does the job and is still substantially faster (and more precise) than Monte Carlo subsampling.

UPDATE (21st Jan 2014): When using sampling without replacement (classic rarefaction), older versions of this function would not work with datasets that have a large number of objects (individuals/sites/species) to be rarefied. This is because phylocurve must calculate the number of possible combinations of each subset of objects in the total and this can be a very large number. It is possible, therefore, to exceed the largest number that R can handle. If that occurred, no warning would be issued but the output would include NAs instead of expected values of Phylogenetic Diversity. I have fixed this problem such that phylocurve should now be able to handle very large numbers (I had to re-learn some high school maths but I got there eventually). It might still be possible to exceed the limits of R but you should have a lot more headroom. If you do still encounter this problem, your options at this point are to: 1) sample with replacement (replace=TRUE); 2) try using a more powerful computer (although this is unlikely to make much difference); or 3) use the much slower phylocurve.perm function instead.

Usage

phylocurve (x, phy, stepm=1, subsampling = “individual”, replace = FALSE)

Arguments

x is a community data table (as in the vegan package) with species/OTUs as columns and samples/sites as rows. Columns are labelled with the names of the species/OTUs. Rows are labelled with the names of the samples/sites. Data can be either abundance or incidence (0/1). Column labels must match tip labels in the phylogenetic tree exactly!

phy is a rooted phylogenetic tree with branch lengths stored as a phylo object (as in the ape package) with terminal nodes labelled with names matching those of the community data table. Note that the function trims away any terminal taxa not present in the community data table, so it is not necessary to do this beforehand.

stepm is the size of the interval in a sequence of numbers of individuals, sites or species to which x is to be rarefied.

subsampling indicates whether the subsampling will be by “individual” (default), “site” or “species”. When there are multiple sites, rarefaction by individuals or species is done by first pooling the sites.

replace is a boolean indicating whether subsampling should be done with or without replacement.

Details

phylocurve takes a community data table and a rooted phylogenetic tree (with branch lengths) and calculates expected mean and variance of Phylogenetic Diversity (PD) for every specified value of m individuals, sites or species. m will range from 1 to the total number of individuals/sites/species in increments given by stepm. Calculations are done using the exact analytical formulae (Nipperess & Matsen, 2013) generalised from the classic equation of Hurlbert (1971). When there are multiple sites in the community table and rarefaction is by individuals or species, sites are first pooled.

Value

phylocurve returns a matrix object of three columns giving the expected PD values (mean and variance) for each value of m.

References

Hurlbert (1971) The nonconcept of Species Diversity: a critique and alternative parameters. Ecology 52: 577-586.
Nipperess & Matsen (2013) The mean and variance of phylogenetic diversity under rarefaction. Methods in Ecology & Evolution 4: 566-572. Manuscript also available from ArXiv.org

>> Download phylocurve.R

Data management and manipulation in R

Here's the materials I used for teaching a module on R for the Genes to Geoscience Research Enrichment Program.

Files for my R module

Zip archive

phylorare: an R function for calculating the rarefied Phylogenetic Diversity of ecological samples

Here's a function I have written for the R statistical environment that calculates expected phylogenetic diversity (under rarefaction) of multiple samples. The function uses an exact analytical formula (Nipperess & Matsen 2013). I am providing it for free and without warranty under the GNU General Public License. You need to be familiar with R to use this function. The function also requires that the ape package be installed. To load the function, place the file in your working folder and type ‘source(“phylorare.R”)’.

UPDATE (21st Jan 2014): When using sampling without replacement (classic rarefaction), older versions of this function would not work with datasets that have a large number of objects (individuals/sites/species) to be rarefied. This is because phylorare must calculate the number of possible combinations of a subset of objects in the total and this can be a very large number. It is possible, therefore, to exceed the largest number that R can handle. If that occurred, no warning would be issued but the output would include NAs instead of expected values of Phylogenetic Diversity. I have fixed this problem such that phylocurve should now be able to handle very large numbers (I had to re-learn some high school maths but I got there eventually). It might still be possible to exceed the limits of R but you should have a lot more headroom. If you do still encounter this problem, your options at this point are to: 1) sample with replacement (replace=T); 2) try using a more powerful computer (although this is unlikely to make much difference); or 3) use a much slower algorithmic solution (such as phylocurve.perm).

Usage

phylorare (x, phy, m=1, subsampling = “individual”, replace =F)

Arguments

x is a community data table (as in the vegan package) with species/OTUs as columns and samples/sites as rows. Columns are labelled with the names of the species/OTUs. Rows are labelled with the names of the samples/sites. Data can be either abundance or incidence (0/1).

phy is a rooted phylogenetic tree with branch lengths stored as a phylo object (as in the ape package) with terminal nodes labelled with names matching those of the community data table. Note that the function trims away any terminal taxa not present in the community data table, so it is not necessary to do this beforehand.

m is the number of individuals, sites or species to which x is to be rarefied.

subsampling indicates whether the subsampling will be by “individual” (default), “site” or “species”. When there are multiple sites, rarefaction by individuals or species is done separately for each site.

replace is a boolean indicating whether subsampling should be done with or without replacement.

Details

phylorare takes a community data table and a rooted phylogenetic tree (with branch lengths) and calculates expected Phylogenetic Diversity (PD) for a given sample size of m individuals, sites or species. Calculations are done using an exact analytical formula generalised from the classic equation of Hurlbert (1971). When there are multiple sites in the community table and rarefaction is either individual or species-based, rarefied PD values are calculated for each site separately. In this case, all PD values will include the root of phy even if all the species in a site share a more recent common ancestor.

Value

phylorare returns a vector object giving the expected PD values of all sample/sites in x.

References

Hurlbert. 1971. The nonconcept of Species Diversity: a critique and alternative parameters. Ecology 52: 577-586.
Nipperess & Matsen. 2013. The mean and variance of phylogenetic diversity under rarefaction. Methods in Ecology & Evolution 4: 566-572. Manuscript also available from ArXiv.org

>> Download phylorare.R

phylodiv: an R function for calculating the phylogenetic diversity of ecological samples

Here's a function I have written for the R statistical environment that calculates the Phylogenetic Diversity (PD) of multiple samples. I am providing it for free and without warranty under the GNU General Public License. You need to be familiar with R to use this function. The function also requires that the ape package be installed. To load the function, place the file in your working folder and type ‘source(“phylodiv.R”)’.

Latest version: 13th December 2012. A couple of small changes to improve efficiency.

Usage

phylodiv (x, phy)

Arguments

x is a community data table (as in the vegan package) with species/OTUs as columns and samples/sites as rows. Columns are labelled with the names of the species/OTUs. Rows are labelled with the names of the samples/sites. Data can be either abundance or incidence (0/1).

phy is a rooted phylogenetic tree with branch lengths stored as a phylo object (as in the ape package) with terminal nodes labelled with names matching those of the community data table. Note that the function trims away any terminal taxa not present in the community data table, so it is not necessary to do this beforehand.

Details

phylodiv takes a community data table and a phylogenetic tree (rooted and with branch lengths) and calculates the Phylogenetic Diversity (PD) of all samples/sites. PD is defined as the total length of all branches spanning a set of terminal taxa representing an ecological sample (Faith, 1992). Please note that, if the common ancestor (node) of the set of taxa of a sample is not the root of the tree, then the set of branches connecting this node to the root are also included in the calculation. Calculations are achieved using the efficient matrix algebra solution of Rodrigues & Gaston (2002).

Value

phylodiv returns a vector giving the Phylogenetic Diversity (PD) of each sample/site in x.

References

Faith DP. 1992. Conservation evaluation and phylogenetic diversity. Biological Conservation 61: 1-10.

Rodrigues A & Gaston KJ. 2002. Maximising phylogenetic diversity in the selection of networks of conservation areas. Biological Conservation 105: 103-111.

>> Download phylodiv.R

The mean and variance of phylogenetic diversity under rarefaction

Erick Matsen and I have just submitted a manuscript to Methods in Ecology and Evolution on the rarefaction of phylogenetic diversity. The manuscript is available to view on ArXiv.org.

phylocurve.perm: an R function for generating a rarefaction curve of Phylogenetic Diversity by randomisation

Here's a function I have written for the R statistical environment that calculates a rarefaction curve giving expected Phylogenetic Diversity (PD) for multiple values of sampling effort. Sampling effort can be defined in terms of the number of individuals, sites or species. This function is a variant of “phylocurve” that uses randomisation rather than an exact analytical formula to calculate an estimate of expected PD.  In addition, this version allows the calculation of estimates for the standard deviation, 95% confidence limits and minimum and maximum values of PD. The accuracy of these estimates is of course dependent on the number of permutations. I am providing it for free and without warranty under the GNU General Public License. You need to be familiar with R to use this function. The function also requires that the ape package be installed. To load the function, place the file in your working folder and type ‘source(“phylocurve_perm.R”)’.

WARNING: this function is much slower (about 10x to 100x) than “phylocurve”. For large datasets, I recommend testing with relatively large values of stepm and small values of perm to gauge the amount of computer time required.

Latest version: 17th August 2011. This version fixes a bug where rarefaction values were calculated incorrectly with community data tables containing species with total abundances of 0.

Usage

phylocurve.perm (x, phy, stepm=1, subsampling = “individual”, replacement = F, perm=1000)

Arguments

x is a community data table (as in the vegan package) with species/OTUs as columns and samples/sites as rows. Columns are labelled with the names of the species/OTUs. Rows are labelled with the names of the samples/sites. Data can be either abundance or incidence (0/1).

phy is a rooted phylogenetic tree with branch lengths stored as a phylo object (as in the ape package) with terminal nodes labelled with names matching those of the community data table. Note that the function trims away any terminal taxa not present in the community data table, so it is not necessary to do this beforehand.

stepm is the size of the interval in a sequence of numbers of individuals, sites or species to which x is to be rarefied.

subsampling indicates whether the subsampling will be by “individual” (default), “site” or “species”. When there are multiple sites, rarefaction by individuals or species is done by first pooling the sites.

replacement is a logical indicating whether subsampling should be done with or without replacement.

perm is the number of iterations of the subsampling routine used to calculate expected PD.

Details

phylocurve.perm takes a community data table and a rooted phylogenetic tree (with branch lengths) and calculates expected Phylogenetic Diversity (PD) for every specified value of m individuals, sites or species. m will range from 1 to the total number of individuals/sites/species in increments given by stepm. Estimates of expected PD as well as variance about that mean are determined by monte carlo randomisations (see Gotelli and Colwell 2001 for a more detailed explanation of the procedure as applied to species richness). When there are multiple sites in the community table and rarefaction is by individuals or species, sites are first pooled.

Value

phylocurve.perm returns a matrix object of seven columns with each row corresponding to a particular value of m. The columns are: 1) the values of m; 2) the mean (expected) values of PD for all randomisations; 3) the standard deviation of PD for all randomisations; 4) the 0.025 quantile (lower 95% confidence limit); 5) the 0.975 quantile (upper 95% confidence limit); 6) the minimum PD; and 7) the maximum PD.

References

Gotelli and Colwell. 2001. Quantifying biodiversity: procedures and pitfalls in the measurement and comparison of species richness. Ecology Letters 4: 379–391.

>> Download phylocurve_perm.R