Package 'SpliceImpactR'

Description

Add one part to an existing SpliceImpactResult

Usage

add_splice_part(
  obj,
  data = NULL,
  res = NULL,
  res_di = NULL,
  matched = NULL,
  sample_frame = NULL,
  hits_final = NULL
)

Arguments

Value

Updated 'SpliceImpactResult'

Examples

obj <- as_splice_impact_result()
res <- data.table::data.table(
  event_id = c("E1", "E1"),
  form = c("inc", "exc"),
  inc = c("100-110", "120-130"),
  exc = c("120-130", "100-110"),
  chr = c("chr1", "chr1"),
  strand = c("+", "+"),
  gene_id = c("ENSG000001", "ENSG000001"),
  padj = c(0.01, 0.01),
  delta_psi = c(0.25, -0.25)
)
obj <- add_splice_part(obj, res = res)
print(as_dt_from_s4(obj, "di_events"))

Description

Convert S4 slots back to data.table

Usage

as_dt_from_s4(
  x,
  slot = c("raw_events", "di_events", "res_di", "matched", "sample_frame",
    "hits_sequences", "paired_hits"),
  keep_internal_keys = FALSE
)

Arguments

Value

'data.table'

Examples

res <- data.table::data.table(
  event_id = c("E1", "E1"),
  form = c("inc", "exc"),
  inc = c("100-110", "120-130"),
  exc = c("120-130", "100-110"),
  chr = c("chr1", "chr1"),
  strand = c("+", "+"),
  gene_id = c("ENSG000001", "ENSG000001"),
  padj = c(0.01, 0.01),
  delta_psi = c(0.20, -0.20)
)
obj <- as_splice_impact_result(res = res)
print(as_dt_from_s4(obj, "di_events"))

Description

Accepts any subset of 'data', 'res', and 'hits_final'. Missing pieces are stored as empty valid slots and can be added later with [add_splice_part()].

Usage

as_splice_impact_result(
  data = NULL,
  res = NULL,
  res_di = NULL,
  matched = NULL,
  sample_frame = NULL,
  hits_final = NULL,
  metadata = list()
)

Arguments

Value

'SpliceImpactResult'

Examples

raw <- data.table::data.table(
  event_id = c("E1", "E1"),
  form = c("inc", "exc"),
  sample = c("S1", "S1"),
  chr = c("chr1", "chr1"),
  strand = c("+", "+"),
  inc = c("100-110", "120-130"),
  exc = c("120-130", "100-110"),
  psi = c(0.70, 0.30),
  inclusion_reads = c(70, 30),
  exclusion_reads = c(30, 70)
)
res <- data.table::data.table(
  event_id = c("E1", "E1"),
  form = c("inc", "exc"),
  inc = c("100-110", "120-130"),
  exc = c("120-130", "100-110"),
  chr = c("chr1", "chr1"),
  strand = c("+", "+"),
  gene_id = c("ENSG000001", "ENSG000001"),
  padj = c(0.01, 0.01),
  delta_psi = c(0.20, -0.20)
)
hits <- data.table::data.table(
  event_id = "E1",
  event_type = "SE",
  gene_id = "ENSG000001",
  chr = "chr1",
  strand = "+",
  transcript_id_control = "TX1",
  transcript_id_case = "TX2",
  inc_control = "100-110",
  inc_case = "100-115",
  exc_control = "120-130",
  exc_case = "121-130",
  n_ppi = 1L,
  diff_n = 1L
)
obj <- as_splice_impact_result(data = raw, res = res, hits_final = hits)
obj

Description

Merges sequence data (transcript and protein sequences) onto an input table of splicing events or transcript annotations, by matching on transcript_id. When multiple sequences share the same transcript_id, the function keeps the entry with a non-missing protein_id and the longest protein_seq.

Usage

attach_sequences(x, sequences, return_class = c("auto", "data.table", "S4"))

Arguments

Details

The join is left-sided: all rows from x are preserved. Duplicate transcript_ids in sequences are resolved internally based on protein presence and sequence length.

Value

A data.table with the same rows as x and appended sequence columns (transcript_seq, protein_seq, etc.).

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annots <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annots$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annots$sequences)
print(x_seq)

Description

Computes per-event differences in the **HIT index** (Hybrid Intron Tolerance) between control and test conditions, performs t test, adjusts for multiple testing, and produces summary visualizations.

Usage

compare_hit_index(
  sample_df,
  condition_map = c(control = "control", test = "case"),
  minimum_proportion = 0.5,
  top_n_heatmap = 5000,
  sig_delta = 0.2,
  fdr_thresh = 0.05
)

Arguments

Details

For each exon-level HIT index event:

• Computes mean HIT per group ('control', 'test')

• Performs a t test if both groups have >50

• Adjusts p-values via Benjamini-Hochberg FDR

• Computes signed ('diff_HIT') and absolute ('delta_HIT') changes

The function returns both the full per-event results table and a combined 2x2 patchwork plot containing:

1. Control vs Test mean HIT scatter

2. Top |deltaHIT| heatmap (control/test columns)

3. Volcano plot of |deltaHIT| vs -log10(FDR)

4. Density of deltaHIT distribution

Value

A named list with:

'results'

Data.table of per-event statistics (means, p, FDR, deltaHIT)

'plot'

Patchwork object with 4-panel summary visualization

See Also

[.getHITindex()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_compare <- compare_hit_index(sample_frame, condition_map = c(control = "control", test = "case"))
print(hit_compare)

Description

Wrapper function that performs sequence alignment (via [compare_sequences_alignment()]) and frame-shift analysis (via [compare_frames()]) for a complete set of inclusion/exclusion transcript pairs. It then summarizes each event by a high-level classification label.

Usage

compare_sequence_frame(
  complete_hits,
  ann,
  return_class = c("auto", "data.table", "S4")
)

Arguments

Details

'compare_sequence_frame()' is a convenience function that integrates sequence and frame comparison stages in one call, producing an annotated table suitable for downstream summarization or visualization.

The summary label 'summary_classification' follows this precedence: 1. '"NMD"' - a member is annotated as nonsense-mediated decay ('transcript_type'). It is taken to make no protein; 'frame_call' still records the frameshift that usually causes the decay. 2. '"Match"' - identical protein sequences. 3. '"Rescue"' - frame restored downstream. 4. '"FrameShift"' - frame disrupted. 5. Otherwise, inherited from 'pc_class'.

Biotypes are taken from 'transcript_type_case'/'transcript_type_control' when the pairs carry them (both matchers add them), otherwise from 'ann'.

Value

A 'data.table' (or updated 'SpliceImpactResult' when 'return_class' resolves to S4) containing all columns from 'complete_hits', plus:

frame_call

Result from [compare_frames()].

rescue

Rescue classification.

summary_classification

One of '"NMD"', '"FrameShift"', '"Rescue"', '"Match"', or the original 'pc_class'.

transcript_type_case, transcript_type_control

Annotated biotypes, added from 'ann' when the input lacks them.

See Also

[compare_frames()], [compare_sequences_alignment()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annots <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annots$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annots$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annots$annotations)
print(seq_compare)

Description

Performs transcript- and protein-level global alignments between included and excluded isoform sequences to quantify sequence similarity, coding differences, and exon coverage differences.

Usage

compare_sequences_alignment(
  hits,
  annotations,
  include_sequences = FALSE,
  verbose = TRUE
)

Arguments

Details

This function wraps internal helpers for alignment (.align_dna(), .align_aa()), exon length summarization (.sum_exon_lengths()), and protein-coding classification (.pc_class()).

Value

A data.table with added columns describing sequence length, coding length, and alignment similarity metrics, including:

• pc_class: protein-coding status ("protein_coding", "onePC", "noPC")

• prot_len_case/control: protein sequence lengths

• tx_len_case/control: transcript sequence lengths

• exon_cds_len_*, exon_len_*: summed exon/CDS lengths

• dna_pid, dna_score, prot_pid, prot_score: alignment metrics

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annots <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annots$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annots$sequences)
pairs <- get_pairs(x_seq, source="multi")
aligned <- compare_sequences_alignment(pairs, annots$annotations)
print(aligned)

Description

Builds a matched-like table for pairs of transcripts by extracting all coding exon coordinates from annotations.

Usage

compare_transcript_pairs(transcript_pairs, annotations)

Arguments

Value

data.table mimicking 'matched' structure, ready for downstream comparison.

Examples

annotation_df <- load_example_data("annotation_df")$annotation_df
pairs <- data.frame(
    transcript1 = c("ENST00000337907", "ENST00000426559"),
    transcript2 = c("ENST00000400908", "ENST00000399728")
)
matched <- compare_transcript_pairs(pairs, annotation_df$annotations)
print(matched)

Description

Convenience wrapper for [enrich_domains_hypergeo()] that runs the enrichment test separately for each database (e.g. '"Pfam"', '"SMART"') and combines the results.

Usage

enrich_by_db(hits, background, dbs, ...)

Arguments

Value

Combined 'data.table' with an added 'database' column.

See Also

