Bio Workflows Rnaseq To De

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.

• CLI + R: fastp -> (salmon | STAR + featureCounts) -> tximport -> DESeq2/edgeR/limma-voom -> lfcShrink -> VST/volcano

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.

1. 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.
2. 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.
3. 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.
4. 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.

The canonical order and why

Each step assumes the previous; two reorderings silently produce wrong results.

1. QC/trim before quantification — adapter/quality tails corrupt pseudo-mapping and duplicate structure.
2. 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.
3. Import via tximport — the tx->gene collapse happens HERE, carrying the length offset (order-trap: summing NumReads biases genes under isoform shift).
4. 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.
5. DESeq() on raw counts with batch in the design — size factors, dispersion, Wald/LRT.
6. 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.
7. 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.