How to Run Isoform-Level Analysis on Your Own RNA-Seq Data

Oak
rna-seq isoforms salmon splicing transcriptomics
An iridescent engineered salmon on a dark studio stand, its body splitting into several translucent overlapping variants with differing fins and tails.

By the end of this guide you will have a transcript-level expression matrix for your own RNA-seq sample or samples, expressed in transcripts per million (TPM) and estimated counts, plus a per-gene table of isoform fractions showing which transcripts carry the expression for each gene. If you have several conditions or timepoints, you will also have a differential transcript usage (DTU) result: a list of genes where the mix of isoforms shifts even when total gene expression does not. What you need is a FASTQ file pair per sample (paired-end 100 or 150 bp reads, ideally 30 million or more read pairs for isoform work), a reference transcriptome and genome FASTA from GENCODE or Ensembl, roughly 16 GB of RAM and 8 cores, and R 4.3 or newer with Bioconductor. Everything below runs on a laptop for a handful of samples and on a single cloud VM for dozens.

A word on why isoform-level analysis is worth the extra effort. Gene-level counts collapse every transcript from a locus into one number, which discards the information that most human genes produce several mRNAs with different exons, promoters, and 3’ ends. Those transcripts often encode proteins with different domain content, localization, or stability, so a locus can look unchanged at the gene level while the functional output changes completely 1. The transcript-level view is also the layer that connects to protein: matching peptide evidence to a sample-specific set of expressed transcripts substantially improves isoform-level interpretation of proteomics data, because peptides that are ambiguous against the full reference become assignable once you know which isoforms are present 2.

1. Choose your reference and understand what it constrains

Short-read isoform quantification is a mapping problem against a fixed catalog of transcripts, which means your answer can only be as good as that catalog. Start with GENCODE for human, taking both the transcript FASTA and the matching genome FASTA for the same release, since mixing releases will silently break your annotation joins later.

REL=46
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.transcripts.fa.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/GRCh38.primary_assembly.genome.fa.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.annotation.gtf.gz

Decide now whether you want the full annotation (about 250,000 transcripts, including retained-intron and processed-transcript entries) or only the protein-coding and lncRNA set. We use the full set for discovery and then filter at the analysis stage, because pre-filtering the index makes reads from an excluded isoform pile onto its nearest included relative and inflate that estimate. The tradeoff is that full GENCODE contains many lowly supported transcripts whose estimates are noisy, so the filtering has to happen somewhere.

Be aware of the ceiling this imposes. Reference catalogs remain incomplete and biased toward well-studied loci, and annotation of even a single isoform of a common gene can lag for years, as the case of the ITGA6 X1X2 variant illustrates 3. If your question concerns a novel junction or an unannotated 3’ end, no short-read quantifier will find it, and you should skip ahead to step 7.

2. Build a decoy-aware Salmon index

We recommend Salmon with selective alignment over alignment-free quantification without decoys, and over a full alignment pipeline (STAR plus RSEM) for most personal-scale work. Selective alignment scores candidate mappings against the actual sequence rather than trusting k-mer compatibility, and adding the genome as a decoy prevents reads from intronic or intergenic origin from being forced onto a transcript. The cost is index build time and about 15 GB of RAM.

# concatenate transcriptome + genome, record decoy names
gunzip -c GRCh38.primary_assembly.genome.fa.gz | grep "^>" | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v46.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index \
  -t gentrome.fa.gz \
  -d decoys.txt \
  -i salmon_idx_v46 \
  -k 31 \
  --gencode \
  -p 8

Two flags matter. --gencode strips the pipe-delimited GENCODE headers down to the bare ENST identifier, which saves you a string-parsing step in R. -k 31 is the default and is right for 100 bp or longer reads; drop to -k 23 only if your reads are 50 bp or shorter, and expect worse isoform resolution either way.

3. Quantify with bias correction and bootstraps

Run Salmon per sample. The flags below are the ones we would not omit.

