R/differentialPlots.R

In ZarnackGroup/BindingSiteFinder: Binding site defintion based on iCLIP data

#' Diagnostic plots for the differential binding background
#'
#' To perform differential binding analysis between two conditions the
#' \code{\link{calculateBsBackground}} function groups crosslinks per gene
#' into those from binding sites and those from background regions.
#' The \code{\link{filterBsBackground}} function perfroms certain
#' filtering operations on that background to ensure that it's suitable
#' for differential testing. This function visually displays the effect
#' of these filtering operations.
#'
#'
#' @param object a \code{\link{BSFDataSet}} object with background counts filtered
#' by \code{\link{filterBsBackground}}
#' @param filter character; which filter to display in the plot (one of:
#' 'minCounts', 'balanceBackground', 'balanceCondition')
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{calculateBsBackground}}
#' @seealso \code{\link{filterBsBackground}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#'
#' # make testset
#' bds = makeBindingSites(bds, bsSize = 7)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = imputeBsDifferencesForTestdata(bds)
#' bds = calculateBsBackground(bds, anno.genes = gns, use.offset = FALSE)
#'
#' # use all filters and remove binding sites that fail (default settings)
#' bds = filterBsBackground(bds)
#'
#' # display minCount filter
#' plotBsBackgroundFilter(bds, filter = "minCounts")
#'
#' # display balance background filter
#' plotBsBackgroundFilter(bds, filter = "balanceBackground")
#'
#' # display balance condition filter
#' plotBsBackgroundFilter(bds, filter = "balanceCondition")
#'
#' @export
plotBsBackgroundFilter <- function(object, filter = c("minCounts", "balanceBackground" ,"balanceCondition")) {

    # bind locally used variables
    s <- r <- name <- value <- ratio.ref <- ratio.comp <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))
    if (is.null(object@params$filterBsBackground)) {
        msg0 = paste0("Global filter was not applied yet. Run filterBsBackground() to compute background first. \n")
        stop(msg0)
    }

    # handle display plot options
    filter = match.arg(filter, choices = c("minCounts", "balanceBackground", "balanceCondition"))

    # check for plot specific data
    if (is.null(object@plotData$filterBsBackground$data.minCounts) & filter == "minCounts") {
        msg0 = paste0("Plot selection: '", filter, "' is not possible.")
        msg1 = paste0("Make sure to run filterBsBackground with option '", filter, "' first.")
        stop(c(msg0,msg1))
    }
    if (is.null(object@plotData$filterBsBackground$data.balanceBackground) & filter == "balanceBackground") {
        msg0 = paste0("Plot selection: '", filter, "' is not possible.")
        msg1 = paste0("Make sure to run filterBsBackground with option '", filter, "' first.")
        stop(c(msg0,msg1))
    }
    if (is.null(object@plotData$filterBsBackground$data.balanceCondition) & filter == "balanceCondition") {
        msg0 = paste0("Plot selection: '", filter, "' is not possible.")
        msg1 = paste0("Make sure to run filterBsBackground with option '", filter, "' first.")
        stop(c(msg0,msg1))
    }


    # MAIN
    # --------------------------------------------------------------------------

    # make plot for minCounts filter
    if (filter == "minCounts") {
        # get options used
        optstr = object@params$filterBsBackground
        optstrNice = paste0("Cutoff=", optstr$minCounts.cutoff)

        # get data
        df = object@plotData$filterBsBackground$data.minCounts %>%
            mutate(r = rank(s))

        p = ggplot(df, aes(x = r, y = s)) +
            ggrastr::rasterise(ggpointdensity::geom_pointdensity(size = 2), dpi = 300) +
            theme_bw() +
            geom_hline(yintercept = optstr$minCounts.cutoff, linetype = "dashed") +
            scale_y_log10() +
            viridis::scale_color_viridis(option = "rocket") +
            theme(legend.position = "top") +
            guides(color = guide_colorbar(title.position = 'top', title.hjust = 0.5,
                                          barwidth = unit(20, 'lines'), barheight = unit(.5, 'lines'))) +
            annotate("text", x = max(df$r)/2, y = optstr$minCounts.cutoff,
                     label = paste0("Cutoff=", optstr$minCounts.cutoff),
                     vjust = -1.1) +
            labs(
                title = "filterBsBackground() - minCounts",
                subtitle = optstrNice,
                x = "Rank",
                y = "#Crosslinks per gene",
                color = "#Genes"
            )
    }

    # make plot for balanceBackground filter
    if (filter == "balanceBackground") {
        # get options used
        optstr = object@params$filterBsBackground
        optstrNice = paste0("Cutoff.bs=", optstr$balanceBackground.cutoff.bs,
                            " cutoff.bg=", optstr$balanceBackground.cutoff.bg)

        # get data
        df = object@plotData$filterBsBackground$data.balanceBackground %>%
            mutate(name = ifelse(name == "ratio.bg", "Background", "Binding site"))

        df = df[!is.na(df$value),]

        p = ggplot(df, aes(x = name, y = value)) +
            ggdist::stat_halfeye(aes(fill = name), adjust = 0.5, width = .6, .width = 0, justification = -.2, point_colour = NA) +
            geom_boxplot(aes(color = name), width = .12, outlier.size = 0.5) +
            geom_point(data = subset(df, name == "ratio.bg" & value < optstr$balanceBackground.cutoff.bg), aes(fill = name), shape = 21, alpha = 0.3) +
            geom_point(data = subset(df, name == "ratio.bs" & value > optstr$balanceBackground.cutoff.bs), aes(fill = name), shape = 21, alpha = 0.3) +
            annotate("rect", xmin = 1.75, xmax = 2.25, ymin = optstr$balanceBackground.cutoff.bs, ymax = 1, alpha = .3, fill = "#435B66") +
            annotate("rect", xmin = 0.75, xmax = 1.25, ymin = 0, ymax = optstr$balanceBackground.cutoff.bg, alpha = .3, fill = "#A76F6F") +
            annotate("text", x = 2.25, y = 0.5, vjust = -1, label = paste0("Remove genes with bs ratio above: ", optstr$balanceBackground.cutoff.bs)) +
            annotate("text", x = 1.25, y = 0.5, vjust = -1, label = paste0("Remove genes with bg ratio below: ", optstr$balanceBackground.cutoff.bg)) +
            coord_flip() +
            theme_bw() +
            scale_fill_manual(values = c("#A76F6F", "#435B66")) +
            scale_color_manual(values = c("#A76F6F", "#435B66")) +
            theme(legend.position = "none") +
            labs(
                title = "filterBsBackground() - balanceBackground",
                subtitle = optstrNice,
                x = NULL,
                y = "Background to binding site count ratio"
            )
    }

    # make plot for balanceCondition filter
    if (filter == "balanceCondition") {
        # get options used
        optstr = object@params$filterBsBackground
        optstrNice = paste0("Cutoff=", optstr$balanceCondition.cutoff)

        # get reference condition
        this.meta = getMeta(object)
        this.condition.reference = levels(this.meta$condition)[1]
        this.condition.comp = levels(this.meta$condition)[2]

        # get data
        df = object@plotData$filterBsBackground$data.balanceCondition %>% as.data.frame()
        df$ratio.ref = df[,3] / df[,2]
        df$ratio.comp = df[,4] / df[,2]
        df = df %>% select(ratio.ref, ratio.comp) %>% pivot_longer(everything())

        p = ggplot(df, aes(x = name, y = value)) +
            ggdist::stat_halfeye(aes(fill = name), adjust = 0.5, width = .6, .width = 0, justification = -.2, point_colour = NA) +
            geom_boxplot(aes(color = name), width = .12, outlier.size = 0.5) +
            geom_point(data = subset(df, name == "ratio.ref" & value < optstr$balanceCondition.cutoff), aes(fill = name), shape = 21, alpha = 0.3) +
            geom_point(data = subset(df, name == "ratio.comp" & value < optstr$balanceCondition.cutoff), aes(fill = name), shape = 21, alpha = 0.3) +
            annotate("rect", xmin = 0.05, xmax = 2.95, ymin = 0, ymax = optstr$balanceCondition.cutoff, alpha = .3, fill = "darkgrey") +
            coord_flip() +
            theme_bw() +
            scale_fill_manual(values = c("#A76F6F", "#435B66")) +
            scale_color_manual(values = c("#A76F6F", "#435B66")) +
            theme(legend.position = "none") +
            labs(
                title = "filterBsBackground() - balanceCondition",
                subtitle = optstrNice,
                x = NULL,
                y = "Condition count ratio"
            ) +
            scale_x_discrete(labels = c(paste0(this.condition.comp, "\n (comparison)"),
                                        paste0(this.condition.reference, "\n (reference)")))


    }

    return(p)
}


