BindingSiteFinder: Binding site defintion based on iCLIP data

Documented in bindingSiteDefinednessPlot duplicatedSitesPlot estimateBsWidthPlot geneOverlapsPlot globalScorePlot makeBsSummaryPlot mergeCrosslinkDiagnosticsPlot mergeSummaryPlot processingStepsFlowChart pureClipGlobalFilterPlot quickFigure rangeCoveragePlot reproducibilityCutoffPlot reproducibilityFilterPlot reproducibilitySamplesPlot reproducibilityScatterPlot supportRatioPlot targetGeneSpectrumPlot transcriptRegionOverlapsPlot transcriptRegionSpectrumPlot

#' Plot the PureCLIP score distribution with global cutoff indicator
#'
#' A diagnostic function that plots the PureCLIP score distribution on a log2
#' scale. The function \code{\link{pureClipGlobalFilter}} is expected to be
#' executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{pureClipGlobalFilter}}
#'
#' @import ggplot2
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # apply 5% filter
#' bds = pureClipGlobalFilter(object = bds, cutoff = 0.05)
#' pureClipGlobalFilterPlot(bds)
#'
#' @export
pureClipGlobalFilterPlot <- function(object) {

    # bind locally used variables
    x <- y <- quant <- name <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))
    if (is.null(object@params$pureClipGlobalFilter)) {
        msg0 = paste0("Global filter was not applied yet. Run BSFind() or pureClipGlobalFilter() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$pureClipGlobalFilter$data) |
        is.null(object@plotData$pureClipGlobalFilter$cutoffPoint)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run pureClipGlobalFilter() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$pureClipGlobalFilter
    optstrNice = paste0("Cutoff=", optstr$cutoff, ", scoreMatchCol=", optstr$match.score)

    # get stored plotting data
    dt = data.frame(s = log2(object@plotData$pureClipGlobalFilter$data))
    cutpoint = log2(object@plotData$pureClipGlobalFilter$cutoffPoint)

    # calc density
    df = density(dt$s)
    df = data.frame(x = df$x, y = df$y)

    # set quantiles to show
    probs = seq(from = 0, to = 1, by = 0.1)
    quants = quantile(dt$s, prob = probs)
    df$quant = factor(findInterval(df$x, quants))

    # set cutoff to show
    cutoffNice = paste0(round(optstr$cutoff, 2) * 100, "%")
    dfCutoff = data.frame(pos = cutpoint, name = paste0(round(optstr$cutoff, 2) * 100, "%"))

    # Setup color pal
    colScaleFun = colorRampPalette(RColorBrewer::brewer.pal(9, "Blues"))

    p = ggplot(df, aes(x = x, y = y)) +
        geom_line() +
        geom_ribbon(aes(ymin = 0, ymax = y, fill = quant)) +
        scale_x_continuous(breaks = quants) +
        scale_fill_manual(values = colScaleFun(length(probs)+1)) +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 6), legend.position = "none") +
        geom_vline(xintercept = cutpoint, linetype = "dashed", size = 1, color = "#b35900") +
        labs(
            title = "pureClipGlobalFilterPlot()",
            subtitle = optstrNice,
            x = "Score [log2]", y = "Density"
        ) +
        geom_label(data = dfCutoff, aes(x = pos, y = Inf, label = name), vjust = 1.2, color = "#b35900")
    return(p)
}


