"""Snakemake workflow: Reproduce MacDonald et al. (2022) GRCh38 LD blocks.

Reproduces the EUR / deCODE analysis and the published pyrho superpopulation
analyses from:
    MacDonald et al. (2022), https://github.com/jmacdon/LDblocks_GRCh38

Pipeline steps:
    0  Data download (VCFs, pyrho maps, deCODE interpolated maps, centromeres)
    1  VCF filtering with bcftools (MAF ≥ 0.01, SNPs only)
    2  Individual list preparation (1000G population dirs ∩ VCF samples)
    2a-2c deCODE map conversion/interpolation -- diagnostic only (regression-
       checks our own ldetect interpolate-maps against MacDonald's published
       maps via validate_maps/compare_maps); the EUR block set itself uses
       MacDonald's published deCODE_interpolated_maps directly, same as pyrho.
       compare_maps/compare_maps_point_mode together are an interpolation
       equivalence test: ldetect interpolate-maps --mode interval is a
       direct port of MacDonald's R interpolation scripts and should match
       their published maps almost exactly, while --mode point (the prior
       default) reproduces the interval-shift bug those scripts fixed -- see
       notes/logs/macdonald2022-interpolation-port.md
    3  LD block detection per chromosome (ldetect run)
    4  Post-processing: centromere removal + small block merging
    5  Genome-wide BED assembly (cat all chromosomes)

Usage:
    cd examples/MacDonald2022/
    uv run snakemake -n                         # dry-run, all chromosomes
    uv run snakemake --cores 8                  # run with 8 parallel jobs
    uv run snakemake --config chromosomes=[22]  # smoke-test on chr22 only
"""

import re
from pathlib import Path

# ---------------------------------------------------------------------------
# Config
# ---------------------------------------------------------------------------


configfile: "config.yaml"


CHROMS = config["chromosomes"]
BLOCK_SETS = list(config["block_sets"].keys())
ACTIVE_BLOCK_SETS = config.get("active_block_sets", BLOCK_SETS)
BLOCK_SET_PATTERN = "|".join(re.escape(block_set) for block_set in BLOCK_SETS)
RAW_DIR = Path(config["raw_vcf_dir"])
FILT_DIR = Path(config["filtered_vcf_dir"])
MAPS_DIR = Path(config["maps_dir"])
RES_DIR = Path(config["resources_dir"])
RESULTS = Path(config["results_dir"])
VCF_TMPL = config["vcf_filename_template"]


wildcard_constraints:
    reference_bed=r"(deCODE_EUR|pyrho_(AFR|EAS|EUR|SAS))_LD_blocks\.bed",


def pop_ne(population):
    return config["populations"][population]["ne"]


def pop_subpops(population):
    return config["populations"][population]["subpops"]


def block_set_population(block_set):
    return config["block_sets"][block_set]["population"]


def block_set_ne(block_set):
    return pop_ne(block_set_population(block_set))


def block_set_reference_bed(block_set):
    return config["block_sets"][block_set]["reference_bed"]


def diagnostic_resource(filename):
    primary = RES_DIR / filename
    nested = RES_DIR / RES_DIR.name / filename
    if not primary.exists() and nested.exists():
        return str(nested)
    return str(primary)


def block_set_remove_centromeres(block_set):
    return config["block_sets"][block_set].get(
        "remove_centromere_blocks",
        config.get("remove_centromere_blocks", False),
    )


def block_set_maf_filter(block_set):
    """'nref' (default, matches MacDonald's literal bcftools command) or
    'minor' (true minor-allele-frequency filtering, matching MacDonald's
    README prose rather than their script -- see
    notes/logs/macdonald2022-pyrho-handoff.md for why this is worth testing).
    """
    return config["block_sets"][block_set].get("maf_filter", "nref")


def block_set_filtered_vcf(wc):
    if block_set_maf_filter(wc.block_set) == "minor":
        return str(FILT_DIR / "minor" / f"chr{wc.chrom}.vcf.gz")
    return str(FILT_DIR / f"chr{wc.chrom}.vcf.gz")


