alterlab-nf-core-sarek

Files (7)

• SKILL.md
• evals.json · evals/
• caller_accuracy.md · references/
• manual_gatk4.md · references/
• samplesheet_schema.md · references/
• usage_3.10.0.md · references/
• make_samplesheet.py · scripts/

SKILL.md (first lines)

# nf-core/sarek — FASTQ-to-VCF Variant Calling

The workflow-runner entry point for raw-reads-to-variants: drive the
**Nextflow [nf-core/sarek](https://nf-co.re/sarek/3.10.0/) pipeline (pinned `-r 3.10.0`)**
to take germline or somatic short-read FASTQ through alignment, GATK4 duplicate
marking and base-quality recalibration, and SNV/indel calling, then hand the
resulting VCFs to the suite's database and parsing skills for interpretation.

This skill is the **command-line / workflow** counterpart to the suite's
Python-library bioinformatics skills. Use it for the *raw-data-to-VCF* leg;
use the library skills (`alterlab-pysam`, `alterlab-tiledbvcf`) once you hold a VCF.

## When to Use This Skill

Trigger this skill when the user wants to:

- Go from **FASTQ to VCF** — call variants on whole-genome (WGS) or whole-exome
  (WES) short reads.
- Run **germline** SNV/indel calling (one or many normal samples).
- Run **somatic / tumor-normal** calling (matched tumor + normal, or tumor-only).
- Use **nf-core/sarek** specifically, or want a reproducible "GATK
  best-practices alignment-to-VCF" pipeline without hand-writing every step.
- Resume a run from an intermediate **`--step`** (already have BAM/CRAM, only need
  recalibration or variant calling).

### Does NOT Trigger — route adjacent requests here

| The request is really about… | Route to |
|---|---|
| Parsing / filtering / reading an **existing** VCF/BAM in Python (pysam/htslib) | `alterlab-pysam` |
| **Storing / querying** large multi-sample variant stores (TileDB-VCF arrays) | `alterlab-tiledbvcf` |
| Clinical significance of a called variant (pathogenic/benign) | `alterlab-clinvar` |
| Population allele frequencies for a called variant | `alterlab-gnomad` |
| Somatic mutation catalogue / cancer census lookup | `alterlab-cosmic` |
| **RNA-seq** transcript/gene quantification (salmon/kallisto), not DNA variants | `alterlab-rnaseq-quant` |
| 16S/ITS **amplicon / microbiome** FASTQ → feature table | `alterlab-qiime2-amplicon` |
| Sequence **homology / similarity search** (BLAST+, DIAMOND) | `alterlab-blast` |
| Spatial transcriptomics neighborhood/SVG analysis | `alterlab-squidpy-spatial` |
| Differential **expression** stats from counts | `alterlab-pydeseq2` |

If the user has no workflow engine and cannot install Nextflow + containers,
do **not** refuse — fall back to the **manual GATK4 recipe** (below /
`references/manual_gatk4.md`).

## The #1 Correctness Trap: sarek's default caller is Strelka

Per the [3.10.0 usage docs](https://nf-co.re/sarek/3.10.0/docs/usage/), **when
`--tools` is not set, sarek runs preprocessing and then Strelka only.** It does
**not** default to GATK HaplotypeCaller or DeepVariant. Always set `--tools`
explicitly to match the user's intent:

| Intent | Pass |
|---|---|
| GATK4 best-practice germline | `--tools haplotypecaller` |
| Highest germline F1 (CNN) | `--tools deepvariant` |
| Somatic, matched tumor/normal | `--tools mutect2` (often `mutect2,strelka`) |
| Joint germline genotyping across a cohort | `--tools haplotypecaller --joint_germline` |
| GPU-accelerated germline (needs an NVIDIA GPU profile) | `--tools parabricks_haplotypecaller` |

`--tools` accepts (per the 3.10.0 schema): `deepvariant`, `freebayes`,
`haplotypecaller`, `parabricks_haplotypecaller`, `mutect2`, `lofreq`, `mpileup`,
`muse`, `strelka`, `sentieon_*`, structural-variant callers (`manta`, `tiddit`,
`indexcov`), CNV/purity tools (`ascat`, `cnvkit`, `controlfreec`), QC
(`ngscheckmate`, `msisensorpro`), `varlociraptor`, and the annotation tools
(`snpeff`, `vep`, `snpsift`, `bcfann`). Caller choice materially changes precision/recall — see
`references/caller_accuracy.md` for the nf-core benchmark (Hanssen et al., 2024).

## Pipeline (how to run it)

### 1. Build the samplesheet

sarek's input is a CSV. Required columns for `--step mapping`:
`patient`, `sample`, `lane`, `fastq_1`, `fastq_2`. Optional: `sex` (XX/XY,
default NA) and `status` (**`0` = normal, `1` = tumor**, default 0) — `status`
is what tells sarek a pair is somatic.

Use the helper to generate a valid sheet from a FASTQ directory (it pairs R1/R2,
fills `lane`, and validates the schema before you burn compute):

```bash
uv run python skills/bioinformatics/alterlab-nf-core-sarek/scripts/make_samplesheet.py \
    --fastq-dir ./fastq --patient PATIENT_01 --sample TUMOR_01 \
    --status 1 --sex XY --out samplesheet.csv
```

Append more rows (e.g. the matched normal with `--status 0 --append`) before
running. See `references/samplesheet_schema.md` for every column, BAM/CRAM
re-entry rows, and a tumor-normal example.

### 2. Run the pipeline (pinned)

```bash
nextflow run nf-core/sarek -r 3.10.0 \
    -profile docker \
    --input samplesheet.csv \
    --outdir ./results \
    --genome GATK.GRCh38 \
    --tools haplotypecaller \
    --aligner bwa-mem2
```

- **Always keep `-r 3.10.0`** — unpinned runs drift to a different pipeline version.
- `-profile` is **mandatory**: `docker`, `singularity`, `apptainer`, or `conda`
  for the local environment (clusters add `test`, institutional configs, etc.).
- `--genome GATK.GRCh38` selects the iGenomes/GATK GRCh38 reference and its
  bundled BQSR known-sites (dbSNP, Mills/1000G indels) automatically.
- `--aligner` options: `bwa-mem` (default), `bwa-mem2`, `dragmap`, `sentieon-bwamem`,
  `parabricks` (GPU).
- For **WES/panel**, pass **`--wes`** *and* `--intervals targets.bed`: `--intervals`
  restricts where calling happens, while `--wes` switches the tools to
  targeted-sequencing settings. Exome data run without `--wes` completes happily with
  WGS-tuned thresholds.
- Resume mid-pipeline with `--step` (`mapping` default, then `markduplicates`,
  `prepare_recalibration`, `recalibrate`, `variant_calling`, `annotate`) and
  Nextflow's `-resume`.

Preprocessing follows GATK best practice: align → **MarkDuplicates** →
**BaseRecalibrator/ApplyBQSR** (BQSR) → variant calling. Details and every flag:
`references/usage_3.10.0.md`.

From AlterLab-IEU/AlterLab-Academic-Skills, MIT.

More from AlterLab-IEU/AlterLab-Academic-Skills