ORF-aware transcript-pair selection

Purpose and scope

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.

Step 1: validate the input and establish the comparison

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.

Step 2: build the reference index once

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.

Step 3: require the candidate to represent its assigned form

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.

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.

Step 4: construct candidate pairs

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.

Step 5: identify the event-specific region

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.

Step 6: compare the coding region outside the event

Exon-level cds_gen_start and cds_gen_stop define the coding bases. Each pair is compared within a coding window:

• When both transcripts are coding, the window is the span coded by either of them (context_window = "pair_cds_span").
• When a member is noncoding, every such pair of the comparison shares one window, the span coded by any candidate ("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:

• the plain skip isoform has the same exons as the inclusion isoform outside the event, but its annotated CDS stops early at the premature stop;
• another isoform also skips the exon but restores the frame with a second change, such as a 2-nt shorter acceptor downstream, so its CDS stays long.

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:

1. A sequence whose length equals the annotated CDS length is used as CDS.
2. An additional terminal TAA/TAG/TGA is removed when it accounts for exactly three extra bases. This accommodates the package’s GENCODE FASTA workflow.
3. A sequence whose length equals the total exon length is sliced to its CDS using exon coordinates and strand.
4. Other length mismatches are reported; the sequence is not silently trimmed or assumed to represent the ORF.

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.

Step 7: choose the pair transparently

Sort by the following rules, in this exact order:

1. More coding transcripts, meaning transcripts with an annotated CDS: two, then one, then zero.
2. Better worst-member TSL tier: both within 1–3, then a member at 4–5, then a member with unknown TSL.
3. Stronger evidence: cds_sequences, cds_annotation, exon_annotation, then no_outside_event_context.
4. Higher context_similarity.
5. Higher orf_similarity.
6. Fewer members that are not exact structural matches (n_inexact): exact outer boundaries win ties.
7. Better worst numeric TSL, then lower summed TSL. Unknown TSL is ranked as 6.
8. Stable case and control transcript IDs for reproducible 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.

Step 7b: fall back to an approximate pair

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.

Step 8: inspect and analyze the result

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.

Differences from the legacy protocol

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.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: Etc/UTC
#> 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] buildtools_1.0.0            RColorBrewer_1.1-3         
#>  [29] abind_1.4-8                 GenomicRanges_1.65.4       
#>  [31] purrr_1.2.2                 R.utils_2.13.0             
#>  [33] BiocGenerics_0.59.12        msigdbr_26.1.1             
#>  [35] RCurl_1.98-1.20             yulab.utils_0.2.5          
#>  [37] tweenr_2.0.3                rappdirs_0.3.4             
#>  [39] aisdk_1.4.12                gdtools_0.5.1              
#>  [41] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [43] enrichplot_1.99.7           ggrepel_0.9.8              
#>  [45] tidytree_0.4.8              maketools_1.3.2            
#>  [47] DelayedArray_0.39.8         codetools_0.2-20           
#>  [49] DOSE_4.7.3                  ggforce_0.5.0              
#>  [51] tidyselect_1.2.1            aplot_0.3.2                
#>  [53] farver_2.1.2                matrixStats_1.5.0          
#>  [55] stats4_4.6.1                Seqinfo_1.3.2              
#>  [57] GenomicAlignments_1.49.2    jsonlite_2.0.0             
#>  [59] Formula_1.2-6               systemfonts_1.3.2          
#>  [61] tools_4.6.1                 ggnewscale_0.5.2           
#>  [63] treeio_1.37.1               PFAM.db_3.22.0             
#>  [65] Rcpp_1.1.2                  glue_1.8.1                 
#>  [67] SparseArray_1.13.4          gridExtra_2.3.1            
#>  [69] BiocBaseUtils_1.15.1        xfun_0.61                  
#>  [71] MatrixGenerics_1.25.0       qvalue_2.45.0              
#>  [73] dplyr_1.2.1                 withr_3.0.3                
#>  [75] BiocManager_1.30.27         fastmap_1.2.0              
#>  [77] callr_3.8.0                 digest_0.6.39              
#>  [79] R6_2.6.1                    gridGraphics_0.5-1         
#>  [81] GO.db_3.23.1                RSQLite_3.53.3             
#>  [83] cigarillo_1.3.1             R.methodsS3_1.8.2          
#>  [85] tidyr_1.3.2                 generics_0.1.4             
#>  [87] fontLiberation_0.1.0        data.table_1.18.6.1        
#>  [89] rtracklayer_1.73.0          S4Arrays_1.13.2            
#>  [91] httr_1.4.9                  htmlwidgets_1.6.4          
#>  [93] scatterpie_0.2.6            pkgconfig_2.0.3            
#>  [95] gtable_0.3.6                blob_1.3.0                 
#>  [97] S7_0.2.2                    XVector_0.53.0             
#>  [99] sys_3.4.3                   clusterProfiler_4.21.3     
#> [101] htmltools_0.5.9             fontBitstreamVera_0.1.1    
#> [103] carData_3.0-6               scales_1.4.0               
#> [105] Biobase_2.73.2              png_0.1-9                  
#> [107] ggfun_0.2.1                 knitr_1.52                 
#> [109] reshape2_1.4.5              rjson_0.2.23               
#> [111] nlme_3.1-171                curl_8.0.0                 
#> [113] org.Hs.eg.db_3.23.1         cachem_1.1.0               
#> [115] stringr_1.6.0               parallel_4.6.1             
#> [117] AnnotationDbi_1.75.2        restfulr_0.0.17            
#> [119] pillar_1.11.1               grid_4.6.1                 
#> [121] vctrs_0.7.3                 ggpubr_1.0.0               
#> [123] car_3.1-5                   tidydr_0.0.6               
#> [125] cluster_2.1.8.3             evaluate_1.0.5             
#> [127] cli_3.6.6                   compiler_4.6.1             
#> [129] Rsamtools_2.29.0            rlang_1.3.0                
#> [131] crayon_1.5.3                ggsignif_0.6.4             
#> [133] labeling_0.4.3              ps_1.9.3                   
#> [135] plyr_1.8.9                  fs_2.1.0                   
#> [137] ggiraph_0.9.6               stringi_1.8.9              
#> [139] assertthat_0.2.1            Biostrings_2.81.9          
#> [141] lazyeval_0.2.3              GOSemSim_2.39.3            
#> [143] fontquiver_0.2.1            Matrix_1.7-6               
#> [145] patchwork_1.3.2             bit64_4.8.6                
#> [147] ggplot2_4.0.3               KEGGREST_1.53.6            
#> [149] SummarizedExperiment_1.43.0 igraph_2.3.4               
#> [151] broom_1.0.13                memoise_2.0.1              
#> [153] bslib_0.12.0                ggtree_4.3.0               
#> [155] bit_4.6.0                   ape_5.8-1                  
#> [157] gson_0.2.1