{
  "skill": "alterlab-nf-core-sarek",
  "evals": [
    {
      "id": "germline-fastq-to-vcf",
      "prompt": "I have paired-end whole-genome FASTQ files for a single patient and I want to call germline SNVs and indels following GATK best practices. What's the cleanest way to go from FASTQ to a VCF?",
      "expected_output": "Invokes alterlab-nf-core-sarek: builds a samplesheet.csv (patient, sample, lane, fastq_1, fastq_2) via scripts/make_samplesheet.py, then runs `nextflow run nf-core/sarek -r 3.10.0 -profile docker --input samplesheet.csv --outdir ./results --genome GATK.GRCh38 --tools haplotypecaller`. Pins the revision to 3.10.0, sets --tools haplotypecaller explicitly for GATK best practice (rather than letting it default to Strelka), describes the MarkDuplicates + BQSR preprocessing against the GRCh38 known-sites bundle, and routes the resulting VCF to alterlab-pysam / annotation skills.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "haplotypecaller" },
        { "type": "behavior", "value": "Pins -r 3.10.0 and sets --tools explicitly to haplotypecaller for GATK best practice instead of relying on the Strelka default." }
      ]
    },
    {
      "id": "somatic-tumor-normal",
      "prompt": "We sequenced a tumor and its matched normal blood sample for one patient. I need to call somatic point mutations. Can you set up the variant-calling pipeline?",
      "expected_output": "Invokes alterlab-nf-core-sarek in somatic mode: constructs a samplesheet with the tumor row (status 1) and the matched normal row (status 0) under the SAME patient ID, runs nf-core/sarek -r 3.10.0 with --tools mutect2 (optionally mutect2,strelka) and --genome GATK.GRCh38. Emphasizes that tumor and normal must share one patient identifier for somatic pairing, and hands the somatic VCF to alterlab-cosmic / alterlab-gnomad for annotation.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "mutect2" },
        { "type": "behavior", "value": "Sets status 1 (tumor) and status 0 (normal) under the same patient and selects a somatic caller (mutect2); does not treat the two samples as unrelated germline samples." }
      ]
    },
    {
      "id": "default-caller-trap-deepvariant",
      "prompt": "I ran nf-core/sarek without setting --tools and I'm confused — I expected GATK HaplotypeCaller output but I got something else. Also, what's the most accurate germline caller I could use instead?",
      "expected_output": "Invokes alterlab-nf-core-sarek: explains the #1 trap — with --tools unset, sarek 3.10.0 runs preprocessing and Strelka only, NOT HaplotypeCaller, so the user got Strelka output. Recommends re-running with an explicit --tools (haplotypecaller for GATK best practice, or deepvariant for highest germline F1 per the nf-core/Hanssen et al. 2024 benchmark). Cites that the default is Strelka and that --tools must be set explicitly.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "Strelka" },
        { "type": "behavior", "value": "Correctly identifies that an unset --tools makes sarek default to Strelka (not HaplotypeCaller) and recommends deepvariant or haplotypecaller explicitly." }
      ]
    },
    {
      "id": "resume-from-recalibrated-bam",
      "prompt": "I already have recalibrated BAM files from a previous sarek run on these samples — I don't want to re-align everything. Can I just re-run the variant calling step with DeepVariant?",
      "expected_output": "Invokes alterlab-nf-core-sarek: directs the user to re-enter the pipeline at --step variant_calling with a samplesheet carrying the bam/bai (CRAM/CRAI) re-entry columns rather than fastq_1/fastq_2, sets --tools deepvariant, keeps -r 3.10.0, and uses Nextflow -resume to reuse cached work. Avoids re-running mapping/markduplicates/BQSR.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "variant_calling" },
        { "type": "behavior", "value": "Uses --step variant_calling with BAM/CRAM re-entry rows instead of starting from FASTQ, so alignment and BQSR are not repeated." }
      ]
    },
    {
      "id": "near-miss-pysam-filter-vcf",
      "prompt": "I already have a multi-sample VCF from a variant-calling run. I just need to filter it in Python to keep variants with QUAL > 30 and DP > 10, and pull out all the variants in a specific gene region. How do I do that with pysam?",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-pysam. The user already holds a VCF and wants to read/filter/region-query it programmatically with pysam/htslib — that is in-memory VCF parsing, not running the FASTQ-to-VCF pipeline. alterlab-nf-core-sarek is for producing variants from raw reads, not post-hoc filtering of an existing VCF.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-pysam" }
      ]
    },
    {
      "id": "near-miss-tiledbvcf-store-cohort",
      "prompt": "I have thousands of per-sample VCFs and flat-file storage is becoming unmanageable for region/sample queries. I want to ingest them into a compressed, incrementally-updatable variant store I can query fast. What should I use?",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-tiledbvcf. The user wants scalable storage and querying of already-called variants in TileDB-VCF arrays, not to run alignment-to-VCF variant calling. alterlab-nf-core-sarek produces the VCFs; storing and querying them at population scale is alterlab-tiledbvcf's job.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-tiledbvcf" }
      ]
    },
    {
      "id": "near-miss-rnaseq-quant",
      "prompt": "I have bulk RNA-seq FASTQ files and I want to quantify transcript and gene-level expression with salmon before differential expression. How do I build the index and quantify?",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-rnaseq-quant. The user wants RNA-seq transcript/gene quantification (salmon/kallisto), not DNA short-variant calling. alterlab-nf-core-sarek is a FASTQ-to-VCF germline/somatic variant pipeline and has nothing to do with expression quantification.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-rnaseq-quant" }
      ]
    }
  ]
}

