script:
    "../scripts/amptk_trim_paired.py"

script:
    "../scripts/amptk_cluster_paired.py"

shell:
    """
    exec 1> "{output}"
    echo "Read merging and primer trimming"
    echo "================================"
    cat {input.merge_trim:q}
    printf "\n\n\n"
    echo "Denoising"
    echo "=========="
    cat {input.cluster:q}
    """

shell:
    """
    gzip -dc {input.otutab} > {output.tmp_tab}
    biom convert -i {input.otutab} \
      -o {output.biom} \
      --table-type 'OTU table' --to-json &> {log}
    """

shell:
    """
    biom convert -i {input.biom}  \
      -o {output.biom_hdf5} \
      --table-type "OTU table" --to-hdf5 &> {log}
    """

shell:
    "rm -Rf results/*/workflow_*/*/*/cmp"

script:
    "../scripts/taxonomy/idtaxa_train.R"

script:
    "../scripts/taxonomy/idtaxa_assign.R"

shell:
    "rm -Rf results/*/workflow_*/*/*/ITSx"

script:
    "../scripts/combine_sample_reports.py"

shell:
    """
    exec &> {log}
    set -xeuo pipefail
    if [ "{wildcards.layout}" = "single.rev" ]; then
        vsearch --fastx_revcomp "{input.otus}" --fastaout "{output.otus}"
    else
        ln -srf "{input.otus}" "{output.otus}"
    fi        
    """

script:
    "../scripts/collect_sample_lists.py"

script:
    "../scripts/dump_config.py"

script:
    "../scripts/prepare_primers.py"

script:
    "../scripts/collect_sample_files.py"

script:
    "../scripts/make_pooling_list.py"

script:
    "../scripts/pool_raw.py"

shell:
    """
    mkdir -p {output}
    ln -sr {input.run_dir} {output}/nested
    """

script:
    "../scripts/make_new_sample_tab.py"

script:
    "../scripts/list_samples_yaml.py"

script:
    "../scripts/collect_unique_files.py"

script:
    "../scripts/link_data_dir.py"

shell:
    """
    ln -srf {input.existing} {output} 2> {log}
    """

shell:
    """
    outdir="$(dirname "{output}")"
    multiqc -fm fastqc -o "$outdir" "{input.fastqc}" &> {log}
    (cd "$outdir" && zip -FSqr -rm multiqc_data.zip multiqc_data) 2> {log}
    """

script:
    "../scripts/make_new_sample_tab.py"

script:
    "../scripts/qiime_trim_single.py"

script:
    "../scripts/qiime_trim_paired.py"

script:
    "../scripts/qiime_dada2.py"

shell:
    """
    exec &> {log}
    # export table
    mkdir -p {output.tmp}
    tab={output.tab}
    tab=${{tab%.gz}}
    qiime tools export \
        --input-path {input.tab} \
        --output-path {output.tmp}
    mv {output.tmp}/feature-table.biom {output.biom_hdf5}
    biom convert -i {output.biom_hdf5}  \
        -o {output.biom_json} \
        --to-json --table-type "OTU table"
    biom convert -i {output.biom_hdf5}  \
        -o $tab \
        --to-tsv --table-type "OTU table"
    sed -i '1,1d' $tab
    gzip -nf $tab

    # export seqs
    qiime tools export \
        --input-path {input.clustered} \
        --output-path {output.tmp}
    cat {output.tmp}/dna-sequences.fasta > {output.clustered}

    # export stats
    qiime metadata tabulate \
        --m-input-file {input.stats} \
        --o-visualization {output.tmp}/stats.qzv
    qiime tools export \
        --input-path {output.tmp}/stats.qzv \
        --output-path {output.tmp}/stats
    statfile={output.tmp}/stats/metadata.tsv
    head -n1 $statfile > {output.stats}
    tail -n+3 $statfile >> {output.stats}
    """

script:
    "../scripts/combine_sample_reports.py"

shell:
    """
    exec 1> "{output}"
    echo "Primer trimming"
    echo "==============="
    cat {input.trim:q}
    printf "\n\n\n"
    echo "Denoising"
    echo "=========="
    cat {input.cluster:q}
    """