#' Plot the number of overlaps when assigning crosslink sites to genes
#'
#' A diagnostic function that plots the number of crosslink sites with their
#' respective overlapping rate. The function \code{\link{pureClipGeneWiseFilter}}
#' is expected to be executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{pureClipGeneWiseFilter}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # Load GRanges with genes
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' # apply 5% gene-wise filter
#' bds = pureClipGeneWiseFilter(object = bds, anno.genes = gns, cutoff = 0.5,
#'  overlaps = "keepSingle")
#' duplicatedSitesPlot(bds)
#'
#' @export
duplicatedSitesPlot <- function(object) {

    # bind locally used variables
    `#N overlaps` <- Freq <- FreqNice <- nRes <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))
    if (is.null(object@params$geneWiseFilter)) {
        msg0 = paste0("Gene-wise filter was not applied yet. Run BSFind() or pureClipGeneWiseFilter() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$geneWiseFilter$data) ) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run pureClipGeneWiseFilter() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$geneWiseFilter
    optstrNice = paste0("Cutoff=", optstr$cutoff, ", overlaps=", optstr$overlaps)

    df = object@plotData$geneWiseFilter$data
    df = df %>% mutate(FreqNice = format(Freq, big.mark = ",", small.mark = "."))

    # Setup color pal
    cols = c("#576F72", rep("lightgrey", nrow(df)-1))

    # Make Plot
    # --------------------------------------------------------------------------
    # Standard plot
    p = ggplot(df, aes(x = `#N overlaps`, y = Freq, fill = `#N overlaps`, label = FreqNice)) +
        geom_col(color = "black") +
        scale_y_log10() +
        scale_fill_manual(values = cols) +
        theme_bw() +
        theme(legend.position = "none") +
        labs(
            title = "dupliatedSitesPlot()",
            subtitle = optstrNice,
            x = "#N overlaps",
            y = "#N (log10)"
        )

    # Add custom options based on overlappingLociOption
    if (optstr$overlaps == "keepSingle") {
        dfResolved = df %>% mutate(nRes = as.numeric(Freq) / as.numeric(`#N overlaps`))
        p = p +
            geom_col(data = dfResolved, aes(x = `#N overlaps`, y = nRes, fill = `#N overlaps`), fill = "#576F72", color = "black") +
            geom_text(data = df, aes(x = `#N overlaps`, y = Freq, label = FreqNice),
                      color = "#b35900", vjust = 1, hjust = 0.8, fontface = "bold", angle = 45)
    }
    if (optstr$overlaps == "removeAll") {
        dfResolved = df %>% mutate(nRes = ifelse(df$`#N overlaps` == 1, TRUE, FALSE))
        p = p +
            geom_col(data = subset(dfResolved, dfResolved$nRes == FALSE),
                     aes(x = `#N overlaps`, y = Freq, fill = `#N overlaps`),
                     fill = "#ECB390", color = "black") +
            geom_text(data = df, aes(x = `#N overlaps`, y = Freq, label = FreqNice),
                      color = "#b35900", vjust = 1, hjust = 0.5, fontface = "bold", angle = 45)
    }
    if (optstr$overlaps == "keepAll") {
        p = p +
            geom_col(data = df, aes(x = `#N overlaps`, y = Freq, fill = `#N overlaps`),
                     fill = "#576F72", color = "black") +
            geom_text(data = df, aes(x = `#N overlaps`, y = Freq, label = FreqNice),
                      color = "#b35900", vjust = 1, hjust = 0.5, fontface = "bold", angle = 45)
    }

    return(p)
}


#' Plot binding site merging diagnostics
#'
#' A diagnostic function that plots the number of regions to merge over the width
#' of these regions for each merging iteration calculated in
#' \code{\link{makeBindingSites}}. The function \code{\link{makeBindingSites}}
#' is expected to be executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{makeBindingSites}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' mergeCrosslinkDiagnosticsPlot(bds)
#'
#' @export
mergeCrosslinkDiagnosticsPlot <- function(object) {

    # bind locally used variables
    w <- s <- iteration <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$makeBindingSites)) {
        msg0 = paste0("Reproducibility filter was not applied yet. Run BSFind() or makeBindingSites() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$makeBindingSites$mergeCsData) |
        is.null(object@plotData$makeBindingSites$data)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run makeBindingSites() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$makeBindingSites
    optstrNice = paste0("bsSize=", optstr$bsSize, ", minWidth=", optstr$minWidth, ", minCrosslinks=", optstr$minCrosslinks)

    dfPlot = object@plotData$makeBindingSites$mergeCsData
    dfPlot = dfPlot %>% mutate(iteration = as.factor(iteration), w = as.factor(w))

    # Make Plot
    # --------------------------------------------------------------------------
    # Standard plot
    cols = c("#716F81","#B97A95","#F6AE99","#F2E1C1","#BEAEE2","#F7DBF0","#CDF0EA","#F9F9F9",rep("grey",5))

    p = ggplot(dfPlot, aes(x = w, y = s, color = iteration, group = iteration)) +
        geom_point(size = 3) +
        geom_line(size = 1) +
        scale_color_manual(values = cols) +
        theme_bw() +
        labs(
            title = "mergeCrosslinkDiagnosticsPlot()",
            subtitle = optstrNice,
            x = "Region-to-fit width",
            y = "Region-to-fit #N",
            color = "Fitting iteration"
        ) +
        theme(legend.position = "bottom") +
        guides(colour = guide_legend(nrow = 1))

    return(p)
}

#' Plot binding site filter diagnostics
#'
#' A diagnostic function that plots the number of binding sites retained
#' after each filtering step calculated in \code{\link{makeBindingSites}}.
#' The function \code{\link{makeBindingSites}} is expected to be executed prior
#' to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{makeBindingSites}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' makeBsSummaryPlot(bds)
#'
#' @export
makeBsSummaryPlot <- function(object) {

    # bind locally used variables
    Option <- nRanges <- slice <- nNice <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$makeBindingSites)) {
        msg0 = paste0("Reproducibility filter was not applied yet. Run BSFind() or makeBindingSites() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$makeBindingSites$mergeCsData) |
        is.null(object@plotData$makeBindingSites$data)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run makeBindingSites() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$makeBindingSites
    optstrNice = paste0("bsSize=", optstr$bsSize, ", minWidth=", optstr$minWidth, ", minCrosslinks=", optstr$minCrosslinks)

    dfPlot = object@plotData$makeBindingSites$data
    dfPlot = dfPlot %>%
        mutate(Option = factor(Option, levels = c(Option))) %>%
        slice(-1) %>%
        mutate(nNice = paste0(format(nRanges, big.mark = ",", decimal.mark = ".")))

    # Make Plot
    # --------------------------------------------------------------------------
    # Standard plot
    cols = c("#262A56","#8D7B68","#A4907C","#C8B6A6","#F1DEC9")

    p = ggplot(dfPlot, aes(x = Option, y = nRanges, fill = Option)) +
        geom_col() +
        scale_fill_manual(values = cols) +
        theme_bw() +
        theme(legend.position = "none") +
        labs(
            title = "makeBsSummaryPlot()",
            subtitle = optstrNice,
            x = "Options",
            y = "Ranges #N"
        ) +
        geom_text(data = dfPlot, aes(x = Option, y = nRanges, label = nNice),
                  color = "#b35900", vjust = 1, hjust = 0.8, fontface = "bold", angle = 45)

    return(p)
}


#' Plot to that shows the crosslink site distribution per replicate
#'
#' A diagnostic function that plots the number of crosslinks sites over the number
#' of crosslink in these sites and highlights the replicate specific
#' reproducibility cutoff that is derived from that distribution.
#' The function \code{\link{reproducibilityFilter}} is expected to be executed prior
#' to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param plotRange numeric; number of crosslinks per sites that should be shown
#' before summing them up
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{reproducibilityFilter}}
#'
#' @import ggplot2
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = reproducibilityFilter(bds)
#' reproducibilityFilterPlot(bds)
#'
#' @export
reproducibilityFilterPlot <- function(object, plotRange = 20) {

    # bind locally used variables
    value <- n <- condtion <- name <- condition <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$reproducibilityFilter)) {
        msg0 = paste0("Reproducibility filter was not applied yet. Run BSFind() or reproducibilityFilter() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$reproducibilityFilterPlot$data) |
        is.null(object@plotData$reproducibilityFilterPlot$cutoffs) |
        is.null(object@plotData$reproducibilitySamplesPlot$data)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run reproducibilityFilter() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$reproducibilityFilter
    optstrNice = paste0("Cutoff=", optstr$cutoff, ", nReps=", optstr$nReps, ", minCrosslinks=", optstr$minCrosslinks)

    dfPlot = object@plotData$reproducibilityFilterPlot$data
    dfPlotCutoffs = object@plotData$reproducibilityFilterPlot$cutoffs

    # Apply plot range zoom in
    dfPlot$value[dfPlot$value > plotRange] = plotRange

    # Make Plot
    # --------------------------------------------------------------------------
    # Standard plot
    p = ggplot(dfPlot, aes(x = value, y = n)) +
        geom_col(color = "#2d5986", aes(fill = condition)) +
        facet_wrap(~name, scales = "free_x") +
        scale_fill_brewer(palette = "Set2") +
        theme_classic() +
        labs(
            title = "reproducibilityFilterPlot()",
            subtitle = optstrNice,
            x = "Crosslinks per site",
            y = "Number of sites",
            fill = "Condition"
        ) +
        geom_vline(data = dfPlotCutoffs, aes(xintercept = value, group = name), color = "#b35900", size = 1, linetype = "dashed") +
        geom_label(data = dfPlotCutoffs, aes(x = value, y = max(dfPlot$n), label = value), alpha = 0.5)
        # geom_label(data = dfPlotCutoffs, aes(x = value, y = 1, label = value), nudge_y = 2)

    return(p)
}


#' UpSet-plot to that shows how each replicate supports binding sites
#'
#' A diagnostic function that plots the set sizes for each replicate, indicating
#' how many binding site the specific replicate supports given its specific
#' threshold. The function \code{\link{reproducibilityFilter}} is expected to be
#' executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param nIntersections numeric; number of intersection to be shown
#' @param show.title logical; if plot title should be visible
#' @param text.size numeric; fontsize of all numbers on axis
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{reproducibilityFilter}},
#' \code{\link{reproducibilityFilterPlot}}
#'
#' @import ComplexHeatmap
#' @importFrom dplyr desc
#' @importFrom utils head tail
#' @importFrom grid gpar
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = reproducibilityFilter(bds)
#' reproducibilitySamplesPlot(bds)
#'
#' @export
reproducibilitySamplesPlot <- function(object,
                                       nIntersections = 20,
                                       show.title = TRUE,
                                       text.size = NULL) {

    # bind locally used variables
    title <- rowTitle <- Freq <- `.` <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$reproducibilityFilter)) {
        msg0 = paste0("Reproducibility filter was not applied yet. Run BSFind() or reproducibilityFilter() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$reproducibilitySamplesPlot$data)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run reproducibilityFilter() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$reproducibilityFilter
    optstrNice = paste0("Cutoff=", paste(optstr$cutoff, collapse = ","),
                        ", nReps=", paste(optstr$nReps, collapse = ","),
                        ", minCrosslinks=", paste(optstr$minCrosslinks, collapse = ","))

    # Get full combination matrix
    df = object@plotData$reproducibilitySamplesPlot$data
    # colnames(df) = colnames(df) %>% as.data.frame() %>% separate(col = ".", sep = "_", into = c(NA, "Sample")) %>% pull(Sample)
    m = make_comb_mat(df)

    rowTitle = "All intersections"

    # Reduce combination matrix for plotting
    if (length(comb_size(m)) > nIntersections) {
        plottingCutoff = comb_size(m) %>%
            as.data.frame() %>% rename("Freq" = ".") %>% arrange(desc(Freq)) %>%
            head(.,nIntersections) %>% tail(1) %>% pull(Freq)
        m = m[comb_size(m) >= plottingCutoff]
        rowTitle = paste0(nIntersections, " largest intersections")
    }

    title = paste0("reproducibilitySamplesPlot()\n", optstrNice)

    # set fontsizes
    if (!is.null(text.size)) {
        fontsize.top = text.size
        fontsize.right = text.size
        fontsize.row.title = text.size
        fontsize.column.title = text.size
        fontsize.row.names = text.size
        fontsize.column.names = text.size
    } else {
        fontsize.top = 8
        fontsize.right = 8
        fontsize.row.title = 12
        fontsize.column.title = 12
        fontsize.row.names = 10
        fontsize.column.names = 10
    }

    UpSet(m, column_title = ifelse(isTRUE(show.title), title, ""),
          row_title = ifelse(isTRUE(show.title), rowTitle, ""),
          comb_order = order(comb_size(m), decreasing = TRUE),
          top_annotation = upset_top_annotation(m, add_numbers = TRUE, annotation_name_gp = gpar(fontsize = 0),
                                                numbers_gp = gpar(col = "#4C6793", fontsize = fontsize.top, fontface = "italic")),
          right_annotation = upset_right_annotation(m, add_numbers = TRUE, annotation_name_gp = gpar(fontsize = 0),
                                                    gp = gpar(fill = "#0B2447"),
                                                    numbers_gp = gpar(col = "#0B2447", fontsize = fontsize.right, fontface = "italic"),
                                                    numbers_rot = 25),
          comb_col = "#4C6793", bg_col = "white", pt_size = unit(.5, "cm"),
          border = T, lwd = 2, bg_pt_col = "grey",
          row_title_gp = gpar(fontsize = fontsize.row.title),
          column_title_gp = gpar(fontsize = fontsize.column.title),
          row_names_gp = gpar(fontsize = fontsize.row.names),
          column_names_gp = gpar(fontsize = fontsize.column.names),
          set_order = seq_along(set_size(m)))

}



#' UpSet-plot to that shows the gene type overlaps
#'
#' A diagnostic function that plots the gene types of binding sites on overlapping
#' loci genes. The function \code{\link{assignToGenes}} is expected to be
#' executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param show.title logical; if plot title should be visible
#' @param text.size numeric; fontsize of all numbers on axis
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{assignToGenes}} \code{\link{targetGeneSpectrumPlot}}
#'
#' @import ComplexHeatmap
#' @importFrom grid gpar
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' # Load GRanges with genes
#' load(list.files(files, pattern = ".rds$", full.names = TRUE)[1])
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#' geneOverlapsPlot(bds)
#'
#' @export
geneOverlapsPlot <- function(object,
                             text.size = NULL,
                             show.title = TRUE) {

    # bind locally used variables
    title <- rowTitle <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$assignToGenes)) {
        msg0 = paste0("Assignment to anno.genes was not applied yet. Run BSFind() or assignToGenes() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$assignToGenes$dataOverlaps) |
        is.null(object@plotData$assignToGenes$dataSpectrum) ) {
        msg1 = paste0("It seems like you used options `remove` or `keep` in assignToGenes(). The plot you try to make is only present for options `frequency` and `hierarchy`. \n")
        stop(msg1)
    }

    optstr = object@params$assignToGenes
    optstrNice = paste0("Source=", optstr$source, ", Overlaps=", optstr$overlaps,
                        ifelse(!rlang::is_empty(optstr$rule), paste0(", rule=", paste(optstr$rule, collapse = ">")), ""))

    # Get full combination matrix
    df = object@plotData$assignToGenes$dataOverlaps
    m = make_comb_mat(df)

    rowTitle = "All intersections"
    title = paste0("geneOverlapsPlot()\n", optstrNice)


    # set fontsizes
    if (!is.null(text.size)) {
        fontsize.top = text.size
        fontsize.right = text.size
        fontsize.row.title = text.size
        fontsize.column.title = text.size
        fontsize.row.names = text.size
        fontsize.column.names = text.size
    } else {
        fontsize.top = 8
        fontsize.right = 8
        fontsize.row.title = 12
        fontsize.column.title = 12
        fontsize.row.names = 10
        fontsize.column.names = 10
    }

    UpSet(m, column_title = ifelse(isTRUE(show.title), title, ""),
          row_title = ifelse(isTRUE(show.title), rowTitle, ""),
          comb_order = order(comb_size(m), decreasing = TRUE),
          top_annotation = upset_top_annotation(m, add_numbers = TRUE, annotation_name_gp = gpar(fontsize = 0),
                                                numbers_gp = gpar(col = "#4C6793", fontsize = fontsize.top, fontface = "italic")),
          right_annotation = upset_right_annotation(m, add_numbers = TRUE, gp = gpar(fill = "#0B2447"), annotation_name_gp = gpar(fontsize = 0),
                                                    numbers_gp = gpar(col = "#0B2447", fontsize = fontsize.right, fontface = "italic"),
                                                    numbers_rot = 25),
          comb_col = "#4C6793", bg_col = "white", pt_size = unit(.5, "cm"),
          row_title_gp = gpar(fontsize = fontsize.row.title),
          column_title_gp = gpar(fontsize = fontsize.column.title),
          row_names_gp = gpar(fontsize = fontsize.row.names),
          column_names_gp = gpar(fontsize = fontsize.column.names),
          border = T, lwd = 2, bg_pt_col = "grey",
          set_order = seq_along(set_size(m)))
}

