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.8.1/) pipeline (pinned `-r 3.8.1`) 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

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.8.1 usage docs, 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:

--tools accepts (per the docs' tool matrix): deepvariant, freebayes, haplotypecaller, mutect2, lofreq, mpileup, strelka (and annotation tools). 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):

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)

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

• Always keep `-r 3.8.1` — 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.
• For WES, pass --intervals targets.bed with the capture-kit BED (there is

no --wes flag in 3.8.1; restrict to target regions via --intervals).

• 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.8.1.md.

3. Interpret the output VCFs

Per-caller VCFs land under results/variant_calling/<tool>/. Then:

• Parse / filter with alterlab-pysam.
• Store / query at scale (multi-sample) with alterlab-tiledbvcf.
• Annotate clinical significance → alterlab-clinvar; population frequency →

alterlab-gnomad; somatic catalogue → alterlab-cosmic.

Fallback: manual GATK4 (no Nextflow)

If the user cannot run Nextflow + containers, run the equivalent GATK4 best-practices chain by hand: bwa-mem2 mem → gatk MarkDuplicates → gatk BaseRecalibrator + gatk ApplyBQSR (with dbSNP + Mills/1000G known sites) → gatk HaplotypeCaller -ERC GVCF → gatk GenotypeGVCFs. Full command sequence and the resource-bundle paths are in references/manual_gatk4.md.

Self-Check Before Reporting

• Is --tools set explicitly? Never let a run fall through to the Strelka

default unless the user truly wants Strelka.

• Is the version pinned (-r 3.8.1) and a -profile chosen?
• For somatic asks, does the samplesheet carry a status 1 tumor and a

status 0 normal under the same `patient`?

• For WES, was --intervals supplied with the capture BED?
• After the run, did you route VCF interpretation to the correct sibling skill

rather than re-deriving variant meaning here?

References

• references/usage_3.8.1.md — pinned run command, profiles, --step/--aligner

options, BQSR preprocessing, sourced from the 3.8.1 usage docs.

• references/samplesheet_schema.md — full CSV column spec, BAM/CRAM re-entry,

tumor-normal worked example.

• references/caller_accuracy.md — choosing --tools, summarizing the nf-core

benchmark (Hanssen et al., 2024, NAR Genomics & Bioinformatics).

• references/manual_gatk4.md — the non-Nextflow GATK4 best-practices fallback.

Part of the AlterLab Academic Skills suite.