Motivation

spliceTerrain was designed to fit into R/Bioconductor workflows. Candidate regions, exon models, and features of interest are often already represented as GRanges or GRangesList objects from GenomicRanges (Lawrence et al. 2013). spliceTerrain allows reuse of these objects for regions, transcript annotation, and highlights alongside indexed BAM files, within the same R session. For example, exon models retrieved from an EnsDb object can become annotation rows without first exporting them to an annotation file.

spliceTerrain brings these inputs together in a focused, sample-resolved sashimi workflow. It returns a patchwork plot or an inspectable genomic context. Its genomic summaries and plotting settings can be edited and replotted.

spliceTerrain displays one panel per BAM for visual inspection of local splicing patterns. It does not perform differential-splicing testing, estimate event-level PSI, or infer transcript isoforms.

Related tools

Several established tools support related visualisation tasks:

• Gviz provides genomic tracks, including BAM-backed coverage and sashimi plots (Hahne and Ivanek 2016).
• ggbio provides genomic visualisation based on ggplot2, with support for Bioconductor objects and transcript-based splice summaries (Yin et al. 2012).
• SpliceWiz provides coverage and sashimi plots within an alternative-splicing analysis workflow. Its event-normalised coverage plots offer an alternative view of splicing differences, summarising replicates by condition rather than relying on individual-sample sashimi plots (Wong et al. 2024).
• ggsashimi provides a command-line sashimi plotting workflow with sample overlays and aggregation (Garrido-Martín et al. 2018).

spliceTerrain focuses on sample-resolved sashimi plots within R/Bioconductor workflows. It is particularly useful when regions and annotations are already represented as genomic objects and plots need iterative refinement. It combines:

• Editable genomic summaries: inspect or modify coverage, junctions, and settings in a returned context, then replot without importing the BAM files again.
• Coverage-aware compression: preserve selected blocks defined by coverage, annotation, and junction anchors while shortening gaps. Coverage retained after filtering can preserve regions within annotated introns at their original width.
• Fragment-level support: count overlapping mates once per covered base and once per distinct junction, while retaining mapped singletons.
• Stable local junction fractions: calculate usage against the pre-display junction reference, so hiding arcs does not change the remaining percentages.
• Integrated sample normalisation: apply supplied library sizes and optional normalisation factors consistently to coverage and junction support, while retaining raw fragment counts for filtering.

These choices offer specific contrasts with existing workflows. When compression is enabled, ggsashimi uses junction-derived, length-dependent shrinking, whereas spliceTerrain also uses retained coverage and annotation to define preserved blocks. By default, Gviz and ggsashimi count support at the read level, so both mates can contribute to coverage and junction counts. In contrast, spliceTerrain uses fragment-level counts. spliceTerrain displays transcript annotation alongside observed junction support without assigning that support to transcript models. This complements the transcript-model-based summaries provided by ggbio’s plotSpliceSum().

spliceTerrain exposes sample normalisation through dedicated arguments. Gviz supports custom transformations that can implement similar scaling, whereas the ggsashimi command-line interface has no equivalent library-size normalisation option.

Versions compared: Gviz 1.57.4, ggbio 1.61.0, SpliceWiz 1.15.0, and ggsashimi 1.1.5.

Input BAM files and regions

The bam argument is a character vector of indexed BAM file paths. If the vector is named, those names are used as sample labels in the plot. If it is unnamed, labels are derived from the BAM file names. BAM files must be indexed for region-restricted import. Index files should be discoverable by Rsamtools using standard same-directory naming, typically file.bam.bai.

bams <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[c(7, 1)]
names(bams) <- c("HNRNPC knockdown", "Control")

The region argument can be a character string, a GRanges, or a GRangesList. Character regions may use comma-separated coordinates, such as "chr14:70,233,810-70,238,690". For Bioconductor workflows, using a GRanges object is often more convenient.

Because the example BAM files use UCSC style chromosome names, we also convert the region to UCSC style before plotting.

ensDb <- EnsDb.Hsapiens.v75
region <- transcripts(ensDb, filter = TxIdFilter("ENST00000394366"))
seqlevelsStyle(region) <- "UCSC"
region
#> GRanges object with 1 range and 6 metadata columns:
#>                   seqnames            ranges strand |           tx_id     tx_biotype tx_cds_seq_start tx_cds_seq_end         gene_id         tx_name
#>                      <Rle>         <IRanges>  <Rle> |     <character>    <character>        <integer>      <integer>     <character>     <character>
#>   ENST00000394366    chr14 70233810-70238690      + | ENST00000394366 protein_coding         70234874       70238178 ENSG00000100650 ENST00000394366
#>   -------
#>   seqinfo: 1 sequence from hg19 genome