salmon quant \
  -i salmon_idx_v46 \
  -l A \
  -1 sample01_R1.fastq.gz \
  -2 sample01_R2.fastq.gz \
  --validateMappings \
  --seqBias \
  --gcBias \
  --posBias \
  --numGibbsSamples 50 \
  --threads 8 \
  -o quant/sample01

-l A lets Salmon infer library type, which you should then check in lib_format_counts.json; a standard dUTP stranded kit reports ISR, and anything else means you should look at your protocol before trusting the numbers. --seqBias and --gcBias correct random-hexamer priming bias and fragment GC bias, both of which distort isoform ratios when isoforms differ in GC content. --numGibbsSamples 50 produces posterior samples of the transcript counts, which is how you get per-transcript uncertainty. That uncertainty is not optional information at the isoform level: two transcripts differing by a single 30 bp cassette exon will have nearly identical read compatibility, and the split between them may be almost arbitrary. Gibbs samples tell you when that is happening.

Check the mapping rate in logs/salmon_quant.log. For a good human poly(A) library against decoy-aware GENCODE, expect 70 to 85 percent. Below 60 percent, stop and diagnose (see Common problems).

4. Import into R and inspect isoform fractions

Use tximport to read the quantifications. Import twice: once at transcript level for isoform work, once summarized to gene level for context.

library(tximport); library(tximeta); library(dplyr); library(readr)
library(GenomicFeatures)

files <- file.path("quant", samples$id, "quant.sf")
names(files) <- samples$id

txdb <- makeTxDbFromGFF("gencode.v46.annotation.gtf.gz")
tx2gene <- AnnotationDbi::select(txdb, keys = keys(txdb, "TXNAME"),
                                keytype = "TXNAME", columns = "GENEID")

# transcript level, counts scaled for DTU-compatible use
txi.tx <- tximport(files, type = "salmon", txOut = TRUE,
                   countsFromAbundance = "scaledTPM")

# gene level, length-corrected offsets for DGE
txi.gene <- tximport(files, type = "salmon", tx2gene = tx2gene,
                     countsFromAbundance = "no")

The countsFromAbundance = "scaledTPM" choice at transcript level is deliberate. Raw transcript counts are not comparable across samples when effective length changes, and the downstream DTU models do not accept length offsets, so scaled TPM is the correct input. For gene-level differential expression you want "no" and let DESeq2 use the average transcript length offsets.

Before any statistics, look at the isoform fractions for a handful of genes you care about. This is often where the interesting result already is.

tpm <- txi.tx$abundance
frac <- as.data.frame(tpm) |>
  tibble::rownames_to_column("tx") |>
  left_join(tx2gene, by = c("tx" = "TXNAME")) |>
  group_by(GENEID) |>
  mutate(across(where(is.numeric), ~ .x / sum(.x))) |>
  ungroup()

Filter to genes with total TPM above roughly 5 before reading fractions, because a fraction computed on a gene at 0.4 TPM is noise. Transcript-level interpretation frameworks such as Isoformic formalize this step, pairing differential expression with isoform-fraction changes so you can distinguish a gene that is simply up from a gene that has switched its dominant transcript 4.

5. Test differential transcript usage

Differential expression and differential usage answer different questions. DTU asks whether the proportions of a gene’s transcripts differ between conditions, and it will find switches that are invisible to gene-level testing. The DRIMSeq-plus-DEXSeq route, with stageR for two-stage testing, is the approach we use for bulk short-read data because it models the Dirichlet-multinomial nature of proportions and controls error at both the gene and transcript levels.

library(DRIMSeq); library(DEXSeq); library(stageR)

cts <- txi.tx$counts
cts <- cts[rowSums(cts) > 0, ]

counts_df <- data.frame(
  gene_id    = tx2gene$GENEID[match(rownames(cts), tx2gene$TXNAME)],
  feature_id = rownames(cts),
  cts, check.names = FALSE)

d <- dmDSdata(counts = counts_df, samples = samples)

n <- nrow(samples); n.small <- min(table(samples$condition))
d <- dmFilter(d,
  min_samps_feature_expr = n.small, min_feature_expr = 10,
  min_samps_feature_prop = n.small, min_feature_prop = 0.10,
  min_samps_gene_expr    = n,       min_gene_expr    = 10)

