Bio Applied Isoform Analysis

Isoform Analysis with Long Reads

When to Use

• Aligning ONT direct-RNA/cDNA or PacBio IsoSeq (HiFi) reads to a genome with splice-aware alignment.
• Discovering novel transcript isoforms (novel exons, splice sites, exon combinations) from long-read RNA-seq.
• Quantifying transcript-level (not just gene-level) expression from full-length reads.
• Testing whether a gene's isoform usage proportions shift between two conditions (differential isoform/transcript usage, DIU/DTU) — distinct from a simple gene-level DE test.
• Investigating oncogenic splice variants (e.g., MET exon 14 skipping) or novel tissue-specific isoforms absent from GENCODE/RefSeq.

Version Compatibility

• minimap2 ≥2.26, samtools ≥1.19
• R ≥4.3, Bioconductor ≥3.18: bambu ≥3.4, DRIMSeq ≥1.30
• Python ≥3.10: pandas ≥2.0, numpy ≥1.26, statsmodels ≥0.14, matplotlib ≥3.8

Prerequisites

• Tools installed: minimap2, samtools, R with Bioconductor packages bambu and DRIMSeq (BiocManager::install(c("bambu","DRIMSeq"))).
• A reference genome FASTA and a GTF annotation (e.g., GENCODE) matching the same genome build as the reads.
• Basic-quality-checked long reads (see long-read-sequencing skill) — chimeric/adapter-contaminated reads should be removed first.
• Familiarity with compositional-data testing (Dirichlet-multinomial) helps interpret DRIMSeq output.

Technology Comparison

Isoform Categories (bambu)

Splice-Aware Alignment

Goal: map full-length long reads to the genome so introns are represented as CIGAR N operations that bambu can parse into exon chains.

Approach: use Minimap2's splice preset (splice:hq for PacBio HiFi), suppress secondary alignments so each read maps once, and sort/index with samtools.

# ONT cDNA
minimap2 -ax splice --secondary=no -C5 --cs \
    hg38.fa cdna_reads.fastq.gz \
    | samtools sort -o cdna_aligned.bam -@ 8
samtools index cdna_aligned.bam

# PacBio HiFi / IsoSeq (splice:hq for high-accuracy reads)
minimap2 -ax splice:hq --secondary=no -C5 \
    hg38.fa isoseq_reads.fastq.gz \
    | samtools sort -o isoseq_aligned.bam -@ 8
samtools index isoseq_aligned.bam

Key flags: -ax splice = long-read spliced alignment preset; --secondary=no = one alignment per read; -C5 = extra cost for non-canonical splice sites (GT-AG = 0, others penalized); --cs = output the alignment-difference string.

Isoform Discovery & Quantification — bambu (R)

Goal: assign reads to known and novel transcript models across all samples jointly, and produce a transcript-level count matrix.

Approach: run bambu() once across all BAMs (multi-sample mode keeps isoform models consistent between conditions), controlling novel-isoform sensitivity with NDR.

library(bambu)

annotations <- prepareAnnotations('gencode.v44.annotation.gtf')

se <- bambu(
  reads       = c('ctrl_rep1.bam', 'ctrl_rep2.bam', 'treat_rep1.bam', 'treat_rep2.bam'),
  annotations = annotations,
  genome      = 'hg38.fa',
  NDR         = 0.1   # novel discovery rate; lower = stricter (0.1 = 90% confidence)
)

writeBambuOutput(se, path = 'bambu_output/')
# counts_transcript.txt:     transcript-level counts (rows=transcripts, cols=samples)
# counts_gene.txt:           gene-level counts
# extended_annotations.gtf:  all transcripts including novel ones

Differential Isoform Usage — DRIMSeq (R)

Goal: test whether the relative proportions of a gene's isoforms change between conditions, independent of the gene's total expression level.

Approach: DRIMSeq models isoform counts per gene as a Dirichlet-multinomial distribution (counts are compositional — they sum to the gene total) and runs a likelihood-ratio test per gene, then per transcript.

library(DRIMSeq)

counts_tx <- read.table('bambu_output/counts_transcript.txt', header = TRUE)
sample_info <- data.frame(
  sample_id = colnames(counts_tx)[-c(1, 2)],
  condition = c('Control', 'Control', 'Treatment', 'Treatment')
)

