R/MAplot.R

In systemPipeR/systemPipeTools: Tools for data visualization

#############
## MAplot ##
#############

#' @title MAplot
#' @description This function plots log2 fold changes (y-axis) versus the mean
#' of normalized counts (on the x-axis). Statistically significant features
#' are colored.

#' @param degseqDF object of class `data.frame` generated by
#' [systemPipeR::run_edgeR()] or [systemPipeR::run_DESeq2()].
#' @param FDR.cutoff filter cutoffs for the p-value adjusted.
#' @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 vecto`r 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
#' MAplot(degseqDF,
#'     comparison = "M12-A12", filter = c(Fold = 1, FDR = 20),
#'     genes = "ATCG00280"
#' )
#' @export
#' @importFrom DESeq2 results DESeq lfcShrink
#' @importFrom ggplot2 ggplot aes aes_string geom_point scale_colour_manual
#' @importFrom ggplot2 scale_x_continuous geom_smooth ggsave .data
#' @importFrom plotly ggplotly
#' @importFrom stats setNames
#' @keywords visualization
MAplot <- function(degseqDF, FDR.cutoff = 0.05, comparison,
                   filter = c(Fold = 2, FDR = 10),
                   genes = "NULL", plotly = FALSE, savePlot = FALSE,
                   filePlot = NULL) {
    padj <- NULL
    if (!is.data.frame(degseqDF)) {
        stop("'degseqDF' needs to be assignes an object of class 'data.frame'.
                                For more information check 'help(run_DESeq2)'")}
    table <- degseqDF %>%
        dplyr::select(
            paste0(comparison, "_baseMean"), paste0(comparison, "_logFC"),
            paste0(comparison, "_FDR")) %>%
        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, "_baseMean")]],
            y = .data[[paste0(comparison, "_logFC")]],
            label = names)) +
        ggplot2::geom_point(ggplot2::aes(colour = 
                .data[[paste0(comparison, "_FDR")]] < FDR.cutoff), size = 0.5) +
        ggplot2::scale_colour_manual(
            name = paste0("FDR < ", FDR.cutoff),
            values = stats::setNames(c("red", "grey"), c(TRUE, FALSE))) +
        ggplot2::scale_x_continuous(trans = "log10", limits = c(0.1, 300000)) +
        ggplot2::ggtitle(comparison) +
        ggrepel::geom_text_repel(data = genes)
    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))))
}