#' MA style plot
#'
#' Wrapper that plots differential binding results as MA plot. For each binding
#' site the estimated baseMean (log2) is shown on X and the fold-change (log2)
#' is shown on Y.
#'
#' @param object a \code{\link{BSFDataSet}} object with results calculated by
#' \code{\link{calculateBsFoldChange}}
#' @param what character; whether to show results for binding sites or the
#' background (one of: 'bs', 'bg')
#' @param sig.threshold numeric; what P value significance level to use
#' (default = 0.05)
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{calculateBsFoldChange}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#'
#' # make testset
#' bds = makeBindingSites(bds, bsSize = 7)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = imputeBsDifferencesForTestdata(bds)
#' bds = calculateBsBackground(bds, anno.genes = gns, use.offset = FALSE)
#'
#' # use all filters and remove binding sites that fail (default settings)
#' bds = filterBsBackground(bds)
#'
#' # calculate fold-changes
#' bds = calculateBsFoldChange(bds)
#'
#' # make MA plot
#' plotBsMA(bds)
#'
#' @export
plotBsMA <- function(object,
                     what = c("bs", "bg"),
                     sig.threshold = 0.05) {

    # bind locally used variables
    bs.padj <- bs.log2FoldChange <- sig <- bs.baseMean <- bg.padj <- NULL
    bg.log2FoldChange <- bg.baseMean <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))
    if (is.null(object@params$calculateBsFoldChange)) {
        msg0 = paste0("Fold-changes were not calculated yet. Run calculateBsFoldChange() first. \n")
        stop(msg0)
    }

    # handle display plot options
    what = match.arg(what, choices = c("bs", "bg"))

    this.ranges = getRanges(object)
    expected.cols = c("bs.padj", "bs.log2FoldChange", "bs.baseMean",
                      "bg.padj", "bg.log2FoldChange", "bg.baseMean")
    if (!all(expected.cols %in% colnames(as.data.frame(mcols(this.ranges))))) {
        msg0 = paste0("MA plot not possible, results columns: ", paste(expected.cols, collapse = ','), " not present.\n")
        msg1 = paste0("Make sure to run calculateBsFoldChange() first.\n")
        stop(c(msg0, msg1))
    }

    # MAIN
    # --------------------------------------------------------------------------

    optstr = object@params$calculateBsFoldChange
    optstrNice = paste0("alpha=", optstr$alpha)

    bright_up_down_not = c("#999999", "#68b1a5", "#874C62")
    dark_up_down_not = c("#4d4d4d", "#2b544d", "#623747")

    if (what == "bs") {
        df = as.data.frame(mcols(this.ranges)) %>%
            mutate(sig = ifelse(bs.padj < sig.threshold & bs.log2FoldChange > 0, "Up",
                                ifelse(bs.padj < sig.threshold & bs.log2FoldChange < 0, "Down", "Not"))) %>%
            mutate(sig = factor(sig, levels = c("Not", "Up", "Down"))) %>%
            arrange(sig)

        p = ggplot(df, aes(x = log2(bs.baseMean), y = bs.log2FoldChange, color = sig, fill = sig)) +
            ggrastr::rasterise(geom_point(shape = 21, stroke = 0.5, size = 1.5), dpi = 300) +
            geom_hline(yintercept = 0, color = "black", alpha = .5) +
            theme_bw() +
            scale_fill_manual(values = bright_up_down_not) +
            scale_color_manual(values = dark_up_down_not) +
            guides(color = guide_legend(override.aes = list(size = 4))) +
            theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
            labs(
                title = "plotBsMA() - binding sites",
                subtitle = optstrNice,
                x = "Base mean",
                y = "Fold-change (log2)",
                color = "Regulation",
                fill = "Regulation")
    }
    if (what == "bg") {
        df = as.data.frame(mcols(this.ranges)) %>%
            mutate(sig = ifelse(bg.padj < sig.threshold & bg.log2FoldChange > 0, "Up",
                                ifelse(bg.padj < sig.threshold & bg.log2FoldChange < 0, "Down", "Not"))) %>%
            mutate(sig = factor(sig, levels = c("Not", "Up", "Down"))) %>%
            arrange(sig)

        p = ggplot(df, aes(x = log2(bg.baseMean), y = bg.log2FoldChange, color = sig, fill = sig)) +
            ggrastr::rasterise(geom_point(shape = 21, stroke = 0.5, size = 1.5), dpi = 300) +
            geom_hline(yintercept = 0, color = "black", alpha = .5) +
            theme_bw() +
            scale_fill_manual(values = bright_up_down_not) +
            scale_color_manual(values = dark_up_down_not) +
            guides(color = guide_legend(override.aes = list(size = 4))) +
            theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
            labs(
                title = "plotBsMA() - background",
                subtitle = optstrNice,
                x = "Base mean",
                y = "Fold-change (log2)",
                color = "Regulation",
                fill = "Regulation")
    }


    return(p)
}

