new file mode 100644

@@ -0,0 +1,1046 @@

+#' @include DataTrack-class.R

+#' @include ReferenceTrack-class.R

+#' @include StackedTrack-class.R

+#' @include SequenceTrack-class.R

+NULL

+

+

+setClassUnion("SequenceTrackOrNULL", c("SequenceTrack", "NULL"))

+

+## AlignmentsTrack Class -----------------------------------------------------

+

+#' AlignmentsTrack class and methods

+#'

+#'

+#' A class to represent short sequences that have been aligned to a reference

+#' genome as they are typically generated in next generation sequencing

+#' experiments.

+#'

+#' @name AlignmentsTrack-class

+#'

+#' @param range An optional meta argument to handle the different input types.

+#' If the \code{range} argument is missing, all the relevant information to

+#' create the object has to be provided as individual function arguments

+#' (see below).

+#'

+#' The different input options for \code{range} are:

+#'

+#' \describe{

+#'

+#' \item{}{A \code{character} string: the path to a \code{BAM} file containing

+#' the read alignments. To be precise, this will result in the instantiation of

+#' a \code{ReferenceAlignmentsTrack} object, but for the user this

+#' implementation detail should be of no concern.}

+#'

+#' \item{}{A \code{GRanges} object: the genomic ranges of the individual reads

+#' as well as the optional additional metadata columns \code{id}, \code{cigar},

+#' \code{mapq}, \code{flag}, \code{isize}, \code{groupid}, \code{status},

+#' \code{md} and \code{seqs} (see description of the individual function

+#' parameters below for details). Calling the constructor on a \code{GRanges}

+#' object without further arguments, e.g. \code{AlignmentsTrack(range=obj)} is

+#' equivalent to calling the coerce method \code{as(obj, "AlignmentsTrack")}.}

+#'

+#' \item{}{An \code{\linkS4class{IRanges}} object: almost identical to the

+#' \code{GRanges} case, except that the chromosome and strand information as

+#' well as all additional metadata has to be provided in the separate

+#' \code{chromosome}, \code{strand}, \code{feature}, \code{group} or \code{id}

+#' arguments, because it can not be directly encoded in an \code{IRanges}

+#' object. Note that none of those inputs are mandatory, and if not provided

+#' explicitely the more or less reasonable default values \code{chromosome=NA}

+#' and \code{strand="*"} are used. }

+#'

+#' \item{}{A \code{data.frame} object: the \code{data.frame} needs to contain

+#' at least the two mandatory columns \code{start} and \code{end} with the

+#' range coordinates. It may also contain a \code{chromosome} and a

+#' \code{strand} column with the chromosome and strand information for each

+#' range. If missing it will be drawn from the separate \code{chromosome} or

+#' \code{strand} arguments. In addition, the \code{id}, \code{cigar},

+#' \code{mapq}, \code{flag}, \code{isize}, \code{groupid}, \code{status},

+#' \code{md} and \code{seqs} data can be provided as additional columns. The

+#' above comments about potential default values also apply here.}

+#'

+#' }

+#'

+#' @template AlignmentsTrack-class_param

+#'

+#' @return

+#'

+#' The return value of the constructor function is a new object of class

+#' \code{AlignmentsTrack} or \code{ReferenceAlignmentsTrack}.

+#' @section Objects from the Class:

+#'

+#' Objects can be created using the constructor function

+#' \code{AlignmentsTrack}.

+#'

+#' @author Florian Hahne

+#'

+#' @inherit GdObject-class seealso

+#'

+#' @examples

+#' ## Creating objects

+#' afrom <- 2960000

+#' ato <- 3160000

+#' alTrack <- AlignmentsTrack(system.file(

+#'     package = "Gviz", "extdata",

+#'     "gapped.bam"

+#' ), isPaired = TRUE)

+#' plotTracks(alTrack, from = afrom, to = ato, chromosome = "chr12")

+#'

+#' ## Omit the coverage or the pile-ups part

+#' plotTracks(alTrack,

+#'     from = afrom, to = ato, chromosome = "chr12",

+#'     type = "coverage"

+#' )

+#' plotTracks(alTrack,

+#'     from = afrom, to = ato, chromosome = "chr12",

+#'     type = "pileup"

+#' )

+#'

+#' ## Including sequence information with the constructor

+#' if (require(BSgenome.Hsapiens.UCSC.hg19)) {

+#'     strack <- SequenceTrack(Hsapiens, chromosome = "chr21")

+#'     afrom <- 44945200

+#'     ato <- 44947200

+#'     alTrack <- AlignmentsTrack(system.file(

+#'         package = "Gviz", "extdata",

+#'         "snps.bam"

+#'     ), isPaired = TRUE, referenceSequence = strack)

+#'     plotTracks(alTrack, chromosome = "chr21", from = afrom, to = ato)

+#'

+#'     ## Including sequence information in the track list

+#'     alTrack <- AlignmentsTrack(system.file(

+#'         package = "Gviz", "extdata",

+#'         "snps.bam"

+#'     ), isPaired = TRUE)

+#'     plotTracks(c(alTrack, strack),

+#'         chromosome = "chr21", from = 44946590,

+#'         to = 44946660

+#'     )

+#' }

+#' @importFrom Rsamtools scanBamFlag scanBamHeader scanBam ScanBamParam scanFaIndex scanFa BamFile scanBamWhat bamWhich

+#'

+#' @exportClass AlignmentsTrack

