From 08467102c88a7f2478cef243a9ead0ebfc2b043a Mon Sep 17 00:00:00 2001 From: Tom Kasper Date: Thu, 17 Sep 2026 22:16:42 +0100 Subject: [PATCH] intial version --- .gitignore | 6 ++ README.md | 24 ++++++ create_species_yaml.py | 13 +++ data/.keep | 0 dorado_models/.keep | 0 reference/dorado_batchsize_files/.keep | 0 reference/species_list.txt | 18 ++++ reference/species_list.yaml | 19 +++++ workflow/.env_example | 1 + workflow/Snakefile | 32 +++++++ workflow/binaries/.keep | 0 workflow/config/dependencies.yaml | 3 + workflow/config/main.yaml | 16 ++++ workflow/envs/minimap.yaml | 4 + workflow/envs/pod5_manipulation.yaml | 7 ++ workflow/envs/samtools.yaml | 5 ++ workflow/rules/baseline_pipeline.smk | 68 +++++++++++++++ workflow/rules/preparation.smk | 84 +++++++++++++++++++ .../scripts/download_reference_sequences.py | 11 +++ 19 files changed, 311 insertions(+) create mode 100644 create_species_yaml.py create mode 100644 data/.keep create mode 100644 dorado_models/.keep create mode 100644 reference/dorado_batchsize_files/.keep create mode 100644 reference/species_list.txt create mode 100644 reference/species_list.yaml create mode 100644 workflow/.env_example create mode 100644 workflow/Snakefile create mode 100644 workflow/binaries/.keep create mode 100644 workflow/config/dependencies.yaml create mode 100644 workflow/config/main.yaml create mode 100644 workflow/envs/minimap.yaml create mode 100644 workflow/envs/pod5_manipulation.yaml create mode 100644 workflow/envs/samtools.yaml create mode 100644 workflow/rules/baseline_pipeline.smk create mode 100644 workflow/rules/preparation.smk create mode 100644 workflow/scripts/download_reference_sequences.py diff --git a/.gitignore b/.gitignore index cac109b..5aa01a7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,9 @@ +# custom +.snakemake/ +data/*/* +dorado_models/*/ + + # ---> Python # Byte-compiled / optimized / DLL files __pycache__/ diff --git a/README.md b/README.md index 732bc79..06ef8ff 100644 --- a/README.md +++ b/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 --resources gpu= +``` \ No newline at end of file diff --git a/create_species_yaml.py b/create_species_yaml.py new file mode 100644 index 0000000..a820e34 --- /dev/null +++ b/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') diff --git a/data/.keep b/data/.keep new file mode 100644 index 0000000..e69de29 diff --git a/dorado_models/.keep b/dorado_models/.keep new file mode 100644 index 0000000..e69de29 diff --git a/reference/dorado_batchsize_files/.keep b/reference/dorado_batchsize_files/.keep new file mode 100644 index 0000000..e69de29 diff --git a/reference/species_list.txt b/reference/species_list.txt new file mode 100644 index 0000000..c886f95 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/reference/species_list.yaml b/reference/species_list.yaml new file mode 100644 index 0000000..4bb7dd0 --- /dev/null +++ b/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 diff --git a/workflow/.env_example b/workflow/.env_example new file mode 100644 index 0000000..e069bc4 --- /dev/null +++ b/workflow/.env_example @@ -0,0 +1 @@ +NCBI_API_KEY= \ No newline at end of file diff --git a/workflow/Snakefile b/workflow/Snakefile new file mode 100644 index 0000000..72ead89 --- /dev/null +++ b/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 + diff --git a/workflow/binaries/.keep b/workflow/binaries/.keep new file mode 100644 index 0000000..e69de29 diff --git a/workflow/config/dependencies.yaml b/workflow/config/dependencies.yaml new file mode 100644 index 0000000..8f9f26d --- /dev/null +++ b/workflow/config/dependencies.yaml @@ -0,0 +1,3 @@ +arch: "amd64" + +datasets_binary : "binaries/datasets" diff --git a/workflow/config/main.yaml b/workflow/config/main.yaml new file mode 100644 index 0000000..1ebb825 --- /dev/null +++ b/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 diff --git a/workflow/envs/minimap.yaml b/workflow/envs/minimap.yaml new file mode 100644 index 0000000..baa87a9 --- /dev/null +++ b/workflow/envs/minimap.yaml @@ -0,0 +1,4 @@ +channels: + - bioconda +dependencies: + - minimap2 \ No newline at end of file diff --git a/workflow/envs/pod5_manipulation.yaml b/workflow/envs/pod5_manipulation.yaml new file mode 100644 index 0000000..22fa4dd --- /dev/null +++ b/workflow/envs/pod5_manipulation.yaml @@ -0,0 +1,7 @@ +channels: + - conda-forge +dependencies: + - python=3.14 + - pip + - pip: + - pod5 \ No newline at end of file diff --git a/workflow/envs/samtools.yaml b/workflow/envs/samtools.yaml new file mode 100644 index 0000000..9d46417 --- /dev/null +++ b/workflow/envs/samtools.yaml @@ -0,0 +1,5 @@ +channels: + - bioconda +dependencies: + - htslib + - samtools \ No newline at end of file diff --git a/workflow/rules/baseline_pipeline.smk b/workflow/rules/baseline_pipeline.smk new file mode 100644 index 0000000..f46d762 --- /dev/null +++ b/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} + """ \ No newline at end of file diff --git a/workflow/rules/preparation.smk b/workflow/rules/preparation.smk new file mode 100644 index 0000000..803b989 --- /dev/null +++ b/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} + """ diff --git a/workflow/scripts/download_reference_sequences.py b/workflow/scripts/download_reference_sequences.py new file mode 100644 index 0000000..e7f26fe --- /dev/null +++ b/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' + \ No newline at end of file