R/phylodiv_expected.R

In davidnipperess/PDcalc: An implementation of the Phylogenetic Diversity (PD) calculus in R

#' Expected Phylogenetic Diversity of ecological samples
#' 
#' Calculates the expected Phylogenetic Diversity (PD) of multiple samples 
#' simultaneously based on probabilities of occurrence (or extinction) for each 
#' combination of species and site. Note that this function uses that version of
#' PD that always includes the path to the root of the tree.
#' @param x is the community data given as a \code{data.frame} or \code{matrix} 
#'   with species/OTUs as columns and samples/sites as rows (like in the 
#'   \code{vegan} package). Columns are labelled with the names of the 
#'   species/OTUs. Rows are labelled with the names of the samples/sites. Data 
#'   are expressed as probabilities (ranging 0-1) of either occurrence or 
#'   extinction. Column labels must match tip labels in the phylogenetic tree 
#'   exactly!
#' @param phy is a rooted phylogenetic tree with branch lengths stored as a 
#'   phylo object (as in the \code{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.
#' @param occur is a \code{logical} indicating whether \code{x} is a matrix of 
#'   probabilities of occurrence (\code{TRUE} - default) or extinction 
#'   (\code{FALSE})
#' @details \code{phylodiv.expected} takes probability data and a phylogenetic 
#'   tree (rooted and with branch lengths) and calculates the expected 
#'   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). Expected PD is the summed branch length
#'   expected for each site given the probability of occurrence (or extinction) 
#'   of species (and ancestral branches) at that site (Witting & Loeschcke 1995;
#'   Faith 2013). 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.
#' @return A vector of expected Phylogenetic Diversity (PD) values for each 
#'   sample/site in \code{x}.
#' @importFrom ape drop.tip
#' @references \itemize{\item{Faith DP. 1992. Conservation evaluation and 
#'   phylogenetic diversity. \emph{Biological Conservation} 61: 1-10.} 
#'   \item{Faith DP. 2013. Biodiversity and evolutionary history: useful 
#'   extensions of the PD phylogenetic diversity assessment framework.
#'   \emph{Annals of the New York Academy of Sciences} 1289: 69–89.} 
#'   \item{Witting L. & Loeschcke V. 1995. The optimization of biodiversity
#'   conservation. \emph{Biological Conservation} 71: 205–207.}}
#' @export
#' 


phylodiv.expected <- function (x, phy, occur=TRUE) {
  
  ### step 1: matching taxa and trimming the tree to match the community data
  ### table thus creating a "community tree" (sensu Cam Webb).
  
  taxon_check <- phylomatchr(x,phy)
  
  if(length(taxon_check[[2]])>0) {
    stop('The following taxa in the community data table were not found on the',
         ' tips of the tree: ', paste(taxon_check[[2]], collapse=', '),
         '.\nPlease fix your taxonomy!',
         call.=FALSE)
  }
  
  if(length(taxon_check[[1]])>0) {
    warning(length(taxon_check[[1]]), " taxa were not in x, and were trimmed",
            " from the tree prior to calculation of PD.")
    phy <- drop.tip (phy, taxon_check[[1]])
  }
  
  ### step 2: converting a community tree into a MRP matrix
  
  # A MRP matrix, used in supertree methods, is where the membership of an OTU
  # in a clade spanned by a branch is indicated by a 0 (no) or 1 (yes). Unlike 
  # supertree MRP matrices, our matrix includes terminal branches. using Peter 
  # Wilson's function
  
  phylomatrix <- FastXtreePhylo(phy)$H1
  
  # this script fills a matrix of 0's with 1's corresponding to the incidence of
  # an OTU on a particular branch.
  
  ### step 3: re-ordering the OTUs of the occurrence table and MRP matrix to
  ### match.
  
  phylomatrix <- phylomatrix[order(phy$tip.label), ]
  x <- x[ ,order(colnames(x))]
  
  ### step 4: occurrence or extinction?
  
  x <- as.matrix(x)
  if(occur==FALSE) {
    x <- 1-x
  }
  
  ### step 5: determining branch probabilities and PD
  
  site_pd <- function(x_vector,phylomatrix) {
    probs <- phylomatrix*x_vector
    probs <- 1-apply(1-probs,MARGIN=2,FUN=prod)
    sum(probs*phy$edge.length)
  }
  
  pd <- apply(x,MARGIN=1,FUN=site_pd,phylomatrix)
  
  return (pd)
  
}