+setClass("AlignmentsTrack",

+    representation = representation(

+        stackRanges = "GRanges",

+        sequences = "DNAStringSet",

+        referenceSequence = "SequenceTrackOrNULL"

+    ),

+    contains = "StackedTrack",

+    prototype = prototype(

+        stacking = "squish",

+        name = "AlignmentsTrack",

+        coverageOnly = FALSE,

+        stackRanges = GRanges(),

+        sequences = Biostrings::DNAStringSet(),

+        referenceSequence = NULL,

+        dp = DisplayPars(

+            alpha.reads = 0.5,

+            alpha.mismatch = 1,

+            cex = 0.7,

+            cex.mismatch = NULL,

+            col.coverage = NULL,

+            col.gap = .DEFAULT_SHADED_COL,

+            col.mates = .DEFAULT_BRIGHT_SHADED_COL,

+            col.deletion = "#000000",

+            col.insertion = "#984EA3",

+            col.mismatch = .DEFAULT_SHADED_COL,

+            col.reads = NULL,

+            col.sashimi = NULL,

+            col = .DEFAULT_SHADED_COL,

+            collapse = FALSE,

+            coverageHeight = 0.1,

+            fill.coverage = NULL,

+            fill.reads = NULL,

+            fill = "#BABABA",

+            fontface.mismatch = 2,

+            lty.coverage = NULL,

+            lty.gap = NULL,

+            lty.mates = NULL,

+            lty.deletion = NULL,

+            lty.insertion = NULL,

+            lty.mismatch = NULL,

+            lty.reads = NULL,

+            lty = 1,

+            lwd.coverage = NULL,

+            lwd.gap = NULL,

+            lwd.mates = NULL,

+            lwd.deletion = NULL,

+            lwd.insertion = NULL,

+            lwd.mismatch = NULL,

+            lwd.reads = NULL,

+            lwd.sashimiMax = 10,

+            lwd = 1,

+            max.height = 10,

+            min.height = 5,

+            minCoverageHeight = 50,

+            minSashimiHeight = 50,

+            noLetters = FALSE,

+            sashimiFilter = NULL,

+            sashimiFilterTolerance = 0L,

+            sashimiHeight = 0.1,

+            sashimiScore = 1,

+            sashimiStrand = "*",

+            sashimiTransformation = NULL,

+            showIndels = FALSE,

+            showMismatches = TRUE,

+            size = NULL,

+            transformation = NULL,

+            type = c("coverage", "pileup")

+        )

+    )

+)

+

+## Initialize ----------------------------------------------------------------

+

+#' @describeIn AlignmentsTrack-class Initialize.

+#' @export

+setMethod("initialize", "AlignmentsTrack", function(.Object, stackRanges = GRanges(), stacks = numeric(), sequences = DNAStringSet(),

+                                                    referenceSequence = NULL, ...) {

+    ## the diplay parameter defaults

+    .makeParMapping()

+    .Object <- .updatePars(.Object, "AlignmentsTrack")

+    .Object@stackRanges <- stackRanges

+    .Object <- callNextMethod(.Object, ...)

+    .Object@stacks <- stacks

+    .Object@sequences <- sequences

+    .Object@referenceSequence <- referenceSequence

+    return(.Object)

+})

+

+## ReferenceAlignmentsTrack Class --------------------------------------------

+

+#' @describeIn AlignmentsTrack-class The file-based version of the `AlignmentsTrack-class`.

+#' @exportClass ReferenceAlignmentsTrack

+setClass("ReferenceAlignmentsTrack", contains = c("AlignmentsTrack", "ReferenceTrack"))

+

+## Initialize ----------------------------------------------------------------

+

+## This just needs to set the appropriate slots that are being inherited

+## from ReferenceTrack because the multiple inheritance has some strange

+## features with regards to method selection

+

+#' @importClassesFrom Biostrings DNAStringSet

+#' @importFrom Biostrings DNAStringSet

+#' @describeIn AlignmentsTrack-class Initialize.

+#' @export

+setMethod("initialize", "ReferenceAlignmentsTrack", function(.Object, stream, reference, mapping = list(),

+                                                             args = list(), defaults = list(), stacks = numeric(),

+                                                             stackRanges = GRanges(), sequences = Biostrings::DNAStringSet(),

+                                                             referenceSequence = NULL, ...) {

+    .Object <- selectMethod("initialize", "ReferenceTrack")(.Object = .Object, reference = reference, stream = stream,

+        mapping = mapping, args = args, defaults = defaults)

+    .Object <- callNextMethod(.Object, ...)

+    .Object@referenceSequence <- referenceSequence

+    return(.Object)

+})

+

+## Constructor ---------------------------------------------------------------

+

+## Constructor

+#' @describeIn AlignmentsTrack-class Constructor for `AlignmentsTrack-class`.

+#' @export