#' Volcano style plot
#'
#' Wrapper that plots differential binding results as volcano plot. For each binding
#' site the estimated fold-change (log2) is shown on X and the adjusted
#' P value (-log10) is shown on Y.
#'
#' @param object a \code{\link{BSFDataSet}} object with results calculated by
#' \code{\link{calculateBsFoldChange}}
#' @param what character; whether to show results for binding sites or the
#' background (one of: 'bs', 'bg')
#' @param sig.threshold numeric; what P value significance level to use
#' (default = 0.05)
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{calculateBsFoldChange}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#'
#' # make testset
#' bds = makeBindingSites(bds, bsSize = 7)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = imputeBsDifferencesForTestdata(bds)
#' bds = calculateBsBackground(bds, anno.genes = gns, use.offset = FALSE)
#'
#' # use all filters and remove binding sites that fail (default settings)
#' bds = filterBsBackground(bds)
#'
#' # calculate fold-changes
#' bds = calculateBsFoldChange(bds)
#'
#' # make volcano plot
#' plotBsVolcano(bds)
#'
#' @export
plotBsVolcano <- function(object,
                          what = c("bs", "bg"),
                          sig.threshold = 0.05) {

    # bind locally used variables
    bs.padj <- bs.log2FoldChange <- sig <- bg.padj <- bg.log2FoldChange <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))
    if (is.null(object@params$calculateBsFoldChange)) {
        msg0 = paste0("Fold-changes were not calculated yet. Run calculateBsFoldChange() first. \n")
        stop(msg0)
    }

    # handle display plot options
    what = match.arg(what, choices = c("bs", "bg"))

    this.ranges = getRanges(object)
    expected.cols = c("bs.padj", "bs.log2FoldChange", "bs.baseMean",
                      "bg.padj", "bg.log2FoldChange", "bg.baseMean")
    if (!all(expected.cols %in% colnames(as.data.frame(mcols(this.ranges))))) {
        msg0 = paste0("MA plot not possible, results columns: ", paste(expected.cols, collapse = ','), " not present.\n")
        msg1 = paste0("Make sure to run calculateBsFoldChange() first.\n")
        stop(c(msg0, msg1))
    }

    # MAIN
    # --------------------------------------------------------------------------

    optstr = object@params$calculateBsFoldChange
    optstrNice = paste0("alpha=", optstr$alpha)

    bright_up_down_not = c("#999999", "#68b1a5", "#874C62")
    dark_up_down_not = c("#4d4d4d", "#2b544d", "#623747")

    if (what == "bs") {
        df = as.data.frame(mcols(this.ranges)) %>%
            mutate(sig = ifelse(bs.padj < sig.threshold & bs.log2FoldChange > 0, "Up",
                                ifelse(bs.padj < sig.threshold & bs.log2FoldChange < 0, "Down", "Not"))) %>%
            mutate(sig = factor(sig, levels = c("Not", "Up", "Down"))) %>%
            arrange(sig)

        p = ggplot(df, aes(x = bs.log2FoldChange, y = -log10(bs.padj), color = sig, fill = sig)) +
            ggrastr::rasterise(geom_point(shape = 21, stroke = 0.5, size = 1.5), dpi = 300) +
            geom_vline(xintercept = 0, color = "black", alpha = .5) +
            theme_bw() +
            scale_fill_manual(values = bright_up_down_not) +
            scale_color_manual(values = dark_up_down_not) +
            guides(color = guide_legend(override.aes = list(size = 4))) +
            theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
            labs(
                title = "plotBsVolcano() - binding sites",
                subtitle = optstrNice,
                x = "Fold-change (log2)",
                y = "Adjusted P value (-log10)",
                color = "Regulation",
                fill = "Regulation")
    }
    if (what == "bg") {
        df = as.data.frame(mcols(this.ranges)) %>%
            mutate(sig = ifelse(bg.padj < sig.threshold & bg.log2FoldChange > 0, "Up",
                                ifelse(bg.padj < sig.threshold & bg.log2FoldChange < 0, "Down", "Not"))) %>%
            mutate(sig = factor(sig, levels = c("Not", "Up", "Down"))) %>%
            arrange(sig)

        p = ggplot(df, aes(x = bg.log2FoldChange, y = -log10(bg.padj), color = sig, fill = sig)) +
            ggrastr::rasterise(geom_point(shape = 21, stroke = 0.5, size = 1.5), dpi = 300) +
            geom_vline(xintercept = 0, color = "black", alpha = .5) +
            theme_bw() +
            scale_fill_manual(values = bright_up_down_not) +
            scale_color_manual(values = dark_up_down_not) +
            guides(color = guide_legend(override.aes = list(size = 4))) +
            theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
            labs(
                title = "plotBsVolcano() - background",
                subtitle = optstrNice,
                x = "Fold-change (log2)",
                y = "Adjusted P value (-log10)",
                color = "Regulation",
                fill = "Regulation")
    }

    return(p)
}