def block_set_genetic_map(wc):
    block_config = config["block_sets"][wc.block_set]
    if block_config["map_source"] == "decode":
        return str(MAPS_DIR / "decode_interpolated_maps" / f"chr{wc.chrom}.tab.gz")
    if block_config["map_source"] == "pyrho":
        return str(
            MAPS_DIR
            / "pyrho_interpolated_maps"
            / block_config["map_pop"]
            / f"chr{wc.chrom}.tab.gz"
        )
    raise ValueError(f"Unknown map source: {block_config['map_source']}")


def expected_individual_count_param(population):
    expected = config["populations"][population].get("expected_count")
    if expected is None:
        return ""
    return f"--expected-count {expected}"


def remove_centromeres_param(wc):
    if block_set_remove_centromeres(wc.block_set):
        return "--remove-centromeres"
    return ""


# ---------------------------------------------------------------------------
# Target rule
# ---------------------------------------------------------------------------


rule all:
    input:
        expand(
            str(RESULTS / "{block_set}_LD_blocks.bed"),
            block_set=ACTIVE_BLOCK_SETS,
        ),
        expand(
            str(RESULTS / "compare" / "{block_set}_block_comparison.tsv"),
            block_set=ACTIVE_BLOCK_SETS,
        ),
        expand(
            str(
                RESULTS / "compare" / "boundaries" / "{block_set}_boundary_offsets.tsv"
            ),
            block_set=ACTIVE_BLOCK_SETS,
        ),
        expand(
            str(RESULTS / "compare" / "raw" / "{block_set}_block_comparison.tsv"),
            block_set=ACTIVE_BLOCK_SETS,
        ),
        expand(
            str(RESULTS / "compare" / "raw" / "{block_set}_boundary_offsets.tsv"),
            block_set=ACTIVE_BLOCK_SETS,
        ),
        str(RESULTS / "compare" / "map_summary.tsv"),
        str(RESULTS / "compare" / "map_ref_comparison.tsv"),
        str(RESULTS / "compare" / "map_ref_comparison_point_mode.tsv"),


# ---------------------------------------------------------------------------
# Step 0: Download
# ---------------------------------------------------------------------------


rule download_vcf:
    output:
        vcf=str(RAW_DIR / VCF_TMPL).replace("{chrom}", "{chrom}"),
        tbi=str(RAW_DIR / VCF_TMPL).replace("{chrom}", "{chrom}") + ".tbi",
    params:
        base_url=config["vcf_base_url"],
        fname=lambda wc: VCF_TMPL.format(chrom=wc.chrom),
    shell:
        """
        mkdir -p {RAW_DIR}
        wget -q -O {output.vcf} {params.base_url}/{params.fname}
        wget -q -O {output.tbi} {params.base_url}/{params.fname}.tbi
        """


rule download_decode_map:
    output:
        str(MAPS_DIR / "decode_raw" / "aau1043_datas3.gz"),
    params:
        url=config["decode_map_url"],
    shell:
        """
        mkdir -p {MAPS_DIR}/decode_raw
        wget -q -O {output} '{params.url}' || {{
            echo "ERROR: Could not download deCODE map."
            echo "Download manually from: {params.url}"
            echo "Place at: {output}"
            exit 1
        }}
        """


rule download_centromeres:
    output:
        str(RES_DIR / "hg38_centromeres.txt.gz"),
    params:
        url=config["centromere_url"],
    shell:
        "mkdir -p {RES_DIR} && wget -q -O {output} {params.url}"


rule download_pyrho_map:
    output:
        str(MAPS_DIR / "pyrho_interpolated_maps" / "{map_pop}" / "chr{chrom}.tab.gz"),
    params:
        base_url=config["pyrho_map_base_url"],
    shell:
        """
        mkdir -p {MAPS_DIR}/pyrho_interpolated_maps/{wildcards.map_pop}
        wget -q -O {output} \
            {params.base_url}/{wildcards.map_pop}/chr{wildcards.chrom}.tab.gz
        """


