Orchestrates the end-to-end bulk RNA-seq differential-expression pipeline from FASTQ to an annotated DE gene table, chaining fastp QC/trim, Salmon (decoy-aware) or STAR+featureCounts quantification, tximport gene-level collapse, DESeq2/edgeR/limma-voom testing, apeglm shrinkage, and VST-based visualization. Use when committing the reference release and gene-ID namespace once for the whole run, sequencing steps in the defensible order (tximport before DE, raw counts into the model, VST only for viz/clustering), choosing alignment-free vs align-then-count and the DE engine, setting strandedness correctly, keeping batch in the design instead of correcting-then-testing, or handing the signed ranking statistic to downstream enrichment. Hands mechanism to the component skills; not a re-teach of any single step.
Orchestrates the end-to-end bulk RNA-seq differential-expression pipeline from FASTQ to an annotated DE gene table, chaining fastp QC/trim, Salmon (decoy-aware) or STAR+featureCounts quantification, tximport gene-level collapse, DESeq2/edgeR/limma-voom testing, apeglm shrinkage, and VST-based visualization. Use when committing the reference release and gene-ID namespace once for the whole run, sequencing steps in the defensible order (tximport before DE, raw counts into the model, VST only for viz/clustering), choosing alignment-free vs align-then-count and the DE engine, setting strandedness correctly, keeping batch in the design instead of correcting-then-testing, or handing the signed ranking statistic to downstream enrichment. Hands mechanism to the component skills; not a re-teach of any single step.
[{"after_qc":"Q30 >80%, adapter content <5% (RNA has a lower quality floor than DNA)"},{"after_quant":"Mapping rate >70%, >10M reads mapped, flat gene-body coverage, low rRNA/intronic"},{"after_import":"tx2gene release matches the Salmon index; ID conversion loses few transcripts"},{"after_de":"Dispersion trend sane, PCA separates condition not batch, no Cook's outliers"}]
Version Compatibility
Reference examples tested with: DESeq2 1.42+, tximport 1.30+, apeglm 1.24+, STAR 2.7.11+, Salmon 1.10+, Subread/featureCounts 2.0.2+ (--countReadPairs added in 2.0.2), fastp 0.23+, ggplot2 3.5+ (kallisto 0.50+ as a Salmon alternative)
Before using code patterns, verify installed versions match. If versions differ:
R: packageVersion('<pkg>') then ?function_name to verify parameters
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.
Note: Salmon selective alignment is default since 1.0 (the historical --validateMappings is now a no-op); DESeqDataSetFromTximport carries the average-transcript-length offset automatically; lfcShrink(type='apeglm') requires coef to name a resultsNames(dds) coefficient and DROPS the stat column. Confirm these in-tool before quoting.
RNA-seq to Differential Expression Workflow
"Find differentially expressed genes from my RNA-seq FASTQ files" -> Chain QC/trim, decoy-aware quantification, tximport gene-level collapse, a count-based DE test, shrinkage, and visualization into one annotated DE table.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. Every step below cross-references the component skill that teaches its mechanism.
The governing principle
A bulk RNA-seq result is decided at four seams between steps, not inside any one tool.
The reference RELEASE + transcriptome/GTF pair is a pipeline-wide commitment made once at quantification and inherited by everything downstream. The transcriptome FASTA that builds the Salmon index and the GTF that builds the tx2gene map (and drives featureCounts) must be the SAME Ensembl/GENCODE release. Mixing an index built on release 104 with a tx2gene from release 110 silently drops renamed/removed transcripts — no error, just missing genes. This choice also fixes the gene-ID namespace (ENSG is the safe backbone; convert to symbol/Entrez only at the reporting/enrichment seam). Changing the release later forces re-quantification.
Raw counts flow forward; normalized/transformed values are terminal. Integer counts (or tximport count-scale output) are the ONLY valid input to DESeq2/edgeR/limma-voom. TPM/CPM are for within-sample ranking only; a VST/rlog matrix is for PCA, clustering, heatmaps, and ML — never for re-running a count-based test. Feeding the wrong scale across a join is the single most common silent corruption.
Collapse transcript->gene through tximport, not by summing counts. tximport carries the average-transcript-length offset that corrects for isoform-usage shifts; naively summing Salmon NumReads biases gene counts whenever isoform usage changes across conditions.
Batch belongs in the design, not "corrected" then tested. Put known batch in the formula (~ batch + condition). Running removeBatchEffect/ComBat and then testing on the corrected matrix exaggerates confidence (Nygaard 2016); the corrected matrix is for visualization/clustering/ML input only.
Pipeline map
FASTQ (paired)
| [1] QC & trim ----------------> fastp (read-qc/fastp-workflow)
v
| [2] Quantify -----------------> salmon (decoy-aware) (rna-quantification/alignment-free-quant)
v | OR STAR + featureCounts (need a BAM?) (read-alignment/star-alignment)
v ^-- commitment: reference RELEASE + tx/GTF pair, gene-ID namespace
| [3] Import & collapse tx->gene -> tximport (rna-quantification/tximport-workflow)
v ^-- carries the length offset; NEVER sum NumReads
| [4] Pre-DE QC ----------------> PCA / dispersion / outliers (rna-quantification/count-matrix-qc)
v
| [5] DE test ------------------> DESeq2 | edgeR-QL | limma-voom (differential-expression/deseq2-basics)
v ^-- RAW counts in; batch in the design, not corrected-then-tested
| [6] Shrink & extract ---------> lfcShrink(apeglm); pull Wald `stat` for ranking (differential-expression/de-results)
v
| [7] Visualize ----------------> VST heatmap/PCA, volcano (differential-expression/de-visualization)
v
Annotated DE table (ENSG + symbol + biotype, log2FC, stat, pvalue, padj, baseMean)
Reference, IDs, and quantification target: the made-once commitments
Decided before the first salmon quant; everything downstream inherits them. Mechanism lives in the component skills; the reasoning below is what a reviewer expects justified.
Commitment
Options
Consequence inherited downstream
Reference release
One Ensembl/GENCODE release for BOTH the transcriptome FASTA (index) and the GTF (tx2gene / featureCounts)
Any mismatch silently drops renamed transcripts; fixes DE row names and the pathway-DB key space
Gene-ID namespace
ENSG backbone (convert to symbol/Entrez only at reporting)
Symbol space is lossy (aliases, many-ENSG-one-symbol merges genes); stripping the ENSG .version with \..* also destroys the GENCODE _PAR_Y tag, collapsing chrY-PAR onto chrX (rna-quantification/tximport-workflow)
Quantification target
Gene-level DGE (countsFromAbundance="no", offset carried) vs transcript-level DTU
DTU needs a DIFFERENT import (txOut=TRUE + dtuScaledTPM); switching later is a re-import, not a filter — see workflows/splicing-pipeline
3'-tagged vs full-length
3'-tagged (QuantSeq/bulk-10x): countsFromAbundance="no", no length offset
Length-bias correction does not apply to 3'-tagged libraries
The canonical order and why
Each step assumes the previous; two reorderings silently produce wrong results.
QC/trim before quantification — adapter/quality tails corrupt pseudo-mapping and duplicate structure.
Quantify to the committed reference — Salmon decoy-aware (genome as decoy) so intron/pseudogene reads are not misassigned to transcripts; STAR only when a genome BAM is also needed downstream.
Import via tximport — the tx->gene collapse happens HERE, carrying the length offset (order-trap: summing NumReads biases genes under isoform shift).
Pre-filter low-count genes (rowSums(counts) >= 10) — this is a speed/memory step, NOT the FDR filter. Order-trap: it does not replace the baseMean independent filtering that results() applies at the FDR step; filterByExpr(y, design) (edgeR) is the design-aware version and must run once BEFORE dispersion, never after.
DESeq() on raw counts with batch in the design — size factors, dispersion, Wald/LRT.
results() then lfcShrink() — independent filtering happens inside results() on baseMean; shrink LFC for effect sizes/ranking, but p-values stay from the unshrunken test. Order-trap: apeglm/ashr objects DROP the stat column — pull the Wald stat from unshrunk results() if a signed ranking metric is needed for GSEA.
Visualize on VST (heatmaps/PCA); volcano uses shrunken LFC + unshrunken p.
Choosing the quantifier and the DE engine
Pipeline-level selection only; mechanism lives in the component skills.
Fork
Lean toward
Hand off to
Alignment-free (Salmon/kallisto) vs align-then-count (STAR+featureCounts)
Salmon for gene-level DGE (decoy-aware, GC/seq-bias correction, no BAM); STAR when a genome BAM is also needed (splicing, coverage, novel junctions, variants)
limma-voom when library sizes vary >3x or outliers dominate; edgeR-QL for tight finite-sample type-I control; DESeq2 for the apeglm/downstream ecosystem (70-90% top-gene overlap on well-designed data)
no (gene DGE via DESeqDataSetFromTximport) / lengthScaledTPM (DGE when the tool can't take offsets) / dtuScaledTPM+txOut (DTU)
rna-quantification/tximport-workflow
Strandedness -s
Confirm, never assume: STAR ReadsPerGene.out.tab cols 3 vs 4, or RSeQC infer_experiment.py; dUTP/TruSeq is reverse (-s 2)
read-qc/rnaseq-qc
Primary path: Salmon + tximport + DESeq2
Goal: turn trimmed FASTQ into a shrunken, annotated gene-level DE table.
Approach: build a decoy-aware index once, quantify each sample, collapse to genes via tximport (release-matched tx2gene), test raw counts with batch in the design, shrink for ranking. Full runnable script: examples/salmon_deseq2_workflow.R.
# Index once: decoy-aware (genome as decoy) so intron/pseudogene reads are not misassigned
grep "^>" genome.fa | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat transcriptome.fa genome.fa > gentrome.fa
salmon index -t gentrome.fa -d decoys.txt -i salmon_index -k 31 -p 8
# Quantify (selective alignment is default since 1.0; --gcBias/--seqBias correct known biases)
salmon quant -i salmon_index -l A -1 trimmed/${s}_R1.fq.gz -2 trimmed/${s}_R2.fq.gz \
-o quants/${s} --gcBias --seqBias -p 8
library(tximport); library(DESeq2)# tx2gene MUST come from the same release as the index (else renamed transcripts drop silently)
txi <- tximport(files, type ='salmon', tx2gene = tx2gene, ignoreTxVersion =TRUE)
dds <- DESeqDataSetFromTximport(txi, colData = coldata, design =~ batch + condition)# batch in design
dds <- dds[rowSums(counts(dds))>=10,]# speed filter, NOT the FDR filter
dds$condition <- relevel(dds$condition, ref ='control')
dds <- DESeq(dds)# RAW counts in
res <- lfcShrink(dds, coef ='condition_treated_vs_control', type ='apeglm')# ranking/effect size# For GSEA ranking, pull the Wald stat from the UNSHRUNK results (apeglm drops `stat`):
res_stat <- results(dds, name ='condition_treated_vs_control')$stat
Alternative path: STAR + featureCounts + DESeq2
Goal: produce a genome BAM (reused by splicing/coverage/variant steps) alongside gene counts.
Approach: align with --sjdbOverhang = readlen-1, count with the verified strandedness, then DESeqDataSetFromMatrix. Full script: examples/star_deseq2_workflow.sh.
STAR --runMode genomeGenerate --genomeDir star_index --genomeFastaFiles genome.fa \
--sjdbGTFfile genes.gtf --sjdbOverhang 149 --runThreadN 8 # 149 for 2x150, not a blanket 100
STAR --genomeDir star_index --readFilesIn trimmed/${s}_R1.fq.gz trimmed/${s}_R2.fq.gz \
--readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts \
--outFileNamePrefix aligned/${s}_ --runThreadN 8
# -s from ReadsPerGene.out.tab cols 3 vs 4 (or infer_experiment.py); -s 2 = dUTP/TruSeq reverse
featureCounts -T 8 -p --countReadPairs -s 2 -a genes.gtf -o counts.txt aligned/*_Aligned.sortedByCoord.out.bam
3' bias = degradation/oligo-dT; high intronic = pre-mRNA/gDNA; high intergenic = gDNA/annotation gap — all compromise DE BEFORE it runs (read-qc/rnaseq-qc)
Import
tx2gene release == index release; few transcripts dropped; report the ID-conversion rate (<0.85 => wrong ID type or organism)
alternative-splicing/isoform-switching - Transcript-level DTU when gene-level is not enough (splicing fork)
pathway-analysis/go-enrichment - Next step: functional enrichment (workflows/expression-to-pathways)
References
Soneson C, Love MI, Robinson MD (2015) Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4:1521. DOI 10.12688/f1000research.7563.1. (tximport; the tx->gene length-offset seam.)
Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550. DOI 10.1186/s13059-014-0550-8.
Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C (2017) Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods 14:417-419. DOI 10.1038/nmeth.4197.
Nygaard V, Rødland EA, Hovig E (2016) Methods that remove batch effects while retaining group differences may lead to exaggerated confidence in downstream analyses. Biostatistics 17:29-39. DOI 10.1093/biostatistics/kxv027. (batch belongs in the design.)
Ewels PA, Peltzer A, Fillinger S, et al (2020) The nf-core framework for community-curated bioinformatics pipelines. Nature Biotechnology 38:276-278. DOI 10.1038/s41587-020-0439-x. (nf-core/rnaseq: the reproducible reference orchestration.)
