R/volcanoplot.R

In systemPipeR/systemPipeTools: Tools for data visualization

##################
## volcanoplot ##
##################
#' @title Volcano plot with `volcanoplot`
#' @description A simple function that shows statistical significance
#' (`p-value`) versus magnitude of change (`log2 fold change`).
#' @param degseqDF object of class `data.frame` generated by
#' [systemPipeR::run_edgeR()] or [systemPipeR::run_DESeq2()].
#' @param comparison `character vector` specifying the factor names for 
#' comparison.
#' @param filter Named vector with filter cutoffs of format c(Fold=2, FDR=1)
#' where Fold refers to the fold change cutoff (unlogged) and FDR to the 
#' p-value cutoff.
#' @param genes `character vector` of genes names to show on the plot.
#' @param plotly logical: when `FALSE` (default), the `ggplot2` plot will be 
#' returned.
#' `TRUE` option returns the `plotly` version of the plot.
#' @param savePlot logical: when `FALSE` (default), the plot will not be saved.
#' If `TRUE` the plot will be saved, and requires the `filePlot` argument.
#' @param filePlot file name where the plot will be saved. For more information,
#' please consult the [ggplot2::ggsave()] function.
#'
#' @return returns an object of `ggplot` or `plotly` class.
#'
#' @examples
#' ## Load targets file and count reads dataframe
#' targetspath <- system.file("extdata", "targets.txt", package = "systemPipeR")
#' targets <- read.delim(targetspath, comment = "#")
#' cmp <- systemPipeR::readComp(
#'     file = targetspath, format = "matrix",
#'     delim = "-"
#' )
#' countMatrixPath <- system.file("extdata", "countDFeByg.xls",
#'     package = "systemPipeR"
#' )
#' countMatrix <- read.delim(countMatrixPath, row.names = 1)
#' ### DEG analysis with `systemPipeR`
#' degseqDF <- systemPipeR::run_DESeq2(
#'     countDF = countMatrix,
#'     targets = targets, cmp = cmp[[1]], independent = FALSE)
#' DEG_list <- systemPipeR::filterDEGs(
#'     degDF = degseqDF,
#'     filter = c(Fold = 2, FDR = 10))
#' ## Plot
#' volcanoplot(degseqDF,
#'     comparison = "M12-A12", filter = c(Fold = 1, FDR = 20),
#'     genes = "ATCG00280")
#' @export
#' @importFrom dplyr select mutate filter
#' @importFrom ggplot2 ggplot aes aes_string geom_point geom_vline ggtitle xlab 
#' ylab theme element_text rel scale_color_manual theme_bw ggsave
#' @importFrom ggrepel geom_text_repel
#' @importFrom plotly ggplotly
#' @importFrom tibble rownames_to_column
volcanoplot <- function(degseqDF, comparison,
                        filter = c(Fold = 2, FDR = 10), genes = "NULL",
                        plotly = FALSE, savePlot = FALSE, filePlot = NULL) {
    if (!inherits(degseqDF, "data.frame")) {
          stop("degseqDF needs to be assigned an object of class 'data.frame'")}
    table <- degseqDF %>%
        dplyr::select(
            paste0(comparison, "_FDR"),
            paste0(comparison, "_logFC")
        ) %>%
        tibble::rownames_to_column(var = "names") %>%
    dplyr::mutate(significant = 
                    .data[[paste0(comparison, "_FDR")]] <= filter["FDR"] / 100 &
        abs(as.numeric(.data[[paste0(comparison, "_logFC")]])) > filter["Fold"])
    table$significant[is.na(table$significant)] <- FALSE
    if (!is.null(genes)) {
        genes <- table %>%
            dplyr::filter(names %in% genes)}
    plot <- ggplot2::ggplot(table, ggplot2::aes(
        x = .data[[paste0(comparison, "_logFC")]],
        y = -log10(as.numeric(.data[[paste0(comparison, "_FDR")]])),
        label = names)) +
        ggplot2::geom_point(ggplot2::aes_string(color = "significant")) +
        ggplot2::geom_vline(
            xintercept = c(-filter["Fold"], filter["Fold"]),
            linetype = 2) +
        ggplot2::ggtitle(comparison) +
        ggplot2::xlab("log2 fold change") +
        ggplot2::ylab("-log10(p-value)") +
        ggrepel::geom_text_repel(data = genes) +
        ggplot2::theme(
            legend.position = "none",
            plot.title = ggplot2::element_text(
                size = ggplot2::rel(1.5), hjust = 0.5),
            axis.title = ggplot2::element_text(size = ggplot2::rel(1.25)) ) +
        ggplot2::scale_color_manual(values = c("#524b4b", "#e11f28")) +
        ggplot2::theme_bw()
    if (savePlot == TRUE) {
        if (is.null(filePlot)) {
            stop("Argument 'filePlot' is missing, please provide file name.") }
        ggplot2::ggsave(plot = plot, filename = filePlot)}
    if (plotly == TRUE) {
        return(suppressWarnings(suppressMessages(plotly::ggplotly(plot))))}
    return(suppressWarnings(suppressMessages(print(plot))))
}