#' Bar-chart to show the hosting gene types of binding sites
#'
#' A diagnostic function that plots the gene type of the hosting gene for
#' each binding site. The function \code{\link{assignToGenes}} is expected to be
#' executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param showNGroups numeric; the number of different gene types to show
#' @param text.size numeric; the size of the text elments on the plot
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{assignToGenes}} \code{\link{geneOverlapsPlot}}
#'
#' @import ggplot2
#' @importFrom dplyr lead desc mutate arrange
#' @importFrom utils head
#' @importFrom tibble add_row
#'
#' @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])
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#' targetGeneSpectrumPlot(bds)
#'
#' @export
targetGeneSpectrumPlot <- function(object,
                                   showNGroups = 5,
                                   text.size = 4) {

    # bind locally used variables
    GeneType <- Freq <- FreqNice <- nGenes <- nBs <- geneType <- nice <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$assignToGenes)) {
        msg0 = paste0("Assignment to genes was not applied yet. Run BSFind() or assigneToGenes() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$assignToGenes$dataOverlaps) |
        is.null(object@plotData$assignToGenes$dataSpectrum) ) {
        msg1 = paste0("It seems like you used options `remove` or `keep` in assignToGenes(). The plot you try to make is only present for options `frequency` and `hierarchy`. \n")
        stop(msg1)
    }

    optstr = object@params$assignToGenes
    optstrNice = paste0("Source=", optstr$source, ", Overlaps=", optstr$overlaps)

    df = object@plotData$assignToGenes$dataSpectrum

    # ----
    # Apply topN cutoff
    if (nrow(df) < showNGroups) {
        msg = paste0("Only ", nrow(df), " different gene types, displaying ", nrow(df), " gene types, not ", showNGroups, ". \n")
        showNGroups = nrow(df)
        message(msg)
    }
    if (nrow(df) > showNGroups) {
        otherFreq = df %>% arrange(desc(nGenes), desc(nBs)) %>% lead(showNGroups-1)
        dfOther = data.frame(geneType = "Other",
                             nBs = sum(otherFreq$nBs, na.rm = TRUE),
                             nGenes = sum(otherFreq$nGenes, na.rm = TRUE))
        df = df %>% as.data.frame() %>% arrange(desc(nGenes)) %>% head(showNGroups-1)
        df = rbind.data.frame(df, dfOther)
        df = df %>% mutate(nice = paste0(
            "#G:",
            format(nGenes, big.mark = ",", small.mark = "."),
            " (#BS:",
            format(nBs, big.mark = ",", small.mark = "."),
            ")")) %>%
            mutate(geneType = forcats::fct_reorder(geneType, nGenes))
    }
    if (nrow(df) == showNGroups) {
        df = df %>% arrange(desc(nGenes), desc(nBs)) %>%
            mutate(nice = paste0(
                "#G:",
                format(nGenes, big.mark = ",", small.mark = "."),
                " (#BS:",
                format(nBs, big.mark = ",", small.mark = "."),
                ")")) %>%
            mutate(geneType = forcats::fct_reorder(geneType, nGenes))

    }
    # ----
    # Setup color pal
    cols = RColorBrewer::brewer.pal(n = showNGroups+2, name = "Greys")
    cols = cols[1:showNGroups+1]

    # Make Plot
    # --------------------------------------------------------------------------
    p = ggplot(df, aes(x = geneType, y = nGenes, fill = geneType, label = nice)) +
        geom_col(color = "black") +
        scale_y_log10() +
        coord_flip() +
        scale_fill_manual(values = cols) +
        theme_bw() +
        theme(legend.position = "none") +
        geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
        labs(
            title = "targetGeneSpectrumPlot()",
            subtitle = optstrNice,
            x = "",
            y = "#N (log10)"
        )

    return(p)
}



#' UpSet-plot to that shows the transcript region overlaps
#'
#' A diagnostic function that plots the transcript regions of binding sites
#' on overlapping loci. The function \code{\link{assignToTranscriptRegions}} is expected to be
#' executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param show.title logical; if plot title should be visible
#' @param text.size numeric; fontsize of all numbers on axis
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{assignToTranscriptRegions}} \code{\link{transcriptRegionSpectrumPlot}}
#'
#' @import ComplexHeatmap
#' @importFrom grid gpar
#' @importFrom dplyr desc
#'
#' @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])
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = assignToTranscriptRegions(object = bds, anno.transcriptRegionList = regions)
#' transcriptRegionOverlapsPlot(bds)
#'
#' @export
transcriptRegionOverlapsPlot <- function(object,
                                         text.size = NULL,
                                         show.title = TRUE) {

    # bind locally used variables
    title <- rowTitle <- Freq <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$assignToTranscriptRegions)) {
        msg0 = paste0("Reproducibility filter was not applied yet. Run BSFind() or assignToTranscriptRegions() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$assignToTranscriptRegions$dataOverlaps) |
        is.null(object@plotData$assignToTranscriptRegions$dataSpectrum)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run assignToTranscriptRegions() BSFind() before plotting. \n")
        stop(msg1)
    }
    optstr = object@params$assignToTranscriptRegions
    optstrNice = paste0("Source=", optstr$source, ", overlaps=", optstr$overlaps,
                        ifelse(!rlang::is_empty(optstr$rule), paste0(", rule=", paste(optstr$rule, collapse = ">")), ""))

    # Get full combination matrix
    df = object@plotData$assignToTranscriptRegions$dataOverlaps
    m = make_comb_mat(df, mode = "distinct")

    rowTitle = "All intersections"
    title = paste0("transcriptRegionOverlapsPlot()\n", optstrNice)

    # set fontsizes
    if (!is.null(text.size)) {
        fontsize.top = text.size
        fontsize.right = text.size
        fontsize.row.title = text.size
        fontsize.column.title = text.size
        fontsize.row.names = text.size
        fontsize.column.names = text.size
    } else {
        fontsize.top = 8
        fontsize.right = 8
        fontsize.row.title = 12
        fontsize.column.title = 12
        fontsize.row.names = 10
        fontsize.column.names = 10
    }

    UpSet(m, column_title = ifelse(isTRUE(show.title), title, ""),
          row_title = ifelse(isTRUE(show.title), rowTitle, ""),
          comb_order = order(comb_size(m), decreasing = TRUE),
          top_annotation = upset_top_annotation(m, add_numbers = TRUE, annotation_name_gp = gpar(fontsize = 0),,
                                                numbers_gp = gpar(col = "#4C6793", fontsize = fontsize.top, fontface = "italic")),
          right_annotation = upset_right_annotation(m, add_numbers = TRUE, gp = gpar(fill = "#0B2447"), annotation_name_gp = gpar(fontsize = 0),,
                                                    numbers_gp = gpar(col = "#0B2447", fontsize = fontsize.right, fontface = "italic"),
                                                    numbers_rot = 25),
          comb_col = "#4C6793", bg_col = "white", pt_size = unit(.5, "cm"),
          row_title_gp = gpar(fontsize = fontsize.row.title),
          column_title_gp = gpar(fontsize = fontsize.column.title),
          row_names_gp = gpar(fontsize = fontsize.row.names),
          column_names_gp = gpar(fontsize = fontsize.column.names),
          border = T, lwd = 2, bg_pt_col = "grey",
          set_order = seq_along(set_size(m)))
}