[enrich_by_event()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "annotated",
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)

enriched_domains <- enrich_by_db(hits_domain, bg, dbs = 'interpro')
print(enriched_domains)

Description

Convenience wrapper for [enrich_domains_hypergeo()] that runs the enrichment test separately for each event type and combines results.

Usage

enrich_by_event(hits, background, events, ...)

Arguments

Value

Combined 'data.table' with an added 'event_type' column.

See Also

[enrich_by_db()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "annotated",
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)

enriched_domains <- enrich_by_event(hits_domain, bg, events = 'AFE', db_filter = 'interpro')
print(enriched_domains)

Description

Tests whether particular protein domains are overrepresented among inclusion/exclusion transcript pairs (foreground) relative to a matched background set, using the hypergeometric test.

Usage

enrich_domains_hypergeo(
  hits,
  background,
  domain_col_fg = "either_domains_list",
  domain_col_bg = "total_sd_domains",
  event_col = "event_type",
  event_filter = NULL,
  db_filter = NULL,
  min_fg_count = 2,
  delim = "\\|"
)

Arguments

Details

Each domain identifier is counted once per transcript pair based on its presence in a list column (e.g. 'either_domains_list'). The probability of observing at least 'k' such pairs is computed under the hypergeometric distribution

$P(X \ge k), \quad X \sim \mathrm{Hypergeom}(M, B-M, K)$

where:

• K = number of foreground pairs,

• B = number of background pairs,

• M = background count of pairs containing the domain,

• k = foreground count of pairs containing the domain.

P-values are Benjamini-Hochberg adjusted ('padj'). Foreground and background pairs without a tested domain change are excluded, so both are drawn from the same conditional population. With 'db_filter', that population is pairs with a change in the selected databases. When transcript IDs are supplied in both tables, repeated unordered pairs are counted once. Foreground pairs and their changed domains must belong to the corresponding background pairs. Odds ratios compare foreground pairs with the rest of the background. Without transcript IDs, rows must already be distinct pair observations and the caller is responsible for membership; impossible counts are rejected. Counts are per transcript pair, not per event: a multi-site event (for example an AFE with several first exons) can contribute one pair per comparison, which 'n_event_comparisons' in the hits identifies.

Optionally, analyses can be restricted by event type or database prefix (e.g. '"Pfam"', '"SMART"') and domains with fewer than 'min_fg_count' foreground occurrences are skipped.

Value

A 'data.table' with one row per domain, including:

'domain_id'

Domain identifier (database prefix removed).

'db'

Database prefix (e.g. '"Pfam"').

'k', 'K'

Foreground domain count and total foreground pairs.

'M', 'B'

Background domain count and total background pairs.

'fg_prop', 'bg_prop'

Proportion of pairs with the domain.

'OR'

Odds ratio (Haldane-Anscombe corrected).

'pval', 'padj'

Raw and BH-adjusted p-values.

'events'

Event IDs contributing to the domain count.

See Also

* [enrich_by_event()] run per event type * [enrich_by_db()] run per database * [add_domain_columns()] attach domain lists to hits

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "annotated",
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)

enriched_domains <- enrich_domains_hypergeo(hits_domain, bg, db_filter = 'interpro')
print(enriched_domains)

Description

Filters 'paired_hits' using one or more logical expressions evaluated in the paired-hit table, then synchronizes all event-linked slots ('segments', 'res_di', 'di_events', 'matched', 'raw_events').

Usage

filter_spliceimpact_hits(obj, ..., keep_sample_frame = TRUE)

Arguments

Value

Filtered [SpliceImpactResult].

Examples

hits <- data.table::data.table(
  event_id = c("E1", "E2"),
  event_type = c("SE", "A3SS"),
  gene_id = c("ENSG000001", "ENSG000002"),
  chr = c("chr1", "chr2"),
  strand = c("+", "-"),
  transcript_id_control = c("TX1", "TX3"),
  transcript_id_case = c("TX2", "TX4"),
  inc_control = c("100-110", "200-210"),
  inc_case = c("100-115", "205-215"),
  exc_control = c("120-130", "220-230"),
  exc_case = c("121-130", "225-235"),
  n_ppi = c(1L, 0L),
  diff_n = c(1L, 0L),
  frame_call = c("Match", "Frameshift")
)
obj <- as_splice_impact_result(hits_final = hits)
obj_keep <- filter_spliceimpact_hits(obj, n_ppi > 0L)
print(as_dt_from_s4(obj_keep, "paired_hits"))

Description

This function loads GENCODE gene models (GTF), processes exon annotations, extracts transcript and protein sequences, identifies hybrid exons, and optionally caches the processed objects for future fast access.

Usage

get_annotation(
  load = c("link", "path", "cached", "test"),
  base_dir = NULL,
  species = c("human", "mouse"),
  release = 45,
  gtf_path = NULL,
  transcript_path = NULL,
  translation_path = NULL,
  filter_tsl = c("1", "2", "3")
)

Arguments

Details

Four loading modes are supported:

'test'

Load small internal test data shipped with the package.

'cached'

Load previously processed objects from package 'BiocFileCache'.

'path'

Read local GTF and FASTA files, process, then cache processed objects in package 'BiocFileCache'.

'link'

Download GENCODE files from URLs, process, then cache processed objects in package 'BiocFileCache'. Optional 'gtf_path', 'transcript_path', and 'translation_path' are only used as overrides when they are existing local files or valid URLs; otherwise downloaded GENCODE assets are used.

Processed objects are cached in package 'BiocFileCache' entries:

annotation/{species}/v{release}/tsl-{...}/annotations.rds
annotation/{species}/v{release}/tsl-{...}/sequences.rds
annotation/{species}/v{release}/tsl-{...}/hybrids.rds

Setting 'filter_tsl = NULL' disables transcript-support filtering and uses separate 'tsl-none' processed cache entries. Use the same 'filter_tsl' when creating and loading the cache. Gene-biotype and incomplete-CDS tag filters still apply. The 'test' mode loads its bundled subset as stored.

Value

A list with:

'annotations'

Processed long-format GTF as 'data.table'

'sequences'

List with elements 'transcripts' and 'proteins' (or 'NULL' if not loaded)

'hybrids'

Hybrid exon annotation list

Examples

# Load bundled test data
ann <- load_example_data("annotation_df")$annotation_df
print(ann)

# Load from local files and cache processed objects
# ann <- get_annotation(
#   load = "path",
#   gtf_path = "/downloaded_gtf_directory/gencode.v45.annotation.gtf.gz",
#   transcript_path = "/downloaded_gtf_directory/gencode.v45.pc_transcripts.fa.gz",
#   translation_path = "/downloaded_gtf_directory/gencode.v45.pc_translations.fa.gz"
# )

# Download files, process, and cache to a custom cache root
# ann <- get_annotation(
#   load = "link",
#   base_dir = "/project/annotation_cache/"
# )

# Load from cached RDS (fast)
# ann <- get_annotation(
#   load = "cached",
#   base_dir = "/project/annotation_cache/"
# )

# Disable TSL filtering, including for transcripts without a TSL value
# ann <- get_annotation(load = "link", filter_tsl = NULL)
# ann <- get_annotation(load = "cached", filter_tsl = NULL)

Description

Constructs a background dataset of transcript pairs suitable for domain-level or exon-level enrichment analyses. Depending on the source parameter, the function can derive the background from HIT index results, annotated transcripts, or a user-provided list of transcript IDs. The returned table has a 'gene_universe' attribute containing eligible genes before filtering for domain-changing pairs. Gene-level enrichment should intersect this universe with the genes actually evaluated for its foreground. 'feature_gene_universe' contains genes with at least one eligible pair whose features differ, i.e. genes that could appear in a domain or PPI foreground; use it for domain and PPI gene enrichment.

Usage

get_background(
  source = c("annotated", "hit_index", "user-given"),
  input = NULL,
  annotations,
  protein_features,
  keep_annotated_first_last = TRUE,
  minOverlap = 0.8,
  BPPARAM = BiocParallel::bpparam()
)

Arguments

Details

Use source = "annotated" (the default) for enrich_domains_hypergeo(). A "hit_index" background keeps one annotated transcript per exon observed in the samples, so it usually lacks the transcript pairs the matcher chose and domain enrichment rejects it; its gene_universe and feature_gene_universe still serve gene-level enrichment.

Value

A data.table in which each row represents a transcript pair annotated with gene ID, CDS/exon length differences, and protein domain differences. Columns include: gene_id, transcript_id_1, transcript_id_2, prot_aa_1, prot_aa_2, domains_1, domains_2, and related summary metrics.

Examples

annots <- load_example_data("annotation_df")$annotation_df
interpro_features <- get_protein_features(c("interpro"), annots$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

# The default, source = "annotated", pairs every annotated transcript (see
# enrich_domains_hypergeo()). A HIT-index background serves gene-level
# universes:
ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annots$annotations,
                     protein_features = protein_feature_total)

