R/LAPP.plot.PCA.R

In ash129/LAPP: Ligand Analysis with Physicochemical Properties

#' Plot LAPP PCA
#'
#' \code{LAPP.plot.PCA()} takes in a modified PCA object returned by LAPP.PCA() and plots it.
#'Uses ggplot2.
#'
#' @param pca modified PCA object returned by LAPP.PCA(). Can take up to 8 different HLA alleles.
#'
#' @param type Describes the type of plot to make.
#' Either 'points' (default) for scatter plo,t or 'density' for a 2D KDE density graph.
#'
#' @param dim1 Index of PC to plot along the x-axis. Defaults to 1
#'
#' @param dim2 Index of PC to plot along the y-axis. Defaults to 2
#'
#' @param type How the PCA should be plotted. Either 'points' for scatterplot or
#' 'density' for an automated 2D KDE plot (default).
#'
#' @export

LAPP.plot.PCA = function(pca, dim1 = 1, dim2 = 2, type = "density", 
                         colors = c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")) {
  
  # names for principal components
  PC.labels = paste0("PC.", 1:ncol(pca$ind$coord))
  
  # data for PCA coordinates
  pca.df = data.frame(pca$ind$coord)
  colnames(pca.df) = PC.labels
  pca.df$sequence = as.character(pca$sequence)
  pca.df$HLA = pca$HLA
  
  # default colors from:
  # colorblind palette from http://jfly.iam.u-tokyo.ac.jp/color/
  # colors = c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
  
  # set color
  palette = scale_color_manual(values = colors)
  
  # empty plot
  pca.plot= ggplot()
  
  # labels and design settings for PCA plot
  pca.settings = ggplot() +
    labs(title = paste0("PCA for Peptides by HLA (", length(levels(pca$HLA)), 
                        " HLA, ", length(PC.labels), " features)"),
         x = paste0("PC", dim1, ": ", round(pca$eig$per[dim1], digits = 4), "%"),
         y = paste0("PC", dim2, ": ", round(pca$eig$per[dim2], digits = 4), "%")) +
    theme(title = element_text(face = "bold", angle = 0),
          axis.text = element_text(face = "bold", angle = 0),
          axis.line= element_blank(),
          panel.grid = element_blank(), panel.background = element_blank()) +
    palette +
    guides(color = guide_legend(label.position = "bottom",
                                title.theme = element_text(face = "bold", angle = 0),
                                label.theme = element_text(, angle = 0),
                                override.aes = list(size = 3)))
  
  # show points
  if (type == "points") { 
    
    # scramble rows of PC table so that plotting them looks nicer
    set.seed(129)
    pca.df = pca.df[sample(1:nrow(pca.df), nrow(pca.df)), ]
    
    pca.plot= pca.settings +
      geom_point(data = pca.df, mapping = aes(x = pca.df[, dim1], y = pca.df[, dim2], color = HLA),
                 size = 4, alpha = 0.5)
    
    # show densities generated by automated 2D kernel density estimation
  } else if (type == "density") {
    
    pca.plot= pca.settings +
      geom_density2d(data = pca.df, mapping = aes(x = pca.df[, dim1], y = pca.df[, dim2], color = HLA),
                     size = 1, alpha = 0.8)
    
  }
  return(pca.plot)
}