+AlignmentsTrack <- function(range = NULL, start = NULL, end = NULL, width = NULL, strand, chromosome, genome,

+                            stacking = "squish", id, cigar, mapq, flag = scanBamFlag(isUnmappedQuery = FALSE), isize, groupid, status, md, seqs,

+                            name = "AlignmentsTrack", isPaired = TRUE, importFunction, referenceSequence, ...) {

+    ## Some defaults

+    if (missing(importFunction)) {

+        importFunction <- .import.bam.alignments

+    }

+    covars <- .getCovars(range)

+    isStream <- FALSE

+    if (!is.character(range)) {

+        n <- max(c(length(start), length(end), length(width)), nrow(covars))

+        id <- .covDefault(id, covars[["id"]], paste("read", seq_len(n), sep = "_"))

+        cigar <- .covDefault(cigar, covars[["cigar"]], paste(if (is(range, "GRangesOrIRanges")) width(range) else width, "M", sep = ""))

+        mapq <- .covDefault(mapq, covars[["mapq"]], rep(as.integer(NA), n))

+        flag <- .covDefault(flag, covars[["flag"]], rep(as.integer(NA), n))

+        isize <- .covDefault(isize, covars[["isize"]], rep(as.integer(NA), n))

+        groupid <- .covDefault(groupid, covars[["groupid"]], seq_len(n))

+        md <- .covDefault(md, covars[["md"]], rep(as.character(NA), n))

+        status <- .covDefault(status, covars[["status"]], ifelse(groupid %in% groupid[duplicated(groupid)], "mated", "unmated"))

+    }

+    ## Build a GRanges object from the inputs

+    .missingToNull(c(

+        "strand", "chromosome", "importFunction", "genome", "id", "cigar", "mapq", "flag", "isize", "groupid", "status",

+        "md", "seqs", "referenceSequence"

+    ))

+    args <- list(

+        id = id, cigar = cigar, mapq = mapq, flag = flag, isize = isize, groupid = groupid, status = status, strand = strand, md = md,

+        chromosome = chromosome, genome = genome

+    )

+    defs <- list(

+        strand = "*", chromosome = "chrNA", genome = NA, id = as.character(NA), cigar = as.character(NA), mapq = as.integer(NA),

+        flag = as.integer(NA), isize = as.integer(NA), groupid = as.character(NA), status = as.character(NA), md = as.character(NA)

+    )

+    range <- .buildRange(

+        range = range, start = start, end = end, width = width,

+        args = args, defaults = defs, chromosome = chromosome, trackType = "AlignmentsTrack",

+        importFun = importFunction, stream = TRUE, autodetect = TRUE

+    )

+    ## This is going to be a list if we have to stream data from a file, otherwise we can compute some additional values

+    if (is.list(range)) {

+        isStream <- TRUE

+        slist <- range

+        range <- GRanges()

+        stackRanges <- GRanges()

+        stacks <- NULL

+        seqs <- DNAStringSet()

+    } else {

+        if (is.null(seqs)) {

+            seqs <- DNAStringSet(vapply(width(range), function(x) paste(rep("N", x), collapse = ""), character(1)))

+        }

+        addArgs <- list(...)

+        if ("showIndels" %in% names(addArgs)) {

+            showIndels <- addArgs$showIndels

+        } else {

+            showIndels <- FALSE

+        }

+        tmp <- .computeAlignments(range, drop.D.ranges = showIndels)

+        range <- tmp$range

+        stackRanges <- tmp$stackRange

+        stacks <- tmp$stacks

+    }

+    ## If no chromosome was explicitly asked for we just take the first one in the GRanges object

+    if (missing(chromosome) || is.null(chromosome)) {

+        chromosome <- if (length(range) > 0) .chrName(as.character(seqnames(range)[1])) else "chrNA"

+    }

+    ## And finally the object instantiation

+    genome <- .getGenomeFromGRange(range, ifelse(is.null(genome), character(), genome[1]))

+    if (!isStream) {

+        return(new("AlignmentsTrack",

+            chromosome = chromosome[1], range = range, stacks = stacks,

+            name = name, genome = genome, stacking = stacking, stackRanges = stackRanges, sequences = seqs,

+            referenceSequence = referenceSequence, ...

+        ))

+    } else {

+        ## A bit hackish but for some functions we may want to know which track type we need but at the

+        ## same time we do not want to enforce this as an additional argument

+        e <- new.env()

+        e[["._trackType"]] <- "AlignmentsTrack"

+        e[["._isPaired"]] <- isPaired

+        e[["._flag"]] <- flag

+        environment(slist[["stream"]]) <- e

+        return(new("ReferenceAlignmentsTrack",

+            chromosome = chromosome[1], range = range, stackRanges = stackRanges,

+            name = name, genome = genome, stacking = stacking, stream = slist[["stream"]], reference = slist[["reference"]],

+            mapping = slist[["mapping"]], args = args, defaults = defs, stacks = stacks, referenceSequence = referenceSequence, ...

+        ))

+    }

+}

+

+## General accessors ---------------------------------------------------------

+

+#' @describeIn AlignmentsTrack-class Return all additional annotation

+#' information except for the genomic coordinates for the track items as

+#' a `data.frame`.

+#' @export

+setMethod("values", "AlignmentsTrack", function(x) .dpOrDefault(x, ".__coverage"))

+

+#' @describeIn AlignmentsTrack-class replace the value of the track's chromosome.

+#' This has to be a valid UCSC chromosome identifier or an integer or character

+#' scalar that can be reasonably coerced into one.

+#' @export

+setReplaceMethod("chromosome", "AlignmentsTrack", function(GdObject, value) {

+    GdObject <- callNextMethod()

+    if (!is.null(GdObject@referenceSequence)) {

+        chromosome(GdObject@referenceSequence) <- value[1]

+    }

+    return(GdObject)

+})

+

+

+

+## Annotation Accessors ------------------------------------------------------

+## Stacking ------------------------------------------------------------------

+

+#' @describeIn AlignmentsTrack-class return the stack indices for each track item.

+#' @export

+setMethod(

+    "stacks", "AlignmentsTrack",

+    function(GdObject) if (length(GdObject)) ranges(GdObject)$stack else 0

+)

+

+#' @describeIn AlignmentsTrack-class recompute the stacks based on the available

+#' space and on the object's track items and stacking settings.

+#' @export

+setMethod("setStacks", "AlignmentsTrack", function(GdObject, ...) {

+    if (length(GdObject)) {

+        bins <- if (!.needsStacking(GdObject)) rep(1, length(GdObject)) else disjointBins(GdObject@stackRanges)

+        GdObject@stacks <- bins

+        ranges(GdObject)$stack <- bins[match(ranges(GdObject)$groupid, names(GdObject@stackRanges))]

+    }

+    return(GdObject)

+})

+

+

+## Consolidate ---------------------------------------------------------------

+## Collapse  -----------------------------------------------------------------

+## Subset --------------------------------------------------------------------

+## AlignmentTracks can be subset by using the information in the stackRanges slot, but for the actual reads we need to make sure that

+## we keep all the bits that belong to a given group. We still want to record the requested ranges in the internal '.__plottingRange'

+## display parameter.

+

+#' @describeIn AlignmentsTrack-class Subset a `AlignmentsTrack` by coordinates

+#' and sort if necessary.

+#' @export