# Or from user-supplied transcript IDs
tx_ids <- c("ENST00000466994","ENST00000484435")
bg_user <- get_background(source = "user-given",
                          input = tx_ids,
                          annotations = annots$annotations,
                          protein_features = protein_feature_total)

Description

Aggregates all available feature tables (e.g., biomaRt + manual) into a single unified long-format annotation table.

Usage

get_comprehensive_annotations(
  protein_feature_list,
  load_path_list = NULL,
  save_path = NULL
)

Arguments

Value

A combined data.table containing all unique feature rows.

Examples

annotation_df <- load_example_data("annotation_df")$annotation_df
user_df <- data.frame(
  ensembl_transcript_id = c(
    "ENST00000511072","ENST00000374900","ENST00000373020","ENST00000456328",
    "ENST00000367770","ENST00000331789","ENST00000335137","ENST00000361567",
    NA,                    "ENST00000380152"
  ),
  ensembl_peptide_id = c(
    "ENSP00000426975", NA,                   "ENSP00000362048","ENSP00000407743",
    "ENSP00000356802","ENSP00000326734", NA,                  "ENSP00000354587",
    "ENSP00000364035", NA
  ),
  name = c(
    "Low complexity","Transmembrane helix","Coiled-coil","Signal peptide",
    "Transmembrane helix","Low complexity","Coiled-coil","Transmembrane helix",
    "Signal peptide","Low complexity"
  ),
  start = c(80L, 201L, 35L, 1L, 410L, 150L, 220L, 30L, 1L, 300L),
  stop  = c(120L,223L, 80L, 20L, 430L, 190L, 260L, 55L, 24L, 360L),
  database   = c("seg","tmhmm","ncoils","signalp","tmhmm","seg","ncoils","tmhmm","signalp", NA),
  alt_name   = c(NA,"TMhelix",NA,"SignalP-noTM", "TMhelix", NA, NA, "TMhelix", "SignalP-TAT", NA),
  feature_id = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)
)
user_features <- get_manual_features(user_df, annotation_df$annotations)
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(user_features, interpro_features))
print(protein_feature_total)

Description

Select genes whose splicing passes significance and effect-size thresholds.

Usage

get_di_gene_enrichment(hits, padj_threshold, delta_psi_threshold)

Arguments

Value

Character vector of gene IDs.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)

annotation_df <- load_example_data("annotation_df")$annotation_df
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
enrichment <- get_enrichment(get_di_gene_enrichment(res, .05, .1), intersect(attr(bg, "gene_universe"), res$gene_id[is.finite(res$padj)]), species = 'human', 'ensembl', 'MSigDB:H')
print(enrichment)

Description

Performs per-site differential inclusion testing from a hit-index or junction-form table. Each site is modeled with a quasi-binomial GLM ('psi_adj ~ condition') to estimate deltaPSI and significance, optionally using parallel processing.

Usage

get_differential_inclusion(
  DT,
  min_total_reads = 10L,
  minimum_proportion_containing_event = 0.5,
  terminal_fill = "event_max",
  cooks_cutoff = "Inf",
  adjust_method = "fdr",
  verbose = TRUE,
  parallel_glm = TRUE,
  chunk_size_glm = 1000,
  BPPARAM = BiocParallel::SerialParam(),
  return_class = c("auto", "data.table", "S4")
)

Arguments

Details

The function filters out low-coverage and low-presence events, optionally fills AFE/ALE sites with zero counts where necessary, and then applies site-level GLMs. Parallelization uses [BiocParallel] back-ends for reproducibility across platforms. To run in parallel, supply 'BPPARAM' (for example [BiocParallel::MulticoreParam()] on Linux/macOS or [BiocParallel::SnowParam()] on Windows).

Value

If 'return_class' resolves to '"data.table"', a 'data.table' with one row per site containing:

• Metadata columns ('site_id', 'event_type', 'event_id', 'gene_id', ...)

• Sample counts ('n_samples', 'n_control', 'n_case')

• Mean PSI per group ('mean_psi_ctrl', 'mean_psi_case')

• deltaPSI ('delta_psi')

• Raw and adjusted p-values ('p.value', 'padj')

• Maximum Cook's distance ('cooks_max')

See Also

[fit_sites_parallel()], [cutoff_num()], [stats::glm()], [BiocParallel::bplapply()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
head(res)

Description

Return genes whose inclusion/exclusion isoforms show unique protein domain gain or loss.

Usage

get_domain_gene_for_enrichment(hits)

Arguments

Value

Character vector of gene IDs.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
enrichment <- get_enrichment(get_domain_gene_for_enrichment(hits_domain), intersect(attr(bg, "feature_gene_universe"), res$gene_id[is.finite(res$padj)]), species = 'human', 'ensembl', 'MSigDB:H')
print(enrichment)

Description

Annotates each splicing event with protein domains that are gained, lost, or uniquely present in inclusion or exclusion isoforms.

Usage

get_domains(
  hits,
  exon_features,
  show_protein_domains = FALSE,
  return_class = c("auto", "data.table", "S4")
)

Arguments

Details

Internally, this function builds a domain lookup table from an exon feature annotation (e.g. InterPro, Pfam) and extracts per-exon and per-transcript domain lists for each isoform in 'hits'. Differences between the inclusion ('*_case') and exclusion ('*_control') isoforms are then summarized as:

* 'case_only_domains': domains unique to the inclusion isoform * 'control_only_domains': domains unique to the exclusion isoform * 'diff_n': total number of non-shared domains

If 'show_protein_domains = TRUE', additional columns report full domain sets across the entire inclusion/exclusion proteins.

Value

The input 'hits' table with added columns (or updated 'SpliceImpactResult' when 'return_class' resolves to S4):

• 'domains_exons_case', 'domains_exons_control' domains found on event exons

• 'case_only_domains', 'control_only_domains' domains unique to each isoform

• 'case_only_domains_list', 'control_only_domains_list', 'either_domains_list' list-columns

• 'case_only_n', 'control_only_n', 'diff_n' domain counts

• optionally, 'domains_protein_case' / 'domains_protein_control'

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)
print(hits_domain)

Description

Perform over-representation analysis for a foreground gene set, optionally against a background universe, using GO, MSigDB, and Reactome categories.

Usage

get_enrichment(
  foreground,
  background = NULL,
  species = c("human", "mouse"),
  gene_id_type = c("symbol", "ensembl"),
  sources = c("GO:BP", "GO:MF", "GO:CC", "MSigDB:H", "MSigDB:C2:CP:REACTOME"),
  min_size = 10,
  max_size = 2000,
  p_adjust_cutoff = 0.05,
  simplify_go = TRUE,
  top_n_plot = 20,
  plot_type = c("dot", "bar")
)

Arguments

Details

This is a convenience wrapper around 'clusterProfiler', 'msigdbr', and optionally 'ReactomePA', producing a combined enrichment table and a quick visualization of top significant terms.

Input genes are internally mapped to Entrez IDs. Enrichment tests are performed using:

* 'clusterProfiler::enrichGO' * 'clusterProfiler::enricher' (MSigDB) * 'ReactomePA::enrichPathway' (optional)

Value

A list with:

results_per_source

List of enrichment result tables per source

results_combined

Combined enrichment table

results_signif

Filtered table by FDR cutoff

plot

A 'ggplot2' object visualizing top terms

Note

Requires these packages installed: 'clusterProfiler', 'msigdbr', 'data.table', 'AnnotationDbi', 'ggplot2', and 'org.Hs.eg.db' or 'org.Mm.eg.db'. For Reactome analysis you must also install 'ReactomePA'.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)

annotation_df <- load_example_data("annotation_df")$annotation_df
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
enrichment <- get_enrichment(res$gene_id, attr(bg, "gene_universe"), species = 'human', 'ensembl', 'MSigDB:H')
print(enrichment)

Description

Overlaps amino acid-based protein features (e.g., InterPro, TMHMM) with exon coding spans defined by transcript-relative CDS coordinates.

Usage

get_exon_features(gtf_dt, feat, inclusive = TRUE)

Arguments

Value

A data.table of overlapping feature-exon pairs with amino acid coordinates (overlap_aa_start, overlap_aa_end, overlap_aa_len) and associated exon and feature metadata.

Examples

annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)
print(exon_features)

Description

Selects a foreground gene vector for enrichment from DI results ('res') or hits-final-like tables ('hits') using a single wrapper.

Usage

get_gene_enrichment(
  mode = c("di", "ppi", "domain"),
  x = NULL,
  res = NULL,
  hits = NULL,
  padj_threshold = 0.05,
  delta_psi_threshold = 0.1
)

Arguments

Details

- 'mode = "di"' uses differential inclusion columns ('gene_id', 'padj', 'delta_psi') - 'mode = "ppi"' uses hits-final-like columns ('gene_id', 'n_ppi') - 'mode = "domain"' uses hits-final-like columns ('gene_id', 'diff_n')

Inputs can be a 'data.table'/'data.frame' or a 'SpliceImpactResult' S4 object.

Value

Character vector of unique gene IDs.

Examples

