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

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

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

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