1. Overview

Regardless of whether your splicing quantification comes from VAST-TOOLS/VastDB or rMATS, the pipeline has the same three stages:

1. Splicing-genome construction: take a table of alternative-splicing events (exon-skipping or intron-retention) and, for each event, define a fixed set of genomic sub-regions around it (flanking exon ends, flanking intron ends, the alternative exon itself, etc; see the CHARM paper's Methods, "Splicing-genome creation"). This produces one row per event with genomic coordinates for each sub-region.
2. Binding metagenome assembly: intersect a table of eCLIP peaks (RBP, chromosome, start, end) against every region of every event, recording where and how much of each region a given RBP's peaks cover. This produces a long-format, per-RBP-per-event overlap table.
3. RNA binding maps: reshape the overlap table into a wide "position matrix" (metagenomic coordinates 1–1000 for exon-skipping, 1–500 for intron-retention), split events by direction of splicing change (increased/decreased/maintained inclusion), and run a chi-squared test at every position to test whether the RBP's binding density differs between event classes. This is what eCLIPSE_full() / eCLIPSE_full_IR() do, and is exactly what the CHARM paper's figures show.

Two independent event sources feed stage 1: VastDB/VAST-TOOLS (event IDs like HsaEX0024923, HsaINT...) and rMATS (SE.MATS.JC.txt, numeric event IDs). Both converge on the same stage-2/stage-3 scripts. That convergence is the whole point of rMATS compatibility, and is what makes this pipeline useful if you don't want to use VAST-TOOLS.

2. Environment & Dependencies

R

Packages used throughout the notebook: data.table, tidyr, dplyr, stringr, ggplot2, ggpubr, scales. The betAS (ASfile) step (Section 6) additionally needs betAS itself.

Python 3

Used by eCLIPSE_Exon.py / eCLIPSE_Intron.py: pandas, numpy.

External tools

VAST-TOOLS/VastDB path only: vast-tools (align/combine). rMATS path only: rMATS-turbo.

3. Path A: VAST-TOOLS / VastDB Splicing Genome VastDB

3.1 Align and combine

Scripts: vt_align.sh, vt_combine.sh

vt_align.sh is a thin wrapper around vast-tools align. Run it from the directory containing your paired-end FASTQs (named <sample>_1.paired.fastq.gz / <sample>_2.paired.fastq.gz):

./vt_align.sh /path/to/fastqs/

Then combine all aligned samples into one inclusion table:

./vt_combine.sh

This produces a vast-tools combine output (an INCLUSION_LEVELS_FULL-hg38-*.tab.gz file) — this is the vttable input to the next script.

3.2 Build the splicing genome.

Script: Eclip_preprocessing_vttools.R

This script reads the VAST-TOOLS inclusion table that you just generated, and writes the exon-skipping and intron-retention splicing-genome tables.

4. Path B: rMATS Splicing Genome rMATS

4.1 Run rMATS-turbo yourself

Scripts: rMATSp1.sh, rMATSp2.sh

rMATSp1.sh indexes every BAM file in a directory with samtools index (rMATS-turbo needs the .bai files present):

./rMATSp1.sh /path/to/bams/

rMATSp2.sh then runs rmats.py itself, comparing two sample groups (one BAM path per line. in this example samples_stressed.txt/samples_unstressed.txt):

rmats.py --b2 samples_unstressed.txt \
      --b1 samples_stressed.txt \
      --gtf ~/Genomes/Human/gencode.v42.primary_assembly.basic.annotation.gtf \
      --bi ~/Genomes/Human/STARindex/ \
      -t single \
      --readLength 76 \
      --nthread 4 \
      --od final_output_rMATS \
      --tmp intermediate_output_rMATS > rMATS_GSE171009.txt

A single rmats.py run writes one *.MATS.JC.txt file per event type into --od: SE.MATS.JC.txt (exon skipping), RI.MATS.JC.txt (intron retention), A5SS.MATS.JC.txt, A3SS.MATS.JC.txt, and MXE.MATS.JC.txt.

4.2 Build the splicing genome

Script: Eclip_preprocessing_rmats.R

This script converts SE.MATS.JC.txt into the same 20-column splicing-genome schema that Eclip_preprocessing_vttools.R produces, so it can be fed straight into eCLIPSE_Exon.py instead of the VastDB path.

5. Binding Metagenome Assembly

Scripts: eCLIPSE_Exon.py, eCLIPSE_Intron.py

These scripts intersect your eCLIP peaks against the splicing-genome regions from Section 3/Section 4, per RBP per event, and write a long-format overlap table.

# Exon-skipping (works for both the VastDB and rMATS splicing-genome tables)
python3 eCLIPSE_Exon.py \
  --events Final_SmallExon_Table.txt \
  --peaks  all_data_combined_38.txt \
  --output SmallExon_events.txt

# Intron retention (VastDB path only. No rMATS IR table exists yet)
python3 eCLIPSE_Intron.py \
  --events Final_IR_Table.txt \
  --peaks  all_data_combined_38.txt \
  --output Intron_events.txt

--peaks expects a space-delimited table (as produced by R's write.table(..., row.names = F)) with at minimum the columns Chrom, StartCord, EndCord, RBP (the combined eCLIP peaks file, all_data_combined_38.txt, already has exactly this format. See Section 9).

6. Producing the betAS (ASfile) Table

This is the step between splicing quantification (Section 3/Section 4) and the RNA binding maps (Section 7/Section 8), where we turn a VAST-TOOLS inclusion table into the ASfile that eCLIPSE_full()/eCLIPSE_full_IR() expect (an Event.ID column plus a dPSI column, Section 8.1). It uses the betAS package to go from raw per-sample PSI values to a statistically-supported per-event delta-PSI.

6.1 Loading and filtering the dataset

The example below frames this as a typical ENCODE-style eCLIP experiment: a shRNA knockdown of the target RBP compared against control samples.

library(betAS)
library(data.table)

dataset <- getDataset(pathTables = "INCLUSION_LEVELS_FULL-hg38-25-v251.tab.gz", tool = "vast-tools")
dataset <- getEvents(dataset, tool = "vast-tools")

# Keep exon-skipping (S) and intron-retention (IR) events
dataset_filtered <- filterEvents(dataset, types = c("S", "IR"), N = 0)
dataset_filtered$PSI  <- na.omit(dataset_filtered$PSI)
dataset_filtered$Qual <- dataset_filtered$Qual[dataset_filtered$Qual$EVENT %in% dataset_filtered$PSI$EVENT, ]

# Two groups: control and shRNA-knockdown samples
samples_control <- c("control_rep1", "control_rep2")
samples_shRNA   <- c("shRNA_rep1", "shRNA_rep2")

dataset_control <- dataset_filtered
dataset_control[["Samples"]] <- samples_control
dataset_control[["PSI"]]  <- dataset_control$PSI[,  c(1:6, which(colnames(dataset_control$PSI)  %in% samples_control))]
dataset_control[["Qual"]] <- dataset_control$Qual[, c(1:6, which(colnames(dataset_control$Qual) %in% paste0(samples_control, ".Q")))]

dataset_shRNA <- dataset_filtered
dataset_shRNA[["Samples"]] <- samples_shRNA
dataset_shRNA[["PSI"]]  <- dataset_shRNA$PSI[,  c(1:6, which(colnames(dataset_shRNA$PSI)  %in% samples_shRNA))]
dataset_shRNA[["Qual"]] <- dataset_shRNA$Qual[, c(1:6, which(colnames(dataset_shRNA$Qual) %in% paste0(samples_shRNA, ".Q")))]

6.2 Filtering to informative events

There are two ways to filter dataset_control/dataset_shRNA down to the events worth testing. Pick one, not both, since both produce a same-named _filtered object.

Option A: betAS's own alternativeEvents(), which filters events to a PSI range:

dataset_control_filtered <- alternativeEvents(dataset_control, minPsi = 1, maxPsi = 99)
dataset_shRNA_filtered   <- alternativeEvents(dataset_shRNA,   minPsi = 1, maxPsi = 99)

Option B: the custom NewFilterFunction(), which is not part of the betAS package. It's defined in eCLIPSE-Supplementary.Rmd and filters differently: instead of a PSI range, it keeps only events that appear in a separately-precomputed event list, FullCassetteEvents.txt, that contains only events with substantial evidence for true alternative splicing. This is the version that is used throughout CHARM:

FullTrial <- fread("FullCassetteEvents.txt")

dataset_control_filtered <- NewFilterFunction(dataset_control, FullTrial)
dataset_shRNA_filtered   <- NewFilterFunction(dataset_shRNA,   FullTrial)

6.3 Computing dPSI and building the ASfile

cols_control <- convertCols(dataset_control_filtered$PSI, samples_control)
cols_shRNA   <- convertCols(dataset_shRNA_filtered$PSI,   samples_shRNA)

volcanoTable <- prepareTableVolcano(
  psitable       = dataset_control_filtered$PSI,   # control/shRNA filtered tables share the same events by construction
  qualtable      = dataset_control_filtered$Qual,
  npoints        = 500,
  colsA          = cols_control,
  colsB          = cols_shRNA,
  labA           = "control",
  labB           = "shRNA",
  basalColor     = "#89C0AE",
  interestColor  = "#E69A9C",
  maxDevTable    = maxDevSimulationN100,     # ships with the betAS package itself
  seed           = TRUE,
  CoverageWeight = FALSE
)

plotVolcano(betasTable = volcanoTable, labA = "control", labB = "shRNA",
            basalColor = "#89C0AE", interestColor = "#E69A9C")   # optional sanity-check plot

ASfile <- volcanoTable[, c("EVENT", "GENE", "deltapsi", "Pdiff")]
colnames(ASfile) <- c("Event.ID", "Gene", "dPSI", "Pdiff")
write.table(ASfile, "my_betAS_table.txt", row.names = FALSE, col.names = TRUE)

prepareTableVolcano() is betAS's core statistical step: it computes delta-PSI (deltapsi) between control and shRNA and a significance-like statistic (Pdiff), via a permutation/simulation approach.

7. Reshaping eCLIPSE_Exon.py Output into a Position Matrix

Notebook: Eclip_position_matrix.Rmd

This is the final pre-processing step of the binding metagenome. Do not take the long-format overlap table from Section 5 directly. You need a wide rnamapfile, with one row per (event, RBP), one column per metagenomic position (1–1000 for exon-skipping, 1–500 for intron retention), plus EVENTS and RBP columns. Converting an event's raw, per-region overlap coordinates into that shared metagenome coordinate system is exactly what this sextion does.

7.1 Parsing the raw overlap output

Near the top of the notebook, before any of the three path-specific sections: parse_overlap_columns(df, regions) turns each "(start,end)"-style column eCLIPSE_Exon.py/eCLIPSE_Intron.py produced (one per biological region, e.g. scaled_values_UPS_EX) into numeric ..._ST/..._EN columns, and drops any row where parsing failed for any region:

exon_regions <- c("UPS_EX", "UPS_IN_ST", "UPS_IN_EN", "AE_ST", "AE_EN",
                   "DOW_IN_ST", "DOW_IN_EN", "DOW_EX")
mapped_plus  <- parse_overlap_columns(subset(raw, STRAND == "+"), exon_regions)
mapped_minus <- parse_overlap_columns(subset(raw, STRAND == "-"), exon_regions)

7.2 One recipe, three paths

Each path's section in the notebook has the same four-part shape: Loading and parsing (read that path's raw overlap file, split by strand), Positive Stranded and Negative Stranded (build one 0/1/NA density block per region, cbind() them together), and Assembling the position matrix (rbind() both strands, attach EVENTS/RBP, saveRDS() the result). Within a "Stranded" section, the same shape of code repeats once per region: Pull that region's raw start/end/length columns; zero-initialise a width-column matrix, one row per event; loop over rows, skipping any row with missing data or a degenerate region, clamping the coordinate into the window, marking the covered columns 1, and padding the columns past the mapped length with NA; then label the columns with the metagenome positions this region occupies.

density_UPS_EX <- matrix(0, nrow = nrow(mapped_plus), ncol = 50)
for (row in seq_len(nrow(mapped_plus))) {
  st <- mapped_plus$scaled_values_UPS_EX_ST[row]
  en <- mapped_plus$scaled_values_UPS_EX_EN[row]
  if (is.na(st) || is.na(en)) next
  if ((en - st) <= 3) next   # degenerate: region too short to trust
  if (en > 50) en <- 50      # clamp into the window
  density_UPS_EX[row, st:en] <- 1
}
colnames(density_UPS_EX) <- as.character(1:50)

Every other region/strand block in the notebook follows this same shape, just with a different width and an NA-pad step added in:

7.3 Running it

Open Eclip_position_matrix.Rmd, point the placeholder path in that path's "Loading and parsing" chunk at your own overlap file, then run that path's chunks top to bottom: Loading and parsing → Positive Stranded → Negative Stranded → Assembling the position matrix. The last chunk of each path ends with a saveRDS() call:

# Exon Skipping — BigExon section, "Assembling the position matrix" chunk:
rnamapfile <- as.data.frame(rbind(plus_mat, minus_mat))
rnamapfile$EVENTS <- c(mapped_plus$EVENT, mapped_minus$EVENT)
rnamapfile$RBP    <- c(mapped_plus$RBP, mapped_minus$RBP)
saveRDS(rnamapfile, "eCLIPSE_position_matrix.rds")

8. RNA Binding Map Functions

Notebook: eCLIPSE_plotting_functions.Rmd

The RNA-binding-map functions themselves.

8.1 Calling the plotting function

plot <- eCLIPSE_full(
  rnamapfile = my_position_matrix,   # wide matrix, see Section 7
  ASfile     = my_splicing_table,    # needs columns Event.ID, dPSI
  rnaBP      = "RBFOX2",              # RBP name to plot, must match a value in rnamapfile$RBP
  metric     = "FDR",                # or "EffectSize"
  title      = "My experiment"
)
print(plot)

ASfile needs a Event.ID column (VastDB HsaEX.../HsaINT... IDs or rMATS numeric IDs, and a dPSI column (used to classify events as increased/decreased/maintained inclusion against PSIthreshold, default 0.05).

Use eCLIPSE_full() for exon-skipping data (VastDB or rMATS) and eCLIPSE_full_IR() for intron-retention data (VastDB only).

8.2 Running it

Open eCLIPSE_plotting_functions.Rmd, run the function-definition chunks (eCLIPSE_full(), eCLIPSE_full_IR()), then the "Running it" chunk with your own data:

rnamapfile <- readRDS("eCLIPSE_position_matrix.rds")  # from Section 7
ASfile     <- fread("my_betAS_table.txt")             # from Section 6

eCLIPSE_full(rnamapfile = rnamapfile, ASfile = ASfile, rnaBP = "RBM39",
             metric = "FDR", plot = TRUE, title = "control vs shRNA")

# Intron-retention path (VastDB only):
rnamapfile_ir <- readRDS("eCLIPSE_position_matrix_IR.rds")
eCLIPSE_full_IR(rnamapfile = rnamapfile_ir, ASfile = ASfile, rnaBP = "RBM39",
                 metric = "FDR", plot = TRUE, title = "control vs shRNA, intron retention")

8.3 Per-RBP binding heatmaps

Alternatively, you can run eCLIPSE for all RBPs for a single perturbation:

eCLIPSE_heatmap(rnamapfile, ASfile,
                rbps = c("SRSF1", "U2AF2"),
                metric = "EffectSize")   # or "FDR"

# Intron-retention path (VastDB only):
eCLIPSE_heatmap_IR(rnamapfile_ir, ASfile,
                    rbps = c("SRSF1", "U2AF2"),
                    metric = "EffectSize")   # or "FDR"

metric = "EffectSize" colours each tile by its chi-squared statistic, metric = "FDR" by -log10(FDR) both categories (decreased and increased inclusion) are drawn on one shared, symmetric colour scale, sized to the largest absolute value across both, so the same colour means the same thing in both panels rather than each panel scaling itself independently.

9. The Large Data Files (Downloadable Example Data)

Several files in the ECLIPSE folder are large precomputed datasets rather than code:

These are included as downloadable example data so you can try the pipeline (or just stages 5 and 7) without first needing your own eCLIP peaks or splicing quantification.