script:
    "../scripts/python/bedtools/functional/filter_cds.py"

shell:
    """
    complementBed -i {input.bed} -g {input.genome} | \
    intersectBed -a stdin -b {input.gapless} -sorted | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/gc/intersect_ranges.py"

shell:
    """
    gunzip -c {input.ref} | \
    sed 's/^[A-Za-z]\\+$/\\U&/' | \
    {input.bin} {wildcards.unit_len} {wildcards.total_len} - 2> {log} \
    > {output} 
    """

shell:
    "sed -n '/unit=[{wildcards.bases}]\+/p' {input} > {output}"

shell:
    "subtractBed -a {input.a} -b {input.b} > {output}"

shell:
    "subtractBed -a {input.a} -b {input.b} > {output}"

shell:
    """
    slopBed -i {input} -b 5 -g {params.genome} | \
    cut -f1-3 | \
    mergeBed -i stdin | \
    intersectBed -a stdin -b {params.gapless} -sorted -g {params.genome} | \
    bgzip -c \
    > {output}
    """

shell:
    """
    grep 'unit={wildcards.base}' {input} | \
    mergeBed -i stdin -d 1 | \
    awk '$3-$2>={wildcards.merged_len}' > \
    {output}
    """

shell:
    """
    multiIntersectBed -i {input} | \
    slopBed -i stdin -b 5 -g {params.genome} | \
    mergeBed -i stdin | \
    intersectBed -a stdin -b {params.gapless} -sorted | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/low_complexity/filter_sort_trf.py"

shell:
    """
    gunzip -c {input} | \
    mergeBed -i stdin | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/low_complexity/filter_sort_rmsk.py"

shell:
    """
    gunzip -c {input.rmsk} | \
    grep {wildcards.rmsk_class} | \
    mergeBed -i stdin | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/low_complexity/filter_sort_censat.py"

shell:
    "mergeBed -i {input} | bgzip -c > {output}"

shell:
    """
    slopBed -i {input.bed} -b 5 -g {input.genome} | \
    mergeBed -i stdin | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

shell:
    """
    complementBed -i {input.bed} -g {params.genome} |
    intersectBed -a stdin -b {params.gapless} -sorted | \
    bgzip -c > {output}
    """

shell:
    """
    slopBed -i {input.perfect} -b 5 -g {input.genome} | \
    mergeBed -i stdin | \
    multiIntersectBed -i stdin {input.imperfect} | \
    mergeBed -i stdin | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c \
    > {output}
    """

shell:
    """
    multiIntersectBed -i {input.beds} | \
    slopBed -i stdin -b 5 -g {input.genome} | \
    mergeBed -i stdin > {output}
    """

shell:
    """
    awk '({params.lower}+10)<=($3-$2) && ($3-$2)<({params.upper}+10)' {input.tr} | \
    subtractBed -a stdin -b {input.hp} | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

shell:
    """
    multiIntersectBed -i {input} | \
    mergeBed -i stdin | \
    bgzip -c > {output}
    """

shell:
    "curl -sS -L -o {output} {params.url}"

shell:
    """
    mkdir -p {config.tools_bin_dir} && \
    tar xjf {input} \
    --directory {config.tools_bin_dir} \
    --strip-components=2 \
    {params.bins}
    """

script:
    "../scripts/python/bedtools/mappability/filter_ref.py"

shell:
    """
    PATH={config.tools_bin_dir}:$PATH
    {input.bin} \
    --complement emulate \
    -T {threads} \
    -i {input.fa} \
    -o {params.base} > {log} 2>&1
    """

shell:
    """
    {input.bin} \
    -m {wildcards.m} \
    -e {wildcards.e} \
    -l {wildcards.l} \
    -T {threads} \
    -I {input.fa} \
    -o {params.base} > {log} 2>&1
    """

shell:
    """
    {input.bin} \
    -I {input.idx} \
    -i {input.map} \
    -o {params.base} 2> {log}
    """

shell:
    """
    sed 's/ AC//' {input} | \
    wig2bed -d | \
    awk '$5>0.9' | \
    cut -f1-3 | \
    gzip -c > \
    {output}
    """

script:
    "../scripts/python/bedtools/mappability/merge_nonunique.py"

