diff --git a/rules/ASSEMBLY_ANNOTATION_RULES b/rules/ASSEMBLY_ANNOTATION_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..e8c536b63acec655dae7ad97d9126a1c9aa3c64d
--- /dev/null
+++ b/rules/ASSEMBLY_ANNOTATION_RULES
@@ -0,0 +1,1079 @@
+# For running the ASSEMBLY and ANNOTATION workflow for ONT data
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+RUNS=[config["runs"]["first"],
+    config["runs"]["second"]]
+#    config["runs"]["third"]]
+BARCODES=config["barcodes"]
+SAMPLES=["ONT3_MG_xx_Rashi_S11"]
+sr_sample=["ONT3_MG_xx_Rashi_S11"]
+# SAMPLES_ALL=config["samples"]
+REFERENCES=["igc", "hg38"]
+IGC_URI=config["igc"]["uri"]
+HG38_URI=config["hg38"]["uri"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa_mem", "minimap2"]
+
+
+######
+# RULES
+######       
+# rule unpack:
+#     input: os.path.join(DATA_DIR, "2018-11-08-minion_latestdata.7z")
+#     output: os.path.join(DATA_DIR, "raw/raw_unpacked.done")
+#     threads: config["p7zip"]["threads"]
+#     shell:
+#         """
+#         date
+#         {config[p7zip][bin]} x -o$(dirname {output}) {input} -mmt{threads} &> {output}
+#         date
+#         """
+
+# Take the per-read FAST5 files and create multi-FAST5 files including a batch of reads (typically not all reads, though).
+#rule create_multifast5s:
+#    input: os.path.join(DATA_DIR, "raw/{run}")
+#    output: protected(directory(os.path.join(DATA_DIR, "multifast5/{run}")))
+#    threads: config["ont_fast5_api"]["single_to_multi_fast5"]["threads"]
+#    log: os.path.join(DATA_DIR, "multifast5/{run}/create_multifast5s")
+#    shell:
+#        """
+#        (date &&\
+#        {config[ont_fast5_api][single_to_multi_fast5][bin]} --input_path {input} --save_path {output} --filename_base multifast5 --batch_size {config[ont_fast5_api][single_to_multi_fast5][batch]} --recursive -t {threads} &&\
+#        touch {output} &&\
+#        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#        """
+
+# FLO-MIN107 SQK-LSK108           dna_r9.5_450bps
+#   --verbose_logs                    Enable verbose logs.
+#   --qscore_filtering                Enable filtering of reads into PASS/FAIL
+#                                     folders based on min qscore.
+#   --min_qscore arg                  Minimum acceptable qscore for a read to be
+#                                     filtered into the PASS folder
+#   --disable_pings                   Disable the transmission of telemetry
+#   --num_callers arg                 Number of parallel basecallers to create.
+#   --num_barcode_threads arg         Number of worker threads to use for
+#                                     barcoding.
+#   --barcode_kits arg                Space separated list of barcoding kit(s) or
+#                                     expansion kit(s) to detect against. Must be
+#                                     in double quotes.
+#   --trim_barcodes                   Trim the barcodes from the output sequences
+#                                     in the FastQ files.
+
+# Left in for legacy look-up reasons. GPU-based basecalling strongly encouraged.
+# rule guppy_cpu_basecall:
+#     input: os.path.join(DATA_DIR, "multifast5/{run}")
+#     output: os.path.join(RESULTS_DIR, "basecalled/guppy/cpu/{run}/basecalling.done")
+#     threads: config["guppy_cpu"]["cpu_threads"]
+#     shell:
+#         """
+#         date
+#         {config[guppy_cpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_cpu][config]} --cpu_threads_per_caller {threads} --disable_pings
+#         touch {output}
+#         date
+#         """
+
+    #input: os.path.join(DATA_DIR, "multifast5/{run}")
+# Basecalling the multi-FAST5 files.
+# Currently requires dummy target, but would be nice to have this be a "real" file target. However, it seems that files are filled incrementally, hence there is no check if the whole process was successful or not if checking only for file presence.
+rule guppy_gpu_basecall:
+    input: os.path.join(DATA_DIR, "multifast5/{run}")
+    output: os.path.join(RESULTS_DIR, "basecalled/{run}/basecalling.done")
+    log: os.path.join(RESULTS_DIR, "basecalled/{run}/basecalling")
+    threads: config["guppy_gpu"]["cpu_threads"]
+    shell:
+        """
+        (date &&\
+        {config[guppy_gpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_gpu][config]} --cpu_threads_per_caller {threads} -x {config[guppy_gpu][gpu_device]} --disable_pings --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --chunk_size {config[guppy_gpu][chunk_size]} --chunks_per_runner {config[guppy_gpu][chunks_per_runner]} --gpu_runners_per_device {config[guppy_gpu][runners_per_device]} --num_callers {config[guppy_gpu][num_callers]} && \
+        touch {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule guppy_gpu_basecall_NO_MOD:
+    input: os.path.join(DATA_DIR, "multifast5/{run}")
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/NO_MOD_basecalling.done")
+    log: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/NO_MOD_basecalling")
+    threads: config["guppy_gpu"]["cpu_threads"]
+    shell: 
+        """
+        (date &&\
+        {config[guppy_gpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_gpu][hac_config]} --cpu_threads_per_caller {threads} -x {config[guppy_gpu][gpu_device]} --disable_pings --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --chunk_size {config[guppy_gpu][chunk_size]} --chunks_per_runner {config[guppy_gpu][chunks_per_runner]} --gpu_runners_per_device {config[guppy_gpu][runners_per_device]} --num_callers {config[guppy_gpu][num_callers]} && \
+        touch {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Demultiplex.
+# Currently requires dummy target, but would be nice to have this be a "real" file target. However, it seems that files are filled incrementally, hence there is no check if the whole process was successful or not if checking only for file presence.
+rule guppy_gpu_demux:
+    input: rules.guppy_gpu_basecall.output
+    output: os.path.join(RESULTS_DIR, "basecalled/{run}/demux.done")
+    log: os.path.join(RESULTS_DIR, "basecalled/{run}/demux")
+    threads: config["guppy_barcoder"]["threads"]
+    shell:
+        """
+        (date &&\
+        {config[guppy_barcoder][bin]} --input_path $(dirname {input}) --save_path $(dirname {output}) -t {threads} --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --trim_barcodes && \
+        touch {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+rule guppy_gpu_demux_NO_MOD:
+    input: rules.guppy_gpu_basecall_NO_MOD.output
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/demux_NO_MOD.done")
+    log: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/demux")
+    threads: config["guppy_barcoder"]["threads"]
+    shell:
+        """
+        (date &&\
+        {config[guppy_barcoder][bin]} --input_path $(dirname {input}) --save_path $(dirname {output}) -t {threads} --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --trim_barcodes && \
+        touch {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Merge multiple basecalled files, i.e., FASTQ files per barcode and run. Each multi-FAST5 file will produce a respective FASTQ file.
+# TODO: Simplify output path to `basecalled/{run}/{barcode}.fastq.gz`. This will require to change `$(dirname {output})` to `$(dirname {output})/{wildcards.barcode}`
+rule merge_individual_fastq_per_barcode:
+    input: rules.guppy_gpu_demux.output
+    output: os.path.join(RESULTS_DIR, "basecalled/{run}/{barcode}/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat $(find $(dirname {output}) -name "*.fastq.gz" | sort) > {output}
+        touch {output}
+        date
+        """
+
+rule merge_individual_fastq_per_barcode_NO_MOD:
+    input: rules.guppy_gpu_demux_NO_MOD.output
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/{barcode}/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat $(find $(dirname {output}) -name "*.fastq.gz" | sort) > {output}
+        touch {output}
+        date
+        """
+
+# Helper function to get the input in order to create a *single* fastq file for all runs with the *same* barcode
+# TODO: Adjust according to TODO-merge_individual_fastq_per_barcode
+def get_all_fastq_for_barcode(wildcards):
+    return expand(os.path.join(RESULTS_DIR, "basecalled/{run}/{barcode}/{barcode}.fastq.gz"), run=RUNS, barcode=wildcards.barcode)
+
+def get_all_fastq_for_barcode_NO_MOD(wildcards):
+    return expand(os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/{barcode}/{barcode}.fastq.gz"), run=RUNS, barcode=wildcards.barcode)
+
+# Merge basecalled files (FASTQ) per barcode, i.e., combine reads from multiple runs according to their barcode.
+rule get_single_fastq_per_barcode:
+    input: get_all_fastq_for_barcode 
+    output: os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat {input} > {output}
+        date
+        """
+rule get_single_fastq_per_barcode_NO_MOD:
+    input: get_all_fastq_for_barcode_NO_MOD
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/merged/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat {input} > {output}
+        date
+        """
+
+# Get basic statistics of the ONT data; over all runs
+# TODO: Simplify top of path: replace `lr` by `nanostat`
+# TODO: Write more generic: have patterned `input: os.path.join(RESULTS_DIR, "basecalled/{prefix}/{barcode}.fastq.gz")` and patterned `output: os.path.join(RESULTS_DIR, "qc/nanostat/{prefix}/{barcode}/{barcode}NanoStats.txt")
+# 
+rule nanostat:
+    input: rules.get_single_fastq_per_barcode.output
+    output: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats.txt")
+    log: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        (date &&\
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output}) &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule nanostat_NO_MOD:
+    input: rules.get_single_fastq_per_barcode_NO_MOD.output
+    output: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats_NO_MOD.txt")
+    log: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats_NO_MOD")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        (date &&\
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output}) &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Deprecated if TODO-nanostat is completed
+# Get basic statistics of the ONT data; per run
+rule nanostat_per_run:
+    input: os.path.join(RESULTS_DIR, "basecalled/{run}/{barcode}/{barcode}.fastq.gz")
+    output: os.path.join(RESULTS_DIR, "qc/lr/{run}/{barcode}/{barcode}NanoStats.txt")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        date
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output})
+        date
+        """
+
+rule nanostat_per_run_NO_MOD:
+    input: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/{barcode}/{barcode}.fastq.gz")
+    output: os.path.join(RESULTS_DIR, "qc/lr/{run}/{barcode}/{barcode}NanoStats_NO_MOD.txt")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        date
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output})
+        date
+        """
+
+# Download the IGC (http://www.nature.com/nbt/journal/vaop/ncurrent/full/nbt.2942.html)
+# Downloading using wget-command instead of FTP-remote in snakemake as the former offers progress information
+rule fetch_igc:
+    output: protected(os.path.join(DB_DIR, "igc/igc.fna.gz"))
+    log: os.path.join(DB_DIR, "igc/fetch_igc")
+    shell:
+        """
+        (date &&\
+        wget ftp://{IGC_URI} -O {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Download HG38 human reference genome
+# Downloading using wget-command instead of FTP-remote in snakemake as the former offers progress information
+rule fetch_hg38:
+    output: protected(os.path.join(DB_DIR, "hg38/GCF_000001405.38_GRCh38.p12_genomic.fna.gz"))
+    log: os.path.join(DB_DIR, "hg38/fetch_hg38")
+    shell:
+        """
+        (date &&\
+        wget ftp://{HG38_URI} -O {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule create_link_to_hg38:
+    input: os.path.join(DB_DIR, "hg38/GCF_000001405.38_GRCh38.p12_genomic.fna")
+    output: protected(os.path.join(DB_DIR, "hg38/hg38.fna"))
+    shell:
+        """
+        date
+        ln -s $(basename -s "" {input}) {output}
+        date
+        """
+
+rule unpack_reference_data:
+    input: os.path.join(DB_DIR, "{prefix}.fna.gz")
+    output: protected(os.path.join(DB_DIR, "{prefix}.fna"))
+    shell:
+        """
+        date
+        gzip -dc {input} > {output}
+        date
+        """
+
+# Build a minimap2 index to map *ONT* reads against
+rule build_minimap2_ont_index:
+    input: "{prefix}.fna"
+    output: "{prefix}-ont.mmi"
+    log: "{prefix}.minimap2_ont_index"
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -x map-ont -d {output} {input} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Build a minimap2 short read index in order to map short reads to a reference/long read-based contigs
+# NOTE: Unclear whether this would be also the right option for mapping short reads to long reads (i.e., reads-to-reads), due to the increased error rate in long reads. For more discussions on this, see https://github.com/lh3/minimap2/issues/407#issuecomment-493823000
+# TODO: Combine this rule and `build_minimap2_ont_index` using a conditional on the suffix (ont vs. sr). This makes the code more readable as it removes largely redundant `input` and `output` keywords
+rule build_minimap2_sr_index:
+    input: "{prefix}.fna"
+    output: "{prefix}-sr.mmi"
+    log: "{prefix}.minimap2_sr_index"
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -x sr -d {output} {input} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Build a bwa index (may be used to map short reads against it, or long reads by using the "-x ont2d" flag)
+rule build_bwa_index:
+    input: "{prefix}.fna"
+    output: index_amb="{prefix}.amb",
+       index_ann="{prefix}.ann",
+       index_bwt="{prefix}.bwt",
+       index_pac="{prefix}.pac",
+       index_sa="{prefix}.sa"
+    log: "{prefix}.bwa_index"
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa index {input} -p {wildcards.prefix} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Simple utility function
+#rule convert_fastq_to_fasta:
+#    input: "{prefix}.fastq.gz"
+#    output: "{prefix}.fna"
+#    shell:
+#        """
+#        date
+#        zcat {input} | sed -n '1~4s/^@/>/p;2~4p' > {output}
+#        date
+#        """
+#
+#rule convert_fq_to_fasta:
+#    input: "{prefix}.fq.gz"
+#    output: "{prefix}.fna"
+#    shell:
+#        """
+#        date
+#        zcat {input} | sed -n '1~4s/^@/>/p;2~4p' > {output}
+#        date
+#        """
+
+# Map long reads to a reference (IGC/HG38/other). The resulting BAM file can be use for many things, among other to estimate which long reads are "known" or "novel", or to estimate the coverage of reference sequences (regions thereof).
+rule map_long_reads_to_reference_w_minimap2:
+    input: 
+        query=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        index=os.path.join(DB_DIR, "{reference}/{reference}-ont.mmi"),
+        index_fa=os.path.join(DB_DIR, "{reference}/{reference}.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}-x-{reference}.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}-x-{reference}")
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}-x-{reference}.minimap2_samtools")
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query} | samtools view -bT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Preprocess the short reads using fastp
+rule run_fastp_on_short_reads:
+    input:
+        r1=os.path.join(config["short_reads_prefix"], "ONT3_MG_xx_Rashi_S11_R1_001.fastq.gz"),
+        r2=os.path.join(config["short_reads_prefix"], "ONT3_MG_xx_Rashi_S11_R2_001.fastq.gz")
+    output:
+        r1=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11_R1_001.fastp.fq.gz"),
+        r2=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11_R2_001.fastp.fq.gz"),
+        html=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11.fastp.html"),
+        json=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11.fastp.json")
+    log: os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11.fastp")
+    conda: "../envs/fastp.yaml"
+    shell:
+        """
+        (date &&\
+        fastp -l {config[fastp][min_length]} -i {input.r1} -I {input.r2} -o {output.r1} -O {output.r2} -h {output.html} -j {output.json} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Simplify paths by removing the `sr`.
+rule run_fastqc_on_preprocessed_short_reads:
+    input: os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_{suffix}.fq.gz"),
+    output:
+        html=os.path.join(RESULTS_DIR, "qc/sr/fastqc/{sr_sample}/{sr_sample}_{suffix}.fastqc.html"),
+        zip=os.path.join(RESULTS_DIR, "qc/sr/fastqc/{sr_sample}/{sr_sample}_{suffix}.fastqc.zip") # the suffix _fastqc.zip is necessary for multiqc to find the file. If not using multiqc, you are free to choose an arbitrary filename
+    conda: "../envs/fastqc.yaml"
+    script: "../scripts/run_fastqc.py" # Use a separate script rather than a simple call to FASTQC to overcome possible race conditions when running in parallel (https://bitbucket.org/snakemake/snakemake-wrappers/raw/7f5e01706728e32d52aa413f213e950f50bf4129/bio/fastqc/wrapper.py)
+
+rule map_short_reads_to_reference_w_minimap2:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index=os.path.join(DB_DIR, "{reference}/{reference}-sr.mmi"),
+        index_fa=os.path.join(DB_DIR, "{reference}/{reference}.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}-x-{reference}.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}-x-{reference}")
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}-x-{reference}.minimap2_samtools")
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -bT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_short_reads_to_reference_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(DB_DIR, "{reference}/{reference}.amb"),
+        index_ann=os.path.join(DB_DIR, "{reference}/{reference}.ann"),
+        index_bwt=os.path.join(DB_DIR, "{reference}/{reference}.bwt"),
+        index_pac=os.path.join(DB_DIR, "{reference}/{reference}.pac"),
+        index_sa=os.path.join(DB_DIR, "{reference}/{reference}.sa"),
+        index_fa=os.path.join(DB_DIR, "{reference}/{reference}.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}-x-{reference}.bam")
+    params:
+        index_prefix=os.path.join(DB_DIR, "{reference}/{reference}"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}-x-{reference}")
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}-x-{reference}.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -bT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Compute the genome coverage histogram - https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html
+# The resulting histogram is used to compute the expected, i.e., average, coverage of the underyling reference (which may be de novo contigs).
+rule run_genomecov:
+    input: os.path.join(RESULTS_DIR, "mapping/{prefix}.bam")
+    output: os.path.join(RESULTS_DIR, "genomecov/{prefix}.cov.txt")
+    conda: "../envs/bedtools.yaml"
+    shell:
+        """
+        date
+        bedtools genomecov -ibam {input} > {output}
+        date
+        """
+
+# Compute the expected, i.e., average, coverage given the coverage histogram
+rule compute_avg_genomecov:
+    input: os.path.join(RESULTS_DIR, "genomecov/{prefix}.cov.txt")
+    output: os.path.join(RESULTS_DIR, "genomecov/{prefix}.avg_cov.txt")
+    shell:
+        """
+        date
+        cat {input} | awk -f {config[compute_avg_coverage][bin]} | tail -n+2 > {output}
+        date
+        """
+
+# TODO: All above rules that build an index build it in the same path as the reference/target file
+# In fact, it is redundant with `build_bwa_index`. The only advantage is that this replication allows more fine-grained control of the resources (slurm, etc.). This is, however, probably not needed as `build_bwa_index` can be quite resource-intensive, e.g., for the IGC.
+# Similar to building a minimap2 index, but this time for bwa
+rule build_long_read_bwa_index:
+    input: rules.get_single_fastq_per_barcode.output
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.amb"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.ann"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.bwt"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.pac"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.sa")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}")
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.bwa_index")
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa index {input} -p {params.prefix} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+# TODO: If TODO-build_long_read_bwa_index is completed, need to adjust the `input.index` here too
+# config[bwa][long_reads_index][opts] according to https://github.com/sfu-compbio/colormap/blob/baa6802b43da8640c305032b7f7fcadd7c496a03/run_corr.sh#L45 and https://academic.oup.com/bioinformatics/article/32/17/i545/2450793#112481513 -> "2.2.3 Mapping parameters"
+rule map_short_reads_to_long_reads:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.amb"),
+        index_fa=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fna"),
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{sr_sample}-x-{barcode}.bam")
+    threads: config["bwa"]["threads"]
+    wildcard_constraints:
+        barcode='^\-' 
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{sr_sample}-x-{barcode}.")
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{sr_sample}-x-{barcode}.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa mem {config[bwa][long_reads_index][opts]} -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -S -b -u -T {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Add extra rules for other assemblers instead of a generic rule as this offers then greater flexibility, e.g., in terms of resource-reservation.
+rule assemble_long_reads_w_flye:
+    input: rules.get_single_fastq_per_barcode.output
+    output: os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna")
+    threads: config["flye"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.flye")
+#    conda: "../envs/flye.yaml"
+    conda: "../envs/flye_v2_7.yaml"
+    shell:
+        """
+        (date &&\
+        flye --nano-raw {input} --meta --out-dir $(dirname {output}) --genome-size {config[flye][genome_size]} --threads {threads} &&\
+        ln -s assembly.fasta {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# "r941_min_high" is the MinION model, high accuarcy
+rule polish_long_read_contigs_w_medaka:
+    input: 
+        basecalls=rules.get_single_fastq_per_barcode.output,
+        contigs=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/medaka/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    threads: config["medaka"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/medaka/round_{round}/lr/merged/{barcode}/assembly_polished.medaka")
+    conda: "../envs/medaka.yaml"
+    shell:
+        """
+        (date &&\
+        medaka_consensus -i {input.basecalls} -d {input.contigs} -o $(dirname {output}) -t {threads} -m r941_min_high &&\
+        ln -s consensus.fasta {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Having the "initial" explicitly in the rule name is not ideal. However, the path to the original contigs is somewhat variable and hence a discriminative rule name was chosen.
+# As piping is no longer supported by racon since v. 1.0.0 (https://github.com/isovic/racon/issues/46#issuecomment-519380742) and the file extensions are checked for all files including the contigs), some temporary files have to be created and removed afterwards.
+rule polish_long_read_contigs_w_racon_initial_round:
+    input: 
+        basecalls=rules.get_single_fastq_per_barcode.output,
+        contigs=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna"),
+        lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    params:
+        tmp_sam=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.tmp.sam"),
+        tmp_contigs=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/assembly_polished.tmp_draft.fa")
+    threads: config["racon"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/assembly_polished.racon.log")
+    conda: "../envs/racon.yaml"
+    shell:
+        """
+        (date &&\
+        samtools view -h {input.lr_to_lr_contigs_mapping} > {params.tmp_sam} &&\
+        cp {input.contigs} {params.tmp_contigs} &&\
+        racon -t {threads} {input.basecalls} {params.tmp_sam} {params.tmp_contigs} > {output} &&\
+        rm {params.tmp_sam} {params.tmp_contigs} &&\
+        date) 2> >(tee {log}) > >(tee {log})
+        """
+
+#polished_contigs=os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/contigs.polished.fna")
+#ln -s contigs.polished.fasta {output.polished_contigs} &&\
+#rule assemble_long_and_short_reads_w_operams:
+#    input: 
+#        lr=rules.get_single_fastq_per_barcode.output,
+#        sr_r1=rules.run_fastp_on_short_reads.output.r1,
+#        sr_r2=rules.run_fastp_on_short_reads.output.r2,
+#        sr_contigs=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+#    output: 
+#        contigs=os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+#    params:
+#        tmp_lr=os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/tmp_lr_{barcode}.fna")
+#    threads: config["operams"]["threads"]
+#    log: os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/opera_ms.log")
+#    shell:
+#        """
+#        (date &&\
+#        zcat {input.lr} > {params.tmp_lr} &&\
+#        {config[operams][bin]} --long-read-mapper minimap2 --num-processors {threads} --short-read1 {input.sr_r1} --short-read2 {input.sr_r2} --long-read {params.tmp_lr} --contig-file {input.sr_contigs} --out-dir $(dirname {output.contigs}) &&\
+#        ln -s contigs.fasta {output.contigs} &&\
+#        rm {params.tmp_lr} &&\
+#        date) 2> >(tee {log}) > >(tee {log})
+#        """
+
+#zcat {input.lr} > {params.tmp_lr} &&\
+rule assemble_long_and_short_reads_w_metaspades:
+    input: 
+        lr=rules.get_single_fastq_per_barcode.output,
+        sr_r1=rules.run_fastp_on_short_reads.output.r1,
+        sr_r2=rules.run_fastp_on_short_reads.output.r2,
+        sr_contigs=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+    output: 
+        contigs=os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    params:
+        tmp_lr=os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_{barcode}-sr_{sr_sample}/tmp_lr_{barcode}.fna")
+    threads: config["metaspades"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_{barcode}-sr_{sr_sample}/metaspades.log")
+    conda: "../envs/spades.yaml"
+    shell:
+        """
+        (date &&\
+        spades.py --meta -o $(dirname {output.contigs}) -k 21,33,55,77 -t {threads} -1 {input.sr_r1} -2 {input.sr_r2} --nanopore {input.lr} &&\
+        ln -s contigs.fasta {output.contigs} &&\
+        date) &> >(tee {log})
+        """
+
+rule map_short_reads_to_hybrid_contigs_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs")
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_long_reads_to_hybrid_contigs_w_bwa_mem:
+    input: 
+        query_lr=rules.get_single_fastq_per_barcode.output,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -x ont2d -t {threads} {params.index_prefix} {input.query_lr} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# rule build_minimap2_long_read_contigs_ont_index:
+#     input: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}.fna")
+#     output: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}-ont.mmi")
+#     conda: "../envs/minimap2.yaml"
+#     shell:
+#         """
+#         date
+#         minimap2 -x map-ont -d {output} {input}
+#         date
+#         """
+
+rule map_long_reads_to_long_read_contigs:
+    input: 
+        query=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        index=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly-ont.mmi"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam")
+#    wildcard_constraints:
+#        barcode='barcode\d+'
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x_{barcode}_contigs")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.minimap2_samtools")
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_long_reads_to_long_read_polished_contigs:
+    input: 
+        query=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        index=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/{polisher}/round_{round}/lr/merged/{barcode}/assembly_polished-ont.mmi"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/{polisher}/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/polished/{polisher}/round_{round}/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_polished_contigs.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/polished/{polisher}/round_{round}/lr/merged/{barcode}/{barcode}_reads-x_{barcode}_polished_contigs")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/polished/{polisher}/round_{round}/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_polished_contigs.minimap2_samtools")
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# rule build_minimap2_long_read_contigs_sr_index:
+#     input: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}.fna")
+#     output: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}-sr.mmi")
+#     conda: "../envs/minimap2.yaml"
+#     shell:
+#         """
+#         date
+#         minimap2 -x sr -d {output} {input}
+#         date
+#         """
+
+rule map_short_reads_to_long_read_contigs_w_minimap2:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly-sr.mmi"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}_reads-x_{barcode}_contigs")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.minimap2_samtools")
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output}
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Also map the short reads with bwa mem as it might be more sensitive (https://github.com/lh3/minimap2/issues/214#issuecomment-413492006)
+rule map_short_reads_to_long_read_contigs_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x_{barcode}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule assemble_short_reads_w_megahit:
+    input: 
+        sr_r1=rules.run_fastp_on_short_reads.output.r1,
+        sr_r2=rules.run_fastp_on_short_reads.output.r2
+    output: os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+    threads: config["megahit"]["threads"]
+    conda: "../envs/megahit.yaml"
+    params:
+        tmp_dir=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}_tmp")
+    log: os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.megahit")
+    shell:
+        """
+        (date &&\
+        megahit -1 {input.sr_r1} -2 {input.sr_r2} -t {threads} -o {params.tmp_dir} &&\
+        mv {params.tmp_dir}/* $(dirname {output}) &&\
+        rmdir {params.tmp_dir} &&\
+        ln -s final.contigs.fa {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# N.B. A similar function exists above and code duplication is generally discouraged.
+rule map_short_reads_to_short_read_contigs_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.sa"),
+        index_fa=rules.assemble_short_reads_w_megahit.output
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_long_reads_to_short_read_contigs_w_bwa_mem:
+    input: 
+        query_lr=rules.get_single_fastq_per_barcode.output,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.sa"),
+        index_fa=rules.assemble_short_reads_w_megahit.output
+    output: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -x ont2d -t {threads} {params.index_prefix} {input.query_lr} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+#npa, npc, .npl, and .npo
+rule run_nonpareil_on_short_reads:
+    input: os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_R{orientation}_{suffix}.fastp.fq.gz"),
+    output: 
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npa"),
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npc"),
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npl"),
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npo")
+    params:
+        out_prefix=os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp"),
+        tmp_fasta=os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.tmp.fna")
+    threads: config["nonpareil"]["threads"]
+    log: os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.nonpareil")
+    conda: "../envs/nonpareil.yaml"
+    shell:
+        """
+        (date &&\
+        zcat {input} | sed -n '1~4s/^@/>/p;2~4p' > {params.tmp_fasta} &&\
+        nonpareil -s {params.tmp_fasta} -R {config[nonpareil][memory]} -t {threads} -T kmer -f fasta -b {params.out_prefix} &&\
+        rm {params.tmp_fasta} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Cleanup how the input files are defined. Currently explicitly specified, which is discouraged to do.
+rule create_nanopolish_index:
+    input:
+        multifast5_dir_01=os.path.join(DATA_DIR, "multifast5/20181108_0827_test"),
+        multifast5_dir_02=os.path.join(DATA_DIR, "multifast5/20181107_0906_same"),
+        multifast5_dir_03=os.path.join(DATA_DIR, "multifast5/20181106_1450_noselection_sizeselection"),
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz")
+    output:
+        index=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index"),
+        fai=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.fai"),
+        index_gzi=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.gzi"),
+        readdb=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.readdb")
+    log: os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.nanopolish_index")
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        (date &&\
+        nanopolish index -d {input.multifast5_dir_01}/ -d {input.multifast5_dir_02}/ -d {input.multifast5_dir_03}/ {input.reads} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule samtools_index:
+    input: "{prefix}.bam"
+    output: "{prefix}.bam.bai"
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        date
+        samtools index {input} {output}
+        date
+        """
+
+# # Legacy version before snakemake's `checkpoint` was used
+# rule polish_long_read_contigs_w_nanopolish:
+#     input:
+#         lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+#         lr_to_lr_contigs_mapping_index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam.bai"),
+#         reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+#         draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+#         nanopolish_index=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index"),
+#         nanopolish_fai=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.fai"),
+#         nanopolish_index_gzi=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.gzi"),
+#         nanopolish_readdb=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.readdb")
+#     output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+#     log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.nanopolish")
+#     threads: 4
+#     conda: "../envs/nanopolish.yaml"
+#     shell:
+#         """
+#         (date &&\
+#         nanopolish_makerange.py {input.draft_contigs} | parallel --results $(dirname {output}) -P 7 \
+#         nanopolish variants --consensus -o $(dirname {output})/polished.{{1}}.vcf -w {{1}} -r {input.reads} -b {input.lr_to_lr_contigs_mapping} -g {input.draft_contigs} -t {threads} --min-candidate-frequency 0.1 &&\
+#         nanopolish vcf2fasta -g {input.draft_contigs} $(dirname {output})/polished.*.vcf > {output} &&\
+# 
+#         date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#         """
+
+# Define the individual windows that nanopolish will parallelize over
+checkpoint nanopolish_makerange:
+    input: os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna")
+    output: directory(os.path.join(RESULTS_DIR, "assembly/flye/polished/nanopolish/round_{round}/lr/merged/{barcode}/windows"))
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        mkdir {output}
+        nanopolish_makerange.py {input} > {output}/ranges.txt
+        for i in $(head -n 20 {output}/ranges.txt); do echo $i > {output}/"$i".coord; done
+        date
+        """
+
+# Perform variant calling on single window
+rule nanopolish_variant_call:
+    input: 
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+        lr_to_lr_contigs_mapping_index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam.bai"),
+        draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+        window=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/windows/{window}.coord")
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/vcfs/{window}.vcf")
+    threads: 4
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        nanopolish variants --consensus -o {output} -w {wildcards.window} -r {input.reads} -b {input.lr_to_lr_contigs_mapping} -g {input.draft_contigs} -t {threads} --min-candidate-frequency 0.1
+        """
+
+# Gather all the individual VCF files that were produced.
+# The "overview" comes from the rule that created the "split-up", i.e., from the `nanopolish_makerange` rule.
+# But the actual files that are aggregated here are those that should eventually be created and used, i.e., the VCF files
+# N.B. snakemake threw an error when the `window=` line below was wrong, but the error was supposedly coming from `nanopolish_makerange`. So much about tracking errors, which is not always straightforward as one might wish.
+def aggregate_nanopolish_vcf(wildcards):
+    checkpoint_output = checkpoints.nanopolish_makerange.get(**wildcards).output[0]
+    return expand(os.path.join(RESULTS_DIR, "assembly/flye/polished/nanopolish/round_{round}/lr/merged/{barcode}/vcfs/{window}.vcf"),
+            round=wildcards.round,
+            barcode=wildcards.barcode,
+            window=glob_wildcards(os.path.join(checkpoint_output, "{window}.coord")).window)
+
+# Combine the information in the individual VCF files to produce a polished assembly
+rule polish_long_read_contigs_w_nanopolish:
+    input: 
+        vcf_calls=aggregate_nanopolish_vcf,
+        draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.nanopolish")
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        (date &&\
+        nanopolish vcf2fasta -g {input.draft_contigs} {input.vcf_calls} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule nanopolish_methylation_call:
+    input: 
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+        lr_to_lr_contigs_mapping_index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam.bai"),
+        draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+        window=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/windows/{window}.coord")
+    output: os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_calls/{window}.methylation_calls.tsv")
+    threads: 4
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        nanopolish call-methylation -t {threads} -r {input.reads} -b {input.lr_to_lr_contigs_mapping} -g {input.draft_contigs} -w {wildcards.window} > {output}
+        """
+
+def aggregate_nanopolish_methylation(wildcards):
+    checkpoint_output = checkpoints.nanopolish_makerange.get(**wildcards).output[0]
+    return expand(os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_calls/{window}.methylation_calls.tsv"),
+            assembler=wildcards.assembler,
+            barcode=wildcards.barcode,
+            round=wildcards.round,
+            window=glob_wildcards(os.path.join(checkpoint_output, "{window}.coord")).window)
+
+# Combine the information in the individual methylation call files to get the overall methylation state
+rule merge_methylation_calls_from_long_read_contigs_w_nanopolish:
+    input: aggregate_nanopolish_methylation,
+    output: os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_frequencies.tsv")
+    log: os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_frequencies")
+    shell:
+        """
+        (date &&\
+        ../scripts/calculate_methylation_frequency.py {input} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Works until the final shred and polish step, then crashes with a Key error
+# rule run_rebaler_on_short_read_contigs:
+#     input:
+#         lr=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+#         sr_contigs=rules.assemble_short_reads_w_megahit.output
+#     output: os.path.join(RESULTS_DIR, "assembly/rebaler/lr/merged/{barcode}/lr_{barcode}-sr_contigs_{sr_sample}.fna")
+#     threads: config["rebaler"]["threads"]
+#     params:
+#         tmp_sr_contigs=os.path.join(RESULTS_DIR, "assembly/rebaler/lr/merged/{barcode}/tmp_sr_contigs.fna")
+#     conda: "../envs/rebaler.yaml"
+#     shell:
+#         """
+#         date
+#         awk 'BEGIN{{counter=0}} FNR==NR{{if($1 ~ /^>/) {{a[$1]=counter; counter=counter+1}};next}}{{if($1 ~/^>/) print ">"a[$1]; else print}}' {input.sr_contigs} {input.sr_contigs} > {params.tmp_sr_contigs} # Rename the contigs from 0 to num. contigs to ensure rebaler has no "Key error"
+#         rebaler -t {threads} {params.tmp_sr_contigs} {input.lr} > {output}
+#         rm {params.tmp_sr_contigs}
+#         date
+#         """
+
+# For info on "--longs-reads" option see https://github.com/bbuchfink/diamond/releases/tag/v0.9.23
+rule run_diamond_on_long_reads:
+    input:
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        db=config["diamond"]["db"]
+    output: os.path.join(RESULTS_DIR, "annotation/diamond/lr/merged/{barcode}.diamond.daa")
+    threads: config["diamond"]["threads"]
+    conda: "../envs/diamond.yaml"
+    shell:
+        """
+        date
+        diamond blastx -q {input.reads} --db {input.db} -p {threads} --long-reads --sensitive --outfmt 100 --out {output}
+        date
+        """
+
+rule call_genes_w_prodigal:
+    input: os.path.join(RESULTS_DIR, "assembly/{prefix}.fna")
+    output: os.path.join(RESULTS_DIR, "annotation/proteins/{prefix}.faa")
+    log: os.path.join(RESULTS_DIR, "annotation/proteins/{prefix}.prodigal.log")
+    conda: "../envs/prodigal.yaml"
+    shell: 
+        """
+        (date &&\
+        prodigal -a {output} -p meta -i {input} &&\
+        date) &> >(tee {log})
+        """
+
+# Write to DAA format to be able to use `diamond view` later to define custom formats.
+# N.B. `diamond view` seems not to support compressed input AND seems to require the *prefix* of the DAA file.
+# TODO: This could maybe be parametrized by the respectively used DB
+rule run_diamond_on_proteins:
+    input: 
+        proteins=rules.call_genes_w_prodigal.output,
+        db=config["diamond"]["db"]
+    output: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.daa")
+    params:
+        outfmt="100"
+    log: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.diamond.log")
+    conda: "../envs/diamond.yaml"
+    threads: config["diamond"]["threads"]
+    shell: 
+        """
+        (date &&\
+        diamond blastp -q {input.proteins} --db {input.db} -p {threads} --outfmt {params.outfmt} --out {output} &&\
+        date) &> >(tee {log})
+        """
+
+# Not sure why, but with the conda version of diamond (and `diamond view`) it was necessary to specify the entire column layout (outfmt) explicitly
+rule format_diamond_output_for_ideel:
+    input: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.daa")
+    output: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.tsv")
+    params:
+        outfmt="6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen"
+    log: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.diamond_view.log")
+    conda: "../envs/diamond.yaml"
+    shell:
+        """
+        (date &&\
+        diamond view --daa $(echo {input} | sed 's/\.daa$//') --max-target-seqs 1 --outfmt {params.outfmt} --out {output} &&\
+        date) &> >(tee {log})
+        """
diff --git a/rules/BINNING_RULES b/rules/BINNING_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..a8ce0377b16e1b1466ceb1019c1b7cad6ba48fea
--- /dev/null
+++ b/rules/BINNING_RULES
@@ -0,0 +1,167 @@
+# Rules for BINNING workflow, i.e. generating MAGs from assemblies
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mm"]
+#SAMPLES=config["samples"]
+SAMPLES=["flye", "megahit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+###################
+# Preparing files #
+###################
+rule prepare_assembly_files:
+    input:
+         ass1=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/no_barcode/assembly.fasta"),
+         ass2=os.path.join(RESULTS_DIR, "assembly/megahit/ONT3_MG_xx_Rashi_S11/final.contigs.fa"),
+         ass3=os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_no_barcode-sr_ONT3_MG_xx_Rashi_S11/contigs.fasta"),
+         bam1=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/ONT3_MG_xx_Rashi_S11_reads-x-ONT3_MG_xx_Rashi_S11-megahit_contigs.bam"),
+         bam2=os.path.join(RESULTS_DIR, "mapping/lr/merged/no_barcode/no_barcode_reads-x-no_barcode-flye_contigs.bam")
+    output:
+         fa1=os.path.join(RESULTS_DIR, "assembly/flye.fa"),
+         fa2=os.path.join(RESULTS_DIR, "assembly/megahit.fa"),
+         fa3=os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid.fa"),
+#         fa4=expand(os.path.join(RESULTS_DIR, "assembly/{mapper}_{reads}_{hybrid_assembler}.fa"), mapper=["bwa", "mmi"], reads=["sr", "lr", "merged"], hybrid_assembler=HYBRID_ASSEMBLER),
+         fa4=os.path.join(RESULTS_DIR, "assembly/bwa_sr_metaspades_hybrid.fa"),
+         fa5=os.path.join(RESULTS_DIR, "assembly/bwa_lr_metaspades_hybrid.fa"),
+         fa6=os.path.join(RESULTS_DIR, "assembly/bwa_merged_metaspades_hybrid.fa"),
+         fa7=os.path.join(RESULTS_DIR, "assembly/mmi_sr_metaspades_hybrid.fa"),
+         fa8=os.path.join(RESULTS_DIR, "assembly/mmi_lr_metaspades_hybrid.fa"),
+         fa9=os.path.join(RESULTS_DIR, "assembly/mmi_merged_metaspades_hybrid.fa"),
+         bout1=os.path.join(RESULTS_DIR, "mapping/megahit/megahit.bam"),
+         bout2=os.path.join(RESULTS_DIR, "mapping/flye/flye.bam")
+    shell:
+         """
+         (date &&\
+         ln {input.ass1} {output.fa1}
+         ln {input.ass2} {output.fa2}
+         ln {input.ass3} {output.fa3}
+         ln {input.ass3} {output.fa4}
+         ln {input.ass3} {output.fa5}
+         ln {input.ass3} {output.fa6}
+         ln {input.ass3} {output.fa7}
+         ln {input.ass3} {output.fa8}
+         ln {input.ass3} {output.fa9}
+         ln {input.bam1} {output.bout1}
+         ln {input.bam2} {output.bout2}
+         date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+         """
+
+
+#####################################
+########### MAXBIN2 #################
+#####################################
+rule maxbin2:
+    input:
+          bbref=os.path.join(RESULTS_DIR, "assembly/{sample}.fa"),
+          bbam=os.path.join(RESULTS_DIR, "mapping/{sample}/{sample}.bam")
+    threads:config["threads"]
+    conda:"../envs/binning.yaml"
+    output:
+           map=os.path.join(RESULTS_DIR, "Binning/{sample}/{sample}.sam"),
+           coverage=os.path.join(RESULTS_DIR, "Binning/{sample}/{sample}_cov.txt"),
+           abund=os.path.join(RESULTS_DIR, "Binning/{sample}/{sample}_abundance.txt"),
+           file=os.path.join(RESULTS_DIR, "Binning/{sample}/maxbin_output/maxbin_output.001.fasta")
+#        mx=directory(os.path.join(RESULTS_DIR, "Binning/{sample}/maxbin_output/"))
+    shell:
+           """
+            (date &&\
+            samtools view -h -o {output.map} {input.bbam} &&\
+            pileup.sh in={output.map} out={output.coverage} &&\
+            awk '{{print $1"\\t"$2}}' {output.coverage} | grep -v '^#' > {output.abund} &&\
+            run_MaxBin.pl -thread {threads} -contig {input.bbref} -out results/Binning/{wildcards.sample}/maxbin_output/maxbin_output -abund {output.abund} &&\
+            touch maxbin.done &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+#####################################
+########### METABAT2 ################
+#####################################
+
+rule depth_files:
+    input:
+          dep1=os.path.join(RESULTS_DIR, "mapping/{sample}/{sample}.bam")
+    threads:config["threads"]
+    conda:"../envs/binning.yaml"
+    output:
+          depout=os.path.join(RESULTS_DIR, "Binning/{sample}/{sample}_depth.txt")
+    shell:
+           """
+            (date &&\
+            jgi_summarize_bam_contig_depths --outputDepth {output.depout} {input.dep1} &&\
+            touch depth_file.done &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+rule metabat2:
+    input:
+           met1=os.path.join(RESULTS_DIR, "assembly/{sample}.fa"),
+           met2=os.path.join(RESULTS_DIR, "Binning/{sample}/{sample}_depth.txt")
+    threads:config["threads"]
+    conda:"../envs/binning.yaml"
+    output:
+           outdir=directory(os.path.join(RESULTS_DIR, "Binning/{sample}/metabat_output/"))
+    shell:
+           """
+            (date &&\
+            metabat2 -i {input.met1} -a {input.met2} -o results/Binning/{wildcards.sample}/metabat_output/{wildcards.sample}.metabat -t {threads} -m 1500 -v --unbinned &&\
+            touch metabat2.done &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+#####################################
+########### DAS_Tool ################
+#####################################
+rule scaffold_list:
+    input:
+        sca1=os.path.join(RESULTS_DIR, "Binning/{sample}/metabat_output/"),
+#       sca2=os.path.join(RESULTS_DIR, "Binning/{sample}/maxbin_output/")
+        file=os.path.join(RESULTS_DIR, "Binning/{sample}/maxbin_output/maxbin_output.001.fasta")
+    threads:config["threads"]
+#    conda:"/home/users/sbusi/apps/environments/base.yml"
+    conda: "../envs/dastool.yaml"
+    output:
+            scout1=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_metabat.scaffolds2bin.tsv"),
+            scout2=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_maxbin.scaffolds2bin.tsv")
+    shell:
+            """
+            (date &&\
+            export Rscript={config[Rscript]} &&\
+            Fasta_to_Scaffolds2Bin.sh -i {input.sca1} -e fa > {output.scout1} &&\
+            Fasta_to_Scaffolds2Bin.sh -i results/Binning/{wildcards.sample}/maxbin_output -e fasta > {output.scout2} &&\
+            touch scaffold_list.done &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+rule DAS_Tool:
+    input:
+            da1=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_metabat.scaffolds2bin.tsv"),
+            da2=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_maxbin.scaffolds2bin.tsv"),
+            da3=os.path.join(RESULTS_DIR, "assembly/{sample}.fa"),
+            db=config["dastool_database"]
+    threads:config["threads"]
+#    conda:"/home/users/sbusi/apps/environments/base.yml"
+    conda: "../envs/dastool.yaml"
+    params:
+            basename=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}")
+    output:
+            daout=directory(os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_DASTool_bins")),
+            dafile=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_proteins.faa"),
+            damfile=touch(os.path.join(RESULTS_DIR, "Binning/{sample}_das_tool.done"))
+    shell:
+            """
+            (date &&\
+            export Rscript={config[Rscript]} &&\
+            DAS_Tool -i {input.da1},{input.da2} -c {input.da3} -o {params.basename} --search_engine diamond -l maxbin2,metabat2 --write_bins 1 --write_bin_evals 1 --threads {threads} --db_directory {input.db} --create_plots 1 &&\
+            touch {output.damfile} &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
diff --git a/rules/MAPPING_RULES b/rules/MAPPING_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..0827a17e3002a22c9942fee6e9baeee8892a7282
--- /dev/null
+++ b/rules/MAPPING_RULES
@@ -0,0 +1,154 @@
+# Rule for running the mapping "workflow" for ONT analsyes
+
+#shell.executable("/bin/bash")
+#shell.prefix("source ~/.bashrc; ")
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mmi"]
+#SAMPLES=config["samples"]
+SAMPLES=["flye", "megahit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+#######################
+## Mapping to hybrid ##
+#######################
+rule build_hybrid_bwa_index:
+    input: expand(os.path.join(RESULTS_DIR, "assembly/{hybrid_assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna"), sr_sample="ONT3_MG_xx_Rashi_S11", barcode=BARCODES, hybrid_assembler=HYBRID_ASSEMBLER)
+    output:
+#           done=touch("bwa_index.done")
+           index_amb=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}.amb"),
+           index_ann=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}.ann"),
+           index_bwt=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}.bwt"),
+           index_pac=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}.pac"),
+           index_sa=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}.sa")
+#    wildcard_constraints:
+#          hybrid_assembler='^\/'
+    params:
+           prefix=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}")
+    log: os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{prefix}.bwa_index")
+    conda: "../envs/bwa.yaml"
+    shell: """
+           (date &&\
+           bwa index {input} -p {params.prefix} &&\
+           date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+           """
+
+rule hybrid_map_short_reads_to_hybrid_contigs_w_bwa_mem:
+    input:
+          query_r1=expand(os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_R1_001.fastp.fq.gz"), sr_sample="ONT3_MG_xx_Rashi_S11"),
+          query_r2=expand(os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_R2_001.fastp.fq.gz"), sr_sample="ONT3_MG_xx_Rashi_S11"),
+          index_amb=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.amb"),          
+          index_ann=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.ann"),
+          index_bwt=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.bwt"),
+          index_pac=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.pac"),
+          index_sa=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.sa"),
+          index_fa=expand(os.path.join(RESULTS_DIR, "assembly/{hybrid_assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna"), sr_sample="ONT3_MG_xx_Rashi_S11", barcode=BARCODES, hybrid_assembler=HYBRID_ASSEMBLER)
+    output:
+          os.path.join(RESULTS_DIR, "mapping/bwa_sr_{hybrid_assembler}/bwa_sr_{hybrid_assembler}.bam")
+    params:
+          index_prefix=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades"),
+          out_prefix=os.path.join(RESULTS_DIR, "mapping/bwa_sr_{hybrid_assembler}/bwa_sr_{hybrid_assembler}")
+    log: os.path.join(RESULTS_DIR, "mapping/bwa_sr_{hybrid_assembler}/sr_contigs.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    threads:
+          config["bwa"]["threads"]
+#    wildcard_constraints:
+#          sample='\w+'
+    shell: """
+           (date &&\
+           bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+           date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+           """
+
+rule hybrid_map_long_reads_to_hybrid_contigs_w_bwa_mem:
+    input:
+          query_lr=expand(os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"), barcode=BARCODES),
+          index_amb=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.amb"),
+          index_ann=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.ann"),
+          index_bwt=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.bwt"),
+          index_pac=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.pac"),
+          index_sa=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades.sa"),
+          index_fa=expand(os.path.join(RESULTS_DIR, "assembly/{hybrid_assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna"), sr_sample="ONT3_MG_xx_Rashi_S11", barcode=BARCODES, hybrid_assembler=HYBRID_ASSEMBLER)
+    output:
+          os.path.join(RESULTS_DIR, "mapping/bwa_lr_{hybrid_assembler}/bwa_lr_{hybrid_assembler}.bam")
+    params:
+          index_prefix=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/metaspades"),
+          out_prefix=os.path.join(RESULTS_DIR, "mapping/bwa_lr_{hybrid_assembler}/bwa_lr_{hybrid_assembler}"),
+    log:  os.path.join(RESULTS_DIR, "mapping/bwa_lr_{hybrid_assembler}/lr_contigs.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+#    wildcard_constraints:
+#          sample='\w+'
+    shell: """
+          (date &&\
+          bwa mem -x ont2d -t {threads} {params.index_prefix} {input.query_lr} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+          date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+          """
+
+rule hybrid_build_minimap2_long_read_contigs_ont_index:
+    input: expand(os.path.join(RESULTS_DIR, "assembly/{hybrid_assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna"), sr_sample="ONT3_MG_xx_Rashi_S11", barcode=BARCODES, hybrid_assembler=HYBRID_ASSEMBLER)
+    output: os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{hybrid_assembler}.mmi")
+    conda: "../envs/minimap2.yaml"
+    shell:
+           """
+            date
+            minimap2 -x map-ont -d {output} {input}
+            date
+            """
+
+rule hybrid_map_long_reads_to_hybrid_contigs_w_minimap2:
+    input:
+          query=expand(os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"), barcode=BARCODES),
+          index=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{hybrid_assembler}.mmi"),
+          index_fa=expand(os.path.join(RESULTS_DIR, "assembly/{hybrid_assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna"), sr_sample="ONT3_MG_xx_Rashi_S11", barcode=BARCODES, hybrid_assembler=HYBRID_ASSEMBLER)
+    output: os.path.join(RESULTS_DIR, "mapping/mmi_lr_{hybrid_assembler}/mmi_lr_{hybrid_assembler}.bam")
+    params:
+          prefix=os.path.join(RESULTS_DIR, "mapping/mmi_lr_{hybrid_assembler}/mmi_lr_{hybrid_assembler}")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/mmi_lr_{hybrid_assembler}/metaspades.minimap2_samtools")
+    shell: """
+            (date &&\
+            minimap2 -a -t {threads} {input.index} {input.query} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+rule hybrid_map_short_reads_to_hybrid_contigs_w_minimap2:
+    input:
+          query_r1=expand(os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_R1_001.fastp.fq.gz"), sr_sample="ONT3_MG_xx_Rashi_S11"),
+          query_r2=expand(os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_R2_001.fastp.fq.gz"), sr_sample="ONT3_MG_xx_Rashi_S11"),
+          index=os.path.join(RESULTS_DIR, "mapping/{hybrid_assembler}/{hybrid_assembler}.mmi"),
+          index_fa=expand(os.path.join(RESULTS_DIR, "assembly/{hybrid_assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna"), sr_sample="ONT3_MG_xx_Rashi_S11", barcode=BARCODES, hybrid_assembler=HYBRID_ASSEMBLER)
+    output: os.path.join(RESULTS_DIR, "mapping/mmi_sr_{hybrid_assembler}/mmi_sr_{hybrid_assembler}.bam")
+    params:
+          prefix=os.path.join(RESULTS_DIR, "mapping/mmi_sr_{hybrid_assembler}/mmi_sr_{hybrid_assembler}")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/mmi_sr_{hybrid_assembler}/metaspades.minimap2_samtools")
+    shell:
+           """
+            (date &&\
+           minimap2 -a -t {threads} {input.index} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output}
+           date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+rule hybrid_merge_bam_files_for_metaspades:
+    input:
+          bam_sr=os.path.join(RESULTS_DIR, "mapping/{mapper}_sr_{hybrid_assembler}/{mapper}_sr_{hybrid_assembler}.bam"),
+          bam_lr=os.path.join(RESULTS_DIR, "mapping/{mapper}_lr_{hybrid_assembler}/{mapper}_lr_{hybrid_assembler}.bam")
+    output: os.path.join(RESULTS_DIR, "mapping/{mapper}_merged_{hybrid_assembler}/{mapper}_merged_{hybrid_assembler}.bam")
+    conda: "../envs/binning.yaml"
+    shell:
+          """
+          (date &&\
+          samtools merge {output} {input.bam_sr} {input.bam_lr} &&\
+          date)
+          """
diff --git a/rules/METAT_RULES b/rules/METAT_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..636ea21c9649bfa82c31e132d14e56d9b5dca24f
--- /dev/null
+++ b/rules/METAT_RULES
@@ -0,0 +1,121 @@
+# Rules for mapping the metaT reads to the different assemblies
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mmi"]
+# SAMPLES=config["samples"]
+SAMPLES=["flye", "megathit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+###########
+## metaT ##
+###########
+# Preprocess the metaT reads using fastp
+rule run_fastp_on_metaT_reads:
+    input:
+        r1=os.path.join(config["metaT_prefix"], "{metaT_sample}_R1_001.fastq.gz"),
+        r2=os.path.join(config["metaT_prefix"], "{metaT_sample}_R2_001.fastq.gz")
+    output:
+        r1=os.path.join(RESULTS_DIR, "preprocessing/metaT/{metaT_sample}_R1_001.fastp.fq.gz"),
+        r2=os.path.join(RESULTS_DIR, "preprocessing/metaT/{metaT_sample}_R2_001.fastp.fq.gz"),
+        html=os.path.join(RESULTS_DIR, "preprocessing/metaT/{metaT_sample}.fastp.html"),
+        json=os.path.join(RESULTS_DIR, "preprocessing/metaT/{metaT_sample}.fastp.json")
+    log: os.path.join(RESULTS_DIR, "preprocessing/metaT/{metaT_sample}.fastp")
+    conda: "../envs/fastp.yaml"
+    shell:
+        """
+        (date &&\
+        fastp -l {config[fastp][min_length]} -i {input.r1} -I {input.r2} -o {output.r1} -O {output.r2} -h {output.html} -j {output.json} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule run_fastqc_on_preprocessed_metaT_reads:
+    input: os.path.join(RESULTS_DIR, "preprocessing/metaT/{metaT_sample}_{suffix}.fq.gz"),
+    output:
+        html=os.path.join(RESULTS_DIR, "qc/metaT/fastqc/{metaT_sample}/{metaT_sample}_{suffix}.fastqc.html"),
+        zip=os.path.join(RESULTS_DIR, "qc/metaT/fastqc/{metaT_sample}/{metaT_sample}_{suffix}.fastqc.zip") # the suffix _fastqc.zip is necessary for multiqc to find the file. If not using multiqc, you are free to choose an arbitrary filename
+    conda: "../envs/fastqc.yaml"
+    script: "../scripts/run_fastqc.py" # Use a separate script rather than a simple call to FASTQC to overcome possible race conditions when running in parallel (https://bitbucket.org/snakemake/snakemake-wrappers/raw/7f5e01706728e32d52aa413f213e950f50bf4129/bio/fastqc/wrapper.py)
+
+
+rule map_metaT_reads_to_long_read_contigs_w_bwa_mem:
+    input:
+        query_r1=rules.run_fastp_on_metaT_reads.output.r1,
+        query_r2=rules.run_fastp_on_metaT_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/metaT/lr/{metaT_sample}_reads-x-{barcode}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/metaT/lr/{metaT_sample}_reads-x_{barcode}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/metaT/lr/{metaT_sample}_reads-x-{barcode}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+
+rule map_metaT_reads_to_short_read_contigs_w_bwa_mem:
+    input:
+        query_r1=rules.run_fastp_on_metaT_reads.output.r1,
+        query_r2=rules.run_fastp_on_metaT_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/metaT/sr/{metaT_sample}_reads-x-{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/metaT/sr/{metaT_sample}_reads-x-{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/metaT/sr/{metaT_sample}_reads-x-{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+
+rule map_metaT_reads_to_hybrid_contigs_w_bwa_mem:
+    input:
+        query_r1=rules.run_fastp_on_metaT_reads.output.r1,
+        query_r2=rules.run_fastp_on_metaT_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/metaT/sr/{metaT_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/metaT/sr/{metaT_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs")
+    log: os.path.join(RESULTS_DIR, "mapping/metaT/sr/{metaT_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
diff --git a/rules/MMSEQ_RULES b/rules/MMSEQ_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..6e805348a14e05a953e0c2c102325d3209428f16
--- /dev/null
+++ b/rules/MMSEQ_RULES
@@ -0,0 +1,114 @@
+# For running the MMSEQ2 comparison of proteins after assemblies are run through prokka/prodigal
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mm"]
+# SAMPLES=config["samples"]
+SAMPLES=["flye", "megahit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+#############################
+######## MMSEQ2  ############
+#############################
+rule mmseq2_database:
+    input:
+        int1=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/lr/merged/{barcode}/assembly.faa"), assembler="flye", barcode=BARCODES),
+        int2=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/{sr_sample}/final.contigs.faa"), assembler="megahit", sr_sample="ONT3_MG_xx_Rashi_S11"),
+        int3=expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.faa"), assembler="metaspades_hybrid", barcode=BARCODES, sr_sample="ONT3_MG_xx_Rashi_S11")
+#        expand(os.path.join(RESULTS_DIR, "annotation/proteins/{assembler}/{prefix}.faa"), assembler=["flye", "megahit", "metaspades_hybrid"])
+    output:
+        dat1=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="flye"),
+        dat2=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="megahit"),
+        dat3=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="metaspades_hybrid")
+    log: os.path.join(RESULTS_DIR, "annotation/mmseq2/database.mmseq2.log")
+#    conda: "../envs/cd-hit.yaml"
+    shell:
+        """
+        (date &&\
+        {config[mmseqs][createdb]} {input.int1} {output.dat1} &&\
+        {config[mmseqs][createdb]} {input.int2} {output.dat2} &&\
+        {config[mmseqs][createdb]} {input.int3} {output.dat3} &&\
+        date) &> >(tee {log})
+        """
+
+rule mmseq2_compare:
+    input:
+        mm1=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="flye"),
+        mm2=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="megahit"),
+        mm3=expand(os.path.join(RESULTS_DIR, "annotation/mmseq2/{assembler}_db"), assembler="metaspades_hybrid")
+    output:
+        mo1=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_megahit_rbh"),
+        mo2=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_metaspades_hybrid_rbh"),
+        mo3=os.path.join(RESULTS_DIR, "annotation/mmseq2/megahit_metaspades_hybrid_rbh")
+    log: os.path.join(RESULTS_DIR, "annotation/mmseq2/compare.mmseq2.log")
+    threads: config["mmseq2"]["threads"]
+#    conda: "../envs/cd-hit.yaml"
+    shell:
+        """
+        (date &&\
+        {config[mmseqs][rbh]} {input.mm1} {input.mm2} {output.mo1} --min-seq-id 0.9 mmseq2_tmp --threads {threads} &&\
+        {config[mmseqs][rbh]} {input.mm1} {input.mm3} {output.mo2} --min-seq-id 0.9 mmseq2_tmp --threads {threads} &&\
+        {config[mmseqs][rbh]} {input.mm2} {input.mm3} {output.mo3} --min-seq-id 0.9 mmseq2_tmp --threads {threads} &&\
+        date) &> >(tee {log})
+        """
+
+rule mmseq2_m8_format:
+    input:
+        fo1=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_db"),
+        fo2=os.path.join(RESULTS_DIR, "annotation/mmseq2/megahit_db"),
+        fo3=os.path.join(RESULTS_DIR, "annotation/mmseq2/metaspades_hybrid_db"),
+        fo4=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_megahit_rbh"),
+        fo5=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_metaspades_hybrid_rbh"),
+        fo6=os.path.join(RESULTS_DIR, "annotation/mmseq2/megahit_metaspades_hybrid_rbh")
+    output:
+        form1=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_megahit.m8"),
+        form2=os.path.join(RESULTS_DIR, "annotation/mmseq2/flye_metaspades_hybrid.m8"),
+        form3=os.path.join(RESULTS_DIR, "annotation/mmseq2/megahit_metaspades_hybrid.m8")
+    log: os.path.join(RESULTS_DIR, "annotation/mmseq2/convertalis.mmseq2.log")
+#    conda: "../envs/cd-hit.yaml"
+    shell:
+        """
+        (date &&\
+        {config[mmseqs][convertalis]} {input.fo1} {input.fo2} {input.fo4} {output.form1}  &&\
+        {config[mmseqs][convertalis]} {input.fo1} {input.fo3} {input.fo5} {output.form2} &&\
+        {config[mmseqs][convertalis]} {input.fo2} {input.fo3} {input.fo6} {output.form3} &&\
+        date) &> >(tee {log})
+        """
+
+rule prepare_plot_files:
+    input: rules.mmseq2_m8_format.output
+    output:
+        touch(os.path.join(RESULTS_DIR, "annotation/mmseq2/plot_files_ready.done")),
+        prep1=os.path.join(RESULTS_DIR, "annotation/mmseq2/total.txt"), 
+        prep2=os.path.join(RESULTS_DIR, "annotation/mmseq2/overlap_sizes.txt")
+    log: os.path.join(RESULTS_DIR, "annotation/mmseq2/files_ready.mmseq2.log")
+    shell:
+        """
+        (date &&\
+        scripts/prepare_plot_files.sh &&\
+        date) &> >(tee {log})
+        """
+
+rule mmseq2_plots:
+    input:
+        plot1=os.path.join(RESULTS_DIR, "annotation/mmseq2/overlap_sizes.txt"),
+        plot2=os.path.join(RESULTS_DIR, "annotation/mmseq2/total.txt")
+    output:
+        touch(os.path.join(RESULTS_DIR, "annotation/mmseq2/upset_plots.done"))
+    log: os.path.join(RESULTS_DIR, "annotation/mmseq2/plot.mmseq2.log")
+    conda: "../envs/renv.yaml"
+    shell:
+        """
+        (date &&\
+        Rscript scripts/mmseq_plots.R {input.plot1} {input.plot2} &&\
+        date) &> >(tee {log})
+        """
diff --git a/rules/TAXONOMY_RULES b/rules/TAXONOMY_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..4be8354f58be5ff5229b19d6490f78980434c752
--- /dev/null
+++ b/rules/TAXONOMY_RULES
@@ -0,0 +1,73 @@
+# For running the contamination check and taxonomy on the generated MAGs
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+BARCODES=config["barcodes"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa", "mm"]
+#SAMPLES=config["samples"]
+SAMPLES=["flye", "megahit", "metaspades_hybrid"]
+BINNING_SAMPLES=config["binning_samples"]
+HYBRID_ASSEMBLER=config["hybrid_assembler"]
+
+
+#####################################
+######## GTDBTK & CHECKM ############
+#####################################
+rule gtdbtk:
+    input:
+           gen=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_DASTool_bins")
+    threads:config["threads"]
+    conda:"../envs/gtdbtk.yaml"
+    params:
+           outdir=os.path.join(RESULTS_DIR, "Binning/gtdbtk_output/{sample}")
+    output:
+#           gtdb=directory(os.path.join(RESULTS_DIR, "Binning/gtdbtk_output/{sample}")),
+           gtdbfile=os.path.join(RESULTS_DIR, "Binning/gtdbtk_output/{sample}/gtdbtk.bac120.summary.tsv")
+    wildcard_constraints:
+           sample='\w+'
+    shell:
+            """
+            (date &&\
+            export GTDBTK_DATA_PATH={config[GTDBTK][DATA]} &&\
+            gtdbtk classify_wf --cpus {threads} -x fa --genome_dir {input.gen} --out_dir {params.outdir} &&\
+            touch gtdbtk.done &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+rule checkm:
+    input:
+           chkm=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_DASTool_bins")
+    threads:config["threads"]
+    conda:"../envs/checkm.yaml"
+    output:
+           chkmdir=directory(os.path.join(RESULTS_DIR, "Binning/checkm_output/{sample}/lineage.ms")),
+           chkout=os.path.join(RESULTS_DIR, "Binning/checkm_output/{sample}_output.txt")
+    shell:
+            """
+            (date &&\
+            checkm lineage_wf -r -t {threads} -x fa {input.chkm} {output.chkmdir}  | tee {output.chkout} &&\
+            touch checkm.done &&\
+            date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+            """
+
+rule plot_checkm:
+    input:
+           pl1=os.path.join(RESULTS_DIR, "Binning/checkm_output/{sample}"),
+           pl2=os.path.join(RESULTS_DIR, "Binning/{sample}/dastool_output/{sample}_DASTool_bins")
+    threads:config["threads"]
+    conda:"../envs/checkm.yaml"
+    output:
+           plt=directory(os.path.join(RESULTS_DIR, "Binning/checkm_output/{sample}/plots"))
+    shell:
+         """
+         date
+         checkm bin_qa_plot --image_type png --width 8 --dpi 200 --font_size 8 -x fa -q {input.pl1} {input.pl2} {output.plt}
+         touch plot_checkm.done
+         date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+         """
diff --git a/rules/checkpoint_ASSEMBLY_ANNOTATION_RULES b/rules/checkpoint_ASSEMBLY_ANNOTATION_RULES
new file mode 100755
index 0000000000000000000000000000000000000000..2313f31da921edc67dd669e6aa3dedeb183089ee
--- /dev/null
+++ b/rules/checkpoint_ASSEMBLY_ANNOTATION_RULES
@@ -0,0 +1,1151 @@
+# For running the ASSEMBLY and ANNOTATION workflow for ONT data
+
+import os
+from tempfile import TemporaryDirectory
+
+configfile: "config/CONFIG.yaml"
+DATA_DIR = config["data_dir"]
+RESULTS_DIR = config["results_dir"]
+DB_DIR=config["db_dir"]
+RUNS=[config["runs"]["first"],
+    config["runs"]["second"]]
+#    config["runs"]["third"]]
+BARCODES=config["barcodes"]
+SAMPLES=["ONT3_MG_xx_Rashi_S11"]
+sr_sample=["ONT3_MG_xx_Rashi_S11"]
+# SAMPLES_ALL=config["samples"]
+REFERENCES=["igc", "hg38"]
+IGC_URI=config["igc"]["uri"]
+HG38_URI=config["hg38"]["uri"]
+ASSEMBLERS=config["assemblers"]
+MAPPERS=["bwa_mem", "minimap2"]
+
+
+######
+# RULES
+######       
+# rule unpack:
+#     input: os.path.join(DATA_DIR, "2018-11-08-minion_latestdata.7z")
+#     output: os.path.join(DATA_DIR, "raw/raw_unpacked.done")
+#     threads: config["p7zip"]["threads"]
+#     shell:
+#         """
+#         date
+#         {config[p7zip][bin]} x -o$(dirname {output}) {input} -mmt{threads} &> {output}
+#         date
+#         """
+
+# Take the per-read FAST5 files and create multi-FAST5 files including a batch of reads (typically not all reads, though).
+#rule create_multifast5s:
+#    input: os.path.join(DATA_DIR, "raw/{run}")
+#    output: protected(directory(os.path.join(DATA_DIR, "multifast5/{run}")))
+#    threads: config["ont_fast5_api"]["single_to_multi_fast5"]["threads"]
+#    log: os.path.join(DATA_DIR, "multifast5/{run}/create_multifast5s")
+#    shell:
+#        """
+#        (date &&\
+#        {config[ont_fast5_api][single_to_multi_fast5][bin]} --input_path {input} --save_path {output} --filename_base multifast5 --batch_size {config[ont_fast5_api][single_to_multi_fast5][batch]} --recursive -t {threads} &&\
+#        touch {output} &&\
+#        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#        """
+
+# FLO-MIN107 SQK-LSK108           dna_r9.5_450bps
+#   --verbose_logs                    Enable verbose logs.
+#   --qscore_filtering                Enable filtering of reads into PASS/FAIL
+#                                     folders based on min qscore.
+#   --min_qscore arg                  Minimum acceptable qscore for a read to be
+#                                     filtered into the PASS folder
+#   --disable_pings                   Disable the transmission of telemetry
+#   --num_callers arg                 Number of parallel basecallers to create.
+#   --num_barcode_threads arg         Number of worker threads to use for
+#                                     barcoding.
+#   --barcode_kits arg                Space separated list of barcoding kit(s) or
+#                                     expansion kit(s) to detect against. Must be
+#                                     in double quotes.
+#   --trim_barcodes                   Trim the barcodes from the output sequences
+#                                     in the FastQ files.
+
+# Left in for legacy look-up reasons. GPU-based basecalling strongly encouraged.
+# rule guppy_cpu_basecall:
+#     input: os.path.join(DATA_DIR, "multifast5/{run}")
+#     output: os.path.join(RESULTS_DIR, "basecalled/guppy/cpu/{run}/basecalling.done")
+#     threads: config["guppy_cpu"]["cpu_threads"]
+#     shell:
+#         """
+#         date
+#         {config[guppy_cpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_cpu][config]} --cpu_threads_per_caller {threads} --disable_pings
+#         touch {output}
+#         date
+#         """
+
+    #input: os.path.join(DATA_DIR, "multifast5/{run}")
+# Basecalling the multi-FAST5 files.
+# Currently requires dummy target, but would be nice to have this be a "real" file target. However, it seems that files are filled incrementally, hence there is no check if the whole process was successful or not if checking only for file presence.
+checkpoint guppy_gpu_basecall:
+    input: os.path.join(DATA_DIR, "multifast5/{run}")
+    output: directory(os.path.join(RESULTS_DIR, "basecalled/{run}/checkpoint_bc"))
+    log: os.path.join(RESULTS_DIR, "basecalled/{run}/basecalling")
+    threads: config["guppy_gpu"]["cpu_threads"]
+    shell:
+        """
+        (date &&\
+        {config[guppy_gpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_gpu][config]} --cpu_threads_per_caller {threads} -x {config[guppy_gpu][gpu_device]} --disable_pings --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --chunk_size {config[guppy_gpu][chunk_size]} --chunks_per_runner {config[guppy_gpu][chunks_per_runner]} --gpu_runners_per_device {config[guppy_gpu][runners_per_device]} --num_callers {config[guppy_gpu][num_callers]} && \
+        touch {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+checkpoint guppy_gpu_basecall_NO_MOD:
+    input: os.path.join(DATA_DIR, "multifast5/{run}")
+    output: directory(os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/checkpoint_bc_NM"))
+    log: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/NO_MOD_basecalling")
+    wildcard_constraints:
+#        run='S.*'
+#        run='S.*^\/'
+#        run='\w+'
+#        run="^no.*"
+    threads: config["guppy_gpu"]["cpu_threads"]
+    shell:
+        """
+        (date &&\
+        {config[guppy_gpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_gpu][hac_config]} --cpu_threads_per_caller {threads} -x {config[guppy_gpu][gpu_device]} --disable_pings --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --chunk_size {config[guppy_gpu][chunk_size]} --chunks_per_runner {config[guppy_gpu][chunks_per_runner]} --gpu_runners_per_device {config[guppy_gpu][runners_per_device]} --num_callers {config[guppy_gpu][num_callers]} && \
+        touch {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+#rule guppy_gpu_basecall:
+#    input: os.path.join(DATA_DIR, "multifast5/{run}")
+#    output: os.path.join(RESULTS_DIR, "basecalled/{run}/basecalling.done")
+#    log: os.path.join(RESULTS_DIR, "basecalled/{run}/basecalling")
+#    threads: config["guppy_gpu"]["cpu_threads"]
+#    shell:
+#        """
+#        (date &&\
+#        {config[guppy_gpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_gpu][config]} --cpu_threads_per_caller {threads} -x {config[guppy_gpu][gpu_device]} --disable_pings --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --chunk_size {config[guppy_gpu][chunk_size]} --chunks_per_runner {config[guppy_gpu][chunks_per_runner]} --gpu_runners_per_device {config[guppy_gpu][runners_per_device]} --num_callers {config[guppy_gpu][num_callers]} && \
+#        touch {output} &&\
+#        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#        """
+#
+#rule guppy_gpu_basecall_NO_MOD:
+#    input: os.path.join(DATA_DIR, "multifast5/{run}")
+#    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/NO_MOD_basecalling.done")
+#    log: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/NO_MOD_basecalling")
+#    threads: config["guppy_gpu"]["cpu_threads"]
+#    shell: 
+#        """
+#        (date &&\
+#        {config[guppy_gpu][bin]} --input_path {input} --save_path $(dirname {output}) --config {config[guppy_gpu][hac_config]} --cpu_threads_per_caller {threads} -x {config[guppy_gpu][gpu_device]} --disable_pings --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --chunk_size {config[guppy_gpu][chunk_size]} --chunks_per_runner {config[guppy_gpu][chunks_per_runner]} --gpu_runners_per_device {config[guppy_gpu][runners_per_device]} --num_callers {config[guppy_gpu][num_callers]} && \
+#        touch {output} &&\
+#        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#        """
+
+# Demultiplex.
+# Currently requires dummy target, but would be nice to have this be a "real" file target. However, it seems that files are filled incrementally, hence there is no check if the whole process was successful or not if checking only for file presence.
+#rule guppy_gpu_demux:
+#    input: rules.guppy_gpu_basecall.output
+#    output: os.path.join(RESULTS_DIR, "basecalled/{run}/demux.done")
+#    log: os.path.join(RESULTS_DIR, "basecalled/{run}/demux")
+#    threads: config["guppy_barcoder"]["threads"]
+#    shell:
+#        """
+#        (date &&\
+#        {config[guppy_barcoder][bin]} --input_path $(dirname {input}) --save_path $(dirname {output}) -t {threads} --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --trim_barcodes && \
+#        touch {output} &&\
+#        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#        """
+#rule guppy_gpu_demux_NO_MOD:
+#    input: rules.guppy_gpu_basecall_NO_MOD.output
+#    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/demux_NO_MOD.done")
+#    log: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/demux")
+#    threads: config["guppy_barcoder"]["threads"]
+#    shell:
+#        """
+#        (date &&\
+#        {config[guppy_barcoder][bin]} --input_path $(dirname {input}) --save_path $(dirname {output}) -t {threads} --compress_fastq --records_per_fastq {config[guppy_gpu][records_per_fastq]} --trim_barcodes && \
+#        touch {output} &&\
+#        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#        """
+
+# Modifications (IMPORTANT)
+# Since 2019_GDB sample was "NOT" barcoded the demux folders are obsolete.
+# The below section gathers all the fastq.gz files from basecalling into a "dummy" folder
+# The "dummy" folder will be used downstream similar to the "barcode" folder used in the 2018 analyses
+# The barcode in the CONFIG file will be given the name of this dummy folder
+rule intermediate_basecalling:
+    input: os.path.join(RESULTS_DIR, "basecalled/{run}/checkpoint_bc/{i}.fastq.gz")
+    output: os.path.join(RESULTS_DIR, "basecalled/{run}/no_nobarcode/{i}.fastq.gz")
+    log: os.path.join(RESULTS_DIR, "basecalled/{run}/no_barcode_{i}")
+    shell:
+        """
+        (date &&\
+        ln -s {input} {output}  &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule intermediate_NO_MOD_basecalling:
+    input: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/checkpoint_bc_NM/{i}.fastq.gz")
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/no_nobarcode/{i}.fastq.gz")
+    log: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/no_barcode_{i}")
+    shell:
+        """
+        (date &&\
+        ln -s {input} {output}  &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+def aggregate_dummy_basecalling(wildcards):
+    checkpoint_output = checkpoints.guppy_gpu_basecall.get(**wildcards).output[0]
+    return expand(os.path.join(RESULTS_DIR, "basecalled/{run}/no_nobarcode/{id}.fastq.gz"),
+        run=wildcards.run,
+        i=glob_wildcards(os.path.join(checkpoint_output, "{i}.fastq.gz")).i)
+
+
+def aggregate_dummy_basecalling_NO_MOD(wildcards):
+    dummy_output=checkpoints.guppy_gpu_basecall_NO_MOD.get(**wildcards).output[0]
+    return expand(os.path.join(RESULTS_DIR, "basecalled/{run}/no_nobarcode/{i}.fastq.gz"),
+        run=wildcards.run,
+        i=glob_wildcards(os.path.join(dummy_output, "{i}.fastq.gz")).i)
+
+# Merge multiple basecalled files, i.e., FASTQ files per barcode and run. Each multi-FAST5 file will produce a respective FASTQ file.
+# TODO: Simplify output path to `basecalled/{run}/{barcode}.fastq.gz`. This will require to change `$(dirname {output})` to `$(dirname {output})/{wildcards.barcode}`
+rule merge_individual_fastq_per_barcode:
+#    input: rules.guppy_gpu_demux.output
+    input: aggregate_dummy_basecalling
+    output: os.path.join(RESULTS_DIR, "basecalled/{run}/{barcode}/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat $(find $(dirname {output}) -name "*.fastq.gz" | sort) > {output}
+        touch {output}
+        date
+        """
+
+rule merge_individual_fastq_per_barcode_NO_MOD:
+#    input: rules.guppy_gpu_demux_NO_MOD.output
+    input: aggregate_dummy_basecalling_NO_MOD
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/{barcode}/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat $(find $(dirname {output}) -name "*.fastq.gz" | sort) > {output}
+        touch {output}
+        date
+        """
+
+# Helper function to get the input in order to create a *single* fastq file for all runs with the *same* barcode
+# TODO: Adjust according to TODO-merge_individual_fastq_per_barcode
+def get_all_fastq_for_barcode(wildcards):
+    return expand(os.path.join(RESULTS_DIR, "basecalled/{run}/{barcode}/{barcode}.fastq.gz"), run=RUNS, barcode=wildcards.barcode)
+
+def get_all_fastq_for_barcode_NO_MOD(wildcards):
+    return expand(os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/{barcode}/{barcode}.fastq.gz"), run=RUNS, barcode=wildcards.barcode)
+
+# Merge basecalled files (FASTQ) per barcode, i.e., combine reads from multiple runs according to their barcode.
+rule get_single_fastq_per_barcode:
+    input: get_all_fastq_for_barcode 
+    output: os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat {input} > {output}
+        date
+        """
+rule get_single_fastq_per_barcode_NO_MOD:
+    input: get_all_fastq_for_barcode_NO_MOD
+    output: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/merged/{barcode}.fastq.gz")
+    shell:
+        """
+        date
+        cat {input} > {output}
+        date
+        """
+
+# Get basic statistics of the ONT data; over all runs
+# TODO: Simplify top of path: replace `lr` by `nanostat`
+# TODO: Write more generic: have patterned `input: os.path.join(RESULTS_DIR, "basecalled/{prefix}/{barcode}.fastq.gz")` and patterned `output: os.path.join(RESULTS_DIR, "qc/nanostat/{prefix}/{barcode}/{barcode}NanoStats.txt")
+# 
+rule nanostat:
+    input: rules.get_single_fastq_per_barcode.output
+    output: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats.txt")
+    log: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        (date &&\
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output}) &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule nanostat_NO_MOD:
+    input: rules.get_single_fastq_per_barcode_NO_MOD.output
+    output: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats_NO_MOD.txt")
+    log: os.path.join(RESULTS_DIR, "qc/lr/merged/{barcode}/{barcode}NanoStats_NO_MOD")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        (date &&\
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output}) &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Deprecated if TODO-nanostat is completed
+# Get basic statistics of the ONT data; per run
+rule nanostat_per_run:
+    input: os.path.join(RESULTS_DIR, "basecalled/{run}/{barcode}/{barcode}.fastq.gz")
+    output: os.path.join(RESULTS_DIR, "qc/lr/{run}/{barcode}/{barcode}NanoStats.txt")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        date
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output})
+        date
+        """
+
+rule nanostat_per_run_NO_MOD:
+    input: os.path.join(RESULTS_DIR, "basecalled_NO_MOD/{run}/{barcode}/{barcode}.fastq.gz")
+    output: os.path.join(RESULTS_DIR, "qc/lr/{run}/{barcode}/{barcode}NanoStats_NO_MOD.txt")
+    conda: "../envs/nanostat.yaml"
+    shell:
+        """
+        date
+        NanoStat --fastq {input} --outdir $(dirname {output}) --prefix {wildcards.barcode} -n $(basename -s "" {output})
+        date
+        """
+
+# Download the IGC (http://www.nature.com/nbt/journal/vaop/ncurrent/full/nbt.2942.html)
+# Downloading using wget-command instead of FTP-remote in snakemake as the former offers progress information
+rule fetch_igc:
+    output: protected(os.path.join(DB_DIR, "igc/igc.fna.gz"))
+    log: os.path.join(DB_DIR, "igc/fetch_igc")
+    shell:
+        """
+        (date &&\
+        wget ftp://{IGC_URI} -O {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Download HG38 human reference genome
+# Downloading using wget-command instead of FTP-remote in snakemake as the former offers progress information
+rule fetch_hg38:
+    output: protected(os.path.join(DB_DIR, "hg38/GCF_000001405.38_GRCh38.p12_genomic.fna.gz"))
+    log: os.path.join(DB_DIR, "hg38/fetch_hg38")
+    shell:
+        """
+        (date &&\
+        wget ftp://{HG38_URI} -O {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule create_link_to_hg38:
+    input: os.path.join(DB_DIR, "hg38/GCF_000001405.38_GRCh38.p12_genomic.fna")
+    output: protected(os.path.join(DB_DIR, "hg38/hg38.fna"))
+    shell:
+        """
+        date
+        ln -s $(basename -s "" {input}) {output}
+        date
+        """
+
+rule unpack_reference_data:
+    input: os.path.join(DB_DIR, "{prefix}.fna.gz")
+    output: protected(os.path.join(DB_DIR, "{prefix}.fna"))
+    shell:
+        """
+        date
+        gzip -dc {input} > {output}
+        date
+        """
+
+# Build a minimap2 index to map *ONT* reads against
+rule build_minimap2_ont_index:
+    input: "{prefix}.fna"
+    output: "{prefix}-ont.mmi"
+    log: "{prefix}.minimap2_ont_index"
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -x map-ont -d {output} {input} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Build a minimap2 short read index in order to map short reads to a reference/long read-based contigs
+# NOTE: Unclear whether this would be also the right option for mapping short reads to long reads (i.e., reads-to-reads), due to the increased error rate in long reads. For more discussions on this, see https://github.com/lh3/minimap2/issues/407#issuecomment-493823000
+# TODO: Combine this rule and `build_minimap2_ont_index` using a conditional on the suffix (ont vs. sr). This makes the code more readable as it removes largely redundant `input` and `output` keywords
+rule build_minimap2_sr_index:
+    input: "{prefix}.fna"
+    output: "{prefix}-sr.mmi"
+    log: "{prefix}.minimap2_sr_index"
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -x sr -d {output} {input} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Build a bwa index (may be used to map short reads against it, or long reads by using the "-x ont2d" flag)
+rule build_bwa_index:
+    input: "{prefix}.fna"
+    output: index_amb="{prefix}.amb",
+       index_ann="{prefix}.ann",
+       index_bwt="{prefix}.bwt",
+       index_pac="{prefix}.pac",
+       index_sa="{prefix}.sa"
+    log: "{prefix}.bwa_index"
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa index {input} -p {wildcards.prefix} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Simple utility function
+#rule convert_fastq_to_fasta:
+#    input: "{prefix}.fastq.gz"
+#    output: "{prefix}.fna"
+#    shell:
+#        """
+#        date
+#        zcat {input} | sed -n '1~4s/^@/>/p;2~4p' > {output}
+#        date
+#        """
+#
+#rule convert_fq_to_fasta:
+#    input: "{prefix}.fq.gz"
+#    output: "{prefix}.fna"
+#    shell:
+#        """
+#        date
+#        zcat {input} | sed -n '1~4s/^@/>/p;2~4p' > {output}
+#        date
+#        """
+
+# Map long reads to a reference (IGC/HG38/other). The resulting BAM file can be use for many things, among other to estimate which long reads are "known" or "novel", or to estimate the coverage of reference sequences (regions thereof).
+rule map_long_reads_to_reference_w_minimap2:
+    input: 
+        query=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        index=os.path.join(DB_DIR, "{reference}/{reference}-ont.mmi"),
+        index_fa=os.path.join(DB_DIR, "{reference}/{reference}.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}-x-{reference}.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}-x-{reference}")
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}-x-{reference}.minimap2_samtools")
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query} | samtools view -bT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Preprocess the short reads using fastp
+rule run_fastp_on_short_reads:
+    input:
+        r1=os.path.join(config["short_reads_prefix"], "ONT3_MG_xx_Rashi_S11_R1_001.fastq.gz"),
+        r2=os.path.join(config["short_reads_prefix"], "ONT3_MG_xx_Rashi_S11_R2_001.fastq.gz")
+    output:
+        r1=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11_R1_001.fastp.fq.gz"),
+        r2=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11_R2_001.fastp.fq.gz"),
+        html=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11.fastp.html"),
+        json=os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11.fastp.json")
+    log: os.path.join(RESULTS_DIR, "preprocessing/sr/ONT3_MG_xx_Rashi_S11.fastp")
+    conda: "../envs/fastp.yaml"
+    shell:
+        """
+        (date &&\
+        fastp -l {config[fastp][min_length]} -i {input.r1} -I {input.r2} -o {output.r1} -O {output.r2} -h {output.html} -j {output.json} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Simplify paths by removing the `sr`.
+rule run_fastqc_on_preprocessed_short_reads:
+    input: os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_{suffix}.fq.gz"),
+    output:
+        html=os.path.join(RESULTS_DIR, "qc/sr/fastqc/{sr_sample}/{sr_sample}_{suffix}.fastqc.html"),
+        zip=os.path.join(RESULTS_DIR, "qc/sr/fastqc/{sr_sample}/{sr_sample}_{suffix}.fastqc.zip") # the suffix _fastqc.zip is necessary for multiqc to find the file. If not using multiqc, you are free to choose an arbitrary filename
+    conda: "../envs/fastqc.yaml"
+    script: "../scripts/run_fastqc.py" # Use a separate script rather than a simple call to FASTQC to overcome possible race conditions when running in parallel (https://bitbucket.org/snakemake/snakemake-wrappers/raw/7f5e01706728e32d52aa413f213e950f50bf4129/bio/fastqc/wrapper.py)
+
+rule map_short_reads_to_reference_w_minimap2:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index=os.path.join(DB_DIR, "{reference}/{reference}-sr.mmi"),
+        index_fa=os.path.join(DB_DIR, "{reference}/{reference}.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}-x-{reference}.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}-x-{reference}")
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}-x-{reference}.minimap2_samtools")
+    conda: "../envs/minimap2.yaml"
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -bT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_short_reads_to_reference_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(DB_DIR, "{reference}/{reference}.amb"),
+        index_ann=os.path.join(DB_DIR, "{reference}/{reference}.ann"),
+        index_bwt=os.path.join(DB_DIR, "{reference}/{reference}.bwt"),
+        index_pac=os.path.join(DB_DIR, "{reference}/{reference}.pac"),
+        index_sa=os.path.join(DB_DIR, "{reference}/{reference}.sa"),
+        index_fa=os.path.join(DB_DIR, "{reference}/{reference}.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}-x-{reference}.bam")
+    params:
+        index_prefix=os.path.join(DB_DIR, "{reference}/{reference}"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}-x-{reference}")
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}-x-{reference}.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -bT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Compute the genome coverage histogram - https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html
+# The resulting histogram is used to compute the expected, i.e., average, coverage of the underyling reference (which may be de novo contigs).
+rule run_genomecov:
+    input: os.path.join(RESULTS_DIR, "mapping/{prefix}.bam")
+    output: os.path.join(RESULTS_DIR, "genomecov/{prefix}.cov.txt")
+    conda: "../envs/bedtools.yaml"
+    shell:
+        """
+        date
+        bedtools genomecov -ibam {input} > {output}
+        date
+        """
+
+# Compute the expected, i.e., average, coverage given the coverage histogram
+rule compute_avg_genomecov:
+    input: os.path.join(RESULTS_DIR, "genomecov/{prefix}.cov.txt")
+    output: os.path.join(RESULTS_DIR, "genomecov/{prefix}.avg_cov.txt")
+    shell:
+        """
+        date
+        cat {input} | awk -f {config[compute_avg_coverage][bin]} | tail -n+2 > {output}
+        date
+        """
+
+# TODO: All above rules that build an index build it in the same path as the reference/target file
+# In fact, it is redundant with `build_bwa_index`. The only advantage is that this replication allows more fine-grained control of the resources (slurm, etc.). This is, however, probably not needed as `build_bwa_index` can be quite resource-intensive, e.g., for the IGC.
+# Similar to building a minimap2 index, but this time for bwa
+rule build_long_read_bwa_index:
+    input: rules.get_single_fastq_per_barcode.output
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.amb"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.ann"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.bwt"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.pac"),
+            os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.sa")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}")
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.bwa_index")
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa index {input} -p {params.prefix} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+# TODO: If TODO-build_long_read_bwa_index is completed, need to adjust the `input.index` here too
+# config[bwa][long_reads_index][opts] according to https://github.com/sfu-compbio/colormap/blob/baa6802b43da8640c305032b7f7fcadd7c496a03/run_corr.sh#L45 and https://academic.oup.com/bioinformatics/article/32/17/i545/2450793#112481513 -> "2.2.3 Mapping parameters"
+rule map_short_reads_to_long_reads:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}.amb"),
+        index_fa=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fna"),
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{sr_sample}-x-{barcode}.bam")
+    threads: config["bwa"]["threads"]
+    wildcard_constraints:
+        barcode='^\-' 
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{sr_sample}-x-{barcode}.")
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{sr_sample}-x-{barcode}.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        (date &&\
+        bwa mem {config[bwa][long_reads_index][opts]} -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -S -b -u -T {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Add extra rules for other assemblers instead of a generic rule as this offers then greater flexibility, e.g., in terms of resource-reservation.
+rule assemble_long_reads_w_flye:
+    input: rules.get_single_fastq_per_barcode.output
+    output: os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna")
+    threads: config["flye"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.flye")
+#    conda: "../envs/flye.yaml"
+    conda: "../envs/flye_v2_7.yaml"
+    shell:
+        """
+        (date &&\
+        flye --nano-raw {input} --meta --out-dir $(dirname {output}) --genome-size {config[flye][genome_size]} --threads {threads} &&\
+        ln -s assembly.fasta {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# "r941_min_high" is the MinION model, high accuarcy
+rule polish_long_read_contigs_w_medaka:
+    input: 
+        basecalls=rules.get_single_fastq_per_barcode.output,
+        contigs=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/medaka/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    threads: config["medaka"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/medaka/round_{round}/lr/merged/{barcode}/assembly_polished.medaka")
+    conda: "../envs/medaka.yaml"
+    shell:
+        """
+        (date &&\
+        medaka_consensus -i {input.basecalls} -d {input.contigs} -o $(dirname {output}) -t {threads} -m r941_min_high &&\
+        ln -s consensus.fasta {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Having the "initial" explicitly in the rule name is not ideal. However, the path to the original contigs is somewhat variable and hence a discriminative rule name was chosen.
+# As piping is no longer supported by racon since v. 1.0.0 (https://github.com/isovic/racon/issues/46#issuecomment-519380742) and the file extensions are checked for all files including the contigs), some temporary files have to be created and removed afterwards.
+rule polish_long_read_contigs_w_racon_initial_round:
+    input: 
+        basecalls=rules.get_single_fastq_per_barcode.output,
+        contigs=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna"),
+        lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    params:
+        tmp_sam=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.tmp.sam"),
+        tmp_contigs=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/assembly_polished.tmp_draft.fa")
+    threads: config["racon"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/racon/round_{round}/lr/merged/{barcode}/assembly_polished.racon.log")
+    conda: "../envs/racon.yaml"
+    shell:
+        """
+        (date &&\
+        samtools view -h {input.lr_to_lr_contigs_mapping} > {params.tmp_sam} &&\
+        cp {input.contigs} {params.tmp_contigs} &&\
+        racon -t {threads} {input.basecalls} {params.tmp_sam} {params.tmp_contigs} > {output} &&\
+        rm {params.tmp_sam} {params.tmp_contigs} &&\
+        date) 2> >(tee {log}) > >(tee {log})
+        """
+
+#polished_contigs=os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/contigs.polished.fna")
+#ln -s contigs.polished.fasta {output.polished_contigs} &&\
+#rule assemble_long_and_short_reads_w_operams:
+#    input: 
+#        lr=rules.get_single_fastq_per_barcode.output,
+#        sr_r1=rules.run_fastp_on_short_reads.output.r1,
+#        sr_r2=rules.run_fastp_on_short_reads.output.r2,
+#        sr_contigs=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+#    output: 
+#        contigs=os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+#    params:
+#        tmp_lr=os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/tmp_lr_{barcode}.fna")
+#    threads: config["operams"]["threads"]
+#    log: os.path.join(RESULTS_DIR, "assembly/operams/lr_{barcode}-sr_{sr_sample}/opera_ms.log")
+#    shell:
+#        """
+#        (date &&\
+#        zcat {input.lr} > {params.tmp_lr} &&\
+#        {config[operams][bin]} --long-read-mapper minimap2 --num-processors {threads} --short-read1 {input.sr_r1} --short-read2 {input.sr_r2} --long-read {params.tmp_lr} --contig-file {input.sr_contigs} --out-dir $(dirname {output.contigs}) &&\
+#        ln -s contigs.fasta {output.contigs} &&\
+#        rm {params.tmp_lr} &&\
+#        date) 2> >(tee {log}) > >(tee {log})
+#        """
+
+#zcat {input.lr} > {params.tmp_lr} &&\
+rule assemble_long_and_short_reads_w_metaspades:
+    input: 
+        lr=rules.get_single_fastq_per_barcode.output,
+        sr_r1=rules.run_fastp_on_short_reads.output.r1,
+        sr_r2=rules.run_fastp_on_short_reads.output.r2,
+        sr_contigs=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+    output: 
+        contigs=os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    params:
+        tmp_lr=os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_{barcode}-sr_{sr_sample}/tmp_lr_{barcode}.fna")
+    threads: config["metaspades"]["threads"]
+    log: os.path.join(RESULTS_DIR, "assembly/metaspades_hybrid/lr_{barcode}-sr_{sr_sample}/metaspades.log")
+    conda: "../envs/spades.yaml"
+    shell:
+        """
+        (date &&\
+        spades.py --meta -o $(dirname {output.contigs}) -k 21,33,55,77 -t {threads} -1 {input.sr_r1} -2 {input.sr_r2} --nanopore {input.lr} &&\
+        ln -s contigs.fasta {output.contigs} &&\
+        date) &> >(tee {log})
+        """
+
+rule map_short_reads_to_hybrid_contigs_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs")
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_long_reads_to_hybrid_contigs_w_bwa_mem:
+    input: 
+        query_lr=rules.get_single_fastq_per_barcode.output,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr_{barcode}-sr_{sr_sample}/contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-lr_{barcode}_sr_{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -x ont2d -t {threads} {params.index_prefix} {input.query_lr} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# rule build_minimap2_long_read_contigs_ont_index:
+#     input: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}.fna")
+#     output: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}-ont.mmi")
+#     conda: "../envs/minimap2.yaml"
+#     shell:
+#         """
+#         date
+#         minimap2 -x map-ont -d {output} {input}
+#         date
+#         """
+
+rule map_long_reads_to_long_read_contigs:
+    input: 
+        query=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        index=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly-ont.mmi"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam")
+#    wildcard_constraints:
+#        barcode='barcode\d+'
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x_{barcode}_contigs")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.minimap2_samtools")
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_long_reads_to_long_read_polished_contigs:
+    input: 
+        query=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        index=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/{polisher}/round_{round}/lr/merged/{barcode}/assembly_polished-ont.mmi"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/{polisher}/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/polished/{polisher}/round_{round}/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_polished_contigs.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/polished/{polisher}/round_{round}/lr/merged/{barcode}/{barcode}_reads-x_{barcode}_polished_contigs")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/polished/{polisher}/round_{round}/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_polished_contigs.minimap2_samtools")
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# rule build_minimap2_long_read_contigs_sr_index:
+#     input: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}.fna")
+#     output: os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/{prefix}-sr.mmi")
+#     conda: "../envs/minimap2.yaml"
+#     shell:
+#         """
+#         date
+#         minimap2 -x sr -d {output} {input}
+#         date
+#         """
+
+rule map_short_reads_to_long_read_contigs_w_minimap2:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly-sr.mmi"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.bam")
+    params:
+        prefix=os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}_reads-x_{barcode}_contigs")
+    conda: "../envs/minimap2.yaml"
+    threads: config["minimap2"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/minimap2/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.minimap2_samtools")
+    shell:
+        """
+        (date &&\
+        minimap2 -a -t {threads} {input.index} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.prefix} > {output}
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Also map the short reads with bwa mem as it might be more sensitive (https://github.com/lh3/minimap2/issues/214#issuecomment-413492006)
+rule map_short_reads_to_long_read_contigs_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.sa"),
+        index_fa=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly.fna")
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/lr/merged/{barcode}/assembly"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x_{barcode}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{barcode}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule assemble_short_reads_w_megahit:
+    input: 
+        sr_r1=rules.run_fastp_on_short_reads.output.r1,
+        sr_r2=rules.run_fastp_on_short_reads.output.r2
+    output: os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.fna")
+    threads: config["megahit"]["threads"]
+    conda: "../envs/megahit.yaml"
+    params:
+        tmp_dir=os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}_tmp")
+    log: os.path.join(RESULTS_DIR, "assembly/megahit/{sr_sample}/final.contigs.megahit")
+    shell:
+        """
+        (date &&\
+        megahit -1 {input.sr_r1} -2 {input.sr_r2} -t {threads} -o {params.tmp_dir} &&\
+        mv {params.tmp_dir}/* $(dirname {output}) &&\
+        rmdir {params.tmp_dir} &&\
+        ln -s final.contigs.fa {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# N.B. A similar function exists above and code duplication is generally discouraged.
+rule map_short_reads_to_short_read_contigs_w_bwa_mem:
+    input: 
+        query_r1=rules.run_fastp_on_short_reads.output.r1,
+        query_r2=rules.run_fastp_on_short_reads.output.r2,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.sa"),
+        index_fa=rules.assemble_short_reads_w_megahit.output
+    output: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/sr/bwa_mem/{sr_sample}_reads-x-{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -t {threads} {params.index_prefix} {input.query_r1} {input.query_r2} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule map_long_reads_to_short_read_contigs_w_bwa_mem:
+    input: 
+        query_lr=rules.get_single_fastq_per_barcode.output,
+        index_amb=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.amb"),
+        index_ann=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.ann"),
+        index_bwt=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.bwt"),
+        index_pac=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.pac"),
+        index_sa=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs.sa"),
+        index_fa=rules.assemble_short_reads_w_megahit.output
+    output: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-{sr_sample}-{assembler}_contigs.bam")
+    params:
+        index_prefix=os.path.join(RESULTS_DIR, "assembly/{assembler}/{sr_sample}/final.contigs"),
+        out_prefix=os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-{sr_sample}-{assembler}_contigs")
+    conda: "../envs/bwa.yaml"
+    threads: config["bwa"]["threads"]
+    log: os.path.join(RESULTS_DIR, "mapping/lr/bwa_mem/{barcode}_reads-x-{sr_sample}-{assembler}_contigs.bwa_mem_samtools")
+    shell:
+        """
+        (date &&\
+        bwa mem -x ont2d -t {threads} {params.index_prefix} {input.query_lr} | samtools view -@ {config[samtools][view][threads]} -SbT {input.index_fa} | samtools sort -@ {config[samtools][sort][threads]} -m {config[samtools][sort][chunk_size]} -T {params.out_prefix} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+#npa, npc, .npl, and .npo
+rule run_nonpareil_on_short_reads:
+    input: os.path.join(RESULTS_DIR, "preprocessing/sr/{sr_sample}_R{orientation}_{suffix}.fastp.fq.gz"),
+    output: 
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npa"),
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npc"),
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npl"),
+        os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.npo")
+    params:
+        out_prefix=os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp"),
+        tmp_fasta=os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.tmp.fna")
+    threads: config["nonpareil"]["threads"]
+    log: os.path.join(RESULTS_DIR, "annotation/nonpareil/{sr_sample}_R{orientation}_{suffix}.fastp.nonpareil")
+    conda: "../envs/nonpareil.yaml"
+    shell:
+        """
+        (date &&\
+        zcat {input} | sed -n '1~4s/^@/>/p;2~4p' > {params.tmp_fasta} &&\
+        nonpareil -s {params.tmp_fasta} -R {config[nonpareil][memory]} -t {threads} -T kmer -f fasta -b {params.out_prefix} &&\
+        rm {params.tmp_fasta} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# TODO: Cleanup how the input files are defined. Currently explicitly specified, which is discouraged to do.
+rule create_nanopolish_index:
+    input:
+        multifast5_dir_01=os.path.join(DATA_DIR, "multifast5/20181108_0827_test"),
+        multifast5_dir_02=os.path.join(DATA_DIR, "multifast5/20181107_0906_same"),
+        multifast5_dir_03=os.path.join(DATA_DIR, "multifast5/20181106_1450_noselection_sizeselection"),
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz")
+    output:
+        index=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index"),
+        fai=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.fai"),
+        index_gzi=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.gzi"),
+        readdb=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.readdb")
+    log: os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.nanopolish_index")
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        (date &&\
+        nanopolish index -d {input.multifast5_dir_01}/ -d {input.multifast5_dir_02}/ -d {input.multifast5_dir_03}/ {input.reads} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule samtools_index:
+    input: "{prefix}.bam"
+    output: "{prefix}.bam.bai"
+    conda: "../envs/bwa.yaml"
+    shell:
+        """
+        date
+        samtools index {input} {output}
+        date
+        """
+
+# # Legacy version before snakemake's `checkpoint` was used
+# rule polish_long_read_contigs_w_nanopolish:
+#     input:
+#         lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+#         lr_to_lr_contigs_mapping_index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam.bai"),
+#         reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+#         draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+#         nanopolish_index=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index"),
+#         nanopolish_fai=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.fai"),
+#         nanopolish_index_gzi=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.gzi"),
+#         nanopolish_readdb=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz.index.readdb")
+#     output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+#     log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.nanopolish")
+#     threads: 4
+#     conda: "../envs/nanopolish.yaml"
+#     shell:
+#         """
+#         (date &&\
+#         nanopolish_makerange.py {input.draft_contigs} | parallel --results $(dirname {output}) -P 7 \
+#         nanopolish variants --consensus -o $(dirname {output})/polished.{{1}}.vcf -w {{1}} -r {input.reads} -b {input.lr_to_lr_contigs_mapping} -g {input.draft_contigs} -t {threads} --min-candidate-frequency 0.1 &&\
+#         nanopolish vcf2fasta -g {input.draft_contigs} $(dirname {output})/polished.*.vcf > {output} &&\
+# 
+#         date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+#         """
+
+# Define the individual windows that nanopolish will parallelize over
+checkpoint nanopolish_makerange:
+    input: os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna")
+    output: directory(os.path.join(RESULTS_DIR, "assembly/flye/polished/nanopolish/round_{round}/lr/merged/{barcode}/windows"))
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        mkdir {output}
+        nanopolish_makerange.py {input} > {output}/ranges.txt
+        for i in $(head -n 20 {output}/ranges.txt); do echo $i > {output}/"$i".coord; done
+        date
+        """
+
+# Perform variant calling on single window
+rule nanopolish_variant_call:
+    input: 
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+        lr_to_lr_contigs_mapping_index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam.bai"),
+        draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+        window=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/windows/{window}.coord")
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/vcfs/{window}.vcf")
+    threads: 4
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        nanopolish variants --consensus -o {output} -w {wildcards.window} -r {input.reads} -b {input.lr_to_lr_contigs_mapping} -g {input.draft_contigs} -t {threads} --min-candidate-frequency 0.1
+        """
+
+# Gather all the individual VCF files that were produced.
+# The "overview" comes from the rule that created the "split-up", i.e., from the `nanopolish_makerange` rule.
+# But the actual files that are aggregated here are those that should eventually be created and used, i.e., the VCF files
+# N.B. snakemake threw an error when the `window=` line below was wrong, but the error was supposedly coming from `nanopolish_makerange`. So much about tracking errors, which is not always straightforward as one might wish.
+def aggregate_nanopolish_vcf(wildcards):
+    checkpoint_output = checkpoints.nanopolish_makerange.get(**wildcards).output[0]
+    return expand(os.path.join(RESULTS_DIR, "assembly/flye/polished/nanopolish/round_{round}/lr/merged/{barcode}/vcfs/{window}.vcf"),
+            round=wildcards.round,
+            barcode=wildcards.barcode,
+            window=glob_wildcards(os.path.join(checkpoint_output, "{window}.coord")).window)
+
+# Combine the information in the individual VCF files to produce a polished assembly
+rule polish_long_read_contigs_w_nanopolish:
+    input: 
+        vcf_calls=aggregate_nanopolish_vcf,
+        draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+    output: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.fna")
+    log: os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/assembly_polished.nanopolish")
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        (date &&\
+        nanopolish vcf2fasta -g {input.draft_contigs} {input.vcf_calls} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+rule nanopolish_methylation_call:
+    input: 
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        lr_to_lr_contigs_mapping=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam"),
+        lr_to_lr_contigs_mapping_index=os.path.join(RESULTS_DIR, "mapping/lr/merged/{barcode}/{barcode}_reads-x-{barcode}-{assembler}_contigs.bam.bai"),
+        draft_contigs=os.path.join(RESULTS_DIR, "assembly/flye/lr/merged/{barcode}/assembly.fna"),
+        window=os.path.join(RESULTS_DIR, "assembly/{assembler}/polished/nanopolish/round_{round}/lr/merged/{barcode}/windows/{window}.coord")
+    output: os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_calls/{window}.methylation_calls.tsv")
+    threads: 4
+    conda: "../envs/nanopolish.yaml"
+    shell:
+        """
+        nanopolish call-methylation -t {threads} -r {input.reads} -b {input.lr_to_lr_contigs_mapping} -g {input.draft_contigs} -w {wildcards.window} > {output}
+        """
+
+def aggregate_nanopolish_methylation(wildcards):
+    checkpoint_output = checkpoints.nanopolish_makerange.get(**wildcards).output[0]
+    return expand(os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_calls/{window}.methylation_calls.tsv"),
+            assembler=wildcards.assembler,
+            barcode=wildcards.barcode,
+            round=wildcards.round,
+            window=glob_wildcards(os.path.join(checkpoint_output, "{window}.coord")).window)
+
+# Combine the information in the individual methylation call files to get the overall methylation state
+rule merge_methylation_calls_from_long_read_contigs_w_nanopolish:
+    input: aggregate_nanopolish_methylation,
+    output: os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_frequencies.tsv")
+    log: os.path.join(RESULTS_DIR, "annotation/methylation/{assembler}/nanopolish/round_{round}/lr/merged/{barcode}/methylation_frequencies")
+    shell:
+        """
+        (date &&\
+        ../scripts/calculate_methylation_frequency.py {input} > {output} &&\
+        date) 2> >(tee {log}.stderr) > >(tee {log}.stdout)
+        """
+
+# Works until the final shred and polish step, then crashes with a Key error
+# rule run_rebaler_on_short_read_contigs:
+#     input:
+#         lr=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+#         sr_contigs=rules.assemble_short_reads_w_megahit.output
+#     output: os.path.join(RESULTS_DIR, "assembly/rebaler/lr/merged/{barcode}/lr_{barcode}-sr_contigs_{sr_sample}.fna")
+#     threads: config["rebaler"]["threads"]
+#     params:
+#         tmp_sr_contigs=os.path.join(RESULTS_DIR, "assembly/rebaler/lr/merged/{barcode}/tmp_sr_contigs.fna")
+#     conda: "../envs/rebaler.yaml"
+#     shell:
+#         """
+#         date
+#         awk 'BEGIN{{counter=0}} FNR==NR{{if($1 ~ /^>/) {{a[$1]=counter; counter=counter+1}};next}}{{if($1 ~/^>/) print ">"a[$1]; else print}}' {input.sr_contigs} {input.sr_contigs} > {params.tmp_sr_contigs} # Rename the contigs from 0 to num. contigs to ensure rebaler has no "Key error"
+#         rebaler -t {threads} {params.tmp_sr_contigs} {input.lr} > {output}
+#         rm {params.tmp_sr_contigs}
+#         date
+#         """
+
+# For info on "--longs-reads" option see https://github.com/bbuchfink/diamond/releases/tag/v0.9.23
+rule run_diamond_on_long_reads:
+    input:
+        reads=os.path.join(RESULTS_DIR, "basecalled/merged/{barcode}.fastq.gz"),
+        db=config["diamond"]["db"]
+    output: os.path.join(RESULTS_DIR, "annotation/diamond/lr/merged/{barcode}.diamond.daa")
+    threads: config["diamond"]["threads"]
+    conda: "../envs/diamond.yaml"
+    shell:
+        """
+        date
+        diamond blastx -q {input.reads} --db {input.db} -p {threads} --long-reads --sensitive --outfmt 100 --out {output}
+        date
+        """
+
+rule call_genes_w_prodigal:
+    input: os.path.join(RESULTS_DIR, "assembly/{prefix}.fna")
+    output: os.path.join(RESULTS_DIR, "annotation/proteins/{prefix}.faa")
+    log: os.path.join(RESULTS_DIR, "annotation/proteins/{prefix}.prodigal.log")
+    conda: "../envs/prodigal.yaml"
+    shell: 
+        """
+        (date &&\
+        prodigal -a {output} -p meta -i {input} &&\
+        date) &> >(tee {log})
+        """
+
+# Write to DAA format to be able to use `diamond view` later to define custom formats.
+# N.B. `diamond view` seems not to support compressed input AND seems to require the *prefix* of the DAA file.
+# TODO: This could maybe be parametrized by the respectively used DB
+rule run_diamond_on_proteins:
+    input: 
+        proteins=rules.call_genes_w_prodigal.output,
+        db=config["diamond"]["db"]
+    output: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.daa")
+    params:
+        outfmt="100"
+    log: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.diamond.log")
+    conda: "../envs/diamond.yaml"
+    threads: config["diamond"]["threads"]
+    shell: 
+        """
+        (date &&\
+        diamond blastp -q {input.proteins} --db {input.db} -p {threads} --outfmt {params.outfmt} --out {output} &&\
+        date) &> >(tee {log})
+        """
+
+# Not sure why, but with the conda version of diamond (and `diamond view`) it was necessary to specify the entire column layout (outfmt) explicitly
+rule format_diamond_output_for_ideel:
+    input: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.daa")
+    output: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.tsv")
+    params:
+        outfmt="6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen"
+    log: os.path.join(RESULTS_DIR, "annotation/diamond/{prefix}.diamond_view.log")
+    conda: "../envs/diamond.yaml"
+    shell:
+        """
+        (date &&\
+        diamond view --daa $(echo {input} | sed 's/\.daa$//') --max-target-seqs 1 --outfmt {params.outfmt} --out {output} &&\
+        date) &> >(tee {log})
+        """