#' Gene Regulation Plot
#'
#' Display the fold-change of all binding sites from a given gene on a relative
#' per-nucleotide scale. Binding sites are displayed as dots and with increasing
#' log2 fold-change, they deviate stronger from the center line.
#'
#' For this function to work, binding sites must be assigned to hosting genes
#' using \code{\link{assignToGenes}}. It is also recommended to assing
#' binding sites to transcript regions with \code{\link{assignToTranscriptRegions}}.
#'
#' It is also necessary to calculate the log2 fold-change of binding sites between
#' two conditions using the differential binding workflow \code{\link{calculateBsFoldChange}}.
#'
#' If in addition the transcript regions of the binding sites are given, then
#' shapes are changed accordingly. An edge case can arise from the merging of two
#' \code{\link{BSFDataSet}} objects. If binding sites are overlapping and
#' slightly offset close to the end of a particular transcript region annotation,
#' they might be assigned to different regions in both objects. This results in
#' some ambiguity after the merge, where for instance a binding site can be
#' assigned to CDS and 3'UTR. To handle how such edge cases are displayed, the
#' \code{transcript.regions.outlier.handle} exists. As default, simply the
#' region of the object that was merged first is shown. If one is interested in
#' showing all regions, then the options \code{both} displays both annotations
#' at the same time and labels them accordingly.
#'
#' @param object object; a \code{\link{BSFDataSet}} object
#' @param plot.geneID character;  the gene id of the gene to display. The id must
#' match with the gene ids given in the annotation object.
#' @param anno.annoDB object; an object of class \code{OrganismDbi} that contains
#' the gene annotation (!!! Experimental !!!).
#' @param anno.genes object; an object of class \code{\link{GenomicRanges}}
#' that represents the gene ranges directly.
#' @param match.geneID character; meta column name of the gene ID
#' @param match.geneName character; meta column name of the gene name
#' @param plot.gene.n.tiles numeric; number of tiles the gene should be split in
#' @param alpha numeric; the alpha value to show significantly regulated
#' binding sites. This should match the alpha value used in \code{\link{calculateBsFoldChange}}.
#' @param lfc.cutoff numeric; log2 fold-change cutoff to show significantly
#' regulated binding sites. This should match the lfc.cutoff value used in
#' \code{\link{calculateBsFoldChange}}.
#' @param transcript.regions.outlier.handle character; the option how to handle
#' multiple transcript region annotations being present for the same binding
#' site.
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{ggplot2}
#'
#' @seealso \code{\link{BSFind}},
#' \code{\link{calculateBsFoldChange}}
#' \code{\link{assignToGenes}}
#' \code{\link{assignToTranscriptRegions}}
#'
#' @import ggplot2
#'
#' @examples
#'
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[2])
#'
#' # make testset
#' bds = makeBindingSites(bds, bsSize = 7)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = assignToTranscriptRegions(object = bds, anno.transcriptRegionList = regions)
#' bds = imputeBsDifferencesForTestdata(bds)
#' bds = calculateBsBackground(bds, anno.genes = gns, use.offset = FALSE)
#'
#' # use all filters and remove binding sites that fail (default settings)
#' bds = filterBsBackground(bds)
#'
#' # calculate fold-changes
#' bds = calculateBsFoldChange(bds)
#'
#' # make example plot
#' exampleGeneId = "ENSG00000253352.10"
#' geneRegulationPlot(bds, plot.geneID = exampleGeneId, anno.genes = gns)
#'
#' @export
geneRegulationPlot <- function(object,
                               plot.geneID = NULL,
                               anno.annoDB = NULL,
                               anno.genes = NULL,
                               match.geneID = "gene_id",
                               match.geneName = "gene_name",
                               plot.gene.n.tiles = 100,
                               alpha = 0.05,
                               lfc.cutoff = 2,
                               transcript.regions.outlier.handle = c("first", "second", "both", "remove"),
                               quiet = FALSE) {

    # bind locally used variables
    bs.padj <- bs.log2FoldChange <- sig <- sigDir <- plot.position <- transcriptRegion <-  NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    # type checks
    stopifnot(is(object, "BSFDataSet"))
    stopifnot(is.logical(quiet))

    # check presence of fold-changes -> must be present
    if (is.null(object@params$calculateBsFoldChange)) {
        msg0 = paste0("Fold-changes were not calculated yet. Run calculateBsFoldChange() first. \n")
        stop(msg0)
    }
    # check presence of gene assignment -> must be present
    if (is.null(getRanges(object)$geneID)) {
        msg0 = paste0("Gene assignment was not calculated. Run assignToGenes() first. \n")
        stop(msg0)
        # TODO  ideally this should test for the param present in object@params$assignToGenes,
        #       but this gets overwritten by the mergeing of two bds objects. Change this in the future.
    }
    # check presence of transcript regions -> can be missing
    if (is.null(getRanges(object)$transcriptRegion)) {
        msg0 = paste0("Transcript region assignment was not calculated. Run assignToTranscriptRegions() first to link binding sites to transcript regions. \n")
        warning(msg0)
        # TODO  ideally this should test for the param present in object@params$assignToTranscriptRegions,
        #       but this gets overwritten by the mergeing of two bds objects. Change this in the future.
        # TODO make a version of the plot that does not encode the transcript regions as shape information
    }

    # check presence of correct annotation resource -> one of both must be present
    # Check if none is specified
    if (is.null(anno.annoDB) & is.null(anno.genes)) {
        msg = paste0("None of the required annotation sources anno.annoDB or anno.genes was specified. ")
        stop(msg)
    }
    # Check if both are specified
    if (!is.null(anno.annoDB) & !is.null(anno.genes)) {
        msg = paste0("Both of the required annotation sources anno.annoDB or anno.genes are specified. Please provide only one of the two. ")
        stop(msg)
    }
    # Checks if anno.annoDB should be used
    if (!is.null(anno.annoDB) & is.null(anno.genes)) {
        stopifnot(is(anno.annoDB, "OrganismDb"))
        if (!is.null(anno.genes)) {
            msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
            stop(msg)
        } else {
            datasource = "anno.annoDB"
            # extract relevant annotation
            anno.genes = genes(anno.annoDB, columns=c("ENSEMBL", "SYMBOL", "GENEID"))
            # Create matching vectors for columns from input annotation
            # --------------------------------------------------------------------------
            selectID = as.character(anno.genes$GENEID)
            selectName = as.character(anno.genes$SYMBOL)
        }
    }
    # Checks if anno.genes should be used
    if (is.null(anno.annoDB) & !is.null(anno.genes)) {
        stopifnot(is(anno.genes, "GenomicRanges"))
        if (!is.null(anno.annoDB)) {
            msg = paste0("Parameter anno.annoDB and anno.genes are specified at the same time. Use only one of them.")
            stop(msg)
        } else {
            datasource = "anno.genes"
            # extract relevant annotation
            # check for duplicated annotation columns
            if (sum(colnames(mcols(anno.genes)) == match.geneID) > 1) {
                msg = paste0("The names of multiple columns of the annotation match the match.geneID parameter. Please use a unique column name for match.geneID in anno.genes.")
                stop(msg)
            }
            if (sum(colnames(mcols(anno.genes)) == match.geneName) > 1) {
                msg = paste0("The names of multiple columns of the annotation match the match.geneName parameter. Please use a unique column name for match.geneName in anno.genes.")
                stop(msg)
            }
            # check correct annotation columns
            inNames = c(match.geneID, match.geneName)
            annoColNames = colnames(mcols(anno.genes))
            presentNames = inNames[inNames %in% annoColNames]

            # extract columns from annotation based on available names
            present.cols = lapply(presentNames, function(x){
                match.col = mcols(anno.genes)[match(x, colnames(mcols(anno.genes)))][[1]]
                return(match.col)
            })
            names(present.cols) = presentNames

            # extract values from annotation for all present columns
            if (match.geneID %in% names(present.cols)) {
                selectID = present.cols[[match(match.geneID, names(present.cols))]]
            } else {
                msg = paste0("No meta column for ", match.geneID, " present. Creating custom ID.\n")
                if (!quiet) message(msg)
                selectID = paste0("CUSTOM", seq_along(anno.genes))
            }
            if (match.geneName %in% names(present.cols)) {
                selectName = present.cols[[match(match.geneName, names(present.cols))]]
            } else {
                msg = paste0("No meta column for ", match.geneName, " present.\n")
                if (!quiet) message(msg)
            }
        }
    }

    # check if gene to plot is in bds object -> must be present
    if (!(plot.geneID %in% getRanges(object)$geneID)) {
        msg0 = stop(paste0("Selected geneID: ", plot.geneID, " is not present in the meta columns of the given bds object. \n"))
        stop(msg0)
    }
    # check if gene to plot is in annotation object -> must be present
    if (!(plot.geneID %in% selectID)) {
        msg0 = stop(paste0("Selected geneID: ", plot.geneID, " is not present in the meta columns of the given gene annotation. \n"))
        stop(msg0)
    }

    transcript.regions.outlier.handle = match.arg(transcript.regions.outlier.handle, choices = c("first", "second", "both", "remove"))

    # --------------------------------------------------------------------------
    # MAIN
    # --------------------------------------------------------------------------

    # get ranges of binding sites for selected gene
    plot.bs = getRanges(object)
    plot.bs = plot.bs[plot.bs$geneID == plot.geneID]
    plot.bs = sort(plot.bs)
    plot.bs = resize(plot.bs, width = 1, fix = "center")

    # get gene ranges for selected gene
    plot.gene = anno.genes[anno.genes$gene_id == plot.geneID]
    plot.gene.name = plot.gene$gene_name

    # claculate relative position of binding sites in gene
    plot.gene.tiles = unlist(tile(plot.gene, n = plot.gene.n.tiles))
    plot.gene.tiles$position = 1:length(plot.gene.tiles)

    ols = findOverlaps(plot.bs, plot.gene.tiles)

    gene.strand = unique(as.character(strand(plot.gene)))
    if(gene.strand == "-") {
        position = plot.gene.n.tiles-(subjectHits(ols))
    } else {
        position = subjectHits(ols)
    }
    plot.bs$plot.position = position

    # --------------------------------------------------------------------------
    # PLOTS
    # --------------------------------------------------------------------------

    # make plots without transcript region annotation
    if (is.null(getRanges(object)$transcriptRegion)) {
        df.plot = mcols(plot.bs) %>% as.data.frame() %>%
            mutate(sig = (ifelse(bs.padj < alpha & abs(bs.log2FoldChange) > log2(lfc.cutoff), TRUE, FALSE))) %>%
            mutate(dir = factor(ifelse(bs.log2FoldChange > 0, "Up", "Down"), levels = c("Up", "Down"))) %>%
            mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
            mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down")))

        p = ggplot(df.plot, aes(x = plot.position, y = bs.log2FoldChange, color = sigDir, fill = sigDir)) +
            geom_hline(yintercept = 0, color = "#4d4d4d") +
            geom_segment(aes(x=plot.position, xend=plot.position, y=0, yend=bs.log2FoldChange), size = 1, color = "#4d4d4d") +
            geom_point(size = 4, stroke = 1, shape = 21) +
            scale_fill_manual(values = c("Not" = "#999999", "Up" = "#68b1a5", "Down" = "#874C62")) +
            scale_color_manual(values = c("Not" = "#4d4d4d", "Up" = "#2b544d", "Down" = "#623747")) +
            theme_minimal() +
            guides(color = guide_legend(override.aes = list(size = 4))) +
            theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
            ylim(-max(abs(df.plot$bs.log2FoldChange)), max(abs(df.plot$bs.log2FoldChange))) +
            xlim(0,plot.gene.n.tiles) +
            labs(
                title = paste0(plot.gene.name, " (", plot.geneID, ")"),
                x = "Relative position on the gene",
                y = "Fold-change (log2)",
                fill = "Regulation",
                color = "Regulation"
            )
        return(p)
    }

    # make plots with transcript region annotation
    if (!is.null(getRanges(object)$transcriptRegion)) {

        expected.regions = c("INTRON", "CDS", "UTR3", "UTR5", "OTHER")
        all.regions.present = unique(plot.bs$transcriptRegion)

        # plot cases where no region duplication is present
        if (all(all.regions.present %in% expected.regions)) {
            df.plot = mcols(plot.bs) %>% as.data.frame() %>%
                mutate(sig = (ifelse(bs.padj < alpha & abs(bs.log2FoldChange) > log2(lfc.cutoff), TRUE, FALSE))) %>%
                mutate(dir = factor(ifelse(bs.log2FoldChange > 0, "Up", "Down"), levels = c("Up", "Down"))) %>%
                mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
                mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
                mutate(duplicated = "no")

            p = ggplot(df.plot, aes(x = plot.position, y = bs.log2FoldChange, color = sigDir, fill = sigDir, shape = transcriptRegion)) +
                geom_hline(yintercept = 0, color = "#4d4d4d") +
                geom_segment(aes(x=plot.position, xend=plot.position, y=0, yend=bs.log2FoldChange), size = 1, color = "#4d4d4d") +
                geom_point(size = 4, stroke = 1) +
                scale_shape_manual(values=c("INTRON" = 23, "CDS" = 22, "UTR3" = 25, "UTR5" = 24, "OTHER" = 21)) +
                scale_fill_manual(values = c("Not" = "#999999", "Up" = "#68b1a5", "Down" = "#874C62")) +
                scale_color_manual(values = c("Not" = "#4d4d4d", "Up" = "#2b544d", "Down" = "#623747")) +
                theme_minimal() +
                guides(color = guide_legend(override.aes = list(size = 4))) +
                theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
                ylim(-max(abs(df.plot$bs.log2FoldChange)), max(abs(df.plot$bs.log2FoldChange))) +
                xlim(0,plot.gene.n.tiles) +
                labs(
                    title = paste0(plot.gene.name, " (", plot.geneID, ")"),
                    x = "Relative position on the gene",
                    y = "Fold-change (log2)",
                    fill = "Regulation",
                    color = "Regulation",
                    shape = "Transcript Region"
                )
            return(p)
        }

        # plot cases where region duplication is present
        if (any(all.regions.present %in% expected.regions)) {
            msg0 = paste0("At least one additional transcript region type detected compared to the currently supported expected regions: ",
                          paste(expected.regions, collapse = ", "),
                          ". Using transcript.regions.outlier.handle=", transcript.regions.outlier.handle," to solve. \n")
            if (!quiet) {
                warning(msg0)
            }

            if (transcript.regions.outlier.handle == "remove") {
                df.plot = mcols(plot.bs) %>% as.data.frame() %>%
                    filter(transcriptRegion %in% expected.regions) %>%
                    mutate(sig = (ifelse(bs.padj < alpha & abs(bs.log2FoldChange) > log2(lfc.cutoff), TRUE, FALSE))) %>%
                    mutate(dir = factor(ifelse(bs.log2FoldChange > 0, "Up", "Down"), levels = c("Up", "Down"))) %>%
                    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
                    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
                    mutate(duplicated = "no")

                p = ggplot(df.plot, aes(x = plot.position, y = bs.log2FoldChange, color = sigDir, fill = sigDir, shape = transcriptRegion)) +
                    geom_hline(yintercept = 0, color = "#4d4d4d") +
                    geom_segment(aes(x=plot.position, xend=plot.position, y=0, yend=bs.log2FoldChange), size = 1, color = "#4d4d4d") +
                    geom_point(size = 4, stroke = 1) +
                    scale_shape_manual(values=c("INTRON" = 23, "CDS" = 22, "UTR3" = 25, "UTR5" = 24, "OTHER" = 21)) +
                    scale_fill_manual(values = c("Not" = "#999999", "Up" = "#68b1a5", "Down" = "#874C62")) +
                    scale_color_manual(values = c("Not" = "#4d4d4d", "Up" = "#2b544d", "Down" = "#623747")) +
                    theme_minimal() +
                    guides(color = guide_legend(override.aes = list(size = 4))) +
                    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
                    ylim(-max(abs(df.plot$bs.log2FoldChange)), max(abs(df.plot$bs.log2FoldChange))) +
                    xlim(0,plot.gene.n.tiles) +
                    labs(
                        title = paste0(plot.gene.name, " (", plot.geneID, ")"),
                        x = "Relative position on the gene",
                        y = "Fold-change (log2)",
                        fill = "Regulation",
                        color = "Regulation",
                        shape = "Transcript Region",
                        caption = paste0("Not all binding sites might be shown due to: transcript.regions.outlier.handle=", transcript.regions.outlier.handle)
                    )
                return(p)
            }

            if (transcript.regions.outlier.handle == "first") {
                df.duplicated = mcols(plot.bs) %>% as.data.frame() %>% filter(grepl(",", transcriptRegion)) %>% mutate(duplicated = "yes") %>%
                    mutate(transcriptRegion = sapply(strsplit(transcriptRegion,","), `[`, 1))
                df.single = mcols(plot.bs) %>% as.data.frame() %>% filter(!grepl(",", transcriptRegion)) %>% mutate(duplicated = "no")
                df.plot = rbind.data.frame(df.duplicated, df.single) %>%
                    mutate(sig = (ifelse(bs.padj < alpha & abs(bs.log2FoldChange) > log2(lfc.cutoff), TRUE, FALSE))) %>%
                    mutate(dir = factor(ifelse(bs.log2FoldChange > 0, "Up", "Down"), levels = c("Up", "Down"))) %>%
                    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
                    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
                    mutate(transcriptRegion = stringr::str_trim(transcriptRegion, "left")) %>%
                    mutate(transcriptRegion = stringr::str_trim(transcriptRegion, "right"))

                p = ggplot(df.plot, aes(x = plot.position, y = bs.log2FoldChange, color = sigDir, fill = sigDir, shape = transcriptRegion)) +
                    geom_hline(yintercept = 0, color = "#4d4d4d") +
                    geom_segment(aes(x=plot.position, xend=plot.position, y=0, yend=bs.log2FoldChange), size = 1, color = "#4d4d4d") +
                    geom_point(size = 4, stroke = 1) +
                    scale_shape_manual(values=c("INTRON" = 23, "CDS" = 22, "UTR3" = 25, "UTR5" = 24, "OTHER" = 21)) +
                    scale_fill_manual(values = c("Not" = "#999999", "Up" = "#68b1a5", "Down" = "#874C62")) +
                    scale_color_manual(values = c("Not" = "#4d4d4d", "Up" = "#2b544d", "Down" = "#623747")) +
                    theme_minimal() +
                    guides(color = guide_legend(override.aes = list(size = 4))) +
                    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
                    ylim(-max(abs(df.plot$bs.log2FoldChange)), max(abs(df.plot$bs.log2FoldChange))) +
                    xlim(0,plot.gene.n.tiles) +
                    labs(
                        title = paste0(plot.gene.name, " (", plot.geneID, ")"),
                        x = "Relative position on the gene",
                        y = "Fold-change (log2)",
                        fill = "Regulation",
                        color = "Regulation",
                        shape = "Transcript Region",
                        caption = paste0("Not all binding sites might be shown due to: transcript.regions.outlier.handle=", transcript.regions.outlier.handle)
                    )
                return(p)
            }

            if (transcript.regions.outlier.handle == "second") {
                df.duplicated = mcols(plot.bs) %>% as.data.frame() %>% filter(grepl(",", transcriptRegion)) %>% mutate(duplicated = "yes") %>%
                    mutate(transcriptRegion = sapply(strsplit(transcriptRegion,","), `[`, 2))
                df.single = mcols(plot.bs) %>% as.data.frame() %>% filter(!grepl(",", transcriptRegion)) %>% mutate(duplicated = "no")
                df.plot = rbind.data.frame(df.duplicated, df.single) %>%
                    mutate(sig = (ifelse(bs.padj < alpha & abs(bs.log2FoldChange) > log2(lfc.cutoff), TRUE, FALSE))) %>%
                    mutate(dir = factor(ifelse(bs.log2FoldChange > 0, "Up", "Down"), levels = c("Up", "Down"))) %>%
                    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
                    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
                    mutate(transcriptRegion = stringr::str_trim(transcriptRegion, "left")) %>%
                    mutate(transcriptRegion = stringr::str_trim(transcriptRegion, "right"))

                p = ggplot(df.plot, aes(x = plot.position, y = bs.log2FoldChange, color = sigDir, fill = sigDir, shape = transcriptRegion)) +
                    geom_hline(yintercept = 0, color = "#4d4d4d") +
                    geom_segment(aes(x=plot.position, xend=plot.position, y=0, yend=bs.log2FoldChange), size = 1, color = "#4d4d4d") +
                    geom_point(size = 4, stroke = 1) +
                    scale_shape_manual(values=c("INTRON" = 23, "CDS" = 22, "UTR3" = 25, "UTR5" = 24, "OTHER" = 21)) +
                    scale_fill_manual(values = c("Not" = "#999999", "Up" = "#68b1a5", "Down" = "#874C62")) +
                    scale_color_manual(values = c("Not" = "#4d4d4d", "Up" = "#2b544d", "Down" = "#623747")) +
                    theme_minimal() +
                    guides(color = guide_legend(override.aes = list(size = 4))) +
                    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
                    ylim(-max(abs(df.plot$bs.log2FoldChange)), max(abs(df.plot$bs.log2FoldChange))) +
                    xlim(0,plot.gene.n.tiles) +
                    labs(
                        title = paste0(plot.gene.name, " (", plot.geneID, ")"),
                        x = "Relative position on the gene",
                        y = "Fold-change (log2)",
                        fill = "Regulation",
                        color = "Regulation",
                        shape = "Transcript Region",
                        caption = paste0("Not all binding sites might be shown due to: transcript.regions.outlier.handle=", transcript.regions.outlier.handle)
                    )
                return(p)
            }

            if (transcript.regions.outlier.handle == "both") {
                df.duplicated = mcols(plot.bs) %>% as.data.frame() %>% filter(grepl(",", transcriptRegion)) %>% mutate(duplicated = "yes") %>%
                    separate_longer_delim(transcriptRegion, delim = ",")
                df.single = mcols(plot.bs) %>% as.data.frame() %>% filter(!grepl(",", transcriptRegion)) %>% mutate(duplicated = "no")
                df.plot = rbind.data.frame(df.duplicated, df.single) %>%
                    mutate(sig = (ifelse(bs.padj < alpha & abs(bs.log2FoldChange) > log2(lfc.cutoff), TRUE, FALSE))) %>%
                    mutate(dir = factor(ifelse(bs.log2FoldChange > 0, "Up", "Down"), levels = c("Up", "Down"))) %>%
                    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
                    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
                    mutate(transcriptRegion = stringr::str_trim(transcriptRegion, "left")) %>%
                    mutate(transcriptRegion = stringr::str_trim(transcriptRegion, "right"))

                p = ggplot(df.plot, aes(x = plot.position, y = bs.log2FoldChange, color = sigDir, fill = sigDir, shape = transcriptRegion)) +
                    geom_hline(yintercept = 0, color = "#4d4d4d") +
                    geom_segment(aes(x=plot.position, xend=plot.position, y=0, yend=bs.log2FoldChange), size = 1, color = "#4d4d4d") +
                    geom_point(size = 4, stroke = 1) +
                    ggrepel::geom_label_repel(data = df.plot %>% filter(duplicated == "yes"), aes(label = transcriptRegion)) +
                    scale_shape_manual(values=c("INTRON" = 23, "CDS" = 22, "UTR3" = 25, "UTR5" = 24, "OTHER" = 21)) +
                    scale_fill_manual(values = c("Not" = "#999999", "Up" = "#68b1a5", "Down" = "#874C62")) +
                    scale_color_manual(values = c("Not" = "#4d4d4d", "Up" = "#2b544d", "Down" = "#623747")) +
                    theme_minimal() +
                    guides(color = guide_legend(override.aes = list(size = 4))) +
                    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
                    ylim(-max(abs(df.plot$bs.log2FoldChange)), max(abs(df.plot$bs.log2FoldChange))) +
                    xlim(0,plot.gene.n.tiles) +
                    labs(
                        title = paste0(plot.gene.name, " (", plot.geneID, ")"),
                        x = "Relative position on the gene",
                        y = "Fold-change (log2)",
                        fill = "Regulation",
                        color = "Regulation",
                        shape = "Transcript Region",
                        size = "Duplicated Region",
                        caption = paste0("Some binding sites might not be visible due to overplotting: transcript.regions.outlier.handle=", transcript.regions.outlier.handle)
                    )
                return(p)
            }

        }
    }
}