How to Analyze Alternative Splicing in Your Own RNA-Seq Data

By the end of this guide you will have four artifacts from a single person’s RNA-seq: a filtered splice junction table with read support per junction, per-event PSI (percent spliced in, the fraction of transcripts that include a given exon or junction) for annotated and novel events, transcript-level abundances with bootstrap uncertainty, and a shortlist of rare variants from your whole genome that a splice prediction model flags and that your own RNA-seq either supports or contradicts. You need paired-end RNA-seq at 75 bp or longer with roughly 50–100 million read pairs per sample (junction-level work is read-hungry in a way that gene-level expression is not), a machine with at least 64 GB of RAM and 500 GB of disk, GRCh38 with GENCODE annotation, and, for the last step, a VCF from whole-genome sequencing. Everything below runs on Linux with conda-installed tools.

A word on what splicing is, since the vocabulary matters for reading output. Human genes are transcribed as pre-mRNA containing exons (kept) and introns (removed). The spliceosome cuts at the 5’ donor and 3’ acceptor sites and ligates the exons. Alternative splicing is the use of different donor/acceptor combinations in the same gene, producing different mature transcripts. The event types every tool reports are exon skipping (SE), alternative 5’ splice site (A5SS), alternative 3’ splice site (A3SS), mutually exclusive exons (MXE), and intron retention (RI). Most human multi-exon genes are alternatively spliced, and splicing changes are a recurring feature of both normal tissue identity and disease biology ¹.

1. Build a splice-aware genome index

Use the primary assembly FASTA and a matching GENCODE GTF. Do not mix annotation versions between the index, the event caller, and the quantifier: junction coordinates will line up but transcript IDs will not.

STAR --runMode genomeGenerate --runThreadN 16 \
  --genomeDir star_gencode44_sjdb100 \
  --genomeFastaFiles GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile gencode.v44.primary_assembly.annotation.gtf \
  --sjdbOverhang 100

Set --sjdbOverhang to read length minus one. Index construction needs about 32 GB of RAM for the human genome and takes under an hour on 16 threads.

We use STAR rather than a lightweight pseudoaligner here because splicing analysis needs actual spliced alignments: you want to know how many reads crossed a specific donor-acceptor pair with a specific overhang on each side. Pseudoalignment gives transcript abundances but discards the per-junction evidence you will use to check individual events later. Run both, for different purposes.

2. Align with two-pass STAR and keep the junction attributes

STAR --genomeDir star_gencode44_sjdb100 \
  --readFilesIn sample_R1.fq.gz sample_R2.fq.gz --readFilesCommand zcat \
  --twopassMode Basic \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS nM NM MD jM jI MC \
  --outSAMstrandField intronMotif \
  --outFilterType BySJout \
  --alignSJoverhangMin 8 --alignSJDBoverhangMin 3 \
  --outFilterMismatchNoverReadLmax 0.04 \
  --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 \
  --quantMode GeneCounts TranscriptomeSAM \
  --runThreadN 16 --outFileNamePrefix sample_

Two-pass mode matters more for splicing than for expression. In the first pass STAR collects novel junctions, then re-indexes with them inserted, so reads that only marginally overlap a novel junction in pass one get properly spliced alignments in pass two. Without it you systematically undercount novel and low-usage junctions, which are exactly the ones worth looking at in one person’s data.

--outFilterType BySJout removes reads whose junctions did not survive STAR’s own filtering, which cuts a large share of spurious junctions before you ever see them. --alignSJoverhangMin 8 requires eight bases on the short side of a novel junction; lowering it inflates false positives from short repeats near splice sites. Keep jM and jI in the SAM attributes so you can trace individual reads back to junctions during manual review.

3. Extract and filter junctions

STAR writes sample_SJ.out.tab with columns for chromosome, intron start and end, strand, motif, annotation status, uniquely mapped reads, multi-mapped reads, and maximum overhang. That file is the raw material for most of what follows. For a cleaner, strand-aware version keyed to the BAM, use regtools:

regtools junctions extract -a 8 -m 50 -M 500000 -s XS \
  sample_Aligned.sortedByCoord.out.bam -o sample.junc

(-s takes RF, FR, or XS depending on library strandedness and regtools version; check regtools junctions extract -h for yours.)

Filter before interpreting anything. Our defaults: at least six uniquely mapping reads supporting the junction, maximum overhang of at least 12 bp, canonical GT-AG or GC-AG or AT-AC motif unless you have a specific reason to keep a non-canonical junction, and intron length between 50 and 500,000 bp. In a single 60-million-pair library you will typically get 150,000–250,000 junctions passing these filters, of which the large majority are annotated.

Once junctions are filtered, the simplest useful quantity is splice-site usage: for each donor site, the fraction of junction reads leaving that donor that go to each acceptor, and symmetrically for acceptors. This sidesteps the need to classify events into SE/A5SS/A3SS categories at all, which is an advantage because the classification schemes disagree with each other and force ambiguous events into boxes ². We compute both: splice-site usage for screening, event-level PSI for anything we plan to report.