#' Bar-chart to show the hosting transcript regions of binding sites
#'
#' A diagnostic function that plots the transcript regions of the hosting gene
#' for each binding site. The function \code{\link{assignToTranscriptRegions}}
#' is expected to be executed prior to calling this plot function.
#'
#' Count frequencies can be normalized to the length of the hosting region with
#' option \code{normalize}. The specific factor how the hosting region length is
#' used is given by \code{normalize.factor}. In the case of
#' \code{normalize.factor = "sum"} binding site frequencies are divided by the
#' summed length of all regions that host the specific binding site.
#'
#' Further with option \code{values} once can indicate whether raw or normalized
#' frequencies should be shown 'as-is' or normalized to 'percentages'.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param values character; if values should be presented 'as-is', that means
#' for example as frequencies in case \code{normalization = FALSE}, or as
#' percentages
#' @param normalize logical; whether to normalize values
#' @param normalize.factor character; indicate by what factor values should be
#' normalized to region length by
#' @param show.others logical; whether to show 'others' category. Has to be false
#' if \code{normalize = TRUE}
#' @param text.size numeric; font size of the numbers to be displayed on each bar
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{assignToTranscriptRegions}} \code{\link{transcriptRegionOverlapsPlot}}
#'
#' @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])
#' bds = makeBindingSites(object = bds, bsSize = 9)
#' bds = assignToGenes(bds, anno.genes = gns)
#' bds = assignToTranscriptRegions(object = bds, anno.transcriptRegionList = regions)
#' transcriptRegionSpectrumPlot(bds)
#'
#' @export
transcriptRegionSpectrumPlot <- function(object,
                                         values = c("asis", "percentage"),
                                         normalize = FALSE,
                                         normalize.factor = c("sum", "median", "mean"),
                                         show.others = FALSE,
                                         text.size = 4) {

    # bind locally used variables
    TranscriptRegion <- FreqNice <- Freq <-  NULL
    per <- perNice <- value <- name <- norm.value <- norm.per <- freq.per <- . <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$assignToTranscriptRegions)) {
        msg0 = paste0("Transcript region assignment was not computed yet. Run BSFind() or assignToTranscriptRegions() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$assignToTranscriptRegions$dataOverlaps) |
        is.null(object@plotData$assignToTranscriptRegions$dataSpectrum)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run assignToTranscriptRegions() or BSFind() before plotting. \n")
        stop(msg1)
    }

    # handle options
    values = match.arg(values, choices = c("asis", "percentage"))
    normalize.factor = match.arg(normalize.factor, choices = c("sum", "median", "mean"))

    # make option string
    optstr = object@params$assignToTranscriptRegions
    optstrNice = paste0("values=", values, ", normalize=", normalize, ifelse(isTRUE(normalize), paste0(", normalize.factor=", normalize.factor), ""),
                        ", overlaps=", optstr$overlaps,
                        ifelse(! rlang::is_empty(optstr$rule), paste0(", rule=", paste(optstr$rule, collapse = ">")), ""))

    # Get full combination matrix
    df = object@plotData$assignToTranscriptRegions$dataSpectrum

    # Nice formatting of standard names
    df = df %>% mutate(TranscriptRegion = ifelse(toupper(TranscriptRegion) == "UTR3", "3'UTR",
                                                 ifelse(toupper(TranscriptRegion) == "UTR5", "5'UTR",
                                                        ifelse(toupper(TranscriptRegion) == "CDS", "CDS",
                                                               ifelse(toupper(TranscriptRegion) == "INTRON", "Intron", TranscriptRegion)))))

    df = df %>% as.data.frame() %>% arrange(desc(Freq)) %>%
        mutate(FreqNice = format(Freq, big.mark = ",", small.mark = ".")) %>%
        mutate(TranscriptRegion = forcats::fct_reorder(TranscriptRegion, Freq))

    # set color scheme
    cols = c("#A7D2CB", "#F2D388", "#C98474", "#874C62", "#576F72", "#B4CDE6", "#7D6E83", "#D0B8A8")
    cols = cols[1:nrow(df)]

    # Make plot variations
    # --------------------------------------------------------------------------
    if (!isTRUE(normalize)) {
        # don't normalize

        if (!isTRUE(show.others)) {
            # remove other category to match raw plot with normalized plot
            df = df %>% drop_na() %>%
                mutate(TranscriptRegion = forcats::fct_reorder(TranscriptRegion, Freq))
        }

        if (values == "asis") {
            # show raw numbers
            p = ggplot(df, aes(x = TranscriptRegion, y = Freq, fill = TranscriptRegion, label = FreqNice)) +
                geom_col(color = "darkgrey", linewidth = 0.5) +
                theme_bw() +
                theme(legend.position = "none") +
                geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                scale_fill_manual(values = cols) +
                coord_flip(clip = "on", expand = TRUE) +
                labs(title = "transcriptRegionSpectrumPlot()",
                     subtitle = optstrNice,
                     x = "Transcript regions",
                     y = "#N Binding sites")
        }
        if (values == "percentage") {
            # show percentages
            df = df %>%
                mutate(per = round(Freq / sum(Freq) * 100, digits = 1)) %>%
                mutate(perNice = paste0(per, "%"))

            p = ggplot(df, aes(x = TranscriptRegion, y = per, fill = TranscriptRegion, label = perNice)) +
                geom_col(color = "darkgrey", linewidth = 0.5) +
                theme_bw() +
                theme(legend.position = "none") +
                geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                scale_fill_manual(values = cols) +
                coord_flip(clip = "on", expand = TRUE) +
                labs(title = "transcriptRegionSpectrumPlot()",
                     subtitle = optstrNice,
                     x = "Transcript regions",
                     y = "% Binding sites")

        }

    }

    if (isTRUE(normalize)) {
        # normalize counts
        df = df %>%
            drop_na() %>%
            pivot_longer(-c(TranscriptRegion, Freq, FreqNice)) %>%
            mutate(norm.value = Freq / value) %>%
            group_by(name) %>%
            mutate(norm.per = round(norm.value / sum(norm.value) * 100, digits = 1),
                   freq.per = round(Freq / sum(Freq) * 100, digits = 1)) %>%
            mutate(diff.per = round(norm.per - freq.per, digits = 1) ) %>%
            mutate(FreqNice = format(Freq, big.mark = ","))

        if (values == "asis") {
            # show raw numbers
            if (normalize.factor == "mean") {
                # normalize by mean length
                p = df %>% filter(name == "norm.hosting.mean") %>%
                    ggplot(., aes(x = TranscriptRegion, y = norm.value, fill = TranscriptRegion, label = round(norm.value, digits = 2))) +
                    geom_col(color = "darkgrey", linewidth = 0.5) +
                    theme_bw() +
                    theme(legend.position = "none") +
                    geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                    scale_fill_manual(values = cols) +
                    coord_flip(clip = "on", expand = TRUE) +
                    labs(title = "transcriptRegionSpectrumPlot()",
                         subtitle = optstrNice,
                         x = "Transcript regions",
                         y = "% Binding sites")
            }
            if (normalize.factor == "median") {
                # normalize by mean length
                p = df %>% filter(name == "norm.hosting.median") %>%
                    ggplot(., aes(x = TranscriptRegion, y = norm.value, fill = TranscriptRegion, label = round(norm.value, digits = 2))) +
                    geom_col(color = "darkgrey", linewidth = 0.5) +
                    theme_bw() +
                    theme(legend.position = "none") +
                    geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                    scale_fill_manual(values = cols) +
                    coord_flip(clip = "on", expand = TRUE) +
                    labs(title = "transcriptRegionSpectrumPlot()",
                         subtitle = optstrNice,
                         x = "Transcript regions",
                         y = "% Binding sites")
            }
            if (normalize.factor == "sum") {
                # normalize by mean length
                p = df %>% filter(name == "norm.hosting.sum") %>%
                    ggplot(., aes(x = TranscriptRegion, y = norm.value, fill = TranscriptRegion, label = round(norm.value, digits = 2))) +
                    geom_col(color = "darkgrey", linewidth = 0.5) +
                    theme_bw() +
                    theme(legend.position = "none") +
                    geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                    scale_fill_manual(values = cols) +
                    coord_flip(clip = "on", expand = TRUE) +
                    labs(title = "transcriptRegionSpectrumPlot()",
                         subtitle = optstrNice,
                         x = "Transcript regions",
                         y = "% Binding sites")
            }

        }
        if (values == "percentage") {
            # show percentages
            if (normalize.factor == "mean") {
                # normalize by mean length
                p = df %>% filter(name == "norm.hosting.mean") %>%
                    ggplot(., aes(x = TranscriptRegion, y = norm.per, fill = TranscriptRegion, label = paste0(round(norm.per, digits = 1), "%"))) +
                    geom_col(color = "darkgrey", linewidth = 0.5) +
                    theme_bw() +
                    theme(legend.position = "none") +
                    geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                    scale_fill_manual(values = cols) +
                    coord_flip(clip = "on", expand = TRUE) +
                    labs(title = "transcriptRegionSpectrumPlot()",
                         subtitle = optstrNice,
                         x = "Transcript regions",
                         y = "% Binding sites")
            }
            if (normalize.factor == "median") {
                # normalize by mean length
                p = df %>% filter(name == "norm.hosting.median") %>%
                    ggplot(., aes(x = TranscriptRegion, y = norm.per, fill = TranscriptRegion, label = paste0(round(norm.per, digits = 1), "%"))) +
                    geom_col(color = "darkgrey", linewidth = 0.5) +
                    theme_bw() +
                    theme(legend.position = "none") +
                    geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                    scale_fill_manual(values = cols) +
                    coord_flip(clip = "on", expand = TRUE) +
                    labs(title = "transcriptRegionSpectrumPlot()",
                         subtitle = optstrNice,
                         x = "Transcript regions",
                         y = "% Binding sites")
            }
            if (normalize.factor == "sum") {
                # normalize by mean length
                p = df %>% filter(name == "norm.hosting.sum") %>%
                    ggplot(., aes(x = TranscriptRegion, y = norm.per, fill = TranscriptRegion, label = paste0(round(norm.per, digits = 1), "%"))) +
                    geom_col(color = "darkgrey", linewidth = 0.5) +
                    theme_bw() +
                    theme(legend.position = "none") +
                    geom_text(aes(y = Inf), color = "#b35900", hjust = 1, fontface = "bold", size = text.size) +
                    scale_fill_manual(values = cols) +
                    coord_flip(clip = "on", expand = TRUE) +
                    labs(title = "transcriptRegionSpectrumPlot()",
                         subtitle = optstrNice,
                         x = "Transcript regions",
                         y = "% Binding sites")
            }
        }

    }

    return(p)
}


#' Plot the PureCLIP score distribution after re-assignment
#'
#' A diagnostic function that plots the PureCLIP score distribution on a log2
#' scale after the re-assignment on binding site level.
#' The function \code{\link{annotateWithScore}} is expected to be executed
#' prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{annotateWithScore}}
#'
#' @import ggplot2
#' @importFrom stats density
#' @importFrom grDevices colorRampPalette
#'
#' @examples
#' # load clip data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' bds1 = makeBindingSites(object = bds, bsSize = 9)
#' bds1 = annotateWithScore(bds1, match.ranges = getRanges(bds))
#' globalScorePlot(bds1)
#'
#' @export
globalScorePlot <- function(object) {

    # bind locally used variables
    x <- y <- quant <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))
    if (is.null(object@params$annotateWithScore)) {
        msg0 = paste0("Score transfer was not applied yet. Run BSFind() or annotateWithScore() to compute values first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$annotateWithScore$data) ) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run annotateWithScore() or BSFind() before plotting. \n")
        stop(msg1)
    }

    optstr = object@params$annotateWithScore
    optstrNice = paste0("match.score=", optstr$match.score, ", match.option=", optstr$match.option)

    # get stored plotting data
    dt = data.frame(s = log2(object@plotData$annotateWithScore$data))

    # calc density
    df = density(dt$s)
    df = data.frame(x = df$x, y = df$y)

    # set quantiles to show
    probs = seq(from = 0, to = 1, by = 0.1)
    quants = quantile(dt$s, prob = probs)
    df$quant = factor(findInterval(df$x, quants))

    # Setup color pal
    colScaleFun = colorRampPalette(RColorBrewer::brewer.pal(9, "Blues"))

    p = ggplot(df, aes(x = x, y = y)) +
        geom_line() +
        geom_ribbon(aes(ymin = 0, ymax = y, fill = quant)) +
        scale_x_continuous(breaks = quants) +
        scale_fill_manual(values = colScaleFun(length(probs)+1)) +
        theme_bw() +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), legend.position = "none") +
        labs(
            title = "globalScorePlot()",
            subtitle = optstrNice,
            x = "Score [log2]", y = "Density"
        )
    return(p)
}


#' Plot the signal-to-flank score for varying gene-wise filter and binding
#' site width
#'
#' A diagnostic function that plots the the signal-to-flank score as a mean
#' for each binding site width and gene-wise filter as indicated when executing
#' \code{\link{estimateBsWidth}}. Additionally a mean of means visualizes the
#' overall trend and a red line indicates the suggested optimal binding site
#' width. The function \code{\link{estimateBsWidth}} is expected to be
#' executed prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{estimateBsWidth}}
#'
#' @import ggplot2
#' @importFrom stats sd
#'
#' @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])
#' bds = estimateBsWidth(bds, anno.genes = gns, est.maxBsWidth = 19,
#'  geneResolution = "coarse", bsResolution = "coarse", est.subsetChromosome = "chr22")
#' estimateBsWidthPlot(bds)
#' @export
estimateBsWidthPlot <- function(object) {

    # bind locally used variables
    bsSize <- ms <- geneWiseFilter <- signalToFlankRatio <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$estimateBsWidth) | is.null(object@params$bsSize) |
        is.null(object@params$geneFilter)) {
        msg0 = paste0("Binding site width was not estimated. Run estimateBsWidth() or BSFind() first. \n")
        stop(msg0)
    }
    if (is.null(object@plotData$estimateBsWidth$data)) {
        msg1 = paste0("It seems like someting went wrong with your data. Please check the input and make sure to run estimateBsWidth() or BSFind() before plotting. \n")
    }

    optstr = object@params$estimateBsWidth
    optstrNice = paste0("bsResolution=", optstr$bsResolution,
                        ", geneResolution=", optstr$geneResolution,
                        ",\noption=", optstr$option,
                        ", subsetChromosome=", ifelse(!is.null(optstr$est.subsetChromosome),
                                                      ifelse(length(optstr$est.subsetChromosome) == 1,
                                                             optstr$est.subsetChromosome, paste(optstr$est.subsetChromosome, collapse = ",")
                                                             ), "None"
                                                      )
                        )

    est.bsSize = object@params$bsSize
    est.GeneFilter = object@params$geneFilter
    df = object@plotData$estimateBsWidth$data

    dfPlot = df %>%
        mutate(bsSize = as.factor(bsSize)) %>%
        mutate(geneWiseFilter = paste0((geneWiseFilter)*100, "%")) %>%
        mutate(geneWiseFilter = factor(geneWiseFilter, levels = c(paste0(seq(0,95, by = 5), "%"))))

    dfMean = dfPlot %>%
        group_by(bsSize) %>%
        summarise(ms = mean(signalToFlankRatio), sd = sd(signalToFlankRatio)) %>%
        mutate(geneWiseFilter = "mean")


    dfAnno = dfPlot %>%
        mutate(bsSize = as.numeric(as.character(bsSize))) %>%
        group_by(geneWiseFilter) %>%
        filter(bsSize == max(bsSize)) %>%
        mutate(bsSize = as.factor(bsSize))

    colfunc = colorRampPalette(c("#dceaef", "#4C6793"))
    cols = colfunc(n = length(unique(dfPlot$geneWiseFilter)))

    p = ggplot(dfMean, aes(x = bsSize, y = ms, group = geneWiseFilter)) +
        geom_line(data = dfPlot, aes(x = bsSize, y = signalToFlankRatio, group = geneWiseFilter, color = geneWiseFilter), linewidth = 0.7) +
        geom_point(data = dfPlot, aes(x = bsSize, y = signalToFlankRatio, group = geneWiseFilter, color = geneWiseFilter), size = 2) +
        scale_colour_manual(values = cols) +
        geom_errorbar(aes(ymin = ms-sd, ymax = ms+sd), color = "#874C62", width = 0.2, linewidth = 1.1) +
        geom_line(color = "#874C62", linewidth = 1.1) +
        geom_point(size = 4, shape = 21, color = "#874C62", fill = "#874C62", stroke = 1) +
        theme_bw() +
        theme(legend.position = "bottom") +
        guides(color=guide_legend(nrow=2,byrow=TRUE)) +
        labs(
            title = "estimateBsWidthPlot()",
            subtitle = optstrNice,
            y = "Signal-to-flank score",
            x = "Binding site width",
            color = "Gene-wise filter"
        ) +
        geom_vline(xintercept = dfMean$bsSize[dfMean$bsSize == est.bsSize], color = "#b35900", linetype = "dashed", size = 1) +
        geom_label(data = dfMean[dfMean$bsSize == est.bsSize,], aes(x = Inf, y = Inf),
                   label = paste0("Size=", est.bsSize, "/ filter=", est.GeneFilter), vjust = 1, size = 3, hjust = 1, alpha = 0.3)
    return(p)
}