shell:
    """
    complementBed -i {input} -g {params.genome} | \
    intersectBed -a stdin -b {params.gapless} -sorted -g {params.genome} | \
    bgzip -c > {output}
    """

shell:
    """
    complementBed -i {input.gapless} -g {input.genome} | \
    slopBed -i stdin -b 15000 -g {input.genome} | \
    mergeBed -i stdin | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/otherdifficult/filter_vdj.py"

shell:
    """
    mergeBed -i {input.bed} | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/other/filter_sort_other.py"

script:
    "../scripts/python/bedtools/postprocess/list_strats.py"

script:
    "../scripts/python/bedtools/postprocess/list_md5.py"

script:
    "../scripts/python/bedtools/postprocess/generate_tsv.py"

script:
    "../scripts/python/bedtools/postprocess/run_unit_tests.py"

shell:
    """
    curl -fsSq -o {params.tar} {params.url} && \
    mkdir {output} && \
    tar xzf {params.tar} --directory {output} --strip-components=1
    """

script:
    "../scripts/python/bedtools/postprocess/diff_previous_strats.py"

run:
    with open(output[0], "w") as f:
        for ref_key, build_key in zip(rks, bks):
            for i in config.buildkey_to_chr_indices(ref_key, build_key):
                pattern = config.refkey_to_final_chr_pattern(ref_key)
                line = [
                    str(i.value),
                    f"{ref_key}@{build_key}",
                    i.chr_name_full(pattern),
                ]
                f.write("\t".join(line) + "\n")

script:
    "../scripts/rmarkdown/rmarkdown/validate.Rmd"

shell:
    "gunzip -c {input} > {output}"

script:
    "../scripts/python/bedtools/postprocess/run_happy.py"

script:
    "../scripts/rmarkdown/rmarkdown/benchmark.Rmd"

shell:
    """
    cp {input.main} {output.main}
    cp {input.validation} {output.validation}
    """

shell:
    """
    tar czf {output} -C {params.parent} {params.target}
    """

shell:
    """
    find {params.root} -type f -exec md5sum {{}} + | \
    sed 's|{params.root}|\.|' | \
    grep -v "\./genome-stratifications-md5s\.txt" | \
    sort -k2,2 > {output} && \
    cd {params.root} && \
    md5sum -c --quiet genome-stratifications-md5s.txt
    """

script:
    "../scripts/python/bedtools/xy/write_par.py"

script:
    "../scripts/python/bedtools/misc/get_file.py"

shell:
    """
    gunzip -c {input} | \
    samtools faidx - -o - \
    2> {log} > {output}
    """

script:
    "../scripts/python/bedtools/ref/get_genome.py"

shell:
    """
    samtools faidx {input.fa} $(cut -f1 {input.genome} | tr '\n' ' ') 2> {log} | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/ref/get_gapless.py"

script:
    "../scripts/python/bedtools/ref/get_ftbl_mapper.py"

shell:
    """
    gunzip -c {input} | \
    grep -v '^#' | \
    awk -F "\t" 'OFS="\t" {{ print $1,$4,$5,$2,$3,$9 }}' | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/ref/filter_sort_bench_bed.py"

shell:
    "curl -sS -L -o {output} {params.url}"

shell:
    """
    mkdir {output} && \
    tar xzf {input} --directory {output} --strip-components=1
    """

shell:
    "make -C {input} 2>&1 > {log} && mv {input}/repseq {output}"

script:
    "../scripts/python/bedtools/segdups/filter_sort_superdups.py"

shell:
    """
    mergeBed -i {input.bed} -d 100 | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

shell:
    """
    gunzip -c {input} | \
    awk '($3-$2 > 10000)' | \
    bgzip -c > {output}
    """

shell:
    """
    complementBed -i {input.bed} -g {input.genome} | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

shell:
    """
    seqtk telo {input.ref} 2> {log} | \
    cut -f1-3 | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c \
    > {output}
    """

shell:
    """
    multiIntersectBed -i {input} | \
    mergeBed -i stdin | \
    bgzip -c \
    > {output}
    """

shell:
    """
    complementBed -i {input.bed} -g {params.genome} |
    intersectBed -a stdin -b {params.gapless} -sorted -g {params.genome} | \
    bgzip -c \
    > {output}
    """

