alterlab-nf-core-sarek

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

The workflow-runner entry point for raw-reads-to-variants: drive the Nextflow nf-core/sarek 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, 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):

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.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.

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.10.0), a -profile chosen, and is Nextflow >=25.10.4 (the version 3.10.0 requires)?
• For somatic asks, does the samplesheet carry a status 1 tumor and a status 0 normal under the same patient?
• For WES/panel, were both --wes and --intervals <capture.bed> supplied?
• After the run, did you route VCF interpretation to the correct sibling skill rather than re-deriving variant meaning here?

References

• references/usage_3.10.0.md — pinned run command, profiles, --step/--aligner options, the --wes + --intervals pairing, BQSR preprocessing, sourced from the 3.10.0 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.