Path B — Standalone tools (reads → quant)

Run each stage yourself when you want transparency, have only a few samples, or can't use Nextflow/containers. Do not assume numerical equivalence to Path A: inspect all assignment, filtering and quantification options. Quantify every sample identically.

Install (bioconda): conda create -n rnaseq -c conda-forge -c bioconda --strict-channel-priority fastqc fastp trim-galore "star=2.7.11b" "salmon=2.8.0" subread multiqc rseqc.

These commands are illustrative except for the Salmon 2.8.0 index/quant flags, which were exercised on a tiny synthetic paired-end fixture. Pin the solved environment. Salmon 2.x uses a different index format from C++ Salmon <=1.12.0; rebuild indexes and use the same binary/reference for every sample. See the current migration guide.

0. Reference data

You need, for your organism and a pinned annotation release:

• genome FASTA (genome.fa) and matching annotation GTF (annotation.gtf) — for STAR/featureCounts.
• transcriptome FASTA (transcripts.fa, cDNA) — for Salmon.

Fetch genome, transcript FASTA and GTF from one explicit Ensembl/GENCODE release (or use the gget skill to locate that release). Keep genome and GTF from the same release.

1. QC raw reads — FastQC

mkdir -p qc/raw
fastqc -t 8 -o qc/raw reads/*.fastq.gz

Inspect per-base quality, adapter content, duplication, and over-represented sequences. Interpretation/thresholds: design-and-qc.md.

2. Trim — fastp (recommended) or Trim Galore

fastp is fast and emits a JSON/HTML report MultiQC understands:

mkdir -p trimmed
fastp \
  -i reads/s1_R1.fastq.gz -I reads/s1_R2.fastq.gz \
  -o trimmed/s1_R1.fq.gz  -O trimmed/s1_R2.fq.gz \
  --detect_adapter_for_pe --qualified_quality_phred 20 --length_required 36 \
  --thread 4 --json qc/s1.fastp.json --html qc/s1.fastp.html

Trim Galore (wraps Cutadapt + FastQC; auto-detects adapters):

trim_galore --paired --cores 4 --fastqc -o trimmed reads/s1_R1.fastq.gz reads/s1_R2.fastq.gz

Aggressive quality trimming is usually unnecessary for modern data and for STAR (which soft-clips); adapter removal is the main goal. Re-run FastQC on trimmed reads to confirm.

3a. STAR — genome alignment + gene counts

Build the index once per genome+annotation+read-length. --sjdbOverhang = read length − 1 (100 is a safe default). Human needs ~30 GB RAM.

mkdir -p star_index
STAR --runMode genomeGenerate --runThreadN 12 \
  --genomeDir star_index \
  --genomeFastaFiles genome.fa \
  --sjdbGTFfile annotation.gtf \
  --sjdbOverhang 100

Align each sample, asking STAR to also count reads per gene:

mkdir -p star
STAR --runThreadN 12 --genomeDir star_index \
  --readFilesIn trimmed/s1_R1.fq.gz trimmed/s1_R2.fq.gz --readFilesCommand gunzip -c \
  --outSAMtype BAM SortedByCoordinate \
  --quantMode GeneCounts \
  --outFileNamePrefix star/s1.

--quantMode GeneCounts writes s1.ReadsPerGene.out.tab (4 columns, see strandedness below). The sorted BAM is useful for QC (RSeQC, IGV) and for featureCounts.

3b. Salmon — decoy-aware selective alignment

Build a decoy-aware index (genome as decoy) so reads from unannotated/genomic regions don't get miscounted against transcripts:

# 1. Decoys = all genome sequence names; gentrome = transcriptome THEN genome (order matters)
rg '^>' genome.fa | sed 's/^>//; s/ .*//' > decoys.txt
cat transcripts.fa genome.fa | gzip > gentrome.fa.gz

# 2. Index (k=31 works for reads >=75 bp; use smaller k for shorter reads)
salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_index -k 31 -p 12

Quantify each sample (-l A auto-detects library type/strandedness; enable bias correction):

salmon quant -i salmon_index -l A \
  -1 trimmed/s1_R1.fq.gz -2 trimmed/s1_R2.fq.gz \
  --gcBias --seqBias -p 8 \
  -o quant/s1

Selective alignment is the default; --validateMappings is an accepted no-op in 2.8.0. For GENCODE bar-delimited transcript FASTA names, add --gencode at indexing so IDs match the GTF. Record any version stripping separately. Do not use this length-corrected bridge for 3′ tag assays.

Each quant/<sample>/quant.sf is transcript-level; aggregate to gene level with scripts/build_counts_matrix.py --from salmon (needs a tx2gene map — see counts-and-handoff.md). Always check the reported mapping rate (quant/<sample>/logs/salmon_quant.log).

3c. featureCounts — counts from a STAR BAM (alternative to STAR GeneCounts)

mkdir -p counts
# Strandedness: 0 unstranded, 1 forward, 2 reverse; verify before choosing.
featureCounts -T 8 -p --countReadPairs \
  -a annotation.gtf -g gene_id \
  -s 2 \
  -o counts/featurecounts.txt \
  star/s1.Aligned.sortedByCoord.out.bam star/s2.Aligned.sortedByCoord.out.bam

Pass all sample BAMs at once to get one matrix. Parse it with scripts/build_counts_matrix.py --from featurecounts.

Strandedness — get this right

The wrong setting can discard most assigned reads; there is no fixed loss fraction. Determine it once, then apply consistently:

• Salmon -l A auto-detects and reports the library type in quant/<sample>/lib_format_counts.json (ISR = reverse-stranded paired, ISF = forward, IU = unstranded; IS is not a valid unstranded code).
• Or run RSeQC infer_experiment.py -r genes.bed -i s1.bam on a STAR BAM.

Map the result to each tool:

Illumina TruSeq Stranded mRNA is reverse (-s 2, STAR col 4). When in doubt, let Salmon auto-detect and match the others to it.

4. Aggregate QC — MultiQC

multiqc qc/ star/ quant/ counts/ -o qc/multiqc

MultiQC collates FastQC, fastp/Trim Galore, STAR, Salmon, and featureCounts logs into one report — your QC narrative for the methods section. Then build the counts matrix (counts-and-handoff.md) and hand off to pydeseq2.

Official references reviewed

• STAR manual: GeneCounts columns, index and alignment flags.
• Salmon 2.8.0 and its installed CLI help: index/quant flags and migration.
• featureCounts CLI source: -p declares paired reads; --countReadPairs counts fragments. --fraction produces fractional counts and is not supported by this bridge.
• fastp, Trim Galore, FastQC, MultiQC.