# MacDonald et al. also publish their own already-interpolated deCODE maps
# (same chr/pos/cM format as the pyrho maps above), produced by their own R
# interpolation script. Use that directly for the EUR/deCODE block set instead
# of recomputing our own interpolation from the raw (paywalled) deCODE
# supplementary file -- avoids an interpolation-method mismatch, mirroring how
# the pyrho block sets already avoid it.
rule download_decode_interpolated_map:
    output:
        str(MAPS_DIR / "decode_interpolated_maps" / "chr{chrom}.tab.gz"),
    params:
        base_url=config["decode_interpolated_map_base_url"],
    shell:
        """
        mkdir -p {MAPS_DIR}/decode_interpolated_maps
        wget -q -O {output} {params.base_url}/chr{wildcards.chrom}.tab.gz
        """


# ---------------------------------------------------------------------------
# Step 1: Filter VCFs
# ---------------------------------------------------------------------------


rule filter_vcf:
    input:
        vcf=str(RAW_DIR / VCF_TMPL).replace("{chrom}", "{chrom}"),
        tbi=str(RAW_DIR / VCF_TMPL).replace("{chrom}", "{chrom}") + ".tbi",
    output:
        vcf=str(FILT_DIR / "chr{chrom}.vcf.gz"),
        tbi=str(FILT_DIR / "chr{chrom}.vcf.gz.tbi"),
    shell:
        """
        mkdir -p {FILT_DIR}
        bcftools view --min-af 0.01 --types snps -Oz -o {output.vcf} {input.vcf}
        bcftools index --tbi {output.vcf}
        """


# Experimental variant: MacDonald's README says they filter by "MAF>0.01",
# but their actual subsetVcf.sh runs plain `--min-af 0.01` with no `:minor`
# suffix, which bcftools defaults to non-reference (nref, effectively ALT)
# allele frequency, not minor allele frequency. This rule tests the prose's
# literal meaning instead. Kept as a fully separate output path (not a
# parameterized version of filter_vcf above) so it never invalidates or
# forces recomputation of the already-validated default-filtered pipeline --
# opt-in only via a block_set with `maf_filter: minor` in config.yaml, never
# part of `active_block_sets`. `:minor` and `:nonmajor` are numerically
# identical for this biallelic-SNV 1000G release, so only one is needed.
rule filter_vcf_minor:
    input:
        vcf=str(RAW_DIR / VCF_TMPL).replace("{chrom}", "{chrom}"),
        tbi=str(RAW_DIR / VCF_TMPL).replace("{chrom}", "{chrom}") + ".tbi",
    output:
        vcf=str(FILT_DIR / "minor" / "chr{chrom}.vcf.gz"),
        tbi=str(FILT_DIR / "minor" / "chr{chrom}.vcf.gz.tbi"),
    shell:
        """
        mkdir -p {FILT_DIR}/minor
        bcftools view --min-af 0.01:minor --types snps -Oz -o {output.vcf} {input.vcf}
        bcftools index --tbi {output.vcf}
        """


# ---------------------------------------------------------------------------
# Step 2: Individual list
# ---------------------------------------------------------------------------


rule prep_individuals:
    input:
        vcf=str(FILT_DIR / "chr22.vcf.gz"),  # any chromosome suffices
    output:
        str(RES_DIR / "{population}_inds.txt"),
    params:
        subpops=lambda wc: " ".join(pop_subpops(wc.population)),
        base_url=config["sample_data_base_url"],
        expected_count=lambda wc: expected_individual_count_param(wc.population),
    shell:
        """
        uv run python scripts/prep_individuals.py \
            --sample-data-base-url {params.base_url} \
            --subpops {params.subpops} \
            --vcf {input.vcf} \
            --output {output} \
            {params.expected_count}
        """


# ---------------------------------------------------------------------------
# Step 2a: Convert deCODE map to 3-column per chromosome
# ---------------------------------------------------------------------------


rule convert_decode_map:
    input:
        rules.download_decode_map.output[0],
    output:
        str(MAPS_DIR / "decode_raw" / "chr{chrom}.tab.gz"),
    shell:
        """
        uv run python scripts/convert_decode_map.py \
            --input {input} \
            --output-dir {MAPS_DIR}/decode_raw \
            --chromosomes {wildcards.chrom}
        """


# ---------------------------------------------------------------------------
# Step 2b: Extract SNP positions from filtered VCF
# ---------------------------------------------------------------------------