res <- data.table::data.table(
  gene_id = c("ENSG000001", "ENSG000001", "ENSG000002"),
  padj = c(0.01, 0.20, 0.03),
  delta_psi = c(0.25, 0.05, -0.30)
)
hits <- data.table::data.table(
  gene_id = c("ENSG000001", "ENSG000002", "ENSG000003"),
  n_ppi = c(1L, 0L, 2L),
  diff_n = c(0L, 1L, 2L)
)
get_gene_enrichment(mode = "di", res = res)
get_gene_enrichment(mode = "ppi", hits = hits)
get_gene_enrichment(mode = "domain", hits = hits)

hits_s4 <- data.table::data.table(
  event_id = c("E1", "E2", "E3"),
  event_type = c("SE", "A3SS", "MXE"),
  gene_id = c("ENSG000001", "ENSG000002", "ENSG000003"),
  chr = c("chr1", "chr1", "chr2"),
  strand = c("+", "+", "-"),
  transcript_id_control = c("TX1", "TX3", "TX5"),
  transcript_id_case = c("TX2", "TX4", "TX6"),
  inc_control = c("100-110", "200-210", "300-310"),
  inc_case = c("100-115", "205-215", "305-315"),
  exc_control = c("120-130", "220-230", "320-330"),
  exc_case = c("121-130", "225-235", "325-335"),
  n_ppi = c(1L, 0L, 2L),
  diff_n = c(0L, 1L, 2L)
)
obj <- as_splice_impact_result(hits_final = hits_s4)
get_gene_enrichment(mode = "ppi", x = obj)
get_gene_enrichment(mode = "domain", x = obj)

Description

Load HIT index PSI files for one or more samples/conditions.

Usage

get_hitindex(paths_df, keep_annotated_first_last = FALSE)

Arguments

Value

A standardized 'data.table' of HIT index PSI values with inclusion/exclusion and metadata.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame, keep_annotated_first_last = TRUE)
print(hit_index)

Description

Convenience accessor for core paired-hit columns

Usage

get_hits_core(x, drop_missing = TRUE, keep_internal_keys = FALSE)

Arguments

Value

'data.table' with the 'core' subset. Works for both S4 and paired-hits 'data.table' input.

Examples

hits <- data.table::data.table(
  event_id = c("E1", "E2"),
  event_type = c("SE", "A3SS"),
  gene_id = c("ENSG000001", "ENSG000002"),
  chr = c("chr1", "chr2"),
  strand = c("+", "-"),
  transcript_id_control = c("TX1", "TX3"),
  transcript_id_case = c("TX2", "TX4"),
  n_ppi = c(1L, 0L),
  diff_n = c(1L, 0L)
)
print(get_hits_core(hits))

Description

Convenience accessor for domain-focused paired-hit columns

Usage

get_hits_domain(x, drop_missing = TRUE, keep_internal_keys = FALSE)

Arguments

Value

'data.table' with the 'domain' subset. Works for both S4 and paired-hits 'data.table' input.

Examples

hits <- data.table::data.table(
  event_id = c("E1", "E2"),
  event_type = c("SE", "A3SS"),
  gene_id = c("ENSG000001", "ENSG000002"),
  chr = c("chr1", "chr2"),
  strand = c("+", "-"),
  transcript_id_control = c("TX1", "TX3"),
  transcript_id_case = c("TX2", "TX4"),
  case_only_domains = c("IPR0001", ""),
  control_only_domains = c("", "IPR0002"),
  case_only_n = c(1L, 0L),
  control_only_n = c(0L, 1L),
  diff_n = c(1L, 1L)
)
print(get_hits_domain(hits))

Description

Extracts 'paired_hits' columns from a [SpliceImpactResult] (or a paired-hits 'data.table') using predefined subset groups such as 'core', 'domain', 'ppi', and 'sequence'.

Usage

get_hits_final_view(
  x,
  col_subset = c("core"),
  cols = NULL,
  drop_missing = TRUE,
  keep_internal_keys = FALSE
)

Arguments

Value

A 'data.table' containing the selected columns. Row count and row order are preserved from the input ('SpliceImpactResult@paired_hits' or provided 'data.table').

Examples

hits <- data.table::data.table(
  event_id = c("E1", "E2"),
  event_type = c("SE", "A3SS"),
  gene_id = c("ENSG000001", "ENSG000002"),
  chr = c("chr1", "chr2"),
  strand = c("+", "-"),
  transcript_id_control = c("TX1", "TX3"),
  transcript_id_case = c("TX2", "TX4"),
  protein_id_control = c("P1", "P3"),
  protein_id_case = c("P2", "P4"),
  inc_control = c("100-110", "200-210"),
  inc_case = c("100-115", "205-215"),
  exc_control = c("120-130", "220-230"),
  exc_case = c("121-130", "225-235"),
  case_only_domains = c("IPR0001", ""),
  control_only_domains = c("", "IPR0002"),
  case_only_n = c(1L, 0L),
  control_only_n = c(0L, 1L),
  diff_n = c(1L, 1L),
  case_ppi = c("A;B", "C"),
  control_ppi = c("A", "C;D"),
  n_case_ppi = c(2L, 1L),
  n_control_ppi = c(1L, 2L),
  n_ppi = c(1L, 1L),
  dna_pid = c(0.95, 0.90),
  prot_pid = c(0.90, 0.85),
  frame_call = c("Match", "Frameshift")
)
print(get_hits_final_view(hits, col_subset = c("core", "ppi")))

Description

Convenience accessor for PPI-focused paired-hit columns

Usage

get_hits_ppi(x, drop_missing = TRUE, keep_internal_keys = FALSE)

Arguments

Value

'data.table' with the 'ppi' subset. Works for both S4 and paired-hits 'data.table' input.

Examples

hits <- data.table::data.table(
  event_id = c("E1", "E2"),
  event_type = c("SE", "A3SS"),
  gene_id = c("ENSG000001", "ENSG000002"),
  chr = c("chr1", "chr2"),
  strand = c("+", "-"),
  transcript_id_control = c("TX1", "TX3"),
  transcript_id_case = c("TX2", "TX4"),
  case_ppi = c("A;B", "C"),
  control_ppi = c("A", "C;D"),
  n_case_ppi = c(2L, 1L),
  n_control_ppi = c(1L, 2L),
  n_ppi = c(1L, 1L)
)
print(get_hits_ppi(hits))

Description

Convenience accessor for sequence/frame-focused paired-hit columns

Usage

get_hits_sequence(x, drop_missing = TRUE, keep_internal_keys = FALSE)

Arguments

Value

'data.table' with the 'sequence' subset. Works for both S4 and paired-hits 'data.table' input.

Examples

hits <- data.table::data.table(
  event_id = c("E1", "E2"),
  event_type = c("SE", "A3SS"),
  gene_id = c("ENSG000001", "ENSG000002"),
  chr = c("chr1", "chr2"),
  strand = c("+", "-"),
  transcript_id_control = c("TX1", "TX3"),
  transcript_id_case = c("TX2", "TX4"),
  protein_id_control = c("P1", "P3"),
  protein_id_case = c("P2", "P4"),
  dna_pid = c(0.95, 0.90),
  prot_pid = c(0.90, 0.85),
  frame_call = c("Match", "Frameshift")
)
print(get_hits_sequence(hits))

Description

Converts a manual feature table into the standardized long format and optionally merges it with biomaRt-derived features.

Usage

get_manual_features(
  manual_features,
  gtf_df,
  biomaRt_features = NULL,
  load_path = NULL,
  save_path = NULL
)

Arguments

Details

Rows that cannot be placed are dropped and reported in a message that names up to five of them: rows whose transcript or peptide ID is not in gtf_df, and rows on transcripts without an annotated CDS. Features that extend past the end of the CDS are truncated to it.

Value

A data.table of manual (and optionally combined) protein features.

See Also

[add_user_features()]

Examples

annotation_df <- load_example_data("annotation_df")$annotation_df
# Illustrative features on coding transcripts in the bundled annotation
# (PPFIA1, COL1A2, EPS15, NUP62); start and stop are amino-acid positions.
# The NUP62 feature is given by its peptide ID only.
user_df <- data.frame(
  ensembl_transcript_id = c("ENST00000253925", "ENST00000297268",
                            "ENST00000371730", NA),
  ensembl_peptide_id = c("ENSP00000253925", "ENSP00000297268",
                         "ENSP00000360795", "ENSP00000305503"),
  name = c("Coiled-coil", "Signal peptide", "Low complexity",
           "Low complexity"),
  start = c(100L, 1L, 600L, 1L),
  stop  = c(400L, 22L, 650L, 300L),
  database = c("ncoils", "signalp", "seg", "seg")
)
user_features <- get_manual_features(user_df, annotation_df$annotations)
print(user_features)

Description

This is a wrapper around match_events_to_annotations_vec that processes large event tables in manageable chunks to reduce memory usage. It sequentially runs matching per chunk and concatenates the results.

Usage

get_matched_events_chunked(
  events,
  annotations,
  chunk_size = 50000,
  minOverlap = 0.05,
  return_class = c("auto", "data.table", "S4"),
  verbose = TRUE
)