shell:
    """
    exec &> {log}
    par=$(cat {input.params)
    if [ "$par" != "{}" ]; then
        echo "Unknown QIIME naive-bayes classifier config: $par"
        exit 1
    fi
    mkdir -p {output.tmp}
    zstd -dcqf {input.seq} > {output.tmp}/seq.fasta 2> {log}
    # extract taxonomic lineages from FASTA for import in QIIME
    zstd -dcqf {input} | st . --to-tsv id,desc > {output.tmp}/tax.txt
    qiime tools import \
        --type 'FeatureData[Sequence]' \
        --input-path {output.tmp}/seq.fasta \
        --output-path {output.tmp}/seq.qza
    qiime tools import \
        --type 'FeatureData[Taxonomy]' \
        --input-format HeaderlessTSVTaxonomyFormat \
        --input-path {output.tmp}/tax.txt \
        --output-path {output.tmp}/tax.qza
    qiime feature-classifier fit-classifier-naive-bayes \
        --i-reference-reads {output.tmp}/seq.qza \
        --i-reference-taxonomy {output.tmp}/tax.qza \
        --o-classifier {output.tmp}/trained.qza
    zstd -dcq {output.tmp}/trained.qza > {output.trained}
    """

shell:
    """
    exec &> {log}
    mkdir -p {output.tmp}
    # lowercase letters cause problems -> convert to uppercase
    # (cannot use 'st upper' because seqtool is not in conda environment,
    # therefore using awk)
    awk '/^>/ {{print($0)}}; /^[^>]/ {{print(toupper($0))}}' {input.seq} > {output.tmp}/input.fasta
    qiime tools import \
        --input-path {output.tmp}/input.fasta \
        --type 'FeatureData[Sequence]' \
        --input-format DNAFASTAFormat \
        --output-path {output.tmp}/seqs.qza

    qiime feature-classifier classify-sklearn \
        --i-classifier {input.db} \
        --i-reads {output.tmp}/seqs.qza \
        --o-classification {output.tmp}/classified.qza \
        --p-reads-per-batch 1000 \
        --p-confidence {params.par[confidence]} \
        --p-n-jobs {threads} \
        --verbose

    qiime tools export \
    --input-path {output.tmp}/classified.qza \
    --output-path {output.tmp}

    gzip -nc {output.tmp}/taxonomy.tsv > {output.tax}
    """

script:
    "../scripts/taxonomy/write_db_config.py"

script:
    "../scripts/taxonomy/obtain_taxdb.py"

script:
    "../scripts/taxonomy/write_db_config.py"

script:
    "../scripts/taxonomy/filter_taxdb.py"

script:
    "../scripts/taxonomy/write_db_config.py"

script:
    "../scripts/taxonomy/obtain_taxdb_formatted.py"

shell:
    """
    zstd -dcqf "{input.db}" > "{output.db}"
    """

script:
    "../scripts/taxonomy/rank_propagate.py"

shell:
    """
    tax={input.tax}
    st set -ul <(gzip -dc "$tax") -d {{l:2}} "{input.fa}" |
      st replace -d '__' ':' |
      st replace -dr ' *; *' ' ' |
      gzip -nc > {output}
    """

shell:
    """
    exec &> {log}
    set -xeuo pipefail
    mkdir -p $(dirname {output.tax_tmp})
    gzip -dc {input.tax} | 
    sed 's/Taxon/taxonomy/g' |
    sed 's/Feature ID/# Feature ID/g' > {output.tax_tmp}
    if [[ $(wc -l < "{output.tax_tmp}") -ge 2 ]]; then
            biom add-metadata -i {input.biom}  \
                -o /dev/stdout \
                --observation-metadata-fp {output.tax_tmp} \
                --sc-separated taxonomy --float-fields Confidence --output-as-json |
            gzip -nc > {output.biom}
            biom add-metadata -i {input.biom_hdf5}  \
                -o {output.biom_tmp} \
                --observation-metadata-fp {output.tax_tmp} \
                --sc-separated taxonomy --float-fields Confidence
            gzip -nc {output.biom_tmp} > {output.biom_hdf5}
    else
        # no taxa
        echo -n | gzip -nc > {output.biom}
        echo -n | gzip -nc > {output.biom_hdf5}
    fi
    """

shell:
    "rm -Rf results/*/workflow_*/*/*/taxonomy"

script:
    "../scripts/uvsnake_gen_config.py"

script:
    "../scripts/uvsnake_run.py"