4. Call and quantify events

For a single sample, or a small longitudinal series from one person, we use SplAdder to build the event catalog. It augments the annotation graph with junctions observed in your data, enumerates events of each type, and quantifies them, so novel exons and novel splice sites present in your transcriptome appear in the output rather than being silently discarded by an annotation-only quantifier ³.

spladder build -o spladder_out \
  -a gencode.v44.primary_assembly.annotation.gtf \
  -b sample1.bam,sample2.bam,sample3.bam \
  --confidence 3 --merge-strat merge_graphs \
  --output-txt --parallel 16

Confidence level 3 is the strictest setting and the one we default to for personal data, where a false novel exon costs more than a missed one. The output includes per-event PSI and the supporting read counts, which you should always carry alongside PSI: a PSI of 0.5 from 4 reads and from 400 reads are entirely different observations.

If you have two conditions with replicates (fasted versus fed, pre- and post-exercise, two timepoints months apart), add rMATS-turbo for statistics:

rmats.py --b1 cond1_bams.txt --b2 cond2_bams.txt \
  --gtf gencode.v44.primary_assembly.annotation.gtf \
  -t paired --readLength 100 --variable-read-length \
  --libType fr-firststrand --nthread 16 --od rmats_out --tmp rmats_tmp \
  --cstat 0.05 --paired-stats

--cstat 0.05 tests against a 5% PSI difference rather than against zero, which is the right null for personal data: a statistically detectable 1% shift in PSI is not something you can interpret. --paired-stats is appropriate when the same person contributes both conditions, which is the usual case here, and it substantially improves power by modeling the within-subject pairing.

5. Quantify transcripts independently

Event-level and transcript-level views fail in different ways, so run both and look for agreement. Build a Salmon index with decoy sequences so that reads from unannotated genomic regions are not forced onto transcripts:

grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f1 | sed 's/>//g' > decoys.txt
cat gencode.v44.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz
salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx -k 31 --gencode -p 16

salmon quant -i salmon_idx -l A -1 sample_R1.fq.gz -2 sample_R2.fq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numBootstraps 100 -p 16 -o salmon_sample

The 100 bootstrap replicates are the point. Isoform-level estimates for genes with many similar transcripts carry large inferential uncertainty, and the bootstrap distribution tells you which isoform calls are stable. Load them with tximport and inspect the coefficient of variation per transcript before trusting any isoform ratio.

Then compute PSI from the transcript TPMs with SUPPA2, which maps transcript abundances onto event definitions:

suppa.py generateEvents -i gencode.v44.primary_assembly.annotation.gtf \
  -o events -f ioe -e SE SS MX RI FL
suppa.py psiPerEvent -i events_SE_strict.ioe -e sample_tpm.txt -o sample_SE

Where SplAdder’s junction-based PSI and SUPPA2’s abundance-based PSI agree within about 0.1 for an event with good coverage, we treat the value as solid. Where they disagree, the usual cause is an unannotated isoform, and the junction-based number is closer to the truth.

6. Put your values in context

A PSI of 0.62 for an exon means nothing on its own. You need a distribution to compare against, and for one person the practical reference is GTEx v8, which provides per-tissue exon inclusion levels across hundreds of donors. Match the tissue (whole blood RNA-seq compares to whole blood, not to “median tissue”) and note the library difference: GTEx used polyA selection, so if your library is ribo-depleted total RNA your intron retention values will be systematically higher and are not comparable. Restrict cross-study comparisons to exon skipping and alternative splice site events, where the bias is smaller.

The more informative comparison is against yourself. Two or three samples of the same tissue from different months give you a within-person variance estimate, and events that move more than your own noise band are the ones worth attention. Splicing programs are also strongly context-dependent, so a shift may reflect cell composition rather than regulation: a change in neutrophil fraction in whole blood will move hundreds of PSI values at once. Estimate cell proportions (CIBERSORTx or a simple marker-gene approach) before attributing a splicing change to splicing.

7. Connect splicing to your genome

This is where a personal profile does something a public dataset cannot. Take rare variants from your WGS (gnomAD allele frequency below 0.1%) within 500 bp of any annotated splice site and score them with SpliceAI:

bcftools view -i 'INFO/AF<0.001' sample.vcf.gz | \
  bedtools intersect -a - -b splice_regions_500bp.bed -header > rare_splice_region.vcf

spliceai -I rare_splice_region.vcf -O rare_spliceai.vcf \
  -R GRCh38.primary_assembly.genome.fa -A grch38 -D 500

SpliceAI reports delta scores for acceptor gain, acceptor loss, donor gain, and donor loss. We triage at delta score 0.5 and above for candidates worth checking, and 0.2–0.5 as a lower-confidence tier. Then go back to your junction table and ask the direct question: at the predicted position, does your own RNA-seq show a junction that the reference does not have, or loss of a junction that should be there? A prediction with matching junction evidence in your own transcriptome is a far stronger observation than either alone. If the gene is heterozygous and you have phased reads, count junction reads by allele: true splice-disrupting variants usually produce an allele-specific effect rather than a uniform shift.