shell:
    """
    intersectBed -a {input.bed} -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

script:
    "../scripts/python/bedtools/xy/filter_sort_features.py"

shell:
    """
    complementBed -i {input.bed} -g {input.genome} | \
    grep {wildcards.sex_chr} | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

shell:
    """
    awk -v OFS='\t' {{'print $1,\"0\",$2'}} {input.bed} | \
    grep -v \"X\|Y\" | \
    intersectBed -a stdin -b {input.gapless} -sorted -g {input.genome} | \
    bgzip -c > {output}
    """

from pathlib import Path
from typing import Any, NamedTuple, Callable
import subprocess as sp
import common.config as cfg
import json

# use subprocess here because I don't see an easy way to stream a bedtools
# python object to a bgzip file (it only seems to do gzip...oh well)


class GCInput(NamedTuple):
    bed: Path
    fraction: int
    is_range_bound: bool


def write_simple_range_beds(
    final_path: Callable[[str], Path],
    gs: list[GCInput],
    is_low: bool,
) -> list[str]:
    def fmt_out(bigger_frac: int, smaller_frac: int) -> Path:
        lower_frac, upper_frac = (
            (smaller_frac, bigger_frac) if is_low else (bigger_frac, smaller_frac)
        )
        return final_path(f"gc{lower_frac}to{upper_frac}_slop50")

    pairs = zip(gs[1:], gs[:-1]) if is_low else zip(gs[:-1], gs[1:])
    torun = [
        (fmt_out(bigger.fraction, smaller.fraction), bigger, smaller)
        for bigger, smaller in pairs
    ]
    for out, bigger, smaller in torun:
        with open(out, "wb") as f:
            p0 = sp.Popen(
                ["subtractBed", "-a", bigger.bed, "-b", smaller.bed],
                stdout=sp.PIPE,
            )
            p1 = sp.run(["bgzip", "-c"], stdin=p0.stdout, stdout=f)
            p0.wait()
            if not (p0.returncode == p1.returncode == 0):
                exit(1)
    return [str(t[0]) for t in torun]


def write_middle_range_bed(
    final_path: Callable[[str], Path],
    lower: GCInput,
    upper: GCInput,
    genome: Path,
    gapless: Path,
) -> str:
    out = final_path(f"gc{lower.fraction}to{upper.fraction}_slop50")
    with open(out, "wb") as f:
        p0 = sp.Popen(
            ["complementBed", "-i", upper.bed, "-g", genome],
            stdout=sp.PIPE,
        )
        p1 = sp.Popen(
            ["subtractBed", "-a", "stdin", "-b", lower.bed],
            stdin=p0.stdout,
            stdout=sp.PIPE,
        )
        p2 = sp.Popen(
            ["intersectBed", "-a", "stdin", "-b", gapless, "-sorted", "-g", genome],
            stdin=p1.stdout,
            stdout=sp.PIPE,
        )
        p3 = sp.run(["bgzip", "-c"], stdin=p2.stdout, stdout=f)
        p0.wait()
        p1.wait()
        p2.wait()
        if not (p0.returncode == p1.returncode == p2.returncode == p3.returncode == 0):
            exit(1)
    return str(out)


