intial version

.gitignore

@@ -1,3 +1,9 @@
+# custom
+.snakemake/
+data/*/*
+dorado_models/*/
+
+
 # ---> Python
 # Byte-compiled / optimized / DLL files
 __pycache__/

README.md

@@ -1,2 +1,26 @@
 # eDNA_Stream_E01
 
+## Scope
+- Establish whether squiggle-based taxonomy is feasible on the small computational budget that is available
+
+## Dataset:
+- NO-MISS bacterial isolates, v10.4.1 chemistry on PromethION flow cells (https://epi2me.nanoporetech.com/nomiss_96bc_p2i_sup_2026/)
+    - Downsampling at the pod5 level intended for compute/time reasons -> the full dataset is 293 pod5 / 1.6 TB of raw data 
+
+## Dependencies:
+- dorado >= v 2.0.0
+- aws CLI
+- conda/mamba (mamba recommended)
+- python env with:
+    - python 3.14
+    - snakemake
+    - python-dotenv
+
+## Optional:
+- .env file with NCBI API key -> slightly speeds up reference fasta download
+
+## Running:
+```bash
+cd workflow
+snakemake --cores <cores> --resources gpu=<n_gpus>
+```

create_species_yaml.py

@@ -0,0 +1,13 @@
+"""
+Converts a txt file with species of interest and comments (1 per line) into the yaml format expected by the workflow
+"""
+LEADING_SPACE = '  '
+
+if __name__ == '__main__':
+    with open('reference/species_list.txt','rt') as ih, open('reference/species_list.yaml','wt') as oh:
+        oh.write('species_list:\n')
+        for line in ih:
+            if line.startswith('#'):
+                oh.write(LEADING_SPACE+line.strip()+'\n')
+                continue
+            oh.write(LEADING_SPACE+'- '+'_'.join(x.lower() for x in line.strip().split())+'\n')

reference/species_list.txt

@@ -0,0 +1,18 @@
+# Phylum Bacillota / Firmicutes (Represented)
+Bacillus subtilis
+Listeria monocytogenes
+Staphylococcus aureus
+# Phylum Bacillota / Firmicutes (False-Positive Checkpoint - Near Bacillus)
+Paenibacillus polymyxa
+# Phylum Pseudomonadota / Proteobacteria (Represented)
+Cronobacter sakazakii
+Citrobacter freundii
+Enterobacter cloacae
+Klebsiella pneumoniae
+Pseudomonas aeruginosa
+Salmonella enterica
+Shigella flexneri
+Vibrio cholerae
+# Phylum Pseudomonadota / Proteobacteria (False-Positive Checkpoints - Near Enterobacteriaceae & Vibrio)
+Escherichia coli K-12
+Aeromonas hydrophila

reference/species_list.yaml

@@ -0,0 +1,19 @@
+species_list:
+  # Phylum Bacillota / Firmicutes (Represented)
+  - bacillus_subtilis
+  - listeria_monocytogenes
+  - staphylococcus_aureus
+  # Phylum Bacillota / Firmicutes (False-Positive Checkpoint - Near Bacillus)
+  - paenibacillus_polymyxa
+  # Phylum Pseudomonadota / Proteobacteria (Represented)
+  - cronobacter_sakazakii
+  - citrobacter_freundii
+  - enterobacter_cloacae
+  - klebsiella_pneumoniae
+  - pseudomonas_aeruginosa
+  - salmonella_enterica
+  - shigella_flexneri
+  - vibrio_cholerae
+  # Phylum Pseudomonadota / Proteobacteria (False-Positive Checkpoints - Near Enterobacteriaceae & Vibrio)
+  - escherichia_coli_k-12
+  - aeromonas_hydrophila

workflow/.env_example

@@ -0,0 +1 @@
+NCBI_API_KEY=

workflow/Snakefile

@@ -0,0 +1,32 @@
+import os
+from dotenv import load_dotenv
+
+load_dotenv()
+
+configfile: "config/main.yaml"
+configfile: "config/dependencies.yaml"
+configfile: "../reference/species_list.yaml"
+
+module preparation:
+    snakefile: "rules/preparation.smk"
+    config: config
+
+module baseline:
+    snakefile: "rules/baseline_pipeline.smk"
+    config: config
+
+# Pulldown
+rule all:
+    input:
+        #"binaries/datasets",
+        #f"{config['dorado_model_dir']}/{config['fast_model_name']}",
+        #f"{config['dorado_model_dir']}/{config['hac_model_name']}",
+        #'../data/reference_genomes/full_reference.mmi',
+        #'../data/raw_pod5/'
+        #'../data/pod5_files_to_pull',
+        #'../data/basecalled_reads/hac.fastq.gz'
+        expand('../data/raw_pod5/PBK98658_853a956f_57f83f46_{batch}.fastq.gz',batch=range(1,config["pod5_dataset_size"]+1,config["pod5_stride"]))
+
+use rule * from preparation
+use rule * from baseline
+

workflow/config/dependencies.yaml

@@ -0,0 +1,3 @@
+arch: "amd64"
+
+datasets_binary : "binaries/datasets"

workflow/config/main.yaml

@@ -0,0 +1,16 @@
+tmp_dir: "/tmp"
+
+# dorado
+dorado_model_dir : "../dorado_models"
+hac_model_name: "dna_r10.4.1_e8.2_400bps_hac@v6.0.0"
+hac_benchmarks_file: "" # add path if available
+hac_min_q: 9 # empty string or 0 to disable
+fast_model_name: "dna_r10.4.1_e8.2_400bps_fast@v5.2.0"
+fast_benchmarks_file: "" #add path if available
+fast_min_q: 8 # empty string or 0 to disable
+
+
+# squiqqle dataset
+pod5_dataset: "nomiss_96BC_P2I_SUP_2026"
+pod5_dataset_size : 293
+pod5_stride : 12 # Use 1 in x pod5 files from the ONT NO-MISS dataset to reduce dataset size / compute requirements

workflow/envs/minimap.yaml

@@ -0,0 +1,4 @@
+channels:
+  - bioconda
+dependencies:
+  - minimap2

workflow/envs/pod5_manipulation.yaml

@@ -0,0 +1,7 @@
+channels:
+  - conda-forge
+dependencies:
+  - python=3.14
+  - pip
+  - pip:
+    - pod5

workflow/envs/samtools.yaml

@@ -0,0 +1,5 @@
+channels:
+  - bioconda
+dependencies:
+  - htslib
+  - samtools

workflow/rules/baseline_pipeline.smk

@@ -0,0 +1,68 @@
+def get_model_name(wildcards):
+    config_param=f'{wildcards.model}_model_name'
+    return config[config_param]
+
+def get_model_requirement(wildcards):
+    return f"{config['dorado_model_dir']}/{config[f'{wildcards.model}_model_name']}
+
+def get_benchmarking_file(wildcards):
+    config_param = f'{wildcards.model}_benchmarks_file'
+    if config[config_param]:
+        return f'--batchsize-benchmarks-file {config[config_param]}'
+    return ''
+
+def get_qscore(wildcards):
+    config_param = f'{wildcards.model}_min_q'
+    if config[config_param]:
+        return f'--min-qscore {config[config_param]}'
+    return ''
+
+def get_batch_names(input_file):
+    with open(input_file,'rt') as ih:
+        return [line.strip().split('/')[-1].split('.')[0] for line in ih]
+
+rule download_pod5:
+    output:
+        '../data/raw_pod5/PBK98658_853a956f_57f83f46_{batch}.pod5'
+    threads: 1
+    wildcard_constraints:
+        batch="\d+"
+    shell:
+        """
+        aws s3 cp --no-sign-request s3://ont-open-data/nomiss_96BC_P2I_SUP_2026/raw/pod5/PBK98658_853a956f_57f83f46_{wildcards.batch}.pod5 {output}
+        """
+
+rule basecall_pod5:
+    input:
+        pod5='../data/raw_pod5/PBK98658_853a956f_57f83f46_{batch}.pod5',
+        model=get_model_requirement
+    output:
+        '../data/basecalled_reads/{model}/PBK98658_853a956f_57f83f46_{batch}.fastq.gz'
+    threads:
+        32
+    resources:
+        gpu=1
+    params:
+        benchmarking=get_benchmarking_file,
+        min_qscore=get_qscore,
+        dorado_model=get_model_name
+    wildcard_constraints:
+        batch="\d+",
+        model="hac|fast"
+    shell:
+        """
+        dorado basecaller --models-directory {config[dorado_model_dir]} --emit-fastq {params.benchmarking} {params.min_qscore} {params.dorado_model} {input.pod5} > {output}
+        """
+
+rule concatenate_basecalled_fastq:
+    input:
+        expand('../data/basecalled_reads/{{model}}/PBK98658_853a956f_57f83f46_{batch}.fastq.gz',batch=range(1,config["pod5_dataset_size"]+1,config["pod5_stride"]))
+    output:
+        '../data/basecalled_reads/{model}.fastq.gz'
+    threads: 1
+    wildcard_constraints:
+        model="hac|fast"
+    shell:
+        """
+        zcat {input} > {output}
+        """

workflow/rules/preparation.smk

@@ -0,0 +1,84 @@
+import os
+import hashlib
+
+# Prepwork
+
+def get_genome_file_name(wildcards):
+    return '../data/reference_genomes/'+'_'.join(x.lower() for x in wildcards.species.split())+'.fasta'
+
+def get_unique_download_dir(wildcards):
+    unique_id = hashlib.md5(wildcards.species.encode()).hexdigest()[:8]
+    return config["tmp_dir"]+'/'+unique_id
+
+def format_cli_arg(wildcards):
+    # Replaces underscores with spaces for the CLI command
+    return wildcards.species.replace("_", " ")
+
+rule pull_ncbi_datasets_cli:
+    output:
+        config["datasets_binary"]
+    params:
+        arch=config["arch"]
+    threads: 1
+    shell:
+        """
+        curl https://ftp.ncbi.nlm.nih.gov/pub/datasets/command-line/v2/linux-{params.arch}/datasets -o {output}
+        chmod a+x {output}
+        """
+
+rule pull_dorado_models:
+    output:
+        directory(f"{config['dorado_model_dir']}/{config['fast_model_name']}"),
+        directory(f"{config['dorado_model_dir']}/{config['hac_model_name']}")
+    threads: 1
+    shell:
+        """
+        dorado download --models-directory {config[dorado_model_dir]} --model {config[fast_model_name]}
+        dorado download --models-directory {config[dorado_model_dir]} --model {config[hac_model_name]}
+        """
+
+rule download_reference_genome:
+    input:
+        config["datasets_binary"]
+    output:
+        temp('../data/reference_genomes/{species}.fasta')
+    params:
+        api_key=os.environ['NCBI_API_KEY'],
+        wd=get_unique_download_dir,
+        ncbi_tax_name=lambda w: format_cli_arg(w)
+    threads: 1
+    shell:
+        """
+        echo {wildcards.species}
+        rm -rf {params.wd}
+        mkdir {params.wd}
+        {config[datasets_binary]} download genome taxon "{params.ncbi_tax_name}" --reference --no-progressbar --filename {params.wd}/{wildcards.species}.zip
+        unzip -d {params.wd} -o {params.wd}/{wildcards.species}.zip 
+        mv {params.wd}/ncbi_dataset/data/*/*fna {output}
+        rm -rf {params.wd}
+        """
+
+rule concatenate_reference_genomes:
+    input:
+        expand("../data/reference_genomes/{species}.fasta",species=config["species_list"])
+    output:
+        "../data/reference_genomes/full_reference.fasta"
+    threads: 1
+    shell:
+        """
+        cat {input} > {output}
+        """
+
+rule create_minimap_index:
+    input: 
+        "../data/reference_genomes/full_reference.fasta"
+    output:
+        "../data/reference_genomes/full_reference.mmi"
+    threads:
+        32
+    conda:
+        "../envs/minimap.yaml"
+    shell:
+        """
+        minimap2 -x map-ont -t {threads} -d {output} {input}
+        """

workflow/scripts/download_reference_sequences.py

@@ -0,0 +1,11 @@
+from typing import TYPE_CHECKING
+from pathlib import Path
+if TYPE_CHECKING:
+    from snakemake.iocontainers import snakemake
+
+
+
+def download_species_fasta(species_name:str, ncbi_datasets_binary: str, tmp_dir : Path) -> Path:
+    """Download the reference genome fasta for a given species"""
+    output_fasta = tmp_dir / f'{species_name}.fasta'
+