Skip to contents

Prepare Sashimi plot data

Usage

prepareSashimi(
  sjenv = NULL,
  flatExonsByGene = NULL,
  filesDF = NULL,
  gene,
  sample_id = NULL,
  minJunctionScore = 10,
  gapWidth = getOption("splicejam.gap", 200),
  addGaps = TRUE,
  baseline = 0,
  compressGR = TRUE,
  compress_introns = TRUE,
  ref2c = NULL,
  gap_feature_type = "intron",
  default_feature_type = "exon",
  feature_type_colname = "feature_type",
  exon_label_type = c("none", "repel", "mark"),
  junc_label_type = c("repel", "mark", "none"),
  return_data = c("df", "ref2c"),
  include_strand = c("both", "+", "-"),
  junc_color = jamba::alpha2col("goldenrod3", 0.7),
  junc_fill = jamba::alpha2col("goldenrod1", 0.4),
  doStackJunctions = TRUE,
  coord_method = c("coord", "scale", "none"),
  scoreFactor = 1,
  scoreArcFactor = 0.2,
  scoreArcMinimum = 100,
  covGR = NULL,
  juncGR = NULL,
  use_memoise = FALSE,
  memoise_coverage_path = "coverage_memoise",
  memoise_junction_path = "junctions_memoise",
  do_shiny_progress = NULL,
  verbose = FALSE,
  ...
)

Arguments

sjenv

environment (or list) containing argument values indicated below. This input is intended to be compatible with output from sashimiDataConstants(). It currently recognizes the following, although function arguments take priority when provided.

  • 'flatExonsByGene'

  • 'ref2c' (optional)

  • 'filesDF'

flatExonsByGene

GRangesList named by gene, specifically named by 'gene_name' as in the 'tx2geneDF' data.frame. GRanges elements are disjoint, non-overlapping genomic ranges per gene and should represent only exon regions. All non-exon regions are assumed to be introns or gaps.

filesDF

data.frame with columns:

  • url: any valid path compatible with data.table::fread()

  • sample_id: character string representing a biological sample, and often may be used for sample grouping, to combine coverage and/or junction read counts for all samples in the group.

  • type: character string with 'bw' for bigwig coverage, 'junction' for BED12 splice junctions, or 'coverageGR' when coverage data are supplied as GRanges objects named by 'sample_id', using argument 'covGR'.

gene

character string of the gene to prepare, which must be present in names(flatExonsByGene).

sample_id

character default NULL, used to subset entries in 'filesDF' to include only specific matching sample_id entries. When NULL, it uses all entries.

minJunctionScore

numeric default 10, to require junctions to have at least this total score per sample_id to be displayed. This filter it useful to hide spurious junction counts which may only have very few supporting read counts, and is typically adjusted relative to the total coverage of the gene. For example, 'GAPDH' gene often has extremely high coverage, and may be associated with numerous spurious junctions with 10 to 50 reads, which are not useful to display when coverage is in the 10,000 order of magnitude.

gapWidth

integer value, default 200 from options options('splicejam.gap'), to set a fixed gap width, or when NULL the gap width is defined as 3 times the median feature width.

addGaps

logical default TRUE, whether to include gap regions in the coverage plot, for example including introns or intergenic regions. When compressGR=TRUE then gaps regions are down-sampled using running maximum signal with roughly the same x-axis resolution as uncompressed regions.

baseline

numeric vector named by names(flatExonsByGene) where baseline is used to adjust the y-axis baseline position above or below zero. Default zero 0 for all practical purposes.

compressGR

logical default TRUE, whether to compress GRanges coordinates in the output data, where gaps/introns are set to a fixed width. When ref2c is not supplied, and compressGR=TRUE, then ref2c is created using make_ref2compressed().

compress_introns

logicaldefault TRUE, whether to compress the coverage polygon coordinates to approximately the same number of pixels per inch as the exon polygons. This option greatly reduces the size of the polygon, since introns are already about 50 to 100 times wider than exons, and when compressGR is TRUE, the introns are visibly compressed to a fixed width on the x-axis. The data has many more x-axis coordinates than the data visualization, this argument is intended to reduce the intron coordinates accordingly.

ref2c

list object output from make_ref2compressed() used to compress axis coordinates, to compress polygon coverage data in compressed regions, and to adjust splice junction arcs using compressed coordinates.

gap_feature_type

character string, default 'intron', with the feature_type representing gaps when addGaps=TRUE.

default_feature_type

character string, default 'exon' with the default feature_type representing non-gap regions.

feature_type_colname

character string with the column name of values() to represent the feature type (exon/intron).

exon_label_type

character string indicating the type of label for exons:

  • 'none' (default): no exon label is displayed

  • 'repel': use ggrepel to display the exon label

  • 'mark': use ggforce::geom_mark_rect() to display exon label

junc_label_type

character string indicating the type of label for junction counts:

  • 'repel' (default): use ggrepel to display the exon label

  • 'mark': use ggforce::geom_mark_rect() to display exon label

  • 'none': no exon label is displayed

return_data

character string, default 'df', with format of data to return: 'df' is data.frame and 'ref2c' returns only the 'ref2c' axis compression as list. Typically only 'df' is used.

include_strand

character string, default 'both' with strandedness of data to include. Coverage strandedness is determined by filename, assuming some substring which matches 'plus|pos|+' represents positive strand. By convention in splicejam, scores are negative for negative strand coverage, and negative strand junction counts, however not required.

  • "both": include coverage and junctions from both strands,

  • "+": include coverage and junctions from only the '+' strand,

  • "-": include coverage and junctions from only the '-' strand.

junc_color, junc_fill