+setMethod("subset", signature(x = "AlignmentsTrack"), function(x, from = NULL, to = NULL, stacks = FALSE, use.defaults = TRUE, ...) {

+    ## Subset to a single chromosome first

+    lx <- length(x)

+    csel <- seqnames(x) != chromosome(x)

+    if (any(csel)) {

+        x <- x[!csel]

+    }

+    if (length(x)) {

+        ## Nothing to do if everything is within the range, otherwise we subset the stackRanges and keep all group items

+        ranges <- if (use.defaults) {

+            .defaultRange(x, from = from, to = to)

+        } else {

+            c(

+                from = ifelse(is.null(from), min(start(x@stackRanges)) - 1, from),

+                to = ifelse(is.null(to), max(end(x@stackRanges)) + 1, to)

+            )

+        }

+        displayPars(x) <- list(".__plottingRange" = ranges)

+        sr <- subsetByOverlaps(x@stackRanges, GRanges(seqnames = chromosome(x)[1], ranges = IRanges(start = ranges["from"], end = ranges["to"])))

+        if (length(sr) < length(x@stackRanges)) {

+            x@stackRanges <- sr

+            x@range <- x@range[x@range$groupid %in% names(sr)]

+            if (stacks) {

+                x <- setStacks(x)

+            }

+        }

+    }

+    return(x)

+})

+

+## ReferenceAlignmentsTracks need to stream the data from file and then pass the results on to the next method

+

+#' @describeIn AlignmentsTrack-class Subset a `ReferenceAlignmentsTrack` by coordinates

+#' and sort if necessary.

+#' @export

+setMethod("subset", signature(x = "ReferenceAlignmentsTrack"), function(x, from, to, chromosome, ...) {

+    ## We only need to reach out into the referenced file once if the range is already contained in the object

+    if (missing(from) || is.null(from) || missing(to) || is.null(to)) {

+        stop("Need both start and end location to subset a ReferenceAnnotationTrack")

+    }

+    if (missing(chromosome) || is.null(chromosome)) {

+        chromosome <- Gviz::chromosome(x)

+    } else {

+        chromosome <- .chrName(chromosome)

+    }

+    subRegion <- GRanges(seqnames = chromosome[1], ranges = IRanges(start = from, end = to))

+    oldRange <- .dpOrDefault(x, ".__plottingRange")

+    isIn <- length(ranges(x)) != 0 && !is.null(oldRange) && ranges(subRegion) %within% IRanges(oldRange["from"], oldRange["to"])

+    if (!isIn) {

+        cMap <- .resolveColMapping(x@stream(x@reference, subRegion), x@args, x@mapping)

+        seqs <- cMap$data$seq

+        cMap$data$seq <- NULL

+        range <- .computeAlignments(.buildRange(cMap$data, args = cMap$args, defaults = x@defaults, trackType = "AnnotationTrack"), drop.D.ranges = .dpOrDefault(x, "showIndels", FALSE))

+        ranges(x) <- range$range

+        x@stackRanges <- range$stackRanges

+        x@stacks <- range$stacks

+        x@sequences <- seqs

+    } else {

+        x@sequences <- subseq(x@sequences, start = from - (oldRange["from"] - 1), width = to - from + 1)

+    }

+    chromosome(x) <- chromosome[1]

+    return(callNextMethod(x = x, from = from, to = to, drop = FALSE, ...))

+})

+

+

+## Position ------------------------------------------------------------------

+## DrawGrid ------------------------------------------------------------------

+

+setMethod("drawGrid", signature(GdObject = "AlignmentsTrack"), function(GdObject, from, to) {

+    if (.dpOrDefault(GdObject, "grid", FALSE)) {

+        yvals <- values(GdObject)

+        ylim <- .dpOrDefault(GdObject, "ylim", if (!is.null(yvals) && length(yvals)) {

+            range(yvals, na.rm = TRUE, finite = TRUE)

+        } else {

+            c(-1, 1)

+        })

+        if (diff(ylim) == 0) {

+            ylim <- ylim + c(-1, 1)

+        }

+        yscale <- c(if (is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(ylim), max(ylim) + diff(range(ylim)) * 0.05)

+        covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", 0)

+        pushViewport(viewport(y = 1 - covHeight["npc"], height = covHeight["npc"], just = c(0.5, 0), yscale = yscale, clip = TRUE))

+        panel.grid(

+            v = 0, h = .dpOrDefault(GdObject, "h", -1),

+            col = .dpOrDefault(GdObject, "col.grid", "#e6e6e6"), lty = .dpOrDefault(GdObject, "lty.grid", 1),

+            lwd = .dpOrDefault(GdObject, "lwd.grid", 1)

+        )

+        popViewport(1)

+    }

+})

+

+## DrawAxis ------------------------------------------------------------------

+

+#' @describeIn AlignmentsTrack-class add a y-axis to the title panel of a track.

+#' @export

+setMethod("drawAxis", signature(GdObject = "AlignmentsTrack"), function(GdObject, ...) {

+    type <- match.arg(.dpOrDefault(GdObject, "type", .ALIGNMENT_TYPES), .ALIGNMENT_TYPES, several.ok = TRUE)

+    if ("coverage" %in% type) {

+        yvals <- values(GdObject)

+        ylim <- .dpOrDefault(GdObject, "ylim", if (!is.null(yvals) && length(yvals)) {

+            range(yvals, na.rm = TRUE, finite = TRUE)

+        } else {

+            c(-1, 1)

+        })

+        if (diff(ylim) == 0) {

+            ylim <- ylim + c(-1, 1)

+        }

+        hSpaceAvail <- vpLocation()$isize["width"] / 6

+        yscale <- c(if (is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(ylim), max(ylim) + diff(range(ylim)) * 0.05)

+        col <- .dpOrDefault(GdObject, "col.axis", "white")

+        acex <- .dpOrDefault(GdObject, "cex.axis")

+        acol <- .dpOrDefault(GdObject, "col.axis", "white")

+        at <- pretty(yscale)

+        at <- at[at >= sort(ylim)[1] & at <= sort(ylim)[2]]

+        covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", c(npc = 0, points = 0))

+        covSpace <- .dpOrDefault(GdObject, ".__coverageSpace", 0)

+        pushViewport(viewport(y = 1 - (covHeight["npc"] + covSpace), height = covHeight["npc"], just = c(0.5, 0)))

+        if (is.null(acex)) {

+            vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(at), "inches"))) * length(at) * 1.5

+            hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(at), "inches")))

