R/covPlotObject-methods.R
In alexchwong/SpliceWiz: interactive analysis and visualization of alternative splicing in R
Defines functions
.gcd_stack_anno
.gpo_filter_anno
.gpo_highlight_anno
.pV_mod_ggplot
.pV_assemble_ggplot
.pV_assemble_covPlotly2
.plotly_makeRectData
.plotly_makeLineData
.plotly_makeJuncCurveData
.pV_trackList_from_list
.pV_trackList_from_vector
.gPO_getJunc
.gPO_PSIjuncStats
.pV_diffStat_ttest
.gPO_normJuncStats
.gPO_covStats
.gPO_getCoverage
.gPO_check_strand
.gPO_check_tracks
.pV_sashimiArc
.pV_seed_plotly_HLexons
.pV_seed_plotly_annoTrack
.pV_seed_instant_plotly_HLexons
.pV_seed_instant_plotly_annoTrack
.pV_seed_ggplot_HLexons
.pV_seed_ggplot_annoTrack
.pV_highlight_to_colors
.pV_plotExons
.pV_plotRef
.pV_seed_plotly_diffTrack
.pV_seed_ggplot_diffTrack
.pV_seed_plotly_covTrack
.pV_seed_ggplot_covTrack
.pV_plotDiffTracks
.pV_plotCoverageTracks
.pV_getJuncCoords
.pV_filterTranscripts
.pV_getEventCoords
.pV_getAllowedCoords
plotView
plotAnnoTrack
getPlotObject
covPlotObject
Documented in
getPlotObject
plotAnnoTrack
plotView
#' Versatile coverage plots for SpliceWiz
#'
#' Here, we implement fast and versatile ggplot and plotly based
#' coverage and sashimi plots. Users can plot with unlimited number
#' of individual, individual-normalized, or group-normalized tracks.
#' Also implemented is user-defined group-comparison differential plots
#' (including t-test plots). Additionally, users can generate ggplots
#' subsetted by exon groups. See details below.
#'
#' @details
#' The typical pipeline for plotting versatile coverage plots is as follows:
#'
#' * A `covDataObject` is generated by calling `getCoverageData()` using an
#' input `NxtSE` object. This step retrieves coverage, junction counts and
#' normalization data for the relevant genomic region being queried. A new
#' `covDataObject` is necessary when querying a new genomic region.
#' * A `covPlotObject` is generated by calling `getPlotObject()` using an
#' input `covDataObject`. This step retrieves alternative splicing event
#' specific data, such as normalized coverages, or group combined coverages.
#' A new `covPlotObject` is required when changing `condition`, `Event`,
#' `strand`, or when querying using a different set of `tracks`.
#' * Plots can be generated by calling `plotView()` using a `covPlotObject`.
#' Interactive plotly plots can be generated by setting `usePlotly = TRUE`,
#' otherwise, static plots are generated. For interactive plots, a `covPlotly`
#' object is returned, which contains raw data which is downsampled by
#' pixel resolution prior to plotting for performance reasons. A new
#' `covPlotly` is required unless one only wishes to downsample the resolution
#' - see [setResolution] for `covPlotly` objects.
#'
#' Tracks are now versatile (unlimited). Samples are retrieved by individual
#' sample names at `getCoverageData()`. If `condition` is set in
#' `getPlotObject()`, track names are defined by their condition categorical
#' names; otherwise, tracks are named by individual samples when retrieved using
#' `getPlotObject()`.
#'
#' * When calling `plotView()`, `trackList` by default displays all tracks as
#' ordered in the `covPlotObject`. Users can supply a vector containing either
#' the track names (or numbers, as ordered in the `covPlotObject`).
#' Alternatively, multiple traces can be stacked in a single track by
#' using a list, e.g. `trackList = list(A = c(1,2,3), B = c(4,5,6))`.
#' * For differential comparisons, `diffList` takes a list of pairs of samples.
#' For example, if `trackList = list("A", "B")`, then setting
#' `diffList = list(c("A", "B"))` will compare groups "A" and "B". This is only
#' activated by setting `diff_stat` to anything other than `none`. For now,
#' only `t-test` is supported.
#'
#' `plotView()` supports plotting by exon ranges, for which only static plots
#' are currently supported. The workflow for generating such a plot is as
#' follows:
#'
#' * A GRanges object is returned by the `plotView()` function and setting
#' `showExonRanges = TRUE`. `plotView()` will simultaneously show an
#' annotation plot of exons labelled by their "exon names", which is the
#' transcript name appended with "-E" followed by the exon number.
#' * If `plotView()` is called and `usePlotly = TRUE` is set, a `covPlotly`
#' object is returned. Calling `showExons()` on this object will display
#' a plotly plot showing exon names, and returning a GRanges object of exon
#' ranges.
#' * Exon ranges can be supplied to the `plotView()` function by setting
#' the `plotRanges` parameter as a GRanges object. This will generate a static
#' plot showing coverage plots segmented by exons.
#'
#' @param object
#' For `getPlotObject()`, a `covDataObject` created using
#' `getCoverageData`. For `plotView()`, a `covPlotObject` created using
#' `getPlotObject()`.
#' @param Event The EventName of the alternative splicing event
#' which will be highlighted and used for normalization
#' @param strand The strand for coverage / junction plotting. Options are `"+"`,
#' `"-"`, or `"*"` (unstranded - default)
#' @param tracks Sample names or condition categories
#' @param condition For condition-based group plots, the name of the condition.
#' @param view_start,view_end The start and end coordinates for plotting
#' @param oldP (Optional) If plotting the same tracks and track-widths,
#' supplying the old `covPlotly` object (returned from a previous call to
#' `plotView()`) results in faster run-time (as plotly::subplot is a time-
#' consuming function)
#' @param centerByEvent (default `FALSE`) If true, centers the view to the
#' specified `Event`
#' @param EventZoomFactor If `centerByEvent = TRUE`, the zoom-out factor to plot
#' the view. Zooms out in exponents of 3
#' (i.e., zoom of 1 means 3x, 2 means 9x, and 0 means 1x)
#' @param EventBasesFlanking (default `100`)
#' If `centerByEvent = TRUE`, includes how many bases flanking the event.
#' @param resolution The number of horizontal "pixels" or data-points to plot.
#' This is calculated per sub-plot. Smaller numbers lead to lower resolution
#' but faster plots. Default is `5000`
#' @param trackList A list, with each element being a vector of 1 or more track
#' names or indices to plot. If a vector is supplied it will be coerced to a
#' list
#' @param diffList A list, with each element being a vector of size 2,
#' containing names or indices of which tracks to contrast.
#' @param diff_stat (default "t-test") Which statistical method to perform
#' differential comparisons.
#' @param reverseGenomeCoords If `TRUE`, the genomic coordinate axis will be
#' reversed to plot negative stranded genes
#' @param ribbon_mode The statistic to represent variance. Options are "sd" -
#' standard deviation, "sem" - standard error of the mean, "ci" -
#' 95% confidence interval, or "none"
#' @param normalizeCoverage If `TRUE`, coverages and junctions of individual
#' samples will be normalized by the given `Event`.
#' @param plotAnnotations Whether the main annotation track should be plotted
#' @param plotAnnoSubTrack If plotting by exon ranges (using `plotRanges`),
#' whether a separate sub-track showing zoomed-in exons should be shown
#' above the main annotation track (and below the coverage plots)
#' @param showExonRanges (only applies if `usePlotly = FALSE`) Whether the
#' main annotation track should be replaced by labeled exon names. If `TRUE`
#' the returned value of `plotView()` is a named GRanges object containing
#' the exon ranges
#' @param verticalLayout A vector (of length 4) containing relative heights of
#' the following elements: (1) main block of coverage tracks, (2) differential
#' track, (3) annotation sub-track, and (4) main annotation track. Default
#' `c(4,1,1,2)`
#' @param horizontalLayout A vector containing relative widths of coverage
#' tracks. Only used alongside `plotRanges` with more than 1 range to plot.
#' If omitted, `plotView` will attempt to scale widths to the widths of the
#' exon ranges.
#' @param filterByTranscripts (default `""`) One or more named transcripts to
#' filter the annotation track.
#' @param filterByEventTranscripts (default `FALSE`) If `TRUE`, only transcripts
#' involved in the given `Event` will be plotted, if any
#' @param filterByExpressedTranscripts (default `TRUE`) Only transcripts with
#' supported junctions will be plotted on the annotation axis. An expressed
#' junction is that which contains more than the minimum `junctionThreshold`
#' in at least 1 track
#' @param condenseTranscripts Whether to plot by genes `TRUE` or transcripts
#' `FALSE`
#' @param plotJunctions Whether to plot junction counts as numbered arcs. Plots
#' normalized junctions if `normalizeCoverage = TRUE`.
#' @param plotJuncPSI If plotting group coverage plots, whether to plot
#' mean +/- sd of normalized junction counts `FALSE`, or estimated junction
#' PSI based on SpliceOver metric applied to each junction `TRUE`.
#' @param junctionThreshold (default `0.01`) Junctions with expressions below
#' this threshold will not be plotted. For raw counts, this is a fraction of
#' maximum coverage value of the track.
#' @param plotRanges A GRanges object containing one or more exon ranges to
#' plot. If given, `view_start` and `view_end` will be ignored. Typical use
#' is to use the output of the `plotView(..., usePlotly = FALSE)`, which
#' returns a named GRanges object, then subset this output by exon name.
#' @param rangesBasesFlanking (default `100`) How many flanking bases to add to
#' each of `plotRanges`. Ignored if only 1 range given (or using
#' `view_start` and `view_end`)
#' @param usePlotly If `TRUE`, returns a `covPlotly` object containing the
#' plotly-based interactive plot. If `FALSE`, returns a ggplot object.
#' @param ... Ignored / not used
#' @return
#' For `getPlotObject()`: A `covPlotObject` object containing Event-based
#' data to create coverage plots using `plotView()`.
#'
#' For `plotView()`:
#'
#' * If `usePlotly = TRUE`, returns a `covPlotly` object containing plotly-based
#' interactive plot
#' * If `usePlotly = FALSE`, returns a patchwork-assembled static plot, unless
#' `showExonRanges = TRUE` in which it shows the plot and returns a named
#' GRanges object containing exon ranges.
#' @examples
#' se <- SpliceWiz_example_NxtSE(novelSplicing = TRUE)
#'
#' # Assign annotation of the experimental conditions
#' colData(se)$treatment <- rep(c("A", "B"), each = 3)
#'
#' # Retrieve coverage data for all samples for the gene "SRSF3" (and surrounds)
#'
#' dataObj <- getCoverageData(
#' se,
#' Gene = "SRSF3",
#' tracks = colnames(se)
#' )
#'
#' # Retrieves raw / normalized coverage / junction data for the
#' # specified SRSF3 skipped exon event:
#'
#' plotObj_samples <- getPlotObject(
#' dataObj,
#' Event = "SE:SRSF3-203-exon4;SRSF3-202-int3"
#' )
#'
#' # Retrieves data for samples grouped by the specified condition
#'
#' plotObj_group <- getPlotObject(
#' dataObj,
#' Event = "SE:SRSF3-203-exon4;SRSF3-202-int3",
#' condition = "treatment",
#' tracks = c("A", "B")
#' )
#'
#' # Display tracks and conditions of covPlotObject
#'
#' tracks(plotObj_group)
#' condition(plotObj_group)
#'
#' # Show static ggplots
#'
#' plotView(plotObj_samples)
#'
#' plotView(plotObj_group, centerByEvent = TRUE)
#'
#' # Plot junctions using PSI estimates based on individual junction SpliceOver
#' # metrics
#'
#' plotView(plotObj_group, centerByEvent = TRUE, plotJuncPSI = TRUE)
#'
#'
#' # Show normalized coverages, individual samples stacked in grouped tracks
#'
#' plotView(
#' plotObj_samples,
#' normalizeCoverage = TRUE,
#' trackList = list(A = c(1,2,3), B = c(4,5,6))
#' )
#'
#' # Show stacked group comparisons with t-test
#'
#' plotView(
#' plotObj_group,
#' trackList = list(c(1,2)),
#' diffList = list(c("A", "B")),
#' diff_stat = "t-test"
#' )
#'
#' # Show interactive plotly:
#'
#' if(interactive()) {
#' p <- plotView(plotObj_samples, usePlotly = TRUE)
#' show(p)
#' }
#'
#' # Show exons with coverage plot
#'
#' # static:
#' gr <- plotView(plotObj_samples, showExonRanges = TRUE)
#'
#' # interactive:
#' if(interactive()) {
#' p <- plotView(plotObj_samples, usePlotly = TRUE)
#' gr <- showExons(p)
#' }
#'
#' # Plot coverage by exons
#'
#' p <- plotView(plotObj_samples,
#' plotRanges = gr[c("SRSF3-203-E3", "SRSF3-203-E4", "SRSF3-203-E5")],
#' horizontalLayout = c(1,1,1)
#' )
#'
#' @name covPlotObject-class
#' @aliases
#' tracks tracks,covPlotObject-method
#' condition condition,covPlotObject-method
#' @seealso [getCoverageData] [covPlotly-class]
#' @md
NULL
covPlotObject <- function(
args = list(),
cov = list(),
norm_cov = list(),
junc = list(),
norm_junc = list(),
junc_PSI = list(),
junc_stats = list(),
cov_stats = list(),
annotation = list()
) {
obj <- new("covPlotObject",
args = args,
cov = cov,
norm_cov = norm_cov,
junc = junc,
norm_junc = norm_junc,
junc_PSI = junc_PSI,
junc_stats = junc_stats,
cov_stats = cov_stats,
annotation = annotation
)
obj
}
#' @describeIn covPlotObject-class Generates a covPlotObject object from a
#' covDataObject. Allows users to change parameters such as viewing window,
#' conditions, tracks, and other parameters, for customizing plot parameters
#' @export
getPlotObject <- function(
object,
Event, # To specify normalization event
strand = c("*", "+", "-"),
tracks, # can be a list of iterables
condition,
...
) {
args <- object@args
if(!missing(tracks)) args[["tracks"]] <- tracks
if(!missing(condition)) args[["condition"]] <- condition
if(!missing(Event)) args[["Event"]] <- Event
if(!("Event" %in% names(args)) && "condition" %in% names(args)) {
.log(paste("In getPlotObject,",
"for plotting by condition, normalizing `Event` must be supplied"
))
}
# Check tracks and conditions in args are legit, modify as necessary
args <- .gPO_check_tracks(object, args)
# Junction strand is always same as requested strand
# But need to check strandedness of RNA-seq data before deciding coverage
if(
!missing(strand) &&
length(strand) == 1 &&
strand %in% c("*", "+", "-")
) {
args[["strand"]] <- strand
} else if(!("strand" %in% names(args))) {
args[["strand"]] <- "*" # default strand
}
# Check strand against stranded-ness of RNA-seq
args <- .gPO_check_strand(object, args)
# Always get these
cov <- .gPO_getCoverage(object, args)
junc <- .gPO_getJunc(object, args, cov)
norm_cov <- norm_junc <- list()
if("Event" %in% names(args)) {
norm_cov <- .gPO_getCoverage(object, args, normalize = TRUE)
# if single sample
norm_junc <- .gPO_getJunc(object, args, norm_cov, normalized = TRUE)
}
# Plot stats if applicable
cov_stats <- junc_PSI <- junc_stats <- list()
if("condition" %in% names(args)) {
# compile junction data
cov_stats <- .gPO_covStats(norm_cov)
junc_stats <- .gPO_normJuncStats(object, args, cov_stats)
junc_PSI <- .gPO_PSIjuncStats(object, args, norm_cov)
}
# plot relevant annotations
highlight_gr <- list()
if(
"Event" %in% names(args) &&
args[["Event"]] %in% rownames(object@normData$rowData)
) {
row <- object@normData$rowData[args[["Event"]],]
args[["EventRegion"]] <- row$EventRegion
if (row$EventType %in% c("MXE", "SE")) {
highlight_gr[[1]] <- coord2GR(c(row$Event1a, row$Event2a))
} else {
highlight_gr[[1]] <- coord2GR(row$Event1a)
}
if (row$EventType %in% c("MXE")) {
highlight_gr[[2]] <- coord2GR(c(row$Event1b, row$Event2b))
} else if (row$EventType %in% c("SE", "A3SS", "A5SS", "ALE", "AFE")) {
highlight_gr[[2]] <- coord2GR(row$Event1b)
}
}
args[["highlight_gr"]] <- highlight_gr
DTlist <- object@annotation
DTlist$reduced.DT <- .gpo_highlight_anno(
DTlist$reduced.DT, args[["highlight_gr"]]
)
return(covPlotObject(
args = args,
cov = cov,
norm_cov = norm_cov,
junc = junc,
norm_junc = norm_junc,
junc_PSI = junc_PSI,
junc_stats = junc_stats,
cov_stats = cov_stats,
annotation = DTlist
))
}
#' @describeIn covPlotObject-class Returns the tracks contained in the
#' covPlotObject object
#' @export
setMethod("tracks", c(object = "covPlotObject"), function(
object
) {
return(object@args[["tracks"]])
})
#' @describeIn covPlotObject-class Returns the condition value set in the
#' covPlotObject object
#' @export
setMethod("condition", c(object = "covPlotObject"), function(
object
) {
return(object@args[["condition"]])
})
#' @describeIn covDataObject-class Directly plots the annotation from a
#' covDataObject.
#' @export
plotAnnoTrack <- function(
object, Event,
view_start, view_end,
reverseGenomeCoords = FALSE,
condensed = FALSE,
selected_transcripts = "",
plot_key_isoforms = FALSE,
usePlotly = FALSE,
...
) {
if(missing(view_start)) view_start <- object@args$view_start
if(missing(view_end)) view_end <- object@args$view_end
if(view_start > view_end) {
reverseGenomeCoords <- !reverseGenomeCoords
view_start <- view_start + view_end
view_end <- view_start - view_end
view_start <- view_start - view_end
}
highlight_gr <- list()
if(!missing(Event) && !("rowData" %in% names(object@normData))) {
.log(paste(
"In plotAnnoTrack(),",
"Plotting ASE from a SpliceWiz reference is currently not supported"
))
}
if(!missing(Event) && Event %in% rownames(object@normData$rowData)) {
row <- object@normData$rowData[Event,]
if (row$EventType %in% c("MXE", "SE")) {
highlight_gr[[1]] <- coord2GR(c(row$Event1a, row$Event2a))
} else {
highlight_gr[[1]] <- coord2GR(row$Event1a)
}
if (row$EventType %in% c("MXE")) {
highlight_gr[[2]] <- coord2GR(c(row$Event1b, row$Event2b))
} else if (row$EventType %in% c("SE", "A3SS", "A5SS", "ALE", "AFE")) {
highlight_gr[[2]] <- coord2GR(row$Event1b)
}
}
DTlist <- object@annotation
DTlist$reduced.DT <- .gpo_highlight_anno(
DTlist$reduced.DT, highlight_gr
)
DTlist$reduced.DT <- .gpo_filter_anno(
DTlist$reduced.DT, DTlist$transcripts.DT,
selected_transcripts = selected_transcripts,
plot_key_isoforms = plot_key_isoforms
)
DTplotlist <- .gcd_stack_anno(
DTlist, view_start, view_end, reverseGenomeCoords
)
view_gr <- GRanges(
object@args[["view_chr"]],
IRanges(object@args[["view_start"]], object@args[["view_end"]])
)
p <- .pV_plotRef(
DTplotlist, view_gr,
reverseGenomeCoords = reverseGenomeCoords,
add_information = TRUE,
usePlotly = usePlotly,
instantPlotly = usePlotly
)
return(p)
}
#' @describeIn covPlotObject-class Creates a coverage plot using the stored
#' data in the covPlotObject
#' @export
plotView <- function(
object,
# for single window view (to simplify things)
view_start,
view_end,
# debug = FALSE,
oldP = covPlotly(),
# Event-centric plotting
centerByEvent = FALSE,
EventZoomFactor = 0.2,
EventBasesFlanking = 100,
resolution = 5000,
# specify tracks
trackList = list(),
# specify differential comparisons
diff_stat = c("t-test", "none"),
diffList = list(),
reverseGenomeCoords = FALSE,
ribbon_mode = c("sd", "sem", "ci", "none"),
# plot raw or normalized coverage (only works for sample tracks)
normalizeCoverage = FALSE,
# specify annotation track
plotAnnotations = TRUE,
plotAnnoSubTrack = TRUE,
showExonRanges = FALSE,
# Plot layout
verticalLayout = c(4,1,1,2),
horizontalLayout = c(),
# filter annotations
filterByTranscripts = "",
# plot only by transcripts with highlights
filterByEventTranscripts = FALSE,
# filter annotations
filterByExpressedTranscripts = TRUE,
# whether to condense transcripts by gene
condenseTranscripts = FALSE,
# whether to plot sashimi arcs
plotJunctions = TRUE,
plotJuncPSI = FALSE,
junctionThreshold = 0.01,
# whether to plot by exon windows
plotRanges = GRanges(),
rangesBasesFlanking = 100,
# ggplot or plotly object?
usePlotly = FALSE,
...
) {
use_DT <- TRUE
if(!is(object, "covPlotObject")) .log(paste(
"In plotView,", "object must be a covPlotObject"
))
args <- object@args
# inject variables into args
if(!missing(view_start)) args[["view_start"]] <- view_start
if(!missing(view_end)) args[["view_end"]] <- view_end
args[["reverseGenomeCoords"]] <- reverseGenomeCoords
ribbon_mode <- match.arg(ribbon_mode)
if(!is_valid(ribbon_mode)) ribbon_mode <- "none"
args[["ribbon_mode"]] <- ribbon_mode
if(
(missing(view_start) | missing(view_end)) &&
centerByEvent &
"EventRegion" %in% names(args)
) {
args <- .pV_getEventCoords(
object, args, EventZoomFactor,
EventBasesFlanking
)
}
diff_stat <- match.arg(diff_stat)
if(!is_valid(diff_stat)) diff_stat <- "none"
# plotRanges checking
if(!is(plotRanges, "GRanges")) .log(paste(
"In plot() for class covPlotObject,",
"`plotRanges` must be a GRanges object"
))
# Plot single window if no plotRanges given
if(length(plotRanges) == 0) plotRanges <- GRanges(
args[["view_chr"]],IRanges(args[["view_start"]], args[["view_end"]])
) # strand is ignored
if(
any(as.character(seqnames(plotRanges)) != args[["view_chr"]])
) {
seqs <- unique(as.character(seqnames(plotRanges)))
.log(paste(
"In plotView() for class covPlotObject,",
"Some elements in `plotRanges` have seqnames that do not match that",
"of covPlotObject", seqs[!(seqs %in% args[["view_chr"]])]
))
}
# Make all plotRanges strand-agnostic
strand(plotRanges) <- "*"
# disable interactive if multi plot
if(usePlotly == TRUE & length(plotRanges) > 1) {
.log(paste(
"In plotView,",
"interactive plots are not supported for multi-view plotting"
), "warning")
usePlotly <- FALSE
}
# Expand plotRanges using rangesBasesFlanking
# - only used if more than 1 range
rbf <- rangesBasesFlanking
if(length(plotRanges) > 1) {
for(i in seq_len(length(plotRanges))) {
if(start(plotRanges[i]) - rbf > args[["limit_start"]]) {
start(plotRanges[i]) <- start(plotRanges[i]) - rbf
} else {
start(plotRanges[i]) <- args[["limit_start"]]
}
if(end(plotRanges[i]) + rbf < args[["limit_end"]]) {
end(plotRanges[i]) <- end(plotRanges[i]) + rbf
} else {
end(plotRanges[i]) <- args[["limit_end"]]
}
}
}
# What is the full range
fullRange <- range(plotRanges)
if(!missing(view_start) & !missing(view_end)) {
# whatever's larger
fullRange <- range(c(fullRange,
GRanges(args[["view_chr"]], IRanges(view_start, view_end))
))
}
if(
start(fullRange) < args[["limit_start"]] |
end(fullRange) > args[["limit_end"]]
) {
.log(paste(
"Given range is outside that supported by parent covPlotObject.",
"Suggest regenerating covPlotObject by calling getPlotObject()"
))
}
fetchRange <- fullRange
if(usePlotly) {
start(fetchRange) <- max(args[["limit_start"]],
start(fullRange) - width(fullRange))
end(fetchRange) <- min(args[["limit_end"]],
end(fullRange) + width(fullRange))
}
# Sort plotRanges
plotRanges <- sort(plotRanges, decreasing = args[["reverseGenomeCoords"]])
# Track list checking - should be a list of indices
# If no trackList given, plot all tracks individually
if(length(trackList) == 0) trackList <- lapply(
seq_len(length(args[["tracks"]])),
function(i) i
)
if(!is(trackList, "list")) {
args[["trackList"]] <- .pV_trackList_from_vector(trackList, args)
} else {
args[["trackList"]] <- .pV_trackList_from_list(trackList, args)
}
# Structure of plot
# | 1a || 1b || 1c |
# | 2a || 2b || 2c |
# | 3a || 3b || 3c |
# | pa || pb || pc | - diff track
# | A || B || C | - annotation subtrack
# | annotation | - annotation full track
covTrack <- list()
diffTrack <- list()
annoSubTrack <- list()
annoFullTrack <- list() # list of 1
exonFullTrack <- list() # list of 1
# Proportion all viewing frames
widthSum <- sum(width(plotRanges))
widthFrac <- width(plotRanges) / widthSum
if(use_DT && !usePlotly) use_DT <- FALSE
# For filtering transcripts by expression later
aggJuncCoords <- c()
reservedCoords <- c()
for(j in seq_len(length(plotRanges))) {
range_gr <- plotRanges[j]
# Work out which coords to plot
juncCoords <- .pV_getJuncCoords(
object, args, normalizeCoverage,
range_gr = range_gr,
junctionThreshold
)
aggJuncCoords <- c(aggJuncCoords, juncCoords)
}
reservedCoords <- sort(unique(c(
aggJuncCoords[["jc_start"]] - 1,
aggJuncCoords[["jc_start"]],
aggJuncCoords[["jc_end"]],
aggJuncCoords[["jc_end"]] + 1
)))
if(usePlotly) {
covTrack[[j]] <- .pV_plotCoverageTracks(
object, args, normalizeCoverage,
range_gr = fullRange,
plotJunctions = plotJunctions, plotJuncPSI = plotJuncPSI,
junctionThreshold = junctionThreshold,
usePlotly = TRUE,
use_DT = use_DT
)
if(diff_stat != "none" && length(diffList) > 0) {
diffTrack[[j]] <- .pV_plotDiffTracks(
object, args, diffList,
diff_stat, fullRange,
usePlotly = usePlotly
)
}
} else {
# Generate full plotly first, then segment later
tmpCov <- .pV_plotCoverageTracks(
object, args, normalizeCoverage,
range_gr = fullRange,
plotJunctions = plotJunctions, plotJuncPSI = plotJuncPSI,
junctionThreshold = junctionThreshold,
usePlotly = FALSE,
use_DT = use_DT
)
tmpCovPlot <- tmpCov$ggplot
tmpCovJuncAnno <- as.data.table(tmpCov$dfJn_all)
tmpDiffPlot <- NULL
if(diff_stat != "none" && length(diffList) > 0) {
tmpDiffPlot <- .pV_plotDiffTracks(
object, args, diffList,
diff_stat, fullRange,
usePlotly = FALSE
)
}
for(j in seq_len(length(plotRanges))) {
range_gr <- plotRanges[j]
okCoords <- .pV_getAllowedCoords(
start(range_gr), end(range_gr),
reservedCoords, resolution
)
if(args[["reverseGenomeCoords"]]) {
plotViewStart <- end(range_gr)
plotViewEnd <- start(range_gr)
} else {
plotViewStart <- start(range_gr)
plotViewEnd <- end(range_gr)
}
covTrack[[j]] <- list()
dfJnSum_list <- list()
if(plotJunctions) {
dtJn_all <- as.data.table(tmpCovJuncAnno)
allGroups <- unique(dtJn_all$group)
for(elem in allGroups) {
dtJn <- copy(dtJn_all[get("group") == elem])
dtJn <- dtJn[
get("x") >= start(range_gr) &
get("x") <= end(range_gr)
]
dtJn[, c("xlabel", "ylabel") := list(
mean(get("x")), mean(get("ytmp"))),
by = "info"
]
dtJn <- unique(
dtJn[,
c("info", "value", "xlabel",
"ylabel", "group", "covTrack"),
with = FALSE
]
)
dfJnSum_list[[elem]] <- dtJn
}
}
dfJnSum_all <- rbindlist(dfJnSum_list)
for(i in seq_len(length(tmpCovPlot))) {
curTrack <- names(args[["trackList"]])[i]
covTrack[[j]][[i]] <- tmpCovPlot[[i]]
ymax <- max(layer_scales(tmpCovPlot[[i]])$y$range$range)
if(nrow(dfJnSum_all) > 0) {
tmpCovSubset <- dfJnSum_all[get("covTrack") == curTrack]
# disable junction plots for multiple samples per track
if(
nrow(tmpCovSubset) > 0 &&
length(unique(tmpCovSubset$group)) == 1
) {
covTrack[[j]][[i]] <- covTrack[[j]][[i]] +
geom_text(
data = tmpCovSubset,
aes(
x = get("xlabel"),
y = get("ylabel"),
label = get("value")
)
)
}
}
if(!is.null(okCoords)) {
covTrack[[j]][[i]]$data <- covTrack[[j]][[i]]$data[
covTrack[[j]][[i]]$data$x %in% okCoords,]
}
covTrack[[j]][[i]] <- covTrack[[j]][[i]] +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(0, 1.2 * ymax),
expand = FALSE
)
}
if(!is.null(tmpDiffPlot)) {
diffTrack[[j]] <- list()
for(i in seq_len(length(tmpDiffPlot))) {
diffTrack[[j]][[i]] <- tmpDiffPlot[[i]]
ymax <- max(layer_scales(tmpDiffPlot[[i]])$y$range$range)
if(!is.null(okCoords)) {
diffTrack[[j]][[i]]$data <- diffTrack[[j]][[i]]$data[
diffTrack[[j]][[i]]$data$x %in% okCoords,]
}
diffTrack[[j]][[i]] <- diffTrack[[j]][[i]] +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(0, 1.2 * ymax),
expand = FALSE
)
}
}
}
}
# bench <- system.time({
DTlist <- .pV_filterTranscripts(
object, args,
filterByTranscripts = filterByTranscripts,
filterByEventTranscripts = filterByEventTranscripts,
filterByExpressedTranscripts = filterByExpressedTranscripts,
aggJuncCoords = aggJuncCoords
)
DTplotlist <- .gcd_stack_anno(
DTlist,
start(fetchRange), end(fetchRange),
reverseGenomeCoords = args[["reverseGenomeCoords"]],
condensed = condenseTranscripts
)
DTplotlist[["exonRanges"]] <- plotRanges
exons_gr <- GRanges()
if(length(plotRanges) > 1) {
if(plotAnnoSubTrack) {
for(j in seq_len(length(plotRanges))) {
range_gr <- plotRanges[j]
annoSubTrack[[j]] <- .pV_plotRef(
DTplotlist, range_gr,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = FALSE,
usePlotly = usePlotly
)
}
if(plotAnnotations) {
annoFullTrack[[1]] <- .pV_plotRef(
DTplotlist, fullRange,
showPlotRanges = plotRanges,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
}
} else if(plotAnnotations) {
annoFullTrack[[1]] <- .pV_plotRef(
DTplotlist, fullRange,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
}
} else {
# Annotate and highlight exon regions
transcripts <- copy(DTlist$transcripts.DT)
reduced <- copy(DTplotlist$reduced.DT)
reduced <- reduced[!is.na(reduced$plot_level)]
exonsOnly <- reduced[get("type") == "exon"]
exonsOnly$transcript_name <- transcripts$transcript_name[match(
exonsOnly$transcript_id, transcripts$transcript_id)]
exonsOnly[, c("exon_name") := paste0(
get("transcript_name"), "-E", get("aux_id")
)]
exonsOnly[, c("feature_id") := NULL]
exons_gr <- .grDT(exonsOnly, keep.extra.columns = TRUE)
names(exons_gr) <- exons_gr$exon_name
# Remove exons that lie outside view range
OL <- findOverlaps(exons_gr, fetchRange, type = "within")
exons_gr <- exons_gr[sort(unique(OL@from))]
if(!plotAnnotations & showExonRanges) showExonRanges <- FALSE
if(plotAnnotations) {
exonFullTrack[[1]] <- .pV_plotExons(
DTplotlist, fullRange, exons_gr,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
annoFullTrack[[1]] <- .pV_plotRef(
DTplotlist, fullRange,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
}
}
# })
# message("Annotation data generation time:")
# print(bench)
keepVLayout <- c(
length(covTrack) > 0,
length(diffTrack) > 0,
length(annoSubTrack) > 0,
length(annoFullTrack) > 0
)
originalVerticalLayout <- verticalLayout
if(length(originalVerticalLayout) == 4) {
originalVerticalLayout <- originalVerticalLayout[-3]
}
if(length(verticalLayout) == 4) {
verticalLayout <- verticalLayout[keepVLayout]
} else if(length(verticalLayout) == sum(keepVLayout)) {
# do nothing
} else {
# revert to default
verticalLayout <- c(4,1,1,2)
verticalLayout <- verticalLayout[keepVLayout]
}
if(length(horizontalLayout) == 0) {
horizontalLayout <- ceiling(10 * widthFrac / sum(widthFrac))
} else if(length(horizontalLayout) == length(covTrack)) {
horizontalLayout <- ceiling(horizontalLayout)
} else {
horizontalLayout <- ceiling(10 * widthFrac / sum(widthFrac))
}
if(!args[["reverseGenomeCoords"]]) {
plotViewStart <- start(fullRange)
plotViewEnd <- end(fullRange)
} else {
plotViewStart <- end(fullRange)
plotViewEnd <- start(fullRange)
}
# if(debug) {
# return(list(
# plotViewStart = plotViewStart, plotViewEnd = plotViewEnd,
# covTrack = covTrack,
# diffTrack = diffTrack,
# annoSubTrack = annoSubTrack,
# annoFullTrack = annoFullTrack,
# exonFullTrack = exonFullTrack
# ))
# }
# bench <- system.time({
if(usePlotly) {
if(is(oldP, "covPlotly")) {
p <- oldP
p@covTrack <- covTrack
p@diffTrack <- diffTrack
p@annoTrack <- annoFullTrack
p@exonTrack <- exonFullTrack
} else {
p <- covPlotly(
covTrack = covTrack,
diffTrack = diffTrack,
annoTrack = annoFullTrack
)
}
p@args[["resolution"]] <- resolution
p@args[["reservedCoords"]] <- reservedCoords
p@args[["exonRanges"]] <- exons_gr
p <- .pV_assemble_covPlotly2(p, originalVerticalLayout)
p@fig[[2]] <- p@fig[[2]] %>% layout(
dragmode = "pan",
xaxis = list(
range = c(plotViewStart, plotViewEnd)
)
)
return(p)
} else {
if(showExonRanges) {
p <- .pV_assemble_ggplot(
covTrack, diffTrack, annoSubTrack, exonFullTrack,
verticalLayout, horizontalLayout
)
show(p)
return(exons_gr)
} else {
p <- .pV_assemble_ggplot(
covTrack, diffTrack, annoSubTrack, annoFullTrack,
verticalLayout, horizontalLayout
)
return(p)
}
}
# })
# message("Plotly generation time:")
# print(bench)
return(p)
}
################################################################################
# Given a range, resolution, and must-include coords
# - returns a set of coordinates to keep in final plot
.pV_getAllowedCoords <- function(
viewStart, viewEnd, reservedCoords, resolution
) {
rC <- reservedCoords[reservedCoords %in% seq(viewStart, viewEnd)]
okCoords <- NULL
if(resolution > 0 & !is.null(rC)) {
okCoords <- round(seq(
from = min(viewStart, viewEnd),
to = max(viewStart, viewEnd),
length.out = max(2, resolution - length(rC))
))
okCoords <- sort(unique(c(okCoords, rC)))
}
return(okCoords)
}
# Returns view_start / view_end of coords centered around event
.pV_getEventCoords <- function(
x, args,
zoom_factor = 0.2,
bases_flanking = 100
) {
# sanitise parameters
if(!is.numeric(bases_flanking) || bases_flanking < 0) bases_flanking <- 0
if(!is.numeric(zoom_factor) || zoom_factor < 0) zoom_factor <- 0
if(!("EventRegion" %in% names(args))) return(NULL)
EventRegion <- args[["EventRegion"]]
grEvent <- coord2GR(EventRegion)
eventStart <- start(grEvent)
eventEnd <- end(grEvent)
view_length <- eventEnd - eventStart
view_center <- (eventEnd + eventStart) / 2
new_view_length <- view_length * 3^zoom_factor + 2 * bases_flanking
view_start <- round(view_center - new_view_length / 2)
view_end <- round(view_center + new_view_length / 2)
seqInfo <- args[["seqInfo"]][args[["view_chr"]]]
seqmax <- GenomeInfoDb::seqlengths(seqInfo)
if (view_end > seqmax) {
view_end <- seqmax - 1
view_start <- view_end - new_view_length
if(view_start < 1) view_start <- 1
}
args[["view_start"]] <- view_start
args[["view_end"]] <- view_end
return(args)
}
################################################################################
# Returns reduced.DT and transcript.DT after filtering annotations by
# given criteria
.pV_filterTranscripts <- function(
x, args,
filterByTranscripts,
filterByEventTranscripts,
filterByExpressedTranscripts,
aggJuncCoords
) {
# Filter annotations by observed coords
reduced.DT <- copy(x@annotation[["reduced.DT"]])
transcripts.DT <- copy(x@annotation[["transcripts.DT"]])
transcriptFiltered <- FALSE
if(is_valid(filterByTranscripts)) {
if(
all(filterByTranscripts %in%
c(
transcripts.DT$transcript_id,
transcripts.DT$transcript_name
)
)
) {
transcripts.DT <- transcripts.DT[
get("transcript_name") %in% filterByTranscripts |
get("transcript_id") %in% filterByTranscripts
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
transcriptFiltered <- TRUE
}
}
if(!transcriptFiltered & filterByEventTranscripts) {
reduced.DT.HL <- reduced.DT[get("highlight") != "0"]
reduced.DT.nonHL <- reduced.DT[get("highlight") == "0"]
reduced.DT.nonHL <- reduced.DT.nonHL[
grepl("novel", get("transcript_id"))
]
if(nrow(reduced.DT.HL) != 0) {
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% reduced.DT.HL$transcript_id &
# remove non-highlighted novel elements
!(get("transcript_id") %in% reduced.DT.nonHL$transcript_id)
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
transcriptFiltered <- TRUE
}
}
if(!transcriptFiltered && filterByExpressedTranscripts) {
reduced.DT.introns <- reduced.DT[get("type") == "intron"]
expressed.introns <- reduced.DT.introns[
get("start") %in% aggJuncCoords[["jc_start"]] &
get("end") %in% aggJuncCoords[["jc_end"]]
]
nonexpressed.novel <- reduced.DT.introns[
!(get("start") %in% aggJuncCoords[["jc_start"]]) |
!(get("end") %in% aggJuncCoords[["jc_end"]])
]
nonexpressed.novel <- nonexpressed.novel[
grepl("novel", get("transcript_id"))
]
if(nrow(expressed.introns) > 0) {
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% expressed.introns$transcript_id &
# remove non-highlighted novel elements
!(get("transcript_id") %in% nonexpressed.novel$transcript_id)
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
transcriptFiltered <- TRUE
}
}
return(list(
reduced.DT = reduced.DT,
transcripts.DT = transcripts.DT
))
}
################################################################################
# Given a list of tracks (in args)
# - returns a list of junction start/end coordinates of expressed transcripts
.pV_getJuncCoords <- function(
x, args,
normalizeCoverage = FALSE,
range_gr,
junctionThreshold = 0.01
) {
plotMeanCov <- FALSE
if("condition" %in% names(args)) plotMeanCov <- TRUE
if(normalizeCoverage) {
if(length(x@norm_cov) == 0 || plotMeanCov) {
normalizeCoverage <- FALSE
}
}
trackList <- args[["trackList"]]
fetch_gr <- range_gr
start(fetch_gr) <- max(args[["limit_start"]],
start(range_gr) - width(range_gr))
end(fetch_gr) <- min(args[["limit_end"]],
end(range_gr) + width(range_gr))
jc_start <- c()
jc_end <- c()
for(track_id in seq_len(length(trackList))) {
idx <- trackList[[track_id]]
for(elem in args[["tracks"]][idx]) {
if(plotMeanCov) {
junc <- x@junc_PSI[[elem]]
} else if(normalizeCoverage) {
junc <- x@norm_junc[[elem]]
} else {
junc <- x@junc[[elem]]
}
OL <- findOverlaps(junc, fetch_gr)
if(length(from(OL)) > 0) {
junc <- junc[unique(from(OL))]
} else {
junc <- NULL
}
if(!is.null(junc)) {
isMean <- ("mean" %in% names(mcols(junc)))
if(isMean) {
junc <- junc[mcols(junc)$mean >= junctionThreshold]
} else if(normalizeCoverage) {
junc <- junc[mcols(junc)$count >= junctionThreshold]
} else {
cov <- x@cov[[elem]]
junc <- junc[mcols(junc)$count >=
junctionThreshold * max(cov$depth)
]
}
if(length(junc) > 0) {
jc_start <- c(jc_start, start(junc))
jc_end <- c(jc_end, end(junc))
}
}
}
}
return(list(
jc_start = jc_start, jc_end = jc_end
))
}
# Creates coverage plot tracks
.pV_plotCoverageTracks <- function(
x, args,
normalizeCoverage,
range_gr,
plotJunctions = TRUE, plotJuncPSI = FALSE,
junctionThreshold = 0.01,
usePlotly = FALSE,
use_DT = FALSE
) {
plotMeanCov <- FALSE
if("condition" %in% names(args)) plotMeanCov <- TRUE
rb <- "none"
if(plotMeanCov) rb <- args[["ribbon_mode"]]
if(normalizeCoverage) {
if(length(x@norm_cov) == 0 || plotMeanCov) {
normalizeCoverage <- FALSE
}
}
# junction processing is identical
trackName <- "Coverage"
trackList <- args[["trackList"]]
# plot view x cartesian coordinates
if(args[["reverseGenomeCoords"]]) {
plotViewStart <- end(range_gr)
plotViewEnd <- start(range_gr)
} else {
plotViewStart <- start(range_gr)
plotViewEnd <- end(range_gr)
}
# fetchViewStart and fetchViewEnd determine how much coverage to fetch
# best to fetch a bit more than we need
fetch_gr <- range_gr
start(fetch_gr) <- max(args[["limit_start"]],
start(range_gr) - width(range_gr))
end(fetch_gr) <- min(args[["limit_end"]],
end(range_gr) + width(range_gr))
DT_list <- list()
DTJn_list <- list()
for(track_id in seq_len(length(trackList))) {
idx <- trackList[[track_id]]
for(elem in args[["tracks"]][idx]) {
if(plotJunctions) {
if(plotMeanCov) {
if(plotJuncPSI) {
junc <- x@junc_PSI[[elem]]
} else {
junc <- x@junc_stats[[elem]]
}
} else if(normalizeCoverage) {
junc <- x@norm_junc[[elem]]
} else {
junc <- x@junc[[elem]]
}
OL <- findOverlaps(junc, fetch_gr)
if(length(from(OL)) > 0) {
junc <- junc[unique(from(OL))]
} else {
junc <- NULL
}
} else {
junc <- NULL
}
if(plotMeanCov) {
groupName <- paste(x@args[["condition"]], elem)
DT <- x@cov_stats[[elem]]
DT_sub <- DT[,seq_len(2)] # only 2 columns
DTJn <- c()
if(!is.null(junc)) {
DTJn <- .pV_sashimiArc(
junc,
arcHeight = 0.1 * max(DT$mean),
junctionThreshold = junctionThreshold * max(DT$mean)
)
}
} else {
groupName <- elem
if(normalizeCoverage) {
DT <- x@norm_cov[[elem]]
} else {
DT <- x@cov[[elem]]
}
DTJn <- c()
if(!is.null(junc)) {
DTJn <- .pV_sashimiArc(
junc,
arcHeight =
0.1 * max(DT$depth),
junctionThreshold = junctionThreshold *
max(DT$depth)
)
}
}
DT$group <- groupName
DT$covTrack <- names(trackList)[track_id]
DT_list[[groupName]] <- DT
if(!is.null(junc) && nrow(DTJn) > 0) {
DTJn$group <- groupName
DTJn$covTrack <- names(trackList)[track_id]
DTJn_list[[groupName]] <- DTJn
}
}
}
df_all <- rbindlist(DT_list)
dfJn_all <- rbindlist(DTJn_list)
if(!use_DT) {
df_all <- as.data.frame(df_all)
dfJn_all <- as.data.frame(dfJn_all)
}
# group annotation
df_all$info <- ""
# junc processing - common
dfJnSum_list <- list()
if(plotJunctions) {
dtJn_all <- as.data.table(dfJn_all)
allGroups <- unique(dtJn_all$group)
for(elem in allGroups) {
dtJn <- copy(dtJn_all[get("group") == elem])
dtJn <- dtJn[
get("x") >= start(range_gr) &
get("x") <= end(range_gr)
]
dtJn[, c("xlabel", "ylabel") := list(
mean(get("x")), mean(get("ytmp"))),
by = "info"
]
dtJn <- unique(
dtJn[,
c("info", "value", "xlabel", "ylabel", "group", "covTrack"),
with = FALSE
]
)
dfJnSum_list[[elem]] <- dtJn
}
}
if(use_DT) {
dfJnSum_all <- rbindlist(dfJnSum_list)
} else {
dfJnSum_all <- as.data.frame(rbindlist(dfJnSum_list))
}
# factor covTrack
df_all$covTrack <- factor(df_all$covTrack, names(trackList))
if(nrow(dfJn_all) > 0) {
dfJn_all$covTrack <- factor(dfJn_all$covTrack, names(trackList))
}
if(nrow(dfJnSum_all) > 0) {
dfJnSum_all$covTrack <- factor(dfJnSum_all$covTrack, names(trackList))
}
if(!usePlotly) {
p <- .pV_seed_ggplot_covTrack(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
)
return(list(
ggplot = p,
dfJn_all = dfJn_all
))
} else {
p <- .pV_seed_plotly_covTrack(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
)
return(p)
}
}
# Creates differential plot tracks
.pV_plotDiffTracks <- function(
x, args, diffList,
diff_stat,
range_gr,
usePlotly = FALSE
) {
# plot view x cartesian coordinates
if(args[["reverseGenomeCoords"]]) {
plotViewStart <- end(range_gr)
plotViewEnd <- start(range_gr)
} else {
plotViewStart <- start(range_gr)
plotViewEnd <- end(range_gr)
}
# fetchViewStart and fetchViewEnd determine how much coverage to fetch
# best to fetch a bit more than we need
fetch_gr <- range_gr
start(fetch_gr) <- max(args[["limit_start"]],
start(range_gr) - width(range_gr))
end(fetch_gr) <- min(args[["limit_end"]],
end(range_gr) + width(range_gr))
# Assume diffList is a list
if(is.null(names(diffList))) {
listNames <- character(length(diffList))
for(i in seq_len(length(diffList))) {
listNames[i] <- paste(diffList[[i]], collapse = "_vs_")
}
names(diffList) <- listNames
} else {
listNames <- names(diffList)
for(i in seq_len(listNames)) {
if(listNames[i] == "") {
listNames[i] <- paste(diffList[[i]], collapse = "_vs_")
}
}
names(diffList) <- listNames
}
#
df_list <- list()
for(i in seq_len(length(diffList))) {
df <- c()
diffItem <- diffList[[i]]
diffName <- names(diffList)[i]
if(length(diffItem) == 2) {
if(
is.vector(diffItem, "character") &&
all(diffItem %in% names(x@norm_cov))
) {
# Assume these are track names
cov1 <- x@norm_cov[[diffItem[1]]]
cov2 <- x@norm_cov[[diffItem[2]]]
} else if(
is.vector(diffItem, "numeric") &&
all(diffItem %in% seq_len(length(x@norm_cov)))
) {
cov1 <- x@norm_cov[[diffItem[1]]]
cov2 <- x@norm_cov[[diffItem[2]]]
} else {
next # ignore this
}
} else {
next # ignore this
}
cov1 <- cov1[get("x") >= start(fetch_gr) & get("x") <= end(fetch_gr)]
cov2 <- cov2[get("x") >= start(fetch_gr) & get("x") <= end(fetch_gr)]
if(diff_stat == "t-test") {
df <- .pV_diffStat_ttest(cov1, cov2)
} else {
next
}
if(!is.null(df)) {
df$diffTrack <- diffName
df_list[[diffName]] <- df
}
}
df_all <- as.data.frame(rbindlist(df_list))
if(nrow(df_all) == 0) return(NULL)
# group annotation
df_all$info <- ""
# factor covTrack
df_all$diffTrack <- factor(df_all$diffTrack, names(diffList))
if(!usePlotly) {
p <- .pV_seed_ggplot_diffTrack(
plotViewStart, plotViewEnd,
df_all, diff_stat
)
} else {
p <- .pV_seed_plotly_diffTrack(
plotViewStart, plotViewEnd,
df_all, diff_stat
)
}
return(p)
}
################################################################################
# Creates a ggplot of coverage tracks
.pV_seed_ggplot_covTrack <- function(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
) {
covTracks <- levels(df_all$covTrack)
plotList <- list()
viewMin <- min(plotViewStart, plotViewEnd)
viewMax <- max(plotViewStart, plotViewEnd)
for(i in seq_len(length(covTracks))) {
covTrack = covTracks[i]
df <- df_all[df_all$covTrack == covTrack,]
dfJn <- dfJn_all[dfJn_all$covTrack == covTrack,]
nGroups <- length(unique(df$group))
# Filter by [plotViewStart, plotViewEnd]
df <- df[df$x >= viewMin & df$x <= viewMax,]
dfJn <- dfJn[dfJn$x >= viewMin & dfJn$x <= viewMax,]
p <- ggplot(df) + geom_hline(yintercept = 0)
# plot ribbon (group cov only)
if(plotMeanCov && rb %in% c("ci", "sd", "sem")) {
p <- p +
geom_ribbon(data = df, alpha = 0.2,
aes(
x = get("x"), y = get("mean"),
ymin = get("mean") - get(rb),
ymax = get("mean") + get(rb),
group = get("group"),
fill = get("group")
)
)
}
# plot line
if(plotMeanCov) {
p <- p + geom_line(aes(
x = get("x"), y = get("mean"),
group = get("group"), color = get("group")
))
} else {
p <- p + geom_line(aes(
x = get("x"), y = get("depth"),
group = get("group"), color = get("group")
))
}
if(nGroups == 1) {
p <- p + scale_color_manual(values = "black") +
scale_fill_manual(values = "black")
}
# plot junctions
if(plotJunctions & nGroups == 1) {
dfJn <- dfJn_all[dfJn_all$covTrack == covTrack,]
p <- p +
geom_line(
data = dfJn,
aes(x = get("x"), y = get("yarc"), group = get("info")),
color = "darkred"
)
# NB sashimi arc numbers are plotted later
}
if(nGroups > 1) {
p <- p + theme_white_legend +
theme(
axis.title.x = element_blank(),
legend.title = element_blank()
) +
labs(x = "", y = covTrack)
} else {
p <- p + theme_white +
theme(axis.title.x = element_blank()) +
labs(x = "", y = covTrack)
}
plotList[[i]] <- p
}
return(plotList)
}
# Creates plotly (data) for coverage tracks
.pV_seed_plotly_covTrack <- function(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
) {
covTrack <- levels(df_all$covTrack)
# nested list of data
dataList <- list()
# y axis range, title
layoutList <- list()
for(i in seq_len(length(covTrack))) {
track <- covTrack[i]
df <- df_all[df_all$covTrack == track,]
dfJn <- dfJn_all[dfJn_all$covTrack == track,]
dfJnSum <- dfJnSum_all[dfJnSum_all$covTrack == track,]
groupNames <- unique(df$group)
nGroups <- length(groupNames)
if(nrow(df) > 0) {
df$group <- factor(df$group, groupNames)
}
if(nrow(dfJn) > 0) {
dfJn$group <- factor(dfJn$group, groupNames)
}
if(nrow(dfJnSum) > 0) {
dfJnSum$group <- factor(dfJnSum$group, groupNames)
}
dataList[[i]] <- list()
dataList[[i]][[1]] <- list() # always have junction track
for(j in seq_len(nGroups)) {
dataList[[i]][[2*j]] <- list()
dataList[[i]][[2*j + 1]] <- list()
}
if(plotJunctions & nGroups == 1) {
# sashimi arcs
plotJuncData <- .plotly_makeJuncCurveData(
dfJn$x, dfJn$yarc, dfJn$info
)
dataList[[i]][[1]] <- list(
x = plotJuncData$x, y = plotJuncData$y,
text = plotJuncData$text
)
# sashimi numbers
dataList[[i]][[2]] <- list(
x = dfJnSum$xlabel, y = dfJnSum$ylabel,
text = dfJnSum$value
)
}
for(j in seq_len(nGroups)) {
df_sub <- as.data.frame(df[df$group == levels(df$group)[j],])
if("mean" %in% colnames(df)) {
dataList[[i]][[2*j + 1]] <- list(
x = df_sub$x, y = df_sub$mean, text = df_sub$info
)
dataList[[i]][[2*j + 1]][["hovertemplate"]] <- paste0(
"<b>", levels(df$group)[j], "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>Norm-Depth (mean):</i> %{y:.3f}<br>",
"<extra></extra>"
)
if(rb %in% c("ci", "sd", "sem")) {
dataList[[i]][[2*j + 2]] <- list(
x = c(df_sub$x, rev(df_sub$x)),
y = c(
df_sub$mean + unname(unlist(df_sub[,rb])),
rev(df_sub$mean - unname(unlist(df_sub[,rb])))
),
text = c(df_sub$info, rev(df_sub$info))
)
dataList[[i]][[2*j + 2]][["hovertemplate"]] <- paste0(
"<b>", levels(df$group)[j], "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>Norm-Depth (", rb,"):</i> %{y:.3f}<br>",
"<extra></extra>"
)
} else {
dataList[[i]][[2*j + 2]] <- list()
}
} else {
dataList[[i]][[2*j + 1]] <- list(
x = df_sub$x, y = df_sub$depth, text = df_sub$info
)
dataList[[i]][[2*j + 1]][["hovertemplate"]] <- paste0(
"<b>", levels(df$group)[j], "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>Depth:</i> %{y}<br>",
"<extra></extra>"
)
dataList[[i]][[2*j + 2]] <- list()
}
if(nGroups > 1) {
dataList[[i]][[2*j + 1]][["showlegend"]] <- TRUE
dataList[[i]][[2*j + 1]][["name"]] <- levels(df$group)[j]
}
}
# what is yrange?
df_sub <- df[
df$x >= min(plotViewStart, plotViewEnd) &
df$x <= max(plotViewStart, plotViewEnd),
]
if("mean" %in% colnames(df_sub)) {
ymax <- max(df_sub$mean)
} else {
ymax <- max(df_sub$depth)
}
layoutList[[i]] <- list(
yrange = c(0, ymax * 1.2),
title = track
)
}
return(list(
dataList = dataList,
layoutList = layoutList
))
}
# Creates a ggplot of differential tracks
.pV_seed_ggplot_diffTrack <- function(
plotViewStart, plotViewEnd, df_all, diff_stat
) {
diffTracks <- levels(df_all$diffTrack)
plotList <- list()
for(i in seq_len(length(diffTracks))) {
diffTrack = diffTracks[i]
df <- df_all[df_all$diffTrack == diffTrack,]
p <- ggplot(df) + geom_hline(yintercept = 0)
# plot line
p <- p + geom_line(aes(x = get("x"), y = get("stat")),
color = "black")
# plot y axis and format
p <- p + theme_white_legend + labs(y = "-log10 P")
p <- p +
theme(axis.title.x = element_blank()) +
labs(
x = "", y = diffTrack
)
plotList[[i]] <- p
}
return(plotList)
}
# Creates a plotly (data) of diff tracks
.pV_seed_plotly_diffTrack <- function(
plotViewStart, plotViewEnd, df_all, diff_stat
) {
diffTrack <- levels(df_all$diffTrack)
# nested list of data
dataList <- list()
# y axis range, title
layoutList <- list()
for(i in seq_len(length(diffTrack))) {
track <- diffTrack[i]
df <- df_all[df_all$diffTrack == track,]
dataList[[i]] <- list()
dataList[[i]][[1]] <- list(
x = df$x, y = df$stat, text = df$info
)
dataList[[i]][[1]][["hovertemplate"]] <- paste0(
"<b>", track, "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>-log10(P):</i> %{y}<br>",
"<extra></extra>"
)
df_sub <- df[
df$x >= min(plotViewStart, plotViewEnd) &
df$x <= max(plotViewStart, plotViewEnd),
]
ymax <- max(df_sub$stat)
layoutList[[i]] <- list(
yrange = c(0, ymax * 1.2),
title = track
)
}
return(list(
dataList = dataList,
layoutList = layoutList
))
}
################################################################################
# Plots an annotation track (or sub-track)
.pV_plotRef <- function(
DTplotlist, view_gr, showPlotRanges = GRanges(),
reverseGenomeCoords = FALSE,
add_information = TRUE,
usePlotly = FALSE,
instantPlotly = FALSE
) {
view_chr <- as.character(seqnames(view_gr))
view_start <- start(view_gr)
view_end <- end(view_gr)
group.DT <- copy(DTplotlist$group.DT)
# use one paste function
if(add_information) {
group.DT[, c("left_extra", "right_extra") := list("","")]
if(reverseGenomeCoords) {
group.DT[get("strand") == "-", c("right_extra") := " -->"]
group.DT[get("strand") == "+", c("left_extra") := "<-- "]
} else {
group.DT[get("strand") == "+", c("right_extra") := " -->"]
group.DT[get("strand") == "-", c("left_extra") := "<-- "]
}
group.DT[, c("display_name") := paste0(
get("left_extra"),
get("group_name"),
" (", get("strand"), ") ",
get("group_biotype"),
get("right_extra")
)]
group.DT[, c("left_extra", "right_extra") := list(NULL,NULL)]
group.DT[, c("disp_x") := 0.5 * (get("start") + get("end"))]
group.DT[get("start") < view_start & get("end") > view_start,
c("disp_x") := 0.5 * (view_start + get("end"))]
group.DT[get("end") > view_end & get("start") < view_end,
c("disp_x") := 0.5 * (get("start") + view_end)]
group.DT[get("start") < view_start & get("end") > view_end,
c("disp_x") := 0.5 * (view_start + view_end)]
}
reduced <- copy(DTplotlist$reduced.DT)
reduced <- reduced[!is.na(reduced$plot_level)]
condense_this <- DTplotlist$condense_this
exonRanges <- DTplotlist$exonRanges
# Hover Text annotation
reduced[, c("Information") := ""]
reduced[get("type") %in% c("exon", "CDS"), c("Information") := paste(
paste(get("transcript_id"), "exon", get("aux_id")),
paste0("(", get("feature_id"), ")"),
paste0(get("seqnames"), ":", get("start"), "-", get("end"), "/",
get("strand")),
sep = "\n"
)]
reduced[get("type") == "intron", c("Information") := paste(
get("feature_id"),
paste0(get("seqnames"), ":", get("start"), "-", get("end"), "/",
get("strand")),
sep = "\n"
)]
reduced <- as.data.frame(reduced)
reducedIntrons <- data.frame()
if (sum(reduced$type == "intron", na.rm = TRUE) > 0) {
reducedIntrons <- reduced[reduced$type == "intron", ]
reducedIntronsExpanded <- c()
for(i in seq_len(nrow(reducedIntrons))) {
reducedIntronsExpanded <- rbind(reducedIntronsExpanded, data.frame(
start = seq(reducedIntrons$start[i], reducedIntrons$end[i],
length.out = 10),
end = seq(reducedIntrons$start[i], reducedIntrons$end[i],
length.out = 10),
plot_level = reducedIntrons$plot_level[i],
highlight = reducedIntrons$highlight[i],
Information = reducedIntrons$Information[i]
))
}
col_highlights <- sort(unique(reducedIntronsExpanded$highlight))
}
nonIntrons <- reduced[reduced$type != "intron", ]
fill_highlights <- sort(unique(nonIntrons$highlight))
if(!reverseGenomeCoords) {
plotViewStart <- view_start
plotViewEnd <- view_end
} else {
plotViewStart <- view_end
plotViewEnd <- view_start
}
if(!usePlotly) {
p <- .pV_seed_ggplot_annoTrack(
plotViewStart, plotViewEnd, view_chr,
reducedIntronsExpanded, nonIntrons,
group.DT, condense_this, add_information,
showPlotRanges
)
} else if(instantPlotly) {
p <- .pV_seed_instant_plotly_annoTrack(
plotViewStart, plotViewEnd, view_chr,
reducedIntrons, nonIntrons,
group.DT, condense_this, add_information
)
} else {
p <- .pV_seed_plotly_annoTrack(
plotViewStart, plotViewEnd, view_chr,
reducedIntrons, nonIntrons,
group.DT, condense_this, add_information
)
}
return(p)
}
# Plots an annotation track showing only named exons
.pV_plotExons <- function(
DTplotlist, view_gr, exons_gr,
reverseGenomeCoords = FALSE,
add_information = TRUE,
usePlotly = FALSE,
instantPlotly = FALSE
) {
view_chr <- as.character(seqnames(view_gr))
view_start <- start(view_gr)
view_end <- end(view_gr)
exons <- as.data.frame(exons_gr)
if(!reverseGenomeCoords) {
plotViewStart <- view_start
plotViewEnd <- view_end
} else {
plotViewStart <- view_end
plotViewEnd <- view_start
}
if(!usePlotly) {
p <- .pV_seed_ggplot_HLexons(
plotViewStart, plotViewEnd, view_chr, exons
)
} else if(instantPlotly) {
p <- .pV_seed_instant_plotly_HLexons(
plotViewStart, plotViewEnd, view_chr, exons
)
} else {
p <- .pV_seed_plotly_HLexons(
plotViewStart, plotViewEnd, view_chr, exons
)
}
return(p)
}
# Convert highlight codes to colors
.pV_highlight_to_colors <- function(highlights, usePlotly = FALSE) {
if(is.null(highlights)) return(NULL)
if(usePlotly) {
highlights <- sub("0", "rgba(0,0,0,1)", highlights)
highlights <- sub("1", "rgba(0,0,255,1)", highlights)
highlights <- sub("2", "rgba(255,0,0,1)", highlights)
highlights <- sub("3", "rgba(255,0,255,1)", highlights)
} else {
highlights <- sub("0", "black", highlights)
highlights <- sub("1", "blue", highlights)
highlights <- sub("2", "red", highlights)
highlights <- sub("3", "purple", highlights)
}
return(highlights)
}
# Creates a ggplot of annotation track
.pV_seed_ggplot_annoTrack <- function(
plotViewStart, plotViewEnd, view_chr,
introns, exons, group.DT, condense_this, add_information,
showPlotRanges
) {
col_highlights <- sort(unique(introns$highlight))
fill_highlights <- sort(unique(exons$highlight))
col_highlights <- .pV_highlight_to_colors(col_highlights)
fill_highlights <- .pV_highlight_to_colors(fill_highlights)
anno <- NULL
if(add_information) {
if (condense_this == TRUE) {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.5 + 0.3 *
runif(rep(1, nrow(group.DT))),
Information = group.DT$display_name
)
} else {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.4,
Information = group.DT$display_name
)
}
}
if (nrow(group.DT) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(group.DT$plot_level)
}
p <- ggplot()
if(nrow(introns) > 0) {
p <- p + geom_line(data = introns,
aes(x = get("start"), y = get("plot_level"),
color = get("highlight"), group = get("Information"))
)
}
if(nrow(exons) > 0) {
p <- p + geom_rect(data = exons,
aes(xmin = get("start"), xmax = get("end"),
ymin = get("plot_level") - 0.1 -
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
ymax = get("plot_level") + 0.1 +
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
fill = get("highlight")
)
)
}
if(length(showPlotRanges) > 0) {
ref_ymin <- min(layer_scales(p)$y$range$range)
ref_ymax <- max(layer_scales(p)$y$range$range)
dfRanges <- as.data.frame(showPlotRanges)
dfRanges$rangeName <- as.character(seq_len(length(showPlotRanges)))
p <- p + geom_rect(data = dfRanges,
mapping = aes(xmin = get("start"), xmax = get("end"),
ymin = ref_ymin - 0.5, ymax = ref_ymax + 0.5
), fill = NA, color = "red"
) + geom_text(data = dfRanges,
aes(
x = 0.5 * (get("start") + get("end")),
y = ref_ymin - 1, label = get("rangeName")
)
)
}
if(!is.null(col_highlights))
p <- p + scale_color_manual(values = col_highlights)
if(!is.null(fill_highlights))
p <- p + scale_fill_manual(values = fill_highlights)
p <- p + theme_white_legend_plot_track +
theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank()
)
if(!is.null(anno)) {
p <- p + geom_text(
data = anno,
aes(x = get("x"), y = get("y"), label = get("Information"))
)
}
p <- p + theme(legend.position = "none") +
labs(x = paste("Chromosome", view_chr))
ref_ymin <- min(layer_scales(p)$y$range$range)
ref_ymax <- max(layer_scales(p)$y$range$range)
p <- p +
scale_x_continuous(labels = label_number(scale_cut = cut_si(""))) +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(ref_ymin - 1, ref_ymax + 1),
expand = FALSE
)
return(p)
}
# Creates a ggplot of exons track
.pV_seed_ggplot_HLexons <- function(
plotViewStart, plotViewEnd, view_chr, exons
) {
fill_highlights <- sort(unique(exons$highlight))
fill_highlights <- .pV_highlight_to_colors(fill_highlights)
if (nrow(exons) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(exons$plot_level)
}
p <- ggplot()
if(nrow(exons) > 0) {
p <- p + geom_rect(data = exons,
aes(xmin = get("start"), xmax = get("end"),
ymin = get("plot_level") - 0.1 -
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
ymax = get("plot_level") + 0.1 +
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
fill = get("highlight")
)
) + geom_text(
data = exons,
aes(
x = 0.5 * (get("start") + get("end")),
y = get("plot_level") - 0.4,
label = get("exon_name")
)# ,
# hjust=1, angle=45
)
}
p <- p +
scale_fill_manual(values = fill_highlights)
p <- p + theme_white_legend_plot_track +
theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank()
)
p <- p + theme(legend.position = "none") +
labs(x = paste("Chromosome", view_chr))
ref_ymin <- min(layer_scales(p)$y$range$range)
ref_ymax <- max(layer_scales(p)$y$range$range)
p <- p +
scale_x_continuous(labels = label_number(scale_cut = cut_si(""))) +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(ref_ymin - 1, ref_ymax + 1),
expand = FALSE
)
return(p)
}
# Creates an instant plotly of annotation track
.pV_seed_instant_plotly_annoTrack <- function(
plotViewStart, plotViewEnd, view_chr,
introns, exons, group.DT, condense_this, add_information
) {
col_highlights <- sort(unique(introns$highlight))
fill_highlights <- sort(unique(exons$highlight))
col_colors <- .pV_highlight_to_colors(col_highlights)
fill_colors <- .pV_highlight_to_colors(fill_highlights)
anno <- NULL
if(add_information) {
if (condense_this == TRUE) {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.5 + 0.3 *
runif(rep(1, nrow(group.DT))),
Information = group.DT$display_name
)
} else {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.4,
Information = group.DT$display_name
)
}
}
if (nrow(group.DT) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(group.DT$plot_level)
}
fig <- plot_ly()
if(nrow(introns) > 0) {
intronLineData <- .plotly_makeLineData(
introns$start, introns$plot_level, introns$end, introns$plot_level,
introns$Information, introns$highlight
)
for(i in seq_len(length(col_highlights))) {
hl <- col_highlights[i]
color <- col_colors[i]
fig <- fig %>% add_trace(
data = intronLineData[intronLineData$colorInfo == hl,],
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~text"),
hoveron = "points", hoverinfo = 'text',
line = list(
color = color
),
type = 'scatter', mode = 'lines',
showlegend = FALSE
)
}
}
typeThick <- c("CDS", "start_codon", "stop_codon")
if(nrow(exons) > 0) {
exons_thin <- exons[!(exons$type %in% typeThick), ]
exons_thick <- exons[exons$type %in% typeThick, ]
exons_thin$ymin <- exons_thin$plot_level - 0.15
exons_thin$ymax <- exons_thin$plot_level + 0.15
exons_thick$ymin <- exons_thick$plot_level - 0.3
exons_thick$ymax <- exons_thick$plot_level + 0.3
exons_comb <- rbind(exons_thin, exons_thick)
exons_instr <- .plotly_makeRectData(
exons_comb$start, exons_comb$end,
exons_comb$ymin, exons_comb$ymax,
exons_comb$Information,
exons_comb$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
fig <- fig %>% add_trace(
data = rects,
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~text"),
type = 'scatter', mode = 'lines',
hoveron = "points+fills", hoverinfo = 'text',
line = list(
color = "transparent"
),
fill = "toself",
fillcolor = fill_colors[i],
showlegend = FALSE
)
}
}
if(add_information) {
xtitle <- paste("Chromosome/Scaffold", view_chr)
fig <- fig %>% add_trace(
data = anno,
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~Information"),
hoverinfo = 'text',
type = 'scatter', mode = 'text',
textposition = 'middle',
showlegend = FALSE
)
} else {
xtitle <- ""
}
fig <- fig %>% layout(
dragmode = "pan",
xaxis = list(
range = c(plotViewStart, plotViewEnd),
title = xtitle
),
yaxis = list(
title = "",
range = c(0, 1 + max_plot_level)# ,
# fixedrange = TRUE
)
)
return(fig)
}
# Creates a plotly of exons track
.pV_seed_instant_plotly_HLexons <- function(
plotViewStart, plotViewEnd, view_chr, exons
) {
fill_highlights <- sort(unique(exons$highlight))
fill_colors <- .pV_highlight_to_colors(fill_highlights)
if (nrow(exons) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(exons$plot_level)
}
fig <- plot_ly()
if(nrow(exons) > 0) {
exons$ymin <- exons$plot_level - 0.15
exons$ymax <- exons$plot_level + 0.15
exons$xdisp <- 0.5 * (exons$start + exons$end)
exons$ydisp <- exons$plot_level - 0.4
exons_instr <- .plotly_makeRectData(
exons$start, exons$end,
exons$ymin, exons$ymax,
exons$exon_name,
exons$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
fig <- fig %>% add_trace(
data = rects,
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~text"),
type = 'scatter', mode = 'lines',
hoveron = "points+fills", hoverinfo = 'text',
line = list(
color = "transparent"
),
fill = "toself",
fillcolor = fill_colors[i],
showlegend = FALSE
)
fig <- fig %>% add_trace(
data = exons,
x = as.formula("~xdisp"),
y = as.formula("~ydisp"),
text = as.formula("~exon_name"),
hoverinfo = 'text',
type = 'scatter', mode = 'text',
textposition = 'middle',
showlegend = FALSE
)
}
}
xtitle <- paste("Chromosome/Scaffold", view_chr)
fig <- fig %>% layout(
dragmode = "pan",
xaxis = list(
range = c(plotViewStart, plotViewEnd),
title = xtitle
),
yaxis = list(
title = "",
range = c(0, 1 + max_plot_level),
fixedrange = TRUE
)
)
return(fig)
}
# Creates a plotly (data) of annotation track
.pV_seed_plotly_annoTrack <- function(
plotViewStart, plotViewEnd, view_chr,
introns, exons, group.DT, condense_this, add_information
) {
dataList <- list()
# x,y axis range, title
layoutList <- list()
col_highlights <- fill_highlights <- c("0", "1", "2", "3")
if(nrow(introns) > 0) {
intronLineData <- .plotly_makeLineData(
introns$start, introns$plot_level, introns$end, introns$plot_level,
introns$Information, introns$highlight
)
for(i in seq_len(length(col_highlights))) {
dataList[[i]] <- list()
hl <- col_highlights[i]
data_df <- intronLineData[intronLineData$colorInfo == hl,]
if(nrow(data_df) > 0) {
dataList[[i]] <- list(
x = data_df$x, y = data_df$y, text = data_df$text
)
}
}
} else {
for(i in seq_len(length(col_highlights))) {
dataList[[i]] <- list()
}
}
typeThick <- c("CDS", "start_codon", "stop_codon")
if(nrow(exons) > 0) {
exons_thin <- exons[!(exons$type %in% typeThick), ]
exons_thick <- exons[exons$type %in% typeThick, ]
exons_thin$ymin <- exons_thin$plot_level - 0.15
exons_thin$ymax <- exons_thin$plot_level + 0.15
exons_thick$ymin <- exons_thick$plot_level - 0.3
exons_thick$ymax <- exons_thick$plot_level + 0.3
exons_comb <- rbind(exons_thin, exons_thick)
exons_instr <- .plotly_makeRectData(
exons_comb$start, exons_comb$end,
exons_comb$ymin, exons_comb$ymax,
exons_comb$Information,
exons_comb$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
if(nrow(rects) > 0) {
dataList[[i + 4]] <- list(
x = rects$x, y = rects$y, text = rects$text
)
}
}
} else {
for(i in seq_len(length(col_highlights))) {
dataList[[i + 4]] <- list()
}
}
if(add_information) {
if (condense_this == TRUE) {
dataList[[9]] <- list(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.5 + 0.3 *
runif(rep(1, nrow(group.DT))),
text = group.DT$display_name
)
} else {
dataList[[9]] <- list(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.4,
text = group.DT$display_name
)
}
} else {
dataList[[9]] <- list()
}
if(add_information) {
xtitle <- paste("Chromosome/Scaffold", view_chr)
} else {
xtitle <- ""
}
if (nrow(group.DT) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(group.DT$plot_level)
}
layoutList <- list(
xrange = c(plotViewStart, plotViewEnd),
xtitle = xtitle,
yrange = c(0, 1 + max_plot_level)
)
return(list(
dataList = dataList,
layoutList = layoutList
))
}
# Creates a plotly (data) of exons track
.pV_seed_plotly_HLexons <- function(
plotViewStart, plotViewEnd, view_chr, exons
) {
dataList <- list()
# x,y axis range, title
layoutList <- list()
col_highlights <- fill_highlights <- c("0", "1", "2", "3")
for(i in seq_len(length(col_highlights))) {
dataList[[i]] <- list()
}
typeThick <- c("CDS", "start_codon", "stop_codon")
if(nrow(exons) > 0) {
exons$ymin <- exons$plot_level - 0.15
exons$ymax <- exons$plot_level + 0.15
exons$xdisp <- 0.5 * (exons$start + exons$end)
exons$ydisp <- exons$plot_level - 0.4
exons_instr <- .plotly_makeRectData(
exons$start, exons$end,
exons$ymin, exons$ymax,
exons$exon_name,
exons$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
if(nrow(rects) > 0) {
dataList[[i + 4]] <- list(
x = rects$x, y = rects$y, text = rects$text
)
}
}
} else {
for(i in seq_len(length(col_highlights))) {
dataList[[i + 4]] <- list()
}
}
dataList[[9]] <- list(
x = exons$xdisp,
y = exons$ydisp,
text = exons$exon_name
)
xtitle <- paste("Chromosome/Scaffold", view_chr)
if (nrow(exons) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(exons$plot_level)
}
layoutList <- list(
xrange = c(plotViewStart, plotViewEnd),
xtitle = xtitle,
yrange = c(0, 1 + max_plot_level)
)
return(list(
dataList = dataList,
layoutList = layoutList
))
}
################################################################################
.seq2 <- Vectorize(seq.default, vectorize.args = c("from", "to", "length.out"))
# Creates arcs for sashimi plot
.pV_sashimiArc <- function(
junc,
arcHeight = 0,
junctionThreshold = 0
) {
if(is.null(junc)) return(NA)
DT <- data.table(
juncName = paste0(
seqnames(junc), ":",
start(junc), "-", end(junc), "/", strand(junc)
),
leftX = start(junc) - 1,
rightX = end(junc) + 1,
leftY = mcols(junc)$leftCoordHeight,
rightY = mcols(junc)$rightCoordHeight
)
isMean <- ("mean" %in% names(mcols(junc)))
if(isMean) {
DT$count <- mcols(junc)$mean
DT$countsd <- mcols(junc)$sd
DT$value <- paste0(
round(100 * DT$count, 1), "+/-",
round(100 * DT$countsd, 1), " %"
)
DT$info <- paste(
paste0("Junction: ", DT$juncName),
paste0("PSI: ", DT$value),
sep = "\n"
)
} else {
DT$count <- mcols(junc)$count
DT$countsd <- 0
DT$value <- round(DT$count,2)
DT$info <- paste(
paste0("Junction: ", DT$juncName),
paste0("Depth: ", DT$count),
sep = "\n"
)
}
DT <- DT[get("count") >= junctionThreshold]
final_DT <- data.table(
coord = rep(DT$juncName, each = 90),
value = rep(DT$value, each = 90),
info = rep(DT$info, each = 90),
x = as.vector(.seq2(DT$leftX, DT$rightX, length.out = 90)),
ytmp = as.vector(.seq2(DT$leftY, DT$rightY, length.out = 90)),
yarc = as.vector(.seq2(DT$leftY, DT$rightY, length.out = 90)) +
rep(arcHeight * sinpi(seq(0,1,length.out = 90)), nrow(DT))
)
return(final_DT)
}
################################################################################
# Check requested tracks are legitimate in getPlotObject
.gPO_check_tracks <- function(
obj, args
) {
# sanity check condition / tracks against colData
colData <- obj@colData
sampleList <- list()
samples_to_check <- c()
# tracks must be a character vector
if(!is.vector(args[["tracks"]], "character")) {
.log(paste(
"In getPlotObject,",
"`tracks` must be a character vector"
))
}
# remove invalid tracks
invalid_tracks <- c()
if("condition" %in% names(args)) {
if(args[["condition"]] %in% colnames(colData)) {
for(cond_elem in args[["tracks"]]) {
# cond_elem is of length 1
if(
any(colData[, args[["condition"]]] == cond_elem)
) {
sampleList[[cond_elem]] <- rownames(colData)[
colData[, args[["condition"]]] == cond_elem]
samples_to_check <- c(samples_to_check,
sampleList[[cond_elem]])
} else {
# ignore track altogether
invalid_tracks <- c(invalid_tracks, cond_elem)
}
}
} else {
.log(paste("In getPlotObject,",
"condition is not a recognised column in colData"
))
}
} else {
for(i in seq_len(length(args[["tracks"]]))) {
for(samp in args[["tracks"]]) {
if(samp %in% rownames(colData)) {
sampleList[[samp]] <- samp
samples_to_check <- c(samples_to_check, samp)
} else {
invalid_tracks <- c(invalid_tracks, samp)
}
}
}
}
# Remove invalid tracks
if(length(invalid_tracks) == length(args[["tracks"]])) {
.log(paste(
"In getPlotObject,",
"all elements in `tracks` are invalid"
))
} else if(length(invalid_tracks) > 0) {
args[["tracks"]] <- args[["tracks"]][
!(args[["tracks"]] %in% invalid_tracks)
]
}
args[["sampleList"]] <- sampleList
args[["samples_to_check"]] <- samples_to_check
return(args)
}
# Check that getPlotObject knows which stranded data to retrieve
.gPO_check_strand <- function(
obj, args
) {
args[["view_strand_jn"]] <- args[["strand"]]
samples_to_check <- args[["samples_to_check"]]
if(args[["strand"]] != "*") {
if(all(obj@normData[["sampleStrand"]][samples_to_check] == -1)) {
if(args[["strand"]] == "+") {
args[["view_strand"]] <- "-"
} else {
args[["view_strand"]] <- "+"
}
}
} else {
args[["view_strand"]] <- "*"
}
return(args)
}
################################################################################
# Get normalized coverage from obj
.gPO_getCoverage <- function(obj, args, normalize = FALSE) {
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand <- args[["view_strand"]]
Event <- args[["Event"]]
outList <- list()
sampleList <- args[["sampleList"]]
if(view_strand == "+") {
cov <- obj@covData$pos
} else if(view_strand == "-") {
cov <- obj@covData$neg
} else {
cov <- obj@covData$uns
}
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
matList <- list()
for(s in samples) {
view <- IRanges::Views(cov[[s]], view_start, view_end)
matList[[s]] <- as.matrix(view[[1]])
}
mat <- do.call(cbind, matList)
colnames(mat) <- samples
# Normalize coverage
if(normalize) {
normFactor <- unlist(obj@normData$norms[Event,samples])
# discard samples with low normalization factors
which_lowNorm <- which(normFactor < 5)
if(length(which_lowNorm) == length(normFactor)) {
outList[[track]] <- data.frame(x = seq(view_start, view_end))
next
} else if(length(which_lowNorm) > 0) {
which_Norm <- which(normFactor > 5)
mat <- mat[, which_Norm, drop = FALSE]
normFactor <- normFactor[which_Norm]
}
mat <- t(t(mat) / normFactor)
}
outList[[track]] <- data.table(x = seq(view_start, view_end))
outList[[track]] <- cbind(outList[[track]], as.data.table(mat))
if(length(samples) == 1) {
colnames(outList[[track]])[2] <- "depth"
}
}
return(outList)
}
# Calculates mean and se/sem/ci of normalized coverages per group
.gPO_covStats <- function(cov) {
out <- list()
outnames <- names(cov)
conf.int <- 0.95
for(i in seq_len(length(cov))) {
covmat <- as.matrix(cov[[outnames[i]]][,-1])
out[[outnames[i]]] <- data.table(
x = cov[[outnames[i]]]$x,
mean = rowMeans(covmat)
)
n <- ncol(covmat)
out[[outnames[i]]]$sd <- rowSds(covmat)
out[[outnames[i]]]$sem <- out[[outnames[i]]]$sd / sqrt(n)
out[[outnames[i]]]$ci <- qt((1 + conf.int) / 2, df = n - 1) *
out[[outnames[i]]]$sem
}
return(out)
}
# Compiles per-condition junction stats (normalized junction counts)
.gPO_normJuncStats <- function(obj, args, cov_stats) {
view_chr <- args[["view_chr"]]
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand_jn <- args[["view_strand_jn"]]
outList <- list()
sampleList <- args[["sampleList"]]
gr_select <- GRanges(
view_chr,
IRanges(view_start, view_end),
view_strand_jn
)
OL <- findOverlaps(obj@juncData$junc_gr, gr_select)
if(length(unique(from(OL))) == 0) return(list())
hits <- sort(unique(from(OL)))
junc_gr <- obj@juncData$junc_gr[hits]
leftCoords <- start(junc_gr) - 1
rightCoords <- end(junc_gr) + 1
lC_inrange <- leftCoords >= view_start & leftCoords <= view_end
rC_inrange <- rightCoords >= view_start & rightCoords <= view_end
Event <- args[["Event"]]
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
if(view_strand_jn == "*") {
junc <- obj@juncData$junc_counts_uns[
hits, samples, drop = FALSE]
} else {
junc <- obj@juncData$junc_counts[
hits, samples, drop = FALSE]
}
cov <- cov_stats[[track]][, seq_len(2), with = FALSE]
normFactor <- unlist(obj@normData$norms[Event,samples])
junc <- t(t(as.matrix(junc)) / normFactor)
# Returned data
outList[[track]] <- .grDT(cbind(
as.data.frame(junc_gr),
data.frame(
mean = rowMeans(junc),
sd = rowSds(junc),
leftCoordHeight = max(cov[,2]),
rightCoordHeight = max(cov[,2])
)
), keep.extra.columns = TRUE)
if(sum(lC_inrange) > 0) {
covL <- unname(unlist(
cov[leftCoords[lC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$leftCoordHeight[lC_inrange] <- covL
}
if(sum(rC_inrange) > 0) {
covR <- unname(unlist(
cov[rightCoords[rC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$rightCoordHeight[rC_inrange] <- covR
}
}
return(outList)
}
.pV_diffStat_ttest <- function(cov1, cov2) {
coords <- sort(intersect(cov1$x, cov2$x))
mat <- cbind(
as.matrix(cov1[get("x") %in% coords, -1]),
as.matrix(cov2[get("x") %in% coords, -1])
)
fac <- factor(
c(
rep("1", ncol(cov1) - 1),
rep("2", ncol(cov2) - 1)
)
)
t_test <- genefilter::rowttests(mat, fac)
ret <- data.frame(
x = coords,
stat = -log10(t_test$p.value)
)
ret$stat[!is.finite(ret$stat)] <- 0
return(ret)
}
.gPO_PSIjuncStats <- function(obj, args, norm_cov) {
view_chr <- args[["view_chr"]]
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand_jn <- args[["view_strand_jn"]]
outList <- list()
sampleList <- args[["sampleList"]]
gr_select <- GRanges(
view_chr,
IRanges(view_start, view_end),
view_strand_jn
)
OL <- findOverlaps(obj@juncData$junc_gr, gr_select)
if(length(unique(from(OL))) == 0) return(list())
hits <- sort(unique(from(OL)))
junc_gr <- obj@juncData$junc_gr[hits]
leftCoords <- start(junc_gr) - 1
rightCoords <- end(junc_gr) + 1
lC_inrange <- leftCoords >= view_start & leftCoords <= view_end
rC_inrange <- rightCoords >= view_start & rightCoords <= view_end
# normalize by starting coordinates:
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
junc <- obj@juncData$junc_PSI[hits, samples, drop = FALSE]
# Returned data
outList[[track]] <- .grDT(cbind(
as.data.frame(junc_gr),
data.frame(
mean = rowMeans(junc),
sd = rowSds(junc),
leftCoordHeight = 1,
rightCoordHeight = 1
)
), keep.extra.columns = TRUE)
cov <- norm_cov[[track]]
if(sum(lC_inrange) > 0) {
covL <- cov[leftCoords[lC_inrange] - cov$x[1] + 1, -1]
mcols(outList[[track]])$leftCoordHeight[lC_inrange] <- rowMeans(covL)
}
if(sum(rC_inrange) > 0) {
covR <- cov[rightCoords[rC_inrange] - cov$x[1] + 1, -1]
mcols(outList[[track]])$rightCoordHeight[rC_inrange] <- rowMeans(covR)
}
}
return(outList)
}
.gPO_getJunc <- function(obj, args, raw_cov, normalized = FALSE) {
view_chr <- args[["view_chr"]]
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand_jn <- args[["view_strand_jn"]]
outList <- list()
sampleList <- args[["sampleList"]]
gr_select <- GRanges(
view_chr,
IRanges(view_start, view_end),
view_strand_jn
)
OL <- findOverlaps(obj@juncData$junc_gr, gr_select)
if(length(unique(from(OL))) == 0) return(list())
hits <- sort(unique(from(OL)))
junc_gr <- obj@juncData$junc_gr[hits]
leftCoords <- start(junc_gr) - 1
rightCoords <- end(junc_gr) + 1
lC_inrange <- leftCoords >= view_start & leftCoords <= view_end
rC_inrange <- rightCoords >= view_start & rightCoords <= view_end
Event <- NULL
if(normalized) Event <- args[["Event"]]
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
# Only do it if cov is a single track
if(length(samples) == 1) {
if(view_strand_jn == "*") {
junc <- obj@juncData$junc_counts_uns[
hits, samples, drop = FALSE]
} else {
junc <- obj@juncData$junc_counts[
hits, samples, drop = FALSE]
}
cov <- raw_cov[[track]]
# Normalize
if(normalized) {
normFactor <- unlist(obj@normData$norms[Event,samples])
junc <- junc / normFactor
}
# Returned data
outList[[track]] <- .grDT(cbind(
as.data.frame(junc_gr),
data.frame(
count = unname(unlist(junc)),
leftCoordHeight = max(cov[,2]),
rightCoordHeight = max(cov[,2])
)
), keep.extra.columns = TRUE)
if(sum(lC_inrange) > 0) {
covL <- unname(unlist(
cov[leftCoords[lC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$leftCoordHeight[lC_inrange] <- covL
}
if(sum(rC_inrange) > 0) {
covR <- unname(unlist(
cov[rightCoords[rC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$rightCoordHeight[rC_inrange] <- covR
}
} else {
# empty
outList[[track]] <- GRanges()
}
}
return(outList)
}
# list checker for plotViews
.pV_trackList_from_vector <- function(trackVec, args) {
trackList <- list()
if(is.vector(trackVec, "numeric")) {
idx <- trackVec
if(all(idx %in% seq_len(length(args[["tracks"]])))) {
trackList <- lapply(idx, function(i) i)
} else if(any(idx %in% seq_len(length(args[["tracks"]])))) {
idx <- idx[idx %in% seq_len(length(args[["tracks"]]))]
trackList <- lapply(idx, function(i) i)
} else {
.log(paste(
"In plotView,",
"`trackList` contains indices outside the range of tracks",
"found in tracks(covPlotObject)"
))
}
} else if(is.vector(trackVec, "character")) {
idxNames <- trackVec
if(all(idxNames %in% args[["tracks"]])) {
trackList <- lapply(idxNames, function(i) {
which(args[["tracks"]] == i)[1]
})
} else if(any(idxNames %in% args[["tracks"]])) {
idxNames <- idxNames[idxNames %in% args[["tracks"]]]
trackList <- lapply(idxNames, function(i) {
which(args[["tracks"]] == i)[1]
})
} else {
.log(paste(
"In plotView,",
"`trackList` contains tracks",
"not found in tracks(covPlotObject)"
))
}
}
if(!is.null(names(trackVec))) {
names(trackList) <- names(trackVec)
if(any(names(trackList) == "")) {
which_noname <- which(names(trackList) == "")
names(trackList)[which_noname] <- paste(
"Track", as.character(which_noname), sep = "_"
)
}
} else {
for(i in seq_len(length(trackList))) {
elem <- trackList[[i]]
elemNames <- args[["tracks"]][elem]
names(trackList)[i] <- paste(elemNames, collapse = ",")
}
}
return(trackList)
}
.pV_trackList_from_list <- function(trackList, args) {
outList <- list()
for(z in seq_len(length(trackList))) {
trackVec <- trackList[[z]]
if(is.vector(trackVec, "numeric")) {
idx <- trackVec
if(all(idx %in% seq_len(length(args[["tracks"]])))) {
outList[[z]] <- idx
} else if(any(idx %in% seq_len(length(args[["tracks"]])))) {
idx <- idx[idx %in% seq_len(length(args[["tracks"]]))]
outList[[z]] <- idx
} else {
.log(paste(
"In plotView,",
"`trackList` contains indices outside the range of tracks",
"found in tracks(covPlotObject)"
))
}
} else if(is.vector(trackVec, "character")) {
idxNames <- trackVec
if(all(idxNames %in% args[["tracks"]])) {
outList[[z]] <- which(args[["tracks"]] %in% idxNames)
} else if(any(idxNames %in% args[["tracks"]])) {
idxNames <- idxNames[idxNames %in% args[["tracks"]]]
outList[[z]] <- which(args[["tracks"]] %in% idxNames)
} else {
.log(paste(
"In plotView,",
"`trackList` contains tracks",
"not found in tracks(covPlotObject)"
))
}
}
}
# transfer names
if(!is.null(names(trackList))) {
names(outList) <- names(trackList)
if(any(names(outList) == "")) {
which_noname <- which(names(outList) == "")
names(outList)[which_noname] <- paste(
"Track", as.character(which_noname), sep = "_"
)
}
} else {
for(i in seq_len(length(outList))) {
elem <- outList[[i]]
elemNames <- args[["tracks"]][elem]
names(outList)[i] <- paste(elemNames, collapse = ",")
}
}
return(outList)
}
# Make plotly x/y/text data given yarc data
.plotly_makeJuncCurveData <- function(
x, yarc, info
) {
sourceDT <- data.table(
x = x, y = yarc, info = info
)
setorderv(c("info", "x"))
infoVec <- unique(sourceDT$info)
df_list <- list()
for(i in seq_len(length(infoVec))) {
subDT <- sourceDT[get("info") == infoVec[i]]
df_list[[2 * i - 1]] <- data.table(
x = subDT$x,
y = subDT$y,
text = infoVec[i]
)
if(i < length(info)) {
df_list[[2 * i]] <- data.table(
x = NA, y = NA, text = NA
)
}
}
return(as.data.frame(rbindlist(df_list)))
}
# Make plotly line data based on line info
.plotly_makeLineData <- function(
xstart, ystart, xend, yend, info, colorInfo, resolution = 20
) {
df_list <- list()
for(i in seq_len(length(info))) {
df_list[[2 * i - 1]] <- data.table(
x = seq(xstart[i], xend[i], length.out = resolution),
y = seq(ystart[i], yend[i], length.out = resolution),
text = info[i],
colorInfo = colorInfo[i]
)
if(i < length(info)) {
df_list[[2 * i]] <- data.table(
x = NA, y = NA, text = NA, colorInfo = colorInfo[i]
)
}
}
return(as.data.frame(rbindlist(df_list)))
}
# Make plotly rectangle data based on x/y min/max info
.plotly_makeRectData <- function(
xmin, xmax, ymin, ymax, info, colorInfo, resolution = 20
) {
df_list <- list()
if(resolution < 1) resolution <- 1
for(i in seq_len(length(info))) {
df_list[[5 * i - 4]] <- data.table(
x = seq(xmin[i], xmax[i], length.out = resolution + 1)[-1],
y = ymin[i],
text = info[i],
colorInfo = colorInfo[i]
)
df_list[[5 * i - 3]] <- data.table(
x = xmax[i],
y = seq(ymin[i], ymax[i], length.out = resolution + 1)[-1],
text = info[i],
colorInfo = colorInfo[i]
)
df_list[[5 * i - 2]] <- data.table(
x = seq(xmax[i], xmin[i], length.out = resolution + 1)[-1],
y = ymax[i],
text = info[i],
colorInfo = colorInfo[i]
)
df_list[[5 * i - 1]] <- data.table(
x = xmin[i],
y = seq(ymax[i], ymin[i], length.out = resolution + 1),
text = info[i],
colorInfo = colorInfo[i]
)
if(i < length(info)) {
df_list[[5 * i]] <- data.table(
x = NA, y = NA, text = NA, colorInfo = colorInfo[i]
)
}
}
return(as.data.frame(rbindlist(df_list)))
}
# Assemble plotly given seed plotly object and x/y/info data
.pV_assemble_covPlotly2 <- function(p, vLayout) {
numCovTraces <- numDiffTraces <- c()
if(length(p@covTrack) > 0) {
numCovTraces <- vapply(p@covTrack[[1]][["dataList"]], length, 0)
numCovTraces <- (numCovTraces - 2)/2
}
if(length(p@diffTrack) > 0) {
numDiffTraces <- vapply(p@diffTrack[[1]][["dataList"]], length, 0)
}
p <- .knitPlotly(p, numCovTraces, numDiffTraces, vLayout)
# inject plot data
for(i in seq_len(length(numCovTraces))) {
p <- .injectPlotData(p, p@args[["covTrackPos"]][i],
p@covTrack[[1]][["dataList"]][[i]],
p@covTrack[[1]][["layoutList"]][[i]][["title"]]
)
}
for(i in seq_len(length(numDiffTraces))) {
p <- .injectPlotData(p, p@args[["diffTrackPos"]][i],
p@diffTrack[[1]][["dataList"]][[i]],
p@diffTrack[[1]][["layoutList"]][[i]][["title"]]
)
}
p <- .injectPlotData(p, p@args[["annoTrackPos"]],
p@annoTrack[[1]][["dataList"]],
p@annoTrack[[1]][["layoutList"]][["xtitle"]]
)
# Layout modifications
# x axis
if("xtitle" %in% names(p@annoTrack[[1]][["layoutList"]])) {
p <- .adjustXtitle(p, p@annoTrack[[1]][["layoutList"]][["xtitle"]])
}
if("xrange" %in% names(p@annoTrack[[1]][["layoutList"]])) {
p <- .adjustXrange(p, p@annoTrack[[1]][["layoutList"]][["xrange"]])
}
# y axis
axisNum <- 0
for(i in seq_len(length(numCovTraces))) {
axisNum <- axisNum + 1
p <- .fixYrange(p, axisNum)
if("yrange" %in% names(p@covTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYrange(p,
p@covTrack[[1]][["layoutList"]][[i]][["yrange"]],
axisNum
)
}
if("title" %in% names(p@covTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYtitle(p,
p@covTrack[[1]][["layoutList"]][[i]][["title"]],
axisNum
)
}
}
for(i in seq_len(length(numDiffTraces))) {
axisNum <- axisNum + 1
p <- .fixYrange(p, axisNum)
if("yrange" %in% names(p@diffTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYrange(p,
p@diffTrack[[1]][["layoutList"]][[i]][["yrange"]],
axisNum
)
}
if("title" %in% names(p@diffTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYtitle(p,
p@diffTrack[[1]][["layoutList"]][[i]][["title"]],
axisNum
)
}
}
axisNum <- axisNum + 1
p <- .adjustYtitle(p, "", axisNum)
if("yrange" %in% names(p@covTrack[["layoutList"]][[i]])) {
p <- .adjustYrange(p, p@covTrack[["layoutList"]][[i]][["yrange"]],
axisNum)
}
return(p)
}
# Assemble ggplot given sub-track data
.pV_assemble_ggplot <- function(
covTrack, diffTrack, annoSubTrack, annoFullTrack,
verticalLayout, horizontalLayout
) {
# Correct vertical layout for multiple cov/diff tracks
nCovV <- nDiffV <- 0
if(length(covTrack) > 0) nCovV <- length(covTrack[[1]])
if(length(diffTrack) > 0) nDiffV <- length(diffTrack[[1]])
realVL <- c()
Vcursor <- 0
if(nCovV > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL,
rep(
verticalLayout[Vcursor] / nCovV,
nCovV
)
)
}
if(nDiffV > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL,
rep(
verticalLayout[Vcursor] / nDiffV,
nDiffV
)
)
}
if(length(annoSubTrack) > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL, verticalLayout[Vcursor])
}
if(length(annoFullTrack) > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL, verticalLayout[Vcursor])
}
realVL <- round(realVL * 1000) # round to nearest integer
realHL <- horizontalLayout
hCount <- length(realHL)
vCount <- length(realVL)
bVec <- seq_len(length(realVL))
rVec <- seq_len(length(realHL))
if(vCount > 1){
tVec <- c(1, 1 + bVec[-vCount])
} else {
tVec <- 1
}
if(hCount > 1){
lVec <- c(1, 1 + rVec[-hCount])
} else {
lVec <- 1
}
# Set up patchwork
hasAnnoFT <- length(annoFullTrack) == 1
for(v in seq_len(length(realVL))) {
if(v == vCount & hasAnnoFT) {
layout_instr <- c(layout_instr, patchwork::area(
t = tVec[v], b = bVec[v], l = 1, r = hCount
))
}
for(h in seq_len(length(realHL))) {
if(v == 1 & h == 1) {
# initialize
layout_instr <- patchwork::area(
t = tVec[1], b = bVec[1], l = lVec[1], r = rVec[1]
)
} else {
layout_instr <- c(layout_instr, patchwork::area(
t = tVec[v], b = bVec[v],
l = lVec[h], r = rVec[h]
))
}
}
}
fullPlotList <- list()
plotCount <- 0
rowCount <- 0
if(length(covTrack) > 0) {
for(j in seq_len(nCovV)) {
rowCount <- rowCount + 1
for(i in seq_len(length(covTrack))) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- .pV_mod_ggplot(
covTrack[[i]][[j]], i, rowCount, hCount, length(realVL)
)
}
}
}
if(length(diffTrack) > 0) {
for(j in seq_len(nDiffV)) {
rowCount <- rowCount + 1
for(i in seq_len(length(diffTrack))) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- .pV_mod_ggplot(
diffTrack[[i]][[j]], i, rowCount, hCount, length(realVL)
)
}
}
}
if(length(annoSubTrack) > 0) {
rowCount <- rowCount + 1
for(i in seq_len(length(annoSubTrack))) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- .pV_mod_ggplot(
annoSubTrack[[i]], i, rowCount, hCount, length(realVL)
)
}
}
if(hasAnnoFT) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- annoFullTrack[[1]]
}
return(
patchwork::wrap_plots(fullPlotList) +
patchwork::plot_layout(
design = layout_instr,
widths = realHL,
heights = realVL
)
)
}
# Modify gg subplots given their position
# - removes y axis on non-left plots
# - removes legend on non-right plots
# - removes x axis on non-bottom plots
.pV_mod_ggplot <- function(p, x, y, maxX, maxY) {
theme_nonright <- theme(
legend.position = "none"
)
theme_nonleft <- theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
strip.background = element_blank(),
strip.text.x = element_blank()
)
theme_nonbottom <- theme(
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_blank()
)
if(x > 1) p <- p + theme_nonleft
if(x < maxX) p <- p + theme_nonright
if(y < maxY) p <- p + theme_nonbottom
return(p)
}
################################################################################
.gpo_highlight_anno <- function(
reduced.DT, highlight_gr
) {
introns <- reduced.DT[get("type") == "intron"]
introns[, c("highlight") := "0"]
exons <- reduced.DT[get("type") == "exon"]
exons[, c("highlight") := "0"]
misc <- reduced.DT[get("type") == "CDS"]
misc[, c("highlight") := "0"]
tr_filter <- c()
if (length(highlight_gr) == 1) {
# IR / RI only
gr <- highlight_gr[[1]]
introns.gr <- .grDT(introns)
OL <- findOverlaps(gr, introns.gr)
introns[OL@to, c("highlight") := "3"] # purple
OL2 <- findOverlaps(gr, introns.gr, type = "equal")
introns[OL2@to, c("highlight") := "2"] # red
# exons that overlap introns
exons.gr <- .grDT(exons)
OL3 <- findOverlaps(gr, exons.gr, type = "within")
exons[OL3@to, c("highlight") := "1"] # blue
} else if (length(highlight_gr) == 2) {
AS_count <- 1
for (event in highlight_gr) {
tr <- list()
for(i in seq_len(length(event))) {
OL <- findOverlaps(event[i], .grDT(introns), type = "equal")
tr[[i]] <- introns[OL@to]$transcript_id
}
tr_final <- NULL
if(length(event) == 2) {
tr_final <- intersect(tr[[1]], tr[[2]])
} else {
tr_final <- tr[[1]]
}
if(length(tr_final) > 0) {
# Highlight introns that match exact junction
gr <- event
introns.gr <- .grDT(introns)
OL <- findOverlaps(gr, introns.gr, type = "equal")
introns[
seq_len(nrow(introns)) %in% OL@to &
get("transcript_id") %in% tr_final,
c("highlight") := as.character(AS_count)
]
# Remove highlight from novel transcripts
# if not all introns highlighted
introns_novel <- introns[
grepl("novel", get("transcript_id")) &
get("highlight") == as.character(AS_count)
]
introns_novel_noHL <- introns[
get("transcript_id") %in% introns_novel$transcript_id &
get("highlight") == "0"
]
if(nrow(introns_novel_noHL) > 0) {
introns[
get("transcript_id") %in% introns_novel_noHL$transcript_id,
c("highlight") := "0"
]
tr_final <- setdiff(tr_final,
introns_novel_noHL$transcript_id)
}
coord_keys_start <- end(gr[1]) + 1
coord_keys_end <- start(gr[1]) - 1
if (length(gr) == 2) {
coord_keys_start <- c(coord_keys_start, end(gr[2]) + 1)
coord_keys_end <- c(coord_keys_end, start(gr[2]) - 1)
}
exons[get("transcript_id") %in% tr_final &
(get("start") %in% coord_keys_start |
get("end") %in% coord_keys_end),
c("highlight") := as.character(AS_count)]
}
AS_count <- AS_count + 1
}
}
for(hl in c("1", "2")) {
exons_selected <- exons[get("highlight") == hl]
hl_trid <- unique(exons_selected$transcript_id)
OL <- findOverlaps(
.grDT(misc),
.grDT(exons_selected)
)
misc[seq_len(nrow(misc)) %in% unique(OL@from) &
get("transcript_id") %in% hl_trid,
c("highlight") := hl]
}
return(rbind(introns, exons, misc))
}
.gpo_filter_anno <- function(
reduced.DT, transcripts.DT,
selected_transcripts = "",
plot_key_isoforms = FALSE
) {
# Remove by named transcripts only
if(is_valid(selected_transcripts)) {
if (
all(selected_transcripts %in%
c(
transcripts.DT$transcript_id,
transcripts.DT$transcript_name
)
)
) {
transcripts.DT <- transcripts.DT[
get("transcript_name") %in% selected_transcripts |
get("transcript_id") %in% selected_transcripts
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
return(reduced.DT)
}
}
if(plot_key_isoforms) {
reduced.DT.HL <- reduced.DT[get("highlight") != "0"]
reduced.DT.nonHL <- reduced.DT[get("highlight") == "0"]
reduced.DT.nonHL <- reduced.DT.nonHL[
grepl("novel", get("transcript_id"))
]
if(nrow(reduced.DT.HL) != 0) {
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% reduced.DT.HL$transcript_id &
# remove non-highlighted novel elements
!(get("transcript_id") %in% reduced.DT.nonHL$transcript_id)
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
return(reduced.DT)
}
}
# Otherwise - do nothing
return(reduced.DT)
}
# Allows re-stack. Useful for filtered DTlist
.gcd_stack_anno <- function(
DTlist, view_start, view_end,
reverseGenomeCoords = FALSE,
condensed = FALSE
) {
transcripts.DT <- copy(DTlist$transcripts.DT)
reduced.DT <- copy(DTlist$reduced.DT)
condense_this <- condensed
# Filter before condense
reduced.DT <- reduced.DT[get("end") > view_start & get("start") < view_end]
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% reduced.DT$transcript_id]
if (condense_this != TRUE & nrow(transcripts.DT) <= 100) {
condense_this <- FALSE
transcripts.DT[, c("group_id") := get("transcript_id")]
reduced.DT[, c("group_id") := get("transcript_id")]
} else {
condense_this <- TRUE
transcripts.DT[, c("group_id") := get("gene_id")]
reduced.DT[transcripts.DT, on = "transcript_id",
c("group_id") := get("gene_id")]
}
group.grl <- split(.grDT(transcripts.DT), transcripts.DT$group_id)
group.DT <- as.data.table(range(group.grl))
group.DT$group <- NULL
data.table::setnames(group.DT, "group_name", "group_id")
# Filter before stack
# group.DT <- group.DT[get("end") > view_start & get("start") < view_end]
# apply plot_order on transcripts.DT
OL <- findOverlaps(.grDT(group.DT), .grDT(group.DT), ignore.strand = TRUE)
group.DT$plot_level <- 1
cur_level <- 1
while (any(group.DT$plot_level == cur_level)) {
j <- match(cur_level, group.DT$plot_level)
repeat {
bump_up_trs <- unique(OL@to[OL@from == j])
bump_up_trs <- bump_up_trs[bump_up_trs > j]
bump_up_trs <- bump_up_trs[
group.DT$plot_level[bump_up_trs] == cur_level]
if (length(bump_up_trs) > 0)
group.DT[bump_up_trs, c("plot_level") := cur_level + 1]
j <- j + match(cur_level, group.DT$plot_level[-seq_len(j)])
if (is.na(j)) break
}
cur_level <- cur_level + 1
}
if (condense_this == TRUE) {
group.DT[transcripts.DT, on = "group_id",
c("group_name", "group_biotype") :=
list(get("i.gene_name"), get("i.gene_biotype"))]
} else {
group.DT[transcripts.DT, on = "group_id",
c("group_name", "group_biotype") :=
list(get("i.transcript_name"), get("i.transcript_biotype"))]
}
reduced.DT$group_id <- factor(reduced.DT$group_id,
unique(group.DT$group_id), ordered = TRUE)
reduced.DT[group.DT, on = "group_id",
c("plot_level") := get("i.plot_level")]
setorderv(reduced.DT, "highlight")
return(list(
group.DT = group.DT,
reduced.DT = reduced.DT,
condense_this = condense_this
))
}
alexchwong/SpliceWiz documentation built on March 17, 2024, 3:16 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .gcd_stack_anno .gpo_filter_anno .gpo_highlight_anno .pV_mod_ggplot .pV_assemble_ggplot .pV_assemble_covPlotly2 .plotly_makeRectData .plotly_makeLineData .plotly_makeJuncCurveData .pV_trackList_from_list .pV_trackList_from_vector .gPO_getJunc .gPO_PSIjuncStats .pV_diffStat_ttest .gPO_normJuncStats .gPO_covStats .gPO_getCoverage .gPO_check_strand .gPO_check_tracks .pV_sashimiArc .pV_seed_plotly_HLexons .pV_seed_plotly_annoTrack .pV_seed_instant_plotly_HLexons .pV_seed_instant_plotly_annoTrack .pV_seed_ggplot_HLexons .pV_seed_ggplot_annoTrack .pV_highlight_to_colors .pV_plotExons .pV_plotRef .pV_seed_plotly_diffTrack .pV_seed_ggplot_diffTrack .pV_seed_plotly_covTrack .pV_seed_ggplot_covTrack .pV_plotDiffTracks .pV_plotCoverageTracks .pV_getJuncCoords .pV_filterTranscripts .pV_getEventCoords .pV_getAllowedCoords plotView plotAnnoTrack getPlotObject covPlotObject
Documented in getPlotObject plotAnnoTrack plotView
#' Versatile coverage plots for SpliceWiz
#'
#' Here, we implement fast and versatile ggplot and plotly based
#' coverage and sashimi plots. Users can plot with unlimited number
#' of individual, individual-normalized, or group-normalized tracks.
#' Also implemented is user-defined group-comparison differential plots
#' (including t-test plots). Additionally, users can generate ggplots
#' subsetted by exon groups. See details below.
#'
#' @details
#' The typical pipeline for plotting versatile coverage plots is as follows:
#'
#' * A `covDataObject` is generated by calling `getCoverageData()` using an
#' input `NxtSE` object. This step retrieves coverage, junction counts and
#' normalization data for the relevant genomic region being queried. A new
#' `covDataObject` is necessary when querying a new genomic region.
#' * A `covPlotObject` is generated by calling `getPlotObject()` using an
#' input `covDataObject`. This step retrieves alternative splicing event
#' specific data, such as normalized coverages, or group combined coverages.
#' A new `covPlotObject` is required when changing `condition`, `Event`,
#' `strand`, or when querying using a different set of `tracks`.
#' * Plots can be generated by calling `plotView()` using a `covPlotObject`.
#' Interactive plotly plots can be generated by setting `usePlotly = TRUE`,
#' otherwise, static plots are generated. For interactive plots, a `covPlotly`
#' object is returned, which contains raw data which is downsampled by
#' pixel resolution prior to plotting for performance reasons. A new
#' `covPlotly` is required unless one only wishes to downsample the resolution
#' - see [setResolution] for `covPlotly` objects.
#'
#' Tracks are now versatile (unlimited). Samples are retrieved by individual
#' sample names at `getCoverageData()`. If `condition` is set in
#' `getPlotObject()`, track names are defined by their condition categorical
#' names; otherwise, tracks are named by individual samples when retrieved using
#' `getPlotObject()`.
#'
#' * When calling `plotView()`, `trackList` by default displays all tracks as
#' ordered in the `covPlotObject`. Users can supply a vector containing either
#' the track names (or numbers, as ordered in the `covPlotObject`).
#' Alternatively, multiple traces can be stacked in a single track by
#' using a list, e.g. `trackList = list(A = c(1,2,3), B = c(4,5,6))`.
#' * For differential comparisons, `diffList` takes a list of pairs of samples.
#' For example, if `trackList = list("A", "B")`, then setting
#' `diffList = list(c("A", "B"))` will compare groups "A" and "B". This is only
#' activated by setting `diff_stat` to anything other than `none`. For now,
#' only `t-test` is supported.
#'
#' `plotView()` supports plotting by exon ranges, for which only static plots
#' are currently supported. The workflow for generating such a plot is as
#' follows:
#'
#' * A GRanges object is returned by the `plotView()` function and setting
#' `showExonRanges = TRUE`. `plotView()` will simultaneously show an
#' annotation plot of exons labelled by their "exon names", which is the
#' transcript name appended with "-E" followed by the exon number.
#' * If `plotView()` is called and `usePlotly = TRUE` is set, a `covPlotly`
#' object is returned. Calling `showExons()` on this object will display
#' a plotly plot showing exon names, and returning a GRanges object of exon
#' ranges.
#' * Exon ranges can be supplied to the `plotView()` function by setting
#' the `plotRanges` parameter as a GRanges object. This will generate a static
#' plot showing coverage plots segmented by exons.
#'
#' @param object
#' For `getPlotObject()`, a `covDataObject` created using
#' `getCoverageData`. For `plotView()`, a `covPlotObject` created using
#' `getPlotObject()`.
#' @param Event The EventName of the alternative splicing event
#' which will be highlighted and used for normalization
#' @param strand The strand for coverage / junction plotting. Options are `"+"`,
#' `"-"`, or `"*"` (unstranded - default)
#' @param tracks Sample names or condition categories
#' @param condition For condition-based group plots, the name of the condition.
#' @param view_start,view_end The start and end coordinates for plotting
#' @param oldP (Optional) If plotting the same tracks and track-widths,
#' supplying the old `covPlotly` object (returned from a previous call to
#' `plotView()`) results in faster run-time (as plotly::subplot is a time-
#' consuming function)
#' @param centerByEvent (default `FALSE`) If true, centers the view to the
#' specified `Event`
#' @param EventZoomFactor If `centerByEvent = TRUE`, the zoom-out factor to plot
#' the view. Zooms out in exponents of 3
#' (i.e., zoom of 1 means 3x, 2 means 9x, and 0 means 1x)
#' @param EventBasesFlanking (default `100`)
#' If `centerByEvent = TRUE`, includes how many bases flanking the event.
#' @param resolution The number of horizontal "pixels" or data-points to plot.
#' This is calculated per sub-plot. Smaller numbers lead to lower resolution
#' but faster plots. Default is `5000`
#' @param trackList A list, with each element being a vector of 1 or more track
#' names or indices to plot. If a vector is supplied it will be coerced to a
#' list
#' @param diffList A list, with each element being a vector of size 2,
#' containing names or indices of which tracks to contrast.
#' @param diff_stat (default "t-test") Which statistical method to perform
#' differential comparisons.
#' @param reverseGenomeCoords If `TRUE`, the genomic coordinate axis will be
#' reversed to plot negative stranded genes
#' @param ribbon_mode The statistic to represent variance. Options are "sd" -
#' standard deviation, "sem" - standard error of the mean, "ci" -
#' 95% confidence interval, or "none"
#' @param normalizeCoverage If `TRUE`, coverages and junctions of individual
#' samples will be normalized by the given `Event`.
#' @param plotAnnotations Whether the main annotation track should be plotted
#' @param plotAnnoSubTrack If plotting by exon ranges (using `plotRanges`),
#' whether a separate sub-track showing zoomed-in exons should be shown
#' above the main annotation track (and below the coverage plots)
#' @param showExonRanges (only applies if `usePlotly = FALSE`) Whether the
#' main annotation track should be replaced by labeled exon names. If `TRUE`
#' the returned value of `plotView()` is a named GRanges object containing
#' the exon ranges
#' @param verticalLayout A vector (of length 4) containing relative heights of
#' the following elements: (1) main block of coverage tracks, (2) differential
#' track, (3) annotation sub-track, and (4) main annotation track. Default
#' `c(4,1,1,2)`
#' @param horizontalLayout A vector containing relative widths of coverage
#' tracks. Only used alongside `plotRanges` with more than 1 range to plot.
#' If omitted, `plotView` will attempt to scale widths to the widths of the
#' exon ranges.
#' @param filterByTranscripts (default `""`) One or more named transcripts to
#' filter the annotation track.
#' @param filterByEventTranscripts (default `FALSE`) If `TRUE`, only transcripts
#' involved in the given `Event` will be plotted, if any
#' @param filterByExpressedTranscripts (default `TRUE`) Only transcripts with
#' supported junctions will be plotted on the annotation axis. An expressed
#' junction is that which contains more than the minimum `junctionThreshold`
#' in at least 1 track
#' @param condenseTranscripts Whether to plot by genes `TRUE` or transcripts
#' `FALSE`
#' @param plotJunctions Whether to plot junction counts as numbered arcs. Plots
#' normalized junctions if `normalizeCoverage = TRUE`.
#' @param plotJuncPSI If plotting group coverage plots, whether to plot
#' mean +/- sd of normalized junction counts `FALSE`, or estimated junction
#' PSI based on SpliceOver metric applied to each junction `TRUE`.
#' @param junctionThreshold (default `0.01`) Junctions with expressions below
#' this threshold will not be plotted. For raw counts, this is a fraction of
#' maximum coverage value of the track.
#' @param plotRanges A GRanges object containing one or more exon ranges to
#' plot. If given, `view_start` and `view_end` will be ignored. Typical use
#' is to use the output of the `plotView(..., usePlotly = FALSE)`, which
#' returns a named GRanges object, then subset this output by exon name.
#' @param rangesBasesFlanking (default `100`) How many flanking bases to add to
#' each of `plotRanges`. Ignored if only 1 range given (or using
#' `view_start` and `view_end`)
#' @param usePlotly If `TRUE`, returns a `covPlotly` object containing the
#' plotly-based interactive plot. If `FALSE`, returns a ggplot object.
#' @param ... Ignored / not used
#' @return
#' For `getPlotObject()`: A `covPlotObject` object containing Event-based
#' data to create coverage plots using `plotView()`.
#'
#' For `plotView()`:
#'
#' * If `usePlotly = TRUE`, returns a `covPlotly` object containing plotly-based
#' interactive plot
#' * If `usePlotly = FALSE`, returns a patchwork-assembled static plot, unless
#' `showExonRanges = TRUE` in which it shows the plot and returns a named
#' GRanges object containing exon ranges.
#' @examples
#' se <- SpliceWiz_example_NxtSE(novelSplicing = TRUE)
#'
#' # Assign annotation of the experimental conditions
#' colData(se)$treatment <- rep(c("A", "B"), each = 3)
#'
#' # Retrieve coverage data for all samples for the gene "SRSF3" (and surrounds)
#'
#' dataObj <- getCoverageData(
#' se,
#' Gene = "SRSF3",
#' tracks = colnames(se)
#' )
#'
#' # Retrieves raw / normalized coverage / junction data for the
#' # specified SRSF3 skipped exon event:
#'
#' plotObj_samples <- getPlotObject(
#' dataObj,
#' Event = "SE:SRSF3-203-exon4;SRSF3-202-int3"
#' )
#'
#' # Retrieves data for samples grouped by the specified condition
#'
#' plotObj_group <- getPlotObject(
#' dataObj,
#' Event = "SE:SRSF3-203-exon4;SRSF3-202-int3",
#' condition = "treatment",
#' tracks = c("A", "B")
#' )
#'
#' # Display tracks and conditions of covPlotObject
#'
#' tracks(plotObj_group)
#' condition(plotObj_group)
#'
#' # Show static ggplots
#'
#' plotView(plotObj_samples)
#'
#' plotView(plotObj_group, centerByEvent = TRUE)
#'
#' # Plot junctions using PSI estimates based on individual junction SpliceOver
#' # metrics
#'
#' plotView(plotObj_group, centerByEvent = TRUE, plotJuncPSI = TRUE)
#'
#'
#' # Show normalized coverages, individual samples stacked in grouped tracks
#'
#' plotView(
#' plotObj_samples,
#' normalizeCoverage = TRUE,
#' trackList = list(A = c(1,2,3), B = c(4,5,6))
#' )
#'
#' # Show stacked group comparisons with t-test
#'
#' plotView(
#' plotObj_group,
#' trackList = list(c(1,2)),
#' diffList = list(c("A", "B")),
#' diff_stat = "t-test"
#' )
#'
#' # Show interactive plotly:
#'
#' if(interactive()) {
#' p <- plotView(plotObj_samples, usePlotly = TRUE)
#' show(p)
#' }
#'
#' # Show exons with coverage plot
#'
#' # static:
#' gr <- plotView(plotObj_samples, showExonRanges = TRUE)
#'
#' # interactive:
#' if(interactive()) {
#' p <- plotView(plotObj_samples, usePlotly = TRUE)
#' gr <- showExons(p)
#' }
#'
#' # Plot coverage by exons
#'
#' p <- plotView(plotObj_samples,
#' plotRanges = gr[c("SRSF3-203-E3", "SRSF3-203-E4", "SRSF3-203-E5")],
#' horizontalLayout = c(1,1,1)
#' )
#'
#' @name covPlotObject-class
#' @aliases
#' tracks tracks,covPlotObject-method
#' condition condition,covPlotObject-method
#' @seealso [getCoverageData] [covPlotly-class]
#' @md
NULL
covPlotObject <- function(
args = list(),
cov = list(),
norm_cov = list(),
junc = list(),
norm_junc = list(),
junc_PSI = list(),
junc_stats = list(),
cov_stats = list(),
annotation = list()
) {
obj <- new("covPlotObject",
args = args,
cov = cov,
norm_cov = norm_cov,
junc = junc,
norm_junc = norm_junc,
junc_PSI = junc_PSI,
junc_stats = junc_stats,
cov_stats = cov_stats,
annotation = annotation
)
obj
}
#' @describeIn covPlotObject-class Generates a covPlotObject object from a
#' covDataObject. Allows users to change parameters such as viewing window,
#' conditions, tracks, and other parameters, for customizing plot parameters
#' @export
getPlotObject <- function(
object,
Event, # To specify normalization event
strand = c("*", "+", "-"),
tracks, # can be a list of iterables
condition,
...
) {
args <- object@args
if(!missing(tracks)) args[["tracks"]] <- tracks
if(!missing(condition)) args[["condition"]] <- condition
if(!missing(Event)) args[["Event"]] <- Event
if(!("Event" %in% names(args)) && "condition" %in% names(args)) {
.log(paste("In getPlotObject,",
"for plotting by condition, normalizing `Event` must be supplied"
))
}
# Check tracks and conditions in args are legit, modify as necessary
args <- .gPO_check_tracks(object, args)
# Junction strand is always same as requested strand
# But need to check strandedness of RNA-seq data before deciding coverage
if(
!missing(strand) &&
length(strand) == 1 &&
strand %in% c("*", "+", "-")
) {
args[["strand"]] <- strand
} else if(!("strand" %in% names(args))) {
args[["strand"]] <- "*" # default strand
}
# Check strand against stranded-ness of RNA-seq
args <- .gPO_check_strand(object, args)
# Always get these
cov <- .gPO_getCoverage(object, args)
junc <- .gPO_getJunc(object, args, cov)
norm_cov <- norm_junc <- list()
if("Event" %in% names(args)) {
norm_cov <- .gPO_getCoverage(object, args, normalize = TRUE)
# if single sample
norm_junc <- .gPO_getJunc(object, args, norm_cov, normalized = TRUE)
}
# Plot stats if applicable
cov_stats <- junc_PSI <- junc_stats <- list()
if("condition" %in% names(args)) {
# compile junction data
cov_stats <- .gPO_covStats(norm_cov)
junc_stats <- .gPO_normJuncStats(object, args, cov_stats)
junc_PSI <- .gPO_PSIjuncStats(object, args, norm_cov)
}
# plot relevant annotations
highlight_gr <- list()
if(
"Event" %in% names(args) &&
args[["Event"]] %in% rownames(object@normData$rowData)
) {
row <- object@normData$rowData[args[["Event"]],]
args[["EventRegion"]] <- row$EventRegion
if (row$EventType %in% c("MXE", "SE")) {
highlight_gr[[1]] <- coord2GR(c(row$Event1a, row$Event2a))
} else {
highlight_gr[[1]] <- coord2GR(row$Event1a)
}
if (row$EventType %in% c("MXE")) {
highlight_gr[[2]] <- coord2GR(c(row$Event1b, row$Event2b))
} else if (row$EventType %in% c("SE", "A3SS", "A5SS", "ALE", "AFE")) {
highlight_gr[[2]] <- coord2GR(row$Event1b)
}
}
args[["highlight_gr"]] <- highlight_gr
DTlist <- object@annotation
DTlist$reduced.DT <- .gpo_highlight_anno(
DTlist$reduced.DT, args[["highlight_gr"]]
)
return(covPlotObject(
args = args,
cov = cov,
norm_cov = norm_cov,
junc = junc,
norm_junc = norm_junc,
junc_PSI = junc_PSI,
junc_stats = junc_stats,
cov_stats = cov_stats,
annotation = DTlist
))
}
#' @describeIn covPlotObject-class Returns the tracks contained in the
#' covPlotObject object
#' @export
setMethod("tracks", c(object = "covPlotObject"), function(
object
) {
return(object@args[["tracks"]])
})
#' @describeIn covPlotObject-class Returns the condition value set in the
#' covPlotObject object
#' @export
setMethod("condition", c(object = "covPlotObject"), function(
object
) {
return(object@args[["condition"]])
})
#' @describeIn covDataObject-class Directly plots the annotation from a
#' covDataObject.
#' @export
plotAnnoTrack <- function(
object, Event,
view_start, view_end,
reverseGenomeCoords = FALSE,
condensed = FALSE,
selected_transcripts = "",
plot_key_isoforms = FALSE,
usePlotly = FALSE,
...
) {
if(missing(view_start)) view_start <- object@args$view_start
if(missing(view_end)) view_end <- object@args$view_end
if(view_start > view_end) {
reverseGenomeCoords <- !reverseGenomeCoords
view_start <- view_start + view_end
view_end <- view_start - view_end
view_start <- view_start - view_end
}
highlight_gr <- list()
if(!missing(Event) && !("rowData" %in% names(object@normData))) {
.log(paste(
"In plotAnnoTrack(),",
"Plotting ASE from a SpliceWiz reference is currently not supported"
))
}
if(!missing(Event) && Event %in% rownames(object@normData$rowData)) {
row <- object@normData$rowData[Event,]
if (row$EventType %in% c("MXE", "SE")) {
highlight_gr[[1]] <- coord2GR(c(row$Event1a, row$Event2a))
} else {
highlight_gr[[1]] <- coord2GR(row$Event1a)
}
if (row$EventType %in% c("MXE")) {
highlight_gr[[2]] <- coord2GR(c(row$Event1b, row$Event2b))
} else if (row$EventType %in% c("SE", "A3SS", "A5SS", "ALE", "AFE")) {
highlight_gr[[2]] <- coord2GR(row$Event1b)
}
}
DTlist <- object@annotation
DTlist$reduced.DT <- .gpo_highlight_anno(
DTlist$reduced.DT, highlight_gr
)
DTlist$reduced.DT <- .gpo_filter_anno(
DTlist$reduced.DT, DTlist$transcripts.DT,
selected_transcripts = selected_transcripts,
plot_key_isoforms = plot_key_isoforms
)
DTplotlist <- .gcd_stack_anno(
DTlist, view_start, view_end, reverseGenomeCoords
)
view_gr <- GRanges(
object@args[["view_chr"]],
IRanges(object@args[["view_start"]], object@args[["view_end"]])
)
p <- .pV_plotRef(
DTplotlist, view_gr,
reverseGenomeCoords = reverseGenomeCoords,
add_information = TRUE,
usePlotly = usePlotly,
instantPlotly = usePlotly
)
return(p)
}
#' @describeIn covPlotObject-class Creates a coverage plot using the stored
#' data in the covPlotObject
#' @export
plotView <- function(
object,
# for single window view (to simplify things)
view_start,
view_end,
# debug = FALSE,
oldP = covPlotly(),
# Event-centric plotting
centerByEvent = FALSE,
EventZoomFactor = 0.2,
EventBasesFlanking = 100,
resolution = 5000,
# specify tracks
trackList = list(),
# specify differential comparisons
diff_stat = c("t-test", "none"),
diffList = list(),
reverseGenomeCoords = FALSE,
ribbon_mode = c("sd", "sem", "ci", "none"),
# plot raw or normalized coverage (only works for sample tracks)
normalizeCoverage = FALSE,
# specify annotation track
plotAnnotations = TRUE,
plotAnnoSubTrack = TRUE,
showExonRanges = FALSE,
# Plot layout
verticalLayout = c(4,1,1,2),
horizontalLayout = c(),
# filter annotations
filterByTranscripts = "",
# plot only by transcripts with highlights
filterByEventTranscripts = FALSE,
# filter annotations
filterByExpressedTranscripts = TRUE,
# whether to condense transcripts by gene
condenseTranscripts = FALSE,
# whether to plot sashimi arcs
plotJunctions = TRUE,
plotJuncPSI = FALSE,
junctionThreshold = 0.01,
# whether to plot by exon windows
plotRanges = GRanges(),
rangesBasesFlanking = 100,
# ggplot or plotly object?
usePlotly = FALSE,
...
) {
use_DT <- TRUE
if(!is(object, "covPlotObject")) .log(paste(
"In plotView,", "object must be a covPlotObject"
))
args <- object@args
# inject variables into args
if(!missing(view_start)) args[["view_start"]] <- view_start
if(!missing(view_end)) args[["view_end"]] <- view_end
args[["reverseGenomeCoords"]] <- reverseGenomeCoords
ribbon_mode <- match.arg(ribbon_mode)
if(!is_valid(ribbon_mode)) ribbon_mode <- "none"
args[["ribbon_mode"]] <- ribbon_mode
if(
(missing(view_start) | missing(view_end)) &&
centerByEvent &
"EventRegion" %in% names(args)
) {
args <- .pV_getEventCoords(
object, args, EventZoomFactor,
EventBasesFlanking
)
}
diff_stat <- match.arg(diff_stat)
if(!is_valid(diff_stat)) diff_stat <- "none"
# plotRanges checking
if(!is(plotRanges, "GRanges")) .log(paste(
"In plot() for class covPlotObject,",
"`plotRanges` must be a GRanges object"
))
# Plot single window if no plotRanges given
if(length(plotRanges) == 0) plotRanges <- GRanges(
args[["view_chr"]],IRanges(args[["view_start"]], args[["view_end"]])
) # strand is ignored
if(
any(as.character(seqnames(plotRanges)) != args[["view_chr"]])
) {
seqs <- unique(as.character(seqnames(plotRanges)))
.log(paste(
"In plotView() for class covPlotObject,",
"Some elements in `plotRanges` have seqnames that do not match that",
"of covPlotObject", seqs[!(seqs %in% args[["view_chr"]])]
))
}
# Make all plotRanges strand-agnostic
strand(plotRanges) <- "*"
# disable interactive if multi plot
if(usePlotly == TRUE & length(plotRanges) > 1) {
.log(paste(
"In plotView,",
"interactive plots are not supported for multi-view plotting"
), "warning")
usePlotly <- FALSE
}
# Expand plotRanges using rangesBasesFlanking
# - only used if more than 1 range
rbf <- rangesBasesFlanking
if(length(plotRanges) > 1) {
for(i in seq_len(length(plotRanges))) {
if(start(plotRanges[i]) - rbf > args[["limit_start"]]) {
start(plotRanges[i]) <- start(plotRanges[i]) - rbf
} else {
start(plotRanges[i]) <- args[["limit_start"]]
}
if(end(plotRanges[i]) + rbf < args[["limit_end"]]) {
end(plotRanges[i]) <- end(plotRanges[i]) + rbf
} else {
end(plotRanges[i]) <- args[["limit_end"]]
}
}
}
# What is the full range
fullRange <- range(plotRanges)
if(!missing(view_start) & !missing(view_end)) {
# whatever's larger
fullRange <- range(c(fullRange,
GRanges(args[["view_chr"]], IRanges(view_start, view_end))
))
}
if(
start(fullRange) < args[["limit_start"]] |
end(fullRange) > args[["limit_end"]]
) {
.log(paste(
"Given range is outside that supported by parent covPlotObject.",
"Suggest regenerating covPlotObject by calling getPlotObject()"
))
}
fetchRange <- fullRange
if(usePlotly) {
start(fetchRange) <- max(args[["limit_start"]],
start(fullRange) - width(fullRange))
end(fetchRange) <- min(args[["limit_end"]],
end(fullRange) + width(fullRange))
}
# Sort plotRanges
plotRanges <- sort(plotRanges, decreasing = args[["reverseGenomeCoords"]])
# Track list checking - should be a list of indices
# If no trackList given, plot all tracks individually
if(length(trackList) == 0) trackList <- lapply(
seq_len(length(args[["tracks"]])),
function(i) i
)
if(!is(trackList, "list")) {
args[["trackList"]] <- .pV_trackList_from_vector(trackList, args)
} else {
args[["trackList"]] <- .pV_trackList_from_list(trackList, args)
}
# Structure of plot
# | 1a || 1b || 1c |
# | 2a || 2b || 2c |
# | 3a || 3b || 3c |
# | pa || pb || pc | - diff track
# | A || B || C | - annotation subtrack
# | annotation | - annotation full track
covTrack <- list()
diffTrack <- list()
annoSubTrack <- list()
annoFullTrack <- list() # list of 1
exonFullTrack <- list() # list of 1
# Proportion all viewing frames
widthSum <- sum(width(plotRanges))
widthFrac <- width(plotRanges) / widthSum
if(use_DT && !usePlotly) use_DT <- FALSE
# For filtering transcripts by expression later
aggJuncCoords <- c()
reservedCoords <- c()
for(j in seq_len(length(plotRanges))) {
range_gr <- plotRanges[j]
# Work out which coords to plot
juncCoords <- .pV_getJuncCoords(
object, args, normalizeCoverage,
range_gr = range_gr,
junctionThreshold
)
aggJuncCoords <- c(aggJuncCoords, juncCoords)
}
reservedCoords <- sort(unique(c(
aggJuncCoords[["jc_start"]] - 1,
aggJuncCoords[["jc_start"]],
aggJuncCoords[["jc_end"]],
aggJuncCoords[["jc_end"]] + 1
)))
if(usePlotly) {
covTrack[[j]] <- .pV_plotCoverageTracks(
object, args, normalizeCoverage,
range_gr = fullRange,
plotJunctions = plotJunctions, plotJuncPSI = plotJuncPSI,
junctionThreshold = junctionThreshold,
usePlotly = TRUE,
use_DT = use_DT
)
if(diff_stat != "none" && length(diffList) > 0) {
diffTrack[[j]] <- .pV_plotDiffTracks(
object, args, diffList,
diff_stat, fullRange,
usePlotly = usePlotly
)
}
} else {
# Generate full plotly first, then segment later
tmpCov <- .pV_plotCoverageTracks(
object, args, normalizeCoverage,
range_gr = fullRange,
plotJunctions = plotJunctions, plotJuncPSI = plotJuncPSI,
junctionThreshold = junctionThreshold,
usePlotly = FALSE,
use_DT = use_DT
)
tmpCovPlot <- tmpCov$ggplot
tmpCovJuncAnno <- as.data.table(tmpCov$dfJn_all)
tmpDiffPlot <- NULL
if(diff_stat != "none" && length(diffList) > 0) {
tmpDiffPlot <- .pV_plotDiffTracks(
object, args, diffList,
diff_stat, fullRange,
usePlotly = FALSE
)
}
for(j in seq_len(length(plotRanges))) {
range_gr <- plotRanges[j]
okCoords <- .pV_getAllowedCoords(
start(range_gr), end(range_gr),
reservedCoords, resolution
)
if(args[["reverseGenomeCoords"]]) {
plotViewStart <- end(range_gr)
plotViewEnd <- start(range_gr)
} else {
plotViewStart <- start(range_gr)
plotViewEnd <- end(range_gr)
}
covTrack[[j]] <- list()
dfJnSum_list <- list()
if(plotJunctions) {
dtJn_all <- as.data.table(tmpCovJuncAnno)
allGroups <- unique(dtJn_all$group)
for(elem in allGroups) {
dtJn <- copy(dtJn_all[get("group") == elem])
dtJn <- dtJn[
get("x") >= start(range_gr) &
get("x") <= end(range_gr)
]
dtJn[, c("xlabel", "ylabel") := list(
mean(get("x")), mean(get("ytmp"))),
by = "info"
]
dtJn <- unique(
dtJn[,
c("info", "value", "xlabel",
"ylabel", "group", "covTrack"),
with = FALSE
]
)
dfJnSum_list[[elem]] <- dtJn
}
}
dfJnSum_all <- rbindlist(dfJnSum_list)
for(i in seq_len(length(tmpCovPlot))) {
curTrack <- names(args[["trackList"]])[i]
covTrack[[j]][[i]] <- tmpCovPlot[[i]]
ymax <- max(layer_scales(tmpCovPlot[[i]])$y$range$range)
if(nrow(dfJnSum_all) > 0) {
tmpCovSubset <- dfJnSum_all[get("covTrack") == curTrack]
# disable junction plots for multiple samples per track
if(
nrow(tmpCovSubset) > 0 &&
length(unique(tmpCovSubset$group)) == 1
) {
covTrack[[j]][[i]] <- covTrack[[j]][[i]] +
geom_text(
data = tmpCovSubset,
aes(
x = get("xlabel"),
y = get("ylabel"),
label = get("value")
)
)
}
}
if(!is.null(okCoords)) {
covTrack[[j]][[i]]$data <- covTrack[[j]][[i]]$data[
covTrack[[j]][[i]]$data$x %in% okCoords,]
}
covTrack[[j]][[i]] <- covTrack[[j]][[i]] +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(0, 1.2 * ymax),
expand = FALSE
)
}
if(!is.null(tmpDiffPlot)) {
diffTrack[[j]] <- list()
for(i in seq_len(length(tmpDiffPlot))) {
diffTrack[[j]][[i]] <- tmpDiffPlot[[i]]
ymax <- max(layer_scales(tmpDiffPlot[[i]])$y$range$range)
if(!is.null(okCoords)) {
diffTrack[[j]][[i]]$data <- diffTrack[[j]][[i]]$data[
diffTrack[[j]][[i]]$data$x %in% okCoords,]
}
diffTrack[[j]][[i]] <- diffTrack[[j]][[i]] +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(0, 1.2 * ymax),
expand = FALSE
)
}
}
}
}
# bench <- system.time({
DTlist <- .pV_filterTranscripts(
object, args,
filterByTranscripts = filterByTranscripts,
filterByEventTranscripts = filterByEventTranscripts,
filterByExpressedTranscripts = filterByExpressedTranscripts,
aggJuncCoords = aggJuncCoords
)
DTplotlist <- .gcd_stack_anno(
DTlist,
start(fetchRange), end(fetchRange),
reverseGenomeCoords = args[["reverseGenomeCoords"]],
condensed = condenseTranscripts
)
DTplotlist[["exonRanges"]] <- plotRanges
exons_gr <- GRanges()
if(length(plotRanges) > 1) {
if(plotAnnoSubTrack) {
for(j in seq_len(length(plotRanges))) {
range_gr <- plotRanges[j]
annoSubTrack[[j]] <- .pV_plotRef(
DTplotlist, range_gr,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = FALSE,
usePlotly = usePlotly
)
}
if(plotAnnotations) {
annoFullTrack[[1]] <- .pV_plotRef(
DTplotlist, fullRange,
showPlotRanges = plotRanges,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
}
} else if(plotAnnotations) {
annoFullTrack[[1]] <- .pV_plotRef(
DTplotlist, fullRange,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
}
} else {
# Annotate and highlight exon regions
transcripts <- copy(DTlist$transcripts.DT)
reduced <- copy(DTplotlist$reduced.DT)
reduced <- reduced[!is.na(reduced$plot_level)]
exonsOnly <- reduced[get("type") == "exon"]
exonsOnly$transcript_name <- transcripts$transcript_name[match(
exonsOnly$transcript_id, transcripts$transcript_id)]
exonsOnly[, c("exon_name") := paste0(
get("transcript_name"), "-E", get("aux_id")
)]
exonsOnly[, c("feature_id") := NULL]
exons_gr <- .grDT(exonsOnly, keep.extra.columns = TRUE)
names(exons_gr) <- exons_gr$exon_name
# Remove exons that lie outside view range
OL <- findOverlaps(exons_gr, fetchRange, type = "within")
exons_gr <- exons_gr[sort(unique(OL@from))]
if(!plotAnnotations & showExonRanges) showExonRanges <- FALSE
if(plotAnnotations) {
exonFullTrack[[1]] <- .pV_plotExons(
DTplotlist, fullRange, exons_gr,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
annoFullTrack[[1]] <- .pV_plotRef(
DTplotlist, fullRange,
reverseGenomeCoords = args[["reverseGenomeCoords"]],
add_information = TRUE,
usePlotly = usePlotly
)
}
}
# })
# message("Annotation data generation time:")
# print(bench)
keepVLayout <- c(
length(covTrack) > 0,
length(diffTrack) > 0,
length(annoSubTrack) > 0,
length(annoFullTrack) > 0
)
originalVerticalLayout <- verticalLayout
if(length(originalVerticalLayout) == 4) {
originalVerticalLayout <- originalVerticalLayout[-3]
}
if(length(verticalLayout) == 4) {
verticalLayout <- verticalLayout[keepVLayout]
} else if(length(verticalLayout) == sum(keepVLayout)) {
# do nothing
} else {
# revert to default
verticalLayout <- c(4,1,1,2)
verticalLayout <- verticalLayout[keepVLayout]
}
if(length(horizontalLayout) == 0) {
horizontalLayout <- ceiling(10 * widthFrac / sum(widthFrac))
} else if(length(horizontalLayout) == length(covTrack)) {
horizontalLayout <- ceiling(horizontalLayout)
} else {
horizontalLayout <- ceiling(10 * widthFrac / sum(widthFrac))
}
if(!args[["reverseGenomeCoords"]]) {
plotViewStart <- start(fullRange)
plotViewEnd <- end(fullRange)
} else {
plotViewStart <- end(fullRange)
plotViewEnd <- start(fullRange)
}
# if(debug) {
# return(list(
# plotViewStart = plotViewStart, plotViewEnd = plotViewEnd,
# covTrack = covTrack,
# diffTrack = diffTrack,
# annoSubTrack = annoSubTrack,
# annoFullTrack = annoFullTrack,
# exonFullTrack = exonFullTrack
# ))
# }
# bench <- system.time({
if(usePlotly) {
if(is(oldP, "covPlotly")) {
p <- oldP
p@covTrack <- covTrack
p@diffTrack <- diffTrack
p@annoTrack <- annoFullTrack
p@exonTrack <- exonFullTrack
} else {
p <- covPlotly(
covTrack = covTrack,
diffTrack = diffTrack,
annoTrack = annoFullTrack
)
}
p@args[["resolution"]] <- resolution
p@args[["reservedCoords"]] <- reservedCoords
p@args[["exonRanges"]] <- exons_gr
p <- .pV_assemble_covPlotly2(p, originalVerticalLayout)
p@fig[[2]] <- p@fig[[2]] %>% layout(
dragmode = "pan",
xaxis = list(
range = c(plotViewStart, plotViewEnd)
)
)
return(p)
} else {
if(showExonRanges) {
p <- .pV_assemble_ggplot(
covTrack, diffTrack, annoSubTrack, exonFullTrack,
verticalLayout, horizontalLayout
)
show(p)
return(exons_gr)
} else {
p <- .pV_assemble_ggplot(
covTrack, diffTrack, annoSubTrack, annoFullTrack,
verticalLayout, horizontalLayout
)
return(p)
}
}
# })
# message("Plotly generation time:")
# print(bench)
return(p)
}
################################################################################
# Given a range, resolution, and must-include coords
# - returns a set of coordinates to keep in final plot
.pV_getAllowedCoords <- function(
viewStart, viewEnd, reservedCoords, resolution
) {
rC <- reservedCoords[reservedCoords %in% seq(viewStart, viewEnd)]
okCoords <- NULL
if(resolution > 0 & !is.null(rC)) {
okCoords <- round(seq(
from = min(viewStart, viewEnd),
to = max(viewStart, viewEnd),
length.out = max(2, resolution - length(rC))
))
okCoords <- sort(unique(c(okCoords, rC)))
}
return(okCoords)
}
# Returns view_start / view_end of coords centered around event
.pV_getEventCoords <- function(
x, args,
zoom_factor = 0.2,
bases_flanking = 100
) {
# sanitise parameters
if(!is.numeric(bases_flanking) || bases_flanking < 0) bases_flanking <- 0
if(!is.numeric(zoom_factor) || zoom_factor < 0) zoom_factor <- 0
if(!("EventRegion" %in% names(args))) return(NULL)
EventRegion <- args[["EventRegion"]]
grEvent <- coord2GR(EventRegion)
eventStart <- start(grEvent)
eventEnd <- end(grEvent)
view_length <- eventEnd - eventStart
view_center <- (eventEnd + eventStart) / 2
new_view_length <- view_length * 3^zoom_factor + 2 * bases_flanking
view_start <- round(view_center - new_view_length / 2)
view_end <- round(view_center + new_view_length / 2)
seqInfo <- args[["seqInfo"]][args[["view_chr"]]]
seqmax <- GenomeInfoDb::seqlengths(seqInfo)
if (view_end > seqmax) {
view_end <- seqmax - 1
view_start <- view_end - new_view_length
if(view_start < 1) view_start <- 1
}
args[["view_start"]] <- view_start
args[["view_end"]] <- view_end
return(args)
}
################################################################################
# Returns reduced.DT and transcript.DT after filtering annotations by
# given criteria
.pV_filterTranscripts <- function(
x, args,
filterByTranscripts,
filterByEventTranscripts,
filterByExpressedTranscripts,
aggJuncCoords
) {
# Filter annotations by observed coords
reduced.DT <- copy(x@annotation[["reduced.DT"]])
transcripts.DT <- copy(x@annotation[["transcripts.DT"]])
transcriptFiltered <- FALSE
if(is_valid(filterByTranscripts)) {
if(
all(filterByTranscripts %in%
c(
transcripts.DT$transcript_id,
transcripts.DT$transcript_name
)
)
) {
transcripts.DT <- transcripts.DT[
get("transcript_name") %in% filterByTranscripts |
get("transcript_id") %in% filterByTranscripts
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
transcriptFiltered <- TRUE
}
}
if(!transcriptFiltered & filterByEventTranscripts) {
reduced.DT.HL <- reduced.DT[get("highlight") != "0"]
reduced.DT.nonHL <- reduced.DT[get("highlight") == "0"]
reduced.DT.nonHL <- reduced.DT.nonHL[
grepl("novel", get("transcript_id"))
]
if(nrow(reduced.DT.HL) != 0) {
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% reduced.DT.HL$transcript_id &
# remove non-highlighted novel elements
!(get("transcript_id") %in% reduced.DT.nonHL$transcript_id)
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
transcriptFiltered <- TRUE
}
}
if(!transcriptFiltered && filterByExpressedTranscripts) {
reduced.DT.introns <- reduced.DT[get("type") == "intron"]
expressed.introns <- reduced.DT.introns[
get("start") %in% aggJuncCoords[["jc_start"]] &
get("end") %in% aggJuncCoords[["jc_end"]]
]
nonexpressed.novel <- reduced.DT.introns[
!(get("start") %in% aggJuncCoords[["jc_start"]]) |
!(get("end") %in% aggJuncCoords[["jc_end"]])
]
nonexpressed.novel <- nonexpressed.novel[
grepl("novel", get("transcript_id"))
]
if(nrow(expressed.introns) > 0) {
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% expressed.introns$transcript_id &
# remove non-highlighted novel elements
!(get("transcript_id") %in% nonexpressed.novel$transcript_id)
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
transcriptFiltered <- TRUE
}
}
return(list(
reduced.DT = reduced.DT,
transcripts.DT = transcripts.DT
))
}
################################################################################
# Given a list of tracks (in args)
# - returns a list of junction start/end coordinates of expressed transcripts
.pV_getJuncCoords <- function(
x, args,
normalizeCoverage = FALSE,
range_gr,
junctionThreshold = 0.01
) {
plotMeanCov <- FALSE
if("condition" %in% names(args)) plotMeanCov <- TRUE
if(normalizeCoverage) {
if(length(x@norm_cov) == 0 || plotMeanCov) {
normalizeCoverage <- FALSE
}
}
trackList <- args[["trackList"]]
fetch_gr <- range_gr
start(fetch_gr) <- max(args[["limit_start"]],
start(range_gr) - width(range_gr))
end(fetch_gr) <- min(args[["limit_end"]],
end(range_gr) + width(range_gr))
jc_start <- c()
jc_end <- c()
for(track_id in seq_len(length(trackList))) {
idx <- trackList[[track_id]]
for(elem in args[["tracks"]][idx]) {
if(plotMeanCov) {
junc <- x@junc_PSI[[elem]]
} else if(normalizeCoverage) {
junc <- x@norm_junc[[elem]]
} else {
junc <- x@junc[[elem]]
}
OL <- findOverlaps(junc, fetch_gr)
if(length(from(OL)) > 0) {
junc <- junc[unique(from(OL))]
} else {
junc <- NULL
}
if(!is.null(junc)) {
isMean <- ("mean" %in% names(mcols(junc)))
if(isMean) {
junc <- junc[mcols(junc)$mean >= junctionThreshold]
} else if(normalizeCoverage) {
junc <- junc[mcols(junc)$count >= junctionThreshold]
} else {
cov <- x@cov[[elem]]
junc <- junc[mcols(junc)$count >=
junctionThreshold * max(cov$depth)
]
}
if(length(junc) > 0) {
jc_start <- c(jc_start, start(junc))
jc_end <- c(jc_end, end(junc))
}
}
}
}
return(list(
jc_start = jc_start, jc_end = jc_end
))
}
# Creates coverage plot tracks
.pV_plotCoverageTracks <- function(
x, args,
normalizeCoverage,
range_gr,
plotJunctions = TRUE, plotJuncPSI = FALSE,
junctionThreshold = 0.01,
usePlotly = FALSE,
use_DT = FALSE
) {
plotMeanCov <- FALSE
if("condition" %in% names(args)) plotMeanCov <- TRUE
rb <- "none"
if(plotMeanCov) rb <- args[["ribbon_mode"]]
if(normalizeCoverage) {
if(length(x@norm_cov) == 0 || plotMeanCov) {
normalizeCoverage <- FALSE
}
}
# junction processing is identical
trackName <- "Coverage"
trackList <- args[["trackList"]]
# plot view x cartesian coordinates
if(args[["reverseGenomeCoords"]]) {
plotViewStart <- end(range_gr)
plotViewEnd <- start(range_gr)
} else {
plotViewStart <- start(range_gr)
plotViewEnd <- end(range_gr)
}
# fetchViewStart and fetchViewEnd determine how much coverage to fetch
# best to fetch a bit more than we need
fetch_gr <- range_gr
start(fetch_gr) <- max(args[["limit_start"]],
start(range_gr) - width(range_gr))
end(fetch_gr) <- min(args[["limit_end"]],
end(range_gr) + width(range_gr))
DT_list <- list()
DTJn_list <- list()
for(track_id in seq_len(length(trackList))) {
idx <- trackList[[track_id]]
for(elem in args[["tracks"]][idx]) {
if(plotJunctions) {
if(plotMeanCov) {
if(plotJuncPSI) {
junc <- x@junc_PSI[[elem]]
} else {
junc <- x@junc_stats[[elem]]
}
} else if(normalizeCoverage) {
junc <- x@norm_junc[[elem]]
} else {
junc <- x@junc[[elem]]
}
OL <- findOverlaps(junc, fetch_gr)
if(length(from(OL)) > 0) {
junc <- junc[unique(from(OL))]
} else {
junc <- NULL
}
} else {
junc <- NULL
}
if(plotMeanCov) {
groupName <- paste(x@args[["condition"]], elem)
DT <- x@cov_stats[[elem]]
DT_sub <- DT[,seq_len(2)] # only 2 columns
DTJn <- c()
if(!is.null(junc)) {
DTJn <- .pV_sashimiArc(
junc,
arcHeight = 0.1 * max(DT$mean),
junctionThreshold = junctionThreshold * max(DT$mean)
)
}
} else {
groupName <- elem
if(normalizeCoverage) {
DT <- x@norm_cov[[elem]]
} else {
DT <- x@cov[[elem]]
}
DTJn <- c()
if(!is.null(junc)) {
DTJn <- .pV_sashimiArc(
junc,
arcHeight =
0.1 * max(DT$depth),
junctionThreshold = junctionThreshold *
max(DT$depth)
)
}
}
DT$group <- groupName
DT$covTrack <- names(trackList)[track_id]
DT_list[[groupName]] <- DT
if(!is.null(junc) && nrow(DTJn) > 0) {
DTJn$group <- groupName
DTJn$covTrack <- names(trackList)[track_id]
DTJn_list[[groupName]] <- DTJn
}
}
}
df_all <- rbindlist(DT_list)
dfJn_all <- rbindlist(DTJn_list)
if(!use_DT) {
df_all <- as.data.frame(df_all)
dfJn_all <- as.data.frame(dfJn_all)
}
# group annotation
df_all$info <- ""
# junc processing - common
dfJnSum_list <- list()
if(plotJunctions) {
dtJn_all <- as.data.table(dfJn_all)
allGroups <- unique(dtJn_all$group)
for(elem in allGroups) {
dtJn <- copy(dtJn_all[get("group") == elem])
dtJn <- dtJn[
get("x") >= start(range_gr) &
get("x") <= end(range_gr)
]
dtJn[, c("xlabel", "ylabel") := list(
mean(get("x")), mean(get("ytmp"))),
by = "info"
]
dtJn <- unique(
dtJn[,
c("info", "value", "xlabel", "ylabel", "group", "covTrack"),
with = FALSE
]
)
dfJnSum_list[[elem]] <- dtJn
}
}
if(use_DT) {
dfJnSum_all <- rbindlist(dfJnSum_list)
} else {
dfJnSum_all <- as.data.frame(rbindlist(dfJnSum_list))
}
# factor covTrack
df_all$covTrack <- factor(df_all$covTrack, names(trackList))
if(nrow(dfJn_all) > 0) {
dfJn_all$covTrack <- factor(dfJn_all$covTrack, names(trackList))
}
if(nrow(dfJnSum_all) > 0) {
dfJnSum_all$covTrack <- factor(dfJnSum_all$covTrack, names(trackList))
}
if(!usePlotly) {
p <- .pV_seed_ggplot_covTrack(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
)
return(list(
ggplot = p,
dfJn_all = dfJn_all
))
} else {
p <- .pV_seed_plotly_covTrack(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
)
return(p)
}
}
# Creates differential plot tracks
.pV_plotDiffTracks <- function(
x, args, diffList,
diff_stat,
range_gr,
usePlotly = FALSE
) {
# plot view x cartesian coordinates
if(args[["reverseGenomeCoords"]]) {
plotViewStart <- end(range_gr)
plotViewEnd <- start(range_gr)
} else {
plotViewStart <- start(range_gr)
plotViewEnd <- end(range_gr)
}
# fetchViewStart and fetchViewEnd determine how much coverage to fetch
# best to fetch a bit more than we need
fetch_gr <- range_gr
start(fetch_gr) <- max(args[["limit_start"]],
start(range_gr) - width(range_gr))
end(fetch_gr) <- min(args[["limit_end"]],
end(range_gr) + width(range_gr))
# Assume diffList is a list
if(is.null(names(diffList))) {
listNames <- character(length(diffList))
for(i in seq_len(length(diffList))) {
listNames[i] <- paste(diffList[[i]], collapse = "_vs_")
}
names(diffList) <- listNames
} else {
listNames <- names(diffList)
for(i in seq_len(listNames)) {
if(listNames[i] == "") {
listNames[i] <- paste(diffList[[i]], collapse = "_vs_")
}
}
names(diffList) <- listNames
}
#
df_list <- list()
for(i in seq_len(length(diffList))) {
df <- c()
diffItem <- diffList[[i]]
diffName <- names(diffList)[i]
if(length(diffItem) == 2) {
if(
is.vector(diffItem, "character") &&
all(diffItem %in% names(x@norm_cov))
) {
# Assume these are track names
cov1 <- x@norm_cov[[diffItem[1]]]
cov2 <- x@norm_cov[[diffItem[2]]]
} else if(
is.vector(diffItem, "numeric") &&
all(diffItem %in% seq_len(length(x@norm_cov)))
) {
cov1 <- x@norm_cov[[diffItem[1]]]
cov2 <- x@norm_cov[[diffItem[2]]]
} else {
next # ignore this
}
} else {
next # ignore this
}
cov1 <- cov1[get("x") >= start(fetch_gr) & get("x") <= end(fetch_gr)]
cov2 <- cov2[get("x") >= start(fetch_gr) & get("x") <= end(fetch_gr)]
if(diff_stat == "t-test") {
df <- .pV_diffStat_ttest(cov1, cov2)
} else {
next
}
if(!is.null(df)) {
df$diffTrack <- diffName
df_list[[diffName]] <- df
}
}
df_all <- as.data.frame(rbindlist(df_list))
if(nrow(df_all) == 0) return(NULL)
# group annotation
df_all$info <- ""
# factor covTrack
df_all$diffTrack <- factor(df_all$diffTrack, names(diffList))
if(!usePlotly) {
p <- .pV_seed_ggplot_diffTrack(
plotViewStart, plotViewEnd,
df_all, diff_stat
)
} else {
p <- .pV_seed_plotly_diffTrack(
plotViewStart, plotViewEnd,
df_all, diff_stat
)
}
return(p)
}
################################################################################
# Creates a ggplot of coverage tracks
.pV_seed_ggplot_covTrack <- function(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
) {
covTracks <- levels(df_all$covTrack)
plotList <- list()
viewMin <- min(plotViewStart, plotViewEnd)
viewMax <- max(plotViewStart, plotViewEnd)
for(i in seq_len(length(covTracks))) {
covTrack = covTracks[i]
df <- df_all[df_all$covTrack == covTrack,]
dfJn <- dfJn_all[dfJn_all$covTrack == covTrack,]
nGroups <- length(unique(df$group))
# Filter by [plotViewStart, plotViewEnd]
df <- df[df$x >= viewMin & df$x <= viewMax,]
dfJn <- dfJn[dfJn$x >= viewMin & dfJn$x <= viewMax,]
p <- ggplot(df) + geom_hline(yintercept = 0)
# plot ribbon (group cov only)
if(plotMeanCov && rb %in% c("ci", "sd", "sem")) {
p <- p +
geom_ribbon(data = df, alpha = 0.2,
aes(
x = get("x"), y = get("mean"),
ymin = get("mean") - get(rb),
ymax = get("mean") + get(rb),
group = get("group"),
fill = get("group")
)
)
}
# plot line
if(plotMeanCov) {
p <- p + geom_line(aes(
x = get("x"), y = get("mean"),
group = get("group"), color = get("group")
))
} else {
p <- p + geom_line(aes(
x = get("x"), y = get("depth"),
group = get("group"), color = get("group")
))
}
if(nGroups == 1) {
p <- p + scale_color_manual(values = "black") +
scale_fill_manual(values = "black")
}
# plot junctions
if(plotJunctions & nGroups == 1) {
dfJn <- dfJn_all[dfJn_all$covTrack == covTrack,]
p <- p +
geom_line(
data = dfJn,
aes(x = get("x"), y = get("yarc"), group = get("info")),
color = "darkred"
)
# NB sashimi arc numbers are plotted later
}
if(nGroups > 1) {
p <- p + theme_white_legend +
theme(
axis.title.x = element_blank(),
legend.title = element_blank()
) +
labs(x = "", y = covTrack)
} else {
p <- p + theme_white +
theme(axis.title.x = element_blank()) +
labs(x = "", y = covTrack)
}
plotList[[i]] <- p
}
return(plotList)
}
# Creates plotly (data) for coverage tracks
.pV_seed_plotly_covTrack <- function(
plotViewStart, plotViewEnd,
df_all, dfJn_all, dfJnSum_all,
plotJunctions, plotMeanCov, rb
) {
covTrack <- levels(df_all$covTrack)
# nested list of data
dataList <- list()
# y axis range, title
layoutList <- list()
for(i in seq_len(length(covTrack))) {
track <- covTrack[i]
df <- df_all[df_all$covTrack == track,]
dfJn <- dfJn_all[dfJn_all$covTrack == track,]
dfJnSum <- dfJnSum_all[dfJnSum_all$covTrack == track,]
groupNames <- unique(df$group)
nGroups <- length(groupNames)
if(nrow(df) > 0) {
df$group <- factor(df$group, groupNames)
}
if(nrow(dfJn) > 0) {
dfJn$group <- factor(dfJn$group, groupNames)
}
if(nrow(dfJnSum) > 0) {
dfJnSum$group <- factor(dfJnSum$group, groupNames)
}
dataList[[i]] <- list()
dataList[[i]][[1]] <- list() # always have junction track
for(j in seq_len(nGroups)) {
dataList[[i]][[2*j]] <- list()
dataList[[i]][[2*j + 1]] <- list()
}
if(plotJunctions & nGroups == 1) {
# sashimi arcs
plotJuncData <- .plotly_makeJuncCurveData(
dfJn$x, dfJn$yarc, dfJn$info
)
dataList[[i]][[1]] <- list(
x = plotJuncData$x, y = plotJuncData$y,
text = plotJuncData$text
)
# sashimi numbers
dataList[[i]][[2]] <- list(
x = dfJnSum$xlabel, y = dfJnSum$ylabel,
text = dfJnSum$value
)
}
for(j in seq_len(nGroups)) {
df_sub <- as.data.frame(df[df$group == levels(df$group)[j],])
if("mean" %in% colnames(df)) {
dataList[[i]][[2*j + 1]] <- list(
x = df_sub$x, y = df_sub$mean, text = df_sub$info
)
dataList[[i]][[2*j + 1]][["hovertemplate"]] <- paste0(
"<b>", levels(df$group)[j], "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>Norm-Depth (mean):</i> %{y:.3f}<br>",
"<extra></extra>"
)
if(rb %in% c("ci", "sd", "sem")) {
dataList[[i]][[2*j + 2]] <- list(
x = c(df_sub$x, rev(df_sub$x)),
y = c(
df_sub$mean + unname(unlist(df_sub[,rb])),
rev(df_sub$mean - unname(unlist(df_sub[,rb])))
),
text = c(df_sub$info, rev(df_sub$info))
)
dataList[[i]][[2*j + 2]][["hovertemplate"]] <- paste0(
"<b>", levels(df$group)[j], "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>Norm-Depth (", rb,"):</i> %{y:.3f}<br>",
"<extra></extra>"
)
} else {
dataList[[i]][[2*j + 2]] <- list()
}
} else {
dataList[[i]][[2*j + 1]] <- list(
x = df_sub$x, y = df_sub$depth, text = df_sub$info
)
dataList[[i]][[2*j + 1]][["hovertemplate"]] <- paste0(
"<b>", levels(df$group)[j], "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>Depth:</i> %{y}<br>",
"<extra></extra>"
)
dataList[[i]][[2*j + 2]] <- list()
}
if(nGroups > 1) {
dataList[[i]][[2*j + 1]][["showlegend"]] <- TRUE
dataList[[i]][[2*j + 1]][["name"]] <- levels(df$group)[j]
}
}
# what is yrange?
df_sub <- df[
df$x >= min(plotViewStart, plotViewEnd) &
df$x <= max(plotViewStart, plotViewEnd),
]
if("mean" %in% colnames(df_sub)) {
ymax <- max(df_sub$mean)
} else {
ymax <- max(df_sub$depth)
}
layoutList[[i]] <- list(
yrange = c(0, ymax * 1.2),
title = track
)
}
return(list(
dataList = dataList,
layoutList = layoutList
))
}
# Creates a ggplot of differential tracks
.pV_seed_ggplot_diffTrack <- function(
plotViewStart, plotViewEnd, df_all, diff_stat
) {
diffTracks <- levels(df_all$diffTrack)
plotList <- list()
for(i in seq_len(length(diffTracks))) {
diffTrack = diffTracks[i]
df <- df_all[df_all$diffTrack == diffTrack,]
p <- ggplot(df) + geom_hline(yintercept = 0)
# plot line
p <- p + geom_line(aes(x = get("x"), y = get("stat")),
color = "black")
# plot y axis and format
p <- p + theme_white_legend + labs(y = "-log10 P")
p <- p +
theme(axis.title.x = element_blank()) +
labs(
x = "", y = diffTrack
)
plotList[[i]] <- p
}
return(plotList)
}
# Creates a plotly (data) of diff tracks
.pV_seed_plotly_diffTrack <- function(
plotViewStart, plotViewEnd, df_all, diff_stat
) {
diffTrack <- levels(df_all$diffTrack)
# nested list of data
dataList <- list()
# y axis range, title
layoutList <- list()
for(i in seq_len(length(diffTrack))) {
track <- diffTrack[i]
df <- df_all[df_all$diffTrack == track,]
dataList[[i]] <- list()
dataList[[i]][[1]] <- list(
x = df$x, y = df$stat, text = df$info
)
dataList[[i]][[1]][["hovertemplate"]] <- paste0(
"<b>", track, "</b><br>",
"<i>Coordinate:</i> %{x}<br>",
"<i>-log10(P):</i> %{y}<br>",
"<extra></extra>"
)
df_sub <- df[
df$x >= min(plotViewStart, plotViewEnd) &
df$x <= max(plotViewStart, plotViewEnd),
]
ymax <- max(df_sub$stat)
layoutList[[i]] <- list(
yrange = c(0, ymax * 1.2),
title = track
)
}
return(list(
dataList = dataList,
layoutList = layoutList
))
}
################################################################################
# Plots an annotation track (or sub-track)
.pV_plotRef <- function(
DTplotlist, view_gr, showPlotRanges = GRanges(),
reverseGenomeCoords = FALSE,
add_information = TRUE,
usePlotly = FALSE,
instantPlotly = FALSE
) {
view_chr <- as.character(seqnames(view_gr))
view_start <- start(view_gr)
view_end <- end(view_gr)
group.DT <- copy(DTplotlist$group.DT)
# use one paste function
if(add_information) {
group.DT[, c("left_extra", "right_extra") := list("","")]
if(reverseGenomeCoords) {
group.DT[get("strand") == "-", c("right_extra") := " -->"]
group.DT[get("strand") == "+", c("left_extra") := "<-- "]
} else {
group.DT[get("strand") == "+", c("right_extra") := " -->"]
group.DT[get("strand") == "-", c("left_extra") := "<-- "]
}
group.DT[, c("display_name") := paste0(
get("left_extra"),
get("group_name"),
" (", get("strand"), ") ",
get("group_biotype"),
get("right_extra")
)]
group.DT[, c("left_extra", "right_extra") := list(NULL,NULL)]
group.DT[, c("disp_x") := 0.5 * (get("start") + get("end"))]
group.DT[get("start") < view_start & get("end") > view_start,
c("disp_x") := 0.5 * (view_start + get("end"))]
group.DT[get("end") > view_end & get("start") < view_end,
c("disp_x") := 0.5 * (get("start") + view_end)]
group.DT[get("start") < view_start & get("end") > view_end,
c("disp_x") := 0.5 * (view_start + view_end)]
}
reduced <- copy(DTplotlist$reduced.DT)
reduced <- reduced[!is.na(reduced$plot_level)]
condense_this <- DTplotlist$condense_this
exonRanges <- DTplotlist$exonRanges
# Hover Text annotation
reduced[, c("Information") := ""]
reduced[get("type") %in% c("exon", "CDS"), c("Information") := paste(
paste(get("transcript_id"), "exon", get("aux_id")),
paste0("(", get("feature_id"), ")"),
paste0(get("seqnames"), ":", get("start"), "-", get("end"), "/",
get("strand")),
sep = "\n"
)]
reduced[get("type") == "intron", c("Information") := paste(
get("feature_id"),
paste0(get("seqnames"), ":", get("start"), "-", get("end"), "/",
get("strand")),
sep = "\n"
)]
reduced <- as.data.frame(reduced)
reducedIntrons <- data.frame()
if (sum(reduced$type == "intron", na.rm = TRUE) > 0) {
reducedIntrons <- reduced[reduced$type == "intron", ]
reducedIntronsExpanded <- c()
for(i in seq_len(nrow(reducedIntrons))) {
reducedIntronsExpanded <- rbind(reducedIntronsExpanded, data.frame(
start = seq(reducedIntrons$start[i], reducedIntrons$end[i],
length.out = 10),
end = seq(reducedIntrons$start[i], reducedIntrons$end[i],
length.out = 10),
plot_level = reducedIntrons$plot_level[i],
highlight = reducedIntrons$highlight[i],
Information = reducedIntrons$Information[i]
))
}
col_highlights <- sort(unique(reducedIntronsExpanded$highlight))
}
nonIntrons <- reduced[reduced$type != "intron", ]
fill_highlights <- sort(unique(nonIntrons$highlight))
if(!reverseGenomeCoords) {
plotViewStart <- view_start
plotViewEnd <- view_end
} else {
plotViewStart <- view_end
plotViewEnd <- view_start
}
if(!usePlotly) {
p <- .pV_seed_ggplot_annoTrack(
plotViewStart, plotViewEnd, view_chr,
reducedIntronsExpanded, nonIntrons,
group.DT, condense_this, add_information,
showPlotRanges
)
} else if(instantPlotly) {
p <- .pV_seed_instant_plotly_annoTrack(
plotViewStart, plotViewEnd, view_chr,
reducedIntrons, nonIntrons,
group.DT, condense_this, add_information
)
} else {
p <- .pV_seed_plotly_annoTrack(
plotViewStart, plotViewEnd, view_chr,
reducedIntrons, nonIntrons,
group.DT, condense_this, add_information
)
}
return(p)
}
# Plots an annotation track showing only named exons
.pV_plotExons <- function(
DTplotlist, view_gr, exons_gr,
reverseGenomeCoords = FALSE,
add_information = TRUE,
usePlotly = FALSE,
instantPlotly = FALSE
) {
view_chr <- as.character(seqnames(view_gr))
view_start <- start(view_gr)
view_end <- end(view_gr)
exons <- as.data.frame(exons_gr)
if(!reverseGenomeCoords) {
plotViewStart <- view_start
plotViewEnd <- view_end
} else {
plotViewStart <- view_end
plotViewEnd <- view_start
}
if(!usePlotly) {
p <- .pV_seed_ggplot_HLexons(
plotViewStart, plotViewEnd, view_chr, exons
)
} else if(instantPlotly) {
p <- .pV_seed_instant_plotly_HLexons(
plotViewStart, plotViewEnd, view_chr, exons
)
} else {
p <- .pV_seed_plotly_HLexons(
plotViewStart, plotViewEnd, view_chr, exons
)
}
return(p)
}
# Convert highlight codes to colors
.pV_highlight_to_colors <- function(highlights, usePlotly = FALSE) {
if(is.null(highlights)) return(NULL)
if(usePlotly) {
highlights <- sub("0", "rgba(0,0,0,1)", highlights)
highlights <- sub("1", "rgba(0,0,255,1)", highlights)
highlights <- sub("2", "rgba(255,0,0,1)", highlights)
highlights <- sub("3", "rgba(255,0,255,1)", highlights)
} else {
highlights <- sub("0", "black", highlights)
highlights <- sub("1", "blue", highlights)
highlights <- sub("2", "red", highlights)
highlights <- sub("3", "purple", highlights)
}
return(highlights)
}
# Creates a ggplot of annotation track
.pV_seed_ggplot_annoTrack <- function(
plotViewStart, plotViewEnd, view_chr,
introns, exons, group.DT, condense_this, add_information,
showPlotRanges
) {
col_highlights <- sort(unique(introns$highlight))
fill_highlights <- sort(unique(exons$highlight))
col_highlights <- .pV_highlight_to_colors(col_highlights)
fill_highlights <- .pV_highlight_to_colors(fill_highlights)
anno <- NULL
if(add_information) {
if (condense_this == TRUE) {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.5 + 0.3 *
runif(rep(1, nrow(group.DT))),
Information = group.DT$display_name
)
} else {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.4,
Information = group.DT$display_name
)
}
}
if (nrow(group.DT) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(group.DT$plot_level)
}
p <- ggplot()
if(nrow(introns) > 0) {
p <- p + geom_line(data = introns,
aes(x = get("start"), y = get("plot_level"),
color = get("highlight"), group = get("Information"))
)
}
if(nrow(exons) > 0) {
p <- p + geom_rect(data = exons,
aes(xmin = get("start"), xmax = get("end"),
ymin = get("plot_level") - 0.1 -
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
ymax = get("plot_level") + 0.1 +
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
fill = get("highlight")
)
)
}
if(length(showPlotRanges) > 0) {
ref_ymin <- min(layer_scales(p)$y$range$range)
ref_ymax <- max(layer_scales(p)$y$range$range)
dfRanges <- as.data.frame(showPlotRanges)
dfRanges$rangeName <- as.character(seq_len(length(showPlotRanges)))
p <- p + geom_rect(data = dfRanges,
mapping = aes(xmin = get("start"), xmax = get("end"),
ymin = ref_ymin - 0.5, ymax = ref_ymax + 0.5
), fill = NA, color = "red"
) + geom_text(data = dfRanges,
aes(
x = 0.5 * (get("start") + get("end")),
y = ref_ymin - 1, label = get("rangeName")
)
)
}
if(!is.null(col_highlights))
p <- p + scale_color_manual(values = col_highlights)
if(!is.null(fill_highlights))
p <- p + scale_fill_manual(values = fill_highlights)
p <- p + theme_white_legend_plot_track +
theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank()
)
if(!is.null(anno)) {
p <- p + geom_text(
data = anno,
aes(x = get("x"), y = get("y"), label = get("Information"))
)
}
p <- p + theme(legend.position = "none") +
labs(x = paste("Chromosome", view_chr))
ref_ymin <- min(layer_scales(p)$y$range$range)
ref_ymax <- max(layer_scales(p)$y$range$range)
p <- p +
scale_x_continuous(labels = label_number(scale_cut = cut_si(""))) +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(ref_ymin - 1, ref_ymax + 1),
expand = FALSE
)
return(p)
}
# Creates a ggplot of exons track
.pV_seed_ggplot_HLexons <- function(
plotViewStart, plotViewEnd, view_chr, exons
) {
fill_highlights <- sort(unique(exons$highlight))
fill_highlights <- .pV_highlight_to_colors(fill_highlights)
if (nrow(exons) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(exons$plot_level)
}
p <- ggplot()
if(nrow(exons) > 0) {
p <- p + geom_rect(data = exons,
aes(xmin = get("start"), xmax = get("end"),
ymin = get("plot_level") - 0.1 -
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
ymax = get("plot_level") + 0.1 +
ifelse(get("type") %in%
c("CDS", "start_codon", "stop_codon"), 0.1, 0),
fill = get("highlight")
)
) + geom_text(
data = exons,
aes(
x = 0.5 * (get("start") + get("end")),
y = get("plot_level") - 0.4,
label = get("exon_name")
)# ,
# hjust=1, angle=45
)
}
p <- p +
scale_fill_manual(values = fill_highlights)
p <- p + theme_white_legend_plot_track +
theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank()
)
p <- p + theme(legend.position = "none") +
labs(x = paste("Chromosome", view_chr))
ref_ymin <- min(layer_scales(p)$y$range$range)
ref_ymax <- max(layer_scales(p)$y$range$range)
p <- p +
scale_x_continuous(labels = label_number(scale_cut = cut_si(""))) +
coord_cartesian(
xlim = c(plotViewStart, plotViewEnd),
ylim = c(ref_ymin - 1, ref_ymax + 1),
expand = FALSE
)
return(p)
}
# Creates an instant plotly of annotation track
.pV_seed_instant_plotly_annoTrack <- function(
plotViewStart, plotViewEnd, view_chr,
introns, exons, group.DT, condense_this, add_information
) {
col_highlights <- sort(unique(introns$highlight))
fill_highlights <- sort(unique(exons$highlight))
col_colors <- .pV_highlight_to_colors(col_highlights)
fill_colors <- .pV_highlight_to_colors(fill_highlights)
anno <- NULL
if(add_information) {
if (condense_this == TRUE) {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.5 + 0.3 *
runif(rep(1, nrow(group.DT))),
Information = group.DT$display_name
)
} else {
anno <- data.frame(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.4,
Information = group.DT$display_name
)
}
}
if (nrow(group.DT) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(group.DT$plot_level)
}
fig <- plot_ly()
if(nrow(introns) > 0) {
intronLineData <- .plotly_makeLineData(
introns$start, introns$plot_level, introns$end, introns$plot_level,
introns$Information, introns$highlight
)
for(i in seq_len(length(col_highlights))) {
hl <- col_highlights[i]
color <- col_colors[i]
fig <- fig %>% add_trace(
data = intronLineData[intronLineData$colorInfo == hl,],
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~text"),
hoveron = "points", hoverinfo = 'text',
line = list(
color = color
),
type = 'scatter', mode = 'lines',
showlegend = FALSE
)
}
}
typeThick <- c("CDS", "start_codon", "stop_codon")
if(nrow(exons) > 0) {
exons_thin <- exons[!(exons$type %in% typeThick), ]
exons_thick <- exons[exons$type %in% typeThick, ]
exons_thin$ymin <- exons_thin$plot_level - 0.15
exons_thin$ymax <- exons_thin$plot_level + 0.15
exons_thick$ymin <- exons_thick$plot_level - 0.3
exons_thick$ymax <- exons_thick$plot_level + 0.3
exons_comb <- rbind(exons_thin, exons_thick)
exons_instr <- .plotly_makeRectData(
exons_comb$start, exons_comb$end,
exons_comb$ymin, exons_comb$ymax,
exons_comb$Information,
exons_comb$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
fig <- fig %>% add_trace(
data = rects,
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~text"),
type = 'scatter', mode = 'lines',
hoveron = "points+fills", hoverinfo = 'text',
line = list(
color = "transparent"
),
fill = "toself",
fillcolor = fill_colors[i],
showlegend = FALSE
)
}
}
if(add_information) {
xtitle <- paste("Chromosome/Scaffold", view_chr)
fig <- fig %>% add_trace(
data = anno,
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~Information"),
hoverinfo = 'text',
type = 'scatter', mode = 'text',
textposition = 'middle',
showlegend = FALSE
)
} else {
xtitle <- ""
}
fig <- fig %>% layout(
dragmode = "pan",
xaxis = list(
range = c(plotViewStart, plotViewEnd),
title = xtitle
),
yaxis = list(
title = "",
range = c(0, 1 + max_plot_level)# ,
# fixedrange = TRUE
)
)
return(fig)
}
# Creates a plotly of exons track
.pV_seed_instant_plotly_HLexons <- function(
plotViewStart, plotViewEnd, view_chr, exons
) {
fill_highlights <- sort(unique(exons$highlight))
fill_colors <- .pV_highlight_to_colors(fill_highlights)
if (nrow(exons) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(exons$plot_level)
}
fig <- plot_ly()
if(nrow(exons) > 0) {
exons$ymin <- exons$plot_level - 0.15
exons$ymax <- exons$plot_level + 0.15
exons$xdisp <- 0.5 * (exons$start + exons$end)
exons$ydisp <- exons$plot_level - 0.4
exons_instr <- .plotly_makeRectData(
exons$start, exons$end,
exons$ymin, exons$ymax,
exons$exon_name,
exons$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
fig <- fig %>% add_trace(
data = rects,
x = as.formula("~x"),
y = as.formula("~y"),
text = as.formula("~text"),
type = 'scatter', mode = 'lines',
hoveron = "points+fills", hoverinfo = 'text',
line = list(
color = "transparent"
),
fill = "toself",
fillcolor = fill_colors[i],
showlegend = FALSE
)
fig <- fig %>% add_trace(
data = exons,
x = as.formula("~xdisp"),
y = as.formula("~ydisp"),
text = as.formula("~exon_name"),
hoverinfo = 'text',
type = 'scatter', mode = 'text',
textposition = 'middle',
showlegend = FALSE
)
}
}
xtitle <- paste("Chromosome/Scaffold", view_chr)
fig <- fig %>% layout(
dragmode = "pan",
xaxis = list(
range = c(plotViewStart, plotViewEnd),
title = xtitle
),
yaxis = list(
title = "",
range = c(0, 1 + max_plot_level),
fixedrange = TRUE
)
)
return(fig)
}
# Creates a plotly (data) of annotation track
.pV_seed_plotly_annoTrack <- function(
plotViewStart, plotViewEnd, view_chr,
introns, exons, group.DT, condense_this, add_information
) {
dataList <- list()
# x,y axis range, title
layoutList <- list()
col_highlights <- fill_highlights <- c("0", "1", "2", "3")
if(nrow(introns) > 0) {
intronLineData <- .plotly_makeLineData(
introns$start, introns$plot_level, introns$end, introns$plot_level,
introns$Information, introns$highlight
)
for(i in seq_len(length(col_highlights))) {
dataList[[i]] <- list()
hl <- col_highlights[i]
data_df <- intronLineData[intronLineData$colorInfo == hl,]
if(nrow(data_df) > 0) {
dataList[[i]] <- list(
x = data_df$x, y = data_df$y, text = data_df$text
)
}
}
} else {
for(i in seq_len(length(col_highlights))) {
dataList[[i]] <- list()
}
}
typeThick <- c("CDS", "start_codon", "stop_codon")
if(nrow(exons) > 0) {
exons_thin <- exons[!(exons$type %in% typeThick), ]
exons_thick <- exons[exons$type %in% typeThick, ]
exons_thin$ymin <- exons_thin$plot_level - 0.15
exons_thin$ymax <- exons_thin$plot_level + 0.15
exons_thick$ymin <- exons_thick$plot_level - 0.3
exons_thick$ymax <- exons_thick$plot_level + 0.3
exons_comb <- rbind(exons_thin, exons_thick)
exons_instr <- .plotly_makeRectData(
exons_comb$start, exons_comb$end,
exons_comb$ymin, exons_comb$ymax,
exons_comb$Information,
exons_comb$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
if(nrow(rects) > 0) {
dataList[[i + 4]] <- list(
x = rects$x, y = rects$y, text = rects$text
)
}
}
} else {
for(i in seq_len(length(col_highlights))) {
dataList[[i + 4]] <- list()
}
}
if(add_information) {
if (condense_this == TRUE) {
dataList[[9]] <- list(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.5 + 0.3 *
runif(rep(1, nrow(group.DT))),
text = group.DT$display_name
)
} else {
dataList[[9]] <- list(
x = group.DT$disp_x,
y = group.DT$plot_level - 0.4,
text = group.DT$display_name
)
}
} else {
dataList[[9]] <- list()
}
if(add_information) {
xtitle <- paste("Chromosome/Scaffold", view_chr)
} else {
xtitle <- ""
}
if (nrow(group.DT) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(group.DT$plot_level)
}
layoutList <- list(
xrange = c(plotViewStart, plotViewEnd),
xtitle = xtitle,
yrange = c(0, 1 + max_plot_level)
)
return(list(
dataList = dataList,
layoutList = layoutList
))
}
# Creates a plotly (data) of exons track
.pV_seed_plotly_HLexons <- function(
plotViewStart, plotViewEnd, view_chr, exons
) {
dataList <- list()
# x,y axis range, title
layoutList <- list()
col_highlights <- fill_highlights <- c("0", "1", "2", "3")
for(i in seq_len(length(col_highlights))) {
dataList[[i]] <- list()
}
typeThick <- c("CDS", "start_codon", "stop_codon")
if(nrow(exons) > 0) {
exons$ymin <- exons$plot_level - 0.15
exons$ymax <- exons$plot_level + 0.15
exons$xdisp <- 0.5 * (exons$start + exons$end)
exons$ydisp <- exons$plot_level - 0.4
exons_instr <- .plotly_makeRectData(
exons$start, exons$end,
exons$ymin, exons$ymax,
exons$exon_name,
exons$highlight,
resolution = 10
)
for(i in seq_len(length(fill_highlights))) {
rects <- exons_instr[
exons_instr$colorInfo == fill_highlights[i],
]
if(nrow(rects) > 0) {
dataList[[i + 4]] <- list(
x = rects$x, y = rects$y, text = rects$text
)
}
}
} else {
for(i in seq_len(length(col_highlights))) {
dataList[[i + 4]] <- list()
}
}
dataList[[9]] <- list(
x = exons$xdisp,
y = exons$ydisp,
text = exons$exon_name
)
xtitle <- paste("Chromosome/Scaffold", view_chr)
if (nrow(exons) == 0) {
max_plot_level <- 1
} else {
max_plot_level <- max(exons$plot_level)
}
layoutList <- list(
xrange = c(plotViewStart, plotViewEnd),
xtitle = xtitle,
yrange = c(0, 1 + max_plot_level)
)
return(list(
dataList = dataList,
layoutList = layoutList
))
}
################################################################################
.seq2 <- Vectorize(seq.default, vectorize.args = c("from", "to", "length.out"))
# Creates arcs for sashimi plot
.pV_sashimiArc <- function(
junc,
arcHeight = 0,
junctionThreshold = 0
) {
if(is.null(junc)) return(NA)
DT <- data.table(
juncName = paste0(
seqnames(junc), ":",
start(junc), "-", end(junc), "/", strand(junc)
),
leftX = start(junc) - 1,
rightX = end(junc) + 1,
leftY = mcols(junc)$leftCoordHeight,
rightY = mcols(junc)$rightCoordHeight
)
isMean <- ("mean" %in% names(mcols(junc)))
if(isMean) {
DT$count <- mcols(junc)$mean
DT$countsd <- mcols(junc)$sd
DT$value <- paste0(
round(100 * DT$count, 1), "+/-",
round(100 * DT$countsd, 1), " %"
)
DT$info <- paste(
paste0("Junction: ", DT$juncName),
paste0("PSI: ", DT$value),
sep = "\n"
)
} else {
DT$count <- mcols(junc)$count
DT$countsd <- 0
DT$value <- round(DT$count,2)
DT$info <- paste(
paste0("Junction: ", DT$juncName),
paste0("Depth: ", DT$count),
sep = "\n"
)
}
DT <- DT[get("count") >= junctionThreshold]
final_DT <- data.table(
coord = rep(DT$juncName, each = 90),
value = rep(DT$value, each = 90),
info = rep(DT$info, each = 90),
x = as.vector(.seq2(DT$leftX, DT$rightX, length.out = 90)),
ytmp = as.vector(.seq2(DT$leftY, DT$rightY, length.out = 90)),
yarc = as.vector(.seq2(DT$leftY, DT$rightY, length.out = 90)) +
rep(arcHeight * sinpi(seq(0,1,length.out = 90)), nrow(DT))
)
return(final_DT)
}
################################################################################
# Check requested tracks are legitimate in getPlotObject
.gPO_check_tracks <- function(
obj, args
) {
# sanity check condition / tracks against colData
colData <- obj@colData
sampleList <- list()
samples_to_check <- c()
# tracks must be a character vector
if(!is.vector(args[["tracks"]], "character")) {
.log(paste(
"In getPlotObject,",
"`tracks` must be a character vector"
))
}
# remove invalid tracks
invalid_tracks <- c()
if("condition" %in% names(args)) {
if(args[["condition"]] %in% colnames(colData)) {
for(cond_elem in args[["tracks"]]) {
# cond_elem is of length 1
if(
any(colData[, args[["condition"]]] == cond_elem)
) {
sampleList[[cond_elem]] <- rownames(colData)[
colData[, args[["condition"]]] == cond_elem]
samples_to_check <- c(samples_to_check,
sampleList[[cond_elem]])
} else {
# ignore track altogether
invalid_tracks <- c(invalid_tracks, cond_elem)
}
}
} else {
.log(paste("In getPlotObject,",
"condition is not a recognised column in colData"
))
}
} else {
for(i in seq_len(length(args[["tracks"]]))) {
for(samp in args[["tracks"]]) {
if(samp %in% rownames(colData)) {
sampleList[[samp]] <- samp
samples_to_check <- c(samples_to_check, samp)
} else {
invalid_tracks <- c(invalid_tracks, samp)
}
}
}
}
# Remove invalid tracks
if(length(invalid_tracks) == length(args[["tracks"]])) {
.log(paste(
"In getPlotObject,",
"all elements in `tracks` are invalid"
))
} else if(length(invalid_tracks) > 0) {
args[["tracks"]] <- args[["tracks"]][
!(args[["tracks"]] %in% invalid_tracks)
]
}
args[["sampleList"]] <- sampleList
args[["samples_to_check"]] <- samples_to_check
return(args)
}
# Check that getPlotObject knows which stranded data to retrieve
.gPO_check_strand <- function(
obj, args
) {
args[["view_strand_jn"]] <- args[["strand"]]
samples_to_check <- args[["samples_to_check"]]
if(args[["strand"]] != "*") {
if(all(obj@normData[["sampleStrand"]][samples_to_check] == -1)) {
if(args[["strand"]] == "+") {
args[["view_strand"]] <- "-"
} else {
args[["view_strand"]] <- "+"
}
}
} else {
args[["view_strand"]] <- "*"
}
return(args)
}
################################################################################
# Get normalized coverage from obj
.gPO_getCoverage <- function(obj, args, normalize = FALSE) {
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand <- args[["view_strand"]]
Event <- args[["Event"]]
outList <- list()
sampleList <- args[["sampleList"]]
if(view_strand == "+") {
cov <- obj@covData$pos
} else if(view_strand == "-") {
cov <- obj@covData$neg
} else {
cov <- obj@covData$uns
}
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
matList <- list()
for(s in samples) {
view <- IRanges::Views(cov[[s]], view_start, view_end)
matList[[s]] <- as.matrix(view[[1]])
}
mat <- do.call(cbind, matList)
colnames(mat) <- samples
# Normalize coverage
if(normalize) {
normFactor <- unlist(obj@normData$norms[Event,samples])
# discard samples with low normalization factors
which_lowNorm <- which(normFactor < 5)
if(length(which_lowNorm) == length(normFactor)) {
outList[[track]] <- data.frame(x = seq(view_start, view_end))
next
} else if(length(which_lowNorm) > 0) {
which_Norm <- which(normFactor > 5)
mat <- mat[, which_Norm, drop = FALSE]
normFactor <- normFactor[which_Norm]
}
mat <- t(t(mat) / normFactor)
}
outList[[track]] <- data.table(x = seq(view_start, view_end))
outList[[track]] <- cbind(outList[[track]], as.data.table(mat))
if(length(samples) == 1) {
colnames(outList[[track]])[2] <- "depth"
}
}
return(outList)
}
# Calculates mean and se/sem/ci of normalized coverages per group
.gPO_covStats <- function(cov) {
out <- list()
outnames <- names(cov)
conf.int <- 0.95
for(i in seq_len(length(cov))) {
covmat <- as.matrix(cov[[outnames[i]]][,-1])
out[[outnames[i]]] <- data.table(
x = cov[[outnames[i]]]$x,
mean = rowMeans(covmat)
)
n <- ncol(covmat)
out[[outnames[i]]]$sd <- rowSds(covmat)
out[[outnames[i]]]$sem <- out[[outnames[i]]]$sd / sqrt(n)
out[[outnames[i]]]$ci <- qt((1 + conf.int) / 2, df = n - 1) *
out[[outnames[i]]]$sem
}
return(out)
}
# Compiles per-condition junction stats (normalized junction counts)
.gPO_normJuncStats <- function(obj, args, cov_stats) {
view_chr <- args[["view_chr"]]
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand_jn <- args[["view_strand_jn"]]
outList <- list()
sampleList <- args[["sampleList"]]
gr_select <- GRanges(
view_chr,
IRanges(view_start, view_end),
view_strand_jn
)
OL <- findOverlaps(obj@juncData$junc_gr, gr_select)
if(length(unique(from(OL))) == 0) return(list())
hits <- sort(unique(from(OL)))
junc_gr <- obj@juncData$junc_gr[hits]
leftCoords <- start(junc_gr) - 1
rightCoords <- end(junc_gr) + 1
lC_inrange <- leftCoords >= view_start & leftCoords <= view_end
rC_inrange <- rightCoords >= view_start & rightCoords <= view_end
Event <- args[["Event"]]
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
if(view_strand_jn == "*") {
junc <- obj@juncData$junc_counts_uns[
hits, samples, drop = FALSE]
} else {
junc <- obj@juncData$junc_counts[
hits, samples, drop = FALSE]
}
cov <- cov_stats[[track]][, seq_len(2), with = FALSE]
normFactor <- unlist(obj@normData$norms[Event,samples])
junc <- t(t(as.matrix(junc)) / normFactor)
# Returned data
outList[[track]] <- .grDT(cbind(
as.data.frame(junc_gr),
data.frame(
mean = rowMeans(junc),
sd = rowSds(junc),
leftCoordHeight = max(cov[,2]),
rightCoordHeight = max(cov[,2])
)
), keep.extra.columns = TRUE)
if(sum(lC_inrange) > 0) {
covL <- unname(unlist(
cov[leftCoords[lC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$leftCoordHeight[lC_inrange] <- covL
}
if(sum(rC_inrange) > 0) {
covR <- unname(unlist(
cov[rightCoords[rC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$rightCoordHeight[rC_inrange] <- covR
}
}
return(outList)
}
.pV_diffStat_ttest <- function(cov1, cov2) {
coords <- sort(intersect(cov1$x, cov2$x))
mat <- cbind(
as.matrix(cov1[get("x") %in% coords, -1]),
as.matrix(cov2[get("x") %in% coords, -1])
)
fac <- factor(
c(
rep("1", ncol(cov1) - 1),
rep("2", ncol(cov2) - 1)
)
)
t_test <- genefilter::rowttests(mat, fac)
ret <- data.frame(
x = coords,
stat = -log10(t_test$p.value)
)
ret$stat[!is.finite(ret$stat)] <- 0
return(ret)
}
.gPO_PSIjuncStats <- function(obj, args, norm_cov) {
view_chr <- args[["view_chr"]]
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand_jn <- args[["view_strand_jn"]]
outList <- list()
sampleList <- args[["sampleList"]]
gr_select <- GRanges(
view_chr,
IRanges(view_start, view_end),
view_strand_jn
)
OL <- findOverlaps(obj@juncData$junc_gr, gr_select)
if(length(unique(from(OL))) == 0) return(list())
hits <- sort(unique(from(OL)))
junc_gr <- obj@juncData$junc_gr[hits]
leftCoords <- start(junc_gr) - 1
rightCoords <- end(junc_gr) + 1
lC_inrange <- leftCoords >= view_start & leftCoords <= view_end
rC_inrange <- rightCoords >= view_start & rightCoords <= view_end
# normalize by starting coordinates:
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
junc <- obj@juncData$junc_PSI[hits, samples, drop = FALSE]
# Returned data
outList[[track]] <- .grDT(cbind(
as.data.frame(junc_gr),
data.frame(
mean = rowMeans(junc),
sd = rowSds(junc),
leftCoordHeight = 1,
rightCoordHeight = 1
)
), keep.extra.columns = TRUE)
cov <- norm_cov[[track]]
if(sum(lC_inrange) > 0) {
covL <- cov[leftCoords[lC_inrange] - cov$x[1] + 1, -1]
mcols(outList[[track]])$leftCoordHeight[lC_inrange] <- rowMeans(covL)
}
if(sum(rC_inrange) > 0) {
covR <- cov[rightCoords[rC_inrange] - cov$x[1] + 1, -1]
mcols(outList[[track]])$rightCoordHeight[rC_inrange] <- rowMeans(covR)
}
}
return(outList)
}
.gPO_getJunc <- function(obj, args, raw_cov, normalized = FALSE) {
view_chr <- args[["view_chr"]]
view_start <- args[["limit_start"]]
view_end <- args[["limit_end"]]
view_strand_jn <- args[["view_strand_jn"]]
outList <- list()
sampleList <- args[["sampleList"]]
gr_select <- GRanges(
view_chr,
IRanges(view_start, view_end),
view_strand_jn
)
OL <- findOverlaps(obj@juncData$junc_gr, gr_select)
if(length(unique(from(OL))) == 0) return(list())
hits <- sort(unique(from(OL)))
junc_gr <- obj@juncData$junc_gr[hits]
leftCoords <- start(junc_gr) - 1
rightCoords <- end(junc_gr) + 1
lC_inrange <- leftCoords >= view_start & leftCoords <= view_end
rC_inrange <- rightCoords >= view_start & rightCoords <= view_end
Event <- NULL
if(normalized) Event <- args[["Event"]]
for(i in seq_len(length(sampleList))) {
track <- names(sampleList)[i]
samples <- sampleList[[i]]
# Only do it if cov is a single track
if(length(samples) == 1) {
if(view_strand_jn == "*") {
junc <- obj@juncData$junc_counts_uns[
hits, samples, drop = FALSE]
} else {
junc <- obj@juncData$junc_counts[
hits, samples, drop = FALSE]
}
cov <- raw_cov[[track]]
# Normalize
if(normalized) {
normFactor <- unlist(obj@normData$norms[Event,samples])
junc <- junc / normFactor
}
# Returned data
outList[[track]] <- .grDT(cbind(
as.data.frame(junc_gr),
data.frame(
count = unname(unlist(junc)),
leftCoordHeight = max(cov[,2]),
rightCoordHeight = max(cov[,2])
)
), keep.extra.columns = TRUE)
if(sum(lC_inrange) > 0) {
covL <- unname(unlist(
cov[leftCoords[lC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$leftCoordHeight[lC_inrange] <- covL
}
if(sum(rC_inrange) > 0) {
covR <- unname(unlist(
cov[rightCoords[rC_inrange] - cov$x[1] + 1, -1]
))
mcols(outList[[track]])$rightCoordHeight[rC_inrange] <- covR
}
} else {
# empty
outList[[track]] <- GRanges()
}
}
return(outList)
}
# list checker for plotViews
.pV_trackList_from_vector <- function(trackVec, args) {
trackList <- list()
if(is.vector(trackVec, "numeric")) {
idx <- trackVec
if(all(idx %in% seq_len(length(args[["tracks"]])))) {
trackList <- lapply(idx, function(i) i)
} else if(any(idx %in% seq_len(length(args[["tracks"]])))) {
idx <- idx[idx %in% seq_len(length(args[["tracks"]]))]
trackList <- lapply(idx, function(i) i)
} else {
.log(paste(
"In plotView,",
"`trackList` contains indices outside the range of tracks",
"found in tracks(covPlotObject)"
))
}
} else if(is.vector(trackVec, "character")) {
idxNames <- trackVec
if(all(idxNames %in% args[["tracks"]])) {
trackList <- lapply(idxNames, function(i) {
which(args[["tracks"]] == i)[1]
})
} else if(any(idxNames %in% args[["tracks"]])) {
idxNames <- idxNames[idxNames %in% args[["tracks"]]]
trackList <- lapply(idxNames, function(i) {
which(args[["tracks"]] == i)[1]
})
} else {
.log(paste(
"In plotView,",
"`trackList` contains tracks",
"not found in tracks(covPlotObject)"
))
}
}
if(!is.null(names(trackVec))) {
names(trackList) <- names(trackVec)
if(any(names(trackList) == "")) {
which_noname <- which(names(trackList) == "")
names(trackList)[which_noname] <- paste(
"Track", as.character(which_noname), sep = "_"
)
}
} else {
for(i in seq_len(length(trackList))) {
elem <- trackList[[i]]
elemNames <- args[["tracks"]][elem]
names(trackList)[i] <- paste(elemNames, collapse = ",")
}
}
return(trackList)
}
.pV_trackList_from_list <- function(trackList, args) {
outList <- list()
for(z in seq_len(length(trackList))) {
trackVec <- trackList[[z]]
if(is.vector(trackVec, "numeric")) {
idx <- trackVec
if(all(idx %in% seq_len(length(args[["tracks"]])))) {
outList[[z]] <- idx
} else if(any(idx %in% seq_len(length(args[["tracks"]])))) {
idx <- idx[idx %in% seq_len(length(args[["tracks"]]))]
outList[[z]] <- idx
} else {
.log(paste(
"In plotView,",
"`trackList` contains indices outside the range of tracks",
"found in tracks(covPlotObject)"
))
}
} else if(is.vector(trackVec, "character")) {
idxNames <- trackVec
if(all(idxNames %in% args[["tracks"]])) {
outList[[z]] <- which(args[["tracks"]] %in% idxNames)
} else if(any(idxNames %in% args[["tracks"]])) {
idxNames <- idxNames[idxNames %in% args[["tracks"]]]
outList[[z]] <- which(args[["tracks"]] %in% idxNames)
} else {
.log(paste(
"In plotView,",
"`trackList` contains tracks",
"not found in tracks(covPlotObject)"
))
}
}
}
# transfer names
if(!is.null(names(trackList))) {
names(outList) <- names(trackList)
if(any(names(outList) == "")) {
which_noname <- which(names(outList) == "")
names(outList)[which_noname] <- paste(
"Track", as.character(which_noname), sep = "_"
)
}
} else {
for(i in seq_len(length(outList))) {
elem <- outList[[i]]
elemNames <- args[["tracks"]][elem]
names(outList)[i] <- paste(elemNames, collapse = ",")
}
}
return(outList)
}
# Make plotly x/y/text data given yarc data
.plotly_makeJuncCurveData <- function(
x, yarc, info
) {
sourceDT <- data.table(
x = x, y = yarc, info = info
)
setorderv(c("info", "x"))
infoVec <- unique(sourceDT$info)
df_list <- list()
for(i in seq_len(length(infoVec))) {
subDT <- sourceDT[get("info") == infoVec[i]]
df_list[[2 * i - 1]] <- data.table(
x = subDT$x,
y = subDT$y,
text = infoVec[i]
)
if(i < length(info)) {
df_list[[2 * i]] <- data.table(
x = NA, y = NA, text = NA
)
}
}
return(as.data.frame(rbindlist(df_list)))
}
# Make plotly line data based on line info
.plotly_makeLineData <- function(
xstart, ystart, xend, yend, info, colorInfo, resolution = 20
) {
df_list <- list()
for(i in seq_len(length(info))) {
df_list[[2 * i - 1]] <- data.table(
x = seq(xstart[i], xend[i], length.out = resolution),
y = seq(ystart[i], yend[i], length.out = resolution),
text = info[i],
colorInfo = colorInfo[i]
)
if(i < length(info)) {
df_list[[2 * i]] <- data.table(
x = NA, y = NA, text = NA, colorInfo = colorInfo[i]
)
}
}
return(as.data.frame(rbindlist(df_list)))
}
# Make plotly rectangle data based on x/y min/max info
.plotly_makeRectData <- function(
xmin, xmax, ymin, ymax, info, colorInfo, resolution = 20
) {
df_list <- list()
if(resolution < 1) resolution <- 1
for(i in seq_len(length(info))) {
df_list[[5 * i - 4]] <- data.table(
x = seq(xmin[i], xmax[i], length.out = resolution + 1)[-1],
y = ymin[i],
text = info[i],
colorInfo = colorInfo[i]
)
df_list[[5 * i - 3]] <- data.table(
x = xmax[i],
y = seq(ymin[i], ymax[i], length.out = resolution + 1)[-1],
text = info[i],
colorInfo = colorInfo[i]
)
df_list[[5 * i - 2]] <- data.table(
x = seq(xmax[i], xmin[i], length.out = resolution + 1)[-1],
y = ymax[i],
text = info[i],
colorInfo = colorInfo[i]
)
df_list[[5 * i - 1]] <- data.table(
x = xmin[i],
y = seq(ymax[i], ymin[i], length.out = resolution + 1),
text = info[i],
colorInfo = colorInfo[i]
)
if(i < length(info)) {
df_list[[5 * i]] <- data.table(
x = NA, y = NA, text = NA, colorInfo = colorInfo[i]
)
}
}
return(as.data.frame(rbindlist(df_list)))
}
# Assemble plotly given seed plotly object and x/y/info data
.pV_assemble_covPlotly2 <- function(p, vLayout) {
numCovTraces <- numDiffTraces <- c()
if(length(p@covTrack) > 0) {
numCovTraces <- vapply(p@covTrack[[1]][["dataList"]], length, 0)
numCovTraces <- (numCovTraces - 2)/2
}
if(length(p@diffTrack) > 0) {
numDiffTraces <- vapply(p@diffTrack[[1]][["dataList"]], length, 0)
}
p <- .knitPlotly(p, numCovTraces, numDiffTraces, vLayout)
# inject plot data
for(i in seq_len(length(numCovTraces))) {
p <- .injectPlotData(p, p@args[["covTrackPos"]][i],
p@covTrack[[1]][["dataList"]][[i]],
p@covTrack[[1]][["layoutList"]][[i]][["title"]]
)
}
for(i in seq_len(length(numDiffTraces))) {
p <- .injectPlotData(p, p@args[["diffTrackPos"]][i],
p@diffTrack[[1]][["dataList"]][[i]],
p@diffTrack[[1]][["layoutList"]][[i]][["title"]]
)
}
p <- .injectPlotData(p, p@args[["annoTrackPos"]],
p@annoTrack[[1]][["dataList"]],
p@annoTrack[[1]][["layoutList"]][["xtitle"]]
)
# Layout modifications
# x axis
if("xtitle" %in% names(p@annoTrack[[1]][["layoutList"]])) {
p <- .adjustXtitle(p, p@annoTrack[[1]][["layoutList"]][["xtitle"]])
}
if("xrange" %in% names(p@annoTrack[[1]][["layoutList"]])) {
p <- .adjustXrange(p, p@annoTrack[[1]][["layoutList"]][["xrange"]])
}
# y axis
axisNum <- 0
for(i in seq_len(length(numCovTraces))) {
axisNum <- axisNum + 1
p <- .fixYrange(p, axisNum)
if("yrange" %in% names(p@covTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYrange(p,
p@covTrack[[1]][["layoutList"]][[i]][["yrange"]],
axisNum
)
}
if("title" %in% names(p@covTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYtitle(p,
p@covTrack[[1]][["layoutList"]][[i]][["title"]],
axisNum
)
}
}
for(i in seq_len(length(numDiffTraces))) {
axisNum <- axisNum + 1
p <- .fixYrange(p, axisNum)
if("yrange" %in% names(p@diffTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYrange(p,
p@diffTrack[[1]][["layoutList"]][[i]][["yrange"]],
axisNum
)
}
if("title" %in% names(p@diffTrack[[1]][["layoutList"]][[i]])) {
p <- .adjustYtitle(p,
p@diffTrack[[1]][["layoutList"]][[i]][["title"]],
axisNum
)
}
}
axisNum <- axisNum + 1
p <- .adjustYtitle(p, "", axisNum)
if("yrange" %in% names(p@covTrack[["layoutList"]][[i]])) {
p <- .adjustYrange(p, p@covTrack[["layoutList"]][[i]][["yrange"]],
axisNum)
}
return(p)
}
# Assemble ggplot given sub-track data
.pV_assemble_ggplot <- function(
covTrack, diffTrack, annoSubTrack, annoFullTrack,
verticalLayout, horizontalLayout
) {
# Correct vertical layout for multiple cov/diff tracks
nCovV <- nDiffV <- 0
if(length(covTrack) > 0) nCovV <- length(covTrack[[1]])
if(length(diffTrack) > 0) nDiffV <- length(diffTrack[[1]])
realVL <- c()
Vcursor <- 0
if(nCovV > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL,
rep(
verticalLayout[Vcursor] / nCovV,
nCovV
)
)
}
if(nDiffV > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL,
rep(
verticalLayout[Vcursor] / nDiffV,
nDiffV
)
)
}
if(length(annoSubTrack) > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL, verticalLayout[Vcursor])
}
if(length(annoFullTrack) > 0) {
Vcursor <- Vcursor + 1
realVL <- c(realVL, verticalLayout[Vcursor])
}
realVL <- round(realVL * 1000) # round to nearest integer
realHL <- horizontalLayout
hCount <- length(realHL)
vCount <- length(realVL)
bVec <- seq_len(length(realVL))
rVec <- seq_len(length(realHL))
if(vCount > 1){
tVec <- c(1, 1 + bVec[-vCount])
} else {
tVec <- 1
}
if(hCount > 1){
lVec <- c(1, 1 + rVec[-hCount])
} else {
lVec <- 1
}
# Set up patchwork
hasAnnoFT <- length(annoFullTrack) == 1
for(v in seq_len(length(realVL))) {
if(v == vCount & hasAnnoFT) {
layout_instr <- c(layout_instr, patchwork::area(
t = tVec[v], b = bVec[v], l = 1, r = hCount
))
}
for(h in seq_len(length(realHL))) {
if(v == 1 & h == 1) {
# initialize
layout_instr <- patchwork::area(
t = tVec[1], b = bVec[1], l = lVec[1], r = rVec[1]
)
} else {
layout_instr <- c(layout_instr, patchwork::area(
t = tVec[v], b = bVec[v],
l = lVec[h], r = rVec[h]
))
}
}
}
fullPlotList <- list()
plotCount <- 0
rowCount <- 0
if(length(covTrack) > 0) {
for(j in seq_len(nCovV)) {
rowCount <- rowCount + 1
for(i in seq_len(length(covTrack))) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- .pV_mod_ggplot(
covTrack[[i]][[j]], i, rowCount, hCount, length(realVL)
)
}
}
}
if(length(diffTrack) > 0) {
for(j in seq_len(nDiffV)) {
rowCount <- rowCount + 1
for(i in seq_len(length(diffTrack))) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- .pV_mod_ggplot(
diffTrack[[i]][[j]], i, rowCount, hCount, length(realVL)
)
}
}
}
if(length(annoSubTrack) > 0) {
rowCount <- rowCount + 1
for(i in seq_len(length(annoSubTrack))) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- .pV_mod_ggplot(
annoSubTrack[[i]], i, rowCount, hCount, length(realVL)
)
}
}
if(hasAnnoFT) {
plotCount <- plotCount + 1
fullPlotList[[plotCount]] <- annoFullTrack[[1]]
}
return(
patchwork::wrap_plots(fullPlotList) +
patchwork::plot_layout(
design = layout_instr,
widths = realHL,
heights = realVL
)
)
}
# Modify gg subplots given their position
# - removes y axis on non-left plots
# - removes legend on non-right plots
# - removes x axis on non-bottom plots
.pV_mod_ggplot <- function(p, x, y, maxX, maxY) {
theme_nonright <- theme(
legend.position = "none"
)
theme_nonleft <- theme(
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
strip.background = element_blank(),
strip.text.x = element_blank()
)
theme_nonbottom <- theme(
axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.ticks.x = element_blank()
)
if(x > 1) p <- p + theme_nonleft
if(x < maxX) p <- p + theme_nonright
if(y < maxY) p <- p + theme_nonbottom
return(p)
}
################################################################################
.gpo_highlight_anno <- function(
reduced.DT, highlight_gr
) {
introns <- reduced.DT[get("type") == "intron"]
introns[, c("highlight") := "0"]
exons <- reduced.DT[get("type") == "exon"]
exons[, c("highlight") := "0"]
misc <- reduced.DT[get("type") == "CDS"]
misc[, c("highlight") := "0"]
tr_filter <- c()
if (length(highlight_gr) == 1) {
# IR / RI only
gr <- highlight_gr[[1]]
introns.gr <- .grDT(introns)
OL <- findOverlaps(gr, introns.gr)
introns[OL@to, c("highlight") := "3"] # purple
OL2 <- findOverlaps(gr, introns.gr, type = "equal")
introns[OL2@to, c("highlight") := "2"] # red
# exons that overlap introns
exons.gr <- .grDT(exons)
OL3 <- findOverlaps(gr, exons.gr, type = "within")
exons[OL3@to, c("highlight") := "1"] # blue
} else if (length(highlight_gr) == 2) {
AS_count <- 1
for (event in highlight_gr) {
tr <- list()
for(i in seq_len(length(event))) {
OL <- findOverlaps(event[i], .grDT(introns), type = "equal")
tr[[i]] <- introns[OL@to]$transcript_id
}
tr_final <- NULL
if(length(event) == 2) {
tr_final <- intersect(tr[[1]], tr[[2]])
} else {
tr_final <- tr[[1]]
}
if(length(tr_final) > 0) {
# Highlight introns that match exact junction
gr <- event
introns.gr <- .grDT(introns)
OL <- findOverlaps(gr, introns.gr, type = "equal")
introns[
seq_len(nrow(introns)) %in% OL@to &
get("transcript_id") %in% tr_final,
c("highlight") := as.character(AS_count)
]
# Remove highlight from novel transcripts
# if not all introns highlighted
introns_novel <- introns[
grepl("novel", get("transcript_id")) &
get("highlight") == as.character(AS_count)
]
introns_novel_noHL <- introns[
get("transcript_id") %in% introns_novel$transcript_id &
get("highlight") == "0"
]
if(nrow(introns_novel_noHL) > 0) {
introns[
get("transcript_id") %in% introns_novel_noHL$transcript_id,
c("highlight") := "0"
]
tr_final <- setdiff(tr_final,
introns_novel_noHL$transcript_id)
}
coord_keys_start <- end(gr[1]) + 1
coord_keys_end <- start(gr[1]) - 1
if (length(gr) == 2) {
coord_keys_start <- c(coord_keys_start, end(gr[2]) + 1)
coord_keys_end <- c(coord_keys_end, start(gr[2]) - 1)
}
exons[get("transcript_id") %in% tr_final &
(get("start") %in% coord_keys_start |
get("end") %in% coord_keys_end),
c("highlight") := as.character(AS_count)]
}
AS_count <- AS_count + 1
}
}
for(hl in c("1", "2")) {
exons_selected <- exons[get("highlight") == hl]
hl_trid <- unique(exons_selected$transcript_id)
OL <- findOverlaps(
.grDT(misc),
.grDT(exons_selected)
)
misc[seq_len(nrow(misc)) %in% unique(OL@from) &
get("transcript_id") %in% hl_trid,
c("highlight") := hl]
}
return(rbind(introns, exons, misc))
}
.gpo_filter_anno <- function(
reduced.DT, transcripts.DT,
selected_transcripts = "",
plot_key_isoforms = FALSE
) {
# Remove by named transcripts only
if(is_valid(selected_transcripts)) {
if (
all(selected_transcripts %in%
c(
transcripts.DT$transcript_id,
transcripts.DT$transcript_name
)
)
) {
transcripts.DT <- transcripts.DT[
get("transcript_name") %in% selected_transcripts |
get("transcript_id") %in% selected_transcripts
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
return(reduced.DT)
}
}
if(plot_key_isoforms) {
reduced.DT.HL <- reduced.DT[get("highlight") != "0"]
reduced.DT.nonHL <- reduced.DT[get("highlight") == "0"]
reduced.DT.nonHL <- reduced.DT.nonHL[
grepl("novel", get("transcript_id"))
]
if(nrow(reduced.DT.HL) != 0) {
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% reduced.DT.HL$transcript_id &
# remove non-highlighted novel elements
!(get("transcript_id") %in% reduced.DT.nonHL$transcript_id)
]
reduced.DT <- reduced.DT[
get("transcript_id") %in% transcripts.DT$transcript_id
]
return(reduced.DT)
}
}
# Otherwise - do nothing
return(reduced.DT)
}
# Allows re-stack. Useful for filtered DTlist
.gcd_stack_anno <- function(
DTlist, view_start, view_end,
reverseGenomeCoords = FALSE,
condensed = FALSE
) {
transcripts.DT <- copy(DTlist$transcripts.DT)
reduced.DT <- copy(DTlist$reduced.DT)
condense_this <- condensed
# Filter before condense
reduced.DT <- reduced.DT[get("end") > view_start & get("start") < view_end]
transcripts.DT <- transcripts.DT[
get("transcript_id") %in% reduced.DT$transcript_id]
if (condense_this != TRUE & nrow(transcripts.DT) <= 100) {
condense_this <- FALSE
transcripts.DT[, c("group_id") := get("transcript_id")]
reduced.DT[, c("group_id") := get("transcript_id")]
} else {
condense_this <- TRUE
transcripts.DT[, c("group_id") := get("gene_id")]
reduced.DT[transcripts.DT, on = "transcript_id",
c("group_id") := get("gene_id")]
}
group.grl <- split(.grDT(transcripts.DT), transcripts.DT$group_id)
group.DT <- as.data.table(range(group.grl))
group.DT$group <- NULL
data.table::setnames(group.DT, "group_name", "group_id")
# Filter before stack
# group.DT <- group.DT[get("end") > view_start & get("start") < view_end]
# apply plot_order on transcripts.DT
OL <- findOverlaps(.grDT(group.DT), .grDT(group.DT), ignore.strand = TRUE)
group.DT$plot_level <- 1
cur_level <- 1
while (any(group.DT$plot_level == cur_level)) {
j <- match(cur_level, group.DT$plot_level)
repeat {
bump_up_trs <- unique(OL@to[OL@from == j])
bump_up_trs <- bump_up_trs[bump_up_trs > j]
bump_up_trs <- bump_up_trs[
group.DT$plot_level[bump_up_trs] == cur_level]
if (length(bump_up_trs) > 0)
group.DT[bump_up_trs, c("plot_level") := cur_level + 1]
j <- j + match(cur_level, group.DT$plot_level[-seq_len(j)])
if (is.na(j)) break
}
cur_level <- cur_level + 1
}
if (condense_this == TRUE) {
group.DT[transcripts.DT, on = "group_id",
c("group_name", "group_biotype") :=
list(get("i.gene_name"), get("i.gene_biotype"))]
} else {
group.DT[transcripts.DT, on = "group_id",
c("group_name", "group_biotype") :=
list(get("i.transcript_name"), get("i.transcript_biotype"))]
}
reduced.DT$group_id <- factor(reduced.DT$group_id,
unique(group.DT$group_id), ordered = TRUE)
reduced.DT[group.DT, on = "group_id",
c("plot_level") := get("i.plot_level")]
setorderv(reduced.DT, "highlight")
return(list(
group.DT = group.DT,
reduced.DT = reduced.DT,
condense_this = condense_this
))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.