An alternative splicing event describes local RNA processing. It usually does not identify the complete transcripts expressed in either condition. The aim of matching is to select plausible annotated transcripts for the two forms, then compare the proteins they encode. A transcript that only overlaps a small part of the event is not sufficient evidence for that form.
get_ranked_pairs() is the default matcher in get_splicing_impact(); the
legacy independent matcher remains available with matching = "legacy". It
first checks AS-form compatibility, then ranks pairs of compatible
transcripts. Coding status (an annotated CDS) and transcript support remain
priorities. Within a support tier, similarity of coding sequence outside the
event precedes exact TSL.
This protocol ranks plausible representatives. It does not estimate transcript abundance, phase distant splicing events, establish actual isoform expression, or assign probabilities that a particular pair is expressed. Protein features and interaction consequences remain predictions conditional on the selected pair and the supplied reference.
Supply canonical differential rows with event_id, event_type, form,
gene_id, chr, strand, inc, and exc. Coordinates must be one-based,
closed; 101-200 contains 100 bases. Multiple intervals use semicolons.
The rMATS import functions perform their zero-based-start conversion before
matching. Do not convert their standardized output a second time.
Each event ID must refer to one gene, chromosome, strand, and event type. Each event/form/interval definition must identify one differential row. Malformed coordinates or conflicting metadata are errors raised before any matching; the message names the affected events. IDs must already agree with the reference; matching does not guess gene aliases or mix releases.
In source = "multi", each positive-delta-PSI form is compared with each
negative-delta-PSI form in the same event. Positive is case; negative is control.
Thus a terminal event with two increasing and three decreasing sites has six
form comparisons, recorded as n_event_comparisons on each selected pair. The
ORF protocol does not replace those comparisons with a single gene-wide
winner. In source = "paired", one INC and one EXC are
required. If delta PSI is absent entirely, INC is assigned to case, preserving
the existing form-based convention. Rows without an opposing form are reported.
Significance is separate from matching. Filter with keep_sig_pairs() before
calling get_ranked_pairs() if that is the desired analysis population. The
pipeline wrapper performs this filtering using its existing thresholds.
Index the relevant transcript metadata and exon structures. A candidate must have the event’s gene, chromosome and strand. Exons are ordered by genomic position for junction comparisons and by transcript direction when extracting coding sequence. Exon IDs never determine biological order.
Conflicting transcript metadata, overlapping exons within a transcript, and conflicting sequence records for one transcript ID stop the calculation. Identical duplicate records are harmless. A failed match remains visible in the event diagnostics; it is not interpreted as evidence for a noncoding or missing biological isoform.
Any overlap between an explicit forbidden interval and a candidate exon rejects the candidate, and every event-defining boundary in the table below must match exactly. There is no 5% eligibility threshold in this protocol.
A candidate that also contains every required interval completely is an
exact match. The outer ends of flanking exons, the shared boundary of an
A5SS/A3SS exon and the TSS/PAS end of a terminal exon are not part of the
event, so a candidate that differs only there is a relaxed match. The
structural_match column of candidates and the selected pairs records
which applied. Relaxed matches must still contain their form’s alternative
region (otherwise the pair is rejected as event_region_not_covered), and they
lose ties to exact matches. For AFE/HFE/ALE/HLE they are considered only when
no transcript matches the terminal exon exactly, so two sites that share a
splice site remain distinguishable.
| Event | Additional eligibility rule |
|---|---|
| SE inclusion | With three intervals, the upstream, included and downstream exons must be adjacent in the transcript and have the specified inner splice boundaries. The flanks’ outer ends may differ (relaxed). |
| SE exclusion | The two flanking exons must be adjacent, with the specified skip junction; the skipped exon must be absent. The flanks’ outer ends may differ (relaxed). |
| MXE | Each three-interval branch must follow its specified adjacent-exon chain and exclude the other branch’s forbidden exon. The flanks’ outer ends may differ (relaxed). |
| RI inclusion | The intron must lie in a single annotated exon that reaches into both flanks; exact matches contain both flanks completely. Separate exons do not establish retention. |
| RI exclusion | The flanks must form the specified adjacent-exon junction, with the intron absent from exons. The flanks’ outer ends may differ (relaxed). |
| A5SS/A3SS | The variable donor/acceptor boundary must match exactly on the appropriate strand, and the variable exon must be adjacent to the partner exon with the specified junction. The exon’s shared boundary and the partner’s outer end may differ (relaxed). |
| AFE/HFE | The transcript’s first exon, in transcript direction, must match the specified interval exactly; without such a transcript, a first exon with the same donor and another TSS is accepted as relaxed. |
| ALE/HLE | The transcript’s last exon, in transcript direction, must match the specified interval exactly; without such a transcript, a last exon with the same acceptor and another PAS is accepted as relaxed. |
| Generic event | Full required-interval coverage and forbidden-interval exclusion; labelled interval_only. |
The rMATS converters define each A3SS/A5SS form by the variable exon and the
flanking exon it is spliced to (rMATS’s flankingES/flankingEE), so the
junction is verified (junction_chain). A definition with the variable exon
alone, as in DI tables saved by earlier versions or some user imports, is
labelled splice_site_only: the variable boundary can be verified, but the
partner cannot. The matcher does not invent missing context.
For generic imports labelled SITE, SE skip and RI retained/spliced branches
can be recognized from their supplied interval topology. Original form labels
are preserved in the output; event grouping still comes from the supplied IDs.
A single internal exon supplied for SE inclusion or MXE is accepted only with
exact exon boundaries and is labelled partial_event; its omitted flanking
junctions are unverified. RI supplied as one continuous retained span is also
labelled partial_event. A single-exon transcript can satisfy a terminal rule
at either end. Hybrid first/last events use their terminal role for the chosen
form; they do not require the alternative terminal exon to be absent from all
internal positions in the other transcript.
These rules can still reduce the number of matched events. Novel junctions, annotation-version differences, or partial event definitions may leave no representative pair. That is an explicit unresolved result. Relaxation never moves an event-defining boundary, and every relaxed match is labelled.
For each form comparison, cross the compatible case and control candidates.
A transcript is never paired with itself. A pair is also rejected if both
transcripts satisfy both form definitions: those candidates cannot distinguish
the requested forms. A pair is also rejected if either transcript lacks part of
its own form’s alternative region (event_region_not_covered). These cases are
recorded in pair_rejections. Nothing is selected independently before
considering its partner.
max_candidates limits each form to 100 candidates per comparison, so at most
10,000 pairs are scored. A form with more keeps its best-supported
candidates, in the ranking’s own order: annotated CDS, TSL tier, exact
structure, numeric TSL, then transcript ID. Dropped candidates are counted in
candidates_dropped_case and candidates_dropped_control, and a message
reports how many comparisons were affected. In the default human GENCODE v45
annotation (TSL 1–3) the largest gene has 58 transcripts, so nothing is
dropped; the limit matters with filter_tsl = NULL, where some genes have
several hundred transcripts.
Reduce the required intervals for each form to genomic ranges. The event mask is their symmetric difference, plus the explicit forbidden intervals from either form. Shared flanking sequence is retained as context.
For SE this normally masks the skipped exon. For MXE it masks both alternative exons. For RI it masks the retained intron. For A3SS/A5SS it masks the extension between the alternative boundaries. For terminal events it masks the regions unique to the alternative terminal intervals.
Positions in this mask are excluded from the context comparison in Step 6, so the score describes how the two transcripts differ outside the event.
Exon-level cds_gen_start and cds_gen_stop define the coding bases. Each pair
is compared within a coding window:
context_window = "pair_cds_span")."event_cds_span"), or the
candidates’ full span when none is coding ("event_exon_span"). A shared
window keeps those pairs comparable with each other.context_similarity is the base-pair Jaccard overlap of the two transcripts’
exonic positions in the window, outside the event mask:
\[ S = \frac{\text{exonic positions shared by both transcripts}} {\text{exonic positions in either transcript}}. \]
Positions in one genome determine sequence, so this is a direct measure of how much coding-region structure the two transcripts share. It deliberately does not depend on where each annotated reading frame ends. Consider an exon skip that shifts the frame:
Comparing annotated CDS sequences would favour the second isoform and hide the
frameshift. Comparing exonic positions scores the plain skip 1 and the rescue
slightly lower, because of its extra change. orf_similarity is the same score
including the event region and is a later tie-breaker. The diagnostics also give
the shared fraction of each transcript’s compared positions and the number of
positions compared.
The sequence table is still checked, because the selected pairs carry sequences
into the downstream protein analyses. transcript_seq may contain CDS or full
RNA sequence:
Length agreement is a consistency check, not proof that independently supplied
FASTA and GTF files describe the same genome. Use a coherent reference bundle.
The evidence label records what supports each pair: cds_sequences (both CDS
sequences usable downstream), cds_annotation (both coding, a sequence
unusable), exon_annotation (a member without CDS), or
no_outside_event_context (nothing to compare outside the event).
Protein sequence identity is not used to rank pairs. A legitimate event-induced frameshift can replace much of the downstream protein. Ranking on protein identity, or on annotated CDS sequence, could preferentially hide that consequence. The existing downstream protein alignment, frame/rescue classification and feature analyses still run on the selected pair.
Sort by the following rules, in this exact order:
cds_sequences, cds_annotation, exon_annotation,
then no_outside_event_context.context_similarity.orf_similarity.n_inexact): exact
outer boundaries win ties.Coding status comes from the annotated CDS, not the biotype, so a transcript
annotated as nonsense-mediated decay counts as coding. When an event shifts the
frame, the transcript that carries it faithfully is often annotated that way,
like the plain skip in Step 6. Requiring the protein_coding biotype would hand
such comparisons to a frame-restoring isoform before similarity is considered.
The biotype is still reported: both matchers add transcript_type_case and
transcript_type_control, and compare_sequence_frame() labels a pair with an
NMD-annotated member "NMD", since that transcript is taken to make no protein.
The 1–3 grouping is a policy choice for this implementation, consistent with the package’s default reference filter. It deliberately allows a TSL-2 or TSL-3 transcript to beat a TSL-1 transcript when its coding background fits substantially better. It does not equate their underlying evidence. TSL describes support for the transcript model; it is not a measurement of expression in the samples under study. See the GENCODE FAQ.
n_tied_best counts ties across all substantive criteria, before transcript-ID
ordering. context_score_margin compares the winner’s context score with the
next pair in the same coding, support and evidence tiers. A zero margin can
still be resolved by full ORF similarity or TSL. A missing margin means no
comparable runner-up or no usable context score. Neither quantity is a p-value
or calibrated confidence estimate.
If no structural pair exists for a comparison, and fallback = TRUE (the
default; matching_fallback in get_splicing_impact()), the comparison is
retried with the legacy matcher’s gates: each required interval must overlap
one exon by at least 5% of its width, and no exon may overlap a forbidden
interval by 5%. Event boundaries are not checked. A pair must still use two
different transcripts and must differ at the event in the expected direction.
event_agreement adds, for each form, how much more of that form’s
event-specific bases its own transcript covers than its partner does. It ranges
up to 2 and must be positive.
Approximate pairs are ranked by event_agreement first, then by the Step 7
criteria. They are labelled matching_tier = "approximate", and a member whose
structure was not verified is labelled structural_match = "approximate". For
example, when the annotation lacks the exact alternative exon, the fallback can
choose a transcript with a slightly different version of it; the frame call is
then a consequence of that version. Report or filter approximate pairs
separately where exact event structure matters. Comparisons with no approximate
pair remain in unmatched, whose fallback column records why.
This runnable usage example uses the bundled reference and splicing inputs. The comparison against the legacy method is a separate review exercise.
ex <- load_example_data(c("sample_frame", "annotation_df", "exon_features", "ppi"))
#> [INFO] Loading bundled test annotation data
annotation_df <- ex$annotation_df
protein_feature_total <- ex$protein_feature_total
exon_features <- ex$exon_features
ppi <- ex$ppi
data <- get_rmats_hit(ex$sample_frame)
#> [INFO] Filtered 0 unannotated ALE/AFE events (465 -> 465)
res <- get_differential_inclusion(
data, parallel_glm = FALSE,
BPPARAM = BiocParallel::SerialParam(), verbose = FALSE
)
res_di <- keep_sig_pairs(res)ranked <- get_ranked_pairs(
res_di, annotation_df$annotations, annotation_df$sequences,
source = "multi", verbose = FALSE
)
names(ranked)
#> [1] "pairs" "matched" "candidates" "rankings"
#> [5] "pair_rejections" "events" "unmatched" "settings"
# ranked$pairs: selected pairs, ready for downstream functions
# ranked$rankings: every scored pair, all ranking components and ties
# ranked$candidates: every same-locus candidate, including rejection reasons
# ranked$pair_rejections: self-pairs, form-indistinguishable combinations and
# approximate combinations without an event difference
# ranked$events: every input form, including forms without matches
# ranked$unmatched: form comparisons without a pair in any tier
# ranked$settings: protocol version and scoring parameters
seq_compare <- compare_sequence_frame(ranked$pairs, annotation_df$annotations)
#> [1] "[INFO] Processing 33 transcript and protein sequence alignments, this may take a little..."
#> [1] "[Processing] Identifying frame shifts and rescues"
#> [1] "[INFO] 1 frameshifts (0 rescues) and 16 non-frameshifts were identified, "
hits_domain <- get_domains(seq_compare, exon_features)ranked$matched holds the selected form rows qualified by matching_group.
The same form can select different transcript representatives against different
opposing forms. Use ranked$pairs directly for consequences; running ordinary
get_pairs() on all those selected form rows would cross them again.
out <- get_splicing_impact(
res = res,
annotation_df = annotation_df,
protein_feature_total = protein_feature_total,
ppi = ppi,
matching = "orf",
debug_steps = TRUE,
verbose = FALSE
)
#> [1] "[INFO] Processing 33 transcript and protein sequence alignments, this may take a little..."
#> [1] "[Processing] Identifying frame shifts and rescues"
#> [1] "[INFO] 1 frameshifts (0 rescues) and 16 non-frameshifts were identified, "
names(out$matching)
#> [1] "candidates" "rankings" "pair_rejections" "events"
#> [5] "unmatched" "settings"In S4 mode, the audit tables are in out@metadata$matching_diagnostics;
get_ranked_pairs(return_class = "S4") stores its results the same way for a
stepwise S4 workflow. matching_pair_id includes the selected transcript IDs
and is preserved as the S4 pair key. Legacy results keep their established key
convention. The wrapper does not load unused reference resources when no
significant forms remain; the matching diagnostics are then empty tables, or
NULL when no annotation_df is given.
| Aspect | Legacy | ORF |
|---|---|---|
| Eligibility | At least 5% of each required segment overlaps an exon by default; forbidden-overlap veto at the same threshold | Exact event boundaries and structure, any forbidden overlap rejected; outer flank ends may differ (relaxed); legacy gates only as a labelled approximate fallback |
| Selection unit | One best transcript independently for each form | One best transcript pair for each form comparison |
| Exon class | A ranking preference | Defining boundaries and junction/terminal roles are eligibility requirements where context is supplied |
| Coding preference | protein_coding biotype, as a tie-breaker after fit |
An annotated CDS, so NMD-annotated transcripts count |
| Transcript support | Breaks the last ties, after fit and coding status | Worst TSL tier first; similarity can distinguish candidates within TSL 1–3; exact TSL breaks later ties |
| Similarity at selection | Not used | Shared exonic positions in the coding window, outside the event first, including it second |
| Missing sequence | Appears downstream after transcript selection | Explicit evidence label during ranking |
| Alternative candidates | Discarded at transcript selection | Retained with eligibility reasons and complete pair scores |
| Ambiguity | Implicit in ordering | Tied-best counts, comparable score margin and stable IDs |
| Cost | Overlap search and per-form sorting | Additional candidate cross-products, bounded explicitly; no sequence alignment |
The ORF protocol should reduce cross-context comparisons and prevent structurally incompatible isoforms from winning through TSL or overlap preferences. It also makes disagreements inspectable. Its main costs are approximate pairs where the annotation lacks a transcript with the exact event boundaries, greater computation, and dependence on a correct reference/sequence pairing. Structural differences in the coding window outside the event, such as a different coding start exon, still lower context similarity; a premature stop alone does not. The evidence tiers can favour a pair with usable sequences over one with missing sequence data.
There are no expression weights, new transcript assembly, ORF prediction, machine-learned probabilities or changes to differential testing in this implementation. Before adopting it as the default, inspect retained and rejected real events across all supported classes, review ties and missing context, and confirm that the TSL-tier policy fits the intended analysis. Event definitions follow the package’s standardized inputs; the upstream file conventions are described in the rMATS documentation.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] BiocParallel_1.47.0 SpliceImpactR_1.1.4 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.6.1 BiocIO_1.23.3
#> [3] bitops_1.1-0 ggplotify_0.1.3
#> [5] tibble_3.3.1 R.oo_1.27.1
#> [7] polyclip_1.10-7 enrichit_0.2.6
#> [9] XML_3.99-0.25 lifecycle_1.0.5
#> [11] httr2_1.3.0 pwalign_1.9.1
#> [13] rstatix_1.1.0 processx_3.9.0
#> [15] lattice_0.23-1 MASS_7.3-66
#> [17] backports_1.5.1 magrittr_2.0.5
#> [19] sass_0.4.10 rmarkdown_2.32
#> [21] jquerylib_0.1.4 yaml_2.3.12
#> [23] otel_0.2.0 ggtangle_0.1.3
#> [25] cowplot_1.2.0 DBI_1.3.0
#> [27] RColorBrewer_1.1-3 abind_1.4-8
#> [29] GenomicRanges_1.65.4 purrr_1.2.2
#> [31] R.utils_2.13.0 BiocGenerics_0.59.12
#> [33] msigdbr_26.1.1 RCurl_1.98-1.20
#> [35] yulab.utils_0.2.5 tweenr_2.0.3
#> [37] rappdirs_0.3.4 aisdk_1.4.12
#> [39] gdtools_0.5.1 IRanges_2.47.5
#> [41] S4Vectors_0.51.10 enrichplot_1.99.8
#> [43] ggrepel_0.9.8 tidytree_0.4.8
#> [45] DelayedArray_0.39.8 codetools_0.2-20
#> [47] DOSE_4.7.3 ggforce_0.5.0
#> [49] tidyselect_1.2.1 aplot_0.3.2
#> [51] farver_2.1.2 matrixStats_1.5.0
#> [53] stats4_4.6.1 Seqinfo_1.3.2
#> [55] GenomicAlignments_1.49.2 jsonlite_2.0.0
#> [57] Formula_1.2-6 systemfonts_1.3.2
#> [59] tools_4.6.1 ggnewscale_0.5.2
#> [61] treeio_1.37.1 PFAM.db_3.22.0
#> [63] Rcpp_1.1.2 glue_1.8.1
#> [65] SparseArray_1.13.4 gridExtra_2.3.1
#> [67] BiocBaseUtils_1.15.1 xfun_0.61
#> [69] MatrixGenerics_1.25.0 qvalue_2.45.0
#> [71] dplyr_1.2.1 withr_3.0.3
#> [73] BiocManager_1.30.27 fastmap_1.2.0
#> [75] callr_3.8.0 digest_0.6.39
#> [77] R6_2.6.1 gridGraphics_0.5-1
#> [79] GO.db_3.23.1 dichromat_2.0-1
#> [81] RSQLite_3.53.3 cigarillo_1.3.1
#> [83] R.methodsS3_1.8.2 tidyr_1.3.2
#> [85] generics_0.1.4 fontLiberation_0.1.0
#> [87] data.table_1.18.6.1 rtracklayer_1.73.0
#> [89] S4Arrays_1.13.2 httr_1.4.9
#> [91] htmlwidgets_1.6.4 scatterpie_0.2.6
#> [93] pkgconfig_2.0.3 gtable_0.3.6
#> [95] blob_1.3.0 S7_0.2.2
#> [97] XVector_0.53.0 clusterProfiler_4.21.3
#> [99] htmltools_0.5.9 fontBitstreamVera_0.1.1
#> [101] carData_3.0-6 bookdown_0.48
#> [103] scales_1.4.0 Biobase_2.73.2
#> [105] png_0.1-9 ggfun_0.2.2
#> [107] knitr_1.52 reshape2_1.4.5
#> [109] rjson_0.2.23 nlme_3.1-171
#> [111] curl_8.0.0 org.Hs.eg.db_3.23.1
#> [113] cachem_1.1.0 stringr_1.6.0
#> [115] parallel_4.6.1 AnnotationDbi_1.75.2
#> [117] restfulr_0.0.17 pillar_1.11.1
#> [119] grid_4.6.1 vctrs_0.7.3
#> [121] ggpubr_1.0.0 car_3.1-5
#> [123] tidydr_0.0.7 cluster_2.1.8.3
#> [125] evaluate_1.0.5 tinytex_0.61
#> [127] magick_2.9.1 cli_3.6.6
#> [129] compiler_4.6.1 Rsamtools_2.29.0
#> [131] rlang_1.3.0 crayon_1.5.3
#> [133] ggsignif_0.6.4 labeling_0.4.3
#> [135] ps_1.9.3 plyr_1.8.9
#> [137] fs_2.1.0 ggiraph_0.9.6
#> [139] stringi_1.8.9 assertthat_0.2.1
#> [141] Biostrings_2.81.9 lazyeval_0.2.3
#> [143] GOSemSim_2.39.3 fontquiver_0.2.1
#> [145] Matrix_1.7-6 patchwork_1.3.2
#> [147] bit64_4.8.6 ggplot2_4.0.3
#> [149] KEGGREST_1.53.6 SummarizedExperiment_1.43.0
#> [151] igraph_2.3.4 broom_1.0.13
#> [153] memoise_2.0.1 bslib_0.12.0
#> [155] ggtree_4.3.0 bit_4.6.0
#> [157] ape_5.8-1 gson_0.2.2