Arguments

Details

This function is intended for large-scale event matching across many splicing events, where running the full table at once may exceed memory limits. It can later be parallelized using future.apply or similar.

Value

A data.table with matched transcripts and exons for all events, including each matched transcript's annotated 'transcript_type' ('NA' when unmatched). The output order matches the original event order.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
print(sample_frame)

hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annots <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annots$annotations, chunk_size = 2000)
print(matched)

Description

Builds paired tables of inclusion/exclusion forms for splicing events from rMATS-like or HITindex-like inputs. In rMATS mode, events are paired when both INC and EXC forms exist for a given event ID. In HITindex mode, all positive and negative deltaPSI rows within each event are cross-joined.

Usage

get_pairs(
  x,
  source = c("paired", "multi"),
  return_class = c("auto", "data.table", "S4")
)

Arguments

Details

In source="paired" mode, only events with exactly one INC and one EXC row are retained. When delta PSI is supplied, the positive form is the case and the negative form is the control; without it, INC is the case. In source="multi" mode, all positive deltaPSI rows are joined with all negative deltaPSI rows (cartesian join) within each event. Missing sample-count diagnostics are retained as 'NA'. With matched transcript input, both modes return the same columns, including when no pairs remain. A 'transcript_type' column on the input (added by both matchers) becomes 'transcript_type_case' and 'transcript_type_control', for example to show which selected transcripts are annotated as NMD.

An event with several rising or falling sites, such as an AFE with three first exons, gives one comparison per rising x falling pair of sites, so a site can appear in several rows. 'n_event_comparisons' gives the number of comparisons its event defines (1 for a two-form event). After [keep_sig_pairs()], 'site_significant_case' and 'site_significant_control' show whether each form or site met the thresholds itself; the other sites of a passing event are kept as comparators. Counts of pairs, including the domain-enrichment foreground, are therefore counts of comparisons. Forms without a matched transcript ID cannot produce a pair. An unmatched form is not interpreted as an absent or noncoding isoform.

Value

A data.table (or updated 'SpliceImpactResult' when 'return_class' resolves to S4) where each row represents an inclusion-exclusion pair of the same event. S4 output sets 'metadata$matching' to '"legacy"' and removes 'matching_diagnostics' left by [get_ranked_pairs()].

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annots <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annots$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annots$sequences)
pairs <- get_pairs(x_seq, source="multi")
print(pairs)

Description

Return genes whose splicing affects known protein-protein interactions.

Usage

get_ppi_gene_enrichment(hits)

Arguments

Value

Character vector of gene IDs.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)
ppi <- get_ppi_interactions()             
hits_final <- get_ppi_switches(hits_domain, ppi, protein_feature_total)
bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
enrichment <- get_enrichment(get_ppi_gene_enrichment(hits_final), intersect(attr(bg, "feature_gene_universe"), res$gene_id[is.finite(res$padj)]), species = 'human', 'ensembl', 'MSigDB:H')
print(enrichment)

Description

Generation details in inst/scripts

Usage

get_ppi_interactions()

Value

A 'data.table' of interaction edges used by PPI switching utilities.

Examples

ppi_int <- get_ppi_interactions()
print(ppi_int)

Description

Adds list-cols case_ppi/control_ppi (partner genes) plus counts. Also returns (optionally useful) per-event token sets in PFAM + ELM forms.

Usage

get_ppi_switches(
  hits_domain,
  ppi,
  protein_feature_total,
  return_class = c("auto", "data.table", "S4")
)

Arguments

Value

A 'data.table' identical to 'hits_domain' with added columns (or updated 'SpliceImpactResult' when 'return_class' resolves to S4):

'case_ppi', 'control_ppi'

Lists of partner transcripts unique to inclusion or exclusion isoforms.

'n_control_ppi', 'n_control_ppi'

Counts of gained/lost interactions.

'n_ppi'

Total PPI changes (sum of both directions).

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)

annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
ppi <- get_ppi_interactions()             
hits_ppi <- get_ppi_switches(hits_domain, ppi, protein_feature_total)
print(hits_ppi)
hits_ppi[n_ppi > 0, .(event_id, gene_id, n_case_ppi, n_control_ppi, n_ppi, case_ppi, control_ppi)]

Description

Here we also remove any duplicate and overlapping domains We also add the genomic location to the name of the protein feature for downstream safeguarding and precision. This is to prevent different occurences of the same domain being called as the same in domain identification and enrichment.

Usage

get_protein_features(
  biomaRt_databases = c("interpro", "mobidblite", "seg", "ncoils", "tmhmm", "signalp",
    "elm", "gene3d", "pfam"),
  gtf_df,
  sequences = NULL,
  load_path = NULL,
  save_path = NULL,
  base_dir = NULL,
  use_cache = TRUE,
  force_refresh = FALSE,
  timeout = 600,
  ensembl_mirror = NULL,
  species = c("human", "mouse"),
  release = 111,
  test = FALSE,
  combine_overlaps = FALSE,
  ensembl_host = NULL
)

Arguments

Details

Ensembl 116 (June 2026) is the last release with BioMart, and www.ensembl.org no longer serves it, so features are fetched from 'https://e<release>.ensembl.org', which Ensembl forwards to the archive of 'release' (for example 'jan2024.archive.ensembl.org' for 111). This works with every biomaRt version; biomaRt before 2.70 cannot use the archive hosts directly. Ensembl BioMart is sometimes temporarily unavailable (HTTP 500/503 or time-outs); retry later or choose another release. 'elm' motifs come from the ELM database and are mapped to transcripts through the same BioMart host.

Value

A 'data.table' with one row per protein feature and transcript coupling

Examples

annotation_df <- load_example_data("annotation_df")$annotation_df
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, annotation_df$sequences, timeout = 600, test = TRUE)
print(interpro_features)

Description

Determines whether inclusion/exclusion events correspond to proximal or distal terminal exon usage for **AFE** (Alternative First Exon) and **ALE** (Alternative Last Exon) events, based on genomic coordinates and strand orientation.

Usage

get_proximal_shift_from_hits(hits)

Arguments

Details

Only AFE and ALE pairs are used; other event types are dropped. Each pair's case exon ('inc_case') is compared with its control exon ('inc_control') at the end facing the gene body:

• For **AFE** pairs, '"proximal"' means the case first exon lies further downstream (larger end on the '+' strand, smaller start on the '-' strand).

• For **ALE** pairs, '"proximal"' means the case last exon lies further upstream (smaller start on the '+' strand, larger end on the '-' strand).

The reverse order is '"distal"', and equal ends are '"overlap"'. Each 'inc_case' and 'inc_control' must be a single 'start-end' interval.

Value

A list with 'data', a 'data.table' of the AFE/ALE pairs ('event_id', 'event_type', 'strand', 'pos', 'neg', the two 'delta_psi' columns, a row index 'I', and 'V1' specifying '"proximal"', '"distal"' or '"overlap"'), and 'plot', the summary from [plot_prox_dist()]. Without AFE/ALE pairs, 'data' has zero rows and 'plot' is a placeholder.

See Also

[plot_prox_dist()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
proximal_output <- get_proximal_shift_from_hits(pairs)
print(proximal_output)

Description

The default transcript-pair selection in [get_splicing_impact()]; the legacy matcher, which selects each form's transcript independently, remains available with 'matching = "legacy"'. Candidates must represent the supplied AS form before coding status, transcript support and coding-region similarity are considered.

Usage

get_ranked_pairs(
  events,
  annotations,
  sequences,
  source = c("multi", "paired"),
  max_candidates = 100L,
  fallback = TRUE,
  verbose = TRUE,
  return_class = c("list", "S4")
)

Arguments

Details

Ranking is lexicographic: number of coding transcripts (descending), worst TSL tier (1–3, then 4–5, then unknown), evidence (both CDS sequences usable, CDS annotated without a usable sequence, then a member without CDS), context similarity, whole-window similarity ('orf_similarity'), number of members that are not exact structural matches, worst numeric TSL, summed TSL, and transcript IDs. TSL 1–3 is a package policy tier; within it similarity precedes exact TSL. This is not a transcript-expression likelihood model.

A transcript is coding when it has an annotated CDS; its biotype is not used. Transcripts annotated as nonsense-mediated decay therefore count as coding: when an event shifts the frame, the transcript that carries it faithfully is often annotated that way, and requiring the 'protein_coding' biotype would favour a frame-restoring isoform.

The event mask is the symmetric difference between required intervals of the two forms, plus their explicit forbidden intervals. Context similarity is the base-pair Jaccard overlap of the two transcripts' exonic positions in a coding window, outside the mask: the span coded by either transcript when both are coding, otherwise one window shared by every pair of the comparison (the span coded by any candidate, or the candidates' full span without CDS). Positions in one genome determine sequence, so this compares the coding region without depending on where each annotated reading frame ends: a frameshifted isoform keeps credit for its downstream exons and is not outranked by a frame-restoring isoform. 'orf_similarity' is the same score including the event. Protein identity is deliberately not a criterion.