d <- dmDSdata(counts = counts_tx, samples = sample_info)

# Filter: require >=2 isoforms per gene and adequate expression
d <- dmFilter(d,
  min_samps_gene_expr    = 2,
  min_samps_feature_expr = 2,
  min_gene_expr          = 10,
  min_feature_expr       = 5
)

d <- dmPrecision(d)
d <- dmFit(d)
d <- dmTest(d, coef = 'conditionTreatment')

res_gene <- results(d, level = 'gene')     # omnibus test
res_tx   <- results(d, level = 'feature')  # per-transcript test
sig      <- res_gene[res_gene$adj_pvalue < 0.05, ]
write.csv(res_gene, 'drimseq_gene_results.csv', row.names = FALSE)

QC and Downstream Summary — Python

Goal: summarize bambu isoform categories and apply multiple-testing correction to DRIMSeq results outside R (e.g., for reporting or plotting pipelines).

Approach: load the bambu/DRIMSeq CSV exports with pandas and reuse statsmodels for FDR control — never re-derive BH correction by hand.

import pandas as pd
from statsmodels.stats.multitest import multipletests


def summarize_bambu_categories(counts_tx_path: str) -> pd.DataFrame:
    """Load a bambu extended annotation/category table and report the
    fraction of transcripts in each isoform category (annotated, novel_*, intergenic).

    Parameters
    ----------
    counts_tx_path : str
        Path to a table with at least a 'category' column (e.g. derived from
        bambu's extended_annotations.gtf via gffutils, or a custom export).
    """
    df = pd.read_csv(counts_tx_path, sep=None, engine='python')
    if 'category' not in df.columns:
        raise ValueError("expected a 'category' column (annotated/novel_in_catalog/...)")
    breakdown = df['category'].value_counts(normalize=True).rename('fraction').to_frame()
    breakdown['n'] = df['category'].value_counts()
    return breakdown


def diu_fdr_summary(drimseq_gene_csv: str, alpha: float = 0.05) -> pd.DataFrame:
    """Apply Benjamini-Hochberg FDR correction to DRIMSeq gene-level p-values
    and return genes passing the given significance threshold.

    Parameters
    ----------
    drimseq_gene_csv : str
        Path to the CSV written by `write.csv(res_gene, ...)` in R.
    alpha : float
        FDR threshold for calling a gene differentially-used (default 0.05).
    """
    res = pd.read_csv(drimseq_gene_csv)
    pvals = res['pvalue'].fillna(1.0).to_numpy()
    _, adj_pvals, _, _ = multipletests(pvals, alpha=alpha, method='fdr_bh')
    res['adj_pvalue'] = adj_pvals
    return res.loc[res['adj_pvalue'] < alpha].sort_values('adj_pvalue')

Pitfalls

• Coordinate systems: BED uses 0-based half-open; VCF/GFF use 1-based inclusive — mixing them causes off-by-one errors.
• PCR-cDNA amplification bias: PCR amplification distorts isoform frequency ratios — avoid when input is sufficient for PCR-free protocols.
• NDR threshold (bambu): Default NDR=1 accepts all novel isoforms. Use NDR=0.1 for strict filtering (90% confidence a transcript is genuine). Too lenient produces many false positives.
• Isoform quantification is compositional: Counts per gene sum to a total — use Dirichlet-multinomial models (DRIMSeq), not simple t-tests, for differential isoform usage.
• Multi-sample bambu runs: always run bambu() once across all samples/conditions together so isoform models stay consistent; running samples separately then merging breaks DRIMSeq's per-gene compositional assumptions.
• Batch effects: check for batch confounding before interpreting biological signal.
• Multiple testing: apply FDR correction (Benjamini-Hochberg) when testing thousands of genes/transcripts simultaneously.

See Also

• long-read-sequencing — upstream basecalling, QC, and read processing for ONT/PacBio data.
• bio-applied-rna-seq-analysis — short-read RNA-seq and gene-level DE for comparison.
• ai-science-splicing-models — deep-learning splice-site and isoform prediction models.
• scrna-seq-analysis — single-cell isoform/splicing analysis when data is cell-resolved.