R/plotBISEPOutput.R

In chapmandu2/CollateralVulnerability2016: Collateral Vulnerability Pipeline

#' BISEP Plot
#'
#' Plot to visualise BISEP output
#'
#' @param con A \code{SQLiteConnection} object
#' @param highlighted_gene Highlight the provided gene on the plot.  Default NULL - ie don't highlight anything
#' @param pi_value_th Set a PI value threshold.  Default 0.2
#' @param bisep_pval_th Set a BISEP p-value threshold.  Default 0.1
#' @param bi_value_th Set a Bimodal Index threshold.  Default 0
#' @param chromosome_sel A vector of chromosome id's.  Default is 0 which means all chromosomes.
#' @param xval The variable to plot on the x-axis.
#' @param yval The variable to plot on the y-axis
#' @return A ggplot2 object
#' @import ggplot2
#' @export
plotBISEPOutput <- function(con, highlighted_gene=NULL, pi_value_th=0.2, bisep_pval_th=0.1, bi_value_th=0, chromosome_sel=0, xval='pi_value', yval='bisep_pval', colval='selected') {

    #get the BISEP output data
    bisep_output <- dplyr::src_sqlite(con@dbname) %>%
        dplyr::tbl('combined_results') %>%
        dplyr::collect()

    plot_data <- bisep_output %>%
        dplyr::mutate(selected=pi_value <= pi_value_th & bisep_pval <= bisep_pval_th)

    if(bi_value_th > 0) {
        plot_data <- plot_data %>%
            dplyr::mutate(selected=selected & BI >= bi_value_th)
    }

    if(chromosome_sel != 0) {
        plot_data <- plot_data %>%
            dplyr::mutate(selected=selected & as.character(seqnames) %in% chromosome_sel)
    }

    num_filtered_genes <- plot_data %>% dplyr::filter(selected) %>% nrow()

    #internal function to work out what to use as an intercept
    intercept_val <- function(axis_var) {
        switch(axis_var,
               pi_value=pi_value_th,
               bisep_pval=bisep_pval_th,
               BI=bi_value_th)
    }

    #bisep_pval vs pi_val plot
    p1 <- ggplot(plot_data, aes_string(x=xval, y=yval, colour=colval)) +
        geom_point(size=rel(0.3)) +
        geom_vline(xintercept=intercept_val(xval), colour='blue') +
        geom_hline(yintercept=intercept_val(yval), colour='blue') +
        theme_bw() +
        ggtitle(sprintf("%s genes @ pi_value < %s & bisep_pval < %s & bi_value > %s", num_filtered_genes, pi_value_th, bisep_pval_th, bi_value_th))

    #handle the need for the bisep_pval to be on log scale
    if(xval == 'bisep_pval') {
        p1 <- p1 +  scale_x_log10()
    }

    if(yval == 'bisep_pval') {
        p1 <- p1 + scale_y_log10()
    }

    #highlight the selected gene if provided
    if(!is.null(highlighted_gene)) {
        highlight_data <- plot_data %>%
            dplyr::filter(gene_id == highlighted_gene)
        p1 <- p1 + geom_point(data=highlight_data, shape=1, colour='black', size=rel(2)) +
            ggrepel::geom_label_repel(data=highlight_data,
                                      aes(label = gene_name),
                                      size = 5, colour='black',
                                      box.padding = unit(0.35, "lines"),
                                      point.padding = unit(0.3, "lines"))
    }


    return(p1)

}