Forbidden segments cannot overlap any exon, and event-defining boundaries must match exactly: SE/MXE and spliced RI check adjacent exons and junction boundaries, retained RI pieces must lie within one exon, A3SS/A5SS check the variable boundary and adjacent exon existence (an omitted partner interval cannot be checked), and terminal events check the first/last exon. A transcript that also covers every required segment is an 'exact' match. The outer ends of flanking exons, the shared A3SS/A5SS boundary and terminal TSS/PAS ends are not part of an event, so a transcript that differs only there is a 'relaxed' match ('structural_match'). Relaxed matches must still contain their form's alternative region and lose ties to exact matches. For terminal events they are used only when no transcript matches exactly, so sites sharing a splice site stay distinguishable. Partial and generic interval-only definitions are labelled in diagnostics.

With 'fallback = TRUE', a comparison without a structural pair is retried with the legacy matcher's gates: each required interval overlaps an exon by at least 5 two transcripts and 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 (0 to 2), and must be positive. Approximate pairs are ranked by 'event_agreement' first, then by the criteria above, and are labelled 'matching_tier = "approximate"' with 'approximate' structural matches; their event boundaries are unverified.

Input errors (missing metadata, strands other than '+'/'-', unknown forms, malformed intervals, duplicate definitions) stop the run before any matching and name the affected events.

Value