script:
    "../scripts/uvsnake_run.py"

shell:
    """
    indir="$(dirname "{input.results[0]}")"
    outdir="$(dirname "{output.results[0]}")"
    cp -f "$indir/{wildcards.cluster_method}.fasta" "$outdir/clusters.fasta"
    cp -f "$indir/{wildcards.cluster_method}_otutab.txt.gz" "$outdir/otutab.txt.gz"
    cp -f "$indir/{wildcards.cluster_method}.biom" "$outdir/otutab.biom"
    cp -f "{input.stats}" "{output.stats}"
    cp -f "{input.log}" "{output.log}"
    """

script:
    "../scripts/taxonomy/convert_utax.py"

import os
from os.path import exists, abspath, relpath
import shutil
import sys
from typing import *

from utils import file_logging
from utils.sample_list import SampleList


def link_samples(sample_tab, outdir):
    # delete output dir to make sure there are no orphan files
    # if exists(outdir):
    #     shutil.rmtree(outdir)
    os.makedirs(outdir, exist_ok=True)
    # link
    sample_list = SampleList(sample_tab)
    for sample, read_paths in sample_list.samples():
        # R1 and R2 files
        for i, source in enumerate(read_paths):
            target = os.path.join(outdir, f"{sample}_R{i+1}.fastq.gz")
            os.symlink(abspath(source), abspath(target))
            # report
            source = relpath(source, ".")
            target = relpath(target, ".")
            print("{} > {}".format(source, target), file=sys.stdout)


with file_logging(snakemake.log[0]):
    link_samples(snakemake.input.sample_tab, snakemake.output.sample_dir)

from utils import file_logging
from utils.sample_list import SampleList


def collect_sample_lists(run_meta, reserved_chars, path_template):
    for d in run_meta:
        tsv_out = path_template.format(ext="tsv", **d)
        l = SampleList(d["sample_file"], reserved_chars=reserved_chars)
        with open(tsv_out, "w") as o:
            l.write(o, sort_by_sample=True)
        yaml_out = path_template.format(ext="yaml", **d)
        with open(yaml_out, "w") as o:
            l.write_yaml(o)


with file_logging(snakemake.log[0]):
    collect_sample_lists(**snakemake.params)

from collections import defaultdict
import csv
import hashlib
import os
from os.path import exists, abspath, relpath
import shutil
import sys
from typing import *

from utils import file_logging
from utils.sample_list import SampleList


def get_samples(run_meta, path_template):
    # get all necessary metadata
    _run_meta = [
        (
            (d["technology"], d["layout"], d["run"]),
            list(SampleList(path_template.format(**d).samples()))
        )
        for d in run_meta
    ]
    # get list of runs per sample
    sample2run = defaultdict(set)
    for run, samples in _run_meta:
        for sample in samples:
            sample2run[sample].add(run)
    # Which layout/run/technology combinations have duplicate samples?
    # Runs/layouts without duplicates are set to None
    dupes = sorted(set(
        run
        for runs in sample2run.values() if len(runs) > 1
        for run in runs
    ))
    # Define unique suffixes for run/layout combination that have duplicate samples.
    # For consistency, all sample names in these runs will receive a suffix, 
    # even if some of them are not duplicated.
    suffixes = {r: dupes.index(dupes[r]) + 1 for r in dupes}
    sample_dict = {}
    for run, samples in _run_meta:
        try:
            suffix = suffixes[run]
        except KeyError:
            suffix = None
        for sample, reads in samples:
            unique_sample = sample if suffix is None else f"{sample}_{suffix}"
            assert not unique_sample in sample_dict
            paths = [(path, f"{unique_sample}_R{i+1}.fastq.gz")
                     for i, path in enumerate(reads)]
            sample_dict[unique_sample] = (sample, run, paths)
    # convert to flat sorted list
    samples = [(unique_sample, *other) for unique_sample, other in sample_dict.items()]
    return sorted(samples)


def file_md5(filename):
    md5 = hashlib.md5()
    with open(filename, "rb") as f:
        # Read and update hash in chunks of 4K
        for chunk in iter(lambda: f.read(4096), b""):
            md5.update(chunk)
        return md5.hexdigest()


