R/pcHeatmap.R

In blackthornrx/rbnctools: Useful R Functions for RBNC

#' pcHeatmap, plot to quickly find PC correlations
#'
#' Returns a PDF plot of PC by sample metric correlations. Automatically checks
#' all numeric columns of colData(rse) against the top 10 principal components.
#'
#' @param rse Gene RSE, holding phenodata and/or sequencing metrics in colData
#' @param pca Object of prcomp() output, holding PCA information
#' @param plotdir Path of directory to save the plot pdf
#' @param width Integer width of plot, default 10
#' @param height Integer height of plot, default 10
#'
#' @return Creates a PDF plot of PC by sample metric correlations.
#' @export
#'
#' @examples
#' # pcHeatmap(rse_gene, pca, getwd())
#'

pcHeatmap = function(rse, pca, plotdir, width=10, height=10) {

  ## Load RColorBrewer if necessary
  if (!require("RColorBrewer", character.only = TRUE)) {
    install.packages("RColorBrewer", dependencies = TRUE) }
  library(RColorBrewer)

  n = min(10, ncol(pca$x))  # number of PCs to plot

  vars = colnames(colData(rse))[sapply(colData(rse), class) == "numeric"]
  ccPcaMetrics = cor(pca$x[,1:n],
                     as.data.frame(colData(rse)[,vars]))
  tPcaMetrics = ccPcaMetrics/sqrt((1-ccPcaMetrics^2)/(ncol(rse)-2))
  pvalPcaMetrics = 2*pt(-abs(tPcaMetrics), df = ncol(rse)-1)
  colnames(pvalPcaMetrics) = gsub("recount_qc_", "", colnames(pvalPcaMetrics),fixed=TRUE )
  colnames(pvalPcaMetrics) = gsub("recount_seq_qc_", "", colnames(pvalPcaMetrics),fixed=TRUE )

  filepath = file.path(plotdir,"pcs_vs_metrics_heatmap.pdf")
  filepath = gsub("//", "/", filepath)

  pdf(filepath, w=width, h=height)
    logPvalMat = -log10(pvalPcaMetrics)
    logPvalMat[logPvalMat>10] = 10
    theSeq = seq(0,10,by=0.1)
    my.col <- colorRampPalette(c("white","blue"))(length(theSeq))
    print(levelplot(logPvalMat, aspect = "fill",
                    at = theSeq,pretty=TRUE,xlab="",ylab="",
                    scales=list(x=list(rot=90, cex=1.2), y=list(cex=1.2)),
                    panel = panel.levelplot.raster, col.regions = my.col))
    dev.off()
}