#' @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)
            mateRanges <- unlist(range(stTmp))
            mateGaps <- gaps(GRanges(ranges = ranges(readInfo), seqnames = readInfo$entityId))
            rmap <- mateRanges[as.character(seqnames(mateGaps))]
            mateGaps <- mateGaps[start(rmap) <= start(mateGaps) & end(rmap) >= end(mateGaps)]
            gy <- readInfo$stack[match(as.character(seqnames(mateGaps)), readInfo$entityId)]
            lineCoords <- data.frame(
                x1 = start(mateGaps), y1 = gy, x2 = end(mateGaps) + 1, y2 = gy,
                col = .dpOrDefault(GdObject, c("col.gap", "col"), .DEFAULT_SHADED_COL),
                lwd = .dpOrDefault(GdObject, c("lwd.gap", "lwd"), 1),
                lty = .dpOrDefault(GdObject, c("lty.gap", "lty"), 1),
                alpha = .alpha(GdObject, "gap"), stringsAsFactors = FALSE
            )
            lineCoords <- lineCoords[!duplicated(lineCoords), ]
        } else {
            mateRanges <- setNames(readInfo, readInfo$entityId)
        }
        if (any(readInfo$status != "unmated") && anyDuplicated(readInfo$groupid) != 0) {
            pairGaps <- gaps(GRanges(ranges = ranges(mateRanges), seqnames = readInfo$groupid[match(names(mateRanges), readInfo$entityId)]))
            rmap <- GdObject@stackRanges[as.character(seqnames(pairGaps))]
            pairGaps <- pairGaps[start(rmap) <= start(pairGaps) & end(rmap) >= end(pairGaps)]
            gy <- readInfo$stack[match(as.character(seqnames(pairGaps)), readInfo$groupid)]
            if (length(pairGaps)) {
                pairsCoords <- data.frame(
                    x1 = start(pairGaps) - 1, y1 = gy, x2 = end(pairGaps) + 1, y2 = gy,
                    col = .dpOrDefault(GdObject, c("col.mates", "col"), .DEFAULT_BRIGHT_SHADED_COL),
                    lwd = .dpOrDefault(GdObject, c("lwd.mates", "lwd"), 1),
                    lty = .dpOrDefault(GdObject, c("lty.mates", "lty"), 1),
                    alpha = .alpha(GdObject, "mates"), stringsAsFactors = FALSE
                )
            } else {
                pairsCoords <- NULL
            }
            lineCoords <- rbind(lineCoords, pairsCoords[!duplicated(pairsCoords), ])
        }
        ## Consider the indels if needed
        ## - deletion as lines
        ## - insertions as vertical bars
        showIndels <- .dpOrDefault(GdObject, "showIndels", FALSE)
        delCoords <- NULL
        insCoords <- NULL
        if (showIndels) {
            cigarTmp <- DataFrame(cigar = readInfo$cigar, start = start(readInfo), entityId = readInfo$entityId, groupId = readInfo$groupid)
            cigarTmp <- cigarTmp[order(cigarTmp$entityId, cigarTmp$start), ]
            cigarTmp <- cigarTmp[!duplicated(cigarTmp$entityId), ]
            delGaps <- unlist(cigarRangesAlongReferenceSpace(cigarTmp$cigar, pos = cigarTmp$start, ops = "D", f = as.factor(cigarTmp$entityId)))
            gy <- readInfo$stack[match(names(delGaps), readInfo$entityId)]
            if (length(delGaps)) {
                delCoords <- data.frame(
                    x1 = start(delGaps) + 1, y1 = gy, x2 = end(delGaps) + 1, y2 = gy,
                    col = .dpOrDefault(GdObject, c("col.deletion", "col"), .DEFAULT_BRIGHT_SHADED_COL),
                    lwd = .dpOrDefault(GdObject, c("lwd.deletion", "lwd"), 1),
                    lty = .dpOrDefault(GdObject, c("lty.deletion", "lty"), 1),
                    alpha = .alpha(GdObject, "deletions"), stringsAsFactors = FALSE
                )
                lineCoords <- rbind(delCoords, lineCoords)
                lineCoords <- lineCoords[!duplicated(lineCoords[, c("x1", "y1", "x2", "y2")]), ]
            }
            insGaps <- unlist(cigarRangesAlongReferenceSpace(cigarTmp$cigar, pos = cigarTmp$start, ops = "I", f = as.factor(cigarTmp$entityId)))
            gy <- readInfo$stack[match(names(insGaps), readInfo$entityId)]
            if (length(insGaps)) {
                ## should both x coordinates be equal to start
                insCoords <- data.frame(
                    x1 = start(insGaps), y1 = gy - sh, x2 = start(insGaps), y2 = gy + sh,
                    col = .dpOrDefault(GdObject, c("col.insertion", "col"), .DEFAULT_BRIGHT_SHADED_COL),
                    lwd = .dpOrDefault(GdObject, c("lwd.insertion", "lwd"), 1),
                    lty = .dpOrDefault(GdObject, c("lty.insertion", "lty"), 1),
                    alpha = .alpha(GdObject, "insertions"), stringsAsFactors = FALSE
                )
            }
        }
        ## The mismatch information on the reads if needed
        mmLetters <- NULL
        if (!is.null(mm)) {
            fcol <- .dpOrDefault(GdObject@referenceSequence, "fontcolor", getBioColor("DNA_BASES_N"))
            ccol <- ifelse(rgb2hsv(col2rgb(fcol))["s", ] < 0.5, "black", "white")
            vpl <- vpLocation()
            lwidth <- max(as.numeric(convertUnit(stringWidth(DNA_ALPHABET), "inches"))) * 0.9337632
            lheight <- max(as.numeric(convertUnit(stringHeight(DNA_ALPHABET), "inches")))
            perLetterW <- vpl$isize["width"] / (maxBase - minBase + 1)
            perLetterH <- vpl$isize["height"] / abs(diff(current.viewport()$yscale))
            res <- .pxResolution(coord = "x")
            mw <- res * .dpOrDefault(GdObject, "min.width", 1)
            mwy <- max(1, mw)
            if (nrow(mm)) {
                x <- c(mm$position + rep(c(0, mwy, mwy, 0), each = nrow(mm)))
                y <- c(rep(mm$stack - sh, 2), rep(mm$stack + sh, 2))
                id <- c(rep(seq(max(id, na.rm = TRUE) + 1, len = nrow(mm)), 4))
                gps <- data.frame(
                    col = rep(if (!(lwidth < perLetterW && lheight < perLetterH)) {
                        "transparent"
                    } else {
                        .dpOrDefault(GdObject, "col.mismatch", .DEFAULT_SHADED_COL)
                    }, nrow(mm)),
                    fill = rep(fcol[as.character(mm$base)]),
                    lwd = rep(.dpOrDefault(GdObject, "lwd.mismatch", 1), nrow(mm)),
                    lty = rep(.dpOrDefault(GdObject, "lty.mismatch", 1), nrow(mm)),
                    alpha = rep(.dpOrDefault(GdObject, "alpha.mismatch", 1), nrow(mm)),
                    stringsAsFactors = FALSE
                )
                ## Finally we draw mm (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")
                if (!.dpOrDefault(GdObject, "noLetters", FALSE) && lwidth < perLetterW && lheight < perLetterH) {
                    mmLetters <- data.frame(x = mm$position + 0.5, y = mm$stack, label = mm$base, col = ccol[mm$base], stringsAsFactors = FALSE)
                }
            }
        }
        if (!is.null(lineCoords)) {
            grid.segments(lineCoords$x1, lineCoords$y1, lineCoords$x2, lineCoords$y2,
                gp = gpar(
                    col = lineCoords$col, alpha = unique(lineCoords$alpha),
                    lwd = lineCoords$lwd, lty = lineCoords$lty
                ),
                default.units = "native"
            )
        }
        if (!is.null(insCoords)) {
            grid.segments(insCoords$x1, insCoords$y1, insCoords$x2, insCoords$y2,
                gp = gpar(
                    col = insCoords$col, alpha = unique(insCoords$alpha),
                    lwd = insCoords$lwd, lty = insCoords$lty
                ),
                default.units = "native"
            )
        }
        if (!is.null(mmLetters)) {
            grid.text(
                x = mmLetters$x, y = mmLetters$y, label = mmLetters$label,
                gp = gpar(col = mmLetters$col), default.units = "native"
            )
        }
        popViewport(2)
    }
    ## Eventually we set up the image map
    ## imageMap(GdObject) <- im
    return(invisible(GdObject))
})

## SetAs ---------------------------------------------------------------------
## Show ----------------------------------------------------------------------

#' @describeIn AlignmentsTrack-class Show method.
#' @export
setMethod(
    "show", signature(object = "AlignmentsTrack"),
    function(object) {
        cat(sprintf(
            paste("AlignmentsTrack track '%s' \n",
                "| genome: %s\n",
                "| active chromosome: %s\n",
                "| containing %i read%s\n",
                sep = ""
            ),
            names(object), genome(object),
            gsub("^chr", "", chromosome(object)),
            length(object),
            ifelse(length(object) == 1, "", "s")
        ))
    }
)

#' @describeIn AlignmentsTrack-class Show method.
#' @export
setMethod("show", signature(object = "ReferenceAlignmentsTrack"), function(object) {
    .referenceTrackInfo(object, "ReferenceAlignmentsTrack")
})