def do_link(path_iter, outdir):
    for paths in path_iter:
        for source, target in paths:
            target = os.path.join(outdir, target)
            assert not os.path.exists(target)
            os.symlink(abspath(source), abspath(target))
            # report
            source = relpath(source, ".")
            target = relpath(target, ".")
            print("{} > {}".format(source, target), file=sys.stdout)


def dump_tsv(samples, outfile):
    with open(outfile, "w") as o:
        w = csv.writer(o, delimiter="\t")
        header = ["technology", "run", "layout", 
                  "sample", "unique_sample",
                  "source_read_1", "source_read_2",
                  "read_1", "read_2",
                  "read_1_md5", "read_2_md5"]
        w.writerow(header)
        for unique_sample, sample, meta, paths in samples:
            orig_reads, target_files = zip(*paths)
            orig_reads = list(orig_reads)
            target_files = list(target_files)
            md5 = [file_md5(p) for p in orig_reads]
            if len(orig_reads) == 1:
                orig_reads.append("")
                target_files.append("")
                md5.append("")
            row = list(meta) + [sample, unique_sample] + orig_reads + target_files + md5
            w.writerow(row)


def link_samples(run_meta, read_dir, sample_file):
    # delete output dir to make sure there are no orphan files
    if exists(read_dir):
        shutil.rmtree(read_dir)
    if not exists(read_dir):
        os.makedirs(read_dir)
    # get the sample paths
    samples = get_samples(run_meta)
    # do the work
    do_link((paths for _, _, _, paths in samples), read_dir)
    dump_tsv(samples, sample_file)


with file_logging(snakemake.log[0]):
    link_samples(snakemake.params.run_meta, snakemake.input.path_template, 
                 snakemake.output.read_dir, snakemake.output.tsv)

from collections import OrderedDict
from copy import copy
import csv
import re

from utils import file_logging


class SampleReport(object):
    def __init__(self, path, pattern):
        # read data
        params = pattern.parse(path)
        self.group_header = list(params)
        self.groups = list(params.values())
        with open(path) as f:
            rdr = csv.reader(f, delimiter="\t")
            try:
                self.header = next(rdr)
            except StopIteration:
                self.header = []
            self.rows = [row for row in rdr if row]

class PathPattern(object):
    def __init__(self, pattern):
        self._pattern = re.compile(pattern)

    def parse(self, path):
        out = OrderedDict()
        for m in re.finditer(self._pattern, path):
            out.update(m.groupdict())
        return out


def combine_reports(sample_reports, path_pattern, outfile):
    pattern = PathPattern(path_pattern)
    sample_files = [SampleReport(path, pattern) for path in sample_reports]
    assert all(f.group_header == sample_files[0].group_header for f in sample_files[1:])
    # obtain header
    header = []
    for f in sample_files:
        f.field_idx = []
        for field in f.header:
            try:
                f.field_idx.append(header.index(field))
            except ValueError:
                f.field_idx.append(len(header))
                header.append(field)
    # write data
    with open(outfile, "w") as o:
        w = csv.writer(o, delimiter="\t")
        if sample_files:
            w.writerow(sample_files[0].group_header + header)
            empty_row = [""] * len(header)
            for f in sample_files:
                for row in f.rows:
                    out = copy(empty_row)
                    for i, field in zip(f.field_idx, row):
                        out[i] = field
                    w.writerow(f.groups + out)


with file_logging(snakemake.log[0]):
    combine_reports(
        snakemake.input.reports,
        snakemake.params.path_pattern,
        snakemake.output.report
    )

import yaml
from copy import deepcopy

from utils import file_logging


def dump_config(config, outfile):
    with open(outfile, "w") as o:
        c = deepcopy(config)  # make copy since we modify it
        del c["settings"]["input"]
        del c["settings"]["primers"]
        del c["settings"]["taxonomy_db_sources"]
        del c["settings"]["taxonomy_dbs"]
        del c["settings"]["taxonomy_methods"]
        c["taxonomy"] = {
            marker: {"-".join(name): config for name, config in tax.items()}
            for marker, tax in c["taxonomy"].items()
        }
        yaml.safe_dump(c, o)


with file_logging(snakemake.log[0]):
    dump_config(snakemake.params.config, snakemake.output[0])

import os
from os.path import dirname

from utils import file_logging