A list with 'pairs' (one selected pair per positive/negative form combination, compatible with [compare_sequence_frame()]; members' biotypes are in 'transcript_type_case' and 'transcript_type_control', 'NA' without a 'transcript_type' column, and 'n_event_comparisons' counts the comparisons its event defines, resolved or not), 'matched' (selected form rows, qualified by 'matching_group'), 'candidates' (eligibility, structural match, approximate eligibility and rejection reasons), 'rankings' (all scored combinations, by tier), 'events' (all input rows and candidate counts), 'unmatched' (form pairs unresolved by every tier, with the structural 'reason', the 'fallback' outcome and any candidates dropped by 'max_candidates'), 'pair_rejections' (self-pairs, combinations unable to distinguish the forms, and approximate combinations without an event difference), and 'settings'. Ties and context-score margins are descriptive diagnostics, not statistical confidence or probabilities. With 'return_class = "S4"', a [SpliceImpactResult] holding the same tables (see 'return_class').

See Also

[get_matched_events_chunked()], [get_pairs()], [get_splicing_impact()]

Examples

ex <- load_example_data(c("sample_frame", "annotation_df"))
hit_index <- get_hitindex(ex$sample_frame)
res <- get_differential_inclusion(hit_index)
res_di <- keep_sig_pairs(res)
ranked <- get_ranked_pairs(res_di, ex$annotation_df$annotations,
                           ex$annotation_df$sequences)
# One case/control transcript pair per resolved comparison
head(ranked$pairs[, c("event_id", "transcript_id_case",
                      "transcript_id_control", "matching_tier")])
table(ranked$pairs$matching_tier)
# Comparisons no tier resolved, and why
table(ranked$unmatched$reason)

Description

Converts rMATS "event" tables (SE, MXE, A3SS, A5SS, RI) into standardized scalar representations with explicit inclusion/exclusion segments for downstream genomic mapping.

Usage

get_rmats(DT)

Arguments

Details

Handles all five canonical rMATS event types (SE, MXE, A3SS, A5SS, RI), applying strand-aware logic for MXE and coordinate adjustments for A3/A5. A3SS/A5SS forms contain the long or short exon and its partner (rMATS flanking) exon, in genomic order; the excluded segment is the extension. Raw rMATS starts are converted from zero-based to one-based coordinates; ends are unchanged. Output intervals are closed, including one-base segments. Non-standard columns (e.g. IJC_SAMPLE_1) are checked for presence.

Value

A standardized 'data.table' containing:

• event_id unique event identifier.

• event_type rMATS event type (SE, MXE, etc.).

• form inclusion/exclusion form.

• gene_id, chr, strand.

• inc, exc scalar genomic segments (string: e.g. '"100-200;300-400"').

• inclusion_reads, exclusion_reads, psi - numeric metrics.

• condition, sample, source_file - carried forward if present.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
rmats <- get_rmats(load_rmats(sample_frame, use = "JCEC", event_types = c("MXE", "SE", "A3SS", "A5SS", "RI")))
print(rmats)

Description

Wrapper function to get both rmats and hit index cleanly

Usage

get_rmats_hit(
  sample_frame,
  event_types = c("ALE", "AFE", "MXE", "SE", "A3SS", "A5SS", "RI"),
  use = "JCEC",
  keep_annotated_first_last = TRUE
)

Arguments

Value

a 'data.table' for all the event types desired from the paths supplied contains: event_id (unique id for event), event_type (AS event type), form (INC/EXC), gene_id (ensembl id), strand, inc, exc (inclusion and exclusion coords) inclusion reads, exclusion reads, psi, sample, condition, source file

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
data <- get_rmats_hit(sample_frame, event_types = c("ALE", "AFE", "MXE", "SE", "A3SS", "A5SS", "RI"))
print(data)

Description

This function reads post-DI rMATS results and converts them into the standardized SpliceImpactR long format with one row per event x (INC/EXC) form.

Input can be: * a data.frame with columns 'path', 'grp1', 'grp2', 'event_type', in which case each file is read and processed; or * a single rMATS results data.frame, in which case 'event_type' must be supplied.

For each event, the function constructs paired INC and EXC entries: * 'inc' contains genomic segments included in the form; for A3SS/A5SS, the long or short exon and its partner (rMATS flanking) exon * 'exc' contains the excluded segment(s) * 'delta_psi', 'p.value', and 'padj' are assigned using the rMATS-reported values, with 'delta_psi' oriented by 'case_group' Raw rMATS coordinates are converted to one-based, closed intervals. An already standardized SpliceImpactR DI table is returned without conversion.

rMATS reports 'IncLevelDifference' as mean('IncLevel1') - mean('IncLevel2'), where group 1 holds the '–b1' samples and group 2 the '–b2' samples. SpliceImpactR uses 'delta_psi' = case - control, so 'case_group' states which rMATS group holds the case samples. With the default, group 1, 'IncLevelDifference' is used as reported; with group 2 its sign is reversed. Each import reports the group it treated as the case.

Event IDs are automatically generated (event_type:N) if not supplied. With several files, events are numbered across all files, so each ID names one event; identical event definitions are kept once. All files belong to one comparison: an event type listed under different 'grp1'/'grp2' labels is an error.

Usage

get_rmats_post_di(input, event_type = NULL, case_group = 1)

Arguments

Value

A 'data.table' with columns:

event_id

unique event identifier

event_type

splicing event type

form

"INC" or "EXC"

gene_id

gene ID

chr

chromosome

strand

strand

inc

genomic coordinates of included segment(s)

exc

genomic coordinates of excluded segment(s)

p.value

rMATS p-value

padj

FDR

delta_psi

case minus control PSI change; the INC and EXC rows of an event have opposite signs

Examples

# # Multiple files from one rMATS run with the KO samples given as --b2
# input <- data.frame(
#   path = c('/path/A3SS.MATS.JC.txt', '/path2/A5SS.MATS.JC.txt'),
#   grp1 = c("WT","WT"),
#   grp2 = c("KO","KO"),
#   event_type = c("A3SS", "A5SS")
# )
# res <- get_rmats_post_di(input, case_group = "KO")

# Single rMATS table already loaded as df; group 1 (--b1) holds the case
# samples, the default
df <- data.frame(
  ID = 1L,
  GeneID = "ENSG00000182871",
  geneSymbol = "COL18A1",
  chr = "chr21",
  strand = "+",
  longExonStart_0base = 45505834L,
  longExonEnd = 45505966L,
  shortES = 45505837L,
  shortEE = 45505966L,
  flankingES = 45505357L,
  flankingEE = 45505431L,
  ID.2 = 2L,
  IJC_SAMPLE_1 = "4,1,0",
  SJC_SAMPLE_1 = "9,12,3",
  IJC_SAMPLE_2 = "0,4,5",
  SJC_SAMPLE_2 = "11,15,15",
  IncFormLen = 52L,
  SkipFormLen = 49L,
  PValue = 0.6967562,
  FDR = 1,
  IncLevel1 = "0.295,0.073,0.0",
  IncLevel2 = "0.0,0.201,0.239",
  IncLevelDifference = -0.024,
  stringsAsFactors = FALSE
)
res2 <- get_rmats_post_di(df, event_type = "A3SS")
print(res2)

Description

Runs the core SpliceImpactR pipeline from raw event table (or sample paths) through paired domain/PPI calls, then returns either a compact 'data.table' bundle ('data', 'res', 'hits_final') or an S4 [SpliceImpactResult].

Usage

get_splicing_impact(
  sample_frame = NULL,
  data = NULL,
  res = NULL,
  annotation_df = NULL,
  protein_feature_total = NULL,
  exon_features = NULL,
  ppi = NULL,
  source_data = c("rmats_hit", "hitindex", "rmats", "both"),
  event_types = c("ALE", "AFE", "MXE", "SE", "A3SS", "A5SS", "RI", "HFE", "HLE"),
  use = "JCEC",
  keep_annotated_first_last = FALSE,
  min_total_reads = 10L,
  minimum_proportion_containing_event = 0.5,
  terminal_fill = "event_max",
  cooks_cutoff = "Inf",
  adjust_method = "fdr",
  fdr_threshold = 0.05,
  delta_psi_threshold = 0.1,
  parallel_glm = TRUE,
  chunk_size_glm = 1000L,
  BPPARAM = BiocParallel::SerialParam(),
  chunk_size_match = 2000L,
  source_pairs = c("multi", "paired"),
  show_protein_domains = FALSE,
  return_class = c("data.table", "S4"),
  debug_steps = FALSE,
  metadata = list(),
  verbose = TRUE,
  matching = c("orf", "legacy"),
  matching_max_candidates = 100L,
  matching_fallback = TRUE
)

Arguments

Value

If 'return_class = "data.table"', returns a named list with 'data', 'res', and 'hits_final'. A supplied 'res' can be analyzed without raw data. If no events pass the cutoffs or no pairs remain, downstream results are zero-row tables with their usual columns, so plotting and enrichment helpers return empty results; unused reference/feature resources are not loaded. With 'matching = "orf"' (the default), an additional 'matching' list contains the candidate, ranking, event and unmatched diagnostics from [get_ranked_pairs()]. In S4 mode these are stored in 'metadata$matching_diagnostics'.

If 'return_class = "S4"', returns a [SpliceImpactResult] containing 'raw_events', 'di_events', and 'paired_hits' slots.

Examples

ex <- load_example_data(
  c("sample_frame", "annotation_df", "protein_feature_total", "ppi")
)
out <- get_splicing_impact(
  sample_frame = ex$sample_frame,
  annotation_df = ex$annotation_df,
  protein_feature_total = ex$protein_feature_total,
  ppi = ex$ppi,
  source_data = "rmats",
  event_types = c("SE"),
  use = "JCEC",
  parallel_glm = FALSE,
  BPPARAM = BiocParallel::SerialParam(),
  verbose = FALSE
)
print(names(out))

Description

SpliceImpactR accepts user-supplied splicing data in the same structure produced by get_rmats_hit(). Each event must be represented at the event-form and sample level.

**Event representation** - Each splicing event must be split into two forms: - INC: the inclusion isoform - EXC: the exclusion isoform - Alternative first/last exon events (AFE/ALE) or more abstract events may instead provide a single SITE form.

**Coordinates** - inc column: genomic coordinates included in the given form - exc column: genomic coordinates excluded in the given form - Coordinates may be one or multiple ranges (e.g., "100-200" or "100-150;300-350") - User must supply at least an inc and if supplying an exc, accompany with an inc coord

**Counts and PSI** - inclusion_reads and exclusion_reads must be provided per form - psi must be provided per sample (range 0-1) If psi isn't given, it will be extracted through inclusion_reads/exclusion_reads

**Sample structure** - Each event must have INC and EXC rows (or just SITE) - Each event must have >1 sample per condition (e.g., case vs control) - Required sample annotations: sample, condition

**Required columns** event_id, event_type, form, gene_id, chr, strand, inc, exc, inclusion_reads, exclusion_reads, psi, sample, condition

**Defaults** - If event_type not supplied: filled as "unknown" - If source_file not supplied: filled with empty string

This format enables downstream functionality including PSI modeling, annotation integration, and protein consequence prediction.

Usage

get_user_data(df)

Arguments

Value

data.table with cols, detailed above: "event_id","event_type","form", "gene_id","chr", "strand", "inc","exc","inclusion_reads","exclusion_reads", "psi", "sample", "condition","source_file" – designed to match get_rmats_hit output

Examples

example_df <- data.frame(
  event_id = rep("A3SS:1", 8),
  event_type = "A3SS",
  form = rep(c("INC","EXC"), each = 4),
  gene_id = "ENSG00000158286",
  chr = "chrX",
  strand = "-",
  inc = c(rep("149608626-149608834",4), rep("149608626-149608829",4)),
  exc = c(rep("",4), rep("149608830-149608834",4)),
  inclusion_reads = c(30,32,29,31, 2,3,4,3),
  exclusion_reads = c(1,1,2,1, 28,27,26,30),
  sample = c("S1","S2","S3","S4","S1","S2","S3","S4"),
  condition = rep(c("case","case","control","control"), 2),
  stringsAsFactors = FALSE
)
example_df$psi < example_df$inclusion_reads / example_df$exclusion_reads
user_data <- get_user_data(example_df)
print(user_data)

Description

This function converts user-supplied event results into the internal SpliceImpactR DI format. It accepts per-event statistics and ensures each splicing event contains valid inclusion/exclusion structure.

## Requirements **Event representation** - Each splicing event must be split into two forms: - INC: the inclusion isoform - EXC: the exclusion isoform - Alternative first/last exon events (AFE/ALE) or more abstract events may instead provide a single SITE form.

**Coordinates** - inc column: genomic coordinates included in the given form - exc column: genomic coordinates excluded in the given form - Coordinates may be one or multiple ranges (e.g., "100-200" or "100-150;300-350")

- User **must supply 'event_id'** (unique ID per splicing event, event_id = event_type:x for form != SITE. event_id = gene:event_type for form = SITE. Look at example output from get_differential_inclusion using test data for more examples) - Each event must be either: - **INC + EXC** forms (paired isoforms), or - **SITE** (single isoform) - Required columns: 'gene_id, chr, strand, inc, exc, form, event_id'

## Behavior - Does **not** generate event IDs: user must provide them - Constructs 'site_id = event_type|gene_id|chr|inc|exc|form' - 'delta_psi': - If missing: INC = +1, EXC = -1, SITE expands to +1 and -1 each to get all relevant comparisons - 'p.value', 'padj': - If missing: set to 0 - Sample counts ('n_samples, n_control, n_case'): set to 'NA' if missing - Other diagnostic fields ('cooks_max, n, n_used, mean_psi_ctrl, mean_psi_case') set to -1 if missing

## Validation - Throws error if: - 'event_id' is absent, or missing or empty in any row - Any event lacks **INC+EXC** or **SITE** - An event mixes SITE with INC/EXC - Errors name up to five affected rows or events.

## Output Returns a 'data.table' formatted like SpliceImpactR DI output. Ready for annotation, pairing, enrichment, and PPI analysis.

Usage

get_user_data_post_di(df)

Arguments

Value

A data.table formatted like SpliceImpactR DI output (get_differential_inclusion)

Examples

example_user_data <- data.frame(
  event_id = rep("A3SS:1", 8),
  event_type = "A3SS",
  gene_id = "ENSG00000158286",
  chr = "chrX",
  strand = "-",
  form = rep(c("INC","EXC"), each = 4),
  inc = c(
    rep("149608626-149608834", 4),
    rep("149608626-149608829", 4)
  ),
  exc = c(
    rep("", 4),
    rep("149608830-149608834", 4)
  ),
  inclusion_reads = c(30, 28, 25, 32,  2, 3, 4, 3),
  exclusion_reads = c(1, 2, 1, 1, 28, 27, 26, 30),
  sample = c("S1","S2","S3","S4","S1","S2","S3","S4"),
  condition = rep(c("case","case","control","control"), 2),
  stringsAsFactors = FALSE
)

# compute psi if missing, just for demo

post_di_user_data <- get_user_data_post_di(example_user_data)
print(post_di_user_data)

Description

Converts an arbitrary differential inclusion result table into the standardized column format expected by SpliceImpactR functions.

Usage

import_di_table(
  df,
  colmap = list(gene_id = "gene_id", chr = "chr", strand = "strand", inc = "inc", exc =
    "exc", delta_psi = "delta_psi", pvalue = "p.value", event_type = NULL),
  default_event_type = "SITE",
  adjust_method = "fdr",
  add_chr_prefix = FALSE
)

Arguments

Details

This function provides a uniform interface for importing external DI results (e.g. from rMATS, MAJIQ, or SUPPA2) so they can be compared or plotted alongside SpliceImpactR outputs. Strands must be '+' or '-'; other values (such as '*' or '.') are an error naming the affected events, since guessing a strand would mismatch minus-strand events.

Value

A standardized 'data.table' with columns: 'site_id', 'event_type', 'gene_id', 'chr', 'strand', 'inc', 'exc', 'delta_psi', 'p.value', 'padj', 'form', and 'event_id'. Missing sample counts are represented by 'NA'. Without an event ID, each site is its own event.

Examples

df <- data.frame(
  gene_id = "ENSG00000280071",
  chr = "7",
  strand = "+",
  inc = "1940088-1940549",
  exc = "",
  delta_psi = 0.25,
  p.value = 0.01
)
di_std <- import_di_table(df)
head(di_std)

Description

Provides a comprehensive visualization and summary of event classification outcomes (frame shifts, rescues, matches, etc.), alignment quality, and domain change prevalence across alternative splicing event types. Integrates results from ['compare_sequence_frame()'] with pre-filter event tables to display pre- vs post-filter usage, event coordination, and gene/domain overlaps.

Usage

integrated_event_summary(hits, pre_filter_hits)

Arguments

Details

The function integrates multiple layers of event-level characterization:

• Event-type composition and class proportions

• Protein alignment quality distribution

• Domain-change prevalence (case/control/both)

• Event-type coordination heatmap (Jaccard similarity across genes)

• Relative event retention pre- vs post-filtering

Value

A named list with:

'summaries'

List containing per-type tables: by_type, class_counts, score_summary, domain_prevalence, and relative_use.

'plot'

A multi-panel ['patchwork'] composite summarizing classification, alignment, domain changes, and coordination.

With no pairs, the summaries are zero-row tables (relative_use still counts the input events) and 'plot' is a placeholder.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)
ppi <- get_ppi_interactions()             
hits_final <- get_ppi_switches(hits_domain, ppi, protein_feature_total)
int_summary <- integrated_event_summary(hits_final, res)
print(int_summary)