Interpretation of any such finding is a clinical question, not a computational one. Predicted splice-altering variants in genes with known disease associations should go to a clinical geneticist or genetic counselor working from a validated laboratory report, not from your pipeline output.

8. Look at the reads before you believe anything

Every event you intend to report gets a sashimi plot. The number of times a “novel exon” turns out to be a pileup over a repeat, or a “retained intron” turns out to be genomic DNA contamination, is high enough that this step is not optional.

ggsashimi.py -b bams.tsv -c chr17:7668000-7690000 \
  -g gencode.v44.primary_assembly.annotation.gtf \
  -M 10 -C 3 -O 3 -A median_j -P palette.txt -o event_plot

-M 10 hides junctions with fewer than ten supporting reads, which declutters the plot to the junctions you can reason about. What you want to see: coverage that drops to near zero across the skipped exon in the skipping isoform, junction arcs with counts consistent with your PSI, and no coverage in flanking intronic regions that would suggest pre-mRNA contamination.

9. Know the limit of the RNA evidence

A transcript-level PSI change does not entail a protein-level change. Translation efficiency, nonsense-mediated decay, and protein turnover all intervene, and isoform-specific peptides are frequently unobservable because the distinguishing region contains no tryptic peptide of usable mass. Bottom-up mass spectrometry recovers only a minority of predicted splice-variant-specific peptides, and top-down approaches that measure intact proteoforms cover a much narrower slice of the proteome ⁴. Where proteomics and transcript-level splicing data have been compared systematically, agreement is partial and depends heavily on the event type and on whether the affected region is peptide-detectable ⁵. If you have proteomics alongside your RNA-seq, treat a matching isoform-specific peptide as confirmation and its absence as uninformative rather than as refutation.

Common problems

Low junction counts across the board usually mean the library is 3’-biased or degraded. Check the RNA integrity number if you have it, and look at coverage uniformity along transcript bodies: a strong 3’ skew means junction PSI values in 5’ exons rest on almost no reads.

Intron retention values that look uniformly high are nearly always a library artifact, not biology. Ribo-depleted total RNA captures pre-mRNA; polyA-selected libraries do not, to the same degree. Compare intron retention only within one library type, and consider excluding the RI event class entirely if you are comparing to an external reference.

Novel junctions clustering in one genomic region point to a mapping problem: a segmental duplication, an unmasked repeat, or a paralog absorbing reads. Filter junctions overlapping ENCODE blacklist regions and check the multi-mapping read counts in the SJ.out.tab file.

Disagreement between rMATS and SUPPA2 on the same event is expected when the gene has an unannotated dominant isoform, since SUPPA2 can only distribute abundance across transcripts in the GTF. Rerun SplAdder with the junctions from your data included and see whether a novel isoform explains the gap.

Small PSI differences that reach statistical significance in a deep library are a trap. Fix a minimum effect size before you run the test (we use 0.1 absolute PSI) and require at least 20 junction-supporting reads across the event. Splicing shifts that matter biologically, including the ones seen in tumor transcriptomes profiled with matched genomic and proteomic data, are typically much larger than the detection threshold of a deep library ⁶.

Finally, batch effects from library prep dominate splicing metrics more than they dominate gene expression. If your longitudinal samples were prepared in different batches, you cannot separate a genuine splicing change from a prep change, and the reading of the data is that the comparison is confounded. Prepare longitudinal samples together when you can.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Yining Tao, Qi Zhang, Haoyu Wang, et al. Alternative splicing and related RNA binding proteins in human health and disease. Signal Transduction and Targeted Therapy, 2024. https://doi.org/10.1038/s41392-024-01734-2 ↩

2. Craig I Dent, Shilpi Singh, Sourav Mukherjee, et al. Quantifying splice-site usage: a simple yet powerful approach to analyze splicing. NAR Genomics and Bioinformatics, 2021. https://doi.org/10.1093/nargab/lqab041 ↩

3. André Kahles, Cheng Soon Ong, Yi Zhong, et al. SplAdder : identification, quantification and testing of alternative splicing events from RNA-Seq data. Bioinformatics, 2016. https://doi.org/10.1093/bioinformatics/btw076 ↩

4. Taojunfeng Su, Michael A.R. Hollas, Ryan T. Fellers, et al. Identification of Splice Variants and Isoforms in Transcriptomics and Proteomics. Annual Review of Biomedical Data Science, 2023. https://doi.org/10.1146/annurev-biodatasci-020722-044021 ↩

5. Marina Reixachs‐Solé, Eduardo Eyras. Uncovering the impacts of alternative splicing on the proteome with current omics techniques. WIREs RNA, 2022. https://doi.org/10.1002/wrna.1707 ↩

6. Quanyou Wu, Lin Feng, Yaru Wang, et al. Multi-omics analysis reveals RNA splicing alterations and their biological and clinical implications in lung adenocarcinoma. Signal Transduction and Targeted Therapy, 2022. https://doi.org/10.1038/s41392-022-01098-5 ↩