character valid R colors used for junction ribbon arcs. The alpha transparency is maintained, and when junctions are stacked, the actual color is slightly adjusted to vary from light to dark in order of stacked junction read counts. Default uses 'goldenrod3' border, and 'goldenrod1' for the fill color.

doStackJunctions

logical default TRUE, whether to stack junction arcs at each end, this argument is passed to grl2df() which calls stackJunctions(). The purpose of stacking junctions is to place each ribbon junction arc starting above the previous ribbon arc, for each junction start, and junction end. In this way, the total height of junction ribbon widths should represent the total junction read count/score at each junction boundary.

coord_method

character string indicating how the x-axis coordinates are represented for use in ggplot2 plot functions.

  • 'coord' (default): stores coordinates unchanged, without regard to intron/gap compression. It allows the axis to be displayed using ggplot2 techniques, optionally using 'ref2c' to compress introns/gaps visually, keeping numeric values as-is.

  • 'scale': stores numeric coordinates after compression when 'compressGR=TRUE', which means these values are not genome coordinates, instead are display coordinates. The ggplot2 axis labels are expected to be adjusted using 'ref2c' for display purposes.

  • 'none': genome coordinates are stored as-is with no anticipated modifications by 'ref2c'. This approach is effectively equivalent to 'coord'.

scoreFactor

numeric default 1, scalar multiplied by each score, expected to be provided as named numeric vector named by 'sample_id'.

scoreArcFactor

numeric default 0.2, used to adjust the magnitude of splice ribbon arc curvature. Use higher values to increase the height of the midpoint of the ribbon arc relative to each end of the ribbon. This value is a scalar, roughly applied as '1 + scoreArcFactor' such that 0.2 has the effect of 'score * 1.2', a 20% addition.

scoreArcMinimum

numeric default 100, the minimum height of the midpoint of the junction ribbon arc above the flat mean value between the start and end of the ribbon. It is used as minimum, when the scoreArcFactor provide too low a value. Useful when there may be intermediate coverage, or intervening exon, to have the arc curve above the coverage.

covGR

GRanges default NULL, optional coverage data in columns stored as NumericList, where colnames(GenomicRanges::values(covGR)) are present in 'filesDF' for rows with filesDF$type %in% "coverage_gr". The name should match 'url' or 'sample_id' in 'filesDF'.

juncGR

GRanges default NULL, optional splice junctions, 'score' is used for the abundance of splice junction reads, and 'sample_id' is used to define the biological sample as defined in 'filesDF'.

use_memoise

logical default FALSE, passed to getGRcoverageFromBw() and import_juncs_from_bed() as needed, indicating whether to use memoise to cache intermediate data results. Default is TRUE for launchSashimiApp() and FALSE otherwise. In general, TRUE would be a valid default option.

memoise_coverage_path, memoise_junction_path

character string with default folder path to store memoise cache files, used with use_memoise=TRUE.

do_shiny_progress

default NULL uses getOption("splicejam.progress", FALSE) which can be FALSE to show no progress, or a function which updates an existing progress bar, for example the output of progressr::progressor() is such a function.

  • The default is NULL to avoid having this argument be conditional when also wrapping prepareSashimi() as a memoise cached function. This argument remains NULL with or without enabling a progress bar, by setting the option to control the value.

verbose

logical whether to print verbose output.

...

additional arguments are passed to make_ref2compressed(), getGRcoverageFromBw(), exoncov2polygon().

Value

list containing ggSashimi a ggplot2 graphical object containing a full Sashimi plot; ggCov the RNA-seq coverage subset of the Sashimi plot; ggJunc the splice junction subsset of the Sashimi plot; ref2c the output of make_ref2compressed() used for ggplot2 coordinate visualization; covDF, juncDF data.frame objects with the raw data used to create ggplot2 objects; covGR, juncGR the GRanges objects used to create the data.frames; gr the GRanges object representing the exons for the gene of interest; juncLabelDF the data.frame containing exon label coordinates used to add labels to the splice junction arcs.

Details

This function is the workhorse function used to produce Sashimi plots, and is intended to be a convenient wrapper function for several other individual functions.

At a minimum, a Sashimi plot requires three things:

  1. Exons, usually from a gene of interest.

  2. RNA-seq coverage data.

  3. Splice junction data.

There is some required pre-processing before running prepareSashimi():

  • Prepare flattened exons by gene using flattenExonsByGene() and corresponding data, including exonsByGene, cdsByGene, and tx2geneDF. Verify the gene exon model data using gene2gg().

  • Find file paths, or web URLs, for a set of bigWig coverage files, representing RNA-seq coverage for each strand, for the samples of interest. Test the coverage data using getGRcoverageFromBw() for a small set of GRanges data.

  • Find file paths, or web URLs, for a set of BED6 or BED12 format files, note that it cannot currently use bigBed format due to limitations in the rtracklayer package. Test the splice junction data using rtracklayer::import() for a small range of GRanges features, then send the data to spliceGR2junctionDF() to prepare a data.frame summary.

The basic input for coverage and junction data is a data.frame, which defines each file path or url, the type of data "bw" or "junction", and the biological sample "sample_id". Any file path compatible with rtracklayer::import() will work, including web URLs and local files. When using a web URL you may need to use "https://" format to force the use of secure web requests, but this requirement varies by country.

See also

Examples

# The active example below uses sample data
suppressPackageStartupMessages(library(GenomicRanges));

data(test_exon_gr);
data(test_junc_gr);
data(test_cov_gr);
filesDF <- data.frame(url="sample_A",
   type="coverage_gr",
   sample_id="sample_A");
sh1 <- prepareSashimi(
   flatExonsByGene=GRangesList(TestGene1=test_exon_gr),
   filesDF=filesDF,
   gene="TestGene1",
   covGR=test_cov_gr,
   juncGR=test_junc_gr);
plotSashimi(sh1);