R/beta_div_fcm.R

In rprops/Phenoflow_package: Advanced analysis of microbial flow cytometry data

#' Beta diversity function for FCM data
#'
#' This function performs beta diversity analysis on FCM data. 
#' The output can be fed into the plot.beta.fcm() function for visualization.
#' @param x flowbasis object generated by flowBasis()
#' @param d Rounding factor for density values. Defaults to 4.
#' @param INDICES Factor vector indicating the samples to average for the beta-diversity
#' analysis
#' @param dist Distance metric to use in the vegdist() function. Defaults to bray.
#' @param k Number of dimensions to project your samples into. Defaults to 2.
#' @param iter Number of iterations for NMDS analysis. Defaults to 100.
#' @param binary Specify whether the data should be transformed to presence/absence. Defaults to FALSE.
#' @param ord.type Choose between NMDS or PCoA analysis.
#' @keywords betadiversity, fcm
#' @importFrom vegan vegdist metaMDS
#' @importFrom stats cmdscale
#' @examples
#' ## Short example
#' 
#' # Load precomputed fingerprint object
#' data(CoolingTower)
#' 
#' # Calculate diversity values
#' beta <- beta_div_fcm(CoolingTower,ord.type="PCoA")
#' plot_beta_fcm(beta)
#' 
#' ## Full data processing example
#' 
#' # Load raw data (imported using flowCore)
#' data(flowData)
#' # Asinh transform and select parameters of interest (cells were stained with Sybr Green I).
#' flowData_transformed <- flowCore::transform(flowData,`FL1-H`=asinh(`FL1-H`),
#'        `SSC-H`=asinh(`SSC-H`), 
#'        `FL3-H`=asinh(`FL3-H`), 
#'        `FSC-H`=asinh(`FSC-H`))
#' param=c('FL1-H', 'FL3-H','SSC-H','FSC-H')
#' flowData_transformed = flowData_transformed[,param]
#' 
#' # Create a PolygonGate for denoising the dataset
#' # Define coordinates for gate in sqrcut1 in format: c(x,x,x,x,y,y,y,y)
#' sqrcut1 <- matrix(c(8.75,8.75,14,14,3,7.5,14,3),ncol=2, nrow=4)
#' colnames(sqrcut1) <- c('FL1-H','FL3-H')
#' polyGate1 <- flowCore::polygonGate(.gate=sqrcut1, filterId = 'Total Cells')
#' 
#' # Gating quality check
#' flowViz::xyplot(`FL3-H` ~ `FL1-H`, data=flowData_transformed[1], filter=polyGate1,
#'          scales=list(y=list(limits=c(0,14)),
#'          x=list(limits=c(6,16))),
#'          axis = lattice::axis.default, nbin=125, 
#'          par.strip.text=list(col='white', font=2, cex=2), smooth=FALSE)
#'  
#'  # Isolate only the cellular information based on the polyGate1
#'  flowData_transformed <- flowCore::Subset(flowData_transformed, polyGate1)
#'  
#'  # Normalize parameter values to [0,1] interval based on max. value across parameters
#'  summary <- flowCore::fsApply(x=flowData_transformed,FUN=function(x) apply(x,2,max),use.exprs=TRUE)
#'  max = max(summary[,1])
#'  mytrans <- function(x) x/max
#'  flowData_transformed <- flowCore::transform(flowData_transformed,`FL1-H`=mytrans(`FL1-H`),
#'          `FL3-H`=mytrans(`FL3-H`), 
#'          `SSC-H`=mytrans(`SSC-H`),
#'          `FSC-H`=mytrans(`FSC-H`))
#'  
#'  # Calculate fingerprint
#'  fbasis <- flowFDA::flowBasis(flowData_transformed, param, nbin=128, 
#'          bw=0.01, normalize=function(x) x)
#'  
#'  # Calculate diversity
#'  beta <- beta_div_fcm(fbasis, ord.type="PCoA")
#'  plot_beta_fcm(beta)
#' @export

beta_div_fcm <- function(x, d = 4, dist = "bray", k = 2, iter = 100, 
                         ord.type = c("NMDS", "PCoA"), INDICES=NULL,
                         binary = FALSE) {
  x <- x@basis/apply(x@basis, 1, max)
  x <- base::round(x, d)
  if(!is.null(INDICES)){
    x <- by(x, INDICES = INDICES, FUN = colMeans)
    x <- do.call(rbind, x)
  }
    input.dist <- vegan::vegdist(x, method = dist, binary = binary)
    if (ord.type == "NMDS"){ 
      mds.fbasis <- vegan::metaMDS(input.dist, autotransform = FALSE, k = k,
                                   trymax = iter) 
      } else {
        mds.fbasis <- stats::cmdscale(input.dist, k = k, eig = TRUE, add = TRUE)
      }
  return(mds.fbasis)
}