#' @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") })