Browse code

Restructured the code, split according to the class, first version of roxygen2 documentation. Bumped the version.

Robert Ivánek authored on 21/01/2022 15:52:54
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,1046 @@
1
+#' @include DataTrack-class.R
2
+#' @include ReferenceTrack-class.R
3
+#' @include StackedTrack-class.R
4
+#' @include SequenceTrack-class.R
5
+NULL
6
+
7
+
8
+setClassUnion("SequenceTrackOrNULL", c("SequenceTrack", "NULL"))
9
+
10
+## AlignmentsTrack Class -----------------------------------------------------
11
+
12
+#' AlignmentsTrack class and methods
13
+#'
14
+#'
15
+#' A class to represent short sequences that have been aligned to a reference
16
+#' genome as they are typically generated in next generation sequencing
17
+#' experiments.
18
+#'
19
+#' @name AlignmentsTrack-class
20
+#'
21
+#' @param range An optional meta argument to handle the different input types.
22
+#' If the \code{range} argument is missing, all the relevant information to
23
+#' create the object has to be provided as individual function arguments
24
+#' (see below).
25
+#'
26
+#' The different input options for \code{range} are:
27
+#'
28
+#' \describe{
29
+#'
30
+#' \item{}{A \code{character} string: the path to a \code{BAM} file containing
31
+#' the read alignments. To be precise, this will result in the instantiation of
32
+#' a \code{ReferenceAlignmentsTrack} object, but for the user this
33
+#' implementation detail should be of no concern.}
34
+#'
35
+#' \item{}{A \code{GRanges} object: the genomic ranges of the individual reads
36
+#' as well as the optional additional metadata columns \code{id}, \code{cigar},
37
+#' \code{mapq}, \code{flag}, \code{isize}, \code{groupid}, \code{status},
38
+#' \code{md} and \code{seqs} (see description of the individual function
39
+#' parameters below for details). Calling the constructor on a \code{GRanges}
40
+#' object without further arguments, e.g. \code{AlignmentsTrack(range=obj)} is
41
+#' equivalent to calling the coerce method \code{as(obj, "AlignmentsTrack")}.}
42
+#'
43
+#' \item{}{An \code{\linkS4class{IRanges}} object: almost identical to the
44
+#' \code{GRanges} case, except that the chromosome and strand information as
45
+#' well as all additional metadata has to be provided in the separate
46
+#' \code{chromosome}, \code{strand}, \code{feature}, \code{group} or \code{id}
47
+#' arguments, because it can not be directly encoded in an \code{IRanges}
48
+#' object. Note that none of those inputs are mandatory, and if not provided
49
+#' explicitely the more or less reasonable default values \code{chromosome=NA}
50
+#' and \code{strand="*"} are used. }
51
+#'
52
+#' \item{}{A \code{data.frame} object: the \code{data.frame} needs to contain
53
+#' at least the two mandatory columns \code{start} and \code{end} with the
54
+#' range coordinates. It may also contain a \code{chromosome} and a
55
+#' \code{strand} column with the chromosome and strand information for each
56
+#' range. If missing it will be drawn from the separate \code{chromosome} or
57
+#' \code{strand} arguments. In addition, the \code{id}, \code{cigar},
58
+#' \code{mapq}, \code{flag}, \code{isize}, \code{groupid}, \code{status},
59
+#' \code{md} and \code{seqs} data can be provided as additional columns. The
60
+#' above comments about potential default values also apply here.}
61
+#'
62
+#' }
63
+#'
64
+#' @template AlignmentsTrack-class_param
65
+#'
66
+#' @return
67
+#'
68
+#' The return value of the constructor function is a new object of class
69
+#' \code{AlignmentsTrack} or \code{ReferenceAlignmentsTrack}.
70
+#' @section Objects from the Class:
71
+#'
72
+#' Objects can be created using the constructor function
73
+#' \code{AlignmentsTrack}.
74
+#'
75
+#' @author Florian Hahne
76
+#'
77
+#' @inherit GdObject-class seealso
78
+#'
79
+#' @examples
80
+#' ## Creating objects
81
+#' afrom <- 2960000
82
+#' ato <- 3160000
83
+#' alTrack <- AlignmentsTrack(system.file(
84
+#'     package = "Gviz", "extdata",
85
+#'     "gapped.bam"
86
+#' ), isPaired = TRUE)
87
+#' plotTracks(alTrack, from = afrom, to = ato, chromosome = "chr12")
88
+#'
89
+#' ## Omit the coverage or the pile-ups part
90
+#' plotTracks(alTrack,
91
+#'     from = afrom, to = ato, chromosome = "chr12",
92
+#'     type = "coverage"
93
+#' )
94
+#' plotTracks(alTrack,
95
+#'     from = afrom, to = ato, chromosome = "chr12",
96
+#'     type = "pileup"
97
+#' )
98
+#'
99
+#' ## Including sequence information with the constructor
100
+#' if (require(BSgenome.Hsapiens.UCSC.hg19)) {
101
+#'     strack <- SequenceTrack(Hsapiens, chromosome = "chr21")
102
+#'     afrom <- 44945200
103
+#'     ato <- 44947200
104
+#'     alTrack <- AlignmentsTrack(system.file(
105
+#'         package = "Gviz", "extdata",
106
+#'         "snps.bam"
107
+#'     ), isPaired = TRUE, referenceSequence = strack)
108
+#'     plotTracks(alTrack, chromosome = "chr21", from = afrom, to = ato)
109
+#'
110
+#'     ## Including sequence information in the track list
111
+#'     alTrack <- AlignmentsTrack(system.file(
112
+#'         package = "Gviz", "extdata",
113
+#'         "snps.bam"
114
+#'     ), isPaired = TRUE)
115
+#'     plotTracks(c(alTrack, strack),
116
+#'         chromosome = "chr21", from = 44946590,
117
+#'         to = 44946660
118
+#'     )
119
+#' }
120
+#' @importFrom Rsamtools scanBamFlag scanBamHeader scanBam ScanBamParam scanFaIndex scanFa BamFile scanBamWhat bamWhich
121
+#'
122
+#' @exportClass AlignmentsTrack
123
+setClass("AlignmentsTrack",
124
+    representation = representation(
125
+        stackRanges = "GRanges",
126
+        sequences = "DNAStringSet",
127
+        referenceSequence = "SequenceTrackOrNULL"
128
+    ),
129
+    contains = "StackedTrack",
130
+    prototype = prototype(
131
+        stacking = "squish",
132
+        name = "AlignmentsTrack",
133
+        coverageOnly = FALSE,
134
+        stackRanges = GRanges(),
135
+        sequences = Biostrings::DNAStringSet(),
136
+        referenceSequence = NULL,
137
+        dp = DisplayPars(
138
+            alpha.reads = 0.5,
139
+            alpha.mismatch = 1,
140
+            cex = 0.7,
141
+            cex.mismatch = NULL,
142
+            col.coverage = NULL,
143
+            col.gap = .DEFAULT_SHADED_COL,
144
+            col.mates = .DEFAULT_BRIGHT_SHADED_COL,
145
+            col.deletion = "#000000",
146
+            col.insertion = "#984EA3",
147
+            col.mismatch = .DEFAULT_SHADED_COL,
148
+            col.reads = NULL,
149
+            col.sashimi = NULL,
150
+            col = .DEFAULT_SHADED_COL,
151
+            collapse = FALSE,
152
+            coverageHeight = 0.1,
153
+            fill.coverage = NULL,
154
+            fill.reads = NULL,
155
+            fill = "#BABABA",
156
+            fontface.mismatch = 2,
157
+            lty.coverage = NULL,
158
+            lty.gap = NULL,
159
+            lty.mates = NULL,
160
+            lty.deletion = NULL,
161
+            lty.insertion = NULL,
162
+            lty.mismatch = NULL,
163
+            lty.reads = NULL,
164
+            lty = 1,
165
+            lwd.coverage = NULL,
166
+            lwd.gap = NULL,
167
+            lwd.mates = NULL,
168
+            lwd.deletion = NULL,
169
+            lwd.insertion = NULL,
170
+            lwd.mismatch = NULL,
171
+            lwd.reads = NULL,
172
+            lwd.sashimiMax = 10,
173
+            lwd = 1,
174
+            max.height = 10,
175
+            min.height = 5,
176
+            minCoverageHeight = 50,
177
+            minSashimiHeight = 50,
178
+            noLetters = FALSE,
179
+            sashimiFilter = NULL,
180
+            sashimiFilterTolerance = 0L,
181
+            sashimiHeight = 0.1,
182
+            sashimiScore = 1,
183
+            sashimiStrand = "*",
184
+            sashimiTransformation = NULL,
185
+            showIndels = FALSE,
186
+            showMismatches = TRUE,
187
+            size = NULL,
188
+            transformation = NULL,
189
+            type = c("coverage", "pileup")
190
+        )
191
+    )
192
+)
193
+
194
+## Initialize ----------------------------------------------------------------
195
+
196
+#' @describeIn AlignmentsTrack-class Initialize.
197
+#' @export
198
+setMethod("initialize", "AlignmentsTrack", function(.Object, stackRanges = GRanges(), stacks = numeric(), sequences = DNAStringSet(),
199
+                                                    referenceSequence = NULL, ...) {
200
+    ## the diplay parameter defaults
201
+    .makeParMapping()
202
+    .Object <- .updatePars(.Object, "AlignmentsTrack")
203
+    .Object@stackRanges <- stackRanges
204
+    .Object <- callNextMethod(.Object, ...)
205
+    .Object@stacks <- stacks
206
+    .Object@sequences <- sequences
207
+    .Object@referenceSequence <- referenceSequence
208
+    return(.Object)
209
+})
210
+
211
+## ReferenceAlignmentsTrack Class --------------------------------------------
212
+
213
+#' @describeIn AlignmentsTrack-class The file-based version of the `AlignmentsTrack-class`.
214
+#' @exportClass ReferenceAlignmentsTrack
215
+setClass("ReferenceAlignmentsTrack", contains = c("AlignmentsTrack", "ReferenceTrack"))
216
+
217
+## Initialize ----------------------------------------------------------------
218
+
219
+## This just needs to set the appropriate slots that are being inherited
220
+## from ReferenceTrack because the multiple inheritance has some strange
221
+## features with regards to method selection
222
+
223
+#' @importClassesFrom Biostrings DNAStringSet
224
+#' @importFrom Biostrings DNAStringSet
225
+#' @describeIn AlignmentsTrack-class Initialize.
226
+#' @export
227
+setMethod("initialize", "ReferenceAlignmentsTrack", function(.Object, stream, reference, mapping = list(),
228
+                                                             args = list(), defaults = list(), stacks = numeric(),
229
+                                                             stackRanges = GRanges(), sequences = Biostrings::DNAStringSet(),
230
+                                                             referenceSequence = NULL, ...) {
231
+    .Object <- selectMethod("initialize", "ReferenceTrack")(.Object = .Object, reference = reference, stream = stream,
232
+        mapping = mapping, args = args, defaults = defaults)
233
+    .Object <- callNextMethod(.Object, ...)
234
+    .Object@referenceSequence <- referenceSequence
235
+    return(.Object)
236
+})
237
+
238
+## Constructor ---------------------------------------------------------------
239
+
240
+## Constructor
241
+#' @describeIn AlignmentsTrack-class Constructor for `AlignmentsTrack-class`.
242
+#' @export
243
+AlignmentsTrack <- function(range = NULL, start = NULL, end = NULL, width = NULL, strand, chromosome, genome,
244
+                            stacking = "squish", id, cigar, mapq, flag = scanBamFlag(isUnmappedQuery = FALSE), isize, groupid, status, md, seqs,
245
+                            name = "AlignmentsTrack", isPaired = TRUE, importFunction, referenceSequence, ...) {
246
+    ## Some defaults
247
+    if (missing(importFunction)) {
248
+        importFunction <- .import.bam.alignments
249
+    }
250
+    covars <- .getCovars(range)
251
+    isStream <- FALSE
252
+    if (!is.character(range)) {
253
+        n <- max(c(length(start), length(end), length(width)), nrow(covars))
254
+        id <- .covDefault(id, covars[["id"]], paste("read", seq_len(n), sep = "_"))
255
+        cigar <- .covDefault(cigar, covars[["cigar"]], paste(if (is(range, "GRangesOrIRanges")) width(range) else width, "M", sep = ""))
256
+        mapq <- .covDefault(mapq, covars[["mapq"]], rep(as.integer(NA), n))
257
+        flag <- .covDefault(flag, covars[["flag"]], rep(as.integer(NA), n))
258
+        isize <- .covDefault(isize, covars[["isize"]], rep(as.integer(NA), n))
259
+        groupid <- .covDefault(groupid, covars[["groupid"]], seq_len(n))
260
+        md <- .covDefault(md, covars[["md"]], rep(as.character(NA), n))
261
+        status <- .covDefault(status, covars[["status"]], ifelse(groupid %in% groupid[duplicated(groupid)], "mated", "unmated"))
262
+    }
263
+    ## Build a GRanges object from the inputs
264
+    .missingToNull(c(
265
+        "strand", "chromosome", "importFunction", "genome", "id", "cigar", "mapq", "flag", "isize", "groupid", "status",
266
+        "md", "seqs", "referenceSequence"
267
+    ))
268
+    args <- list(
269
+        id = id, cigar = cigar, mapq = mapq, flag = flag, isize = isize, groupid = groupid, status = status, strand = strand, md = md,
270
+        chromosome = chromosome, genome = genome
271
+    )
272
+    defs <- list(
273
+        strand = "*", chromosome = "chrNA", genome = NA, id = as.character(NA), cigar = as.character(NA), mapq = as.integer(NA),
274
+        flag = as.integer(NA), isize = as.integer(NA), groupid = as.character(NA), status = as.character(NA), md = as.character(NA)
275
+    )
276
+    range <- .buildRange(
277
+        range = range, start = start, end = end, width = width,
278
+        args = args, defaults = defs, chromosome = chromosome, trackType = "AlignmentsTrack",
279
+        importFun = importFunction, stream = TRUE, autodetect = TRUE
280
+    )
281
+    ## This is going to be a list if we have to stream data from a file, otherwise we can compute some additional values
282
+    if (is.list(range)) {
283
+        isStream <- TRUE
284
+        slist <- range
285
+        range <- GRanges()
286
+        stackRanges <- GRanges()
287
+        stacks <- NULL
288
+        seqs <- DNAStringSet()
289
+    } else {
290
+        if (is.null(seqs)) {
291
+            seqs <- DNAStringSet(vapply(width(range), function(x) paste(rep("N", x), collapse = ""), character(1)))
292
+        }
293
+        addArgs <- list(...)
294
+        if ("showIndels" %in% names(addArgs)) {
295
+            showIndels <- addArgs$showIndels
296
+        } else {
297
+            showIndels <- FALSE
298
+        }
299
+        tmp <- .computeAlignments(range, drop.D.ranges = showIndels)
300
+        range <- tmp$range
301
+        stackRanges <- tmp$stackRange
302
+        stacks <- tmp$stacks
303
+    }
304
+    ## If no chromosome was explicitly asked for we just take the first one in the GRanges object
305
+    if (missing(chromosome) || is.null(chromosome)) {
306
+        chromosome <- if (length(range) > 0) .chrName(as.character(seqnames(range)[1])) else "chrNA"
307
+    }
308
+    ## And finally the object instantiation
309
+    genome <- .getGenomeFromGRange(range, ifelse(is.null(genome), character(), genome[1]))
310
+    if (!isStream) {
311
+        return(new("AlignmentsTrack",
312
+            chromosome = chromosome[1], range = range, stacks = stacks,
313
+            name = name, genome = genome, stacking = stacking, stackRanges = stackRanges, sequences = seqs,
314
+            referenceSequence = referenceSequence, ...
315
+        ))
316
+    } else {
317
+        ## A bit hackish but for some functions we may want to know which track type we need but at the
318
+        ## same time we do not want to enforce this as an additional argument
319
+        e <- new.env()
320
+        e[["._trackType"]] <- "AlignmentsTrack"
321
+        e[["._isPaired"]] <- isPaired
322
+        e[["._flag"]] <- flag
323
+        environment(slist[["stream"]]) <- e
324
+        return(new("ReferenceAlignmentsTrack",
325
+            chromosome = chromosome[1], range = range, stackRanges = stackRanges,
326
+            name = name, genome = genome, stacking = stacking, stream = slist[["stream"]], reference = slist[["reference"]],
327
+            mapping = slist[["mapping"]], args = args, defaults = defs, stacks = stacks, referenceSequence = referenceSequence, ...
328
+        ))
329
+    }
330
+}
331
+
332
+## General accessors ---------------------------------------------------------
333
+
334
+#' @describeIn AlignmentsTrack-class Return all additional annotation
335
+#' information except for the genomic coordinates for the track items as
336
+#' a `data.frame`.
337
+#' @export
338
+setMethod("values", "AlignmentsTrack", function(x) .dpOrDefault(x, ".__coverage"))
339
+
340
+#' @describeIn AlignmentsTrack-class replace the value of the track's chromosome.
341
+#' This has to be a valid UCSC chromosome identifier or an integer or character
342
+#' scalar that can be reasonably coerced into one.
343
+#' @export
344
+setReplaceMethod("chromosome", "AlignmentsTrack", function(GdObject, value) {
345
+    GdObject <- callNextMethod()
346
+    if (!is.null(GdObject@referenceSequence)) {
347
+        chromosome(GdObject@referenceSequence) <- value[1]
348
+    }
349
+    return(GdObject)
350
+})
351
+
352
+
353
+
354
+## Annotation Accessors ------------------------------------------------------
355
+## Stacking ------------------------------------------------------------------
356
+
357
+#' @describeIn AlignmentsTrack-class return the stack indices for each track item.
358
+#' @export
359
+setMethod(
360
+    "stacks", "AlignmentsTrack",
361
+    function(GdObject) if (length(GdObject)) ranges(GdObject)$stack else 0
362
+)
363
+
364
+#' @describeIn AlignmentsTrack-class recompute the stacks based on the available
365
+#' space and on the object's track items and stacking settings.
366
+#' @export
367
+setMethod("setStacks", "AlignmentsTrack", function(GdObject, ...) {
368
+    if (length(GdObject)) {
369
+        bins <- if (!.needsStacking(GdObject)) rep(1, length(GdObject)) else disjointBins(GdObject@stackRanges)
370
+        GdObject@stacks <- bins
371
+        ranges(GdObject)$stack <- bins[match(ranges(GdObject)$groupid, names(GdObject@stackRanges))]
372
+    }
373
+    return(GdObject)
374
+})
375
+
376
+
377
+## Consolidate ---------------------------------------------------------------
378
+## Collapse  -----------------------------------------------------------------
379
+## Subset --------------------------------------------------------------------
380
+## AlignmentTracks can be subset by using the information in the stackRanges slot, but for the actual reads we need to make sure that
381
+## we keep all the bits that belong to a given group. We still want to record the requested ranges in the internal '.__plottingRange'
382
+## display parameter.
383
+
384
+#' @describeIn AlignmentsTrack-class Subset a `AlignmentsTrack` by coordinates
385
+#' and sort if necessary.
386
+#' @export
387
+setMethod("subset", signature(x = "AlignmentsTrack"), function(x, from = NULL, to = NULL, stacks = FALSE, use.defaults = TRUE, ...) {
388
+    ## Subset to a single chromosome first
389
+    lx <- length(x)
390
+    csel <- seqnames(x) != chromosome(x)
391
+    if (any(csel)) {
392
+        x <- x[!csel]
393
+    }
394
+    if (length(x)) {
395
+        ## Nothing to do if everything is within the range, otherwise we subset the stackRanges and keep all group items
396
+        ranges <- if (use.defaults) {
397
+            .defaultRange(x, from = from, to = to)
398
+        } else {
399
+            c(
400
+                from = ifelse(is.null(from), min(start(x@stackRanges)) - 1, from),
401
+                to = ifelse(is.null(to), max(end(x@stackRanges)) + 1, to)
402
+            )
403
+        }
404
+        displayPars(x) <- list(".__plottingRange" = ranges)
405
+        sr <- subsetByOverlaps(x@stackRanges, GRanges(seqnames = chromosome(x)[1], ranges = IRanges(start = ranges["from"], end = ranges["to"])))
406
+        if (length(sr) < length(x@stackRanges)) {
407
+            x@stackRanges <- sr
408
+            x@range <- x@range[x@range$groupid %in% names(sr)]
409
+            if (stacks) {
410
+                x <- setStacks(x)
411
+            }
412
+        }
413
+    }
414
+    return(x)
415
+})
416
+
417
+## ReferenceAlignmentsTracks need to stream the data from file and then pass the results on to the next method
418
+
419
+#' @describeIn AlignmentsTrack-class Subset a `ReferenceAlignmentsTrack` by coordinates
420
+#' and sort if necessary.
421
+#' @export
422
+setMethod("subset", signature(x = "ReferenceAlignmentsTrack"), function(x, from, to, chromosome, ...) {
423
+    ## We only need to reach out into the referenced file once if the range is already contained in the object
424
+    if (missing(from) || is.null(from) || missing(to) || is.null(to)) {
425
+        stop("Need both start and end location to subset a ReferenceAnnotationTrack")
426
+    }
427
+    if (missing(chromosome) || is.null(chromosome)) {
428
+        chromosome <- Gviz::chromosome(x)
429
+    } else {
430
+        chromosome <- .chrName(chromosome)
431
+    }
432
+    subRegion <- GRanges(seqnames = chromosome[1], ranges = IRanges(start = from, end = to))
433
+    oldRange <- .dpOrDefault(x, ".__plottingRange")
434
+    isIn <- length(ranges(x)) != 0 && !is.null(oldRange) && ranges(subRegion) %within% IRanges(oldRange["from"], oldRange["to"])
435
+    if (!isIn) {
436
+        cMap <- .resolveColMapping(x@stream(x@reference, subRegion), x@args, x@mapping)
437
+        seqs <- cMap$data$seq
438
+        cMap$data$seq <- NULL
439
+        range <- .computeAlignments(.buildRange(cMap$data, args = cMap$args, defaults = x@defaults, trackType = "AnnotationTrack"), drop.D.ranges = .dpOrDefault(x, "showIndels", FALSE))
440
+        ranges(x) <- range$range
441
+        x@stackRanges <- range$stackRanges
442
+        x@stacks <- range$stacks
443
+        x@sequences <- seqs
444
+    } else {
445
+        x@sequences <- subseq(x@sequences, start = from - (oldRange["from"] - 1), width = to - from + 1)
446
+    }
447
+    chromosome(x) <- chromosome[1]
448
+    return(callNextMethod(x = x, from = from, to = to, drop = FALSE, ...))
449
+})
450
+
451
+
452
+## Position ------------------------------------------------------------------
453
+## DrawGrid ------------------------------------------------------------------
454
+
455
+setMethod("drawGrid", signature(GdObject = "AlignmentsTrack"), function(GdObject, from, to) {
456
+    if (.dpOrDefault(GdObject, "grid", FALSE)) {
457
+        yvals <- values(GdObject)
458
+        ylim <- .dpOrDefault(GdObject, "ylim", if (!is.null(yvals) && length(yvals)) {
459
+            range(yvals, na.rm = TRUE, finite = TRUE)
460
+        } else {
461
+            c(-1, 1)
462
+        })
463
+        if (diff(ylim) == 0) {
464
+            ylim <- ylim + c(-1, 1)
465
+        }
466
+        yscale <- c(if (is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(ylim), max(ylim) + diff(range(ylim)) * 0.05)
467
+        covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", 0)
468
+        pushViewport(viewport(y = 1 - covHeight["npc"], height = covHeight["npc"], just = c(0.5, 0), yscale = yscale, clip = TRUE))
469
+        panel.grid(
470
+            v = 0, h = .dpOrDefault(GdObject, "h", -1),
471
+            col = .dpOrDefault(GdObject, "col.grid", "#e6e6e6"), lty = .dpOrDefault(GdObject, "lty.grid", 1),
472
+            lwd = .dpOrDefault(GdObject, "lwd.grid", 1)
473
+        )
474
+        popViewport(1)
475
+    }
476
+})
477
+
478
+## DrawAxis ------------------------------------------------------------------
479
+
480
+#' @describeIn AlignmentsTrack-class add a y-axis to the title panel of a track.
481
+#' @export
482
+setMethod("drawAxis", signature(GdObject = "AlignmentsTrack"), function(GdObject, ...) {
483
+    type <- match.arg(.dpOrDefault(GdObject, "type", .ALIGNMENT_TYPES), .ALIGNMENT_TYPES, several.ok = TRUE)
484
+    if ("coverage" %in% type) {
485
+        yvals <- values(GdObject)
486
+        ylim <- .dpOrDefault(GdObject, "ylim", if (!is.null(yvals) && length(yvals)) {
487
+            range(yvals, na.rm = TRUE, finite = TRUE)
488
+        } else {
489
+            c(-1, 1)
490
+        })
491
+        if (diff(ylim) == 0) {
492
+            ylim <- ylim + c(-1, 1)
493
+        }
494
+        hSpaceAvail <- vpLocation()$isize["width"] / 6
495
+        yscale <- c(if (is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(ylim), max(ylim) + diff(range(ylim)) * 0.05)
496
+        col <- .dpOrDefault(GdObject, "col.axis", "white")
497
+        acex <- .dpOrDefault(GdObject, "cex.axis")
498
+        acol <- .dpOrDefault(GdObject, "col.axis", "white")
499
+        at <- pretty(yscale)
500
+        at <- at[at >= sort(ylim)[1] & at <= sort(ylim)[2]]
501
+        covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", c(npc = 0, points = 0))
502
+        covSpace <- .dpOrDefault(GdObject, ".__coverageSpace", 0)
503
+        pushViewport(viewport(y = 1 - (covHeight["npc"] + covSpace), height = covHeight["npc"], just = c(0.5, 0)))
504
+        if (is.null(acex)) {
505
+            vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(at), "inches"))) * length(at) * 1.5
506
+            hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(at), "inches")))
507
+            vSpaceAvail <- abs(diff(range(at))) / abs(diff(yscale)) * vpLocation()$isize["height"]
508
+            acex <- max(0.6, min(vSpaceAvail / vSpaceNeeded, hSpaceAvail / hSpaceNeeded))
509
+        }
510
+        vpTitleAxis <- viewport(x = 0.95, width = 0.2, yscale = yscale, just = 0)
511
+        pushViewport(vpTitleAxis)
512
+        suppressWarnings(grid.yaxis(gp = gpar(col = acol, cex = acex), at = at))
513
+        grid.lines(x = c(0, 0), y = ylim, gp = gpar(col = acol), default.units = "native")
514
+        popViewport(2)
515
+    } else {
516
+        covHeight <- c(npc = 0, points = 0)
517
+        covSpace <- 0
518
+    }
519
+    if ("sashimi" %in% type) {
520
+        sash <- .dpOrDefault(GdObject, ".__sashimi", list(x = numeric(), y = numeric(), id = integer(), score = numeric()))
521
+        yscale <- if (length(sash$y)) c(-(max(sash$y) + diff(range(sash$y)) * 0.05), 0) else c(-1, 0)
522
+        ylim <- if (length(sash$y)) c(-max(sash$y), yscale[1] + max(sash$y)) else c(-1, 0)
523
+        hSpaceAvail <- vpLocation()$isize["width"] / 6
524
+        col <- .dpOrDefault(GdObject, "col.axis", "white")
525
+        acex <- .dpOrDefault(GdObject, "cex.axis")
526
+        acol <- .dpOrDefault(GdObject, "col.axis", "white")
527
+        labs <- if (length(sash$score)) pretty(c(1, sash$score)) else pretty(c(1, .dpOrDefault(GdObject, ".__sashimiScore", 10)))
528
+        at <- seq(ylim[1], ylim[2], length.out = length(labs))
529
+        sashHeight <- .dpOrDefault(GdObject, ".__sashimiHeight", c(npc = 0, points = 0))
530
+        sashSpace <- .dpOrDefault(GdObject, ".__sashimiSpace", 0)
531
+        pushViewport(viewport(
532
+            y = 1 - (sashHeight["npc"] + sashSpace + covHeight["npc"] + covSpace),
533
+            height = sashHeight["npc"], just = c(0.5, 0)
534
+        ))
535
+        if (is.null(acex)) {
536
+            vSpaceNeeded <- max(as.numeric(convertWidth(stringHeight(labs), "inches"))) * length(at) * 1.5
537
+            hSpaceNeeded <- max(as.numeric(convertWidth(stringWidth(labs), "inches")))
538
+            vSpaceAvail <- abs(diff(range(at))) / abs(diff(yscale)) * vpLocation()$isize["height"]
539
+            acex <- max(0.6, min(vSpaceAvail / vSpaceNeeded, hSpaceAvail / hSpaceNeeded))
540
+        }
541
+        vpTitleAxis <- viewport(x = 0.75, width = 0.2, yscale = yscale, just = 0)
542
+        pushViewport(vpTitleAxis)
543
+        suppressWarnings(grid.yaxis(gp = gpar(col = acol, cex = acex), at = at, label = labs))
544
+        grid.polygon(x = c(0, 0, 1), y = c(ylim[1], ylim[2], ylim[2]), default.units = "native", gp = gpar(col = acol, fill = acol))
545
+        popViewport(2)
546
+    }
547
+    if (.dpOrDefault(GdObject, ".__isCropped", FALSE)) {
548
+        vspacing <- as.numeric(convertHeight(unit(2, "points"), "npc"))
549
+        vsize <- as.numeric(convertHeight(unit(4, "points"), "npc"))
550
+        pushViewport(viewport(height = vsize, y = vspacing, just = c(0.5, 0)))
551
+        .moreInd(direction = "down", lwd = 2)
552
+        popViewport(1)
553
+    }
554
+})
555
+
556
+## DrawGD --------------------------------------------------------------------
557
+
558
+#' @describeIn AlignmentsTrack-class plot the object to a graphics device.
559
+#' The return value of this method is the input object, potentially updated
560
+#' during the plotting operation. Internally, there are two modes in which the
561
+#' method can be called. Either in 'prepare' mode, in which case no plotting is
562
+#' done but the object is preprocessed based on the available space, or in
563
+#' 'plotting' mode, in which case the actual graphical output is created.
564
+#' Since subsetting of the object can be potentially costly, this can be
565
+#' switched off in case subsetting has already been performed before or
566
+#' is not necessary.
567
+#'
568
+#' @importFrom GenomicAlignments cigarRangesAlongReferenceSpace
569
+#' @export
570
+setMethod("drawGD", signature("AlignmentsTrack"), function(GdObject, minBase, maxBase, prepare = FALSE, subset = TRUE, ...) {
571
+    debug <- .dpOrDefault(GdObject, "debug", FALSE)
572
+    imageMap(GdObject) <- NULL
573
+    if (!length(GdObject)) {
574
+        return(invisible(GdObject))
575
+    }
576
+    rev <- .dpOrDefault(GdObject, "reverseStrand", FALSE)
577
+    revs <- !.dpOrDefault(GdObject, "reverseStacking", FALSE)
578
+    vSpacing <- as.numeric(convertHeight(unit(3, "points"), "npc"))
579
+    type <- match.arg(.dpOrDefault(GdObject, "type", .ALIGNMENT_TYPES), .ALIGNMENT_TYPES, several.ok = TRUE)
580
+    ylim <- .dpOrDefault(GdObject, "ylim")
581
+    ## In prepare mode we need to make sure that the stacking information is updated from the optional display parameter (by calling
582
+    ## the StackedTrack drawGD method), and compute the coverage vector which will be needed for the axis
583
+    if (prepare) {
584
+        if ((is.logical(debug) && debug) || "prepare" %in% debug) {
585
+            browser()
586
+        }
587
+        GdObject <- callNextMethod(GdObject, ...)
588
+        ## The mismatched bases need to be extracted from the read sequences and the reference sequence
589
+        if (.dpOrDefault(GdObject, "showMismatches", TRUE) && !is.null(GdObject@referenceSequence)) {
590
+            mm <- .findMismatches(GdObject)
591
+            if (nrow(mm)) {
592
+                displayPars(GdObject) <- list(".__mismatches" = mm)
593
+            }
594
+        }
595
+        ## The coverage calculation and the height of the coverage section
596
+        if ("coverage" %in% type) {
597
+            covSpace <- as.numeric(convertHeight(unit(5, "points"), "npc"))
598
+            if ("pileup" %in% type) {
599
+                covHeight <- .dpOrDefault(GdObject, "coverageHeight", 0.1)
600
+                if (covHeight > 0 && covHeight < 1) {
601
+                    covHeight <- as.numeric(convertHeight(unit(covHeight, "npc"), "points"))
602
+                }
603
+                covHeight <- max(.dpOrDefault(GdObject, "minCoverageHeight", 50), covHeight)
604
+                covHeight <- c(points = covHeight, npc = as.numeric(convertHeight(unit(covHeight, "points"), "npc")))
605
+            } else if ("sashimi" %in% type) {
606
+                covHeight <- c(npc = 0.5 - (covSpace * 2), points = 1)
607
+            } else {
608
+                covHeight <- c(npc = 1 - (covSpace * 2), points = 1)
609
+            }
610
+            coverage <- coverage(range(GdObject), width = maxBase)
611
+            ## apply data transformation if one is set up
612
+            trans <- .dpOrDefault(GdObject, "transformation")
613
+            if (is.list(trans)) {
614
+                trans <- trans[[1]]
615
+            }
616
+            if (!is.null(trans)) {
617
+                if (!is.function(trans) || length(formals(trans)) != 1L) {
618
+                    stop("Display parameter 'transformation' must be a function with a single argument")
619
+                }
620
+                test <- trans(coverage)
621
+                if (!is(test, "Rle") || length(test) != length(coverage)) {
622
+                    stop(
623
+                        "The function in display parameter 'transformation' results in invalid output.\n",
624
+                        "It has to return a numeric matrix with the same dimensions as the input data."
625
+                    )
626
+                }
627
+                coverage <- test
628
+            }
629
+            displayPars(GdObject) <- list(
630
+                ".__coverage" = coverage,
631
+                ".__coverageHeight" = covHeight,
632
+                ".__coverageSpace" = covSpace
633
+            )
634
+        } else {
635
+            covHeight <- c(npc = 0, points = 0)
636
+            covSpace <- 0
637
+        }
638
+        ## The sashimi calculation and the height of the sashimi section
639
+        if ("sashimi" %in% type) {
640
+            sashSpace <- as.numeric(convertHeight(unit(5, "points"), "npc"))
641
+            if ("pileup" %in% type) {
642
+                sashHeight <- .dpOrDefault(GdObject, "sashimiHeight", 0.1)
643
+                if (sashHeight > 0 && sashHeight < 1) {
644
+                    sashHeight <- as.numeric(convertHeight(unit(sashHeight, "npc"), "points"))
645
+                }
646
+                sashHeight <- max(.dpOrDefault(GdObject, "minSashimiHeight", 50), sashHeight)
647
+                sashHeight <- c(points = sashHeight, npc = as.numeric(convertHeight(unit(sashHeight, "points"), "npc")))
648
+            } else if ("coverage" %in% type) {
649
+                sashHeight <- c(npc = 0.5 - (sashSpace * 2), points = 1)
650
+            } else {
651
+                sashHeight <- c(npc = 1 - (sashSpace * 2), points = 1)
652
+            }
653
+            sashScore <- .dpOrDefault(GdObject, "sashimiScore", 1L)
654
+            sashLwdMax <- .dpOrDefault(GdObject, "lwd.sashimiMax", 10)
655
+            sashStrand <- .dpOrDefault(GdObject, "sashimiStrand", "*")
656
+            sashFilter <- .dpOrDefault(GdObject, "sashimiFilter", NULL)
657
+            sashFilterTolerance <- .dpOrDefault(GdObject, "sashimiFilterTolerance", 0L)
658
+            sashNumbers <- .dpOrDefault(GdObject, "sashimiNumbers", FALSE)
659
+            sash <- .dpOrDefault(GdObject, "sashimiJunctions", NULL)
660
+            if (is.null(sash)) {
661
+                sash <- .create.summarizedJunctions.for.sashimi.junctions(ranges(GdObject))
662
+            } else {
663
+                if (!is(sash, "GRanges")) {
664
+                    stop("\"sashimiJunctions\" object must be of \"GRanges\" class!")
665
+                }
666
+                sashMcolName <- if (sashStrand == "+") "plus_score" else if (sashStrand == "-") "minus_score" else "score"
667
+                if (sum(colnames(mcols(sash)) == sashMcolName) != 1) {
668
+                    stop(sprintf("\"mcols\" of \"sashimiJunctions\" object must contain column named \"%s\",\n which matches the specified (%s) \"sashimiStrand\"!", sashMcolName, sashStrand))
669
+                }
670
+            }
671
+            sashTransform <- .dpOrDefault(GdObject, c("sashimiTransformation", "transformation"))
672
+            sash <- .convert.summarizedJunctions.to.sashimi.junctions(
673
+                juns = sash,
674
+                score = sashScore,
675
+                lwd.max = sashLwdMax,
676
+                strand = sashStrand,
677
+                filter = sashFilter,
678
+                filterTolerance = sashFilterTolerance,
679
+                trans = sashTransform
680
+            )
681
+            displayPars(GdObject) <- list(
682
+                ".__sashimi" = sash,
683
+                ".__sashimiHeight" = sashHeight,
684
+                ".__sashimiSpace" = sashSpace,
685
+                ".__sashimiNumbers" = sashNumbers
686
+            )
687
+        } else {
688
+            sashHeight <- c(npc = 0, points = 0)
689
+            sashSpace <- 0
690
+        }
691
+        if ("pileup" %in% type) {
692
+            ## If there are more bins than we can plot we reduce the number until they fit
693
+            pushViewport(viewport(height = 1 - (covHeight["npc"] + covSpace + sashHeight["npc"] + sashSpace) - vSpacing * 2, y = vSpacing, just = c(0.5, 0)))
694
+            bins <- ranges(GdObject)$stack
695
+            curVp <- vpLocation()
696
+            mih <- min(curVp$size["height"], .dpOrDefault(GdObject, c("min.height", "max.height"), 3))
697
+            mah <- min(curVp$size["height"], max(mih, .dpOrDefault(GdObject, c("max.height", "min.height"), 8)))
698
+            bins <- stacks(GdObject)
699
+            if (curVp$size["height"] / max(bins) < mih) {
700
+                maxStack <- curVp$size["height"] %/% mih
701
+                sel <- if (revs) bins <= maxStack else bins > max(bins) - maxStack
702
+                ranges(GdObject) <- ranges(GdObject)[sel]
703
+                displayPars(GdObject) <- list(".__isCropped" = TRUE)
704
+                bins <- stacks(GdObject)
705
+            }
706
+            yrange <- range(bins) + c(-0.5, 0.5)
707
+            add <- max(0, (curVp$size["height"] %/% mah) - max(bins))
708
+            yrange <- if (revs) c(yrange[1], yrange[2] + add) else c(yrange[1] - add, yrange[2])
709
+            displayPars(GdObject) <- list(".__yrange" = yrange)
710
+            popViewport(1)
711
+        }
712
+        return(invisible(GdObject))
713
+    }
714
+    if ((is.logical(debug) && debug) || "draw" %in% debug) {
715
+        browser()
716
+    }
717
+    mm <- .dpOrDefault(GdObject, ".__mismatches")
718
+    ## The coverage plot first
719
+    xscale <- if (!rev) c(minBase, maxBase) else c(maxBase, minBase)
720
+    readInfo <- ranges(GdObject)
721
+    covHeight <- .dpOrDefault(GdObject, ".__coverageHeight", c(npc = 0, points = 0))
722
+    covSpace <- .dpOrDefault(GdObject, ".__coverageSpace", 0)
723
+    if ("coverage" %in% type) {
724
+        cov <- .dpOrDefault(GdObject, ".__coverage", Rle(lengths = maxBase, values = as.integer(0)))
725
+        yscale <- c(if (is.null(.dpOrDefault(GdObject, "transformation"))) 0 else min(cov), max(cov) + diff(range(cov)) * 0.05)
726
+        if (!is.null(ylim)) {
727
+            yscale <- ylim
728
+        }
729
+        vp <- viewport(
730
+            height = covHeight["npc"], y = 1 - (covHeight["npc"] + covSpace), just = c(0.5, 0), xscale = xscale,
731
+            yscale = yscale, clip = TRUE
732
+        )
733
+        pushViewport(vp)
734
+        res <- .pxResolution(coord = "x")
735
+        gp <- gpar(
736
+            col = .dpOrDefault(GdObject, c("col.coverage", "col"), .DEFAULT_SHADED_COL),
737
+            fill = .dpOrDefault(GdObject, c("fill.coverage", "fill"), "#BABABA"),
738
+            lwd = .dpOrDefault(GdObject, c("lwd.coverage", "lwd"), 1),
739
+            lty = .dpOrDefault(GdObject, c("lty.coverage", "lty"), 1),
740
+            alpha = .alpha(GdObject)
741
+        )
742
+        ## We can compact this if the resolution is not sufficient to speed up the drawing.
743
+        mminBase <- max(1, minBase)
744
+        if (res > 2) {
745
+            brks <- ceiling((maxBase - mminBase) / res)
746
+            x <- seq(mminBase, maxBase, len = brks)
747
+            y <- tapply(as.integer(cov[mminBase:maxBase]), cut(mminBase:maxBase, breaks = brks), mean)
748
+        } else {
749
+            x <- mminBase:maxBase
750
+            y <- as.integer(cov[x])
751
+        }
752
+        grid.polygon(c(minBase - max(1, res), 0, x, maxBase + max(1, res)), c(0, 0, y, 0), default.units = "native", gp = gp)
753
+        grid.lines(y = c(0, 0), gp = gpar(col = gp$col, alpha = gp$alpha))
754
+        if (!is.null(mm)) {
755
+            fcol <- .dpOrDefault(GdObject@referenceSequence, "fontcolor", getBioColor("DNA_BASES_N"))
756
+            vpos <- tapply(as.character(mm$base[mm$base != "."]), mm$position[mm$base != "."], table, simplify = FALSE)
757
+            x <- rep(as.integer(names(vpos)), listLen(vpos))
758
+            y <- unlist(lapply(vpos, cumsum), use.names = FALSE)
759
+            col <- fcol[unlist(lapply(vpos, names), use.names = FALSE)]
760
+            grid.rect(
761
+                x = x, y = y, height = unlist(vpos, use.names = FALSE), width = 1, default.units = "native", just = c(0, 1),
762
+                gp = gpar(col = "transparent", fill = col)
763
+            )
764
+        }
765
+        popViewport(1)
766
+        twoPx <- 2 * as.numeric(convertHeight(unit(1, "points"), "npc"))
767
+        vp <- viewport(height = twoPx, y = 1 - (covHeight["npc"] + covSpace + twoPx), just = c(0.5, 0))
768
+        pushViewport(vp)
769
+        grid.rect(gp = gpar(
770
+            fill = .dpOrDefault(GdObject, "background.title"), col = "transparent",
771
+            alpha = .dpOrDefault(GdObject, "alpha.title")
772
+        ), width = 2)
773
+        popViewport(1)
774
+    }
775
+    ## The sashimi plot as second
776
+    xscale <- if (!rev) c(minBase, maxBase) else c(maxBase, minBase)
777
+    sashHeight <- .dpOrDefault(GdObject, ".__sashimiHeight", c(npc = 0, points = 0))
778
+    sashSpace <- .dpOrDefault(GdObject, ".__sashimiSpace", 0)
779
+    if ("sashimi" %in% type) {
780
+        sash <- .dpOrDefault(GdObject, ".__sashimi", list(x = numeric(), y = numeric(), id = integer(), score = numeric(), scaled = numeric()))
781
+        sashNumbers <- .dpOrDefault(GdObject, ".__sashimiNumbers", FALSE)
782
+        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
783
+        vp <- viewport(
784
+            height = sashHeight["npc"], y = 1 - (covHeight["npc"] + covSpace + sashHeight["npc"] + sashSpace), just = c(0.5, 0),
785
+            xscale = xscale, yscale = yscale, clip = TRUE
786
+        )
787
+        pushViewport(vp)
788
+        gp <- gpar(
789
+            col = .dpOrDefault(GdObject, c("col.sashimi", "col"), .DEFAULT_SHADED_COL),
790
+            fill = .dpOrDefault(GdObject, c("fill.sashimi", "fill"), "#FFFFFF"),
791
+            lwd = .dpOrDefault(GdObject, c("lwd.sashimi", "lwd"), 1),
792
+            lty = .dpOrDefault(GdObject, c("lty.sashimi", "lty"), 1),
793
+            alpha = .alpha(GdObject)
794
+        )
795
+        if (length(sash$x)) {
796
+            grid.xspline(sash$x, -sash$y,
797
+                id = sash$id, shape = -1, open = TRUE,
798
+                default.units = "native", gp = gpar(col = gp$col, lwd = sash$scaled)
799
+            )
800
+            ## print the number of reads together with the connecting lines (currently no scaling/resolution)
801
+            if (sashNumbers) {
802
+                grid.rect(sash$x[c(FALSE, TRUE, FALSE)], -sash$y[c(FALSE, TRUE, FALSE)],
803
+                    width = convertUnit(stringWidth(sash$score) * 1.5, "inches"),
804
+                    height = convertUnit(stringHeight(sash$score) * 1.5, "inches"),
805
+                    default.units = "native", gp = gpar(col = gp$col, fill = gp$fill)
806
+                )
807
+                grid.text(
808
+                    label = sash$score, sash$x[c(FALSE, TRUE, FALSE)], -sash$y[c(FALSE, TRUE, FALSE)],
809
+                    default.units = "native", gp = gpar(col = gp$col)
810
+                )
811
+            }
812
+        }
813
+        popViewport(1)
814
+    }
815
+    ## Now the pileups
816
+    if ("pileup" %in% type) {
817
+        cex <- max(0.3, .dpOrDefault(GdObject, c("cex.mismatch", "cex"), 0.7))
818
+        pushViewport(viewport(
819
+            height = 1 - (covHeight["npc"] + covSpace + sashHeight["npc"] + sashSpace) - vSpacing * 4, y = vSpacing, just = c(0.5, 0),
820
+            gp = .fontGp(GdObject, "mismatch", cex = cex)
821
+        ))
822
+        bins <- stacks(GdObject)
823
+        stacks <- max(bins)
824
+        yscale <- .dpOrDefault(GdObject, ".__yrange")
825
+        if (revs) {
826
+            yscale <- rev(yscale)
827
+        }
828
+        pushViewport(dataViewport(xscale = xscale, extension = 0, yscale = yscale, clip = TRUE))
829
+        ## Figuring out resolution and the necessary plotting details
830
+        res <- .pxResolution(coord = "x")
831
+        ylim <- c(0, 1)
832
+        h <- diff(ylim)
833
+        middle <- mean(ylim)
834
+        sh <- max(0, min(h, .dpOrDefault(GdObject, "stackHeight", 0.75))) / 2
835
+        boxOnly <- res > 10
836
+        if (boxOnly) {
837
+            x <- c(start(readInfo), rep(end(readInfo) + 1, 2), start(readInfo))
838
+            y <- c(rep(readInfo$stack + sh, 2), rep(readInfo$stack - sh, 2))
839
+            id <- rep(readInfo$uid, 4)
840
+        } else {
841
+            ## We first precompute the coordinates for all the arrow polygons
842
+            uid2strand <- setNames(as.character(readInfo$readStrand), as.character(readInfo$uid))
843
+            arrowMap <- unlist(setNames(lapply(split(as.character(readInfo$uid), readInfo$entityId), function(x) {
844
+                str <- uid2strand[x][1]
845
+                setNames(str, if (str == "+") tail(x, 1) else head(x, 1))
846
+            }), NULL))
847
+            readInfo$arrow <- as.character(NA)
848
+            readInfo$arrow[match(names(arrowMap), readInfo$uid)] <- arrowMap
849
+            ## The parts that don't need arrow heads
850
+            sel <- is.na(readInfo$arrow) | readInfo$arrow == "*"
851
+            x <- c(start(readInfo)[sel], rep(end(readInfo)[sel] + 1, 2), start(readInfo)[sel])
852
+            y <- c(rep(readInfo$stack[sel] + sh, 2), rep(readInfo$stack[sel] - sh, 2))
853
+            id <- rep(readInfo$uid[sel], 4)
854
+            ## The arrow heads facing right
855
+            w <- Gviz:::.pxResolution(coord = "x", 5)
856
+            sel <- readInfo$arrow == "+"
857
+            ah <- pmax(start(readInfo)[sel], end(readInfo)[sel] + 1 - w)
858
+            x <- c(x, start(readInfo)[sel], ah, end(readInfo)[sel] + 1, ah, start(readInfo)[sel])
859
+            y <- c(y, rep(readInfo$stack[sel] + sh, 2), readInfo$stack[sel], rep(readInfo$stack[sel] - sh, 2))
860
+            id <- c(id, rep(readInfo$uid[sel], 5))
861
+            ## The arrow heads facing left
862
+            sel <- readInfo$arrow == "-"
863
+            ah <- pmin(end(readInfo)[sel], start(readInfo)[sel] + w)
864
+            x <- c(x, start(readInfo)[sel], ah, rep(end(readInfo)[sel] + 1, 2), ah)
865
+            y <- c(y, readInfo$stack[sel], rep(readInfo$stack[sel] + sh, 2), rep(readInfo$stack[sel] - sh, 2))
866
+            id <- c(id, rep(readInfo$uid[sel], 5))
867
+        }
868
+        nn <- length(unique(readInfo$uid))
869
+        gps <- data.frame(
870
+            col = rep(.dpOrDefault(GdObject, c("col.reads", "col"), .DEFAULT_SHADED_COL), nn),
871
+            fill = rep(.dpOrDefault(GdObject, c("fill.reads", "fill"), .DEFAULT_BRIGHT_SHADED_COL), nn),
872
+            lwd = rep(.dpOrDefault(GdObject, c("lwd.reads", "lwd"), 1), nn),
873
+            lty = rep(.dpOrDefault(GdObject, c("lty.reads", "lty"), 1), nn),
874
+            alpha = rep(.alpha(GdObject, "reads"), nn), stringsAsFactors = FALSE
875
+        )
876
+        ## Finally we draw reads (we need to draw in two steps because of the different alpha levels, reads vs. mismatches)
877
+        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")
878
+        ## Now the coordinates for the connecting lines
879
+        lineCoords <- NULL
880
+        if (anyDuplicated(readInfo$entityId) != 0) {
881
+            stTmp <- split(readInfo, readInfo$entityId)
882
+            mateRanges <- unlist(range(stTmp))
883
+            mateGaps <- gaps(GRanges(ranges = ranges(readInfo), seqnames = readInfo$entityId))
884
+            rmap <- mateRanges[as.character(seqnames(mateGaps))]
885
+            mateGaps <- mateGaps[start(rmap) <= start(mateGaps) & end(rmap) >= end(mateGaps)]
886
+            gy <- readInfo$stack[match(as.character(seqnames(mateGaps)), readInfo$entityId)]
887
+            lineCoords <- data.frame(
888
+                x1 = start(mateGaps), y1 = gy, x2 = end(mateGaps) + 1, y2 = gy,
889
+                col = .dpOrDefault(GdObject, c("col.gap", "col"), .DEFAULT_SHADED_COL),
890
+                lwd = .dpOrDefault(GdObject, c("lwd.gap", "lwd"), 1),
891
+                lty = .dpOrDefault(GdObject, c("lty.gap", "lty"), 1),
892
+                alpha = .alpha(GdObject, "gap"), stringsAsFactors = FALSE
893
+            )
894
+            lineCoords <- lineCoords[!duplicated(lineCoords), ]
895
+        } else {
896
+            mateRanges <- setNames(readInfo, readInfo$entityId)
897
+        }
898
+        if (any(readInfo$status != "unmated") && anyDuplicated(readInfo$groupid) != 0) {
899
+            pairGaps <- gaps(GRanges(ranges = ranges(mateRanges), seqnames = readInfo$groupid[match(names(mateRanges), readInfo$entityId)]))
900
+            rmap <- GdObject@stackRanges[as.character(seqnames(pairGaps))]
901
+            pairGaps <- pairGaps[start(rmap) <= start(pairGaps) & end(rmap) >= end(pairGaps)]
902
+            gy <- readInfo$stack[match(as.character(seqnames(pairGaps)), readInfo$groupid)]
903
+            if (length(pairGaps)) {
904
+                pairsCoords <- data.frame(
905
+                    x1 = start(pairGaps) - 1, y1 = gy, x2 = end(pairGaps) + 1, y2 = gy,
906
+                    col = .dpOrDefault(GdObject, c("col.mates", "col"), .DEFAULT_BRIGHT_SHADED_COL),
907
+                    lwd = .dpOrDefault(GdObject, c("lwd.mates", "lwd"), 1),
908
+                    lty = .dpOrDefault(GdObject, c("lty.mates", "lty"), 1),
909
+                    alpha = .alpha(GdObject, "mates"), stringsAsFactors = FALSE
910
+                )
911
+            } else {
912
+                pairsCoords <- NULL
913
+            }
914
+            lineCoords <- rbind(lineCoords, pairsCoords[!duplicated(pairsCoords), ])
915
+        }
916
+        ## Consider the indels if needed
917
+        ## - deletion as lines
918
+        ## - insertions as vertical bars
919
+        showIndels <- .dpOrDefault(GdObject, "showIndels", FALSE)
920
+        delCoords <- NULL
921
+        insCoords <- NULL
922
+        if (showIndels) {
923
+            cigarTmp <- DataFrame(cigar = readInfo$cigar, start = start(readInfo), entityId = readInfo$entityId, groupId = readInfo$groupid)
924
+            cigarTmp <- cigarTmp[order(cigarTmp$entityId, cigarTmp$start), ]
925
+            cigarTmp <- cigarTmp[!duplicated(cigarTmp$entityId), ]
926
+            delGaps <- unlist(cigarRangesAlongReferenceSpace(cigarTmp$cigar, pos = cigarTmp$start, ops = "D", f = as.factor(cigarTmp$entityId)))
927
+            gy <- readInfo$stack[match(names(delGaps), readInfo$entityId)]
928
+            if (length(delGaps)) {
929
+                delCoords <- data.frame(
930
+                    x1 = start(delGaps) + 1, y1 = gy, x2 = end(delGaps) + 1, y2 = gy,
931
+                    col = .dpOrDefault(GdObject, c("col.deletion", "col"), .DEFAULT_BRIGHT_SHADED_COL),
932
+                    lwd = .dpOrDefault(GdObject, c("lwd.deletion", "lwd"), 1),
933
+                    lty = .dpOrDefault(GdObject, c("lty.deletion", "lty"), 1),
934
+                    alpha = .alpha(GdObject, "deletions"), stringsAsFactors = FALSE
935
+                )
936
+                lineCoords <- rbind(delCoords, lineCoords)
937
+                lineCoords <- lineCoords[!duplicated(lineCoords[, c("x1", "y1", "x2", "y2")]), ]
938
+            }
939
+            insGaps <- unlist(cigarRangesAlongReferenceSpace(cigarTmp$cigar, pos = cigarTmp$start, ops = "I", f = as.factor(cigarTmp$entityId)))
940
+            gy <- readInfo$stack[match(names(insGaps), readInfo$entityId)]
941
+            if (length(insGaps)) {
942
+                ## should both x coordinates be equal to start
943
+                insCoords <- data.frame(
944
+                    x1 = start(insGaps), y1 = gy - sh, x2 = start(insGaps), y2 = gy + sh,
945
+                    col = .dpOrDefault(GdObject, c("col.insertion", "col"), .DEFAULT_BRIGHT_SHADED_COL),
946
+                    lwd = .dpOrDefault(GdObject, c("lwd.insertion", "lwd"), 1),
947
+                    lty = .dpOrDefault(GdObject, c("lty.insertion", "lty"), 1),
948
+                    alpha = .alpha(GdObject, "insertions"), stringsAsFactors = FALSE
949
+                )
950
+            }
951
+        }
952
+        ## The mismatch information on the reads if needed
953
+        mmLetters <- NULL
954
+        if (!is.null(mm)) {
955
+            fcol <- .dpOrDefault(GdObject@referenceSequence, "fontcolor", getBioColor("DNA_BASES_N"))
956
+            ccol <- ifelse(rgb2hsv(col2rgb(fcol))["s", ] < 0.5, "black", "white")
957
+            vpl <- vpLocation()
958
+            lwidth <- max(as.numeric(convertUnit(stringWidth(DNA_ALPHABET), "inches"))) * 0.9337632
959
+            lheight <- max(as.numeric(convertUnit(stringHeight(DNA_ALPHABET), "inches")))
960
+            perLetterW <- vpl$isize["width"] / (maxBase - minBase + 1)
961
+            perLetterH <- vpl$isize["height"] / abs(diff(current.viewport()$yscale))
962
+            res <- .pxResolution(coord = "x")
963
+            mw <- res * .dpOrDefault(GdObject, "min.width", 1)
964
+            mwy <- max(1, mw)
965
+            if (nrow(mm)) {
966
+                x <- c(mm$position + rep(c(0, mwy, mwy, 0), each = nrow(mm)))
967
+                y <- c(rep(mm$stack - sh, 2), rep(mm$stack + sh, 2))
968
+                id <- c(rep(seq(max(id, na.rm = TRUE) + 1, len = nrow(mm)), 4))
969
+                gps <- data.frame(
970
+                    col = rep(if (!(lwidth < perLetterW && lheight < perLetterH)) {
971
+                        "transparent"
972
+                    } else {
973
+                        .dpOrDefault(GdObject, "col.mismatch", .DEFAULT_SHADED_COL)
974
+                    }, nrow(mm)),
975
+                    fill = rep(fcol[as.character(mm$base)]),
976
+                    lwd = rep(.dpOrDefault(GdObject, "lwd.mismatch", 1), nrow(mm)),
977
+                    lty = rep(.dpOrDefault(GdObject, "lty.mismatch", 1), nrow(mm)),
978
+                    alpha = rep(.dpOrDefault(GdObject, "alpha.mismatch", 1), nrow(mm)),
979
+                    stringsAsFactors = FALSE
980
+                )
981
+                ## Finally we draw mm (we need to draw in two steps because of the different alpha levels, reads vs. mismatches)
982
+                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")
983
+                if (!.dpOrDefault(GdObject, "noLetters", FALSE) && lwidth < perLetterW && lheight < perLetterH) {
984
+                    mmLetters <- data.frame(x = mm$position + 0.5, y = mm$stack, label = mm$base, col = ccol[mm$base], stringsAsFactors = FALSE)
985
+                }
986
+            }
987
+        }
988
+        if (!is.null(lineCoords)) {
989
+            grid.segments(lineCoords$x1, lineCoords$y1, lineCoords$x2, lineCoords$y2,
990
+                gp = gpar(
991
+                    col = lineCoords$col, alpha = unique(lineCoords$alpha),
992
+                    lwd = lineCoords$lwd, lty = lineCoords$lty
993
+                ),
994
+                default.units = "native"
995
+            )
996
+        }
997
+        if (!is.null(insCoords)) {
998
+            grid.segments(insCoords$x1, insCoords$y1, insCoords$x2, insCoords$y2,
999
+                gp = gpar(
1000
+                    col = insCoords$col, alpha = unique(insCoords$alpha),
1001
+                    lwd = insCoords$lwd, lty = insCoords$lty
1002
+                ),
1003
+                default.units = "native"
1004
+            )
1005
+        }
1006
+        if (!is.null(mmLetters)) {
1007
+            grid.text(
1008
+                x = mmLetters$x, y = mmLetters$y, label = mmLetters$label,
1009
+                gp = gpar(col = mmLetters$col), default.units = "native"
1010
+            )
1011
+        }
1012
+        popViewport(2)
1013
+    }
1014
+    ## Eventually we set up the image map
1015
+    ## imageMap(GdObject) <- im
1016
+    return(invisible(GdObject))
1017
+})
1018
+
1019
+## SetAs ---------------------------------------------------------------------
1020
+## Show ----------------------------------------------------------------------
1021
+
1022
+#' @describeIn AlignmentsTrack-class Show method.
1023
+#' @export
1024
+setMethod(
1025
+    "show", signature(object = "AlignmentsTrack"),
1026
+    function(object) {
1027
+        cat(sprintf(
1028
+            paste("AlignmentsTrack track '%s' \n",
1029
+                "| genome: %s\n",
1030
+                "| active chromosome: %s\n",
1031
+                "| containing %i read%s\n",
1032
+                sep = ""
1033
+            ),
1034
+            names(object), genome(object),
1035
+            gsub("^chr", "", chromosome(object)),
1036
+            length(object),
1037
+            ifelse(length(object) == 1, "", "s")
1038
+        ))
1039
+    }
1040
+)
1041
+
1042
+#' @describeIn AlignmentsTrack-class Show method.
1043
+#' @export
1044
+setMethod("show", signature(object = "ReferenceAlignmentsTrack"), function(object) {
1045
+    .referenceTrackInfo(object, "ReferenceAlignmentsTrack")
1046
+})