R/volcanoPlot.R

In prostarproteomics/DAPAR: Tools for the Differential Analysis of Proteins Abundance with R

#' @title Volcanoplot of the differential analysis
#' 
#' @description 
#' Plots a volcanoplot after the differential analysis.
#' Typically, the log of Fold Change is represented on the X-axis and the
#' log10 of the p-value is drawn on the Y-axis. When the \code{threshold_pVal}
#' and the \code{threshold_logFC} are set, two lines are drawn respectively on
#' the y-axis and the X-axis to visually distinguish between differential and
#' non differential data.
#'
#' @param logFC A vector of the log(fold change) values of the differential
#' analysis.
#' @param pVal A vector of the p-value values returned by the differential
#' analysis.
#' @param threshold_pVal A floating number which represents the p-value that
#' separates differential and non-differential data.
#' @param threshold_logFC A floating number which represents the log of the
#' Fold Change that separates differential and non-differential data.
#' @param conditions A list of the names of condition 1 and 2 used for the
#' differential analysis.
#' @param colors xxx
#' @return A volcanoplot
#' @author Florence Combes, Samuel Wieczorek
#' @examples
#' data(Exp1_R25_prot, package="DAPARdata")
#' obj <- Exp1_R25_prot[seq_len(100)]
#' level <- 'protein'
#' metacell.mask <- match.metacell(GetMetacell(obj), c("Missing POV", "Missing MEC"), level)
#' indices <- GetIndices_WholeMatrix(metacell.mask, op = ">=", th = 1)
#' obj <- MetaCellFiltering(obj, indices, cmd = "delete")
#' qData <- Biobase::exprs(obj$new)
#' sTab <- Biobase::pData(obj$new)
#' limma <- limmaCompleteTest(qData, sTab)
#' diffAnaVolcanoplot(limma$logFC[, 1], limma$P_Value[, 1])
#'
#' @export
#'
#'
diffAnaVolcanoplot <- function(logFC = NULL,
    pVal = NULL,
    threshold_pVal = 1e-60,
    threshold_logFC = 0,
    conditions = NULL,
    colors = NULL) {
    
   pkgs.require('graphics')
    
    xtitle <- paste("log2 ( mean(", conditions[2], ") / mean(", 
        conditions[1], ") )",
        sep = ""
    )

    if (is.null(colors)) {
        colors <- list(In = "orange", Out = "grey")
    }

    if (is.null(logFC) || is.null(pVal)) {
        p <- plot(-1, -1,
            xlab = xtitle,
            ylab = "- log10 ( p-value )",
            xlim = range(0, 1),
            xaxt = "n",
            yaxt = "n"
        )
        return(NULL)
    }

    x <- logFC
    y <- -log10(pVal)

    colorCode <- c(colors$Out, colors$In)
    color <- rep(colorCode[1], length(y))

    for (i in seq_len(length(y))) {
        if ((y[i] >= threshold_pVal) && (abs(x[i]) >= threshold_logFC)) {
            color[i] <- colorCode[2]
        }
    }

    p <- plot(x,
        y,
        xlab = xtitle,
        ylab = "- log10 ( p-value )",
        xlim = c(-max(abs(x)), max(abs(x))),
        ylim = c(0, max(y)),
        col = color,
        pch = 16,
        las = 1,
        cex = 1,
        cex.lab = 1.5,
        cex.axis = 1.5,
        cex.main = 3
    )

    graphics::abline(h = threshold_pVal, col = "gray")
    graphics::abline(v = threshold_logFC, col = "gray")
    graphics::abline(v = -threshold_logFC, col = "gray")

    return(p)
}



