RNA-seq Analysis — Pipeline Skill (Thin Scaffold)

Overview

Adds decision logic for choosing quantification paths (STAR vs Salmon), strandedness determination, and DESeq2 modeling pitfalls that LLMs routinely get wrong.

Usage

• Activate when choosing between alignment-based (STAR+featureCounts) vs alignment-free (Salmon) quantification
• Activate when setting up DESeq2 from either count source
• Activate when debugging unexpected low counts or weak DE signal

Core Concepts

Decision Logic

Need BAMs? (variant calling, fusions, novel transcripts)
├── YES → STAR two-pass + featureCounts
└── NO → Salmon (faster, lower RAM, bias-corrected)

Salmon output → DESeq2:
└── Use DESeqDataSetFromTximport (preserves length offset)
    └── NEVER round txi$counts and pass to DESeqDataSetFromMatrix

Filtering strategy (VQSR equivalent):
├── ≥30 samples → VQSR-like (default DESeq2 shrinkage works well)
└── <30 samples → still works, but pre-filter aggressively

Strandedness determination (MUST verify before counting):

Kit	featureCounts -s	Salmon -l
TruSeq Stranded / dUTP	2	ISR (PE) / SR (SE)
Ligation / forward	1	ISF / SF
Unstranded	0	IU / U

Verify with: infer_experiment.py (RSeQC) or Salmon's lib_format_counts.json.

Critical Parameters

Parameter	Value	When to change
--sjdbOverhang	readLength - 1 (default 100)	Non-100bp reads
--twopassMode Basic	Always for DE	Skip only for speed in QC runs
Salmon --gcBias --seqBias	Always on	Never skip for Illumina
Salmon --validateMappings	Always on	Enables selective alignment
featureCounts -s	Kit-dependent	ALWAYS verify empirically
DESeq2 pre-filter	rowSums(counts >= 10) >= smallest_group	Adjust 10 for shallow sequencing
lfcShrink type	"apeglm"	Use "ashr" for arbitrary contrasts
Design formula	~ batch + condition (interest LAST)	Always put interest term last

Common Mistakes

• Wrong: Guessing strandedness (-s) instead of verifying empirically Right: Check with infer_experiment.py or Salmon lib_format_counts.json before running featureCounts Why: Wrong strandedness causes 2–10× lower counts and weak DE signal

• Wrong: Passing normalized counts (TPM, FPKM, CPM) to DESeq2 Right: Always start from raw integer counts Why: Normalized values break the negative-binomial model

• Wrong: Rounding txi$counts and passing to DESeqDataSetFromMatrix Right: Use DESeqDataSetFromTximport(txi, ...) which preserves the length offset Why: Rounding discards per-gene length correction for transcript-length bias

• Wrong: Placing variable of interest first in design (~ condition + batch) Right: Put interest last: ~ batch + condition Why: DESeq2 tests the last term by default

• Wrong: Assuming DESeq2 reorders samples to match Right: Assert stopifnot(all(colnames(cts) == rownames(coldata))) before creating DESeq object Why: DESeq2 silently mislabels samples if order doesn't match

• Wrong: Skipping pre-filtering of low-count genes Right: Filter with rowSums(counts(dds) >= 10) >= smallest_group_size Why: Near-zero genes inflate runtime and dilute multiple-testing correction

• Wrong: Using STAR ReadsPerGene.out.tab column 2 for stranded libraries Right: Use column 3 (forward) or column 4 (reverse) matching your kit Why: Column 2 is unstranded; using it for stranded data halves effective counts

• Wrong: Including samples with <10M mapped reads or <70% mapping rate Right: Exclude or flag before fitting Why: Low-quality samples distort size factors and dispersion estimates

• Wrong: Modeling when batch is perfectly confounded with condition Right: Check table(coldata$batch, coldata$condition) first Why: No statistical model can separate perfectly correlated effects

Response Format

• Lead with the command or code the user needs — explain after
• Structure as: confirm inputs → working code → key parameters explained → gotchas
• One complete working example per task; do not show every alternative
• Keep code comments minimal and functional (what, not why-it-exists)
• Target: 50-100 lines of code with brief surrounding explanation