+            vSpaceAvail <- abs(diff(range(at))) / abs(diff(yscale)) * vpLocation()$isize["height"]

+            acex <- max(0.6, min(vSpaceAvail / vSpaceNeeded, hSpaceAvail / hSpaceNeeded))

+        }

+        vpTitleAxis <- viewport(x = 0.95, width = 0.2, yscale = yscale, just = 0)

+        pushViewport(vpTitleAxis)

+        suppressWarnings(grid.yaxis(gp = gpar(col = acol, cex = acex), at = at))

+        grid.lines(x = c(0, 0), y = ylim, gp = gpar(col = acol), default.units = "native")

+        popViewport(2)

+    } else {

+        covHeight <- c(npc = 0, points = 0)

+        covSpace <- 0

+    }

+    if ("sashimi" %in% type) {

+        sash <- .dpOrDefault(GdObject, ".__sashimi", list(x = numeric(), y = numeric(), id = integer(), score = numeric()))

+        yscale <- if (length(sash$y)) c(-(max(sash$y) + diff(range(sash$y)) * 0.05), 0) else c(-1, 0)

+        ylim <- if (length(sash$y)) c(-max(sash$y), yscale[1] + max(sash$y)) else c(-1, 0)

+        hSpaceAvail <- vpLocation()$isize["width"] / 6

+        col <- .dpOrDefault(GdObject, "col.axis", "white")

+        acex <- .dpOrDefault(GdObject, "cex.axis")

+        acol <- .dpOrDefault(GdObject, "col.axis", "white")

+        labs <- if (length(sash$score)) pretty(c(1, sash$score)) else pretty(c(1, .dpOrDefault(GdObject, ".__sashimiScore", 10)))

+        at <- seq(ylim[1], ylim[2], length.out = length(labs))

+        sashHeight <- .dpOrDefault(GdObject, ".__sashimiHeight", c(npc = 0, points = 0))

+        sashSpace <- .dpOrDefault(GdObject, ".__sashimiSpace", 0)

+        pushViewport(viewport(

+            y = 1 - (sashHeight["npc"] + sashSpace + covHeight["npc"] + covSpace),

+            height = sashHeight["npc"], just = c(0.5, 0)

+        ))

+        if (is.null(acex)) {

+            vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(labs), "inches"))) * length(at) * 1.5

+            hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(labs), "inches")))

+            vSpaceAvail <- abs(diff(range(at))) / abs(diff(yscale)) * vpLocation()$isize["height"]

+            acex <- max(0.6, min(vSpaceAvail / vSpaceNeeded, hSpaceAvail / hSpaceNeeded))

+        }

+        vpTitleAxis <- viewport(x = 0.75, width = 0.2, yscale = yscale, just = 0)

+        pushViewport(vpTitleAxis)

+        suppressWarnings(grid.yaxis(gp = gpar(col = acol, cex = acex), at = at, label = labs))

+        grid.polygon(x = c(0, 0, 1), y = c(ylim[1], ylim[2], ylim[2]), default.units = "native", gp = gpar(col = acol, fill = acol))

+        popViewport(2)

+    }

+    if (.dpOrDefault(GdObject, ".__isCropped", FALSE)) {

+        vspacing <- as.numeric(convertHeight(unit(2, "points"), "npc"))

+        vsize <- as.numeric(convertHeight(unit(4, "points"), "npc"))

+        pushViewport(viewport(height = vsize, y = vspacing, just = c(0.5, 0)))

+        .moreInd(direction = "down", lwd = 2)

+        popViewport(1)

+    }

+})

+

+## DrawGD --------------------------------------------------------------------

+

+#' @describeIn AlignmentsTrack-class plot the object to a graphics device.

+#' The return value of this method is the input object, potentially updated

+#' during the plotting operation. Internally, there are two modes in which the

+#' method can be called. Either in 'prepare' mode, in which case no plotting is

+#' done but the object is preprocessed based on the available space, or in

+#' 'plotting' mode, in which case the actual graphical output is created.

+#' Since subsetting of the object can be potentially costly, this can be

+#' switched off in case subsetting has already been performed before or

+#' is not necessary.

+#'

+#' @importFrom GenomicAlignments cigarRangesAlongReferenceSpace

+#' @export

