R/results_barplot.R

In NathanSkene/MAGMA_Celltyping: Find Causal Cell-Types Underlying Complex Trait Genetics

#' Plot barplot of enrichment results - plotting -log10 (p-value) but cut-off
#' significance is based on FDR where if either linear or top10% result is <
#' FDR threshold the result will be plotted
#' 
#' Plot of cell-type enrichment results from 
#' \link[MAGMA.Celltyping]{celltype_associations_pipeline}.
#' 
#' @param merged_results Enrichment results generated by 
#' \link[MAGMA.Celltyping]{celltype_associations_pipeline} and merged by 
#' \link[MAGMA.Celltyping]{merge_results}.
#' @param title Plot title.
#' @param x_lab Plot x-axis label.
#' @param y_lab Plot y-axis label.
#' @param fdr_thresh FDR filtering threshold.
#' @param x_var x-axis variable.
#' @param y_var y-axis variable.
#' @param fill_var Fill variable.
#' @param fill_cols The two colours to use for the fill param (Vector).
#' @param horz_line_p add dashed line for a p-value cut-off - default of 0.01
#' @param sort_celltypes_by Sort the order of the cell types in the plot by the
#' enrichment results. Choose to sort by Top 10\% or Linear results or set to
#' NULL to avoid sorting. Default is Top 10\%.
#' @param show_plot Whether to print the plot.
#' @param save_path Path to save plot to.
#' @param height Plot height.
#' @param width Plot width. 
#' @inheritParams ggplot2::ggsave
#' 
#' @return ggplot object.
#' 
#' @examples 
#' MAGMA_results <- MAGMA.Celltyping::enrichment_results
#' merged_results <- MAGMA.Celltyping::merge_results(MAGMA_results)
#' bar <- MAGMA.Celltyping::results_barplot(
#'     merged_results = merged_results,
#'     fdr_thresh = 1)  
#' @export
#' @importFrom ggplot2 ggplot element_blank ggsave geom_vline geom_bar aes_string labs theme_bw theme scale_fill_manual
results_barplot <- function(merged_results,
                            title = NULL,
                            y_lab = NULL,
                            x_lab = "-log10(p-value)",
                            fdr_thresh = .05,
                            y_var = "Celltype",
                            x_var = "-log10p",
                            fill_var = "EnrichmentMode",
                            fill_cols = c("#60A9BD","#3A62A0"),
                            horz_line_p = 0.01,
                            sort_celltypes_by = "Top 10%",
                            show_plot = TRUE,
                            height = 5, 
                            width = 7,
                            dpi = 300,
                            save_path = file.path(
                                tempdir(),
                                "MAGMA_Celltyping.barplot.jpg")
                            ) { 
    #### Avoid confusing checks ####
    FDR <- `-log10p` <- log10p <- NULL;
    #### Check args #####
    if (!"GWAS" %in% colnames(merged_results)) {
        merged_results$GWAS <- merged_results$dataset
    }
    if (is.null(fdr_thresh)) fdr_thresh <- 1
    subtitle <- if (!is.null(fdr_thresh)) paste0("FDR < ", fdr_thresh) else NULL
    #### Filter by FDR ####
    plot_dat <- subset(merged_results, FDR < fdr_thresh)
    messager(formatC(nrow(plot_dat), big.mark = ","),
             "results @ FDR <",fdr_thresh)
    if(nrow(plot_dat)==0) stop("Filtered data must contain >0 rows.")
    if(!sort_celltypes_by %in% unique(plot_dat$EnrichmentMode)){
      if(!is.null(sort_celltypes_by)){
        stop("sort_celltypes_by must be one of: ",
             paste(unique(plot_dat$EnrichmentMode),", "))
      }
    }
    #add -log10p-val col
    plot_dat[,`-log10p`:=-1*log10p]
    #sort by given results
    if(!is.null(sort_celltypes_by)){
      y_var_new <- paste0(y_var,"_fct")
      plot_dat[,eval(y_var_new):= factor(get(y_var),
                                     levels=plot_dat[get(fill_var)==sort_celltypes_by,][order(get(x_var),decreasing=FALSE),][,get(y_var)])]
      y_var <- y_var_new
    }
    #### Plot ####
    bar <- ggplot2::ggplot(
      data = plot_dat,
      ggplot2::aes_string(y = y_var,
                   x = x_var,
                   fill = fill_var)) + 
      ggplot2::geom_bar(stat = "identity", 
                        position = "dodge") +
      ggplot2::labs(x = x_lab, 
                    y = y_lab,
                    title = title,
                    legend = NULL) +
      ggplot2::theme_bw()+
      ggplot2::theme(legend.title=ggplot2::element_blank(),
                     legend.position="top",
                     plot.margin = ggplot2::margin(0, 0, 0, 0, "pt")
                     )+
      ggplot2::scale_fill_manual(values=fill_cols)+
      ggplot2::geom_vline(
        linetype="dotted",
        xintercept = -log(
          as.numeric(horz_line_p),
          10), 
        colour = "darkred", 
        alpha=.8)
    #### Print ####
    if (show_plot) print(bar)
    #### Save ####
    if (!is.null(save_path)){ 
        dir.create(dirname(save_path), showWarnings = FALSE, recursive = TRUE)
        ggplot2::ggsave(filename = save_path, 
                        plot = bar, 
                        dpi = dpi,
                        height = height, width = width)
    }
    return(bar)
}