# Path A — nf-core/rnaseq

`nf-core/rnaseq` is the field-standard, community-audited pipeline for the reads → counts stage. It chains FastQC → trimming (Trim Galore or fastp) → optional contaminant/rRNA removal → alignment + quantification (STAR+Salmon, STAR+RSEM, or HISAT2) → tximport gene/transcript count merging → extensive QC → MultiQC, with reviewed defaults and per-process containers.

This file covers *how to run it and what comes out*. For the Nextflow engine itself — profiles, executors, containers (Docker/Singularity/Conda/Wave), HPC/cloud, `-resume` caching, offline/`nf-core pipelines download` — use the **`nextflow`** skill. Don't duplicate that here.

Current stable revision at time of writing: **3.26.0** (always pin with `-r`).

## Samplesheet

`nf-core/rnaseq` takes a CSV, not loose files. Columns (header required):

```csv
sample,fastq_1,fastq_2,strandedness
CONTROL_REP1,/data/ctrl1_R1.fastq.gz,/data/ctrl1_R2.fastq.gz,auto
CONTROL_REP2,/data/ctrl2_R1.fastq.gz,/data/ctrl2_R2.fastq.gz,auto
TREATED_REP1,/data/trt1_R1.fastq.gz,/data/trt1_R2.fastq.gz,auto
TREATED_REP2,/data/trt1_R1.fastq.gz,,auto
```

- **sample** — sample ID. Rows that share a `sample` value are treated as the same sample sequenced over multiple lanes and are merged.
- **fastq_1 / fastq_2** — paths or URLs to gzipped FASTQ. Leave `fastq_2` empty for single-end.
- **strandedness** — `auto` (recommended; the pipeline infers it with Salmon and warns on mismatch), or `forward` / `reverse` / `unstranded` if you know the kit. TruSeq stranded mRNA is typically `reverse`.

Validate before launching: `python scripts/validate_samplesheet.py --samplesheet samplesheet.csv`.

## Choosing the aligner / quantifier

Set with `--aligner` (genome alignment) or `--pseudo_aligner` (lightweight). Defaults are well chosen.

| Option | What it does | When |
|--------|--------------|------|
| `--aligner star_salmon` (default) | STAR genome alignment, Salmon quantifies against the transcriptome from the BAM | The standard, defensible default — gives counts **and** a genome BAM for QC/IGV |
| `--aligner star_rsem` | STAR + RSEM | You specifically need RSEM estimates |
| `--aligner hisat2` | HISAT2 alignment (no built-in transcript quant) | Lower memory than STAR |
| `--pseudo_aligner salmon` (+ `--skip_alignment`) | Salmon quasi-mapping only, no BAM | Fastest/lightest; you don't need a genome BAM |

`star_salmon` is recommended unless you have a specific reason otherwise. You can also add `--pseudo_aligner salmon` alongside an aligner to get both.

## Running

```bash
# Smoke-test first (tiny bundled data; proves the environment works)
nextflow run nf-core/rnaseq -r 3.26.0 -profile test,docker --outdir test_results

# Real run with an iGenomes reference key
nextflow run nf-core/rnaseq -r 3.26.0 \
  -profile docker \
  --input samplesheet.csv \
  --genome GRCh38 \
  --aligner star_salmon \
  --outdir results \
  -resume

# Or supply your own reference explicitly (more reproducible than iGenomes keys)
nextflow run nf-core/rnaseq -r 3.26.0 \
  -profile singularity \
  --input samplesheet.csv \
  --fasta /ref/genome.fa --gtf /ref/annotation.gtf \
  --aligner star_salmon --outdir results -resume
```

Generate a validated, documented command + params file interactively with `nf-core pipelines launch rnaseq` (see the `nextflow` skill).

### Useful parameters

- `--save_reference` — keep the built STAR/Salmon indices so re-runs and other projects don't rebuild them.
- `--trimmer trimgalore|fastp` — trimming tool (default `trimgalore`).
- `--remove_ribo_rna` — sortmerna rRNA depletion (use if libraries weren't poly-A/ribo-depleted, or to quantify rRNA contamination).
- `--extra_salmon_quant_args='--gcBias'` — pass tool flags through.
- `--skip_*` (e.g. `--skip_markduplicates`, `--skip_stringtie`) — drop stages you don't need.
- `--gencode` — set when using GENCODE (not Ensembl) annotation, so gene IDs/biotypes parse correctly.

Pick the reference release deliberately and record it. iGenomes keys (`--genome GRCh38`) are convenient but version-pinning your own `--fasta`/`--gtf` is more reproducible.

## Outputs

Key paths under `--outdir` (for `star_salmon`):

```
results/
├── multiqc/             # MultiQC report — read this first
├── star_salmon/
│   ├── salmon.merged.gene_counts.tsv                 # raw estimated gene counts (tximport, countsFromAbundance=no)
│   ├── salmon.merged.gene_counts_length_scaled.tsv   # length-scaled counts -> use for DESeq2
│   ├── salmon.merged.gene_tpm.tsv                     # TPM (for visualization, NOT for DESeq2)
│   ├── salmon.merged.gene_counts.rds                  # SummarizedExperiment (R)
│   ├── <SAMPLE>/                                      # per-sample Salmon quant dirs
│   └── deseq2_qc/                                     # PCA + sample-distance plots the pipeline already made
└── pipeline_info/        # execution report, software versions, params
```

**The pipeline already runs tximport**, so on Path A you do **not** need this skill's `build_counts_matrix.py`. Use the merged TSV directly.

## Handoff to PyDESeq2

`salmon.merged.gene_counts_length_scaled.tsv` is genes × samples, with a leading `gene_id` (and usually `gene_name`) column, and non-integer values. Round to integers for PyDESeq2:

```python
import pandas as pd

df = pd.read_csv("results/star_salmon/salmon.merged.gene_counts_length_scaled.tsv", sep="\t")
df = df.drop(columns=[c for c in ["gene_name"] if c in df.columns]).set_index("gene_id")
counts = df.round().astype(int)          # genes x samples, integer
counts.to_csv("counts.csv")              # hand to the pydeseq2 skill (it transposes to samples x genes)
```

Build `metadata.csv` (index = sample IDs matching the count columns; columns = `condition`, `batch`, …) and proceed with the **`pydeseq2`** skill. The pipeline's own `deseq2_qc/` PCA is a good first sanity check before you run your own contrasts. Length-scaled counts are appropriate to round and use directly — see `counts-and-handoff.md` for the reasoning and the alternative offset-based route.