If multiple ranges are supplied, spliceTerrain() reduces them to a single span on one chromosome and uses that span as the BAM query and plotting window. The span is the smallest continuous interval containing all supplied ranges: it starts at the lowest start coordinate and ends at the highest end coordinate.

Adding transcript annotation

Transcript models can be displayed with the annotation argument. This must be a GRangesList, where each list element represents one transcript or feature group. The names of the list elements are used as y-axis labels in the annotation panel. When annotation is present, the final plot contains one panel per BAM file plus one annotation panel. The relative heights of these panels are controlled with panel_heights, where the final value corresponds to the annotation panel. Setting annotated_junctions = TRUE distinguishes junctions that exactly match an intron in the supplied annotation with solid arcs, while unmatched junctions are drawn with dashed arcs. Junction classification is relative to the supplied annotation, so known junctions may appear unmatched when only a subset of transcripts is provided.

For this vignette, the plotting region covers the SRSF5 gene. We retrieve exon ranges for the gene and keep a small number of transcripts so the annotation panel remains readable.

annotation <- exonsBy(
    ensDb, by = "tx", filter = GeneIdFilter("ENSG00000100650")
)
annotation <- annotation[1:3]
seqlevelsStyle(annotation) <- "UCSC"

spliceTerrain(
    bam = bams,
    region = region,
    annotation = annotation,
    annotated_junctions = TRUE,
    panel_heights = c(2, 2, 1)
)

Styling annotation tracks

Annotation features can be coloured by a metadata column with anno_fill_by. Use "group" to colour each transcript or feature group separately. Specific fill colours can be supplied with anno_fill_colours.

Text labels can be drawn inside annotation features with anno_label_by. The label colour and size are controlled by anno_label_colour and anno_label_size. The min_arrow argument controls how wide an annotation intron must be before directional arrows are drawn.

spliceTerrain(
    bam = bams,
    region = region,
    annotation = annotation,
    panel_heights = c(2, 2, 1),
    anno_fill_by = "group",
    anno_fill_colours = c(
        ENST00000553369 = "#66C2A5",
        ENST00000394366 = "#8DA0CB",
        ENST00000451983 = "#FC8D62"
    ),
    anno_label_by = "exon_rank",
    anno_label_colour = "black",
    anno_label_size = 2.5,
    min_arrow = 100
)

Highlighting regions of interest

The highlight argument marks one or more genomic intervals with shaded vertical bands. This is useful for drawing attention to a splice site, variant position, alternative exon, or other feature of interest. The highlighted interval can be supplied as a character string or as a GRanges object.

spliceTerrain(
    bam = bams,
    region = region,
    highlight = "chr14:70234854-70234854",
    highlight_colour = "red"
)

A wider interval can be highlighted with a semi-transparent colour.

spliceTerrain(
    bam = bams,
    region = region,
    highlight = "chr14:70234056-70234097",
    highlight_colour = scales::alpha("red", 0.2)
)

Local junction-usage labels

The junction_usage argument selects junctions by a genomic position or range and labels their relative fragment support within each sample. Junctions are selected when either their start or end anchor overlaps the supplied interval(s). Junction anchors use intronic coordinates: the start anchor is the first base of the intron, and the end anchor is the last base of the intron, not the adjacent exonic bases.

Selecting a single shared splice site gives the cleanest interpretation. For junctions sharing a donor or acceptor, the labels describe their relative support at that site. By default, the junction count and usage percentage are separated by a new line. Use junction_usage_label_sep = " " to display them together on one line.

spliceTerrain(
    bam = bams,
    region = region,
    junction_usage = "chr14:70234854-70234854",
    junction_usage_label_sep = " "
)

A range can also be supplied, including one covering multiple splice sites. It defines a combined set of selected junctions; each junction enters the denominator once, even if both anchors or overlapping intervals select it.

spliceTerrain(
    bam = bams,
    region = region,
    junction_usage = "chr14:70234098-70234854"
)

