Bio Splicing Qc

Version Compatibility

Reference examples tested with: RSeQC 5.0+, STAR 2.7.11+, samtools 1.19+, pysam 0.22+, regtools 1.0+, maxentpy 0.0.1+, spliceai 1.3+, matplotlib 3.8+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Splicing-Specific Quality Control

Splicing analysis is more demanding than DGE on read length, depth, library prep, alignment strategy, and annotation choice. Failures in any of these silently bias PSI estimates and inflate novel-junction false positives. The decision sequence is: experimental design -> library prep -> alignment strategy -> annotation -> diagnostic metrics. Each layer's failure mode is distinct.

QC Layer Taxonomy

Decision Tree by Question

Experimental Design Audit (Before Sequencing)

STAR 2-Pass Alignment

Goal: Maximize novel-junction sensitivity for downstream AS analysis.

Approach: Run STAR once per sample to discover novel junctions (pass 1), merge novel junctions across cohort, then re-align with the augmented junction set (pass 2). Cohort-style 2-pass beats per-sample basic 2-pass for differential splicing because all samples use the same junction reference.

# Pass 1: per-sample
STAR --runMode alignReads \
    --runThreadN 8 \
    --genomeDir genome_index \
    --sjdbGTFfile gencode.v45.basic.gtf \
    --sjdbOverhang 149 \
    --readFilesIn sample_R1.fq.gz sample_R2.fq.gz \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --outFileNamePrefix pass1_${sample}_ \
    --outSJtype Standard \
    --outFilterMultimapNmax 20 \
    --alignSJoverhangMin 8 \
    --alignSJDBoverhangMin 1

# Cohort-style 2-pass: collect all SJ.out.tab from pass 1
cat pass1_*_SJ.out.tab | awk '$5 > 0 && $7 >= 3' | sort -u > cohort_novel_SJ.tab

# Pass 2: re-align with augmented junctions
STAR --runMode alignReads \
    --runThreadN 8 \
    --genomeDir genome_index \
    --sjdbGTFfile gencode.v45.basic.gtf \
    --sjdbFileChrStartEnd cohort_novel_SJ.tab \
    --sjdbOverhang 149 \
    --readFilesIn sample_R1.fq.gz sample_R2.fq.gz \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --outFileNamePrefix pass2_${sample}_ \
    --outSJtype Standard \
    --twopassMode None \
    --quantMode GeneCounts \
    --alignSJoverhangMin 8 \
    --alignSJDBoverhangMin 3

Per-sample 2-pass (--twopassMode Basic) is simpler but produces inconsistent junction sets across samples; for differential splicing the cohort-style version is preferred (Veeneman 2016 Bioinformatics).

The pass-1 filter awk '$5 > 0 && $7 >= 3' keeps junctions with strand info AND >=3 unique reads — adjust threshold to balance discovery vs noise.

Junction Saturation

Goal: Determine whether sequencing depth is sufficient for comprehensive splicing detection.

Approach: Run RSeQC junction saturation; check whether the discovery curve plateaus.

junction_saturation.py \
    -i sample.bam \
    -r gencode_v45.bed \
    -o sample_junc_sat \
    -m 50

import subprocess
import pandas as pd

samples = ['s1.bam', 's2.bam', 's3.bam']
for sample in samples:
    subprocess.run([
        'junction_saturation.py',
        '-i', sample,
        '-r', 'gencode_v45.bed',
        '-o', sample.replace('.bam', '_junc_sat')
    ], check=True)

The output *.junctionSaturation_plot.r plots known + novel junctions vs subsampled reads.

Plateau detection rule: if from 80% to 100% of reads, the junction count rises by <2%, consider it plateaued. Still rising means more sequencing would yield more junctions.

For AS analysis, plateau on the known junction curve is the requirement; novel-junction curves often don't plateau even at deep coverage (which is biologically informative — novel junctions are inherently rarer events).

Novel-vs-Known Junction Ratio

Goal: Detect annotation/mapping issues or biologically interesting cryptic splicing.

Approach: Classify junctions with RSeQC and compute the novel:known ratio.

junction_annotation.py -i sample.bam -r gencode_v45.bed -o sample_junc_annot

import pandas as pd

# RSeQC .junction.xls has a header: chrom, intron_st(0-based), intron_end(1-based), read_count, annotation
junc = pd.read_csv('sample_junc_annot.junction.xls', sep='\t')
total = junc['read_count'].sum()

by_class = junc.groupby('annotation')['read_count'].sum()
known_frac = by_class.get('annotated', 0) / total
novel_frac = (by_class.get('partial_novel', 0) + by_class.get('complete_novel', 0)) / total

print(f'known: {known_frac:.1%}, novel: {novel_frac:.1%}')

High novel-junction rate may be biology, not artifact:

• TDP-43 loss (ALS/FTD post-mortem brain): cryptic exon de-repression in UNC13A, STMN2, ATG4B (Brown 2022 Nature; Klim 2019 Nat Neurosci)
• SF3B1-mutant cancer (MDS, CLL, uveal melanoma): cryptic 3'ss ~10-30nt upstream of canonical (Darman 2015 Cell Rep)
• Non-model organism: GENCODE-grade annotation unavailable; novel junctions reflect annotation gaps not biology
• Microbial / viral contamination: reads aligning to host but with unusual junctions

