R/Evenness.R

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

#' Evenness function for FCM data
#'
#' This function calculates Pareto evenness from FCM data 
#' @param x flowbasis object generated by flowBasis()
#' @param d Rounding factor for density values. Defaults to 4.
#' @param plot Make plot of diversity values? Defaults to FALSE.
#' @param n Number of replicates. Defaults to 1.
#' @keywords evenness, fcm, alpha
#' @examples
#' ## Short example
#' 
#' # Load precomputed fingerprint object
#' data(CoolingTower)
#' 
#' # Calculate diversity values
#' Evenness(CoolingTower, plot=TRUE)
#' 
#' ## 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
#'  Evenness(fbasis, plot=TRUE)
#' @export

Evenness <- function(x, d = 3, n = 1, plot = FALSE) {
  x <- x@basis/apply(x@basis, 1, max)
  AUC = as.numeric(matrix(nrow = length(x[, 1]), ncol = 1))
  for (i in 1:nrow(x)) {
    AUC[i] = MESS::auc(cum_Richness(x[i, ], d = d)[, 1], cum_Richness(x[i, ], d = d)[, 2])
    AUC[i] = 1 - (AUC[i] - 0.5)/(0.5)
  }
  if (n > 1) {
    results = matrix(nrow = length(x[, 1])/3, ncol = 2)
    results[, 1] = trip(AUC, n)[, 1]
    results[, 2] = trip(AUC, n)[, 2]
    results = data.frame(results)
    colnames(results) = c("Evenness", "sdev")
    rownames(results) = trip_col(attr(x, "dimnames")[[1]], n)
  } else {
    results = matrix(nrow = length(x[, 1]), ncol = 1)
    results[, 1] = AUC
    results = data.frame(results)
    colnames(results) = c("Evenness")
    rownames(results) = attr(x, "dimnames")[[1]]
  }
  if (plot == TRUE) {
    graphics::plot(results$Evenness, pch = 21, bg = grDevices::adjustcolor("blue", 0.7), col = grDevices::adjustcolor("blue", 
                                                                                                                      0.7), cex = 1.5, las = 1, ylab = "Evenness", xlab = "Samples")
  }
  return(results)
  cat(paste0("1 = maximum evenness; 0 = minimum evenness"))
}