design <- model.matrix(~ condition, data = DRIMSeq::samples(d))

dxd <- DEXSeqDataSet(
  countData  = round(as.matrix(counts(d)[, -c(1,2)])),
  sampleData = DRIMSeq::samples(d),
  design     = ~ sample + exon + condition:exon,
  featureID  = counts(d)$feature_id,
  groupID    = counts(d)$gene_id)

dxd <- estimateSizeFactors(dxd)
dxd <- estimateDispersions(dxd, quiet = TRUE)
dxd <- testForDEU(dxd, reducedModel = ~ sample + exon)
dxr <- DEXSeqResults(dxd, independentFiltering = FALSE)

The filter above is the part people skip and then regret. Requiring a transcript to hold at least 10 percent of its gene’s expression in the smaller group removes the long tail of near-unidentifiable isoforms whose p-values are driven by quantification uncertainty rather than biology. Then stage the testing: screen genes with a per-gene q-value, and only confirm individual transcripts within genes that pass.

pScreen <- perGeneQValue(dxr)
pConfirm <- matrix(dxr$pvalue, ncol = 1,
                   dimnames = list(dxr$featureID, "transcript"))
stageRObj <- stageRTx(pScreen = pScreen, pConfirmation = pConfirm,
                      pScreenAdjusted = TRUE,
                      tx2gene = as.data.frame(counts(d)[, c("feature_id","gene_id")]))
stageRObj <- stageWiseAdjustment(stageRObj, method = "dtu", alpha = 0.05)
res <- getAdjustedPValues(stageRObj, order = TRUE, onlyReferenceTx = FALSE)

With a single sample and no comparison group, DTU testing does not apply, and the right output is the isoform fraction table from step 4 plus the Gibbs-derived uncertainty. That is still a useful personal baseline: it tells you which transcript of each gene you express, and a second timepoint months later gives you the comparison.

6. Ask what the switch does to the protein

A DTU hit is a statistical statement about mRNA proportions. Turning it into something interpretable means asking what the two transcripts encode, and this is where most analyses stop too early. Pull the coding sequences and compare domain content.

library(IsoformSwitchAnalyzeR)
# after aSwitchList construction from salmon quant dirs + GTF + transcript FASTA
exampleSwitchList <- isoformSwitchAnalysisPart1(
  switchAnalyzeRlist = aSwitchList,
  pathToOutput = "isar_out",
  outputSequences = TRUE,   # writes AA + nt FASTA for CPC2/Pfam/SignalP/IUPred2A
  prepareForWebServers = TRUE)

Run the exported amino acid FASTA through Pfam for domains, SignalP for signal peptides, and CPC2 for coding potential, then re-import the results and let the package annotate consequences: domain loss, intron retention, alternative 3’ UTR, or a premature termination codon predicted to trigger nonsense-mediated decay. This is the step that maps a proportion change onto a plausible functional difference.

Keep the limits in view. RNA and protein isoform evidence agree less often than one would hope, because splicing changes are buffered by translation and degradation, and current mass spectrometry covers only a fraction of isoform-distinguishing peptides 5. Peptides unique to a given isoform are frequently absent from tryptic digests altogether, which is why matched transcript and protein evidence is treated as complementary rather than confirmatory 1. Multi-omics studies that combine splicing calls with protein and clinical data in cancer cohorts make the same point: splicing alterations are widespread and organized, and their consequences become interpretable only when several layers are read together 67. Nextflow pipelines for this kind of integration now exist and are worth using if you have more than two data types to reconcile 8.

If a result concerns a gene relevant to a clinical question, this is where to bring in a physician or genetic counselor. Isoform-level research findings from a personal dataset are not clinical results, and interpreting them for health decisions requires a clinician working with validated assays.

7. Decide whether you need long reads

Short reads cannot tell you which exons co-occur on the same molecule beyond the fragment length, so isoform structure is inferred rather than observed. A benchmark of thirteen short- and long-read isoform detection methods found wide disagreement in recovered structures, with the largest gaps at complex loci and low-expression transcripts. If your question is about full-length structure at a specific locus, add PacBio Iso-Seq or Oxford Nanopore cDNA sequencing. Integrating Iso-Seq structures with short-read quantification is the strongest combination available today: the long reads define the catalog including alternative promoters and polyadenylation sites, and the short reads supply the depth for accurate quantification against it 9.

