R/plotHeatmap.R

In e-myers/rnaseq: Process, Analyze and Visualize RNA-seq Data

#' Read density per region - heatmap
#'
#' Create heatmap with summary profile from matrix generated by computeMatrix()
#' @param matrixFiles Character - List of mat.gz files generated by computeMatrix()
#' @param deeptoolsPath String - Path to deeptools directory
#' @param outDest String - Directory where eps files should be saved
#' @param outSuffix String - will be appended to original filename
#' @param textMatrix Logical - Whether to write matrix to text file
#' @param interpolationMethod String - 'nearest' works for our data so far
#' @param dpi Numeric -
#' @param sortRegions String -
#' @param averageTypeSummaryProfile String -
#' @param missingDataColor String -
#' @param colorMap String -
#' @param zMin Numeric -
#' @param zMax Numeric -
#' @param heatmapHeight Numeric -
#' @param heatmapWidth Numeric -
#' @param startLabel String -
#' @param endLabel String -
#' @param samplesLabel Character -
#' @param regionsLabel Charater -
#' @details Take a list of mat.gz files containing read density for a given set of regions,
#' created by the deeptools computeMatrix tool (called by computeMatrix() in this package),
#' and use deeptools plotHeatmap to create an eps with a heatmap and summary profile for each one.
#' Note that you can't give plotHeatmap a missingDataColor value using #ffffff format when calling
#' it through R.  So far as I can tell the # is getting interpreted as a comment, and it thinks you're
#' not giving it anything.
#' The latest version of deeptools by default uses interpolation method "bilinear".
#' However, it fails for our data (at least the TTX and the small RNA), whereas "nearest" works.
#' TIME: ~10s per matrix file
#' COMMAND LINE EXAMPLE (if you want to play with the parameters while looking at just one file, this might be easiest):
#' plotHeatmap -m samplename_regionsName.mat.gz -out samplename_regionsName.eps
#' IMPROVE: If given sample and region labels, some kind of check to see that they correspond to what's
#' in the matfile somehow.  And check the length of the character vectors to make sure it matches the
#' length of matrixFiles.
#' @examples
#' mats = list.files(pattern='.mat.gz')
#' plotHeatmap(mats, samplesLabels=sub('.*_', '', sub('_Aligned.*', '', matfiles)))
#' @author Emma Myers
#' @export

plotHeatmap = function(matrixFiles, deeptoolsPath="", outDest="./", outSuffix="", textMatrix=FALSE,
                       interpolationMethod=NULL, dpi=NULL, sortRegions=NULL,
                       averageTypeSummaryPlot=NULL, missingDataColor=NULL, colorMap=NULL,
                       zMin=NULL, zMax=NULL, heatmapHeight=NULL, heatmapWidth=NULL,
                       startLabel=NULL, endLabel=NULL, samplesLabels=NULL, regionsLabel=NULL) {

    # Check arguments
    if ( !all(file.exists(matrixFiles)) ) {
        writeLines("Missing input files:")
        writeLines(matrixFiles[which(!file.exists(matrixFiles))])
        stop("Missing input file(s).  See above.")
    }
    if (! deeptoolsPath == "") {deeptoolsPath = dir_check(deeptoolsPath)}
    outDest = dir_check(outDest)

    counter = 1

    # For each matrix file
    for (m in matrixFiles) {

        writeLines("\n")

        # Make sure file exists and ends in .mat.gz
        if ( ! file_checks(m, extension=".mat.gz", verbose=TRUE) ) { next }

        writeLines(paste("Processing file:", m))

        # Define arguments to the plotHeatmap command
        fOut = paste(outDest, sub(".mat.gz", paste(outSuffix, ".eps", sep=""), basename(m)), sep="")
        arguments = c("-m", m, "-out", fOut)
        if (textMatrix) { arguments = c(arguments, "--outFileNameMatrix", sub(".eps", "_heatmap_mat.txt", fOut)) }
        if ( !is.null(interpolationMethod) ) { arguments = c(arguments, "--interpolationMethod", interpolationMethod) }
        if ( !is.null(dpi) ) { arguments = c(arguments, "--dpi", dpi) }
        if ( !is.null(sortRegions) ) { arguments = c(arguments, "--sortRegions", sortRegions) }
        if ( !is.null(averageTypeSummaryPlot) ) { arguments = c(arguments, "--averageTypeSummaryPlot", averageTypeSummaryPlot) }
        if ( !is.null(missingDataColor) ) { arguments = c(arguments, "--missingDataColor", missingDataColor) }
        if ( !is.null(colorMap) ) { arguments = c(arguments, "--colorMap", colorMap) }
        if ( !is.null(zMin) ) { arguments = c(arguments, "-min", zMin) }
        if ( !is.null(zMax) ) { arguments = c(arguments, "-max", zMax) }
        if ( !is.null(heatmapHeight) ) { arguments = c(arguments, "--heatmapHeight", heatmapHeight) }
        if ( !is.null(heatmapWidth) ) { arguments = c(arguments, "--heatmapWidth", heatmapWidth) }
        if ( !is.null(startLabel) ) { arguments = c(arguments, "--startLabel", paste("'", startLabel, "'", sep="")) }
        if ( !is.null(endLabel) ) { arguments = c(arguments, "--endLabel", paste("'", endLabel, "'", sep="")) }
        if ( !is.null(samplesLabels) ) { arguments = c(arguments, "--samplesLabel", samplesLabels[counter]) }
        if ( !is.null(regionsLabel) ) { arguments = c(arguments, "--regionsLabel", paste("'", regionsLabel, "'", collapse=" ")) }


        #################
        # Get output file
        #################
        if ( file_checks(fOut, shouldExist=FALSE, verbose=TRUE) ) {
            writeLines(paste("Creating file ", fOut, "...", sep=""), sep="")
            tStart = proc.time()[3]
            system2(paste(deeptoolsPath, "plotHeatmap", sep=""), args = arguments, stdout = fOut)
            tElapsed = proc.time()[3] - tStart
            writeLines(paste("done (", round(tElapsed/60, digits=2), "m).", sep=""))
            writeLines(paste("Done with file."))
        }

        counter = counter + 1

    }

}