#' Step-wise flowchart plot
#'
#' An overview plot that shows all workflow functions that were executed on the
#' current object, with all input and output binding site numbers and major options
#' that were used. The function can be called at any time in the analysis. Most
#' optimal usage is after a full run of the wrapper function \code{\link{BSFind}}.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param size.all numeric; size of all numbers
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{BSFind}}
#'
#' @import ggplot2
#' @importFrom dplyr n lead slice left_join
#' @importFrom utils tail
#'
#' @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])
#' bds = BSFind(bds, anno.genes = gns, anno.transcriptRegionList = regions,
#'  est.subsetChromosome = "chr22")
#' processingStepsFlowChart(bds)
#'
#' @export
processingStepsFlowChart <- function(object, size.all = 3) {

    # bind locally used variables
    xmin <- ymin <- xmax <- ymax <- Type <- NULL
    y <- x <- Step <- from <- to <- s_e <- Options <- `#N In` <- `#N Out` <- NULL
    y1 <- y2 <- id <- x1 <- x2 <- opt <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    # MAIN COMPUTE
    # --------------------------------------------------------------------------
    # get results
    res = object@results
    res = format(res, big.mark = ",", decimal.mark = ".")
    colnames(res) = c("Step", "Type", "#N In", "#N Out", "%", "Options")

    flowNode = res %>%
        mutate(x = 0, y = n():1) %>%
        mutate(y = y *2) %>%
        mutate(xmin = x - 0.45,
               xmax = x + 0.45,
               ymin = y - 0.65,
               ymax = y + 0.65) %>%
        mutate(Type = as.factor(Type)) %>%
        as.data.frame()

    df = res %>%
        mutate(from = Step, to = lead(Step)) %>%
        slice(-n())

    flowEdges = df %>% dplyr::select(from, to) %>%
        mutate(id = row_number()) %>%
        pivot_longer(cols = c("from", "to"),
                     names_to = "s_e",
                     values_to = "Step") %>%
        left_join(flowNode, by = "Step") %>%
        dplyr::select(-c(Step, Type, y, xmin, xmax)) %>%
        mutate(y = ifelse(s_e == "from", ymin, ymax)) %>%
        dplyr::select(-c(ymin, ymax, Options))

    flowArrowIn = flowNode %>%
        mutate(x2 = xmin -0.01, x1 = xmin -0.25, y2 = ymax -0.45, y1 = ymax +0.4) %>%
        mutate(`#N In` = paste0("#N=", `#N In`),`#N Out` = paste0("#N=", `#N Out`))

    flowArrowOut = flowNode %>%
        mutate(x1 = xmin -0.01, x2 = xmin -0.25, y1 = ymax -0.7, y2 = ymax -1.5) %>%
        mutate(`#N In` = paste0("#N=", `#N In`),`#N Out` = paste0("#N=", `#N Out`))

    flowOutLast = flowArrowIn %>% tail(1) %>%
        mutate(y1 = y1-2, y2 = y2-2)

    suppressWarnings({
        flowOption = flowNode %>%
            separate(Options, sep = ",", into = c("opt")) %>% # throws a warning, that additional columns are dropped
            mutate(x2 = xmax+0.25, x1 = xmax +0.25, y2 = ymax -1, y1 = ymax -0.65)

        flowOptionArrow = flowNode %>%
            separate(Options, sep = ",", into = c("opt")) %>% # throws a warning, that additional columns are dropped
            mutate(x2 = xmax+0.01, x1 = xmax +0.25, y2 = ymax -0.85, y1 = ymax -0.65)
    })

    if (!is.null(object@params$estimateBsWidth) |
        !is.null(object@params$calculateSignalToFlankScore)) {
        cols = c("#F2D7D9", "#D3CEDF", "#ECE5C7", "#9CB4CC")
    } else {
        cols = c("#F2D7D9", "#D3CEDF", "#9CB4CC")
    }

    p = ggplot() +
        geom_rect(data = flowNode, aes(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax, fill = Type), color = "#4d4d4d") +
        geom_text(data = flowNode, aes(x = x, y = y, label = Step), color = "#585c45", size = size.all) +
        geom_path(data = flowEdges,
                  mapping = aes(x = x, y = y, group = id),
                  colour = "#585c45",
                  arrow = arrow(length = unit(0.3, "cm"), type = "closed")) +
        geom_curve(data = flowArrowIn, aes(x = x1, y = y1, xend = x2, yend = y2),
                   arrow = arrow(length = unit(0.2, "cm"), type = "closed"), size = 0.4,
                   color = "gray20", curvature = -0.2) +
        geom_curve(data = flowArrowOut, aes(x = x1, y = y1, xend = x2, yend = y2),
                   arrow = arrow(length = unit(0.2, "cm"), type = "closed"), size = 0.4,
                   color = "gray20", curvature = -0.2) +
        geom_label(data = flowArrowIn, aes(x = x1, y = y1, label = `#N In`), size = size.all, hjust = 1, color = "#4d4d4d", fill = "#EDEDED") +
        geom_label(data = flowOutLast, aes(x = x1, y = y1, label = `#N Out`), size = size.all, hjust = 1, color = "#4d4d4d", fill = "#EDEDED") +
        xlim(-1.3,1.3) +
        scale_fill_manual(values = cols) +
        theme_void() +
        labs(title = "BindingSiteFinder Flowchart",
             subtitle = "Use the function processingStepsTable(), with option=extended to get a table with the full overview of all options that \nwere used. ",
             fill = "") +
        theme(legend.position = "bottom",
              plot.title = element_text(hjust = 0.5, color = "#748DA6", size = 20, face = "bold"),
              plot.subtitle = element_text(color = "#748DA6", size = 8, face = "italic")) +
        geom_label(data = flowOption, aes(x = x1, y = y1, label = opt), size = size.all, hjust = 0, color = "#4d4d4d", fill = "#EDEDED", fontface = "italic") +
        geom_curve(data = flowOptionArrow, aes(x = x1, y = y1, xend = x2, yend = y2),
                   arrow = arrow(length = unit(0.2, "cm"), type = "closed"), size = 0.4,
                   color = "gray20", curvature = -0.2)

    return(p)
}


#' Binding site definedness plot
#'
#' Binding site definedness is given by the percent of crosslinks that fall
#' diretly inside the binding site compare to those around the binding site.
#' This plotting function shows the distribution of those percentage values
#' grouped by what is indicated in the \code{by} argument.
#'
#' If \code{by} = 'all', then all binding site are grouped into one distribution.
#' For options 'transcript_region' and 'gene_type' binding sites are split into
#' groups according to the respective assignment. This requires that the respective
#' assignment function was executed on the dataset prior to calling this plot function.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param by character; the option by which the plot should be grouped by.
#' Options are: "all", "transcript_region", "gene_type"
#' @param showN.genes numeric; if \code{by} is `gene_type`, then this argument
#' set the maximum number of groups to be shown in the plot
#' @param show.others logical; whether to show 'others' category.
#'
#' @return a plot of type \code{\link{ggplot}}
#'
#' @seealso \code{\link{BSFind}}, \code{\link{calculateSignalToFlankScore}}
#'
#' @import ggplot2
#' @importFrom ggdist stat_halfeye
#' @importFrom dplyr tally
#'
#' @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])
#' bds = BSFind(bds, anno.genes = gns, anno.transcriptRegionList = regions,
#'  est.subsetChromosome = "chr22")
#' bds = calculateSignalToFlankScore(bds)
#' bindingSiteDefinednessPlot(bds)
#'
#' @export
bindingSiteDefinednessPlot <- function(
        object,
        by = c("all", "transcript_region", "gene_type"),
        showN.genes = 5,
        show.others = FALSE
) {
    # init local vairables
    transcriptRegion <- signalToFlankRatio <- geneType <- NULL

    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    if (is.null(object@params$calculateSignalToFlankScore)) {
        msg0 = paste0("Signal-to-flank ratio was not calculated yet. Run BSFind() or CalculateSignalToFlankScore() to compute values first. \n")
        stop(msg0)
    }

    # handle by options
    by = match.arg(by, choices = c("all", "transcript_region", "gene_type"))


    optstr = object@params$calculateSignalToFlankScore
    optstrNice = paste0("Flank=", optstr$flank)

    # get main plotting dataframe
    df = as.data.frame(mcols(getRanges(object)), row.names = NULL)

    # filter if needed
    if (!isTRUE(show.others)) {
        # remove other category to match raw plot with normalized plot
        df = subset(df, transcriptRegion != "OTHER")
    }

    # Nice formatting of standard names
    df = df %>% mutate(transcriptRegion = ifelse(toupper(transcriptRegion) == "UTR3", "3'UTR",
                                                 ifelse(toupper(transcriptRegion) == "UTR5", "5'UTR",
                                                        ifelse(toupper(transcriptRegion) == "CDS", "CDS",
                                                               ifelse(toupper(transcriptRegion) == "INTRON", "Intron", transcriptRegion)))))

    # set plotting order based on frequency
    df.order = df %>% group_by(transcriptRegion) %>% tally() %>% arrange(desc(n))
    df$transcriptRegion = factor(df$transcriptRegion, levels = rev(df.order$transcriptRegion), ordered = TRUE)

    # check for what plotting options are available
    if (by == "transcript_region") {
        # check if option is really available
        if (is.null(object@params$assignToTranscriptRegions)) {
            msg0 = paste0("Transcript region assignment was not computed yet. Run BSFind() or assignToTranscriptRegions() to compute values first. \n")
            stop(msg0)
        } else {
            # set color scheme
            cols = c("#A7D2CB", "#F2D388", "#C98474", "#874C62", "#576F72", "#B4CDE6", "#7D6E83", "#D0B8A8")
            cols = cols[1:nrow(df)]
            # make plot
            p = ggplot(df, aes(x = transcriptRegion, y = signalToFlankRatio)) +
                geom_boxplot(aes(color = transcriptRegion), width = .15, outlier.shape = NA) +
                ggdist::stat_halfeye(aes(fill = transcriptRegion), adjust = 1.5, width = .6, .width = 0, justification = -.2, point_colour = NA) +
                geom_boxplot(aes(color = transcriptRegion), width = .15, outlier.shape = NA) +
                scale_fill_manual(values = cols) +
                scale_color_manual(values = cols) +
                theme_bw() +
                theme(legend.position = "none") +
                labs(
                    title = "bindingSiteDefinednessPlot()",
                    subtitle = optstrNice,
                    x = "Transcript region",
                    y = "Percent bound"
                )  +
                ylim(0,1) +
                geom_hline(yintercept = mean(df$signalToFlankRatio), linetype = "dashed", color = "#696969") +
                coord_flip()
        }
    }
    if (by == "gene_type") {
        # check if option is really available
        if (is.null(object@params$assignToGenes)) {
            msg0 = paste0("Assignment to genes was not applied yet. Run BSFind() or assigneToGenes() to compute values first. \n")
            stop(msg0)
        } else {
            selGroup = df %>%
                group_by(geneType) %>%
                summarise(n = n()) %>%
                arrange(desc(n)) %>%
                slice_head(n = showN.genes) %>%
                pull(geneType)

            df = df %>% filter(geneType %in% selGroup) %>%
                mutate(geneType  = factor(geneType, levels = rev(selGroup)))

            #  set color scheme
            colfunc = colorRampPalette(c("#cccccc", "#1a1a1a"))
            cols = colfunc(n = showN.genes)

            # make plot
            p = ggplot(df, aes(x = geneType, y = signalToFlankRatio)) +
                ggdist::stat_halfeye(aes(fill = geneType), adjust = 1.5,
                                     width = .6, .width = 0, justification = -.2,
                                     point_colour = NA, slab_color = "#202020",
                                     slab_linewidth = 0.5) +
                geom_boxplot(aes(fill = geneType), width = .15, outlier.shape = NA, color = "#202020") +
                scale_fill_manual(values = cols) +
                scale_color_manual(values = cols) +
                theme_bw() +
                theme(legend.position = "none") +
                labs(
                    title = "bindingSiteDefinednessPlot()",
                    subtitle = optstrNice,
                    x = "Gene type",
                    y = "Percent bound"
                )  +
                ylim(0,1) +
                geom_hline(yintercept = mean(df$signalToFlankRatio), linetype = "dashed", color = "#696969") +
                coord_flip()

        }
    }
    if (by == "all") {
        # make plot
        p = ggplot(df, aes(x = "Total", y = signalToFlankRatio)) +
            ggdist::stat_halfeye(color = "#696969", fill = "#D3D3D3", adjust = 1.5, width = .6, .width = 0, justification = -.2, point_colour = NA) +
            geom_boxplot(color = "#696969", fill = "#D3D3D3", width = .15, outlier.shape = NA) +
            theme_bw() +
            theme(legend.position = "none") +
            labs(
                title = "bindingSiteDefinednessPlot()",
                subtitle = optstrNice,
                x = "All",
                y = "Percent bound"
            ) +
            ylim(0,1) +
            coord_flip()
    }
    return(p)
}