# Choosing `--tools`: caller selection and accuracy

## The default trap

In nf-core/sarek 3.10.0, **when `--tools` is not specified, the pipeline runs
preprocessing and then Strelka only** (https://nf-co.re/sarek/3.10.0/docs/usage/).
It does not silently run GATK HaplotypeCaller or DeepVariant. Always set
`--tools` to match intent.

## Tool/assay matrix (from the 3.10.0 docs)

Which callers apply to which assay and analysis type:

| Tool (`--tools` value) | WGS | WES | Germline | Tumor-only | Somatic |
|---|:---:|:---:|:---:|:---:|:---:|
| `deepvariant` | x | x | x | | |
| `freebayes` | x | x | x | x | x |
| `haplotypecaller` | x | x | x | | |
| `mutect2` | x | x | | x | x |
| `lofreq` | x | x | | x | |
| `mpileup` | x | x | x | x | |
| `strelka` | x | x | | | x |

(`--tools` also accepts annotation tools such as `snpeff`/`vep`; pass several
comma-separated, e.g. `--tools haplotypecaller,vep`.)

## Recommended choices

| Scenario | `--tools` |
|---|---|
| Germline, GATK4 best practice | `haplotypecaller` |
| Germline, maximize F1 (CNN caller) | `deepvariant` |
| Germline cohort, joint genotyping | `haplotypecaller` + `--joint_germline` |
| Somatic, matched tumor/normal | `mutect2` (or `mutect2,strelka`) |
| Somatic, tumor-only | `mutect2` or `freebayes` |

## Benchmark grounding

nf-core/sarek's own benchmarking study compared callers and infrastructures:

> Hanssen, F., Garcia, M. U., Folkersen, L., Pedersen, A. S., Lescai, F.,
> Jodoin, S., Miller, E., Seybold, M., Wacker, O., Smith, N., Gabernet, G.,
> & Nahnsen, S. (2024). Scalable and efficient DNA sequencing analysis on
> different compute infrastructures aiding variant discovery.
> *NAR Genomics and Bioinformatics, 6*(2), lqae031.
> https://doi.org/10.1093/nargab/lqae031

Reported directional findings (no exact metric values reproduced here — consult
the paper for figures):

- **Germline:** Strelka2 (with Manta) and DeepVariant performed best across the
  evaluated precision/recall/F1 metrics.
- **Somatic:** Mutect2 (filtered) had the highest precision; FreeBayes the
  highest recall; the highest somatic F1 was Mutect2, followed by Strelka2.
- **Aligners:** BWA-MEM and BWA-MEM2 yielded higher recall than DragMap.

Treat these as guidance for `--tools`/`--aligner`, not as a substitute for
benchmarking on your own truth set (e.g. GIAB) when accuracy is critical.

# Manual GATK4 fallback (no Nextflow)

Use this when the user cannot run Nextflow + a container engine. It reproduces
the **same GATK best-practices preprocessing chain** that sarek automates, using
the GATK4 CLI directly. It is more error-prone and you must manage references
yourself — prefer the pipeline when available.

Sources: GATK Best Practices "Data pre-processing for variant discovery"
(https://gatk.broadinstitute.org/hc/en-us/articles/360035535912) and the
germline short-variant discovery best-practice (HaplotypeCaller in GVCF mode →
GenotypeGVCFs). The known-sites resources match those the sarek 3.10.0 GATK genome
key uses (dbSNP + Mills/1000G gold-standard indels).

## Tools needed (bioconda)

`bwa-mem2`, `samtools`, `gatk4`. Install in an isolated env (e.g. conda/mamba);
these run offline once the reference and known-sites are local.

## Reference inputs

- `ref.fasta` — GRCh38 reference (with `.fai` and `.dict`).
- `dbsnp.vcf.gz` — dbSNP known sites.
- `mills_1000G.indels.vcf.gz` — Mills and 1000G gold-standard indels.

All three are in the Broad GATK resource bundle for hg38
(`genomics-public-data/resources/broad/hg38/v0/`).

## Germline single-sample chain

```bash
# 0. Index the reference for bwa-mem2 (once)
bwa-mem2 index ref.fasta

# 1. Align (set a read group; required downstream)
bwa-mem2 mem -t 8 -R '@RG\tID:L001\tSM:NORMAL_01\tPL:ILLUMINA\tLB:lib1' \
    ref.fasta R1.fastq.gz R2.fastq.gz \
  | samtools sort -@ 8 -o aligned.bam -
samtools index aligned.bam

# 2. Mark duplicates
gatk MarkDuplicates -I aligned.bam -O markdup.bam -M markdup.metrics.txt
samtools index markdup.bam

# 3. BQSR — build the recalibration table from known sites
gatk BaseRecalibrator \
    -I markdup.bam -R ref.fasta \
    --known-sites dbsnp.vcf.gz \
    --known-sites mills_1000G.indels.vcf.gz \
    -O recal.table

# 4. BQSR — apply it
gatk ApplyBQSR -I markdup.bam -R ref.fasta --bqsr-recal-file recal.table -O recal.bam
samtools index recal.bam

# 5. Call per-sample variants in GVCF mode
gatk HaplotypeCaller -I recal.bam -R ref.fasta -ERC GVCF -O sample.g.vcf.gz

# 6. Genotype to a final VCF (single-sample shortcut)
gatk GenotypeGVCFs -R ref.fasta -V sample.g.vcf.gz -O sample.vcf.gz
```

## Joint genotyping (cohort)

Produce one `*.g.vcf.gz` per sample (steps 1–5), then combine and joint-genotype:

```bash
gatk CombineGVCFs -R ref.fasta -V a.g.vcf.gz -V b.g.vcf.gz -O cohort.g.vcf.gz
gatk GenotypeGVCFs -R ref.fasta -V cohort.g.vcf.gz -O cohort.vcf.gz
```

(For large cohorts the best-practice path uses `GenomicsDBImport` instead of
`CombineGVCFs`; consult the GATK docs.)

## WES

Restrict every region-aware step with `-L targets.bed` (the capture-kit BED) —
the manual analogue of sarek's `--intervals`.

## After the VCF

Hand off exactly as with the pipeline: `alterlab-pysam` (parse/filter),
`alterlab-tiledbvcf` (store/query), `alterlab-clinvar` / `alterlab-gnomad` /
`alterlab-cosmic` (annotate).

> Caveat: command flags can change between GATK4 minor versions. Verify against
> `gatk <Tool> --help` for the installed version before running on real data.

# nf-core/sarek 3.10.0 — Samplesheet (`--input`) Schema

Source: https://nf-co.re/sarek/3.10.0/docs/usage/. The input is a comma-separated
CSV with a header row. Columns depend on the `--step` you start from.

## Columns for `--step mapping` (FASTQ entry)

| Column | Required | Meaning |
|---|---|---|
| `patient` | yes | Subject identifier. Tumor and matched normal share the **same** `patient`. |
| `sample` | yes | Biological sample identifier (unique per sample). |
| `lane` | yes (mapping) | Sequencing lane/run identifier; lets sarek track read groups and merge lanes. |
| `fastq_1` | yes | Path to the R1 gzipped FASTQ. |
| `fastq_2` | yes (paired-end) | Path to the R2 gzipped FASTQ. |
| `sex` | optional | `XX` / `XY` (default `NA`). |
| `status` | optional | **`0` = normal, `1` = tumor** (default `0`). Drives somatic pairing. |

### Germline example (single normal sample)

```csv
patient,sex,status,sample,lane,fastq_1,fastq_2
PATIENT_01,XY,0,NORMAL_01,L001,/data/N_R1.fastq.gz,/data/N_R2.fastq.gz
```

### Somatic example (matched tumor + normal, same patient)

```csv
patient,sex,status,sample,lane,fastq_1,fastq_2
PATIENT_01,XY,0,NORMAL_01,L001,/data/N_R1.fastq.gz,/data/N_R2.fastq.gz
PATIENT_01,XY,1,TUMOR_01,L001,/data/T_R1.fastq.gz,/data/T_R2.fastq.gz
```

The tumor row (`status 1`) plus a normal row (`status 0`) under one `patient` is
what makes the run somatic. Run with `--tools mutect2` (or `mutect2,strelka`).
For **tumor-only** somatic calling, include only the `status 1` row.

### Multiple lanes

Repeat the sample with different `lane` values; sarek aligns each lane and merges:

```csv
patient,sex,status,sample,lane,fastq_1,fastq_2
PATIENT_01,XY,0,NORMAL_01,L001,/data/N_L001_R1.fastq.gz,/data/N_L001_R2.fastq.gz
PATIENT_01,XY,0,NORMAL_01,L002,/data/N_L002_R1.fastq.gz,/data/N_L002_R2.fastq.gz
```

## Re-entry rows (resuming with `--step`)

When you already have aligned/recalibrated data, the FASTQ columns are replaced
by alignment columns and you pick the matching `--step`:

- **BAM/CRAM re-entry** (`--step markduplicates` / `prepare_recalibration` /
  `recalibrate` / `variant_calling`): provide `bam`/`bai` or `cram`/`crai`
  (and, for the recalibration steps, the BQSR `table`) instead of
  `fastq_1`/`fastq_2`.
- **VCF re-entry** (`--step annotate`): provide a `vcf` column.

Consult the 3.10.0 usage docs for the exact column set required by each step; the
helper script (`scripts/make_samplesheet.py`) writes the FASTQ-entry sheet for
`--step mapping`, which is the common starting point.

## Common mistakes

- Putting tumor and normal under **different** `patient` IDs → sarek treats them
  as unrelated germline samples and never does somatic calling.
- Omitting `lane` for `--step mapping` (it is required there).
- Leaving `--tools` unset and assuming GATK ran — it ran **Strelka** by default
  (see `caller_accuracy.md`).

# nf-core/sarek 3.10.0 — Usage Reference

Source: https://nf-co.re/sarek/3.10.0/docs/usage/ (pinned release `3.10.0`,
"Aktse", released 2026-08-12). Everything below is for that pinned version.
Newer releases may rename or change defaults — keep `-r 3.10.0` unless the user
explicitly asks to upgrade.

`3.10.0` declares `nextflowVersion = '!>=25.10.4'`, so an older Nextflow refuses
to launch it. The release also completed the migration to topic channels and to
Nextflow **strict syntax**, which becomes the parser in Nextflow 26.x — relevant
if you maintain custom local modules or a `modules.config` alongside it.

## Minimal command

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

- `-r 3.10.0` — pins the pipeline revision. Required for reproducibility.
- `-profile` — **mandatory**, reflects the compute/software environment. Common
  values: `docker`, `singularity`, `apptainer`, `conda`. A `test` profile runs a
  tiny built-in dataset. Combine with institutional configs as needed.
- `--input` — path to the samplesheet CSV (see `samplesheet_schema.md`).
- `--outdir` — results directory (required).
- `--genome` — iGenomes/GATK reference key, e.g. `GATK.GRCh38` or `GATK.GRCh37`.
  Selecting the GATK key wires up the BQSR known-sites automatically (below).
- `--tools` — caller selection. **Defaults to Strelka if omitted** (see
  `caller_accuracy.md`). Set it explicitly.

## Aligner (`--aligner`)

- `bwa-mem` — default.
- `bwa-mem2` — faster, same algorithm family.
- `dragmap` — DRAGEN-style mapper.
- `sentieon-bwamem` — Sentieon's BWA implementation (needs a Sentieon licence).
- `parabricks` — GPU-accelerated; requires a GPU profile
  (e.g. `--aligner parabricks -profile docker,gpu`).

The nf-core benchmark (Hanssen et al., 2024) reports BWA-MEM and BWA-MEM2 give
higher recall than DragMap; see `caller_accuracy.md`.

## Steps (`--step`)

Start or resume the pipeline at a point matching the inputs you already have:

| `--step` | Starts from |
|---|---|
| `mapping` (default) | FASTQ → alignment |
| `markduplicates` | mapped BAM/CRAM → duplicate marking |
| `prepare_recalibration` | → BaseRecalibrator (build recal table) |
| `recalibrate` | → ApplyBQSR |
| `variant_calling` | recalibrated BAM/CRAM → calling |
| `annotate` | existing VCF → annotation only |

Pair `--step` with the matching samplesheet columns (BAM/CRAM/VCF re-entry rows;
see `samplesheet_schema.md`) and Nextflow `-resume` to reuse cached work.

## Preprocessing = GATK best practice

With the GATK genome key, sarek runs the GATK data-pre-processing chain:

1. **Alignment** (`--aligner`).
2. **MarkDuplicates** — flag PCR/optical duplicates.
3. **BQSR** — `BaseRecalibrator` builds a recalibration table from known sites,
   then `ApplyBQSR` writes recalibrated reads.

### BQSR known-sites / GRCh38 resource bundle

These reference resources ship with the GATK genome key and feed recalibration
and calling (per the 3.10.0 reference table):

| Resource | Role | Used by (per docs) |
|---|---|---|
| **dbSNP** | known SNP sites | BaseRecalibrator, GenotypeGVCFs, HaplotypeCaller, ControlFREEC |
| **Mills and 1000G gold-standard indels** (`known_indels`) | known indel sites | BaseRecalibrator(Spark), FilterVariantTranches |

Source bundle: the Broad GATK resource bundle for hg38
(`genomics-public-data/resources/broad/hg38/v0/`), referenced by the pipeline as
`GATKBundle`. With `--genome GATK.GRCh38` you do not supply these by hand.

## Whole-exome / panel (WES)

Pass **both**:

- `--wes` — a boolean that flips targeted-sequencing settings in the individual
  tools ("Enable when exome or panel data is provided"), and
- `--intervals targets.bed` — the capture-kit BED, so calling is restricted to the
  targeted regions.

They do different jobs: `--intervals` limits *where* calling happens; `--wes` tells
the callers and QC modules that coverage is targeted rather than uniform. Running
exome data without `--wes` produces results that look fine but carry WGS-tuned
thresholds. (`--intervals` also speeds WGS by parallelizing over interval lists.)

## Joint germline

Add `--joint_germline` to `--tools haplotypecaller` to run joint genotyping
across a GVCF cohort (HaplotypeCaller in GVCF mode → joint GenotypeGVCFs).

## Output layout

Results are written under `--outdir`, including:

- `preprocessing/` — recalibrated CRAM/BAM and BQSR tables.
- `variant_calling/<tool>/` — per-caller VCFs (e.g. `haplotypecaller/`,
  `mutect2/`, `strelka/`, `deepvariant/`).
- `reports/` — MultiQC and per-tool QC.

Hand the VCFs to `alterlab-pysam` (parse), `alterlab-tiledbvcf` (store/query),
and `alterlab-clinvar` / `alterlab-gnomad` / `alterlab-cosmic` (annotate).

#!/usr/bin/env python3
"""Build and validate an nf-core/sarek 3.8.1 `--input` samplesheet from a FASTQ dir.

Generates the FASTQ-entry samplesheet for `--step mapping`, the common starting
point. Pairs R1/R2 files in a directory, fills the required columns
(patient, sample, lane, fastq_1, fastq_2) and the optional sex/status columns,
validates the schema, and writes (or appends to) a CSV ready for:

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

Schema source: https://nf-co.re/sarek/3.8.1/docs/usage/
  Required (mapping): patient, sample, lane, fastq_1, fastq_2
  Optional: sex (XX/XY, default NA), status (0=normal / 1=tumor, default 0)

Stdlib only — no third-party deps; runs in a bare `uv run python` env. It does
NOT call any network service and does NOT run the pipeline; it only prepares the
samplesheet. Review the CSV before launching a (compute-heavy) run.

Examples
--------
    # One tumor sample, auto-pairing R1/R2 in ./fastq
    uv run python make_samplesheet.py --fastq-dir ./fastq \\
        --patient PATIENT_01 --sample TUMOR_01 --status 1 --sex XY \\
        --out samplesheet.csv

    # Append the matched normal under the SAME patient
    uv run python make_samplesheet.py --fastq-dir ./fastq_normal \\
        --patient PATIENT_01 --sample NORMAL_01 --status 0 --sex XY \\
        --out samplesheet.csv --append
"""

from __future__ import annotations

import argparse
import csv
import re
import sys
from pathlib import Path

HEADER = ["patient", "sex", "status", "sample", "lane", "fastq_1", "fastq_2"]

# R1/R2 mate tokens seen in Illumina-style names: _R1 / _R1_001 / _1 before .fastq.gz
R1_PAT = re.compile(r"(.+?)([._-])R?1((?:_\d+)?)(\.f(?:ast)?q\.gz)$", re.IGNORECASE)


def _mate2_name(r1: Path) -> Path:
    """Return the expected R2 path for an R1 filename, or raise if it isn't an R1."""
    m = R1_PAT.search(r1.name)
    if not m:
        raise ValueError(f"Not an R1-style FASTQ name: {r1.name}")
    stem, sep, tail, ext = m.groups()
    # Mirror whichever token form was used (R1 -> R2, or 1 -> 2).
    token = "R2" if re.search(r"R1", r1.name[m.start(2):], re.IGNORECASE) else "2"
    # Reconstruct using the exact separator + optional _001 tail + extension.
    return r1.with_name(f"{stem}{sep}{token}{tail}{ext}")


def pair_fastqs(fastq_dir: Path) -> list[tuple[Path, Path]]:
    """Find R1/R2 pairs in a directory. Returns sorted (r1, r2) absolute-path pairs."""
    if not fastq_dir.is_dir():
        raise NotADirectoryError(f"--fastq-dir not found: {fastq_dir}")
    candidates = sorted(p for p in fastq_dir.iterdir()
                        if p.is_file() and re.search(r"\.f(ast)?q\.gz$", p.name, re.IGNORECASE))
    pairs: list[tuple[Path, Path]] = []
    for r1 in candidates:
        if not R1_PAT.search(r1.name):
            continue  # skip R2 (and unmatched) files; they're picked up via their R1
        r2 = _mate2_name(r1)
        if not r2.exists():
            raise FileNotFoundError(f"R1 {r1.name} has no matching R2 (expected {r2.name})")
        pairs.append((r1.resolve(), r2.resolve()))
    if not pairs:
        raise FileNotFoundError(
            f"No R1/R2 FASTQ pairs found in {fastq_dir} "
            "(expected names like *_R1_001.fastq.gz / *_R1.fq.gz)")
    return pairs


def build_rows(pairs, patient, sample, status, sex, lane_prefix):
    """One CSV row per R1/R2 pair; lanes auto-numbered L001, L002, ... if >1 pair."""
    rows = []
    multi = len(pairs) > 1
    for i, (r1, r2) in enumerate(pairs, start=1):
        lane = f"{lane_prefix}{i:03d}" if multi else f"{lane_prefix}001"
        rows.append({
            "patient": patient, "sex": sex, "status": str(status),
            "sample": sample, "lane": lane,
            "fastq_1": str(r1), "fastq_2": str(r2),
        })
    return rows


def validate_rows(rows):
    """Raise ValueError on any schema violation. Mirrors sarek 3.8.1 requirements."""
    errors = []
    seen_keys = set()
    for n, row in enumerate(rows, start=1):
        for col in ("patient", "sample", "lane", "fastq_1", "fastq_2"):
            if not row.get(col):
                errors.append(f"row {n}: missing required column '{col}'")
        if row.get("status") not in {"0", "1"}:
            errors.append(f"row {n}: status must be 0 (normal) or 1 (tumor), got {row.get('status')!r}")
        if row.get("sex") not in {"XX", "XY", "NA"}:
            errors.append(f"row {n}: sex should be XX/XY/NA, got {row.get('sex')!r}")
        key = (row["patient"], row["sample"], row["lane"])
        if key in seen_keys:
            errors.append(f"row {n}: duplicate patient/sample/lane {key}")
        seen_keys.add(key)
    if errors:
        raise ValueError("Samplesheet validation failed:\n  - " + "\n  - ".join(errors))


def write_csv(rows, out: Path, append: bool):
    existing = []
    if append and out.exists():
        with out.open(newline="") as fh:
            existing = list(csv.DictReader(fh))
    all_rows = existing + rows
    validate_rows(all_rows)  # re-validate the merged sheet (catches cross-row dups)
    with out.open("w", newline="") as fh:
        w = csv.DictWriter(fh, fieldnames=HEADER)
        w.writeheader()
        w.writerows(all_rows)
    return len(all_rows)


def main(argv=None):
    ap = argparse.ArgumentParser(description=__doc__.split("\n")[0],
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument("--fastq-dir", required=True, type=Path, help="Directory of R1/R2 .fastq.gz files")
    ap.add_argument("--patient", required=True, help="Subject ID (tumor+normal share this)")
    ap.add_argument("--sample", required=True, help="Sample ID (unique per sample)")
    ap.add_argument("--status", type=int, choices=(0, 1), default=0,
                    help="0=normal, 1=tumor (default 0)")
    ap.add_argument("--sex", choices=("XX", "XY", "NA"), default="NA", help="default NA")
    ap.add_argument("--lane-prefix", default="L", help="Lane label prefix (default 'L')")
    ap.add_argument("--out", type=Path, default=Path("samplesheet.csv"))
    ap.add_argument("--append", action="store_true",
                    help="Append to an existing sheet (e.g. add the matched normal)")
    args = ap.parse_args(argv)

    try:
        pairs = pair_fastqs(args.fastq_dir)
        rows = build_rows(pairs, args.patient, args.sample, args.status, args.sex, args.lane_prefix)
        total = write_csv(rows, args.out, args.append)
    except (ValueError, FileNotFoundError, NotADirectoryError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        return 1

    print(f"Wrote {len(rows)} row(s) for sample '{args.sample}' "
          f"({'tumor' if args.status == 1 else 'normal'}); sheet now has {total} row(s): {args.out}")
    print("Next: nextflow run nf-core/sarek -r 3.8.1 -profile docker "
          f"--input {args.out} --outdir ./results --genome GATK.GRCh38 --tools haplotypecaller")
    print("Reminder: set --tools explicitly — sarek defaults to Strelka if you omit it.")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())