For each sample, the percentage is the junction’s raw fragment support divided by the summed support of all selected, eligible imported junctions. This is calculated after import, MAPQ, and strand filtering, but before min_junction_reads filtering or plotting-window exclusion. An alternative junction remains in the denominator when its arc is hidden. For example, supports of 90 and 9 give 90.9% and 9.1%; hiding the second arc with a threshold of 10 leaves the first label at 90.9%, not 100%. Visible percentages therefore need not sum to 100%. Library-size normalisation does not change the fractions. The reference is conditional on the imported data: changing the query region, MAPQ threshold, or strand selection can change the available evidence.

cassette <- data.frame(
    Junction = c("Upstream inclusion", "Downstream inclusion", "Skipping"),
    Fragments = c(50L, 50L, 50L)
)
cassette[["Junction usage (%)"]] <- round(
    100 * cassette$Fragments / sum(cassette$Fragments), 1
)
pander(cassette)

Filtering signal

Three arguments control how much signal is imported and displayed. min_mapq filters alignments by mapping quality during BAM import. min_coverage filters positions by per-base fragment coverage before plotting. min_junction_reads sets the minimum fragment support required to display a junction arc. Both mates supporting the same junction contribute one count. Both thresholds are applied to raw fragment counts before library-size normalisation. The latter two arguments also accept percentage strings such as "10%". Percentage thresholds are calculated separately for each BAM relative to its maximum raw coverage or junction count within the plotting region. When supplying one threshold per BAM, whole-number counts and percentages can be mixed in a character vector.

spliceTerrain(
    bam = bams,
    region = region,
    min_mapq = 10,
    min_coverage = 10,
    min_junction_reads = "10%"
)

Stranded RNA-seq data

The strandedness argument describes how read strand should be interpreted. It can be one of "unstranded", "forward", or "reverse", and can be supplied either as a single value for all BAM files or one value per BAM file. The example BAM files are unstranded, so both read orientations contain signal from the transcript.

spliceTerrain(
    bam = bams,
    region = region,
    strandedness = "unstranded"
)

Specifying a stranded protocol for unstranded BAM files filters reads by strand. As a result, either "reverse" or "forward" shows only part of the coverage that is visible when strandedness = "unstranded".

spliceTerrain(
    bam = bams,
    region = region,
    strandedness = "forward"
)

spliceTerrain(
    bam = bams,
    region = region,
    strandedness = "reverse"
)

For paired-end BAM files, spliceTerrain() detects paired-end status from the BAM alignment flags and imports alignments accordingly. Mapped reads are retained even when their mates are unmapped. Read-one/read-two flags determine strand interpretation for both complete pairs and singletons. Unmapped, secondary, supplementary, and below-min_mapq alignments are excluded. If region specifies a strand and strandedness is not "unstranded", alignments are restricted to the resolved region strand.

Coverage and junction support are fragment-level summaries. Each imported mate group, mapped singleton, or single-end read is one fragment. Coverage uses the union of each fragment’s aligned blocks. Overlapping mates contribute at most one count per covered position. Introns and unsequenced gaps between mates remain uncovered. Each fragment supports a distinct junction at most once, even if both mates cross it.

Normalising between samples

When samples have different sequencing depths, coverage and junction counts can be normalised with lib_size. Supply full-library sizes in fragment units before restriction to the plotting region: one per pair or mapped singleton, or one per read for single-end data. Additional normalisation factors, such as trimmed mean of M-values (TMM) factors (Robinson and Oshlack 2010) estimated with edgeR, can be supplied with norm_factors. The effective library size is calculated as lib_size * norm_factors.

By default, counts are normalised to the median effective library size. Use normalise_to to set the target explicitly. The common_y argument is useful when comparing samples because it puts sample panels on the same y-axis scale.

spliceTerrain(
    bam = bams,
    region = region,
    lib_size = c(40e6, 55e6),
    norm_factors = c(0.95, 1.08),
    normalise_to = 50e6,
    common_y = TRUE
)

Compressed genomic layouts

By default, compress_introns = TRUE compacts large uninformative gaps so the plot focuses on observed signal. The intron_width argument controls the approximate width of compacted introns in plot space. The map is built jointly from retained coverage across samples, annotations, and junction/selection anchors. Retained blocks keep their widths; gaps between them receive the chosen width. Consequently, retained intronic coverage or an unannotated covered exon remains part of the displayed signal rather than being compressed simply because it lies inside an annotated intron. Changing coverage filters can change the retained blocks and their map.

spliceTerrain(
    bam = bams,
    region = region,
    compress_introns = TRUE,
    intron_width = 500
)

