PCA plot

pca-plot.r

# Plot two PCA components with optional coloring

plot_pca <- function(x, pcs, annotation,
                     size = 4, alpha = 0.8, legend.position = c(1, 0)) {
  suppressPackageStartupMessages(require(ggplot2))
  require(scales)

  # check pcs
  if (missing(pcs)) pcs <- paste0("PC", 1:2)
  stopifnot(length(pcs) == 2)
  if (!all(grepl("PC", pcs))) stop("pcs must be something like: PC1.")

  pc.df <- data.frame(x$x)
  all.pcs <- colnames(pc.df)

  if (!missing(annotation)) {
    stopifnot(is.data.frame(annotation))
    pc.df <- cbind(pc.df, annotation[rownames(pc.df), , drop = FALSE])
    mapping <- aes_string(pcs[1], pcs[2], color = names(annotation)[1])
  } else {
    mapping <- aes_string(pcs[1], pcs[2])
  }

  if (legend.position %in% c("top", "bottom", "left", "right")) {
    legend.justification <- "center"
  } else {
    legend.justification <- c(1, 0)
  }

  pc.var <- percent(pca_prop_var(x)[pcs])
  

  ggplot(pc.df) + mapping +
    geom_hline(yintercept = 0, color = "grey50") +
    geom_vline(xintercept = 0, color = "grey50") +
    geom_point(size = size, alpha = alpha) +
    labs(x = paste0(pcs[1], " (", pc.var[1], ")"),
         y = paste0(pcs[2], " (", pc.var[2], ")")) +
    theme_bw() +
    theme(legend.justification = legend.justification,
          legend.position = legend.position,
          legend.background = element_rect(fill = alpha("grey90", 0.5)))
}


pca_prop_var <- function(x) {
  sdev <- setNames(x$sdev, colnames(x$x))
  sdev^2 / sum(sdev^2)
}