R/subVolcano.R

In peterawe/CMScaller: CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models

#' volcano plot
#' @export
#' @description Volcano plot of output from \code{\link{subDEG}}.
#' @param deg output from \code{\link{subDEG}} which may be either a list or
#' a data frame.
#' @param geneID \code{deg} column name to use for labeling top hits.
#' @param lfc a numeric absolute \eqn{log2(fold-change)} threshold.
#' @param padj a numeric, adjusted \eqn{p}-value threshold.
#' \code{length(class)==ncol(emat)}.
#' @param ave a numeric, average expression threshold.
#' @param topN an integer, number of genes to label in plot (selected by
#' largest absolute \eqn{fold-change}).
#' @param classCol a character vector specifying class colors.
#' @param cexText a numeric, scaling factor for labels for top hits.
#' @param ... additional arguments passed to \code{\link[graphics]{plot}} or
#' \code{\link[graphics]{smoothScatter}} or if \pkg{KernSmooth} is available
#' (currently only main, xlab and ylab used).
#' @return one or more volcano plots where horizontal dashed line is the maximum
#' \eqn{p}-value below \code{padj} and vertical lines the \code{lfc}-threshold.
subVolcano <- function(deg, geneID = "rownames",
                       lfc = log2(2), padj = .05, ave=0,
                       topN = 10, cexText=1,
                       classCol = getOption("subClassCol"), ...) {

    ddd <- list(...)
    if (!is.null(ddd$main)) main <- ddd$main else main=""
    if (!is.null(ddd$xlab)) xlab <- ddd$xlab else
        xlab=expression(log[2](fold~change))
    if (!is.null(ddd$ylab)) ylab <- ddd$ylab else
        ylab=ylab=expression(-log[10](italic(p)))

    plotVolc <- function(deg, clCol=classCol) {
        x <- deg$logFC
        y <- -log10(deg$P.Value)
        k <- is.infinite(y)
        # for very small p-values, replace with some random large number
        # stochastic to spread points vertically
        if (any(k)) {
            warning("-log10(p-value) infinite")
            y[k] <- stats::runif(sum(k), max(y[!k]), max(y[!k])*1.25)
        }

    if (geneID=="rownames") geneID <- rownames(deg) else geneID <- deg[,geneID]

        # plot background
        xRange <- c(-max(abs(x), na.rm=TRUE), max(abs(x), na.rm=TRUE))*1.2
        yRange <- c(0, max(y, na.rm=TRUE)*1.2)

        if (length(x) < 3e3) {
            graphics::plot(x, y, xlim=xRange, col="gray", cex=.5,
                           main=main, xlab=xlab, ylab=ylab, ...)
        } else {
            if (!packageExists("KernSmooth")) {
                graphics::plot(x, y, xlim=xRange, col="gray", cex=.5,
                               main=main, xlab=xlab, ylab=ylab, ...)
            } else {
                graphics::smoothScatter(x, y, xlim=xRange, nrpoints=0,
                                        main=main, xlab=xlab, ylab=ylab, ...)
            }
        }
        graphics::abline(h=0, lty=1)

        hline <- min(y[which(deg$adj.P.Val < padj)])
        if (length(hline)==0) hline <- 0
        graphics::abline(v=c(-lfc,lfc), h=hline, lty=2)
        # add features crossing threshollds
        ff <- which(abs(x) > lfc & deg$adj.P.Val < padj)
        if (length(ff) >= 1) {
            cc <- ifelse(x > 0, clCol[1], clCol[2])
            graphics::points(x[ff], y[ff], col=cc[ff], cex=.5)

            # label top-DEGs
            if (topN > 0) {
                ff2 <- seq_len(nrow(deg)) %in% ff & abs(deg$AveExpr) > ave
                dn <- rank(x[ff2]) <= topN
                up <- rank(-x[ff2]) <= topN
                top <- dn | up
                if (sum(top) > 0)
                graphics::text(x[ff2][top], y[ff2][top],geneID[ff2][top], cex=cexText)
            }
        }
    return(data.frame(x, y, geneID,
                        row.names=rownames(deg), stringsAsFactors=FALSE))
    }

    mtextFun <- function(...) graphics::mtext(...,side=3, cex=.67, line=0)

    if(class(deg)[1] == "list") {
        K <- length(deg)
        snk <- lapply(seq_len(K), function(k) {
            plotVolc(deg[[k]],
                    clCol=c(classCol[k %% length(classCol)], "gray"))
            mtextFun(paste(names(deg)[k],"up "), adj=1)
            mtextFun(paste(" ",names(deg)[k],"down"), adj=0)
            })
    } else {
        snk <- plotVolc(deg)
        mtextFun(paste(" ",attr(deg, "contrast")[1]), adj=1)
        mtextFun(paste(attr(deg, "contrast")[2], " "), adj=0)
    }
    invisible(snk)
}