# 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*

```{r }
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`.


```{r }
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

```{r }
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)
```