def link_data_dir(outdir_list, clust_files, data_dir):
    # write output dirs
    outdirs = [dirname(f) for f in clust_files]
    with open(outdir_list, "w") as o:
        o.writelines((l + "\n" for l in outdirs))
    if len(outdirs) == 1:
        source_dir = outdirs[0]
        if os.path.exists(data_dir):
            os.remove(data_dir)
        os.symlink(os.path.relpath(source_dir, dirname(data_dir)), data_dir)


with file_logging(snakemake.log[0]):
    link_data_dir(outdir_list=snakemake.output.out_list,
                  clust_files=snakemake.input.clust,
                  data_dir=snakemake.params.data_dir)

from collections import defaultdict
import os

import yaml

from utils import file_logging


def list_samples(sample_files, outfile):
    out = defaultdict(dict)
    for sample_file in sample_files:
        p = sample_file.split(os.sep)
        layout = p[-3]
        run = p[-2]
        with open(sample_file) as f:
            out[run][layout] = yaml.safe_load(f)
    with open(outfile, "w") as o:
        yaml.safe_dump(dict(out), o)


with file_logging(snakemake.log[0]):
    list_samples(sample_files=snakemake.input.sample_files,
                    outfile=snakemake.output.yml)

from os.path import join

from utils import file_logging
from utils.sample_list import SampleList


def make_sample_tab(sample_file, sample_file_out, sample_dir, subdir=None, qiime_style=False):
    if subdir is not None:
        sample_dir = join(sample_dir, subdir)
    sl = SampleList(sample_file)
    with open(sample_file_out, "w") as f:
        sl.write(f, absolute_paths=True, qiime_style=qiime_style,
                 out_pattern=sample_dir+"/{sample}_R{read}.fastq.gz")


with file_logging(snakemake.log[0]):
    make_sample_tab(sample_file=snakemake.input.tab,
                    sample_file_out=snakemake.output.tab,
                    sample_dir=snakemake.input.sample_dir,
                    **snakemake.params)

import yaml

from utils import file_logging
from utils.sample_list import SampleList


def make_pooling_list(sample_files, sample_file_out, info_file_out):
    # read sample files
    # TODO: indexes.tsv for demultiplexing
    sample_lists = [SampleList(f) for f in sample_files]
    # there should be only one layout present (paired or single)
    unique_layout = set(l.layout for l in sample_lists)
    assert len(unique_layout) == 1
    layout = next(iter(unique_layout))
    # get a map of sample -> read files
    sample_dicts = [dict(l.samples()) for l in sample_lists]
    unique_names = sorted(set(s for d in sample_dicts for s in d))
    # write sample lists
    sample_list_out = SampleList(layout=layout)
    sample_pools = {}
    for sample in unique_names:
        # obtain a list of all run files for that sample
        run_files = [d[sample] for d in sample_dicts if sample in d]
        # sort the runs by path to have a consistent output
        run_files.sort()
        # group by R1/R2
        run_files = list(zip(*run_files))
        out_reads = [f"{sample}_R{i+1}.fastq.gz" for i in range(sample_list_out.n_reads)]
        sample_list_out.add(sample, out_reads)
        sample_pools[sample] = {f"R{i+1}": f for i, f in enumerate(run_files)}
    # write new sample list to files
    # TODO: we only list file names, not paths. This sample list is mostly needed
    # for the pipeline to work, but the paths are unimportant since already known
    with open(sample_file_out, "w") as f:
        sample_list_out.write(f)
    # write YAML file with all pooling info
    # note: sample keys are sorted (unless sort_keys=False), but the
    # but the order of read pooling is consistent (files are in a sorted list)
    with open(info_file_out, "w") as f:
        yaml.safe_dump(sample_pools, f)


with file_logging(snakemake.log[0]):
    make_pooling_list(sample_files=snakemake.input.sample_files,
                      sample_file_out=snakemake.output.sample_file,
                      info_file_out=snakemake.output.yml)

import os
from os.path import abspath
from multiprocessing.dummy import Pool   # only threaded version necessary
from subprocess import check_call
import sys
# import gzip
# import shutil

import yaml

from utils import file_logging