rule extract_snp_positions:
    input:
        str(FILT_DIR / "chr{chrom}.vcf.gz"),
    output:
        str(MAPS_DIR / "snps" / "chr{chrom}.bed.gz"),
    shell:
        """
        mkdir -p {MAPS_DIR}/snps
        bcftools query -f '%CHROM\t%POS0\t%POS\t%ID\n' {input} | bgzip > {output}
        """


# ---------------------------------------------------------------------------
# Step 2c: Interpolate deCODE genetic map
# ---------------------------------------------------------------------------


rule interpolate_map:
    input:
        snps=rules.extract_snp_positions.output[0],
        map=rules.convert_decode_map.output[0],
    output:
        str(MAPS_DIR / "interpolated" / "chr{chrom}.tab.gz"),
    shell:
        """
        mkdir -p {MAPS_DIR}/interpolated
        uv run ldetect interpolate-maps \
            --snp-file {input.snps} \
            --genetic-map {input.map} \
            --output {output} \
            --mode interval
        """


# Point-mode counterpart of interpolate_map, kept only so compare_maps_point
# can quantify the interpolation bug --mode interval fixed: feeding
# convert_decode_map.py's interval-rate rows into point-based interpolation
# uses the *next* interval's rate for SNPs in the *current* interval. See
# notes/logs/macdonald2022-interpolation-port.md.
rule interpolate_map_point:
    input:
        snps=rules.extract_snp_positions.output[0],
        map=rules.convert_decode_map.output[0],
    output:
        str(MAPS_DIR / "interpolated_point" / "chr{chrom}.tab.gz"),
    shell:
        """
        mkdir -p {MAPS_DIR}/interpolated_point
        uv run ldetect interpolate-maps \
            --snp-file {input.snps} \
            --genetic-map {input.map} \
            --output {output} \
            --mode point
        """


# ---------------------------------------------------------------------------
# Step 3: Run ldetect-lite per chromosome × block set
#
# The shared ldetect-lite covariance code applies the covariance cutoff before
# adding the diagonal shrinkage term, matching original ldetect behavior and
# preventing population-monomorphic variants from entering the covariance
# output. deCODE block sets use the interpolated deCODE map; pyrho block sets
# use MacDonald's published already-interpolated pyrho maps directly.
# ---------------------------------------------------------------------------


rule run_ldetect:
    input:
        genetic_map=block_set_genetic_map,
        vcf=block_set_filtered_vcf,
        individuals=lambda wc: str(
            RES_DIR / f"{block_set_population(wc.block_set)}_inds.txt"
        ),
    output:
        bed=str(RESULTS / "{block_set}" / "chr{chrom}" / "chr{chrom}-ld-blocks.bed"),
    log:
        run=str(RESULTS / "logs" / "{block_set}" / "chr{chrom}.ldetect.log"),
        timing=str(RESULTS / "logs" / "{block_set}" / "chr{chrom}.timing.log"),
    benchmark:
        str(RESULTS / "logs" / "{block_set}" / "chr{chrom}.benchmark.tsv")
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    threads: config["workers"]
    params:
        output_dir=str(RESULTS / "{block_set}" / "chr{chrom}"),
        log_dir=str(RESULTS / "logs" / "{block_set}"),
        ne=lambda wc: block_set_ne(wc.block_set),
        cutoff=config["cov_cutoff"],
        subset=config["subset"],
        n_snps_bw_bpoints=config["n_snps_bw_bpoints"],
    shell:
        """
        mkdir -p {params.log_dir}
        /usr/bin/time -v -o {log.timing} \
        uv run ldetect run \
            --genetic-map {input.genetic_map} \
            --reference-panel {input.vcf} \
            --individuals {input.individuals} \
            --chromosome chr{wildcards.chrom} \
            --output-dir {params.output_dir} \
            --ne {params.ne} \
            --cov-cutoff {params.cutoff} \
            --n-snps-bw-bpoints {params.n_snps_bw_bpoints} \
            --subset {params.subset} \
            --workers {threads} \
            > {log.run} 2>&1
        """


# ---------------------------------------------------------------------------
# Step 4: Post-processing
# ---------------------------------------------------------------------------


