R/G_PCA.R

In neyhartj/gws:

#' Plot the principal components of a genomic relationship matrix
#'
#' @description
#' This function takes a genomic relationship matrix (i.e. one generated by
#' the \code{A.mat} function in \code{rrBLUP}), performs a singular value
#' decomposition, and plots the designated principal components.
#'
#' @param G A genomic relationship \code{matrix} of dimensions n x n, where n
#' is the number of entries.
#' @param x.PC The principal component to plot on the x-axis.
#' @param y.PC The principal component to plot on the y-axis.
#' @param col A \code{factor} of categories by which to color the plot of
#' principal componts. Must be of length n.
#' @param scree A \code{logical} indicating whether to display a scree plot.
#' @param return.data A \code{logical} indicating whether to return the results
#' of the singular value decomposition.
#'
#'
#'
#' @export
#'
G_PCA <- function(G, x.PC, y.PC, col = NULL, scree = FALSE, return.data = FALSE) {

  # Error checking
  if (nrow(G) != ncol(G)) stop("'G' is not a square matrix.")
  if (is.null(col)) {
    col <- factor(rep(1, nrow(G)))
  } else {
    if (length(col) != nrow(G)) stop("The length of 'color.factor' is not the same as the number of row/columns of 'Gmat'.")
  }

  # Class verification
  G <- as.matrix(G)
  col <- as.factor(col)

  # Perform the svd
  svd.out <- svd(G)

  # Retrieve lambda
  lambda <- svd.out$d
  # Calculation the proportion of variation explained by each PC
  prop.var <- lambda / sum(lambda)

  # If scree is true, display the plot
  if (scree) {
    plot(x = 1:length(lambda),
         y = prop.var,
         xlab = "Principal Component",
         ylab = "Variation Explained by PC")

    # End the function
    return()
  }

  # Subset the principal components
  PC.x <- svd.out$u[,x.PC]
  PC.y <- svd.out$u[,y.PC]

  # Create pretty axis limits
  range.y <- range(pretty(range(PC.y)))
  ylim <- c( min(range.y), (max(range.y) + ( diff(range.y) / 4 )) )
  range.x <- range(pretty(range(PC.x)))
  xlim <- c( (min(range.x) - ( diff(range.x) / 6 )), max(range.x) )

  # Plot the PCA
  plot(x = PC.x,
       y = PC.y,
       xlab = paste("PC", x.PC, " (", round(prop.var[x.PC], 3) * 100, "%)", sep = ""),
       ylab = paste("PC", y.PC, " (", round(prop.var[y.PC], 3) * 100, "%)", sep = ""),
       ylim = ylim,
       xlim = xlim,
       pch = 16,
       col = col)
  # Add a legend
  legend("topleft", legend = unique(col), col = as.numeric(unique(col)), pch = 16)

  # Output data if told
  if (return.data) {
    return(svd.out)
  } else {
    return(NULL)
  }

} # Close the function