def do_pooling(data):
    source_paths, target_path = data
    # make sure that all file names are the same
    # TODO: this assertion fails with 'single.rev' layout
    # assert len(set(os.path.basename(p) for p in list(source_paths) + [target_path])) == 1
    try:
        if len(source_paths[0]) == 1:
            # we only need to link
            assert len(source_paths) == 1
            os.symlink(abspath(source_paths[0]), abspath(target_path))
        else:
            check_call(["zcat {} | gzip -c > {}".format(" ".join(source_paths), target_path)], shell=True)
            # with gzip.open(target_path, 'w') as out:
            #     for path in source_paths:
            #         print("source", path, "to", target_path)
            #         with gzip.open(path) as f:
            #             shutil.copyfileobj(f, out)
            return data
    except Exception as e:
        return e


def pool_raw(pooling_info, outdir, processes=1):
    os.makedirs(outdir, exist_ok=True)
    # list of input -> output paths
    with open(pooling_info) as f:
        pool_files = yaml.safe_load(f)
    # send these lists to the workers
    p = Pool(processes)
    # Create a flat sequence of read files to combine
    # and the corresponding output file name
    args = ((f, os.path.join(outdir, f"{sample}_{read}.fastq.gz"))
            for sample, sample_files in pool_files.items()
            for read, f in sample_files.items())
    for res in p.imap_unordered(do_pooling, args):
        if isinstance(res, Exception):
            raise res
        source_paths, target_path = res
        print("{} > {}".format(" ".join(source_paths), target_path),
              file=sys.stderr)


with file_logging(snakemake.log[0]):
    pool_raw(pooling_info=snakemake.input.yml,
                outdir=snakemake.output.fq,
                processes=snakemake.threads)

import re
import sys
from subprocess import check_call

from utils import file_logging


def denoise_paired(
        input, 
        layout,
        seq_out, 
        tab_out, 
        stats_out,
        trunc_qual=None,
        trunc_len=None,
        max_err=None,
        # the following settings have defaults (can be undefined without error)
        merge_maxdiffs=0,
        chimera_method=None,
        pooling_method=None,
        threads=1,
        **unused
):
    cmd = ["qiime", "dada2"]

    # https://docs.qiime2.org/2023.7/plugins/available/dada2/denoise-paired/
    if layout == "paired":
        cmd += [
            "denoise-paired",
            "--p-trunc-len-f", "{:.0f}".format(trunc_len["fwd"]),
            "--p-trunc-len-r", "{:.0f}".format(trunc_len["rev"]),
            "--p-max-ee-f", str(max_err["fwd"]),
            "--p-max-ee-r", str(max_err["rev"]),
        ]
    else:
        dir_ = "fwd"
        if layout.endswith(".rev"):
            layout = re.sub(r"\.rev$", "", layout)
            dir_ = "rev"
        assert layout == "single"
        # https://docs.qiime2.org/2023.7/plugins/available/dada2/denoise-single
        cmd += [
            "denoise-single",
            "--p-trunc-len", "{:.0f}".format(trunc_len[dir_]),
            "--p-max-ee", str(max_err[dir_]),
        ]

    # common options

    if chimera_method is None:
        chimera_method = "consensus"
    elif chimera_method == "per-sample":
        chimera_method = "none"
    if pooling_method is None:
        pooling_method = "independent"
    elif pooling_method == "pooled":
        print('Warning: "pooling_method = pooled" not possible for QIIME2',
              file=sys.stderr)
    cmd += [
        "--i-demultiplexed-seqs", input,
        "--p-trunc-q", "{:.0f}".format(trunc_qual),
        "--p-n-reads-learn", "1000000",  # TODO: not configurable
        "--p-chimera-method", chimera_method,
        "--verbose",
        "--p-n-threads", str(threads),
        "--o-representative-sequences", seq_out,
        "--o-table", tab_out,
        "--o-denoising-stats", stats_out
    ]

    if merge_maxdiffs > 0:
        if merge_maxdiffs > 1:
            print("WARNING: 'merge_maxdiffs' is > 1, but QIIME only allows up to one difference.",
                  file=sys.stderr)
        cmd.append("--p-allow-one-off")

    print("Call: " + " ".join(cmd), file=sys.stderr)
    check_call(cmd, stdout=sys.stdout, stderr=sys.stderr)


with file_logging(snakemake.log[0]):
    denoise_paired(
        snakemake.input.trim,
        layout=snakemake.wildcards.layout,
        seq_out=snakemake.output.asvs,
        tab_out=snakemake.output.tab,
        stats_out=snakemake.output.stats,
        threads=snakemake.threads,
        **snakemake.params.par
    )