def write_intersected_range_beds(
    final_path: Callable[[str], Path],
    low: list[GCInput],
    high: list[GCInput],
) -> list[str]:
    pairs = zip(
        [x for x in low if x.is_range_bound],
        [x for x in reversed(high) if x.is_range_bound],
    )
    torun = [
        (i, final_path(f"gclt{b1.fraction}orgt{b2.fraction}_slop50"), b1, b2)
        for i, (b1, b2) in enumerate(pairs)
    ]
    for i, bed_out, b1, b2 in torun:
        with open(bed_out, "wb") as f:
            p0 = sp.Popen(
                ["multiIntersectBed", "-i", b1.bed, b2.bed],
                stdout=sp.PIPE,
            )
            p1 = sp.Popen(
                ["mergeBed", "-i", "stdin"],
                stdin=p0.stdout,
                stdout=sp.PIPE,
            )
            p2 = sp.run(["bgzip", "-c"], stdin=p1.stdout, stdout=f)
            p0.wait()
            p1.wait()
            if not (p0.returncode == p1.returncode == p2.returncode == 0):
                exit(1)

    return [str(t[1]) for t in torun]


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards.ref_key)
    bk = cfg.BuildKey(smk.wildcards.build_key)
    gps = sconf.buildkey_to_include(rk, bk).gc
    assert gps is not None, "this should not happen"
    # ASSUME both of these input lists are sorted by GC fraction
    low = [GCInput(p, f, r) for p, (f, r) in zip(smk.input.low, gps.low)]
    high = [GCInput(p, f, r) for p, (f, r) in zip(smk.input.high, gps.high)]
    genome = Path(smk.input.genome)
    gapless = Path(smk.input.gapless)

    def final_path(name: str) -> Path:
        p = Path(str(smk.params.path_pattern).format(name))
        p.parent.mkdir(exist_ok=True, parents=True)
        return p

    low_strats = write_simple_range_beds(final_path, low, True)
    high_strats = write_simple_range_beds(final_path, high, False)
    range_strat = write_middle_range_bed(
        final_path,
        low[-1],
        high[0],
        genome,
        gapless,
    )
    inter_strats = write_intersected_range_beds(
        final_path,
        low,
        high,
    )

    with open(smk.output[0], "w") as f:
        # put the first low and last high input here since these are already
        # in the final directory
        obj = {
            "gc_ranges": [
                low[0][0],
                *low_strats,
                range_strat,
                *high_strats,
                high[-1][0],
            ],
            "widest_extreme": inter_strats[0],
            "other_extremes": inter_strats[1:],
        }
        json.dump(obj, f, indent=4)


main(snakemake, snakemake.config)  # type: ignore

from typing import Any
import common.config as cfg
from common.bed import read_bed, filter_sort_bed, write_bed


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards["ref_key"])
    bk = cfg.BuildKey(smk.wildcards["build_key"])

    bedfile = sconf.refkey_to_strat(rk).low_complexity.satellites
    assert bedfile is not None, "this should not happen"

    ps = bedfile.params
    conv = sconf.buildkey_to_chr_conversion(rk, bk, ps.chr_pattern)

    df = read_bed(smk.input[0], ps, [bedfile.sat_col])
    df = filter_sort_bed(conv, df)
    df = df[~df[3].str.startswith("ct_")]
    write_bed(smk.output[0], df)


main(snakemake, snakemake.config)  # type: ignore

from typing import Any
import common.config as cfg
from common.bed import read_filter_sort_bed


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards["ref_key"])
    bk = cfg.BuildKey(smk.wildcards["build_key"])
    rm = sconf.refkey_to_strat(rk).low_complexity.rmsk
    assert rm is not None, "this should not happen"
    read_filter_sort_bed(
        sconf, smk.input[0], smk.output[0], rm.params, rk, bk, [rm.class_col]
    )


main(snakemake, snakemake.config)  # type: ignore

from typing import Any
import common.config as cfg
from common.bed import read_filter_sort_bed


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards["ref_key"])
    bk = cfg.BuildKey(smk.wildcards["build_key"])
    ss = sconf.refkey_to_strat(rk).low_complexity.simreps
    assert ss is not None, "this should not happen"
    read_filter_sort_bed(sconf, smk.input[0], smk.output[0], ss.params, rk, bk)


main(snakemake, snakemake.config)  # type: ignore

import json
from typing import Any
import common.config as cfg
from common.bed import filter_sort_bed_inner, read_bed, write_bed

VDJ_PAT = "^ID=gene-(IGH|IGK|IGL|TRA|TRB|TRG);"


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards["ref_key"])
    bk = cfg.BuildKey(smk.wildcards["build_key"])
    to_map = sconf.buildkey_to_final_chr_mapping(rk, bk)

    with open(smk.input["mapper"], "r") as f:
        from_map = json.load(f)

    df = read_bed(smk.input["bed"], more=[5])
    df = df[df[3].str.match(VDJ_PAT)][[0, 1, 2]].copy()
    df = filter_sort_bed_inner(from_map, to_map, df)

    write_bed(smk.output[0], df)


main(snakemake, snakemake.config)  # type: ignore

