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. Results are equivalent to Path A when tools, versions, reference, and parameters match. Quantify every sample identically.

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

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 download links with the gget skill (gget ref -w dna,gtf,cdna <species>), or from Ensembl/GENCODE directly. Keep genome and GTF from the same release.

1. QC raw reads — FastQC

bash

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:

bash

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):

bash

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.

bash

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:

bash

STAR --runThreadN 12 --genomeDir star_index \
  --readFilesIn trimmed/s1_R1.fq.gz trimmed/s1_R2.fq.gz --readFilesCommand zcat \
  --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 quasi-mapping

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

bash

# 1. Decoys = all genome sequence names; gentrome = transcriptome THEN genome (order matters)
grep '^>' 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):

bash

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

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)

bash

featureCounts -T 8 -p --countReadPairs \
  -a annotation.gtf -g gene_id \
  -s 2 \                       # strandedness: 0 unstranded, 1 forward, 2 reverse
  -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 silently discards ~half the reads. 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/IS = unstranded).
Or run RSeQC infer_experiment.py -r genes.bed -i s1.bam on a STAR BAM.

Map the result to each tool:

Library	Salmon `-l`	featureCounts `-s`	STAR column (in `ReadsPerGene.out.tab`)
Unstranded	`IU` (auto `A`)	`0`	col 2
Forward (e.g. Ligation)	`ISF` (auto `A`)	`1`	col 3
Reverse (e.g. dUTP/TruSeq stranded)	`ISR` (auto `A`)	`2`	col 4

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

4. Aggregate QC — MultiQC

bash

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.