R/MOTIFBREAKR_plot.R

In RajLabMSSM/echoannot: echoverse module: Annotate fine-mapping results

#' Plot \pkg{motifbreakR} results
#'
#' Plot motif disruption results generated by \link[echoannot]{MOTIFBREAKR}.
#' 
#' \emph{Notes:}\cr
#' \itemize{
#' \item{Saving as a PDF seems to work much better than PNG format 
#' (at least when using \pkg{grDevices}).
#' }
#' }
#' @source
#' \href{https://pubmed.ncbi.nlm.nih.gov/26272984/}{Publication}
#' \href{https://github.com/Simon-Coetzee/MotifBreakR}{GitHub}
#' \href{https://github.com/Simon-Coetzee/motifBreakR/issues/31}{Viewport bug}
#' @param mb_filter \link[GenomicRanges]{GRanges} object 
#' used to filter \code{mb_res} before plotting. 
#' @inheritParams MOTIFBREAKR
#' @inheritParams MOTIFBREAKR_filter
#' @inheritParams grDevices::png
#' @inheritParams motifbreakR::plotMB
#' @returns Named list of motif plot paths.
#' 
#' @export
#' @importFrom grDevices png dev.off
#' @importFrom stats setNames
#' @importFrom echodata dt_to_granges
#' @importFrom BiocGenerics %in%
#' @examples
#' library(BSgenome) ## <-- IMPORTANT!
#' library(BSgenome.Hsapiens.UCSC.hg19) ## <-- IMPORTANT!
#' #### Example fine-mapping results ####
#' merged_DT <- echodata::get_Nalls2019_merged()
#' #### Run motif analyses ####
#' mb_res <- MOTIFBREAKR(rsid_list = c("rs11175620"),
#'                       # limit the number of datasets tested 
#'                       # for demonstration purposes only
#'                       pwmList_max = 5,
#'                       calculate_pvals = FALSE)  
#' \dontrun{
#' plot_paths <- MOTIFBREAKR_plot(mb_res = mb_res) 
#' }
MOTIFBREAKR_plot <- function(mb_res, 
                             mb_filter=NULL,
                             rsid=NULL,
                             effect=c("strong", "weak"),
                             results_dir=file.path(tempdir(),"results"),
                             height=3,
                             width=7,
                             verbose=TRUE){
    
    requireNamespace("motifbreakR")
    requireNamespace("MotifDb")
    requireNamespace("BSgenome.Hsapiens.UCSC.hg19")
    requireNamespace("grDevices")
    id <- NULL;
     
    mb_res <- echodata::dt_to_granges(dat = mb_res, 
                                      style = "UCSC", ## IMPORTANT!
                                      verbose = verbose)
    #### Make unique IDs ####
    mb_res <- MOTIFBREAKR_make_id(mb_res = mb_res) 
    #### Get SNPS ####
    if(is.null(rsid)) rsid <- mb_res$SNP_id
    rsid <- unique(rsid)
    if(length(rsid)>0){
        messager("Plotting",length(rsid),"unique RSID(s).",v=verbose)
    } else {
        stp <- "No valid RSIDs provided."
        stop(stp)
    }
    #### Filter SNPs based on another GRanges object ####
    if(!is.null(mb_filter)){
        mb_filter <- MOTIFBREAKR_make_id(mb_res = mb_filter)
        mb_res <- subset(mb_res, id %in% mb_filter$id)
    }
    #### Set up save path ####
    if(!is.null(results_dir)){
        dir.create(results_dir,showWarnings = FALSE, recursive = TRUE)
    } 
    # par(mfrow=c(2,2))
    plot_paths <- lapply(stats::setNames(rsid,rsid), 
                         function(rs){
        messager("Plotting motif disruption results:",rs,v=verbose)
        save_path <- file.path(results_dir,paste0(rs,'.pdf'))
        if(!is.null(results_dir)){ 
            grDevices::pdf(file = save_path,
                           width = width,
                           height = height)
        }
        motifbreakR::plotMB(results = mb_res,
                            rsid = rs,
                            effect = effect) 
        if(!is.null(results_dir)){
            grDevices::dev.off();    
        }
        return(save_path)
    }) 
    return(plot_paths)
}