#' @title Volcanoplot of the differential analysis
#' 
#' @description 
#' #' Plots an interactive volcanoplot after the differential analysis.
#' Typically, the log of Fold Change is represented on the X-axis and the
#' log10 of the p-value is drawn on the Y-axis. When the \code{threshold_pVal}
#' and the \code{threshold_logFC} are set, two lines are drawn respectively on
#' the y-axis and the X-axis to visually distinguish between differential and
#' non differential data. With the use of the package Highcharter, a
#' customizable tooltip appears when the user put the mouse's pointer over
#' a point of the scatter plot.
#'
#'
#' @param df A dataframe which contains the following slots :
#' x : a vector of the log(fold change) values of the differential analysis,
#' y : a vector of the p-value values returned by the differential analysis.
#' index : a vector of the rowanmes of the data.
#' This dataframe must has been built with the option stringsAsFactors set
#' to FALSE. There may be additional slots which will be used to show
#' informations in the tooltip. The name of these slots must begin with the
#' prefix "tooltip_". It will be automatically removed in the plot.
#' @param threshold_pVal A floating number which represents the p-value that
#' separates differential and non-differential data.
#' @param threshold_logFC A floating number which represents the log of the
#' Fold Change that separates differential and non-differential data.
#' @param conditions A list of the names of condition 1 and 2 used for the
#' differential analysis.
#' @param clickFunction A string that contains a JavaScript function used to
#' show info from slots in df. The variable this.index refers to the slot
#' named index and allows to retrieve the right row to show in the tooltip.
#' @param pal xxx
#' @return An interactive volcanoplot
#' @author Samuel Wieczorek
#' @examples
#' library(highcharter)
#' data(Exp1_R25_prot, package="DAPARdata")
#' obj <- Exp1_R25_prot[seq_len(100)]
#' level <- 'protein'
#' metacell.mask <- match.metacell(GetMetacell(obj), c("Missing POV", "Missing MEC"), level)
#' indices <- GetIndices_WholeMatrix(metacell.mask, op = ">=", th = 1)
#' obj <- MetaCellFiltering(obj, indices, cmd = "delete")$new
#' qData <- Biobase::exprs(obj)
#' sTab <- Biobase::pData(obj)
#' data <- limmaCompleteTest(qData, sTab)
#' df <- data.frame(
#'     x = data$logFC, y = -log10(data$P_Value),
#'     index = as.character(rownames(obj))
#' )
#' colnames(df) <- c("x", "y", "index")
#' tooltipSlot <- c("Fasta_headers", "Sequence_length")
#' df <- cbind(df, Biobase::fData(obj)[, tooltipSlot])
#' colnames(df) <- gsub(".", "_", colnames(df), fixed = TRUE)
#' if (ncol(df) > 3) {
#'     colnames(df)[seq.int(from = 4, to = ncol(df))] <-
#'         paste("tooltip_", colnames(df)[seq.int(from = 4, to = ncol(df))],
#'          sep = "")
#' }
#' hc_clickFunction <- JS("function(event) {
#' Shiny.onInputChange('eventPointClicked',
#' [this.index]+'_'+ [this.series.name]);}")
#' cond <- c("25fmol", "10fmol")
#' diffAnaVolcanoplot_rCharts(df, 2.5, 1, cond, hc_clickFunction)
#'
#' @export
#'
#'
diffAnaVolcanoplot_rCharts <- function(df,
    threshold_pVal = 1e-60,
    threshold_logFC = 0,
    conditions = NULL,
    clickFunction = NULL,
    pal = NULL) {
    xtitle <- paste("log2 ( mean(", conditions[2], ") / mean(", 
        conditions[1], ") )", sep = "")

    if (is.null(clickFunction)) {
        clickFunction <-
            JS("function(event) {
                Shiny.onInputChange(
                'eventPointClicked', 
                [this.index]+'_'+ [this.series.name]);
                }")
    }

    if (is.null(pal)) {
        pal <- list(In = "orange", Out = "gray")
    } else {
        if (length(pal) != 2) {
            warning("The palette must be a list of two items: In and Out. 
                Set to default.")
            pal <- list(In = "orange", Out = "gray")
        }
    }

    df <- cbind(df,
        g = ifelse(df$y >= threshold_pVal & 
                abs(df$x) >= threshold_logFC, "g1", "g2")
    )


    i_tooltip <- which(startsWith(colnames(df), "tooltip"))
    txt_tooltip <- NULL
    for (i in i_tooltip) {
        t <- txt_tooltip <- paste(txt_tooltip, "<b>", gsub("tooltip_", "",
            colnames(df)[i],
            fixed = TRUE
        ),
        " </b>: {point.", colnames(df)[i], "} <br> ",
        sep = ""
        )
    }

    leftBorder <- data.frame(
        x = c(min(df$x), -threshold_logFC, -threshold_logFC),
        y = c(threshold_pVal, threshold_pVal, max(df$y))
    )
    rightBorder <- data.frame(
        x = c(max(df$x), threshold_logFC, threshold_logFC),
        y = c(threshold_pVal, threshold_pVal, max(df$y))
    )

    title <- NULL
    title <- paste0(conditions[1], "_vs_", conditions[2])

    h1 <- highchart() %>%
        hc_add_series(data = df, type = "scatter", hcaes(x, y, group = g)) %>%
        hc_colors(c(pal$In, pal$Out)) %>%
        my_hc_chart(zoomType = "xy", chartType = "scatter") %>%
        hc_legend(enabled = FALSE) %>%
        hc_title(
            text = title,
            margin = 20, align = "center",
            style = list(size = 20, color = "black", useHTML = TRUE)
        ) %>%
        hc_yAxis(title = list(text = "-log10(pValue)")) %>%
        hc_xAxis(
            title = list(text = "logFC"),
            plotLines = list(
                list(
                    color = "grey", 
                    width = 1, 
                    value = 0, 
                    zIndex = 5
                    )
                )
        ) %>%
        hc_tooltip(headerFormat = "", pointFormat = txt_tooltip) %>%
        hc_plotOptions(
            line = list(
                marker = list(enabled = FALSE),
                dashStyle = "Dash"
            ),
            series = list(
                animation = list(duration = 100),
                cursor = "pointer",
                point = list(events = list(
                    click = clickFunction
                ))
            )
        ) %>%
        my_hc_ExportMenu(filename = "volcanoplot") %>%
        hc_add_series(data = leftBorder, type = "line", color = "grey") %>%
        hc_add_series(data = rightBorder, type = "line", color = "grey")

    return(h1)
}