If novel% >40%, drill down: check organism, check spliceosomal mutation status, check known disease signatures.

Junction Read Overhang and Coverage

Goal: Profile per-junction read counts and overhang distribution to identify weakly-supported events.

Approach: Parse CIGAR for N (intron) operations; tally per-junction reads and minimum exon overhangs.

import pysam
from collections import defaultdict

def junction_stats(bam_path):
    bam = pysam.AlignmentFile(bam_path, 'rb')
    counts = defaultdict(int)
    min_overhang = defaultdict(lambda: float('inf'))

    for read in bam.fetch():
        if read.is_unmapped or read.is_secondary:
            continue
        ref_pos = read.reference_start
        cumulative_query = 0
        cigar = read.cigartuples
        for i, (op, length) in enumerate(cigar):
            if op == 3:
                left_match = sum(l for o, l in cigar[:i] if o in (0, 7, 8))
                right_match = sum(l for o, l in cigar[i+1:] if o in (0, 7, 8))
                overhang = min(left_match, right_match)
                key = (read.reference_name, ref_pos, ref_pos + length)
                counts[key] += 1
                min_overhang[key] = min(min_overhang[key], overhang)
            if op in (0, 2, 3, 7, 8):
                ref_pos += length

    bam.close()
    return counts, dict(min_overhang)

counts, overhang = junction_stats('sample.bam')
print(f'total junctions: {len(counts)}')
print(f'>= 10 reads: {sum(1 for c in counts.values() if c >= 10)}')
print(f'overhang >= 8 nt: {sum(1 for k, c in counts.items() if overhang[k] >= 8)}')

Junction reads with overhang <8 nt are common false positives, especially for novel sites. Most callers default to >=8 nt anchor for this reason. Microexon-aware aligners use overhang as low as 6 nt with explicit configuration.

Splice Site Strength (MaxEntScan and SpliceAI)

Goal: Score donor and acceptor splice sites to flag weak / cryptic sites and to predict variant impact on splicing.

Approach: Use MaxEntScan (sequence information content) and SpliceAI (context-aware deep-learning) — they answer different questions.

from maxentpy.maxent import score5, score3

donor = 'CAGGTAAGT'
acceptor = 'TTTTTTTTTTTTTTTTTTTTCAG'
print(f"5'ss MaxEnt: {score5(donor):.2f}")
print(f"3'ss MaxEnt: {score3(acceptor):.2f}")

MaxEntScan vs SpliceAI:

• MaxEntScan scores sequence information content (intrinsic strength). Captures position-wise dependencies at the consensus.
• SpliceAI predicts in-vivo usage probability given full pre-mRNA context (10 kb window).
• A position with high MaxEnt but low SpliceAI is intrinsically strong but contextually silenced (chromatin, trans factors).
• A position with low MaxEnt but high SpliceAI is intrinsically weak but contextually used (enhancer-driven, e.g. weak donors stabilized by ESEs).
• Report both for variant interpretation; for variant impact see splice-variant-prediction.

Picard CollectRnaSeqMetrics and Gene-Body Coverage

Goal: Get integrated RNA-seq QC including intronic / exonic / intergenic mapping rates and gene-body coverage uniformity.

Approach: Run picard CollectRnaSeqMetrics for mapping distribution; RSeQC geneBody_coverage.py for 5'-3' bias.

picard CollectRnaSeqMetrics \
    I=sample.bam \
    O=sample.rna_metrics.txt \
    REF_FLAT=refFlat.txt \
    STRAND_SPECIFICITY=SECOND_READ_TRANSCRIPTION_STRAND \
    RIBOSOMAL_INTERVALS=rRNA_intervals.interval_list

# Strandedness conversion (foot-gun):
# Reverse-stranded (Illumina TruSeq Stranded; NEB Ultra II Directional — both dUTP):
#   rMATS  --libType fr-firststrand
#   featureCounts -s 2
#   Picard STRAND_SPECIFICITY=SECOND_READ_TRANSCRIPTION_STRAND
# Forward-stranded (Lexogen QuantSeq FWD, certain ligation-based kits):
#   rMATS  --libType fr-secondstrand
#   featureCounts -s 1
#   Picard STRAND_SPECIFICITY=FIRST_READ_TRANSCRIPTION_STRAND
# STAR has no library-strand flag; pass --outSAMstrandField intronMotif
# (works for any library) so downstream tools can read XS tags.

geneBody_coverage.py \
    -i sample.bam \
    -r gencode_v45.bed \
    -o sample_geneBody

3' bias (degraded RNA) directly reduces splicing-event detection because junction reads scatter across the gene body; with 3' bias they concentrate near the 3' end and miss CDS junctions.

Strandedness Verification

infer_experiment.py -i sample.bam -r gencode_v45.bed -s 200000

Output reports the fraction of reads consistent with each library type:

Wrong strand setting halves usable junction reads — always verify before quantification. RSeQC infer_experiment.py is fast and authoritative.

Annotation Choice

Comprehensive captures more biology but inflates DTU multiple-testing burden and includes annotation noise. For event-level (rMATS) AS, basic is usually adequate; for transcript-level DTU (DRIMSeq, satuRn), comprehensive may be necessary to capture rare isoforms.