Practically, we would sequence 2 to 3 million Iso-Seq reads for one tissue or blood sample, run isoseq3 cluster then pigeon classify against GENCODE to get a sample-specific GTF, append any novel high-confidence transcripts to the reference FASTA, and rebuild the Salmon index from step 2. Everything downstream is unchanged.

Common problems

A mapping rate below 60 percent almost always has one of four causes. Check for adapter and poly(A) contamination with FastQC and trim if needed. Check that you did not index a cDNA-only FASTA while feeding it total RNA data heavy in pre-mRNA, in which case the decoy genome should have absorbed those reads and its absence is the bug. Check for rRNA carryover, visible as a handful of transcripts holding 30 percent or more of TPM. Finally, confirm species and release: a mouse library against a human index maps at around 20 percent and produces entirely plausible-looking output.

Transcripts with wildly unstable estimates across Gibbs samples are the second recurring issue. Compute the coefficient of variation per transcript from the posterior samples in aux_info/bootstraps and flag anything above roughly 0.5 as unquantifiable in that library rather than reporting it. Two transcripts differing only in a short internal exon or a few bases of 5’ end will often both be flagged, and the correct action is to report their sum as one unit.

Retained-intron transcripts absorbing large fractions are usually a library artifact rather than biology. High intronic signal from incomplete poly(A) selection or degraded RNA inflates those entries. Check the intronic read fraction with a tool such as RSeQC or Qualimap, and if it exceeds 25 to 30 percent, treat retained-intron isoform calls as suspect for that sample.

Finally, resist reading a DTU hit as a mechanism. The test says proportions moved. Whether that reflects a splicing factor change, an alternative promoter, a shift in polyadenylation, or differential decay of one isoform requires the annotation work in step 6 and, for a definitive answer, long reads across the locus. Reporting the switch with its predicted consequence and its uncertainty is the complete and defensible output.

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. 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 ↩ ↩2

  2. Becky C. Carlyle, Robert R. Kitchen, Jing Zhang, et al. Isoform-Level Interpretation of High-Throughput Proteomics Data Enabled by Deep Integration with RNA-seq. Journal of Proteome Research, 2018. https://doi.org/10.1021/acs.jproteome.8b00310 ↩

  3. Ximena Aixa Castro Naser, Alessandro Cestaro, Silvio C. E. Tosatto, et al. Integrative Multi-Omics Characterization and Structural Insights into the Poorly Annotated Integrin ITGA6 X1X2 Isoform in Mammals. Genes, 2025. https://doi.org/10.3390/genes16101134 ↩

  4. Izabela Mamede, Lucio R Queiroz, Carlos Mata-Machado, et al. Isoformic: a workflow for transcript-level RNA-seq interpretation. NAR Genomics and Bioinformatics, 2025. https://doi.org/10.1093/nargab/lqaf176 ↩

  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 ↩

  7. Timothy I. Shaw, Bi Zhao, Yuxin Li, et al. Multi-omics approach to identifying isoform variants as therapeutic targets in cancer patients. Frontiers in Oncology, 2022. https://doi.org/10.3389/fonc.2022.1051487 ↩

  8. Bianka Alexandra Pasat, Eleftherios Pilalis, Katarzyna Mnich, et al. MultiOmicsIntegrator: a nextflow pipeline for integrated omics analyses. Bioinformatics Advances, 2024. https://doi.org/10.1093/bioadv/vbae175 ↩

  9. Haimei Wen, Wei Chen, Yu Chen, et al. Integrative analysis of Iso-Seq and RNA-seq reveals dynamic changes of alternative promoter, alternative splicing and alternative polyadenylation during Angiotensin II-induced senescence in rat primary aortic endothelial cells. Frontiers in Genetics, 2023. https://doi.org/10.3389/fgene.2023.1064624 ↩