demo/custom-PCA.R

In opm: Analysing Phenotype Microarray and Growth Curve Data

#' # Analysing phenotype microarray data: drawing principal components
#'
#' Using **BiodiversityR** and **opm** to generate `PCA` bi-plots of estimated
#' parameters. An advantage of such plots is that the overall differences
#' between the organisms can be plotted together with indications of the effect
#' of certain substrates on these differences. Furthermore, the data are reduced
#' to their independent components before plotting. Depending on the biological
#' question, this might be more interesting than drawing heat maps. Some
#' **BiodiversityR** functions do not always seem to work as intended, see the
#' comments below.
#'
#' Author: *Markus Goeker*, with assistance by *Pia Wuest*


library(BiodiversityR)
library(opm)


#' ## Define auxiliary function


#' ### `PCA` bi-plot of PM data with **BiodiversityR**
#'
#' This function is intended to operate on numeric matrices as generated by
#' `opm::extract`. It uses methods from the **vegan** and **BiodiversityR**
#' packages.
#'
#' * `x` Matrix as exported by `opm::extract`.
#' * `group` Logical scalar. Assign symbols by group? This causes secondary
#'   symbols to be plotted over the primary symbols (which should best be
#'   invisible in that case).
#' * `circle` Logical scalar. Draw equilibrium circle?
#' * `text.col` Used for drawing descriptor names.
#' * `points.col` Used for drawing the primary object symbols.
#' * `arrow.col` Used for drawing descriptor arrows.
#' * `legend.x` Passed to graphics::legend in non-interactive use.
#' * `legend.y` Likewise.
#' * `lwd` Used for the secondary object symbols.
#' * `cex` Used for drawing descriptor names.
#' * `scaling` Passed to `stats::biplot`.
#' * `...` Optional arguments passed to `stats::biplot`.
#'
#' The function invisibly returns an object of class `pca`.
#'
custom_pca <- function(x, group = TRUE, circle = scaling == 1,
    text.col = "darkgrey", points.col = if (group)
      "white"
    else
      "black", arrow.col = "red", legend.x = "topleft", legend.y = NULL,
    lwd = 2, cex = 0.4, scaling = 1, ...) {
  result <- vegan::rda(x) # conducts a PCA in this setting
  result.plot <- biplot(x = result, col = c(points.col, arrow.col),
    scaling = scaling, ...)
  if (group) {
    groups <- data.frame(.N = if (is.null(attr(x, "row.groups")))
        rownames(x)
      else
        attr(x, "row.groups"))
    if (can.click <- dev.capabilities()$locator)
      message("click on the preferred position of the legend")
    # Problem: in interactive use the returned colour names have been incorrect!
    # This did not seem to work on an X11 device! PDF was OK but see below.
    cols.used <- BiodiversityR::ordisymbol(result.plot, groups, ".N",
      legend = can.click, lwd = lwd)
    if (!can.click)
      legend(legend.x, legend.y, legend = levels(groups$.N),
        text.col = cols.used)
  }
  text(result.plot, "species", col = text.col, cex = cex)
  # Another problem: the circle does not appear in the PDF! Only X11 worked.
  if (circle)
    BiodiversityR::ordiequilibriumcircle(result, result.plot)
  invisible(result)
}


#' ## Create `PCA` bi-plot for the `vaas_4` example data

x <- opm::extract(vaas_4, as.labels = list("Strain"),
  as.groups = list("Species"))
x.pca <- custom_pca(x)


detach("package:BiodiversityR", unload = TRUE)
detach("package:vegan", unload = TRUE)