Description

Keeps all rows belonging to events where at least one isoform or site passes adjusted p-value and deltaPSI significance criteria. The 'site_significant' column marks the rows that pass themselves; the other rows are kept as the event's comparators. In a multi-site event (for example an AFE with several first exons) these include sites whose own change does not meet the thresholds.

Usage

keep_sig_pairs(
  DT,
  padj_thr = 0.05,
  dpsi_thr = 0.1,
  return_class = c("auto", "data.table", "S4")
)

Arguments

Value

A 'data.table' (or updated 'SpliceImpactResult' when 'return_class' resolves to S4) containing all rows from event pairs in which at least one row meets the significance criteria, with a logical 'site_significant' column.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
sig_di <- keep_sig_pairs(res)
print(sig_di)

Description

Returns commonly used example objects (sample manifest, test annotations, optional feature mappings, optional PPI table) so man-page examples can stay focused on the documented function rather than setup boilerplate.

Usage

load_example_data(
  include = c("sample_frame", "annotation_df"),
  biomaRt_databases = c("interpro"),
  test = TRUE
)

Arguments

Value

Named list containing the requested objects.

Examples

ex <- load_example_data(c("sample_frame", "annotation_df"))
sample_frame <- ex$sample_frame
annotation_df <- ex$annotation_df

ex2 <- load_example_data(c("annotation_df", "exon_features"))
exon_features <- ex2$exon_features
print(exon_features)

Description

Parses rMATS output (.MATS.JC.txt or .MATS.JCEC.txt) for multiple event types and returns unified event tables ready for downstream inclusion/exclusion processing.

Usage

load_rmats(
  paths,
  use = c("JC", "JCEC"),
  event_types = c("SE", "RI", "A5SS", "A3SS", "MXE")
)

Arguments

Value

A 'data.table' with unified rMATS event annotations, including columns:

• event_type, event_id, gene_id, chr, strand

• sample, condition (if applicable)

• delta_psi, pvalue, fdr (if present)

• inclusion/exclusion read counts

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
rmats <- load_rmats(sample_frame, use = "JCEC", event_types = c("MXE", "SE", "A3SS", "A5SS", "RI"))
print(rmats)

Description

Generates a multi-panel comparison summarizing global splicing characteristics (event counts, events-per-gene, PSI distributions) between experimental and control conditions, normalized by sequencing depth.

Usage

overview_spicing_comparison(
  events,
  sample_df,
  depth_norm = c("exon_files", "user-given"),
  event_type = "AFE",
  conditions = c(control = "control", experimental = "case"),
  minReads = 10,
  output_file = NULL
)

Arguments

Details

The output combines: - **Panel A:** Depth-normalized event counts per sample (Wilcoxon test) - **Panel B:** Mean events per gene per sample (Wilcoxon test) - **Panel D:** PSI cumulative distribution comparison (K-S test)

Size factors are computed internally using [.get_size_factors_from_exons()], enabling robust normalization.

Value

Invisibly returns a combined ['patchwork'] plot object.

See Also

[.get_size_factors_from_exons()], [getSizeFactors()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
ov <- overview_spicing_comparison(hit_index, sample_frame, 'exon_files')
print(ov)

Description

Visualizes alignment scores ('prot_pid' or 'dna_pid') and coding class composition of hits, showing both the categorical composition and histogram of alignment identity values.

Usage

plot_alignment_summary(
  hits,
  mode = c("protein", "transcript"),
  output_file = NULL
)

Arguments

Value

A composite 'ggplot' object (from 'ggpubr::ggarrange') showing stacked bar counts by coding class and histogram of alignment identity scores. With no pairs, a placeholder plot.

See Also

[plot_length_comparison()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
alignment_summary <- plot_alignment_summary(seq_compare)
print(alignment_summary)

Description

Generates a volcano plot of deltaPSI vs. -log10(FDR) highlighting significant events.

Usage

plot_di_volcano_dt(di, padj_thr = 0.05, dpsi_thr = 0.1)

Arguments

Details

Significant sites are colored in 'deeppink4'; nonsignificant sites are shown in light grey. Dashed and dotted lines indicate deltaPSI and FDR thresholds.

Value

A 'ggplot2' object showing differential inclusion significance.

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
plot_di_volcano_dt(res)

Description

Visualizes the top enriched protein domains based on the number of events contributing to each domain's enrichment, optionally coloring by -log10 adjusted p-value.

Usage

plot_enriched_domains_counts(enriched_domains, top_n = 25)

Arguments

Details

The function expects the output of [enrich_domains_hypergeo()], typically a 'data.table' or 'data.frame' containing 'domain_id', 'events', and optionally 'padj' and 'OR'.

Each bar corresponds to one domain, with height proportional to the number of unique 'event_id's contributing to that domain. Bars are ordered by ascending adjusted p-value ('padj'), and colored by -log10(padj) if available. When no 'padj' column is present, the bars are shown in a uniform fill color.

Value

A 'ggplot' object showing bars of domain counts colored by enrichment significance.

See Also

[enrich_domains_hypergeo()], [enrich_by_event()], [enrich_by_db()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)
bg <- get_background(source = "annotated",
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
enriched_domains <- enrich_domains_hypergeo(hits_domain, bg, db_filter = 'interpro')
plot_enriched_domains_counts(enriched_domains, top_n = 20)

Description

Generates a multi-panel summary comparing isoform lengths between inclusion (INC) and exclusion (EXC) events for either protein or transcript modes.

Usage

plot_length_comparison(
  hits,
  phenotypes = c(control = "control", experimental = "case"),
  mode = c("protein", "transcript"),
  output_file = NULL
)

Arguments

Details

Produces three coordinated panels:

1. Paired boxplot showing INC vs EXC lengths per event with Wilcoxon paired test annotation.

2. Density plot of delta length (INC - EXC).

3. Barplot summarizing protein-coding categories.

Value

A composite 'ggplot' object (assembled with 'patchwork') showing paired boxplots, delta-length density, and protein-coding class distribution.

See Also

[plot_alignment_summary()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)
proximal_output <- plot_length_comparison(seq_compare)
print(proximal_output)

Description

Visualizes the frequency and magnitude of gained/lost PPIs per event, using a dual-panel layout: - left: proportion of events with any PPI change - right: histograms of CASE and CONTROL partner counts (non-zero only)

Usage

plot_ppi_summary(
  df,
  bins = 30,
  palette = c(no = "grey80", yes = "deeppink4", CASE = "#2b8cbe", CONTROL = "#e34a33"),
  output_file = NULL,
  width = 9,
  height = 4.8
)

Arguments

Value

A 'ggplot' object combining two panels (using 'patchwork').

See Also

[ppi_switches_for_hits()]

Examples

ex <- load_example_data("sample_frame")
sample_frame <- ex$sample_frame
hit_index <- get_hitindex(sample_frame)
res <- get_differential_inclusion(hit_index)
annotation_df <- load_example_data("annotation_df")$annotation_df
matched <- get_matched_events_chunked(res, annotation_df$annotations, chunk_size = 2000)
x_seq <- attach_sequences(matched, annotation_df$sequences)
pairs <- get_pairs(x_seq, source="multi")
seq_compare <-compare_sequence_frame(pairs, annotation_df$annotations)

annotation_df <- get_annotation(load = 'test')
interpro_features <- get_protein_features(c("interpro"), annotation_df$annotations, timeout = 600, test = TRUE)
protein_feature_total <- get_comprehensive_annotations(list(interpro_features))

exon_features <- get_exon_features(annotation_df$annotations, protein_feature_total)

hits_domain <- get_domains(seq_compare, exon_features)

bg <- get_background(source = "hit_index",
                     input = sample_frame,
                     annotations = annotation_df$annotations,
                     protein_features = protein_feature_total)
ppi <- get_ppi_interactions()             
hits_final <- get_ppi_switches(hits_domain, ppi, protein_feature_total)
ppi_plot <- plot_ppi_summary(hits_final)
print(ppi_plot)