If exact genomic distances are more important than compact display, intron compression can be disabled.

spliceTerrain(
    bam = bams,
    region = region,
    compress_introns = FALSE
)

Junction arc display

Junction arcs can be tuned with arc_height, arc_side, and arc_scale. Increasing arc_height gives stacked arcs more vertical separation. The arc_side argument controls whether arcs are drawn above, below, or on both sides of the coverage track. With arc_scale = TRUE, arc width scales with filtered fragment support. Junction label size is controlled with junc_text_size.

spliceTerrain(
    bam = bams,
    region = region,
    junction_usage = "chr14:70234854-70234854",
    arc_height = 2,
    arc_side = "above",
    arc_scale = TRUE,
    junc_text_size = 2.5
)

Sample colours and panel layout

The colours argument controls the colours used for sample coverage and junction arcs. It can be a single colour used for all samples or one colour per BAM file.

Axis text sizes can be adjusted with axis_title_size and axis_text_size.

spliceTerrain(
    bam = bams,
    region = region,
    annotation = annotation,
    colours = c("darkblue", "darkred"),
    panel_heights = c(2, 2, 1),
    axis_title_size = 20,
    axis_text_size = 12
)

Returning and reusing plot data

Most users can work directly with the final plot returned by spliceTerrain(). For advanced workflows, set return_ctx = TRUE to return the processed context instead. This contains validated settings, imported alignments, summarised coverage, junction counts, annotation, and overlays. The gal field contains one alignment object per sample: GAlignments for single-end BAMs or GAlignmentsList for paired-end BAMs. The paired-end list groups imported mates and includes mapped singleton reads. The cov and juncs fields contain fragment summaries: coverage_raw holds unnormalised fragment counts and coverage holds the plotting values after any library-size normalisation. The juncs_all field retains eligible junction summaries on the plotting chromosome before display filtering or clipping, including out-of-window junctions from imported alignments. It provides a stable reference for the usage_fraction column in juncs. That column is on the 0–1 scale, with NA for unselected junctions or an undefined percentage. All coordinate-bearing fields use genomic coordinates; compressed plot coordinates are generated internally only when the context is plotted.

The context also supports programmatic inspection, including AI-assisted workflows. Scripts or AI agents working within an R session can inspect coverage and junction summaries, check the values underlying plotted labels, and modify annotations or plotting settings before redrawing the figure. This provides access to the data behind the visualisation without requiring those values to be inferred from the plot image.

ctx <- spliceTerrain(
    bam = bams,
    region = region,
    annotation = annotation,
    min_junction_reads = 1,
    panel_heights = c(2, 2, 1),
    min_arrow = 100,
    return_ctx = TRUE
)

names(ctx)
#>  [1] "bam"                      "region"                   "annotation"               "junction_usage"           "junction_usage_label_sep" "highlight"                "strandedness"             "min_mapq"                 "min_coverage"             "min_junction_reads"       "lib_size"                 "norm_factors"             "normalise_to"             "compress_introns"         "intron_width"             "min_arrow"                "common_y"                 "arc_height"               "arc_side"                 "arc_scale"                "colours"                  "highlight_colour"         "anno_fill_by"             "anno_fill_colours"        "anno_label_by"            "anno_label_colour"        "anno_label_size"          "junc_text_size"           "panel_heights"            "axis_title_size"          "axis_text_size"           "annotated_junctions"      "gal"                      "cov"                      "juncs"                    "juncs_all"

The context can be passed back with the ctx argument to redraw the plot without re-importing BAM files or recalculating coverage and junctions. For example, the annotation can be used to retain only coverage overlapping annotated exons and junctions that exactly match annotated introns.

ctx$cov <- ctx$cov[overlapsAny(ctx$cov, ctx$annotation)]
ctx$juncs <- ctx$juncs[ctx$juncs$annotation_match]
spliceTerrain(ctx = ctx)

The annotation_match field records whether each junction exactly matches an intron implied by the supplied annotation. Plotting settings can also be changed directly, for example with ctx$colours <- c("darkblue", "darkred"). Settings that control BAM import, filtering, or normalisation are not reapplied when a context is replotted because its coverage and junctions are already processed.

Filtering ctx$juncs changes which arcs are displayed without changing the usage denominator. Changing ctx$junction_usage to another GRanges selection recalculates the fractions from ctx$juncs_all when plotting. If you intend to change the reference evidence itself, edit ctx$juncs_all, retaining one row per sample/junction.