#' Plot crosslink event coverage over binding site range
#'
#' A diagnostic plot function that allows to check the coverage of crosslink
#' events over different merged regions. The coverage is shown as mean over all
#' replicates and conditions, with a standard deviation corridor.
#'
#' If \code{object} is a single BSFDataObject a single coverage plot will be
#' drawn, whereas if it is a list of BSFDataObjects, then faceting is used to
#' make a plot for each list element.
#'
#' @param object a BSFDataSet, or a list of BSFDataSet
#' @param width numeric; set the plotting range to show (in nt)
#' @param show.samples logical; to show individual samples as lines
#' @param subset.chromosome character; subset by a all ranges on the indicated
#' chromosome. Can also be a vector with multiple chromosomes. If NULL then all
#' ranges are being used.
#' @param quiet logical; whether to print messages
#'
#' @return a plot of type \code{ggplot2} displaying the crosslink coverage over
#' the ranges of the given \code{\link{BSFDataSet}}
#'
#' @seealso \code{\link{BSFDataSet}}, \code{\link{makeBindingSites}}
#'
#' @import ggplot2
#' @importFrom matrixStats colSds
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#'
#' # plotting a single object
#' bds <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 1)
#' rangeCoveragePlot(bds, subset.chromosome = "chr22")
#'
#' # plotting multiple objects
#' bds1 <- makeBindingSites(object = bds, bsSize = 3, minWidth = 2,
#' minCrosslinks = 2, minClSites = 1, sub.chr = "chr22")
#' bds2 <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 1, sub.chr = "chr22")
#' l = list(`1. bsSize = 3` = bds1, `2. bsSize = 9` = bds2)
#' rangeCoveragePlot(l, subset.chromosome = "chr22")
#'
#' @export
rangeCoveragePlot <- function(object,
                              width = 20,
                              show.samples = FALSE,
                              subset.chromosome = "chr1",
                              quiet = TRUE
) {

    # bind locally used variables
    position <- sd <- xmin <- xmax <- ymin <- ymax <- NULL
    Sample <- value <- name <- NULL

    # make plot for individual object
    if (!is(object, "list")) {
        stopifnot(is(object, "BSFDataSet"))
        if (!is.numeric(width)) {
            stop("Plotting width has to be numeric")
        }

        # check if subset is part of the seqnames from ranges
        checkRng = getRanges(object)
        if (!all(subset.chromosome %in% levels(seqnames(checkRng)))){
            msg = paste0("Chromosome to estimate on (", subset.chromosome, "), is not included in the ranges: ", paste(levels(seqnames(checkRng)), collapse = ", "))
            stop(msg)
        }

        # limit estimation to a specific chromosome
        if (!is.null(subset.chromosome)) {
            redObj = .subsetByChr(object, chr = subset.chromosome, quiet = quiet)
        } else {
            redObj = object
        }


        nSamples = nrow(getMeta(redObj))
        nBs = length(getRanges(redObj))
        bsSize = unique(width(getRanges(redObj)))
        rng = getRanges(redObj)
        rngResize = resize(rng, width = 1, fix = "center") + width
        newObj = setRanges(redObj, rngResize)

        optstrNice = paste0("#BS=", format(nBs, big.mark = ","), ", bsSize=", bsSize, ", chr=", subset.chromosome)

        # compute coverage
        cov = clipCoverage(newObj, samples.merge = FALSE, ranges.merge = FALSE,
                           positions.merge = TRUE, out.format = "data.frame", quiet = quiet)

        # main data frame to show
        suppressWarnings({
            df = data.frame(
                mean = colMeans(cov),
                sd = colSds(as.matrix(cov)),
                position = c(rev(-seq_len(width)), 0, seq_len(width))
            )
        })

        # make main plot
        p = ggplot() +
            annotate("rect", xmin = (-floor(bsSize / 2)), xmax = (floor(bsSize / 2)),
                     ymin = -Inf, ymax = Inf, alpha = 0.3, color = "grey") +
            geom_ribbon(data = df, aes(x = position, ymin = mean - sd, ymax = mean + sd),
                        fill = "#8cb3d9") +
            geom_line(data = df, aes(x = position, y = mean), size = 1, color = "#2d5986") +
            theme_classic() +
            labs(
                title = "rangeCoveragePlot()",
                subtitle = optstrNice,
                x = "Position relative to binding site center (nt)",
                y = "Number of crosslink events\n(dark blue = mean of replicates)"
            ) +
            theme(legend.position = "none")

        # add replicates to plot
        if (isTRUE(show.samples)) {
            # add sample
            dfReps = cov %>% tibble::rownames_to_column("Sample") %>%
                pivot_longer(-Sample) %>%
                mutate(position = rep(c(rev(-seq_len(width)), 0, seq_len(width)), nSamples))

            # annotate with sample names
            dfNames = dfReps %>% group_by(Sample) %>%
                filter(position == max(position))
            # set color scheme
            colfunc = colorRampPalette(c("#FFC3A1", "#A75D5D"))
            cols = colfunc(n = nSamples)

            p = p +
                scale_color_manual(values = cols) +
                geom_line(data = dfReps, aes(x = position, y = value, color = Sample)) +
                geom_text(data = dfNames, aes(x = position, y = value, label = Sample, color = Sample), nudge_x = 3) +
                theme(legend.position = "none")
        }

    }
    # make plot for a list of BSFDataSet objects
    if (is(object, "list")) {
        stopifnot(vapply(object, function(x) {
            is(x, "BSFDataSet")},
            FUN.VALUE = logical(1))
        )

        objectNames = names(object)
        df = lapply(seq_along(object), function(x) {

            # check if subset is part of the seqnames from ranges
            checkRng = getRanges(object[[x]])
            if (!all(subset.chromosome %in% levels(seqnames(checkRng)))){
                msg = paste0("Chromosome to estimate on (", subset.chromosome, "), is not included in the ranges: ", paste(levels(seqnames(checkRng)), collapse = ", "))
                stop(msg)
            }

            # limit estimation to a specific chromosome
            if (!is.null(subset.chromosome)) {
                redObj = .subsetByChr(object[[x]], chr = subset.chromosome, quiet = quiet)
            } else {
                redObj = object[[x]]
            }

            # set correct ranges
            rngResize = resize(getRanges(redObj), width = 1, fix = "center") + width
            newObj = setRanges(redObj, rngResize)

            # compute coverage
            cov = clipCoverage(newObj, samples.merge = FALSE, ranges.merge = FALSE, positions.merge = TRUE, out.format = "data.frame")

            # main data frame to show
            suppressWarnings({
                df = data.frame(
                    mean = colMeans(cov),
                    sd = colSds(as.matrix(cov)),
                    position = c(rev(-seq_len(width)), 0, seq_len(width)),
                    name = objectNames[[x]]
                )
            })
            # add sample
            dfReps = cov %>% tibble::rownames_to_column("Sample") %>%
                pivot_longer(-Sample) %>%
                select(-name)

            rownames(df) = NULL
            rownames(dfReps) = NULL

            df = cbind(df, dfReps)

            return(df)
        })
        dfPlot = do.call(rbind, df)
        # annotate with sample names
        dfNames = dfPlot %>% group_by(Sample) %>%
            filter(position == max(position))

        bsSize = vapply(object, function(x) {
            unique(width(getRanges(x)))
        }, FUN.VALUE = integer(1))

        dfFrame = data.frame(
            name = names(bsSize),
            xmin = (-floor(bsSize / 2)),
            xmax = (floor(bsSize / 2)),
            ymin = -Inf,
            ymax = Inf
        )

        p = ggplot() +
            geom_rect(data = dfFrame, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
                      alpha = 0.2, color = "grey") +
            geom_ribbon(data = dfPlot, aes(x = position, ymin = mean - sd, ymax = mean + sd),
                        fill = "#8cb3d9") +
            geom_line(data = dfPlot, aes(x = position, y = mean), size = 1, color = "#2d5986") +
            theme_classic() +
            labs(
                title = "rangeCoveragePlot()",
                x = "Position relative to binding site center (nt)",
                y = "Number of crosslink events\n(dark blue = mean of replicates)"
            ) +
            facet_wrap( ~name) +
            theme(legend.position = "none")


        if(isTRUE(show.samples)) {
            # set color scheme
            colfunc = colorRampPalette(c("#FFC3A1", "#A75D5D"))
            cols = colfunc(n = length(unique(dfPlot$Sample)))

            p = p +
                geom_line(data = dfPlot, aes(x = position, y = value, color = Sample)) +
                geom_text(data = dfNames, aes(x = position, y = value, label = Sample, color = Sample), nudge_x = 3, size = 3) +
                scale_color_manual(values = cols)+
                theme(legend.position = "none")
        }
    }
    return(p)
}