+setMethod("drawGD", signature("AlignmentsTrack"), function(GdObject, minBase, maxBase, prepare = FALSE, subset = TRUE, ...) {

+    debug <- .dpOrDefault(GdObject, "debug", FALSE)

+    imageMap(GdObject) <- NULL

+    if (!length(GdObject)) {

+        return(invisible(GdObject))

+    }

+    rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)

+    revs <- !.dpOrDefault(GdObject, "reverseStacking", FALSE)

+    vSpacing <- as.numeric(convertHeight(unit(3, "points"), "npc"))

+    type <- match.arg(.dpOrDefault(GdObject, "type", .ALIGNMENT_TYPES), .ALIGNMENT_TYPES, several.ok = TRUE)

+    ylim <- .dpOrDefault(GdObject, "ylim")

+    ## In prepare mode we need to make sure that the stacking information is updated from the optional display parameter (by calling

+    ## the StackedTrack drawGD method), and compute the coverage vector which will be needed for the axis

+    if (prepare) {

+        if ((is.logical(debug) && debug) || "prepare" %in% debug) {

+            browser()

+        }

+        GdObject <- callNextMethod(GdObject, ...)

+        ## The mismatched bases need to be extracted from the read sequences and the reference sequence

+        if (.dpOrDefault(GdObject, "showMismatches", TRUE) && !is.null(GdObject@referenceSequence)) {

+            mm <- .findMismatches(GdObject)

+            if (nrow(mm)) {

+                displayPars(GdObject) <- list(".__mismatches" = mm)

+            }

+        }

+        ## The coverage calculation and the height of the coverage section

+        if ("coverage" %in% type) {

+            covSpace <- as.numeric(convertHeight(unit(5, "points"), "npc"))

+            if ("pileup" %in% type) {

+                covHeight <- .dpOrDefault(GdObject, "coverageHeight", 0.1)

+                if (covHeight > 0 && covHeight < 1) {

+                    covHeight <- as.numeric(convertHeight(unit(covHeight, "npc"), "points"))

+                }

+                covHeight <- max(.dpOrDefault(GdObject, "minCoverageHeight", 50), covHeight)

+                covHeight <- c(points = covHeight, npc = as.numeric(convertHeight(unit(covHeight, "points"), "npc")))

+            } else if ("sashimi" %in% type) {

+                covHeight <- c(npc = 0.5 - (covSpace * 2), points = 1)

+            } else {

+                covHeight <- c(npc = 1 - (covSpace * 2), points = 1)

+            }

+            coverage <- coverage(range(GdObject), width = maxBase)

+            ## apply data transformation if one is set up

+            trans <- .dpOrDefault(GdObject, "transformation")

+            if (is.list(trans)) {

+                trans <- trans[[1]]

+            }

+            if (!is.null(trans)) {

+                if (!is.function(trans) || length(formals(trans)) != 1L) {

+                    stop("Display parameter 'transformation' must be a function with a single argument")

+                }

+                test <- trans(coverage)

+                if (!is(test, "Rle") || length(test) != length(coverage)) {

+                    stop(

+                        "The function in display parameter 'transformation' results in invalid output.\n",

+                        "It has to return a numeric matrix with the same dimensions as the input data."

+                    )

+                }

+                coverage <- test

+            }

+            displayPars(GdObject) <- list(

+                ".__coverage" = coverage,

+                ".__coverageHeight" = covHeight,

+                ".__coverageSpace" = covSpace

+            )

+        } else {

+            covHeight <- c(npc = 0, points = 0)

+            covSpace <- 0

+        }

+        ## The sashimi calculation and the height of the sashimi section

+        if ("sashimi" %in% type) {

+            sashSpace <- as.numeric(convertHeight(unit(5, "points"), "npc"))

+            if ("pileup" %in% type) {

+                sashHeight <- .dpOrDefault(GdObject, "sashimiHeight", 0.1)

+                if (sashHeight > 0 && sashHeight < 1) {

+                    sashHeight <- as.numeric(convertHeight(unit(sashHeight, "npc"), "points"))

+                }

+                sashHeight <- max(.dpOrDefault(GdObject, "minSashimiHeight", 50), sashHeight)

+                sashHeight <- c(points = sashHeight, npc = as.numeric(convertHeight(unit(sashHeight, "points"), "npc")))

+            } else if ("coverage" %in% type) {

+                sashHeight <- c(npc = 0.5 - (sashSpace * 2), points = 1)

+            } else {

+                sashHeight <- c(npc = 1 - (sashSpace * 2), points = 1)

+            }

+            sashScore <- .dpOrDefault(GdObject, "sashimiScore", 1L)

+            sashLwdMax <- .dpOrDefault(GdObject, "lwd.sashimiMax", 10)

+            sashStrand <- .dpOrDefault(GdObject, "sashimiStrand", "*")

+            sashFilter <- .dpOrDefault(GdObject, "sashimiFilter", NULL)

+            sashFilterTolerance <- .dpOrDefault(GdObject, "sashimiFilterTolerance", 0L)

+            sashNumbers <- .dpOrDefault(GdObject, "sashimiNumbers", FALSE)

+            sash <- .dpOrDefault(GdObject, "sashimiJunctions", NULL)

+            if (is.null(sash)) {

+                sash <- .create.summarizedJunctions.for.sashimi.junctions(ranges(GdObject))

+            } else {

+                if (!is(sash, "GRanges")) {

+                    stop("\"sashimiJunctions\" object must be of \"GRanges\" class!")

+                }

+                sashMcolName <- if (sashStrand == "+") "plus_score" else if (sashStrand == "-") "minus_score" else "score"

+                if (sum(colnames(mcols(sash)) == sashMcolName) != 1) {

+                    stop(sprintf("\"mcols\" of \"sashimiJunctions\" object must contain column named \"%s\",\n which matches the specified (%s) \"sashimiStrand\"!", sashMcolName, sashStrand))

+                }

+            }

+            sashTransform <- .dpOrDefault(GdObject, c("sashimiTransformation", "transformation"))

+            sash <- .convert.summarizedJunctions.to.sashimi.junctions(

+                juns = sash,

+                score = sashScore,

+                lwd.max = sashLwdMax,

+                strand = sashStrand,

+                filter = sashFilter,

+                filterTolerance = sashFilterTolerance,

+                trans = sashTransform

+            )

+            displayPars(GdObject) <- list(

+                ".__sashimi" = sash,

+                ".__sashimiHeight" = sashHeight,

+                ".__sashimiSpace" = sashSpace,

+                ".__sashimiNumbers" = sashNumbers

+            )

+        } else {

+            sashHeight <- c(npc = 0, points = 0)

+            sashSpace <- 0

+        }

+        if ("pileup" %in% type) {

+            ## If there are more bins than we can plot we reduce the number until they fit

+            pushViewport(viewport(height = 1 - (covHeight["npc"] + covSpace + sashHeight["npc"] + sashSpace) - vSpacing * 2, y = vSpacing, just = c(0.5, 0)))

+            bins <- ranges(GdObject)$stack

+            curVp <- vpLocation()

+            mih <- min(curVp$size["height"], .dpOrDefault(GdObject, c("min.height", "max.height"), 3))

+            mah <- min(curVp$size["height"], max(mih, .dpOrDefault(GdObject, c("max.height", "min.height"), 8)))

+            bins <- stacks(GdObject)

+            if (curVp$size["height"] / max(bins) < mih) {

+                maxStack <- curVp$size["height"] %/% mih

+                sel <- if (revs) bins <= maxStack else bins > max(bins) - maxStack

+                ranges(GdObject) <- ranges(GdObject)[sel]

+                displayPars(GdObject) <- list(".__isCropped" = TRUE)

+                bins <- stacks(GdObject)

+            }

+            yrange <- range(bins) + c(-0.5, 0.5)

+            add <- max(0, (curVp$size["height"] %/% mah) - max(bins))

+            yrange <- if (revs) c(yrange[1], yrange[2] + add) else c(yrange[1] - add, yrange[2])

+            displayPars(GdObject) <- list(".__yrange" = yrange)

+            popViewport(1)

+        }

+        return(invisible(GdObject))

+    }

+    if ((is.logical(debug) && debug) || "draw" %in% debug) {

+        browser()

+    }

+    mm <- .dpOrDefault(GdObject, ".__mismatches")

+    ## The coverage plot first

+    xscale <- if (!rev) c(minBase, maxBase) else c(maxBase, minBase)

+    readInfo <- ranges(GdObject)

+    covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", c(npc = 0, points = 0))

+    covSpace <- .dpOrDefault(GdObject, ".__coverageSpace", 0)

+    if ("coverage" %in% type) {

+        cov <- .dpOrDefault(GdObject, ".__coverage", Rle(lengths = maxBase, values = as.integer(0)))

+        yscale <- c(if (is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(cov), max(cov) + diff(range(cov)) * 0.05)

+        if (!is.null(ylim)) {

+            yscale <- ylim

+        }

+        vp <- viewport(

+            height = covHeight["npc"], y = 1 - (covHeight["npc"] + covSpace), just = c(0.5, 0), xscale = xscale,

+            yscale = yscale, clip = TRUE

+        )

+        pushViewport(vp)

+        res <- .pxResolution(coord = "x")

+        gp <- gpar(

+            col = .dpOrDefault(GdObject, c("col.coverage", "col"), .DEFAULT_SHADED_COL),

+            fill = .dpOrDefault(GdObject, c("fill.coverage", "fill"), "#BABABA"),

+            lwd = .dpOrDefault(GdObject, c("lwd.coverage", "lwd"), 1),

+            lty = .dpOrDefault(GdObject, c("lty.coverage", "lty"), 1),

+            alpha = .alpha(GdObject)

+        )

+        ## We can compact this if the resolution is not sufficient to speed up the drawing.

+        mminBase <- max(1, minBase)

+        if (res > 2) {

+            brks <- ceiling((maxBase - mminBase) / res)

+            x <- seq(mminBase, maxBase, len = brks)

+            y <- tapply(as.integer(cov[mminBase:maxBase]), cut(mminBase:maxBase, breaks = brks), mean)

+        } else {

+            x <- mminBase:maxBase

+            y <- as.integer(cov[x])

+        }

+        grid.polygon(c(minBase - max(1, res), 0, x, maxBase + max(1, res)), c(0, 0, y, 0), default.units = "native", gp = gp)

+        grid.lines(y = c(0, 0), gp = gpar(col = gp$col, alpha = gp$alpha))

+        if (!is.null(mm)) {

+            fcol <- .dpOrDefault(GdObject@referenceSequence, "fontcolor", getBioColor("DNA_BASES_N"))

+            vpos <- tapply(as.character(mm$base[mm$base != "."]), mm$position[mm$base != "."], table, simplify = FALSE)

+            x <- rep(as.integer(names(vpos)), listLen(vpos))

+            y <- unlist(lapply(vpos, cumsum), use.names = FALSE)

+            col <- fcol[unlist(lapply(vpos, names), use.names = FALSE)]

+            grid.rect(

+                x = x, y = y, height = unlist(vpos, use.names = FALSE), width = 1, default.units = "native", just = c(0, 1),

+                gp = gpar(col = "transparent", fill = col)

+            )

+        }

+        popViewport(1)

+        twoPx <- 2 * as.numeric(convertHeight(unit(1, "points"), "npc"))

+        vp <- viewport(height = twoPx, y = 1 - (covHeight["npc"] + covSpace + twoPx), just = c(0.5, 0))

+        pushViewport(vp)

+        grid.rect(gp = gpar(

+            fill = .dpOrDefault(GdObject, "background.title"), col = "transparent",

+            alpha = .dpOrDefault(GdObject, "alpha.title")

+        ), width = 2)

+        popViewport(1)

+    }

+    ## The sashimi plot as second

+    xscale <- if (!rev) c(minBase, maxBase) else c(maxBase, minBase)

+    sashHeight <- .dpOrDefault(GdObject, ".__sashimiHeight", c(npc = 0, points = 0))

+    sashSpace <- .dpOrDefault(GdObject, ".__sashimiSpace", 0)

+    if ("sashimi" %in% type) {

+        sash <- .dpOrDefault(GdObject, ".__sashimi", list(x = numeric(), y = numeric(), id = integer(), score = numeric(), scaled = numeric()))

+        sashNumbers <- .dpOrDefault(GdObject, ".__sashimiNumbers", FALSE)

+        yscale <- if (length(sash$y)) c(-(max(sash$y) + diff(range(sash$y)) * 0.15), 0) else c(-1, 0) # changed from 0.05 to 0.15 to make sure that numbers fit in the viewport

+        vp <- viewport(

+            height = sashHeight["npc"], y = 1 - (covHeight["npc"] + covSpace + sashHeight["npc"] + sashSpace), just = c(0.5, 0),

+            xscale = xscale, yscale = yscale, clip = TRUE

+        )

+        pushViewport(vp)

+        gp <- gpar(

+            col = .dpOrDefault(GdObject, c("col.sashimi", "col"), .DEFAULT_SHADED_COL),

+            fill = .dpOrDefault(GdObject, c("fill.sashimi", "fill"), "#FFFFFF"),

+            lwd = .dpOrDefault(GdObject, c("lwd.sashimi", "lwd"), 1),

+            lty = .dpOrDefault(GdObject, c("lty.sashimi", "lty"), 1),

+            alpha = .alpha(GdObject)

+        )

+        if (length(sash$x)) {

+            grid.xspline(sash$x, -sash$y,

+                id = sash$id, shape = -1, open = TRUE,

+                default.units = "native", gp = gpar(col = gp$col, lwd = sash$scaled)

+            )

+            ## print the number of reads together with the connecting lines (currently no scaling/resolution)

+            if (sashNumbers) {

+                grid.rect(sash$x[c(FALSE, TRUE, FALSE)], -sash$y[c(FALSE, TRUE, FALSE)],

+                    width = convertUnit(stringWidth(sash$score) * 1.5, "inches"),

+                    height = convertUnit(stringHeight(sash$score) * 1.5, "inches"),

+                    default.units = "native", gp = gpar(col = gp$col, fill = gp$fill)

+                )

+                grid.text(

+                    label = sash$score, sash$x[c(FALSE, TRUE, FALSE)], -sash$y[c(FALSE, TRUE, FALSE)],

+                    default.units = "native", gp = gpar(col = gp$col)

+                )

+            }

+        }

+        popViewport(1)

+    }

+    ## Now the pileups

+    if ("pileup" %in% type) {

+        cex <- max(0.3, .dpOrDefault(GdObject, c("cex.mismatch", "cex"), 0.7))

+        pushViewport(viewport(

+            height = 1 - (covHeight["npc"] + covSpace + sashHeight["npc"] + sashSpace) - vSpacing * 4, y = vSpacing, just = c(0.5, 0),

+            gp = .fontGp(GdObject, "mismatch", cex = cex)

+        ))

+        bins <- stacks(GdObject)

+        stacks <- max(bins)

+        yscale <- .dpOrDefault(GdObject, ".__yrange")

+        if (revs) {

+            yscale <- rev(yscale)

+        }

+        pushViewport(dataViewport(xscale = xscale, extension = 0, yscale = yscale, clip = TRUE))

+        ## Figuring out resolution and the necessary plotting details

+        res <- .pxResolution(coord = "x")

+        ylim <- c(0, 1)

+        h <- diff(ylim)

+        middle <- mean(ylim)

+        sh <- max(0, min(h, .dpOrDefault(GdObject, "stackHeight", 0.75))) / 2

+        boxOnly <- res > 10

+        if (boxOnly) {

+            x <- c(start(readInfo), rep(end(readInfo) + 1, 2), start(readInfo))

+            y <- c(rep(readInfo$stack + sh, 2), rep(readInfo$stack - sh, 2))

+            id <- rep(readInfo$uid, 4)

+        } else {

+            ## We first precompute the coordinates for all the arrow polygons

+            uid2strand <- setNames(as.character(readInfo$readStrand), as.character(readInfo$uid))

+            arrowMap <- unlist(setNames(lapply(split(as.character(readInfo$uid), readInfo$entityId), function(x) {

+                str <- uid2strand[x][1]

+                setNames(str, if (str == "+") tail(x, 1) else head(x, 1))

+            }), NULL))

+            readInfo$arrow <- as.character(NA)

+            readInfo$arrow[match(names(arrowMap), readInfo$uid)] <- arrowMap

+            ## The parts that don't need arrow heads

+            sel <- is.na(readInfo$arrow) | readInfo$arrow == "*"

+            x <- c(start(readInfo)[sel], rep(end(readInfo)[sel] + 1, 2), start(readInfo)[sel])

+            y <- c(rep(readInfo$stack[sel] + sh, 2), rep(readInfo$stack[sel] - sh, 2))

+            id <- rep(readInfo$uid[sel], 4)

+            ## The arrow heads facing right

+            w <- Gviz:::.pxResolution(coord = "x", 5)

+            sel <- readInfo$arrow == "+"

+            ah <- pmax(start(readInfo)[sel], end(readInfo)[sel] + 1 - w)

+            x <- c(x, start(readInfo)[sel], ah, end(readInfo)[sel] + 1, ah, start(readInfo)[sel])

+            y <- c(y, rep(readInfo$stack[sel] + sh, 2), readInfo$stack[sel], rep(readInfo$stack[sel] - sh, 2))

+            id <- c(id, rep(readInfo$uid[sel], 5))

+            ## The arrow heads facing left

+            sel <- readInfo$arrow == "-"

+            ah <- pmin(end(readInfo)[sel], start(readInfo)[sel] + w)

+            x <- c(x, start(readInfo)[sel], ah, rep(end(readInfo)[sel] + 1, 2), ah)

+            y <- c(y, readInfo$stack[sel], rep(readInfo$stack[sel] + sh, 2), rep(readInfo$stack[sel] - sh, 2))

+            id <- c(id, rep(readInfo$uid[sel], 5))

+        }

+        nn <- length(unique(readInfo$uid))

+        gps <- data.frame(

+            col = rep(.dpOrDefault(GdObject, c("col.reads", "col"), .DEFAULT_SHADED_COL), nn),

+            fill = rep(.dpOrDefault(GdObject, c("fill.reads", "fill"), .DEFAULT_BRIGHT_SHADED_COL), nn),

+            lwd = rep(.dpOrDefault(GdObject, c("lwd.reads", "lwd"), 1), nn),

+            lty = rep(.dpOrDefault(GdObject, c("lty.reads", "lty"), 1), nn),

+            alpha = rep(.alpha(GdObject, "reads"), nn), stringsAsFactors = FALSE

+        )

+        ## Finally we draw reads (we need to draw in two steps because of the different alpha levels, reads vs. mismatches)

+        grid.polygon(x = x, y = y, id = id, default.units = "native", gp = gpar(col = gps$col, fill = gps$fill, lwd = gps$lwd, lty = gps$lty, alpha = unique(gps$alpha))) # fix for Sys.setenv(`_R_CHECK_LENGTH_1_CONDITION_`="true"); Sys.setenv(`_R_CHECK_LENGTH_1_LOGIC2_`="true")

+        ## Now the coordinates for the connecting lines

+        lineCoords <- NULL

+        if (anyDuplicated(readInfo$entityId) != 0) {

+            stTmp <- split(readInfo, readInfo$entityId)