rule postprocess:
    input:
        bed=rules.run_ldetect.output.bed,
        vcf=block_set_filtered_vcf,
        centromeres=rules.download_centromeres.output[0],
    output:
        str(RESULTS / "{block_set}" / "chr{chrom}" / "chr{chrom}-ld-blocks.final.bed"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        min_snps=config["min_snps_per_block"],
        remove_centromeres=remove_centromeres_param,
    shell:
        """
        uv run python scripts/postprocess.py \
            --bed {input.bed} \
            --vcf {input.vcf} \
            --centromeres {input.centromeres} \
            {params.remove_centromeres} \
            --min-snps {params.min_snps} \
            --output {output}
        """


# ---------------------------------------------------------------------------
# Step 5: Combine chromosomes into genome-wide BED
# ---------------------------------------------------------------------------


rule combine_raw_beds:
    input:
        expand(
            str(RESULTS / "{{block_set}}" / "chr{chrom}" / "chr{chrom}-ld-blocks.bed"),
            chrom=CHROMS,
        ),
    output:
        str(RESULTS / "{block_set}_raw_LD_blocks.bed"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    shell:
        """
        # Write header once, then body lines from all chromosomes (skip per-file headers)
        echo -e "#chr\tstart\tstop" > {output}
        for f in {input}; do
            tail -n +2 "$f" >> {output}
        done
        """


rule combine_beds:
    input:
        expand(
            str(
                RESULTS
                / "{{block_set}}"
                / "chr{chrom}"
                / "chr{chrom}-ld-blocks.final.bed"
            ),
            chrom=CHROMS,
        ),
    output:
        str(RESULTS / "{block_set}_LD_blocks.bed"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    shell:
        """
        # Write header once, then body lines from all chromosomes (skip per-file headers)
        echo -e "#chr\tstart\tstop" > {output}
        for f in {input}; do
            tail -n +2 "$f" >> {output}
        done
        """


# ---------------------------------------------------------------------------
# Step 6a: Validate interpolated genetic maps
# ---------------------------------------------------------------------------


rule validate_maps:
    input:
        expand(
            str(MAPS_DIR / "interpolated" / "chr{chrom}.tab.gz"),
            chrom=CHROMS,
        ),
    output:
        str(RESULTS / "compare" / "map_summary.tsv"),
    shell:
        """
        uv run python scripts/validate_maps.py \
            --interpolated {input} \
            --output {output}
        """


rule compare_maps:
    input:
        expand(
            str(MAPS_DIR / "interpolated" / "chr{chrom}.tab.gz"),
            chrom=CHROMS,
        ),
    output:
        str(RESULTS / "compare" / "map_ref_comparison.tsv"),
    params:
        ref_dir=str(RES_DIR / "macdonalds_maps"),
    shell:
        """
        uv run python scripts/compare_maps.py \
            --ours {input} \
            --ref-dir {params.ref_dir} \
            --output {output}
        """


# Equivalence check: --mode interval (interpolate_map) is a direct port of
# MacDonald's R interpolation scripts, so it should match their published
# deCODE_interpolated_maps almost exactly (near-zero MAE). This rule reruns
# the same comparison against interpolate_map_point's --mode point output so
# the two map_ref_comparison*.tsv files can be diffed directly: interval
# mode's MAE should collapse to ~0 while point mode reproduces the
# previously measured 0.0013-0.0026 cM error. See
# notes/logs/macdonald2022-interpolation-port.md.
rule compare_maps_point_mode:
    input:
        expand(
            str(MAPS_DIR / "interpolated_point" / "chr{chrom}.tab.gz"),
            chrom=CHROMS,
        ),
    output:
        str(RESULTS / "compare" / "map_ref_comparison_point_mode.tsv"),
    params:
        ref_dir=str(RES_DIR / "macdonalds_maps"),
    shell:
        """
        uv run python scripts/compare_maps.py \
            --ours {input} \
            --ref-dir {params.ref_dir} \
            --output {output}
        """


# ---------------------------------------------------------------------------
# Step 6b: Compare LD blocks against MacDonald et al. published output
# ---------------------------------------------------------------------------


rule download_reference_bed:
    output:
        str(RES_DIR / "{reference_bed}"),
    params:
        base_url=config["reference_bed_base_url"],
    shell:
        """
        mkdir -p {RES_DIR}
        wget -q -O {output} \
            {params.base_url}/{wildcards.reference_bed}
        """


rule compare_blocks:
    input:
        ours=str(RESULTS / "{block_set}_LD_blocks.bed"),
        ref=lambda wc: str(RES_DIR / block_set_reference_bed(wc.block_set)),
    output:
        str(RESULTS / "compare" / "{block_set}_block_comparison.tsv"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        tolerance=50_000,
    shell:
        """
        uv run python scripts/compare_blocks.py \
            --ours {input.ours} \
            --ref {input.ref} \
            --output {output} \
            --tolerance {params.tolerance}
        """


rule compare_boundaries:
    input:
        ours=str(RESULTS / "{block_set}_LD_blocks.bed"),
        ref=lambda wc: str(RES_DIR / block_set_reference_bed(wc.block_set)),
    output:
        str(RESULTS / "compare" / "boundaries" / "{block_set}_boundary_offsets.tsv"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        tolerance=50_000,
    shell:
        """
        python scripts/compare_boundaries.py \
            --ours {input.ours} \
            --ref {input.ref} \
            --output {output} \
            --tolerance {params.tolerance}
        """


rule compare_raw_blocks:
    input:
        ours=str(RESULTS / "{block_set}_raw_LD_blocks.bed"),
        ref=lambda wc: str(RES_DIR / block_set_reference_bed(wc.block_set)),
    output:
        str(RESULTS / "compare" / "raw" / "{block_set}_block_comparison.tsv"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        tolerance=50_000,
    shell:
        """
        python scripts/compare_blocks.py \
            --ours {input.ours} \
            --ref {input.ref} \
            --output {output} \
            --tolerance {params.tolerance}
        """


rule compare_raw_boundaries:
    input:
        ours=str(RESULTS / "{block_set}_raw_LD_blocks.bed"),
        ref=lambda wc: str(RES_DIR / block_set_reference_bed(wc.block_set)),
    output:
        str(RESULTS / "compare" / "raw" / "{block_set}_boundary_offsets.tsv"),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        tolerance=50_000,
    shell:
        """
        python scripts/compare_boundaries.py \
            --ours {input.ours} \
            --ref {input.ref} \
            --output {output} \
            --tolerance {params.tolerance}
        """


# These diagnostic rules deliberately treat combined BEDs as snapshots rather
# than workflow inputs. Downloaded aggregate results can therefore be examined
# without Snakemake trying to reconstruct missing per-chromosome intermediates.
rule diagnose_boundaries:
    output:
        str(
            RESULTS
            / "compare"
            / "diagnostics"
            / "{block_set}_boundary_diagnostics.tsv"
        ),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        ours=str(RESULTS / "{block_set}_LD_blocks.bed"),
        ref=lambda wc: diagnostic_resource(
            block_set_reference_bed(wc.block_set)
        ),
        centromeres=diagnostic_resource("hg38_centromeres.txt.gz"),
        tolerance=50_000,
        window=2_000_000,
    shell:
        """
        python scripts/diagnose_boundaries.py \
            --ours {params.ours} \
            --ref {params.ref} \
            --centromeres {params.centromeres} \
            --output {output} \
            --tolerance {params.tolerance} \
            --window {params.window}
        """


rule diagnose_raw_boundaries:
    output:
        str(
            RESULTS
            / "compare"
            / "diagnostics"
            / "raw"
            / "{block_set}_boundary_diagnostics.tsv"
        ),
    wildcard_constraints:
        block_set=BLOCK_SET_PATTERN,
    params:
        ours=str(RESULTS / "{block_set}_raw_LD_blocks.bed"),
        ref=lambda wc: diagnostic_resource(
            block_set_reference_bed(wc.block_set)
        ),
        centromeres=diagnostic_resource("hg38_centromeres.txt.gz"),
        tolerance=50_000,
        window=2_000_000,
    shell:
        """
        python scripts/diagnose_boundaries.py \
            --ours {params.ours} \
            --ref {params.ref} \
            --centromeres {params.centromeres} \
            --output {output} \
            --tolerance {params.tolerance} \
            --window {params.window}
        """