#' Plot summarized results of the different binding site merging and filtering
#' steps
#'
#' Bar charts produced for the different filter steps in the binding site
#' merging routine. Depending on the selected option (\code{select}) all or
#' only a user defined filter can be shown.
#'
#' If \code{object} is a single BSFDataObject a single coverage plot will be
#' drawn, whereas if it is a list of BSFDataObjects, then faceting is used to
#' make a plot for each list element.
#'
#' @param object a BSFDataObject, with the makeBindingSites function already run
#' @param select one of "all", "filter", "inputRanges",
#' "minCLSites", "mergeCrosslinkSites", "minCrosslinks",
#' "centerIsClSite" or "centerIsSummit".
#' Defines which parameter is selected for plotting.
#' @param ... further arguments passed to ggplot
#'
#' @return a plot of type ggplot after the \code{\link{makeBindingSites}}
#' function was run
#'
#' @import ggplot2
#'
#' @seealso \code{\link{makeBindingSites}}
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#'
#' # plotting a single object
#' bds0 <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 1)
#' mergeSummaryPlot(bds0)
#'
#' # plotting mulitple obejcts
#' bds1 <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 1, sub.chr = "chr22")
#' bds2 <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 3, sub.chr = "chr22")
#' l = list(`1. bsSize = 3` = bds1, `2. bsSize = 9` = bds2)
#' mergeSummaryPlot(l, width = 20)
#'
#' @export
mergeSummaryPlot <- function(object,
                             select = c(
                                 "all",
                                 "filter",
                                 "inputRanges",
                                 "minClSites",
                                 "mergeCrosslinkSites",
                                 "minCrosslinks",
                                 "centerIsClSite",
                                 "centerIsSummit"
                             ),
                             ...) {
    # bind locally used variables
    Option <- nRanges <- name <- NULL

    # make plot if only a single summary is given
    if (!is(object, "list")) {
        stopifnot(is(object, "BSFDataSet"))
        smry = getSummary(object)
        # filter not possible if only a single summary is given
        select = "all"
        df = smry
        df$Option = factor(df$Option, levels = c(df$Option))
        p = ggplot(df, aes(
            x = Option,
            y = nRanges,
            fill = Option
        ), ...) +
            geom_col() +
            scale_fill_brewer(palette = "Dark2") +
            theme_classic() +
            theme(axis.text.x = element_text(angle = 90),
                  legend.position = "none") +
            ggtitle("Processing steps")

    }
    # make plot if a list of summaries is give
    if (is(object, "list")) {
        select = match.arg(
            select,
            choices = c(
                "all",
                "filter",
                "inputRanges",
                "minClSites",
                "mergeCrosslinkSites",
                "minCrosslinks",
                "centerIsClSite",
                "centerIsSummit"
            )
        )

        stopifnot(all(vapply(object, function(o) {
            (is(o, "BSFDataSet"))},
            FUN.VALUE = logical(1))))

        smry = lapply(object, getSummary)
        df = do.call("rbind", smry)
        df$name = rep(names(smry), each = 6)
        df$Option = factor(df$Option, levels = c(smry[[1]][[1]]))
        if (select == "all") {
            p = ggplot(df, aes(
                x = Option,
                y = nRanges,
                fill = Option
            ), ...) +
                geom_col() +
                scale_fill_brewer(palette = "Dark2") +
                theme_classic() +
                theme(axis.text.x = element_text(angle = 90),
                      legend.position = "none") +
                facet_wrap( ~ name) +
                ggtitle("Processing steps")
        }
        if (select == "filter") {
            df = df[df$Option != c("inputRanges", "mergeCrosslinkSites"), ]
            p = ggplot(df, aes(
                x = Option,
                y = nRanges,
                fill = Option
            ), ...) +
                geom_col() +
                scale_fill_brewer(palette = "Dark2") +
                theme_classic() +
                theme(axis.text.x = element_text(angle = 90),
                      legend.position = "none") +
                facet_wrap( ~ name) +
                ggtitle("Filtering options")
        }
        if (select == "inputRanges" |
            select == "mergeCrosslinkSites" |
            select == "minCrosslinks" |
            select == "centerIsClSite" |
            select == "centerIsSummit" | select == "minClSites") {
            df = df[df$Option == c(select), ]
            p = ggplot(df, aes(
                x = name,
                y = nRanges,
                fill = name
            ), ...length()) +
                geom_col() +
                scale_fill_brewer(palette = "Dark2") +
                theme_classic() +
                theme(axis.text.x = element_text(angle = 90),
                      legend.position = "none") +
                ggtitle(select)
        }

    }

    return(p)
}

#' Plot to that shows how many replicates support each binding site
#'
#' Plotting function for settings specified in
#' \code{\link{reproducibilityFilter}}.
#'
#' @param object a BSFDataSet object
#' @param cutoff a vector of length = 1, or of length = levels(meta$conditions)
#' with a single number (between 0-1) indicating the quantile cutoff
#' @param min.crosslinks numeric of length = 1, defines the lower boundary for
#' the minimum number of crosslinks a binding site has to be supported by all
#' replicates, regardless of the replicate specific quantile threshold
#' @param max.range maximum number of crosslinks per sites that should be shown
#' @param ... further arguments passed to ggplot
#'
#' @return a plot of type \code{ggplot2} showing the per replicate
#' reproducibility cutoffs based on a given quantile threshold
#'
#' @seealso \code{\link{reproducibilityFilter}}
#'
#' @import ggplot2
#' @importFrom ggforce geom_mark_rect
#'
#' @examples
#'
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#'
#' # merge binding sites
#' bds <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 1)
#'
#' # use the same cutoff for both conditions
#' suppressWarnings(reproducibilityCutoffPlot(bds, max.range = 20, cutoff = c(0.05)))
#'
#' # use different cutoffs for each condition
#' suppressWarnings(reproducibilityCutoffPlot(bds, max.range = 20, cutoff = c(0.1)))
#'
#'
#' @export
reproducibilityCutoffPlot <-
    function(object,
             cutoff = 0.05,
             min.crosslinks = 1,
             max.range = 20,
             ...) {

        # bind locally used variables
        value <- name <- condition <- applyTo <- NULL

        # deprecation notice
        .Deprecated("reproducibilityFilterPlot")

        stopifnot(is(object, "BSFDataSet"))

        cond = getMeta(object)$condition
        df = coverageOverRanges(
            object, returnOptions = "merge_positions_keep_replicates",
            quiet = TRUE)
        df = as.data.frame(mcols(df))

        if (length(cutoff) == 1) {
            if(length(levels(cond)) > 1) {
                stop("Only one cutoff is given for multiple conditions.")
            }
            # calculate sample specific thresholds
            qSel = .selectQuantilesSingleCondtion(
                covDf = df,
                userCond = cond,
                userNreps = 1,
                userCutoff = cutoff
            )
            # apply minimal crosslink threshold
            qSel$value = ifelse(qSel$value < min.crosslinks,
                                qSel$value + min.crosslinks, qSel$value)

            # prepare objects for plotting
            df[df > max.range] = max.range
            df = df %>% pivot_longer(everything())
            df$condition = vapply(strsplit(df$name, "_"), `[`, 2,
                                  FUN.VALUE = character(1))

            qSel$lower = ifelse(qSel$value == min.crosslinks, TRUE, FALSE)

            p = ggplot(df,
                       aes(x = value, group = name)) +
                scale_fill_brewer(palette = "Set1") +
                scale_color_brewer(palette = "Set1") +
                geom_bar(fill = "#2d5986", color = "#2d5986") +
                facet_wrap( ~name, scales = "free_x") +
                geom_vline(data = qSel,
                           aes(xintercept = value, group = name),
                           color = "darkgrey", size = 1.5) +
                theme_classic() +
                geom_label(data = qSel, aes(x = value, y = 1, label = value)) +
                labs(
                    title = paste0(paste(paste0("Cutoff: ", unique(qSel$per), " ",
                                                unique(qSel$sel)), collapse = "; "),
                                   "; min.crosslinks = ", min.crosslinks),
                    x = "Crosslinks per site",
                    y = "Number of sites",
                    fill = "condition"
                )

        }
        if (length(cutoff) > 1) {
            if (length(levels(cond)) == 1) {
                stop("Multiple cutoffs are given but only one condition exists")
            }
            # calculate sample specific thresholds
            qSel = .selectQuantilesMultipleConditions(
                covDf = df,
                userCond = cond,
                userNreps = 1,
                userCutoff = cutoff
            )
            # apply minimal crosslink threshold
            qSel$value = ifelse(qSel$value < min.crosslinks,
                                qSel$value + min.crosslinks, qSel$value)

            # prepare objects for plotting
            df[df > max.range] = max.range
            df = df %>% pivot_longer(everything())
            df$condition = vapply(strsplit(df$name, "_"), `[`, 2,
                                  FUN.VALUE = character(1))
            idx = match(df$condition, unique(cond))
            df$cutoff = cutoff[idx]
            rectDf = data.frame(
                condition = cond,
                name = unique(df$name),
                value = 1
            )
            qSel$lower = ifelse(qSel$value == min.crosslinks, TRUE, FALSE)

            p = ggplot(df,
                       aes(x = value, group = name)) +
                geom_rect(data = rectDf,
                          aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
                              fill = condition), alpha = 0.3) +
                scale_fill_brewer(palette = "Set1") +
                # scale_color_brewer(palette = "Set1") +
                geom_bar(fill = "#2d5986", color = "#2d5986") +
                facet_wrap( ~name, scales = "free_x") +
                geom_vline(data = qSel,
                           aes(xintercept = value, group = name,
                               color = applyTo), size = 1.5) +
                theme_classic() +
                geom_label(data = qSel, aes(x = value, y = 1, label = value)) +
                labs(
                    title = paste0(paste(paste0("Cutoff: ", unique(qSel$per), " ",
                                                unique(qSel$sel)), collapse = "; "),
                                   "; min.crosslinks = ", min.crosslinks),
                    x = "Crosslinks per site",
                    y = "Number of sites",
                    fill = "condition"
                ) +
                scale_color_discrete(guide="none")
        }
        return(p)
    }