from typing import Any
import common.config as cfg
from common.bed import read_bed, filter_sort_bed, write_bed
from pybedtools import BedTool as bt  # type: ignore


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards["ref_key"])
    bk = cfg.BuildKey(smk.wildcards["build_key"])
    lk = cfg.OtherLevelKey(smk.wildcards["other_level_key"])
    sk = cfg.OtherStratKey(smk.wildcards["other_strat_key"])
    bed = sconf.otherkey_to_bed(rk, bk, lk, sk)
    ps = bed.params
    conv = sconf.buildkey_to_chr_conversion(rk, bk, ps.chr_pattern)
    df = filter_sort_bed(conv, read_bed(smk.input["bed"][0], ps))
    if bed.remove_gaps:
        df = (
            bt()
            .from_dataframe(df)
            .intersect(
                smk.input["gapless"],
                sorted=True,
                g=smk.input["genome"][0],
            )
            .to_dataframe()
        )
    write_bed(smk.output[0], df)


main(snakemake, snakemake.config)  # type: ignore

from pathlib import Path
from typing import Any
import common.config as cfg
from common.stratdiff.lib.diff import compare_all
from common.io import setup_logging


log = setup_logging(snakemake.log[0])  # type: ignore


def main(smk: Any, sconf: cfg.GiabStrats) -> None:
    rk = cfg.RefKey(smk.wildcards.ref_key)
    bk = cfg.BuildKey(smk.wildcards.build_key)
    comparison = sconf.buildkey_to_build(rk, bk).comparison
    pattern = sconf.refkey_to_final_chr_pattern(rk)
    ixs = sconf.buildkey_to_chr_indices(rk, bk)

    assert comparison is not None, "this should not happen"

    outdir = Path(smk.output[0]).parent

    logged = compare_all(
        Path(smk.input["new_list"]).parent,
        smk.input["old"],
        outdir,
        comparison.path_mapper,
        comparison.replacements,
        [i.chr_name_full(pattern) for i in ixs],
        comparison.ignore_generated,
        comparison.ignore_other,
    )
    for x in logged:
        log.info(x)


main(snakemake, snakemake.config)  # type: ignore

from pathlib import Path
from typing import Any
import re


def main(smk: Any) -> None:
    with open(smk.input[0], "r") as i, open(smk.output[0], "w") as o:
        for f in i:
            fp = Path(f.strip())
            level = re.sub("^[^_]+_", "", fp.name.replace(".bed.gz", ""))
            o.write(f"{level}\t{Path(fp.parent.name) / fp.name}\n")


main(snakemake)  # type: ignore

from typing import Any
from os.path import dirname
from common.io import get_md5


def strat_files(path: str) -> list[str]:
    with open(path, "r") as f:
        return [s.strip() for s in f]


def main(smk: Any) -> None:
    out = str(smk.output)
    root = dirname(out)
    ss = strat_files(str(smk.input))
    with open(out, "w") as op:
        for s in ss:
            h = get_md5(s)
            p = s.replace(root + "/", "")
            op.write(f"{h}  {p}\n")


main(snakemake)  # type: ignore

from typing import Any


def main(smk: Any) -> None:
    with open(smk.output[0], "w") as f:
        for s in sorted(smk.input):
            f.write(s + "\n")


main(snakemake)  # type: ignore

import gzip
import os
import subprocess as sp
from typing import Any
from pathlib import Path


def gzip_empty(path: str) -> bool:
    with gzip.open(path, "rb") as f:
        return len(f.read(1)) == 0


def main(smk: Any) -> None:
    # ASSUME if filtered benchmark bed is empty then we have nothing to
    # benchmark because the desired chromosomes in this build are not included
    # in the benchmark
    with open(smk.log[0], "w") as f:
        if gzip_empty(smk.input["bench_bed"][0]):
            Path(smk.output[0]).touch()
            f.write("no overlapping chromosomes, writing empty dummy file\n")
        else:
            cmd = [
                "hap.py",
                *["--engine", "vcfeval"],
                "--verbose",
                *["--threads", str(smk.threads)],
                *["--stratification", smk.input["strats"][0]],
                *["-f", smk.input["bench_bed"][0]],
                *["-o", smk.params["prefix"]],
                smk.input["bench_vcf"][0],
                smk.input["query_vcf"][0],
            ]
            res = sp.run(
                cmd,
                stderr=sp.STDOUT,
                stdout=f,
                env={"HGREF": smk.input["ref"][0], **os.environ},
            )
            if res.returncode != 0:
                exit(1)


main(snakemake)  # type: ignore