#' Plot that shows the binding site support ratio
#'
#' Function that shows a ratio to determine how well a given binding
#' site with is supported by the crosslink coverage of the data.
#' Ratios are computed using the \code{\link{supportRatio}} function.
#'
#' The higher the ratio, the more does the given binding site width captures
#' the enrichment of crosslinks compared the the local surrounding. A ratio
#' equal to 1 would mean no enrichemnt at all.
#'
#' @param object a BSFDataSet object
#' @param bsWidths a numeric vector indicating the different binding site
#' width to compute the ratio for
#' @param bsFlank optional; a numeric vector of the same length as
#' \code{bsWidth} used to specify the width of the flanking regions
#' @param ... further arguments passed to \code{makeBindingSites}
#'
#' @return an object of class \code{ggplot2}
#'
#' @import ggplot2
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#'
#' suppressWarnings(supportRatioPlot(bds, bsWidths = c(3,7),
#' minWidth = 1, minClSites = 1, minCrosslinks = 2))
#'
#' @export
supportRatioPlot <- function(object, bsWidths, bsFlank = NA, ...){
    # deprecation notice
    .Deprecated("estimateBsWidthPlot")

    df = supportRatio(object, bsWidths, bsFlank, ...)

    p = ggplot(df, aes(x = bsWidths, y = supportRatio, group = 1)) +
        geom_point() +
        geom_line() +
        labs(
            title = "Binding site widths signal ratios",
            x = "Binding site widths",
            y = "Signal to flank ratio"
        ) +
        theme_bw()
    p
}


#' Plot that shows binding site reproducibility as scatter
#'
#' Function compute the number of crosslinks per binding site on a log2 scale
#' for each sample. Samples are pairwise correlated as a scatter and pairwise
#' pearson correlation is shown.
#'
#' Unlike most plotting functions, this function is actively calculating the
#' values to plot.
#'
#' @param object a BSFDataSet object
#' @param quiet logical; whether to print messages
#'
#' @return an object of class \code{ggplot2}
#'
#' @seealso \code{\link{BSFind}},
#' \code{\link{reproducibilityFilter}}
#'
#' @import ggplot2
#'
#' @examples
#' # load data
#' files <- system.file("extdata", package="BindingSiteFinder")
#' load(list.files(files, pattern = ".rda$", full.names = TRUE))
#' bds <- makeBindingSites(object = bds, bsSize = 9, minWidth = 2,
#' minCrosslinks = 2, minClSites = 3, sub.chr = "chr22")
#' reproducibilityScatterPlot(bds)
#'
#' @export
reproducibilityScatterPlot <- function(object,
                                       quiet = FALSE

){
    # bind locally used variables
    what <- NULL

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

    if (is.null(object@params$reproducibilityFilter)) {
        what = "BEFORE"
        msg1 = paste0("Reproducibility filter was not applied yet.\n")
        msg2 = paste0("Correlations are shown BEFORE reproducibility filtering.\n")
        msg3 = paste0("Run BSFind() or reproducibilityFilter() to show correlations AFTER reproducibility filtering.\n")
        if (!quiet) message(c(msg1, msg2, msg3))
    } else {
        what = "AFTER"
        msg1 = paste0("Reproducibility filter was applied.\n")
        msg2 = paste0("Correlations are shown AFTER reproducibility filtering.\n")
        if (!quiet) message(c(msg1, msg2))
    }
    cov = coverageOverRanges(object, returnOptions = "merge_positions_keep_replicates")
    df = as.data.frame(mcols(cov))
    df = log2(df+1)
    max.value = max(df)
    p = GGally::ggpairs(df,
                        upper = list(continuous = .pairsCorColor),
                        lower = list(continuous = GGally::wrap(.parisCorrelationPlot, max.value = max.value)),
                        diag = list(continuous = .pairsDensity),
                        progress = quiet) +
        labs(title = "reproducibilityScatterPlot()",
             subtitle = paste0("Correlations are ", what, " reproducibility filter was applied."))

    return(p)
}





#' Quick figures
#'
#' Summarize all results in a set of quick figures. Depending on how the
#' function is called a different set of analytic plots are arranged into
#' either a 'main' or 'supplementary' type multi-panel figure.
#'
#' @param object a \code{\link{BSFDataSet}} object
#' @param what character; the plotting option. One of: 'main', 'supp'
#' @param save.filename File name to create on the disc
#' @param save.width numeric; plot size width
#' @param save.height numeric; plot size height
#' @param save.device charcter; Device to use. One of: 'pdf', 'png', ...
#' @param quiet whether to print messages
#' @param ... further arguments passed to \code{\link[ggplot2]{ggsave}}
#'
#' @return a plot
#'
#' @seealso \code{\link{BSFind}}
#'
#' @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])
#' bds = BSFind(bds, anno.genes = gns, anno.transcriptRegionList = regions,
#'  est.subsetChromosome = "chr22")
#' quickFigure(bds)
#'
#' @export
quickFigure <- function(object,
                        what = c("main", "supp"),
                        save.filename = NULL, # if it is NULL, then plot to normal device, if it is a path, then plot to the path
                        save.width = 10,
                        save.height = 12,
                        save.device = "pdf",
                        quiet = TRUE,
                        ...
){
    # INPUT CHECKS
    # --------------------------------------------------------------------------
    stopifnot(is(object, "BSFDataSet"))

    # handle options
    what = match.arg(what, choices = c("main", "supp"))


    # setting ggplot local theme
    theme_clean <- theme(
        plot.title = element_blank(),
        plot.subtitle = element_blank(),
        axis.text = element_text(size = 6),
        axis.title = element_text(size = 6)
    )


    if (what == "main") {

        # overview
        p1 = processingStepsFlowChart(object, size.all = 2) +
            theme(plot.title = element_blank(),
                  plot.subtitle = element_blank())

        # binding sites
        p2 = estimateBsWidthPlot(object) + theme_clean
        p3 = reproducibilitySamplesPlot(object, text.size = 6, show.title = FALSE)
        p3 = ggplotify::as.grob(p3)
        # targets and assignment
        p4 = targetGeneSpectrumPlot(object, text.size = 2) + theme_clean
        p5 = transcriptRegionSpectrumPlot(object, normalize = TRUE, values = "percentage", normalize.factor = "median", text.size = 2) + theme_clean
        p6 = bindingSiteDefinednessPlot(object, by = "transcript_region") + theme_clean

        # stitch plots into patch

        pl = list(p1, p2, p3, p4, p5, p6)
        patch = patchwork::wrap_plots(pl)

        # set layout
        layout <- "
        AAAABBB#
        AAAABBB#
        AAAABBB#
        AAAACCCC
        AAAACCCC
        AAAACCCC
        AAAACCCC
        AAAACCCC
        #DDDEEFF
        #DDDEEFF
        #DDDEEFF
        #DDDEEFF
        #DDDEEFF
        "

        # annotate patches
        patch = patch + patchwork::plot_layout(guides = "collect", design = layout) &
            theme(legend.position = "bottom",
                  legend.title = element_blank(),
                  legend.text = element_text(size = 6))  &
            theme(legend.key.size = unit(3,"mm")) &
            guides(fill = guide_legend(nrow = 4, byrow = TRUE))
        patch = patch +
            patchwork::plot_annotation(
                title = paste0("Binding spectrum of: ", getName(object)),
                subtitle = "Disclaimer: This is an auto generated figure, handle with care.",
                caption = paste0("Binding spectrum of: ", getName(object), "\n",
                                 "A) processingStepsFlowChart, B) estimateBsWidthPlot,
                                 C) reproducibilitySamplesPlot, D) targetGeneSpectrumPlot,
                                 E) transcriptRegionSpectrumPlot, F) bindingSiteDefinednessPlot"),
                tag_levels = "A"
            )

    }
    if (what == "supp") {
        # binding site definition
        p1 = pureClipGlobalFilterPlot(object) + theme_clean
        p2 = duplicatedSitesPlot(object) + theme_clean
        p3 = mergeCrosslinkDiagnosticsPlot(object) + theme_clean
        p4 = makeBsSummaryPlot(object) + theme_clean
        p5 = globalScorePlot(object) + theme_clean
        # targets and assignment
        p6 = reproducibilityFilterPlot(object) + theme_clean
        p7 = reproducibilityScatterPlot(object, quiet = quiet) + theme_clean

        p8 = geneOverlapsPlot(object, text.size = 6, show.title = FALSE)
        p9 = transcriptRegionOverlapsPlot(object, text.size = 6, show.title = FALSE)

        p = reproducibilitySamplesPlot(object, text.size = 6, show.title = FALSE)

        p7 = patchwork::wrap_elements(GGally::ggmatrix_gtable(p7))
        p8 = ggplotify::as.grob(p8)
        p9 = ggplotify::as.grob(p9)

        # stitch plots into patch
        pl = list(p1, p2, p3, p4, p5, p6, p7, p8, p9)
        patch = patchwork::wrap_plots(pl)

        # patch = p1 + p2 + p3 + p4 + p5 +
        #     p6  + patchwork::wrap_elements(GGally::ggmatrix_gtable(p7)) +
        #     ggplotify::as.grob(p8) + ggplotify::as.grob(p9)

        # set layout
        layout <- "
        AAABBBCCCDDDEEE
        AAABBBCCCDDDEEE
        AAABBBCCCDDDEEE
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        FFFFFFGGGGGGGGG
        HHHHHHH#IIIIIII
        HHHHHHH#IIIIIII
        HHHHHHH#IIIIIII
        HHHHHHH#IIIIIII
        "

        patch = patch + patchwork::plot_layout(guides = "collect", design = layout) &
            theme(legend.position = "right",
                  legend.title = element_blank(),
                  legend.text = element_text(size = 6))  &
            theme(legend.key.size = unit(3,"mm")) &
            guides(fill = guide_legend(ncol = 2, byrow = FALSE))
        patch = patch +
            patchwork::plot_annotation(
                title = paste0("Binding spectrum of: ", getName(object)),
                subtitle = "Disclaimer: This is an auto generated figure, handle with care.",
                caption = paste0("Binding spectrum of: ", getName(object), "\n",
                                 "A) pureClipGlobalFilterPlot, B) duplicatedSitesPlot, C) mergeCrosslinkDiagnosticsPlot,
                                 D) makeBsSummaryPlot, E) globalScorePlot, F) reproducibilityFilterPlot,
                                 G) reproducibilityScatterPlot, H) geneOverlapsPlot, I) transcriptRegionOverlapsPlot"),
                tag_levels = "A"
            )

    }

    # manage output
    if (is.null(save.filename)) {
        return(patch)
    } else {
        ggsave(filename = save.filename, plot = patch,
               width = save.width, height = save.height,
               device = save.device, ...)
    }
}
ZarnackGroup/BindingSiteFinder documentation built on May 2, 2024, 12:38 a.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
ZarnackGroup/BindingSiteFinder
Binding site defintion based on iCLIP data

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

R Package Documentation

Browse R Packages

We want your feedback!

ZarnackGroup/BindingSiteFinder Binding site defintion based on iCLIP data

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

R Package Documentation

Browse R Packages

We want your feedback!

ZarnackGroup/BindingSiteFinder
Binding site defintion based on iCLIP data

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