From d939ccfdf987a28ab5728bd548657e1ec87114b5 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Wed, 4 Feb 2026 09:58:57 +0000 Subject: [PATCH 01/15] add gene identification with pprodigal, update CARD version --- config/config.yaml | 6 +++--- workflow/Snakefile | 24 +++++++++++++++--------- workflow/envs/gtdbtk.yaml | 1 + workflow/envs/pprodigal.yaml | 6 ++++++ workflow/rules/assembly.smk | 17 +++++++++++++++++ workflow/rules/classify.smk | 4 ++-- 6 files changed, 44 insertions(+), 14 deletions(-) create mode 100644 workflow/envs/pprodigal.yaml diff --git a/config/config.yaml b/config/config.yaml index a864f70..cbd8a33 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -4,7 +4,7 @@ pepfile: config/pep/config.yaml ## All results can be found under results/project-name/ project-name: "test" # -run-date: "2025-06-26" +run-date: "2026-02-03" data-handling: # path to store data within the workflow @@ -64,10 +64,10 @@ gtdb: use-local: True # if use-local is set to True, please specify the folder where the decompressed database is stored # this path is expected to lay under the data-handling resources folder - db-folder: gtdb/release226/ + db-folder: gtdb/release220/ card: - version: v4.0.0 + version: v4.0.1 dbfile: card.json url: https://card.mcmaster.ca/latest/data diff --git a/workflow/Snakefile b/workflow/Snakefile index f8e0c36..ea4d1ca 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -57,20 +57,15 @@ rule all: project=get_project(), sample=get_samples(), ), - # plasmid analysis + #proteins expand( - "results/{project}/output/plasmids/{sample}/{sample}_plasmid_summary.tsv", - project=get_project(), - sample=get_samples(), - ), - # resistance analysis - expand( - "results/{project}/output/ARGs/reads/{sample}/{sample}_read_ARGs.csv", + "results/{project}/output/proteins/{sample}/{sample}_proteins.faa", project=get_project(), sample=get_samples(), ), + # plasmid analysis expand( - "results/{project}/output/ARGs/assembly/{sample}/{sample}_assembly_ARGs.csv", + "results/{project}/output/plasmids/{sample}/{sample}_plasmid_summary.tsv", project=get_project(), sample=get_samples(), ), @@ -99,6 +94,17 @@ rule all: project=get_project(), sample=get_samples(), ), + # resistance analysis + expand( + "results/{project}/output/ARGs/reads/{sample}/{sample}_read_ARGs.csv", + project=get_project(), + sample=get_samples(), + ), + expand( + "results/{project}/output/ARGs/assembly/{sample}/{sample}_assembly_ARGs.csv", + project=get_project(), + sample=get_samples(), + ), """ diff --git a/workflow/envs/gtdbtk.yaml b/workflow/envs/gtdbtk.yaml index a84bf84..3144b50 100644 --- a/workflow/envs/gtdbtk.yaml +++ b/workflow/envs/gtdbtk.yaml @@ -4,3 +4,4 @@ channels: - nodefaults dependencies: - gtdbtk = 2.4.1 + - numpy = 1.23.1 diff --git a/workflow/envs/pprodigal.yaml b/workflow/envs/pprodigal.yaml new file mode 100644 index 0000000..dad1fe8 --- /dev/null +++ b/workflow/envs/pprodigal.yaml @@ -0,0 +1,6 @@ +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - pprodigal = 1.0.1 diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index 3c9b10a..aba95e9 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -145,3 +145,20 @@ rule cleanup_megahit_output: touch("results/{project}/megahit/{sample}_cleanup.done"), log: "logs/{project}/assembly/{sample}_cleanup.log", + + +rule protein_identification: + input: + contigs=get_assembly, + output: + faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa", + fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna", + gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff", + log: + "logs/{project}/prodigal/{sample}_prodigal_run.log", + threads: 32 + conda: + "../envs/pprodigal.yaml" + shell: + "pprodigal -i {input.contigs} -o {output.gff} -a {output.faa} " + "-d {output.fna} -p meta --tasks {threads} > {log} 2>&1" diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 8c29097..5549053 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -107,14 +107,14 @@ if config["gtdb"]["use-local"]: db_folder=get_gtdb_folder(), clf_outdir=lambda wildcards, output: Path(output.json).parent, json=lambda wildcards, output: Path(output.json).name, - threads: 64 + threads: 32 log: "logs/{project}/gtdbtk/{sample}_classify.log", conda: "../envs/gtdbtk.yaml" shell: "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--mash_db {params.db_folder}mash_db/ --cpus 40 " + "--mash_db {params.db_folder}mash_db/ --cpus {threads} --pplacer_cpus 30 " "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" From 8ecd88610e9172e33d4d8d2e3699c5e66a5d3f60 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Thu, 5 Feb 2026 08:33:04 +0000 Subject: [PATCH 02/15] add gzip of protein files --- workflow/Snakefile | 2 +- workflow/rules/assembly.smk | 24 +++++++++++++++++++++--- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index ea4d1ca..74bdf57 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -59,7 +59,7 @@ rule all: ), #proteins expand( - "results/{project}/output/proteins/{sample}/{sample}_proteins.faa", + "results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", project=get_project(), sample=get_samples(), ), diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index aba95e9..3aea577 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -79,13 +79,13 @@ rule gzip_assembly: contigs=get_assembly, output: "results/{project}/output/fastas/{sample}/{sample}.fa.gz", - threads: 64 + threads: 20 log: "logs/{project}/assembly/{sample}_gzip.log", conda: "../envs/unix.yaml" shell: - "pigz -c {input.contigs} > {output} 2> {log}" + "gzip -c {input.contigs} > {output} 2> {log}" rule assembly_summary: @@ -155,10 +155,28 @@ rule protein_identification: fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna", gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff", log: - "logs/{project}/prodigal/{sample}_prodigal_run.log", + "logs/{project}/proteins/{sample}.log", threads: 32 conda: "../envs/pprodigal.yaml" shell: "pprodigal -i {input.contigs} -o {output.gff} -a {output.faa} " "-d {output.fna} -p meta --tasks {threads} > {log} 2>&1" + + +rule gzip_proteins: + input: + faa=rules.protein_identification.output.faa, + fna=rules.protein_identification.output.fna, + gff=rules.protein_identification.output.gff, + output: + faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", + fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna.gz", + gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff.gz", + threads: 20 + log: + "logs/{project}/proteins/{sample}_gzip.log", + conda: + "../envs/unix.yaml" + shell: + "gzip {input.faa} {input.fna} {input.gff} > {log} 2>&1" From 3b9684b6e6ac64de7a90e9742e54df789fc574b5 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Fri, 6 Feb 2026 13:21:06 +0000 Subject: [PATCH 03/15] update tool and database versions --- config/config.yaml | 2 +- config/multiqc_config.yaml | 5 ++- workflow/Snakefile | 13 ++---- workflow/envs/card.yaml | 3 +- workflow/envs/checkm2.yaml | 2 +- workflow/envs/coverm.yaml | 2 +- workflow/envs/das_tool.yaml | 2 +- workflow/envs/fastqc.yaml | 6 +++ workflow/envs/genomad.yaml | 4 +- workflow/envs/gtdbtk.yaml | 3 +- workflow/envs/metabat.yaml | 2 +- workflow/envs/metacoag.yaml | 2 +- workflow/envs/minimap2.yaml | 6 +-- workflow/envs/python.yaml | 10 ++--- workflow/envs/rbt.yaml | 2 +- workflow/envs/unix.yaml | 2 +- workflow/rules/analysis.smk | 62 +++++++++++++++++++++------ workflow/rules/assembly.smk | 37 +--------------- workflow/rules/classify.smk | 6 +-- workflow/rules/common.smk | 8 ++-- workflow/rules/host_filtering.smk | 14 +++--- workflow/rules/qc.smk | 71 +++++++++++++++++++++++++++++-- workflow/rules/trimming.smk | 48 --------------------- 23 files changed, 163 insertions(+), 149 deletions(-) create mode 100644 workflow/envs/fastqc.yaml delete mode 100644 workflow/rules/trimming.smk diff --git a/config/config.yaml b/config/config.yaml index cbd8a33..cfcc4cb 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -64,7 +64,7 @@ gtdb: use-local: True # if use-local is set to True, please specify the folder where the decompressed database is stored # this path is expected to lay under the data-handling resources folder - db-folder: gtdb/release220/ + db-folder: gtdb/release226/ card: version: v4.0.1 diff --git a/config/multiqc_config.yaml b/config/multiqc_config.yaml index 24f4723..005aa10 100644 --- a/config/multiqc_config.yaml +++ b/config/multiqc_config.yaml @@ -31,7 +31,8 @@ fn_clean_exts: - "_L002_R2_001" - "_S" - ".1" - + - ".2" + # customising general Statistics table_columns_visible: Reads Quality Control: @@ -46,4 +47,4 @@ table_columns_visible: after_filtering_q30_bases: False after_filtering_gc_content: False pct_surviving: True - pct_adapter: False \ No newline at end of file + pct_adapter: False diff --git a/workflow/Snakefile b/workflow/Snakefile index 74bdf57..2497d86 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -20,7 +20,6 @@ report: "report/workflow.rst" include: "rules/common.smk" include: "rules/qc.smk" include: "rules/host_filtering.smk" -include: "rules/trimming.smk" include: "rules/assembly.smk" include: "rules/metabat.smk" include: "rules/metacoag.smk" @@ -57,12 +56,6 @@ rule all: project=get_project(), sample=get_samples(), ), - #proteins - expand( - "results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", - project=get_project(), - sample=get_samples(), - ), # plasmid analysis expand( "results/{project}/output/plasmids/{sample}/{sample}_plasmid_summary.tsv", @@ -96,12 +89,12 @@ rule all: ), # resistance analysis expand( - "results/{project}/output/ARGs/reads/{sample}/{sample}_read_ARGs.csv", + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", project=get_project(), sample=get_samples(), ), expand( - "results/{project}/output/ARGs/assembly/{sample}/{sample}_assembly_ARGs.csv", + "results/{project}/output/resistance/CARD/assembly/{sample}/{sample}_assembly_ARGs.csv", project=get_project(), sample=get_samples(), ), @@ -109,7 +102,7 @@ rule all: """ expand( - "results/{project}/output/ARGs/mags/{sample}/all_mags.done", + "results/{project}/output/resistance/CARD/mags/{sample}/all_mags.done", project=get_project(), sample=get_samples(), ), diff --git a/workflow/envs/card.yaml b/workflow/envs/card.yaml index 88f740f..1c69e2a 100644 --- a/workflow/envs/card.yaml +++ b/workflow/envs/card.yaml @@ -3,5 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - rgi = 6.0.3 - - kma = 1.4.9 \ No newline at end of file + - rgi = 6.0.5 diff --git a/workflow/envs/checkm2.yaml b/workflow/envs/checkm2.yaml index a9d75be..40d8f09 100644 --- a/workflow/envs/checkm2.yaml +++ b/workflow/envs/checkm2.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - checkm2 = 1.0.1 \ No newline at end of file + - checkm2 = 1.1.0 diff --git a/workflow/envs/coverm.yaml b/workflow/envs/coverm.yaml index 4ef59df..3fc7948 100644 --- a/workflow/envs/coverm.yaml +++ b/workflow/envs/coverm.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - coverm = 0.6.1 \ No newline at end of file + - coverm = 0.7.0 diff --git a/workflow/envs/das_tool.yaml b/workflow/envs/das_tool.yaml index bcb2764..3c32a4e 100644 --- a/workflow/envs/das_tool.yaml +++ b/workflow/envs/das_tool.yaml @@ -2,4 +2,4 @@ channels: - conda-forge - bioconda dependencies: - - das_tool=1.1.6 + - das_tool=1.1.7 diff --git a/workflow/envs/fastqc.yaml b/workflow/envs/fastqc.yaml new file mode 100644 index 0000000..2e917d0 --- /dev/null +++ b/workflow/envs/fastqc.yaml @@ -0,0 +1,6 @@ +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - fastqc = 0.12.1 diff --git a/workflow/envs/genomad.yaml b/workflow/envs/genomad.yaml index d6589cf..3424746 100644 --- a/workflow/envs/genomad.yaml +++ b/workflow/envs/genomad.yaml @@ -1,6 +1,6 @@ -channels: +channels: - conda-forge - bioconda - nodefaults dependencies: - - genomad = 1.8.0 + - genomad = 1.11.2 diff --git a/workflow/envs/gtdbtk.yaml b/workflow/envs/gtdbtk.yaml index 3144b50..d95f83b 100644 --- a/workflow/envs/gtdbtk.yaml +++ b/workflow/envs/gtdbtk.yaml @@ -3,5 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - gtdbtk = 2.4.1 - - numpy = 1.23.1 + - gtdbtk = 2.6.1 diff --git a/workflow/envs/metabat.yaml b/workflow/envs/metabat.yaml index 47d94d6..a31cac3 100644 --- a/workflow/envs/metabat.yaml +++ b/workflow/envs/metabat.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - metabat2 = 2.15 \ No newline at end of file + - metabat2 = 2.18 diff --git a/workflow/envs/metacoag.yaml b/workflow/envs/metacoag.yaml index 6a3e0d4..e8c0241 100644 --- a/workflow/envs/metacoag.yaml +++ b/workflow/envs/metacoag.yaml @@ -3,5 +3,5 @@ channels: - bioconda - nodefaults dependencies: - - metacoag = 1.2.1 + - metacoag = 1.2.2 - biopython = 1.81 diff --git a/workflow/envs/minimap2.yaml b/workflow/envs/minimap2.yaml index 30b228a..f625c9c 100644 --- a/workflow/envs/minimap2.yaml +++ b/workflow/envs/minimap2.yaml @@ -3,8 +3,8 @@ channels: - bioconda - nodefaults dependencies: - - minimap2 = 2.24 - - samtools = 1.16.1 + - minimap2 = 2.30 + - samtools = 1.23 - pip - pip: - - Command \ No newline at end of file + - Command diff --git a/workflow/envs/python.yaml b/workflow/envs/python.yaml index bf746a3..1682687 100644 --- a/workflow/envs/python.yaml +++ b/workflow/envs/python.yaml @@ -3,9 +3,9 @@ channels: - bioconda - anaconda dependencies: - - altair = 5.2.0 - - biopython = 1.82 + - altair = 6.0.0 + - biopython = 1.85 - distinctipy = 1.3.4 - - numpy = 1.26.3 - - matplotlib = 3.8.2 - - pandas = 2.1.4 + - numpy = 2.4.1 + - matplotlib = 3.10.8 + - pandas = 2.3.3 diff --git a/workflow/envs/rbt.yaml b/workflow/envs/rbt.yaml index c11823a..ebc035f 100644 --- a/workflow/envs/rbt.yaml +++ b/workflow/envs/rbt.yaml @@ -3,4 +3,4 @@ channels: - bioconda - nodefaults dependencies: - - rust-bio-tools = 0.42.0 \ No newline at end of file + - rust-bio-tools = 0.42.2 diff --git a/workflow/envs/unix.yaml b/workflow/envs/unix.yaml index f148464..fc89a8a 100644 --- a/workflow/envs/unix.yaml +++ b/workflow/envs/unix.yaml @@ -6,4 +6,4 @@ dependencies: - git - sed - curl - - pigz \ No newline at end of file + - gzip diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 1f59532..b8fcec2 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -1,3 +1,39 @@ +#gene identification +rule gene_identification: + input: + contigs=get_assembly, + output: + faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa", + fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna", + gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff", + log: + "logs/{project}/proteins/{sample}.log", + threads: 32 + conda: + "../envs/pprodigal.yaml" + shell: + "pprodigal -i {input.contigs} -o {output.gff} -a {output.faa} " + "-d {output.fna} -p meta --tasks {threads} > {log} 2>&1" + + +rule gzip_proteins: + input: + faa=rules.gene_identification.output.faa, + fna=rules.gene_identification.output.fna, + gff=rules.gene_identification.output.gff, + output: + faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", + fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna.gz", + gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff.gz", + threads: 20 + log: + "logs/{project}/proteins/{sample}_gzip.log", + conda: + "../envs/unix.yaml" + shell: + "gzip {input.faa} {input.fna} {input.gff} > {log} 2>&1" + + # Plasmid analysis rule load_genomad_DB: output: @@ -50,11 +86,11 @@ rule move_genomad_output: sum_folder=lambda wildcards, input: Path(input.plasmid_tsv).parent, log: "logs/{project}/plasmids/{sample}_move_output.log", - threads: 64 + threads: 32 conda: "../envs/unix.yaml" shell: - "scp {params.sum_folder}/* {params.outdir}/ > {log} 2>&1" + "scp {params.sum_folder}/*.tsv {params.sum_folder}/*.json {params.outdir}/ > {log} 2>&1" # Resistance analysis @@ -115,9 +151,9 @@ rule CARD_read_run: fa=get_filtered_gz_fastqs, db=rules.CARD_annotation.output, output: - txt="results/{project}/output/ARGs/reads/{sample}/{sample}.gene_mapping_data.txt", + txt="results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", bam=temp( - "results/{project}/output/ARGs/reads/{sample}/{sample}.sorted.length_100.bam" + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.sorted.length_100.bam" ), params: folder=lambda wildcards, output: Path(output.txt).parent, @@ -136,7 +172,7 @@ rule CARD_read_sample_summary: input: txt=rules.CARD_read_run.output.txt, output: - csv="results/{project}/output/ARGs/reads/{sample}/{sample}_read_ARGs.csv", + csv="results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", params: case="reads", log: @@ -154,7 +190,7 @@ rule CARD_load_DB: input: db=get_card_db_file(), read_args=expand( - "results/{project}/output/ARGs/reads/{sample}/{sample}.gene_mapping_data.txt", + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", sample=get_samples(), project=get_project(), ), @@ -171,11 +207,11 @@ rule CARD_load_DB: rule CARD_assembly_run: input: - fa=rules.gzip_assembly.output, + faa=rules.gzip_proteins.output.faa, db=rules.CARD_load_DB.output, output: - txt="results/{project}/output/ARGs/assembly/{sample}/{sample}.txt", - json="results/{project}/output/ARGs/assembly/{sample}/{sample}.json", + txt="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.txt", + json="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.json", params: path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), log: @@ -184,8 +220,8 @@ rule CARD_assembly_run: conda: "../envs/card.yaml" shell: - "rgi main -i {input.fa} -o {params.path_wo_ext} " - "-t contig -a DIAMOND --low_quality --local " + "rgi main -i {input.faa} -o {params.path_wo_ext} " + "-t protein -a DIAMOND --local " "-n {threads} --clean > {log} 2>&1" @@ -194,7 +230,7 @@ use rule CARD_read_sample_summary as CARD_assembly_sample_summary with: txt=rules.CARD_assembly_run.output.txt, json=rules.CARD_assembly_run.output.json, output: - csv="results/{project}/output/ARGs/assembly/{sample}/{sample}_assembly_ARGs.csv", + csv="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}_assembly_ARGs.csv", params: case="assembly", log: @@ -206,7 +242,7 @@ use rule CARD_assembly_run as CARD_mag_run with: fa="results/{project}/output/fastas/{sample}/mags/{binID}.fa.gz", db=rules.CARD_load_DB.output, output: - txt="results/{project}/output/ARGs/mags/{sample}/{binID}.txt", + txt="results/{project}/output/resistance/CARD/mags/{sample}/{binID}.txt", params: path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), log: diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index 3aea577..5ea7756 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -144,39 +144,4 @@ rule cleanup_megahit_output: output: touch("results/{project}/megahit/{sample}_cleanup.done"), log: - "logs/{project}/assembly/{sample}_cleanup.log", - - -rule protein_identification: - input: - contigs=get_assembly, - output: - faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa", - fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna", - gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff", - log: - "logs/{project}/proteins/{sample}.log", - threads: 32 - conda: - "../envs/pprodigal.yaml" - shell: - "pprodigal -i {input.contigs} -o {output.gff} -a {output.faa} " - "-d {output.fna} -p meta --tasks {threads} > {log} 2>&1" - - -rule gzip_proteins: - input: - faa=rules.protein_identification.output.faa, - fna=rules.protein_identification.output.fna, - gff=rules.protein_identification.output.gff, - output: - faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", - fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna.gz", - gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff.gz", - threads: 20 - log: - "logs/{project}/proteins/{sample}_gzip.log", - conda: - "../envs/unix.yaml" - shell: - "gzip {input.faa} {input.fna} {input.gff} > {log} 2>&1" + "logs/{project}/assembly/{sample}_cleanup.log", \ No newline at end of file diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 5549053..34b6932 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -104,7 +104,6 @@ if config["gtdb"]["use-local"]: ) ), params: - db_folder=get_gtdb_folder(), clf_outdir=lambda wildcards, output: Path(output.json).parent, json=lambda wildcards, output: Path(output.json).name, threads: 32 @@ -114,7 +113,7 @@ if config["gtdb"]["use-local"]: "../envs/gtdbtk.yaml" shell: "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--mash_db {params.db_folder}mash_db/ --cpus {threads} --pplacer_cpus 30 " + "--cpus {threads} --pplacer_cpus 30 " "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" @@ -152,7 +151,8 @@ else: "../envs/gtdbtk.yaml" shell: "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--cpus {threads} --genome_dir {input.bins}/ --out_dir {output.outdir}/ && " + "--cpus {threads} --pplacer_cpus 30 " + "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index abed525..4019493 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -91,15 +91,15 @@ def get_checkm2_db(): def get_filtered_fastqs(wildcards): return [ - "results/{project}/filtered/fastqs/{sample}_R1.fastq", - "results/{project}/filtered/fastqs/{sample}_R2.fastq", + "results/{project}/filtered/{sample}_R1.fastq", + "results/{project}/filtered/{sample}_R2.fastq", ] def get_filtered_gz_fastqs(wildcards): return [ - "results/{project}/filtered/fastqs/{sample}_R1.fastq.gz", - "results/{project}/filtered/fastqs/{sample}_R2.fastq.gz", + "results/{project}/filtered/{sample}_R1.fastq.gz", + "results/{project}/filtered/{sample}_R2.fastq.gz", ] diff --git a/workflow/rules/host_filtering.smk b/workflow/rules/host_filtering.smk index 238d47f..d3a0a9e 100644 --- a/workflow/rules/host_filtering.smk +++ b/workflow/rules/host_filtering.smk @@ -79,7 +79,7 @@ rule filter_human: output: filtered=temp( expand( - "results/{{project}}/filtered/fastqs/{{sample}}_{read}.fastq", + "results/{{project}}/filtered/{{sample}}_{read}.fastq", read=["R1", "R2"], ) ), @@ -97,16 +97,16 @@ rule filter_human: rule gzip_filtered_reads: input: - "results/{project}/filtered/fastqs/{sample}_{read}.fastq", + "results/{project}/filtered/{sample}_{read}.fastq", output: - "results/{project}/filtered/fastqs/{sample}_{read}.fastq.gz", + "results/{project}/filtered/{sample}_{read}.fastq.gz", log: "logs/{project}/human_filtering/gzip_{sample}_{read}.log", - threads: 64 + threads: 20 conda: "../envs/unix.yaml" shell: - "pigz -k {input} > {log} 2>&1" + "gzip -k {input} > {log} 2>&1" if config["host-filtering"]["do-host-filtering"]: @@ -152,9 +152,9 @@ if config["host-filtering"]["do-host-filtering"]: "../envs/minimap2.yaml" shell: "(samtools fastq -F 3584 -f 77 {input.bam} | " - "pigz -c > {output.filtered[0]} && " + "gzip -c > {output.filtered[0]} && " "samtools fastq -F 3584 -f 141 {input.bam} | " - "pigz -c > {output.filtered[1]}) > {log} 2>&1" + "gzip -c > {output.filtered[1]}) > {log} 2>&1" ## TODO diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index b00345e..337bc45 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -1,24 +1,87 @@ ## read QC +# fastp in paired-end mode for Illumina paired-end data +# version in this wrapper: fastp=1.0.1 +rule fastp: + input: + sample=get_local_fastqs, + output: + trimmed=temp( + [ + "results/{project}/trimmed/fastp/{sample}.1.fastq.gz", + "results/{project}/trimmed/fastp/{sample}.2.fastq.gz", + ] + ), + html=temp("results/{project}/trimmed/fastp/{sample}.html"), + json="results/{project}/report_prerequisites/qc/{sample}.fastp.json", + params: + adapters=get_adapters, + extra="--qualified_quality_phred {phred} --length_required {minlen}".format( + phred=(config["quality-criteria"]["min-PHRED"]), + minlen=(config["quality-criteria"]["min-length-reads"]), + ), + log: + "logs/{project}/fastp/{sample}.log", + threads: 16 + wrapper: + "v7.1.0/bio/fastp" + + +"""# version in this wrapper: fastqc=0.12.1 rule fastqc: input: - get_trimmed_fastqs, + rules.fastp.output.trimmed, + #get_trimmed_fastqs, output: html=temp("results/{project}/qc/fastqc/{sample}_trimmed.html"), zip=temp("results/{project}/qc/fastqc/{sample}_trimmed_fastqc.zip"), + threads: 4 + resources: + mem_mb=1024, log: "logs/{project}/fastqc/{sample}.log", wrapper: - "v1.23.5/bio/fastqc" + "v7.6.0/bio/fastqc" +""" + + +rule fastqc: + input: + rules.fastp.output.trimmed, + output: + zip=temp( + expand( + "results/{{project}}/trimmed/fastp/{{sample}}.{read}_fastqc.zip", + read=["1", "2"], + ) + ), + html=temp( + expand( + "results/{{project}}/trimmed/fastp/{{sample}}.{read}_fastqc.html", + read=["1", "2"], + ) + ), + threads: 4 + resources: + mem_mb=1024, + log: + "logs/{project}/fastqc/{sample}.log", + conda: + "../envs/fastqc.yaml" + shell: + "fastqc --memory {resources.mem_mb} --threads {threads} " + "--format fastq --quiet {input} > {log} 2>&1" +# version in this wrapper: multiqc=1.33 rule multiqc: input: expand( [ - "results/{{project}}/qc/fastqc/{sample}_trimmed_fastqc.zip", + "results/{{project}}/trimmed/fastp/{sample}.{read}_fastqc.zip", "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", ], sample=get_samples(), + read=["1", "2"], ), output: report( @@ -38,7 +101,7 @@ rule multiqc: log: "logs/{project}/multiqc.log", wrapper: - "v3.3.1/bio/multiqc" + "v8.1.1/bio/multiqc" rule qc_summary: diff --git a/workflow/rules/trimming.smk b/workflow/rules/trimming.smk deleted file mode 100644 index 2cf0384..0000000 --- a/workflow/rules/trimming.smk +++ /dev/null @@ -1,48 +0,0 @@ -RAW_DATA_PATH = get_data_path() - -""" -# copy files to local -rule copy_fastq: - output: - raw1=f"{RAW_DATA_PATH}{{project}}/{{sample}}_R1.fastq.gz", - raw2=f"{RAW_DATA_PATH}{{project}}/{{sample}}_R2.fastq.gz", - params: - outdir=lambda wildcards, output: Path(output.raw1).parent, - # returns a list of folder and the filenames for R1 and R2 reads - fastqs=get_fastqs, - threads: 20 - log: - "logs/{project}/copy_data/{sample}.log", - conda: - "../envs/unix.yaml" - shell: - "(mkdir -p {params.outdir} && " - "tar cpfz - -C {params.fastqs[0]}/ {params.fastqs[1]} {params.fastqs[2]} | " - "(cd {params.outdir} ; tar xpfz -)) > {log} 2>&1" -""" - - -# fastp in paired-end mode for Illumina paired-end data -rule fastp: - input: - sample=get_local_fastqs, - output: - trimmed=temp( - [ - "results/{project}/trimmed/fastp/{sample}.1.fastq.gz", - "results/{project}/trimmed/fastp/{sample}.2.fastq.gz", - ] - ), - html=temp("results/{project}/trimmed/fastp/{sample}.html"), - json="results/{project}/report_prerequisites/qc/{sample}.fastp.json", - params: - adapters=get_adapters, - extra="--qualified_quality_phred {phred} --length_required {minlen}".format( - phred=(config["quality-criteria"]["min-PHRED"]), - minlen=(config["quality-criteria"]["min-length-reads"]), - ), - log: - "logs/{project}/fastp/{sample}.log", - threads: 16 - wrapper: - "v2.6.0/bio/fastp" From 213f2d18970dddbf695def0338d4d230fc12a9e2 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Tue, 10 Feb 2026 16:50:28 +0000 Subject: [PATCH 04/15] fix gtdbtk, fix and restructure report files --- workflow/envs/gtdbtk.yaml | 2 + workflow/rules/analysis.smk | 5 +- workflow/rules/assembly.smk | 38 +++++---- workflow/rules/bin_qc.smk | 107 ++++++++----------------- workflow/rules/classify.smk | 6 +- workflow/rules/common.smk | 10 +-- workflow/rules/das_tool.smk | 4 +- workflow/rules/host_filtering.smk | 14 ++-- workflow/rules/qc.smk | 37 +++++---- workflow/rules/report.smk | 10 +-- workflow/scripts/assembly_summary.py | 19 ++++- workflow/scripts/bin_summary_sample.py | 25 ++---- workflow/scripts/quality_summary.py | 18 ++--- 13 files changed, 130 insertions(+), 165 deletions(-) diff --git a/workflow/envs/gtdbtk.yaml b/workflow/envs/gtdbtk.yaml index d95f83b..46b47f4 100644 --- a/workflow/envs/gtdbtk.yaml +++ b/workflow/envs/gtdbtk.yaml @@ -4,3 +4,5 @@ channels: - nodefaults dependencies: - gtdbtk = 2.6.1 + # fix because gtdbtk doesn't work with python 3.14 + - python = 3.13 diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index b8fcec2..56f90f6 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -1,4 +1,4 @@ -#gene identification +# gene identification rule gene_identification: input: contigs=get_assembly, @@ -155,6 +155,9 @@ rule CARD_read_run: bam=temp( "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.sorted.length_100.bam" ), + bai=temp( + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.sorted.length_100.bam.bai" + ), params: folder=lambda wildcards, output: Path(output.txt).parent, log: diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index 5ea7756..d4019c2 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -7,7 +7,7 @@ rule megahit: output: contigs=temp("results/{project}/megahit/{sample}/final.contigs.fa"), outdir=temp(directory("results/{project}/megahit/{sample}/")), - log="results/{project}/report_prerequisites/assembly/{sample}_megahit.log", + log="results/{project}/output/report/prerequisites/assembly/{sample}_megahit.log", done=touch("results/{project}/megahit/{sample}.done"), params: threshold=get_contig_length_threshold(), @@ -29,7 +29,7 @@ rule map_to_assembly: fastqs=get_filtered_gz_fastqs, output: bam=temp( - "results/{project}/report_prerequisites/assembly/{sample}_reads_mapped.bam" + "results/{project}/output/report/prerequisites/assembly/{sample}_reads_mapped.bam" ), threads: 64 log: @@ -47,7 +47,7 @@ rule index_assembly_alignment: rules.map_to_assembly.output.bam, output: bai=temp( - "results/{project}/report_prerequisites/assembly/{sample}_reads_mapped.bam.bai" + "results/{project}/output/report/prerequisites/assembly/{sample}_reads_mapped.bam.bai" ), threads: 20 log: @@ -63,7 +63,7 @@ rule reads_mapped_assembly: bam=rules.map_to_assembly.output.bam, bai=rules.index_assembly_alignment.output.bai, output: - bai="results/{project}/report_prerequisites/assembly/{sample}_reads_mapped.txt", + bai="results/{project}/output/report/prerequisites/assembly/{sample}_reads_mapped.txt", threads: 20 log: "logs/{project}/assembly/{sample}_mapping_reads.log", @@ -92,11 +92,19 @@ rule assembly_summary: input: qc_csv=rules.qc_summary.output.csv, asbl=expand( - "results/{{project}}/report_prerequisites/assembly/{sample}_megahit.log", + "results/{{project}}/output/report/prerequisites/assembly/{sample}_megahit.log", sample=get_samples(), ), mapped=expand( - "results/{{project}}/report_prerequisites/assembly/{sample}_reads_mapped.txt", + "results/{{project}}/output/report/prerequisites/assembly/{sample}_reads_mapped.txt", + sample=get_samples(), + ), + csv_bins=expand( + "results/{{project}}/output/report/{sample}/{sample}_bin_summary.csv", + sample=get_samples(), + ), + csv_mags=expand( + "results/{{project}}/output/report/{sample}/{sample}_mags_summary.csv", sample=get_samples(), ), output: @@ -114,13 +122,15 @@ use rule qc_summary_report as assembly_report with: input: rules.assembly_summary.output.vis_csv, output: - report( - directory("results/{project}/output/report/all/assembly/"), - htmlindex="index.html", - category="3. Assembly results", - labels={ - "sample": "all samples", - }, + temp( + report( + directory("results/{project}/output/report/all/assembly/"), + htmlindex="index.html", + category="3. Assembly results", + labels={ + "sample": "all samples", + }, + ) ), params: pin_until="sample", @@ -144,4 +154,4 @@ rule cleanup_megahit_output: output: touch("results/{project}/megahit/{sample}_cleanup.done"), log: - "logs/{project}/assembly/{sample}_cleanup.log", \ No newline at end of file + "logs/{project}/assembly/{sample}_cleanup.log", diff --git a/workflow/rules/bin_qc.smk b/workflow/rules/bin_qc.smk index c2a47be..d9039eb 100644 --- a/workflow/rules/bin_qc.smk +++ b/workflow/rules/bin_qc.smk @@ -37,11 +37,8 @@ rule bin_summary_sample: tool="results/{project}/output/report/{sample}/{sample}_DASTool_summary.tsv", checkm="results/{project}/output/report/{sample}/checkm2_quality_report.tsv", gtdb="results/{project}/output/classification/bins/{sample}/{sample}.summary.tsv", - #args="results/{project}/output/ARGs/bins/{sample}/all_bins.done", output: csv_bins="results/{project}/output/report/{sample}/{sample}_bin_summary.csv", - csv_checkm="results/{project}/output/report/{sample}/{sample}_checkm2_summary.csv", - csv_dastool="results/{project}/output/report/{sample}/{sample}_DASTool_summary.csv", csv_tax="results/{project}/output/report/{sample}/{sample}_bin_taxonomy.csv", csv_mags="results/{project}/output/report/{sample}/{sample}_mags_summary.csv", params: @@ -60,12 +57,14 @@ use rule qc_summary_report as bin_sample_report with: input: "results/{project}/output/report/{sample}/{sample}_bin_summary.csv", output: - report( - directory("results/{project}/output/report/{sample}/bin/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.1 Summary", - labels={"sample": "{sample}"}, + temp( + report( + directory("results/{project}/output/report/{sample}/bin/"), + htmlindex="index.html", + category="4. Binning results", + subcategory="4.1 Summary", + labels={"sample": "{sample}"}, + ) ), params: pin_until="bin", @@ -77,64 +76,18 @@ use rule qc_summary_report as bin_sample_report with: "logs/{project}/report/{sample}/bin_rbt_csv.log", -use rule qc_summary_report as dastool_report with: - input: - "results/{project}/output/report/{sample}/{sample}_DASTool_summary.csv", - output: - report( - directory("results/{project}/output/report/{sample}/dastool/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.2 Quality control", - labels={ - "sample": "{sample}", - "tool": "DAS Tool", - }, - ), - params: - pin_until="bin", - styles="resources/report/tables/", - name="{sample}_DASTool_summary", - header="DAS Tool summary for sample {sample}", - pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/dastool_rbt_csv.log", - - -use rule qc_summary_report as checkm2_report with: - input: - "results/{project}/output/report/{sample}/{sample}_checkm2_summary.csv", - output: - report( - directory("results/{project}/output/report/{sample}/checkm2/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.2 Quality control", - labels={ - "sample": "{sample}", - "tool": "CheckM 2", - }, - ), - params: - pin_until="bin", - styles="resources/report/tables/", - name="{sample}_CheckM2_summary", - header="CheckM2 summary for sample {sample}", - pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/checkm2_rbt_csv.log", - - use rule qc_summary_report as taxonomy_report with: input: "results/{project}/output/report/{sample}/{sample}_bin_taxonomy.csv", output: - report( - directory("results/{project}/output/report/{sample}/taxonomy/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.3 Taxonomy classification", - labels={"sample": "{sample}"}, + temp( + report( + directory("results/{project}/output/report/{sample}/taxonomy/"), + htmlindex="index.html", + category="4. Binning results", + subcategory="4.3 Taxonomy classification", + labels={"sample": "{sample}"}, + ) ), params: pin_until="bin", @@ -150,12 +103,14 @@ use rule qc_summary_report as mag_report with: input: "results/{project}/output/report/{sample}/{sample}_mags_summary.csv", output: - report( - directory("results/{project}/output/report/{sample}/mags/"), - htmlindex="index.html", - category="5. Taxonomic classification", - subcategory="5.1 MAGs classification", - labels={"sample": "{sample}"}, + temp( + report( + directory("results/{project}/output/report/{sample}/mags/"), + htmlindex="index.html", + category="5. Taxonomic classification", + subcategory="5.1 MAGs classification", + labels={"sample": "{sample}"}, + ) ), params: pin_until="MAG", @@ -192,12 +147,14 @@ use rule qc_summary_report as bin_all_report with: input: "results/{project}/output/report/all/binning_summary.csv", output: - report( - directory("results/{project}/output/report/all/binning/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.1 Summary", - labels={"sample": "all"}, + temp( + report( + directory("results/{project}/output/report/all/binning/"), + htmlindex="index.html", + category="4. Binning results", + subcategory="4.1 Summary", + labels={"sample": "all"}, + ) ), params: pin_until="sample", diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 34b6932..6295b22 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -80,7 +80,7 @@ if config["gtdb"]["use-local"]: rule prepare_gtdb: output: - done=touch("results/GTDB_prep.done"), + done=temp(touch("results/GTDB_prep.done")), params: db_folder=get_gtdb_folder(), threads: 1 @@ -106,14 +106,14 @@ if config["gtdb"]["use-local"]: params: clf_outdir=lambda wildcards, output: Path(output.json).parent, json=lambda wildcards, output: Path(output.json).name, - threads: 32 + threads: 20 log: "logs/{project}/gtdbtk/{sample}_classify.log", conda: "../envs/gtdbtk.yaml" shell: "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--cpus {threads} --pplacer_cpus 30 " + "--cpus {threads} --pplacer_cpus {threads} " "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 4019493..4f9a245 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -60,7 +60,7 @@ def get_prefiltered_fastqs(wildcards): def get_host_map_statistics(wildcards): if config["host-filtering"]["do-host-filtering"]: logs = expand( - "results/{{project}}/report_prerequisites/qc/filter_host_{sample}.log", + "results/{{project}}/output/report/prerequisites/qc/filter_host_{sample}.log", sample=get_samples(), ) return logs @@ -91,15 +91,15 @@ def get_checkm2_db(): def get_filtered_fastqs(wildcards): return [ - "results/{project}/filtered/{sample}_R1.fastq", - "results/{project}/filtered/{sample}_R2.fastq", + "results/{project}/output/filtered_reads/{sample}_R1.fastq", + "results/{project}/output/filtered_reads/{sample}_R2.fastq", ] def get_filtered_gz_fastqs(wildcards): return [ - "results/{project}/filtered/{sample}_R1.fastq.gz", - "results/{project}/filtered/{sample}_R2.fastq.gz", + "results/{project}/output/filtered_reads/{sample}_R1.fastq.gz", + "results/{project}/output/filtered_reads/{sample}_R2.fastq.gz", ] diff --git a/workflow/rules/das_tool.smk b/workflow/rules/das_tool.smk index aa1f5ae..93c5032 100644 --- a/workflow/rules/das_tool.smk +++ b/workflow/rules/das_tool.smk @@ -113,13 +113,13 @@ if bins_for_sample: output: bins=directory("results/{project}/output/fastas/{sample}/bins/"), done=touch("results/{project}/binning/das_tool/{sample}_bins.done"), - threads: 64 + threads: 15 log: "logs/{project}/bins/{sample}/gz_bins.log", conda: "../envs/unix.yaml" shell: - "(pigz -k {input.bins}/*.fa && " + "(gzip -k {input.bins}/*.fa && " "mkdir -p {output.bins}/ && " "mv {input.bins}/*.fa.gz {output.bins}/ ) > {log} 2>&1" diff --git a/workflow/rules/host_filtering.smk b/workflow/rules/host_filtering.smk index d3a0a9e..e8729f8 100644 --- a/workflow/rules/host_filtering.smk +++ b/workflow/rules/host_filtering.smk @@ -79,13 +79,13 @@ rule filter_human: output: filtered=temp( expand( - "results/{{project}}/filtered/{{sample}}_{read}.fastq", + "results/{{project}}/output/filtered_reads/{{sample}}_{read}.fastq", read=["R1", "R2"], ) ), threads: 64 log: - "results/{project}/report_prerequisites/qc/filter_human_{sample}.log", + "results/{project}/output/report/prerequisites/qc/{sample}_filter_human.log", conda: "../envs/minimap2.yaml" shell: @@ -97,9 +97,9 @@ rule filter_human: rule gzip_filtered_reads: input: - "results/{project}/filtered/{sample}_{read}.fastq", + "results/{project}/output/filtered_reads/{sample}_{read}.fastq", output: - "results/{project}/filtered/{sample}_{read}.fastq.gz", + "results/{project}/output/filtered_reads/{sample}_{read}.fastq.gz", log: "logs/{project}/human_filtering/gzip_{sample}_{read}.log", threads: 20 @@ -147,7 +147,7 @@ if config["host-filtering"]["do-host-filtering"]: ), threads: 64 log: - "results/{project}/report_prerequisites/qc/filter_host_{sample}.log", + "results/{project}/output/report/prerequisites/qc/{sample}_filter_host.log", conda: "../envs/minimap2.yaml" shell: @@ -167,11 +167,11 @@ if config["host-filtering"]["do-host-filtering"]: sample=get_samples(), ), jsons=expand( - "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", + "results/{{project}}/output/report/prerequisites/qc/{sample}.fastp.json", sample=get_samples(), ), human_logs=expand( - "results/{{project}}/report_prerequisites/qc/filter_human_{sample}.log", + "results/{{project}}/output/report/prerequisites/qc/filter_human_{sample}.log", sample=get_samples(), ), host_logs=get_host_map_statistics, diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index 337bc45..e96359e 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -12,7 +12,7 @@ rule fastp: ] ), html=temp("results/{project}/trimmed/fastp/{sample}.html"), - json="results/{project}/report_prerequisites/qc/{sample}.fastp.json", + json="results/{project}/output/report/prerequisites/qc/{sample}.fastp.json", params: adapters=get_adapters, extra="--qualified_quality_phred {phred} --length_required {minlen}".format( @@ -78,25 +78,22 @@ rule multiqc: expand( [ "results/{{project}}/trimmed/fastp/{sample}.{read}_fastqc.zip", - "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", + "results/{{project}}/output/report/prerequisites/qc/{sample}.fastp.json", ], sample=get_samples(), read=["1", "2"], ), + config="config/multiqc_config.yaml", output: report( - "results/{project}/output/report/all/multiqc.html", - htmlindex="multiqc.html", + "results/{project}/output/report/all/multiqc_{project}.html", + htmlindex="multiqc_{project}.html", category="1. Quality control", labels={"sample": "all samples"}, ), - "results/{project}/output/report/all/multiqc_data.zip", + "results/{project}/output/report/all/multiqc_{project}.zip", params: - extra=( - "--zip-data-dir " - "--config config/multiqc_config.yaml " - "--title 'Results for data from {project} project'" - ), + extra=("--title 'Results for data from {project} project'"), use_input_files_only=True, log: "logs/{project}/multiqc.log", @@ -108,11 +105,11 @@ rule qc_summary: input: # move these to report_prerequistes jsons=expand( - "results/{{project}}/report_prerequisites/qc/{sample}.fastp.json", + "results/{{project}}/output/report/prerequisites/qc/{sample}.fastp.json", sample=get_samples(), ), human_logs=expand( - "results/{{project}}/report_prerequisites/qc/filter_human_{sample}.log", + "results/{{project}}/output/report/prerequisites/qc/{sample}_filter_human.log", sample=get_samples(), ), host_logs=get_host_map_statistics, @@ -134,13 +131,15 @@ rule qc_summary_report: input: rules.qc_summary.output.vis_csv, output: - report( - directory("results/{project}/output/report/all/quality_summary/"), - htmlindex="index.html", - category="1. Quality control", - labels={ - "sample": "all samples", - }, + temp( + report( + directory("results/{project}/output/report/all/quality_summary/"), + htmlindex="index.html", + category="1. Quality control", + labels={ + "sample": "all samples", + }, + ) ), params: pin_until="sample", diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index 149cbc1..6dd4bad 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -7,7 +7,7 @@ rule snakemake_report: input: # 1. Quality control report_input, - "results/{project}/output/report/all/multiqc.html", + "results/{project}/output/report/all/multiqc_{project}.html", # 2. Species diversity "results/{project}/output/report/all/quality_summary/", # 3. Assembly results @@ -18,14 +18,6 @@ rule snakemake_report: "results/{{project}}/output/report/{sample}/bin/", sample=get_samples(), ), - expand( - "results/{{project}}/output/report/{sample}/checkm2/", - sample=get_samples(), - ), - expand( - "results/{{project}}/output/report/{sample}/dastool/", - sample=get_samples(), - ), expand( "results/{{project}}/output/report/{sample}/taxonomy/", sample=get_samples(), diff --git a/workflow/scripts/assembly_summary.py b/workflow/scripts/assembly_summary.py index ee5fbda..f0c5efd 100644 --- a/workflow/scripts/assembly_summary.py +++ b/workflow/scripts/assembly_summary.py @@ -1,11 +1,14 @@ import pandas as pd import sys +import os sys.stderr = open(snakemake.log[0], "w") asbl_logs = snakemake.input.asbl txts = snakemake.input.mapped qc_csv = snakemake.input.qc_csv +csv_mags = snakemake.input.csv_mags +csv_bins = snakemake.input.csv_bins df=pd.read_csv(qc_csv,index_col="sample") samples=df.index.to_list() @@ -16,7 +19,7 @@ summary_sample_dict["#reads_after_filtering"] = df.loc[sample,"#reads_after_filtering"] for asbl_log in asbl_logs: - if asbl_log.rfind(sample) >= 0: + if asbl_log.rfind(f'{sample}_megahit') >= 0: with open(asbl_log, "r") as a_log: for line in a_log: if line.find("total") >= 0: @@ -43,12 +46,22 @@ summary_sample_dict[colname] = value for txt in txts: - if txt.rfind(sample) >= 0: + if txt.rfind(f'{sample}_reads_mapped') >= 0: with open(txt) as t: mapped=int(t.readline()) summary_sample_dict["#assembled_reads"] = mapped - summary_sample_dict["%assembled_reads"] = round(((mapped/summary_sample_dict["#reads_after_filtering"]) * 100),4) + summary_sample_dict["%assembled_reads"] = round(((mapped/summary_sample_dict["#reads_after_filtering"]) * 100),2) break + + for bin_file in csv_bins: + if bin_file.rfind(sample) >= 0: + + bin_df=pd.read_csv(bin_file) + summary_sample_dict["#bins"] = len(bin_df) + + mag_file=[file for file in csv_mags if os.path.basename(os.path.dirname(file)) == sample][0] + mag_df=pd.read_csv(mag_file) + summary_sample_dict["#MAGs"] = len(mag_df) summary_dict[sample] = summary_sample_dict diff --git a/workflow/scripts/bin_summary_sample.py b/workflow/scripts/bin_summary_sample.py index a880292..9876f3f 100644 --- a/workflow/scripts/bin_summary_sample.py +++ b/workflow/scripts/bin_summary_sample.py @@ -6,31 +6,25 @@ ## input files in_dastool = ( snakemake.input.tool -) # "ResMAG/results/autobrewer/das_tool/ABS_24_07/ABS_24_07_DASTool_summary.tsv" +) in_checkm = ( snakemake.input.checkm -) # "ResMAG/results/autobrewer/qc/checkm2/ABS_24_07/quality_report.tsv" +) in_gtdb = ( snakemake.input.gtdb -) # "ResMAG/results/autobrewer/classification/ABS_24_07/ABS_24_07.bac120.summary.tsv" +) ## output files ### csv csv_path_mags = ( snakemake.output.csv_mags -) # "ResMAG/results/autobrewer/report/ABS_24_07/mags_summary.csv" # +) csv_path_bins = ( snakemake.output.csv_bins -) # "ResMAG/results/autobrewer/report/ABS_24_07/bin_summary.csv" +) csv_path_tax = ( snakemake.output.csv_tax -) # "ResMAG/results/autobrewer/report/ABS_24_07/bin_taxonomy.csv" -csv_path_checkm = ( - snakemake.output.csv_checkm -) # "ResMAG/results/autobrewer/report/ABS_24_07/checkm_summary.csv" -csv_path_dastool = ( - snakemake.output.csv_dastool -) # "ResMAG/results/autobrewer/report/ABS_24_07/DASTool_summary.csv" +) ## params max_cont = snakemake.params.max_cont @@ -46,7 +40,7 @@ def save_csv_table(csv_path, summary_df): tool_df = pd.read_table(in_dastool) tool_df.drop(["bin_set"], axis=1, inplace=True) -# tool_df.rename({"bin":"bin"},axis=1,inplace=True) + tool_df.set_index("bin", inplace=True) col_order = [ "bin_score", @@ -64,9 +58,6 @@ def save_csv_table(csv_path, summary_df): del col_order[:2] tool_red_df = tool_df.drop(col_order, axis=1) -save_csv_table(csv_path_dastool, tool_red_df) - - checkm_df = pd.read_table(in_checkm) rm_list = ["Translation_Table_Used", "Additional_Notes"] checkm_df.drop(rm_list, axis=1, inplace=True) @@ -93,8 +84,6 @@ def save_csv_table(csv_path, summary_df): del col_order[:5] checkm_red_df = checkm_df.drop(col_order, axis=1) -save_csv_table(csv_path_checkm, checkm_red_df) - gtdb_df = pd.read_table(in_gtdb) gtdb_df.rename({"user_genome": "bin"}, axis=1, inplace=True) diff --git a/workflow/scripts/quality_summary.py b/workflow/scripts/quality_summary.py index cd9e0f2..6f282aa 100644 --- a/workflow/scripts/quality_summary.py +++ b/workflow/scripts/quality_summary.py @@ -11,8 +11,8 @@ hostname=snakemake.params.hostname host_logs = snakemake.input.host_logs -outfile = snakemake.output.csv #"/local/work/josefa/ResMAG/results/lanuv/output/report/all/seq_summary.csv" -outfile_vis = snakemake.output.vis_csv #"/local/work/josefa/ResMAG/results/lanuv/output/report/all/seq_summary_vis.csv" +outfile = snakemake.output.csv +outfile_vis = snakemake.output.vis_csv results_dict = {} @@ -29,21 +29,21 @@ total_pre_host_filt = int(jdata["summary"]["after_filtering"]["total_reads"]) sample_results_dict["#reads_afterQC"] = total_pre_host_filt - sample_results_dict["%reads_filtered_by_quality"] = round(((1 - (total_pre_host_filt/reads_before)) * 100),4) + sample_results_dict["%reads_filtered_by_quality"] = round(((1 - (total_pre_host_filt/reads_before)) * 100),2) sample_results_dict["#bp_afterQC"] = int(jdata["summary"]["after_filtering"]["total_bases"]) - sample_results_dict["%Q30_reads_afterQC"] = float(jdata["summary"]["after_filtering"]["q30_rate"]) * 100 + sample_results_dict["%Q30_reads_afterQC"] = round((float(jdata["summary"]["after_filtering"]["q30_rate"]) * 100),2) if snakemake.params.other_host: for host_log in host_logs: - if host_log.find(sample) >= 0: + if host_log.find(f'{sample}_filter_host') >= 0: with open(host_log, "r") as file: for line in file.readlines(): if line.find("processed") >= 0: non_host_reads = int(line.split()[2]) * 2 no_reads = total_pre_host_filt - non_host_reads - prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),4) + prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),2) break sample_results_dict[f"#{hostname}_reads"] = no_reads @@ -52,7 +52,7 @@ for human_log in human_logs: - if human_log.find(sample) >= 0: + if human_log.find(f'{sample}_filter_human') >= 0: with open(human_log, "r") as file: for line in file.readlines(): if line.find("processed") >= 0: @@ -61,7 +61,7 @@ no_reads = non_host_reads - non_human_reads else: no_reads = total_pre_host_filt - non_human_reads - prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),4) + prct_reads = round(((int(no_reads) / int(total_pre_host_filt)) * 100),2) break sample_results_dict[f"#human_reads"] = no_reads @@ -69,7 +69,7 @@ break sample_results_dict["#reads_after_filtering"] = non_human_reads - sample_results_dict["%reads_filtered_total"] = round(((1 - (non_human_reads/reads_before)) * 100),4) + sample_results_dict["%reads_filtered_total"] = round(((1 - (non_human_reads/reads_before)) * 100),2) results_dict[sample] = sample_results_dict From 6d07016a931564357f4d3ebece97030dda64955e Mon Sep 17 00:00:00 2001 From: josefawelling Date: Wed, 18 Feb 2026 08:09:12 +0000 Subject: [PATCH 05/15] update to snakemake 9, remove unnecessary rules --- config/config.yaml | 9 -- workflow/Snakefile | 20 ---- workflow/report/host_plot.rst | 2 - workflow/rules/analysis.smk | 28 ++--- workflow/rules/assembly.smk | 3 +- workflow/rules/bin_qc.smk | 7 +- workflow/rules/common.smk | 58 +++-------- workflow/rules/das_tool.smk | 48 ++------- workflow/rules/host_filtering.smk | 77 -------------- workflow/rules/metabat.smk | 11 -- workflow/rules/metacoag.smk | 17 +-- workflow/rules/qc.smk | 8 +- workflow/scripts/binner_control.py | 32 ------ workflow/scripts/diversity_summary.py | 106 ------------------- workflow/scripts/plot_host.py | 142 -------------------------- 15 files changed, 45 insertions(+), 523 deletions(-) delete mode 100644 workflow/report/host_plot.rst delete mode 100644 workflow/scripts/binner_control.py delete mode 100644 workflow/scripts/diversity_summary.py delete mode 100644 workflow/scripts/plot_host.py diff --git a/config/config.yaml b/config/config.yaml index cfcc4cb..15e2b23 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -3,12 +3,8 @@ pepfile: config/pep/config.yaml ## Please change to a name describing your project ## All results can be found under results/project-name/ project-name: "test" -# -run-date: "2026-02-03" data-handling: - # path to store data within the workflow - data: data/ # path where databases and reference genomes are stored resources: resources/ @@ -42,11 +38,6 @@ host-filtering: # minimum contig length used during assembly step min-contig-length: 300 -das-tool: - # cores to use for das_tool - threads: 64 - binner-list: ["metabat", "metacoag"] - MAG-criteria: min-completeness: 50.00 max-contamination: 30.00 diff --git a/workflow/Snakefile b/workflow/Snakefile index 2497d86..6a8051c 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -67,26 +67,6 @@ rule all: "results/{project}/output/report/report_{project}.zip", project=get_project(), ), - expand( - "results/{project}/megahit/{sample}_cleanup.done", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/binning/metabat/{sample}_cleanup.done", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/binning/metacoag/{sample}_cleanup.done", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/binning/das_tool/{sample}_cleanup.done", - project=get_project(), - sample=get_samples(), - ), # resistance analysis expand( "results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", diff --git a/workflow/report/host_plot.rst b/workflow/report/host_plot.rst deleted file mode 100644 index 1388dd2..0000000 --- a/workflow/report/host_plot.rst +++ /dev/null @@ -1,2 +0,0 @@ -A barplot for each sample, showing the percentage of human reads and if given also the percentage of reads of another host. -The reads are classified by mapping against a reference genome using `minimap2 `_ and `samtools `_. \ No newline at end of file diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 56f90f6..caabb88 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -37,23 +37,24 @@ rule gzip_proteins: # Plasmid analysis rule load_genomad_DB: output: + folder=get_genomad_DB_folder(), file=get_genomad_DB_file(), params: - folder=lambda wildcards, output: Path(output.file).parent.parent, + res_folder=lambda wildcards, output: Path(output.folder).parent, log: "logs/load_genomad_DB.log", conda: "../envs/genomad.yaml" shell: - "genomad download-database {params.folder}/ > {log} 2>&1" + "genomad download-database {params.res_folder}/ > {log} 2>&1" rule genomad_run: input: - db=rules.load_genomad_DB.output.file, + db=rules.load_genomad_DB.output.folder, asmbl=rules.gzip_assembly.output, output: - outdir=temp(directory("results/{project}/genomad/{sample}/")), + #outdir=temp(directory("results/{project}/genomad/{sample}/")), plasmid_tsv=temp( "results/{project}/genomad/{sample}/{sample}_summary/{sample}_plasmid_summary.tsv" ), @@ -61,7 +62,7 @@ rule genomad_run: "results/{project}/genomad/{sample}/{sample}_summary/{sample}_virus_summary.tsv" ), params: - db_folder=lambda wildcards, input: Path(input.db).parent, + outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent.parent, log: "logs/{project}/plasmids/{sample}_run.log", threads: 64 @@ -69,13 +70,12 @@ rule genomad_run: "../envs/genomad.yaml" shell: "genomad end-to-end --cleanup -t {threads} " - "{input.asmbl} {output.outdir}/ " - "{params.db_folder}/ > {log} 2>&1" + "{input.asmbl} {params.outdir}/ " + "{input.db}/ > {log} 2>&1" rule move_genomad_output: input: - folder=rules.genomad_run.output.outdir, plasmid_tsv=rules.genomad_run.output.plasmid_tsv, virus_tsv=rules.genomad_run.output.virus_tsv, output: @@ -131,11 +131,11 @@ rule CARD_annotation: json=get_card_db_file(), load=rules.CARD_load_DB_for_reads.output, output: - temp(touch("results/CARD_annotation.done")), + done=temp(touch("results/CARD_annotation.done")), + ann=get_card_annotation_file(), params: folder=lambda wildcards, input: Path(input.json).parent, file=lambda wildcards, input: Path(input.json).name, - ann=get_card_annotation_file(), log: "logs/CARD_annotation.log", conda: @@ -143,7 +143,7 @@ rule CARD_annotation: shell: "(cd {params.folder}/ && " "rgi card_annotation -i {params.file}) && " - "rgi load -i {input.json} --card_annotation {params.ann} --local > {log} 2>&1" + "rgi load -i {input.json} --card_annotation {output.ann} --local > {log} 2>&1" rule CARD_read_run: @@ -152,12 +152,6 @@ rule CARD_read_run: db=rules.CARD_annotation.output, output: txt="results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", - bam=temp( - "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.sorted.length_100.bam" - ), - bai=temp( - "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.sorted.length_100.bam.bai" - ), params: folder=lambda wildcards, output: Path(output.txt).parent, log: diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index d4019c2..90c122c 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -8,7 +8,6 @@ rule megahit: contigs=temp("results/{project}/megahit/{sample}/final.contigs.fa"), outdir=temp(directory("results/{project}/megahit/{sample}/")), log="results/{project}/output/report/prerequisites/assembly/{sample}_megahit.log", - done=touch("results/{project}/megahit/{sample}.done"), params: threshold=get_contig_length_threshold(), threads: 64 @@ -142,6 +141,7 @@ use rule qc_summary_report as assembly_report with: "logs/{project}/report/assembly_rbt_csv.log", +""" # remove megahit intermediate results when all dependent results are produced rule cleanup_megahit_output: input: @@ -155,3 +155,4 @@ rule cleanup_megahit_output: touch("results/{project}/megahit/{sample}_cleanup.done"), log: "logs/{project}/assembly/{sample}_cleanup.log", +""" diff --git a/workflow/rules/bin_qc.smk b/workflow/rules/bin_qc.smk index d9039eb..d2b87e3 100644 --- a/workflow/rules/bin_qc.smk +++ b/workflow/rules/bin_qc.smk @@ -17,10 +17,9 @@ rule checkm2_run: bins="results/{project}/output/fastas/{sample}/bins/", dbfile=get_checkm2_db(), output: - outdir=temp(directory("results/{project}/qc/checkm2/{sample}/")), stats="results/{project}/output/report/{sample}/checkm2_quality_report.tsv", params: - outname="quality_report.tsv", + outdir="results/{project}/qc/checkm2/{sample}/", log: "logs/{project}/checkm2/{sample}.log", threads: 24 @@ -28,8 +27,8 @@ rule checkm2_run: "../envs/checkm2.yaml" shell: "(checkm2 predict -x fa.gz --threads {threads} --force " - "--input {input.bins}/ --output-directory {output.outdir}/ && " - "cp {output.outdir}/{params.outname} {output.stats}) > {log} 2>&1" + "--input {input.bins}/ --output-directory {params.outdir}/ && " + "cp {params.outdir}/quality_report.tsv {output.stats}) > {log} 2>&1" rule bin_summary_sample: diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 4f9a245..41d9ba1 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -4,10 +4,6 @@ import os configfile: "config/config.yaml" -def get_data_path(): - return config["data-handling"]["data"] - - def get_resource_path(): return config["data-handling"]["resources"] @@ -20,6 +16,7 @@ def get_samples(): return list(pep.sample_table["sample_name"].values) +""" def get_fastqs(wildcards): file_r1 = pep.sample_table.loc[wildcards.sample]["fq1"] folder = str(Path(file_r1).parent) @@ -27,12 +24,15 @@ def get_fastqs(wildcards): filename_r2 = Path(pep.sample_table.loc[wildcards.sample]["fq2"]).name return [folder, filename_r1, filename_r2] +""" + -def get_local_fastqs(wildcards): - path = get_data_path() +def get_fastqs(wildcards): + file_r1 = pep.sample_table.loc[wildcards.sample]["fq1"] + file_r2 = pep.sample_table.loc[wildcards.sample]["fq2"] return ( - "{data}{{project}}/{{sample}}_R1.fastq.gz".format(data=path), - "{data}{{project}}/{{sample}}_R2.fastq.gz".format(data=path), + file_r1, + file_r2, ) @@ -125,39 +125,10 @@ def get_contig_length_threshold(): return config["min-contig-length"] -def get_binners(): - return config["das-tool"]["binner-list"] - - -def get_all_contig2bin_files(wildcards): - binners = get_binners() - file_list = [ - "".join( - [ - "results/{project}/output/contig2bins/{sample}/", - binner, - "_contig2bin.tsv", - ] - ) - for binner in binners - ] - return file_list - - -## reads in binner control file and returns list with paths to contig2bin files -## and a list with name of the binners that produced results -def get_paths_binner(wildcards): - file = "results/{}/binning/das_tool/binner_control_{}.csv".format( - wildcards.project, wildcards.sample - ) - lines = open(file).readlines() - paths = str(lines[0].rstrip("\n")) - binner = str(lines[1].rstrip("\n")) - return paths, binner - - def bins_for_sample(wildcards): - if len(get_paths_binner[0]) > 0: + binfolder = "results/{project}/binning/das_tool/{sample}/{sample}_DASTool_bins/" + bin_list = [b for b in os.listdir(binfolder) if b.endswith(".fa")] + if len(bin_list) > 0: return True else: return False @@ -201,8 +172,13 @@ def get_gtdb_folder(): return path +def get_genomad_DB_folder(): + path = "".join([get_resource_path(), "genomad_db/"]) + return path + + def get_genomad_DB_file(): - path = "".join([get_resource_path(), "genomad_db/names.dmp"]) + path = "".join([get_genomad_DB_folder(), "names.dmp"]) return path diff --git a/workflow/rules/das_tool.smk b/workflow/rules/das_tool.smk index 93c5032..d18fb88 100644 --- a/workflow/rules/das_tool.smk +++ b/workflow/rules/das_tool.smk @@ -5,7 +5,7 @@ rule postprocess_metabat: input: outdir=rules.metabat.output.outdir, output: - "results/{project}/output/contig2bins/{sample}/metabat_contig2bin.tsv", + tsv="results/{project}/output/contig2bins/{sample}/metabat_contig2bin.tsv", params: binner="metabat", prefix="bin", @@ -21,9 +21,8 @@ rule postprocess_metabat: rule postprocess_metacoag: input: c2bin="results/{project}/binning/metacoag/{sample}/contig_to_bin.tsv", - folder="results/{project}/binning/metacoag/{sample}/", output: - "results/{project}/output/contig2bins/{sample}/metacoag_contig2bin.tsv", + tsv="results/{project}/output/contig2bins/{sample}/metacoag_contig2bin.tsv", threads: 12 log: "logs/{project}/contig2bins/{sample}/postprocess_metacoag.log", @@ -31,30 +30,14 @@ rule postprocess_metacoag: "../envs/unix.yaml" shell: "(awk '{{print $1 \"\t\" $NF}}' {input.c2bin} | " - "sed 's/len=[0-9]*,/metacoag_/g' > {output}) > {log} 2>&1" - - -## tests if binner created bins and writes file for DAS Tool -rule binner_control: - input: - get_all_contig2bin_files, - output: - temp("results/{project}/binning/das_tool/binner_control_{sample}.csv"), - params: - get_binners(), - log: - "logs/{project}/das_tool/{sample}/binner_control.log", - conda: - "../envs/unix.yaml" - script: - "../scripts/binner_control.py" + "sed 's/len=[0-9]*,/metacoag_/g' > {output.tsv}) > {log} 2>&1" rule dastool_run: input: - binc=rules.binner_control.output, + metacoag=rules.postprocess_metacoag.output.tsv, + metabat=rules.postprocess_metabat.output.tsv, contigs=get_assembly, - asml_folder=rules.megahit.output.outdir, output: summary=temp( "results/{project}/binning/das_tool/{sample}/{sample}_DASTool_summary.tsv" @@ -67,10 +50,8 @@ rule dastool_run: "results/{project}/binning/das_tool/{sample}/{sample}_DASTool_bins/" ) ), - outdir=temp(directory("results/{project}/binning/das_tool/{sample}/")), - done=touch("results/{project}/binning/das_tool/{sample}_run.done"), params: - path_bin_list=get_paths_binner, + outdir=lambda wildcards, output: Path(output.bins).parent, threshold=0.001, threads: 64 log: @@ -79,10 +60,10 @@ rule dastool_run: "../envs/das_tool.yaml" shell: "DAS_Tool --debug --write_bins " - "-i {params.path_bin_list[0]} " - "-l {params.path_bin_list[1]} " + "-i {input.metabat},{input.metacoag} " + "-l metabat,metacoag " "-c {input.contigs} " - "-o {output.outdir}/{wildcards.sample} " + "-o {params.outdir}/{wildcards.sample} " "--score_threshold {params.threshold} " "--threads={threads} " "> {log} 2>&1 " @@ -136,14 +117,3 @@ if bins_for_sample: "../envs/python.yaml" script: "../scripts/move_MAGs.py" - - rule cleanup_dastool_output: - input: - folder=rules.dastool_run.output.outdir, - bins=rules.gzip_bins.output.done, - move=rules.move_dastool_output.output.done, - output: - done=touch("results/{project}/binning/das_tool/{sample}_cleanup.done"), - threads: 2 - log: - "logs/{project}/das_tool/{sample}/cleanup.log", diff --git a/workflow/rules/host_filtering.smk b/workflow/rules/host_filtering.smk index e8729f8..96a9c19 100644 --- a/workflow/rules/host_filtering.smk +++ b/workflow/rules/host_filtering.smk @@ -155,80 +155,3 @@ if config["host-filtering"]["do-host-filtering"]: "gzip -c > {output.filtered[0]} && " "samtools fastq -F 3584 -f 141 {input.bam} | " "gzip -c > {output.filtered[1]}) > {log} 2>&1" - - -## TODO -## change to using kaiju output or just to present human contamination -## all reports are done on human (+ optional other different host) filtered -"""rule diversity_summary: - input: - reports=expand( - "results/{{project}}/output/classification/reads/{sample}/{sample}_kraken2_report.tsv", - sample=get_samples(), - ), - jsons=expand( - "results/{{project}}/output/report/prerequisites/qc/{sample}.fastp.json", - sample=get_samples(), - ), - human_logs=expand( - "results/{{project}}/output/report/prerequisites/qc/filter_human_{sample}.log", - sample=get_samples(), - ), - host_logs=get_host_map_statistics, - output: - csv="results/{project}/output/report/all/diversity_summary.csv", - log: - "logs/{project}/kraken2/summary.log", - params: - other_host=config["host-filtering"]["do-host-filtering"], - hostname=config["host-filtering"]["host-name"], - threads: 2 - conda: - "../envs/python.yaml" - script: - "../scripts/diversity_summary.py" - - -use rule qc_summary_report as diversity_summary_report with: - input: - rules.diversity_summary.output.csv, - output: - report( - directory("results/{project}/output/report/all/diversity_summary/"), - htmlindex="index.html", - caption="../report/kraken.rst", - category="2. Species diversity", - labels={ - "sample": "all samples", - }, - ), - params: - pin_until="sample", - styles="resources/report/tables/", - name="diversity_summary", - header="Diversity summary based on mapping to host genome(s) and Kraken2", - pattern=config["tablular-config"], - log: - "logs/{project}/report/kraken2_rbt_csv.log", - - -rule create_host_plot: - input: - csv=rules.diversity_summary.output.csv, - output: - html=report( - "results/{project}/output/report/all/host_contamination.html", - caption="../report/host_plot.rst", - category="2. Species diversity", - labels={"sample": "all samples"}, - ), - params: - other_host=config["host-filtering"]["do-host-filtering"], - hostname=config["host-filtering"]["host-name"], - log: - "logs/{project}/report/host_plot.log", - conda: - "../envs/python.yaml" - script: - "../scripts/plot_host.py" -""" diff --git a/workflow/rules/metabat.smk b/workflow/rules/metabat.smk index 7d8328c..44e467b 100644 --- a/workflow/rules/metabat.smk +++ b/workflow/rules/metabat.smk @@ -36,14 +36,3 @@ rule metabat: "-m 1500 -s 100000 " "--minCorr 95 --minContigByCorr 300 --minSamples 1 " "-t {threads} -v -l -o {output}/{params.prefix} > {log} 2>&1" - - -rule cleanup_metabat_output: - input: - folder=rules.metabat.output.outdir, - dastool="results/{project}/binning/das_tool/{sample}/{sample}_DASTool_summary.tsv", - output: - done=touch("results/{project}/binning/metabat/{sample}_cleanup.done"), - threads: 2 - log: - "logs/{project}/metabat/{sample}/cleanup.log", diff --git a/workflow/rules/metacoag.smk b/workflow/rules/metacoag.smk index 49319a4..b5fc240 100644 --- a/workflow/rules/metacoag.smk +++ b/workflow/rules/metacoag.smk @@ -63,12 +63,10 @@ rule metacoag_run: contigs=get_assembly, gfa="results/{project}/binning_prep/{sample}/assembly_tree.gfa", abd="results/{project}/binning_prep/{sample}/abundance_metacoag.tsv", + #assembly folder needs to be there + folder=rules.megahit.output.outdir, output: out_tsv=temp("results/{project}/binning/metacoag/{sample}/contig_to_bin.tsv"), - folder=temp(directory("results/{project}/binning/metacoag/{sample}/")), - intermediate=temp( - "results/{project}/binning/metacoag/{sample}final.contigs.fa.normalized_contig_tetramers.pickle" - ), params: outdir=lambda wildcards, output: Path(output.out_tsv).parent, threads: 64 @@ -82,14 +80,3 @@ rule metacoag_run: "--output {params.outdir} --min_length 300 " "--bin_mg_threshold 0.2 --min_bin_size 100000 " "--nthreads {threads} > {log} 2>&1" - - -rule cleanup_metacoag_output: - input: - folder=rules.metacoag_run.output.folder, - dastool="results/{project}/binning/das_tool/{sample}/{sample}_DASTool_summary.tsv", - output: - done=touch("results/{project}/binning/metacoag/{sample}_cleanup.done"), - threads: 2 - log: - "logs/{project}/metacoag/{sample}/cleanup.log", diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index e96359e..a66bc9e 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -3,7 +3,7 @@ # version in this wrapper: fastp=1.0.1 rule fastp: input: - sample=get_local_fastqs, + sample=get_fastqs, output: trimmed=temp( [ @@ -54,12 +54,6 @@ rule fastqc: read=["1", "2"], ) ), - html=temp( - expand( - "results/{{project}}/trimmed/fastp/{{sample}}.{read}_fastqc.html", - read=["1", "2"], - ) - ), threads: 4 resources: mem_mb=1024, diff --git a/workflow/scripts/binner_control.py b/workflow/scripts/binner_control.py deleted file mode 100644 index a0e38cb..0000000 --- a/workflow/scripts/binner_control.py +++ /dev/null @@ -1,32 +0,0 @@ -import os -import sys - -sys.stderr = open(snakemake.log[0], "w") - -def check_files(files, binner_ls): - path_ls = [] - results_binner = binner_ls.copy() - - for binner in binner_ls: - for file_path in files: - if file_path.find(binner) >= 0: - # check if file exists and is not empty - if os.path.isfile(file_path) and os.path.getsize(file_path) == 0: - results_binner.remove(binner) - break - else: - path_ls.append(file_path) - break - - with open(output_csv, 'w') as f_out: - f_out.write("{paths}\n{binners}\n".format(paths=(",".join(path_ls)),binners=(",".join(results_binner)))) - - - -# Get the input and output file paths from command-line arguments -files = snakemake.input -binner_ls = snakemake.params[0] -output_csv = snakemake.output[0] - -# Call the function with the provided file paths -check_files(files, binner_ls) \ No newline at end of file diff --git a/workflow/scripts/diversity_summary.py b/workflow/scripts/diversity_summary.py deleted file mode 100644 index cea0b21..0000000 --- a/workflow/scripts/diversity_summary.py +++ /dev/null @@ -1,106 +0,0 @@ -import pandas as pd -import json -import sys - -sys.stderr = open(snakemake.log[0], "w") - -jsons = snakemake.input.jsons -human_logs = snakemake.input.human_logs -if snakemake.params.other_host: - hostname = snakemake.params.hostname - host_logs = snakemake.input.host_logs -reports = snakemake.input.reports -tax_ids = snakemake.params.taxid_dict -outfile = snakemake.output.csv - -# header for kraken report file -header_ls = ["prct", "total_reads", "lvl_reads", "lvl", "tax_id", "name"] -results_dict = {} - -for report in reports: - sample = report[report.rfind("/") + 1 : report.rfind("_kraken2_report")] - sample_results_dict = {} - - for json_file in jsons: - if json_file.find(sample) >= 0: - with open(json_file, "r") as read_file: - jdata = json.load(read_file) - # fastp json shows reads instead of read pairs as in the kraken reports - total_pre_host_filt = int( - (jdata["summary"]["after_filtering"]["total_reads"]) / 2 - ) - break - sample_results_dict["#total_reads"] = total_pre_host_filt - - if snakemake.params.other_host: - for host_log in host_logs: - if host_log.find(sample) >= 0: - with open(host_log, "r") as file: - for line in file.readlines(): - if line.find("processed") >= 0: - non_host_reads = int(line.split()[2]) - no_reads = total_pre_host_filt - non_host_reads - prct_reads = ( - int(no_reads) / int(total_pre_host_filt) - ) * 100 - break - - sample_results_dict[f"%{hostname}"] = "%.3f" % prct_reads - sample_results_dict[f"#reads_{hostname}"] = no_reads - break - - for human_log in human_logs: - if human_log.find(sample) >= 0: - with open(human_log, "r") as file: - for line in file.readlines(): - if line.find("processed") >= 0: - non_human_reads = int(line.split()[2]) - - if snakemake.params.other_host: - no_reads = ( - total_pre_host_filt - non_host_reads - non_human_reads - ) - - else: - no_reads = total_pre_host_filt - non_human_reads - - prct_reads = (int(no_reads) / int(total_pre_host_filt)) * 100 - break - - - sample_results_dict[f"%human"] = "%.3f" % prct_reads - sample_results_dict[f"#reads_human"] = no_reads - break - - report_df = pd.read_table(report, names=header_ls) - - for ref in tax_ids: - line = report_df.loc[report_df["tax_id"] == tax_ids[ref]] - - no_reads = line.iloc[0]["total_reads"] - prct_reads = (int(no_reads) / int(total_pre_host_filt)) * 100 - - sample_results_dict[f"%{ref}"] = "%.3f" % prct_reads - sample_results_dict[f"#reads_{ref}"] = no_reads - - # add unclassified reads - line = report_df.loc[report_df["tax_id"] == 0] - if not line.empty: - no_reads = line.iloc[0]["total_reads"] - prct_reads = (int(no_reads) / int(total_pre_host_filt)) * 100 - - else: - no_reads = 0 - prct_reads = 0 - - sample_results_dict["%unclassified"] = "%.3f" % prct_reads - sample_results_dict["#reads_unclassified"] = no_reads - - results_dict[sample] = sample_results_dict - - -results_df = pd.DataFrame.from_dict(results_dict, orient="index") -results_df.index.name = "sample" - -results_df = results_df.sort_index() -results_df.to_csv(outfile) diff --git a/workflow/scripts/plot_host.py b/workflow/scripts/plot_host.py deleted file mode 100644 index 87602cc..0000000 --- a/workflow/scripts/plot_host.py +++ /dev/null @@ -1,142 +0,0 @@ -import pandas as pd -import altair as alt -import sys - -## write to log file -sys.stderr = open(snakemake.log[0], "w") - -## input file -csv_in = snakemake.input.csv - -#input parameter -other_host = snakemake.params.other_host - -## output file -host_percentage_html = snakemake.output.html - -## variables -color_red = "#e03e3e" -color_green = "#6aa84f" - -## prepare dataframe -df=pd.read_csv(csv_in) - -human_cont_df=pd.DataFrame() -human_cont_df["sample"]=df["sample"] -human_cont_df["human"]=df["%human"].divide(100) - -if other_host: - hostname=snakemake.params.hostname - human_cont_df["host"]=df[f"%{hostname}"].divide(100) - - -slider = alt.binding_range( - min=0, max=100, step=0.5, name="maximum acceptable:" -) -# selector = alt.param(name='SelectorName', value=50, bind=slider) -selector = alt.selection_point( - name="SelectorName", fields=["max_contamination"], bind=slider, value=50 -) - -# plot for human percentage -if other_host: - human_base = ( - alt.Chart(human_cont_df) - .encode( - alt.X("human:Q") - .axis(format="%", labelFontSize=12, titleFontSize=12) - .title("Percentage of human reads") - .scale(domain=[0, 1]), - alt.Y("sample:N") - .axis(labelFontSize=12, titleFontSize=12), - ) - .add_params(selector) - .properties(width="container") - .interactive() - ) - -else: - title="Share of human reads among all reads (after QC)" - title_object = alt.TitleParams(title, anchor='middle',fontSize=14) - - human_base = ( - alt.Chart(human_cont_df, title=title_object) - .encode( - alt.X("human:Q") - .axis(format="%", labelFontSize=12, titleFontSize=12) - .title("Percentage of human reads") - .scale(domain=[0, 1]), - alt.Y("sample:N") - .axis(labelFontSize=12, titleFontSize=12), - ) - .add_params(selector) - .properties(width="container") - .interactive() - ) - - -human_bars = human_base.mark_bar().encode( - color=alt.condition( - (alt.datum.human * 100) >= selector.max_contamination, - alt.value(color_red), - alt.value(color_green), - ) -) - -human_text = human_base.mark_text( - align="center", - baseline="middle", - dx=30, - fontSize=12, -).encode( - text=alt.Text("human:Q", format=".2%"), -) - -human_full = human_bars + human_text - -#if second host -if other_host: - # plot for other host percentage - host_base = ( - alt.Chart(human_cont_df) - .encode( - alt.X("host:Q") - .axis(format="%", labelFontSize=12, titleFontSize=12) - .title(f"Percentage of {hostname} reads") - .scale(domain=[0, 1]), - alt.Y("sample:N") - .axis(labelFontSize=12, titleFontSize=12), - ) - .add_params(selector) - .properties(width="container") - .interactive() - ) - - host_bars = host_base.mark_bar().encode( - color=alt.condition( - (alt.datum.host * 100) >= selector.max_contamination, - alt.value(color_red), - alt.value(color_green), - ) - ) - - host_text = host_base.mark_text( - align="center", - baseline="middle", - dx=30, - fontSize=12, - ).encode( - text=alt.Text("host:Q", format=".2%"), - ) - host_full = host_bars + host_text - - title=f"Share of human and {hostname} reads among all reads (after QC)" - title_object = alt.TitleParams(title, anchor='middle',fontSize=14) - - all_chart=alt.vconcat(host_full, human_full, title=title_object) - all_chart.save(host_percentage_html) - -else: - human_full.save(host_percentage_html) - - From beb13220063703346d7e4b7c40ec8f8ae264e248 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Fri, 13 Mar 2026 12:48:07 +0000 Subject: [PATCH 06/15] sm9 compatibility of DBs, add priorities, clean up, --- README.md | 38 ++++--- config/config.yaml | 33 +++--- workflow/Snakefile | 23 ++-- workflow/envs/snakemake.yaml | 3 +- workflow/rules/analysis.smk | 70 ++++++------ workflow/rules/assembly.smk | 21 +--- workflow/rules/bin_qc.smk | 71 +++--------- workflow/rules/classify.smk | 146 ++++++++++--------------- workflow/rules/common.smk | 140 +++++++++++------------- workflow/rules/das_tool.smk | 5 +- workflow/rules/host_filtering.smk | 31 +----- workflow/rules/qc.smk | 18 --- workflow/rules/report.smk | 3 - workflow/scripts/bin_summary_all.py | 32 ------ workflow/scripts/bin_summary_sample.py | 44 +------- 15 files changed, 239 insertions(+), 439 deletions(-) delete mode 100644 workflow/scripts/bin_summary_all.py diff --git a/README.md b/README.md index 21adc6a..730ba36 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # ResMAG - name pending -[![Snakemake](https://img.shields.io/badge/snakemake-≥6.3.0-brightgreen.svg)](https://snakemake.github.io) +[![Snakemake](https://img.shields.io/badge/snakemake-≥9.0.0-brightgreen.svg)](https://snakemake.github.io) [![GitHub actions status](https://github.com///workflows/Tests/badge.svg?branch=main)](https://github.com///actions?query=branch%3Amain+workflow%3ATests) @@ -88,24 +88,26 @@ git clone https://github.com/IKIM-Essen/metagenomics_res.git ``` #### Download GTDB -The GTDB needs to be downloaded and decompressed, it requires about 110 Gb. -1. Change to the cloned workflow directory -2. Create a new folder `resources/gtdb/` and change to this directory -3. Download the latest version of GTDB - **or** - if you have already downloaded a version of GTDB move the `gtdbtk_data.tar.gz` file to `resources/gtdb/` -4. Unarchive the downloaded file -5. After successful step 4: the archive can be removed +The GTDB needs to be downloaded and decompressed, it requires about 140 Gb. + +1. Change to the directory where the GTDB should be stored +2. Download the latest or your desired version of GTDB + Please make sure this version is compatible with GTDB-tk version `2.6.1` + ``` + wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz + ``` +3. Decompress the downloaded archive + ``` + tar xzf gtdbtk_data.tar.gz + ``` +4. After successful step 3: the archive can be removed +5. Please specify the path to your decompressed GTDB in the config file (see [Configuring workflow](#configuring-workflow)) -``` -wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz -tar xvzf gtdbtk_data.tar.gz -``` #### Install Snakemake Create a snakemake environment using [mamba](https://mamba.readthedocs.io/en/latest/) via: - ```mamba create -c conda-forge -c bioconda -n snakemake snakemake=7.32.3``` + ```mamba create -c conda-forge -c bioconda -n snakemake snakemake snakemake-storage-plugin-fs``` For installation details, see the [instructions in the Snakemake documentation](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). @@ -114,8 +116,10 @@ For installation details, see the [instructions in the Snakemake documentation]( - Specify a project name (`project-name`) - Specify filtering options for human reads (`human-filtering`) - Specify host filtering options, if you have a non-human host (`host-filtering`) - - Specify options for GTDB database (see [Download GTDB](#Download-GTDB)) -2. Provide sample information in the `config/pep/samples.csv` file while keeping the header and the format: + - Specify options for different databases: + - GTDB database needs to be downloaded before (see [Download GTDB](#Download-GTDB)) + - other databases (kaiju, CheckM2, CARD, genomad) can be given as a local path or downloaded when running the pipeline +2. Provide sample information in the `config/pep/samples.csv` file while keeping the header and the format as: ``` sample_name,fq1,fq2 @@ -131,7 +135,7 @@ Test your configuration by performing a dry-run via ```snakemake --use-conda -n``` Executing the workflow: -```snakemake --use-conda --cores $N --rerun-incomplete``` +```snakemake --use-conda --cores $N -k``` using `$N` cores. It is recommended to use all available cores. diff --git a/config/config.yaml b/config/config.yaml index 15e2b23..04f8769 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -4,10 +4,6 @@ pepfile: config/pep/config.yaml ## All results can be found under results/project-name/ project-name: "test" -data-handling: - # path where databases and reference genomes are stored - resources: resources/ - quality-criteria: # minimal length of acceptable reads min-length-reads: 30 @@ -20,9 +16,9 @@ adapter-seqs: "-a CTGTCTCTTATACACATCT -g AGATGTGTATAAGAGACAG" human-filtering: # if you have a human reference genome on your device, # you can set use-local to True and specify the path to the human reference genome - use-local: False + use-local: True # path to a locally stored human reference genome - local-path: /groups/ds/databases_refGenomes/refGenomes/latest/human/GCA_000001405.29_GRCh38.p14_genomic.fna.gz + local-path: /groups/ds/databases_refGenomes/refGenomes/latest/human/GRCh38_latest_genomic.fna.gz # if use-local = False the reference genome is downloaded via the following url download-path: https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.fna.gz @@ -43,23 +39,30 @@ MAG-criteria: max-contamination: 30.00 ### Handling of different databases +## The GTDB has about 140 Gb +# please specify the folder where the decompressed database is stored +gtdb: /groups/ds/databases_refGenomes/databases/latest/gtdb/release226/ + kaiju: + use-local: True + # path of kaiju DB .fmi file + # if use-local False: DB is downloaded and saved to this path + fmi-file: /groups/ds/databases_refGenomes/databases/latest/kaiju/refseq_nr/unpacked_2024-08-13/kaiju_db_refseq_nr.fmi download: https://kaiju-idx.s3.eu-central-1.amazonaws.com/2024/kaiju_db_refseq_nr_2024-08-13.tgz - fmi-file: kaiju_db/kaiju_db_refseq_nr.fmi -checkm2: CheckM2_database/uniref100.KO.1.dmnd +checkm2: + use-local: True + db-folder: /groups/ds/databases_refGenomes/databases/latest/checkM2/ -## The GTDB has about 100 Gb -## if use-local is set to False this will be downloaded during the pipeline run -gtdb: +# folder of genomad DB +genomad: use-local: True - # if use-local is set to True, please specify the folder where the decompressed database is stored - # this path is expected to lay under the data-handling resources folder - db-folder: gtdb/release226/ + db-folder: /groups/ds/databases_refGenomes/databases/latest/genomad/ card: + use-local: True + db-file: /groups/ds/databases_refGenomes/databases/latest/card/card.json version: v4.0.1 - dbfile: card.json url: https://card.mcmaster.ca/latest/data ## string term used for formatting output tables diff --git a/workflow/Snakefile b/workflow/Snakefile index 6a8051c..9ddbf62 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -5,7 +5,7 @@ from snakemake.utils import min_version -min_version("6.3.0") +min_version("9.0.0") configfile: "config/config.yaml" @@ -45,6 +45,11 @@ rule all: project=get_project(), sample=get_samples(), ), + expand( + "results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", + project=get_project(), + sample=get_samples(), + ), # fasta output expand( "results/{project}/output/fastas/{sample}/bins/", @@ -69,21 +74,7 @@ rule all: ), # resistance analysis expand( - "results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", - project=get_project(), - sample=get_samples(), - ), - expand( - "results/{project}/output/resistance/CARD/assembly/{sample}/{sample}_assembly_ARGs.csv", - project=get_project(), - sample=get_samples(), - ), - - -""" - expand( - "results/{project}/output/resistance/CARD/mags/{sample}/all_mags.done", + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", project=get_project(), sample=get_samples(), ), -""" diff --git a/workflow/envs/snakemake.yaml b/workflow/envs/snakemake.yaml index 612d3a1..9551cc7 100644 --- a/workflow/envs/snakemake.yaml +++ b/workflow/envs/snakemake.yaml @@ -2,4 +2,5 @@ channels: - conda-forge - bioconda dependencies: - - snakemake=7.32.4 \ No newline at end of file + - snakemake=9.16.3 + - snakemake-storage-plugin-fs=1.1.3 diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index caabb88..660094a 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -26,6 +26,7 @@ rule gzip_proteins: fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna.gz", gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff.gz", threads: 20 + priority: 1 log: "logs/{project}/proteins/{sample}_gzip.log", conda: @@ -35,26 +36,28 @@ rule gzip_proteins: # Plasmid analysis -rule load_genomad_DB: - output: - folder=get_genomad_DB_folder(), - file=get_genomad_DB_file(), - params: - res_folder=lambda wildcards, output: Path(output.folder).parent, - log: - "logs/load_genomad_DB.log", - conda: - "../envs/genomad.yaml" - shell: - "genomad download-database {params.res_folder}/ > {log} 2>&1" + +if not config["genomad"]["use-local"]: + + rule load_genomad_DB: + output: + folder=get_genomad_DB_folder(), + file=get_genomad_DB_file(), + params: + res_folder=lambda wildcards, output: Path(output.folder).parent, + log: + "logs/load_genomad_DB.log", + conda: + "../envs/genomad.yaml" + shell: + "genomad download-database {params.res_folder}/ > {log} 2>&1" rule genomad_run: input: - db=rules.load_genomad_DB.output.folder, + db=get_genomad_DB_folder(), asmbl=rules.gzip_assembly.output, output: - #outdir=temp(directory("results/{project}/genomad/{sample}/")), plasmid_tsv=temp( "results/{project}/genomad/{sample}/{sample}_summary/{sample}_plasmid_summary.tsv" ), @@ -94,22 +97,24 @@ rule move_genomad_output: # Resistance analysis -rule download_CARD_data: - output: - json=get_card_db_file(), - params: - download=config["card"]["url"], - folder=lambda wildcards, output: Path(output.json).parent, - filename=lambda wildcards, output: Path(output.json).name, - log: - "logs/CARD_data_download.log", - threads: 30 - conda: - "../envs/unix.yaml" - shell: - "(cd {params.folder} && " - "wget {params.download} && " - "tar -xvf data ./{params.filename}) > {log} 2>&1" +if not config["card"]["use-local"]: + + rule download_CARD_data: + output: + json=get_card_db_file(), + params: + download=config["card"]["url"], + folder=lambda wildcards, output: Path(output.json).parent, + filename=lambda wildcards, output: Path(output.json).name, + log: + "logs/CARD_data_download.log", + threads: 30 + conda: + "../envs/unix.yaml" + shell: + "(cd {params.folder} && " + "wget {params.download} && " + "tar -xvf data ./{params.filename}) > {log} 2>&1" rule CARD_load_DB_for_reads: @@ -131,7 +136,7 @@ rule CARD_annotation: json=get_card_db_file(), load=rules.CARD_load_DB_for_reads.output, output: - done=temp(touch("results/CARD_annotation.done")), + #done=temp(touch("results/CARD_annotation.done")), ann=get_card_annotation_file(), params: folder=lambda wildcards, input: Path(input.json).parent, @@ -165,6 +170,7 @@ rule CARD_read_run: "--clean -n {threads} > {log} 2>&1" +""" rule CARD_read_sample_summary: input: txt=rules.CARD_read_run.output.txt, @@ -181,6 +187,7 @@ rule CARD_read_sample_summary: "../scripts/arg_summary_sample.py" + # updates CARD database to use for contigs instead of reads # read based classification must be finished before rule CARD_load_DB: @@ -253,3 +260,4 @@ rule wrap_mag_ARGs: touch("results/{project}/output/ARGs/mags/{sample}/all_mags.done"), log: "logs/{project}/ARGs/mags/{sample}/all_mags.log", +""" diff --git a/workflow/rules/assembly.smk b/workflow/rules/assembly.smk index 90c122c..6596587 100644 --- a/workflow/rules/assembly.smk +++ b/workflow/rules/assembly.smk @@ -11,6 +11,9 @@ rule megahit: params: threshold=get_contig_length_threshold(), threads: 64 + priority: 2 + resources: + heavy=2, log: "logs/{project}/assembly/{sample}_megahit.log", conda: @@ -79,6 +82,7 @@ rule gzip_assembly: output: "results/{project}/output/fastas/{sample}/{sample}.fa.gz", threads: 20 + priority: 1 log: "logs/{project}/assembly/{sample}_gzip.log", conda: @@ -139,20 +143,3 @@ use rule qc_summary_report as assembly_report with: pattern=config["tablular-config"], log: "logs/{project}/report/assembly_rbt_csv.log", - - -""" -# remove megahit intermediate results when all dependent results are produced -rule cleanup_megahit_output: - input: - # folder to remove - asmbl_folder=rules.megahit.output.outdir, - #dependent results - gz_asmbl=rules.gzip_assembly.output, - asmbl_summary=rules.assembly_summary.output.csv, - binning_done="results/{project}/binning/das_tool/{sample}_run.done", - output: - touch("results/{project}/megahit/{sample}_cleanup.done"), - log: - "logs/{project}/assembly/{sample}_cleanup.log", -""" diff --git a/workflow/rules/bin_qc.smk b/workflow/rules/bin_qc.smk index d2b87e3..54987b4 100644 --- a/workflow/rules/bin_qc.smk +++ b/workflow/rules/bin_qc.smk @@ -1,15 +1,17 @@ ## bin QC -rule checkm2_DB_download: - output: - dbfile=get_checkm2_db(), #"{}/{}".format(config["data-handling"]["resources"], config["checkm2"]), - params: - direct=lambda wildcards, output: Path(output.dbfile).parent.parent, - log: - "logs/checkm2_DB_download.log", - conda: - "../envs/checkm2.yaml" - shell: - "checkm2 database --download --path {params.direct} > {log} 2>&1" +if not config["checkm2"]["use-local"]: + + rule checkm2_DB_download: + output: + dbfile=get_checkm2_db(), + params: + direct=get_checkm2_db_folder(), + log: + "logs/checkm2_DB_download.log", + conda: + "../envs/checkm2.yaml" + shell: + "checkm2 database --download --path {params.direct} > {log} 2>&1" rule checkm2_run: @@ -27,7 +29,8 @@ rule checkm2_run: "../envs/checkm2.yaml" shell: "(checkm2 predict -x fa.gz --threads {threads} --force " - "--input {input.bins}/ --output-directory {params.outdir}/ && " + "--input {input.bins}/ --output-directory {params.outdir}/ " + "--database_path {input.dbfile} && " "cp {params.outdir}/quality_report.tsv {output.stats}) > {log} 2>&1" @@ -119,47 +122,3 @@ use rule qc_summary_report as mag_report with: pattern=config["tablular-config"], log: "logs/{project}/report/{sample}/mag_rbt_csv.log", - - -rule bin_summary_all: - input: - csv_mags=expand( - "results/{{project}}/output/report/{sample}/{sample}_mags_summary.csv", - sample=get_samples(), - ), - csv_bins=expand( - "results/{{project}}/output/report/{sample}/{sample}_bin_summary.csv", - sample=get_samples(), - ), - output: - "results/{project}/output/report/all/binning_summary.csv", - log: - "logs/{project}/bin_summary/all.log", - threads: 4 - conda: - "../envs/python.yaml" - script: - "../scripts/bin_summary_all.py" - - -use rule qc_summary_report as bin_all_report with: - input: - "results/{project}/output/report/all/binning_summary.csv", - output: - temp( - report( - directory("results/{project}/output/report/all/binning/"), - htmlindex="index.html", - category="4. Binning results", - subcategory="4.1 Summary", - labels={"sample": "all"}, - ) - ), - params: - pin_until="sample", - styles="resources/report/tables/", - name="bin_summary", - header="Bin summary for all samples", - pattern=config["tablular-config"], - log: - "logs/{project}/report/all_bin_rbt_csv.log", diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 6295b22..471e5ae 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -1,17 +1,20 @@ -rule download_kaiju: - output: - db_files=get_kaiju_files(), - params: - download=config["kaiju"]["download"], - db_folder=lambda wildcards, output: Path(output.db_files[0]).parent, - log: - "logs/kaiju_DB_download.log", - conda: - "../envs/unix.yaml" - shell: - "(mkdir -p {params.db_folder} && " - "wget -c {params.download} -O - | " - "tar -zxv -C {params.db_folder}) > {log} 2>&1" +# if there is no local database version to use, it is downloaded +if not config["kaiju"]["use-local"]: + + rule download_kaiju: + output: + db_files=get_kaiju_files(), + params: + download=config["kaiju"]["download"], + db_folder=lambda wildcards, output: Path(output.db_files[0]).parent, + log: + "logs/kaiju_DB_download.log", + conda: + "../envs/unix.yaml" + shell: + "(mkdir -p {params.db_folder} && " + "wget -c {params.download} -O - | " + "tar -zxv -C {params.db_folder}) > {log} 2>&1" rule run_kaiju: @@ -76,86 +79,49 @@ rule kaiju2krona: "ktImportText -o {output.html} {output.krona}) > {log} 2>&1" -if config["gtdb"]["use-local"]: - - rule prepare_gtdb: - output: - done=temp(touch("results/GTDB_prep.done")), - params: - db_folder=get_gtdb_folder(), - threads: 1 - log: - "logs/GTDB_prep.log", - conda: - "../envs/gtdbtk.yaml" - shell: - "(conda env config vars set " - "GTDBTK_DATA_PATH='{params.db_folder}') > {log} 2>&1" - - rule gtdbtk_classify_wf: - input: - bins=rules.gzip_bins.output.bins, - db_prep=rules.prepare_gtdb.output.done, - output: - json="results/{project}/output/classification/bins/{sample}/gtdbtk.json", - outdir=temp( - directory( - "results/{project}/output/classification/bins/{sample}/gtdbtk/" - ) - ), - params: - clf_outdir=lambda wildcards, output: Path(output.json).parent, - json=lambda wildcards, output: Path(output.json).name, - threads: 20 - log: - "logs/{project}/gtdbtk/{sample}_classify.log", - conda: - "../envs/gtdbtk.yaml" - shell: - "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--cpus {threads} --pplacer_cpus {threads} " - "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " - "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" - -else: +rule prepare_gtdb: + output: + done=temp(touch("results/GTDB_prep.done")), + params: + db_folder=get_gtdb_folder(), + threads: 1 + log: + "logs/GTDB_prep.log", + conda: + "../envs/gtdbtk.yaml" + shell: + "(conda env config vars set " + "GTDBTK_DATA_PATH='{params.db_folder}') > {log} 2>&1" - rule download_GTDB: - output: - done=touch("results/gtdbtk_download_DB.done"), - log: - "logs/gtdbtk_download_DB.log", - threads: 64 - conda: - "../envs/gtdbtk.yaml" - shell: - "download-db.sh > {log} 2>&1" - rule gtdbtk_classify_wf: - input: - bins=rules.gzip_bins.output.bins, - load=rules.download_GTDB.output.done, - output: - json="results/{project}/output/classification/bins/{sample}/gtdbtk.json", - outdir=temp( - directory( - "results/{project}/output/classification/bins/{sample}/gtdbtk/" - ) - ), - params: - clf_outdir=lambda wildcards, output: Path(output.json).parent, - json=lambda wildcards, output: Path(output.json).name, - threads: 64 - log: - "logs/{project}/gtdbtk/{sample}_classify.log", - conda: - "../envs/gtdbtk.yaml" - shell: - "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " - "--cpus {threads} --pplacer_cpus 30 " - "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " - "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" +rule gtdbtk_classify_wf: + input: + bins=rules.gzip_bins.output.bins, + db_prep=rules.prepare_gtdb.output.done, + output: + json="results/{project}/output/classification/bins/{sample}/gtdbtk.json", + outdir=temp( + directory("results/{project}/output/classification/bins/{sample}/gtdbtk/") + ), + params: + clf_outdir=lambda wildcards, output: Path(output.json).parent, + json=lambda wildcards, output: Path(output.json).name, + threads: 40 + # to ensure only one of these commands are run at the same time + resources: + heavy=1, + log: + "logs/{project}/gtdbtk/{sample}_classify.log", + conda: + "../envs/gtdbtk.yaml" + shell: + "(gtdbtk classify_wf --prefix {wildcards.sample} -x fa.gz " + "--cpus {threads} --pplacer_cpus {threads} " + "--genome_dir {input.bins}/ --out_dir {output.outdir}/ && " + "cp {output.outdir}/{params.json} {params.clf_outdir}/) > {log} 2>&1" +# combines bacterial and archaeal rule gtdb_summary: input: json=rules.gtdbtk_classify_wf.output.json, diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 41d9ba1..5995d50 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -4,10 +4,6 @@ import os configfile: "config/config.yaml" -def get_resource_path(): - return config["data-handling"]["resources"] - - def get_project(): return config["project-name"] @@ -16,17 +12,6 @@ def get_samples(): return list(pep.sample_table["sample_name"].values) -""" -def get_fastqs(wildcards): - file_r1 = pep.sample_table.loc[wildcards.sample]["fq1"] - folder = str(Path(file_r1).parent) - filename_r1 = Path(file_r1).name - filename_r2 = Path(pep.sample_table.loc[wildcards.sample]["fq2"]).name - return [folder, filename_r1, filename_r2] - -""" - - def get_fastqs(wildcards): file_r1 = pep.sample_table.loc[wildcards.sample]["fq1"] file_r2 = pep.sample_table.loc[wildcards.sample]["fq2"] @@ -68,27 +53,6 @@ def get_host_map_statistics(wildcards): return [] -def get_human_ref(): - if config["human-filtering"]["use-local"]: - path = config["human-filtering"]["local-path"] - else: - path = config["human-filtering"]["download-path"] - filename = path.split("/")[-1] - local_ref = "".join([get_resource_path(), "ref_genome/", filename]) - return local_ref - - -def get_human_local_folder(): - path = config["human-filtering"]["local-path"] - folder = Path(path).parent - return folder - - -def get_checkm2_db(): - file = "{}{}".format(get_resource_path(), config["checkm2"]) - return file - - def get_filtered_fastqs(wildcards): return [ "results/{project}/output/filtered_reads/{sample}_R1.fastq", @@ -111,15 +75,6 @@ def get_gz_assembly(wildcards): return "results/{project}/output/fastas/{sample}/{sample}.fa.gz" -def get_kaiju_files(): - file = "".join([get_resource_path(), config["kaiju"]["fmi-file"]]) - path = str(Path(file).parent) - fmi = Path(file).name - names = ["nodes.dmp", fmi, "names.dmp"] - files = ["/".join([path, name]) for name in names] - return files - - ## binning parameters def get_contig_length_threshold(): return config["min-contig-length"] @@ -134,47 +89,49 @@ def bins_for_sample(wildcards): return False -def get_mag_fa(wildcards): - folder = "results/{}/output/fastas/{}/mags/".format( - wildcards.project, wildcards.sample - ) - files = [ - os.path.join(folder, binID) - for binID in os.listdir(folder) - if binID.endswith("fa.gz") - ] +def get_human_ref_download(): + return config["human-filtering"]["download-path"] + + +def get_human_ref(): + if config["human-filtering"]["use-local"]: + path = config["human-filtering"]["local-path"] + return path + else: + path = get_human_ref_download() + filename = path.split("/")[-1] + local_ref = "".join(["resources/ref_genome/", filename]) + return local_ref + + +def get_host_ref(): + return config["host-filtering"]["ref-genome"] + + +def get_kaiju_files(): + file = config["kaiju"]["fmi-file"] + path = str(Path(file).parent) + fmi = Path(file).name + names = ["nodes.dmp", fmi, "names.dmp"] + files = ["/".join([path, name]) for name in names] return files -def get_binIDs_for_sample(wildcards): - folder = "results/{}/output/fastas/{}/mags/".format( - wildcards.project, wildcards.sample - ) - binIDs = [binID for binID in os.listdir(folder) if binID.endswith("fa.gz")] - return binIDs +def get_checkm2_db_folder(): + return config["checkm2"]["db-folder"] -def get_mag_ARGs(wildcards): - bin_fastas = (get_mag_fa(wildcards),) - bin_fastas = list(bin_fastas)[0] - binIDs = [os.path.basename(binID) for binID in bin_fastas] - folder = "results/{}/output/ARGs/mags/{}/".format( - wildcards.project, wildcards.sample - ) - arg_files = [ - os.path.join(folder, binID.replace(".fa.gz", ".txt")) for binID in binIDs - ] - return arg_files +def get_checkm2_db(): + path = "".join([get_checkm2_db_folder(), "CheckM2_database/uniref100.KO.1.dmnd"]) + return path def get_gtdb_folder(): - path = "".join([get_resource_path(), config["gtdb"]["db-folder"]]) - return path + return config["gtdb"] def get_genomad_DB_folder(): - path = "".join([get_resource_path(), "genomad_db/"]) - return path + return config["genomad"]["db-folder"] def get_genomad_DB_file(): @@ -183,11 +140,38 @@ def get_genomad_DB_file(): def get_card_db_file(): - path = "".join([get_resource_path(), "CARD_db/", config["card"]["dbfile"]]) - return path + return config["card"]["db-file"] def get_card_annotation_file(): version = config["card"]["version"] - path = "".join([get_resource_path(), "CARD_db/card_database_", version, ".fasta"]) + folder = str(Path(get_card_db_file()).parent) + path = "".join([folder, "/card_database_", version, ".fasta"]) return path + + +""" +def get_mag_fa(wildcards): + folder = "results/{}/output/fastas/{}/mags/".format( + wildcards.project, wildcards.sample + ) + files = [ + os.path.join(folder, binID) + for binID in os.listdir(folder) + if binID.endswith("fa.gz") + ] + return files + + +def get_mag_ARGs(wildcards): + bin_fastas = (get_mag_fa(wildcards),) + bin_fastas = list(bin_fastas)[0] + binIDs = [os.path.basename(binID) for binID in bin_fastas] + folder = "results/{}/output/ARGs/mags/{}/".format( + wildcards.project, wildcards.sample + ) + arg_files = [ + os.path.join(folder, binID.replace(".fa.gz", ".txt")) for binID in binIDs + ] + return arg_files +""" diff --git a/workflow/rules/das_tool.smk b/workflow/rules/das_tool.smk index d18fb88..81238bd 100644 --- a/workflow/rules/das_tool.smk +++ b/workflow/rules/das_tool.smk @@ -54,6 +54,8 @@ rule dastool_run: outdir=lambda wildcards, output: Path(output.bins).parent, threshold=0.001, threads: 64 + resources: + heavy=2, log: "logs/{project}/das_tool/{sample}/das_tool_run.log", conda: @@ -94,7 +96,8 @@ if bins_for_sample: output: bins=directory("results/{project}/output/fastas/{sample}/bins/"), done=touch("results/{project}/binning/das_tool/{sample}_bins.done"), - threads: 15 + threads: 10 + priority: 1 log: "logs/{project}/bins/{sample}/gz_bins.log", conda: diff --git a/workflow/rules/host_filtering.smk b/workflow/rules/host_filtering.smk index 96a9c19..258e1fd 100644 --- a/workflow/rules/host_filtering.smk +++ b/workflow/rules/host_filtering.smk @@ -1,33 +1,13 @@ from pathlib import Path - -if config["human-filtering"]["use-local"]: - - rule copy_local_human_ref: - output: - fasta=get_human_ref(), - params: - local=get_human_local_folder(), - folder=lambda wildcards, output: Path(output.fasta).parent, - file=lambda wildcards, output: Path(output.fasta).name, - log: - "logs/human_ref_local_copy.log", - group: - "refGenome_depended" - conda: - "../envs/unix.yaml" - shell: - "(mkdir -p {params.folder} && " - "tar cpfz - -C {params.local} {params.file} | " - "(cd {params.folder} ; tar xpfz -)) > {log} 2>&1" - -else: +# Download human reference genome if not local file is given +if not config["human-filtering"]["use-local"]: rule download_human_ref: output: fasta=get_human_ref(), params: - download=config["human-filtering"]["download-path"], + download=get_human_ref_download(), folder=lambda wildcards, output: Path(output.fasta).parent, log: "logs/human_ref_download.log", @@ -103,6 +83,7 @@ rule gzip_filtered_reads: log: "logs/{project}/human_filtering/gzip_{sample}_{read}.log", threads: 20 + priority: 1 conda: "../envs/unix.yaml" shell: @@ -116,11 +97,9 @@ if config["host-filtering"]["do-host-filtering"]: use rule map_to_human as map_to_host with: input: fastqs=get_trimmed_fastqs, - ref=config["host-filtering"]["ref-genome"], + ref=get_host_ref(), output: bam=temp("results/{project}/host_filtering/alignments/{sample}.bam"), - params: - ref=config["host-filtering"]["ref-genome"], threads: 20 log: "logs/{project}/host_filtering/map_to_host_{sample}.log", diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index a66bc9e..4f2d10d 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -26,24 +26,6 @@ rule fastp: "v7.1.0/bio/fastp" -"""# version in this wrapper: fastqc=0.12.1 -rule fastqc: - input: - rules.fastp.output.trimmed, - #get_trimmed_fastqs, - output: - html=temp("results/{project}/qc/fastqc/{sample}_trimmed.html"), - zip=temp("results/{project}/qc/fastqc/{sample}_trimmed_fastqc.zip"), - threads: 4 - resources: - mem_mb=1024, - log: - "logs/{project}/fastqc/{sample}.log", - wrapper: - "v7.6.0/bio/fastqc" -""" - - rule fastqc: input: rules.fastp.output.trimmed, diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index 6dd4bad..e8c401e 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -12,8 +12,6 @@ rule snakemake_report: "results/{project}/output/report/all/quality_summary/", # 3. Assembly results "results/{project}/output/report/all/assembly/", - # 4. Binning results - "results/{project}/output/report/all/binning/", expand( "results/{{project}}/output/report/{sample}/bin/", sample=get_samples(), @@ -42,6 +40,5 @@ rule snakemake_report: "../envs/snakemake.yaml" shell: "snakemake --nolock --report {output} --report-stylesheet {params.style} " - "> {log} 2>&1" #"{params.for_testing} " diff --git a/workflow/scripts/bin_summary_all.py b/workflow/scripts/bin_summary_all.py deleted file mode 100644 index 57ff66d..0000000 --- a/workflow/scripts/bin_summary_all.py +++ /dev/null @@ -1,32 +0,0 @@ -import pandas as pd -import sys -import os - -sys.stderr = open(snakemake.log[0], "w") -csv_mags = snakemake.input.csv_mags -csv_bins = snakemake.input.csv_bins - -summary_dict = {} -for bin_file in csv_bins: - sample = os.path.basename(os.path.dirname(bin_file)) - - summary_sample_dict = {} - - mag_file=[file for file in csv_mags if os.path.basename(os.path.dirname(file)) == sample][0] - mag_df=pd.read_csv(mag_file) - summary_sample_dict["# MAGs"] = len(mag_df) - - bin_df=pd.read_csv(bin_file) - summary_sample_dict["# bins"] = len(bin_df) - - summary_sample_dict["least contigs"] = bin_df["contigs"].min() - - summary_sample_dict["highest N50"] = f'{(bin_df["contig_N50"].max()):,}' - - summary_dict[sample]=summary_sample_dict - -summary_df = pd.DataFrame.from_dict(summary_dict, orient="index") -summary_df.index.name = "sample" -summary_df.sort_index(inplace=True) - -summary_df.to_csv(snakemake.output[0]) diff --git a/workflow/scripts/bin_summary_sample.py b/workflow/scripts/bin_summary_sample.py index 9876f3f..d7acfc9 100644 --- a/workflow/scripts/bin_summary_sample.py +++ b/workflow/scripts/bin_summary_sample.py @@ -4,27 +4,15 @@ sys.stderr = open(snakemake.log[0], "w") ## input files -in_dastool = ( - snakemake.input.tool -) -in_checkm = ( - snakemake.input.checkm -) -in_gtdb = ( - snakemake.input.gtdb -) +in_dastool = snakemake.input.tool +in_checkm = snakemake.input.checkm +in_gtdb = snakemake.input.gtdb ## output files ### csv -csv_path_mags = ( - snakemake.output.csv_mags -) -csv_path_bins = ( - snakemake.output.csv_bins -) -csv_path_tax = ( - snakemake.output.csv_tax -) +csv_path_mags = snakemake.output.csv_mags +csv_path_bins = snakemake.output.csv_bins +csv_path_tax = snakemake.output.csv_tax ## params max_cont = snakemake.params.max_cont @@ -133,23 +121,3 @@ def save_csv_table(csv_path, summary_df): ] mags_df.index.name = "MAG" save_csv_table(csv_path_mags, mags_df) - - -""" -if we want to sort by mean of deviation from min completeness & max contamination - -mean_dict={} -for binid in bins_df.index: - cont=bins_df.at[binid,"contamination"] - cont_per=round(((max_contamination - cont)/max_contamination), 4) - comp=bins_df.at[binid,"completeness"] - comp_per=round((comp/100),4) - mean_dict[binid]=round(((cont_per + comp_per)/2),4) - -to_sort=pd.DataFrame.from_dict(mean_dict,orient="index",columns=['mean']) - -# index list sorted by mean used to reindex original df -new_ind=to_sort.sort_values(["mean"],ascending = [False]).index.to_list() -sorted_all=bins_df.reindex(new_ind) - -""" From abebda491675d01fd709e44cb074381c7ea94ac3 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Mon, 23 Mar 2026 10:42:32 +0000 Subject: [PATCH 07/15] usage of local databases --- config/config.yaml | 5 ++++- workflow/rules/bin_qc.smk | 2 +- workflow/rules/common.smk | 11 +++++------ 3 files changed, 10 insertions(+), 8 deletions(-) diff --git a/config/config.yaml b/config/config.yaml index 04f8769..fa73f93 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -52,12 +52,15 @@ kaiju: checkm2: use-local: True + #if use-local False a folder named CheckM2_database will be created under the given path db-folder: /groups/ds/databases_refGenomes/databases/latest/checkM2/ # folder of genomad DB genomad: use-local: True - db-folder: /groups/ds/databases_refGenomes/databases/latest/genomad/ + # a folder named genomad_db is expected under the given path + # if use-local False this folder will be created + db-folder: /groups/ds/databases_refGenomes/databases/latest/ card: use-local: True diff --git a/workflow/rules/bin_qc.smk b/workflow/rules/bin_qc.smk index 54987b4..c8f1b79 100644 --- a/workflow/rules/bin_qc.smk +++ b/workflow/rules/bin_qc.smk @@ -5,7 +5,7 @@ if not config["checkm2"]["use-local"]: output: dbfile=get_checkm2_db(), params: - direct=get_checkm2_db_folder(), + direct=lambda wildcards, output: Path(output.dbfile).parent.parent, log: "logs/checkm2_DB_download.log", conda: diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 5995d50..097d363 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -117,12 +117,10 @@ def get_kaiju_files(): return files -def get_checkm2_db_folder(): - return config["checkm2"]["db-folder"] - - def get_checkm2_db(): - path = "".join([get_checkm2_db_folder(), "CheckM2_database/uniref100.KO.1.dmnd"]) + path = "".join( + [config["checkm2"]["db-folder"], "CheckM2_database/uniref100.KO.1.dmnd"] + ) return path @@ -131,7 +129,8 @@ def get_gtdb_folder(): def get_genomad_DB_folder(): - return config["genomad"]["db-folder"] + path = "".join([config["genomad"]["db-folder"], "genomad_db/"]) + return path def get_genomad_DB_file(): From d6dbe0ecbe786cff8ce2e387da0529592cf2b2e4 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Tue, 31 Mar 2026 12:39:34 +0000 Subject: [PATCH 08/15] add UniCARD usage for resistance determination --- config/config.yaml | 6 + workflow/Snakefile | 26 ++- workflow/envs/diamond.yaml | 6 + workflow/rules/analysis.smk | 54 ++--- workflow/rules/classify.smk | 114 +++++++++- workflow/rules/common.smk | 40 ++++ workflow/rules/qc.smk | 40 ++-- workflow/rules/report.smk | 6 +- workflow/rules/resistance.smk | 249 ++++++++++++++++++++++ workflow/scripts/contig_classification.py | 75 +++++++ workflow/scripts/uniCARD_filtering.py | 171 +++++++++++++++ 11 files changed, 722 insertions(+), 65 deletions(-) create mode 100644 workflow/envs/diamond.yaml create mode 100644 workflow/rules/resistance.smk create mode 100644 workflow/scripts/contig_classification.py create mode 100644 workflow/scripts/uniCARD_filtering.py diff --git a/config/config.yaml b/config/config.yaml index fa73f93..bf7dfc1 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -68,6 +68,12 @@ card: version: v4.0.1 url: https://card.mcmaster.ca/latest/data +uniCARD: + #path to the .dmnd file of UniCARD database + db-file: /groups/ds/databases_refGenomes/databases/latest/uniCARD/uniCARD.dmnd + # path to CARDs hierarchy file + hierarchy-json: /groups/ds/databases_refGenomes/databases/latest/uniCARD/CARD_hierarchy_v4.0.0.json + ## string term used for formatting output tables tablular-config: '/>github<\/a>/a \\t\t\t\n\t\t\t
  • \n\t\t\t\t' diff --git a/workflow/Snakefile b/workflow/Snakefile index 9ddbf62..d9992ae 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -27,6 +27,7 @@ include: "rules/das_tool.smk" include: "rules/bin_qc.smk" include: "rules/classify.smk" include: "rules/analysis.smk" +include: "rules/resistance.smk" include: "rules/report.smk" @@ -34,7 +35,7 @@ rule all: input: # read classification output expand( - "results/{project}/output/classification/reads/{sample}/{sample}_{level}_kaiju_summary.tsv", + "results/{project}/output/classification/reads/{sample}/{sample}_{level}_abundance.tsv", project=get_project(), sample=get_samples(), level=["species", "genus", "family", "order", "class", "phylum"], @@ -45,11 +46,6 @@ rule all: project=get_project(), sample=get_samples(), ), - expand( - "results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", - project=get_project(), - sample=get_samples(), - ), # fasta output expand( "results/{project}/output/fastas/{sample}/bins/", @@ -78,3 +74,21 @@ rule all: project=get_project(), sample=get_samples(), ), + expand( + "results/{project}/output/resistance/uniCARD/assembly/{sample}.csv", + project=get_project(), + sample=get_samples(), + ), + #contig classification + expand( + "results/{project}/output/classification/assembly/{sample}/{sample}_classified_contigs.csv", + project=get_project(), + sample=get_samples(), + ), + + +"""expand( + "results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", + project=get_project(), + sample=get_samples(), +),""" diff --git a/workflow/envs/diamond.yaml b/workflow/envs/diamond.yaml new file mode 100644 index 0000000..42d39a8 --- /dev/null +++ b/workflow/envs/diamond.yaml @@ -0,0 +1,6 @@ +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - diamond = 2.1.12 diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 660094a..1427fbc 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -8,9 +8,9 @@ rule gene_identification: gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff", log: "logs/{project}/proteins/{sample}.log", - threads: 32 conda: "../envs/pprodigal.yaml" + threads: 32 shell: "pprodigal -i {input.contigs} -o {output.gff} -a {output.faa} " "-d {output.fna} -p meta --tasks {threads} > {log} 2>&1" @@ -25,12 +25,12 @@ rule gzip_proteins: faa="results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", fna="results/{project}/output/proteins/{sample}/{sample}_nucleotides.fna.gz", gff="results/{project}/output/proteins/{sample}/{sample}_annotations.gff.gz", - threads: 20 - priority: 1 log: "logs/{project}/proteins/{sample}_gzip.log", + priority: 1 conda: "../envs/unix.yaml" + threads: 20 shell: "gzip {input.faa} {input.fna} {input.gff} > {log} 2>&1" @@ -43,12 +43,12 @@ if not config["genomad"]["use-local"]: output: folder=get_genomad_DB_folder(), file=get_genomad_DB_file(), - params: - res_folder=lambda wildcards, output: Path(output.folder).parent, log: "logs/load_genomad_DB.log", conda: "../envs/genomad.yaml" + params: + res_folder=lambda wildcards, output: Path(output.folder).parent, shell: "genomad download-database {params.res_folder}/ > {log} 2>&1" @@ -64,13 +64,13 @@ rule genomad_run: virus_tsv=temp( "results/{project}/genomad/{sample}/{sample}_summary/{sample}_virus_summary.tsv" ), - params: - outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent.parent, log: "logs/{project}/plasmids/{sample}_run.log", - threads: 64 conda: "../envs/genomad.yaml" + threads: 64 + params: + outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent.parent, shell: "genomad end-to-end --cleanup -t {threads} " "{input.asmbl} {params.outdir}/ " @@ -84,33 +84,34 @@ rule move_genomad_output: output: plasmid_tsv="results/{project}/output/plasmids/{sample}/{sample}_plasmid_summary.tsv", virus_tsv="results/{project}/output/plasmids/{sample}/{sample}_virus_summary.tsv", - params: - outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent, - sum_folder=lambda wildcards, input: Path(input.plasmid_tsv).parent, log: "logs/{project}/plasmids/{sample}_move_output.log", - threads: 32 conda: "../envs/unix.yaml" + threads: 32 + params: + outdir=lambda wildcards, output: Path(output.plasmid_tsv).parent, + sum_folder=lambda wildcards, input: Path(input.plasmid_tsv).parent, shell: "scp {params.sum_folder}/*.tsv {params.sum_folder}/*.json {params.outdir}/ > {log} 2>&1" +""" # Resistance analysis if not config["card"]["use-local"]: rule download_CARD_data: output: json=get_card_db_file(), - params: - download=config["card"]["url"], - folder=lambda wildcards, output: Path(output.json).parent, - filename=lambda wildcards, output: Path(output.json).name, log: "logs/CARD_data_download.log", - threads: 30 conda: "../envs/unix.yaml" + threads: 30 + params: + download=config["card"]["url"], + folder=lambda wildcards, output: Path(output.json).parent, + filename=lambda wildcards, output: Path(output.json).name, shell: "(cd {params.folder} && " "wget {params.download} && " @@ -127,8 +128,7 @@ rule CARD_load_DB_for_reads: conda: "../envs/card.yaml" shell: - "rgi clean --local && " - "rgi load --card_json {input} --local > {log} 2>&1" + "rgi clean --local && " "rgi load --card_json {input} --local > {log} 2>&1" rule CARD_annotation: @@ -136,15 +136,15 @@ rule CARD_annotation: json=get_card_db_file(), load=rules.CARD_load_DB_for_reads.output, output: - #done=temp(touch("results/CARD_annotation.done")), + done=temp(touch("results/CARD_annotation.done")), ann=get_card_annotation_file(), - params: - folder=lambda wildcards, input: Path(input.json).parent, - file=lambda wildcards, input: Path(input.json).name, log: "logs/CARD_annotation.log", conda: "../envs/card.yaml" + params: + folder=lambda wildcards, input: Path(input.json).parent, + file=lambda wildcards, input: Path(input.json).name, shell: "(cd {params.folder}/ && " "rgi card_annotation -i {params.file}) && " @@ -157,20 +157,20 @@ rule CARD_read_run: db=rules.CARD_annotation.output, output: txt="results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", - params: - folder=lambda wildcards, output: Path(output.txt).parent, log: "logs/{project}/ARGs/reads/{sample}.log", - threads: 64 conda: "../envs/card.yaml" + threads: 64 + params: + folder=lambda wildcards, output: Path(output.txt).parent, shell: "rgi bwt -1 {input.fa[0]} -2 {input.fa[1]} " "-o {params.folder}/{wildcards.sample} --local " "--clean -n {threads} > {log} 2>&1" -""" + rule CARD_read_sample_summary: input: txt=rules.CARD_read_run.output.txt, diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 471e5ae..c90b279 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -16,8 +16,8 @@ if not config["kaiju"]["use-local"]: "wget -c {params.download} -O - | " "tar -zxv -C {params.db_folder}) > {log} 2>&1" - -rule run_kaiju: +# classification of reads +rule run_kaiju_reads: input: db_files=get_kaiju_files(), fastqs=get_filtered_gz_fastqs, @@ -27,7 +27,7 @@ rule run_kaiju: evalue=0.00001, threads: 64 log: - "logs/{project}/kaiju/{sample}_run.log", + "logs/{project}/kaiju/{sample}_reads_run.log", conda: "../envs/kaiju.yaml" shell: @@ -36,15 +36,15 @@ rule run_kaiju: "-j {input.fastqs[1]} -o {output.report} > {log} 2>&1" -rule kaiju_table: +rule kaiju_table_reads: input: - report=rules.run_kaiju.output.report, + report=rules.run_kaiju_reads.output.report, db_files=get_kaiju_files(), output: - summary="results/{project}/output/classification/reads/{sample}/{sample}_{level}_kaiju_summary.tsv", + summary="results/{project}/output/classification/reads/{sample}/{sample}_{level}_abundance.tsv", threads: 2 log: - "logs/{project}/kaiju/{sample}_{level}_table.log", + "logs/{project}/kaiju/{sample}_reads_{level}_table.log", conda: "../envs/kaiju.yaml" shell: @@ -54,23 +54,23 @@ rule kaiju_table: # -p to print full taxonomy -rule kaiju2krona: +rule kaiju2krona_reads: input: db_files=get_kaiju_files(), - report=rules.run_kaiju.output.report, + report=rules.run_kaiju_reads.output.report, output: krona=temp( "results/{project}/output/classification/reads/{sample}/kaiju.out.krona" ), html=report( - "results/{project}/output/report/{sample}/{sample}_kaiju.out.html", + "results/{project}/output/report/{sample}/{sample}_reads_kaiju.out.html", htmlindex="index.html", category="5. Taxonomic classification", subcategory="5.2 Read classification", labels={"sample": "{sample}"}, ), log: - "logs/{project}/kaiju/{sample}_2krona.log", + "logs/{project}/kaiju/{sample}_reads_2krona.log", conda: "../envs/kaiju.yaml" shell: @@ -79,6 +79,98 @@ rule kaiju2krona: "ktImportText -o {output.html} {output.krona}) > {log} 2>&1" +# classification of contigs +rule run_kaiju_contigs: + input: + db_files=get_kaiju_files(), + assembly=get_gz_assembly, + output: + report=temp("results/{project}/output/classification/assembly/{sample}/kaiju.out"), + threads: 60 + log: + "logs/{project}/contig_classification/{sample}_kaiju_run.log", + conda: + "../envs/kaiju.yaml" + shell: + "kaiju -z {threads} -t {input.db_files[0]} -f {input.db_files[1]} " + "-i {input.assembly} -o {output.report} > {log} 2>&1" + + +rule gzip_kaiju_contigs: + input: + report=rules.run_kaiju_contigs.output.report, + output: + report="results/{project}/output/classification/assembly/{sample}/kaiju.out.gz", + threads: 20 + priority: 1 + log: + "logs/{project}/contig_classification/{sample}_kaiju_gzip.log", + conda: + "../envs/unix.yaml" + shell: + "gzip {input.report} > {log} 2>&1" + + +use rule kaiju_table_reads as kaiju_table_contigs with: + input: + report=rules.run_kaiju_contigs.output.report, + db_files=get_kaiju_files(), + output: + summary="results/{project}/output/classification/assembly/{sample}/{sample}_{level}_abundance.tsv", + threads: 2 + log: + "logs/{project}/contig_classification/{sample}_{level}_table.log", + + +use rule kaiju2krona_reads as kaiju2krona_contigs with: + input: + db_files=get_kaiju_files(), + report=rules.run_kaiju_contigs.output.report, + output: + krona=temp( + "results/{project}/output/classification/assembly/{sample}/kaiju.out.krona" + ), + html=report("results/{project}/output/report/{sample}/{sample}_contigs_kaiju.out.html", + htmlindex="index.html", + category="5. Taxonomic classification", + subcategory="5.3 Contig classification", + labels={"sample": "{sample}"}, + ), + log: + "logs/{project}/contig_classification/{sample}_2krona.log", + + +rule extract_contig_headers: + input: + assembly=get_gz_assembly, + output: + headers="results/{project}/output/classification/assembly/{sample}/headers.txt", + threads: 10 + log: + "logs/{project}/contig_classification/{sample}_header.log", + conda: + "../envs/unix.yaml" + shell: + "(zgrep '>' {input.assembly} > {output.headers}) > {log} 2>&1" + + +rule contig_classification: + input: + out=rules.run_kaiju_contigs.output.report, + genus="results/{project}/output/classification/assembly/{sample}/{sample}_genus_abundance.tsv", + species="results/{project}/output/classification/assembly/{sample}/{sample}_species_abundance.tsv", + headers=rules.extract_contig_headers.output.headers, + output: + csv="results/{project}/output/classification/assembly/{sample}/{sample}_classified_contigs.csv", + log: + "logs/{project}/contig_classification/{sample}_classified_contigs.log", + threads: 5 + conda: + "../envs/python.yaml" + script: + "../scripts/contig_classification.py" + + rule prepare_gtdb: output: done=temp(touch("results/GTDB_prep.done")), diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 097d363..cd38d08 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -75,6 +75,10 @@ def get_gz_assembly(wildcards): return "results/{project}/output/fastas/{sample}/{sample}.fa.gz" +def get_proteins(wildcards): + return ("results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz",) + + ## binning parameters def get_contig_length_threshold(): return config["min-contig-length"] @@ -149,7 +153,43 @@ def get_card_annotation_file(): return path +def get_uniCARD_db(): + return config["uniCARD"]["db-file"] + + +def get_CARD_hierarchy(): + return config["uniCARD"]["hierarchy-json"] + + +def get_uniCARD_db_wo_ext(): + db_wo_ext = str(Path(get_uniCARD_db()).with_suffix("")) + return db_wo_ext + + """ +def is_uniCARD_fasta(): + filename = str(Path(get_uniCARD_db()).name) + ext = "".join(Path(filename).suffixes) + if ext in ["fasta", "fasta.gz"]: + return True + else: + return False + + +def get_uniCARD_db_wo_ext(): + db_wo_ext = str(Path(get_uniCARD_db()).with_suffix("")) + if is_uniCARD_fasta(): + db_wo_ext = db_wo_ext.with_suffix("") + return db_wo_ext + + +def get_unicard_dmnd(): + if is_uniCARD_fasta(): + return "".join([get_uniCARD_db_wo_ext(), ".dmnd"]) + else: + return get_uniCARD_db() + + def get_mag_fa(wildcards): folder = "results/{}/output/fastas/{}/mags/".format( wildcards.project, wildcards.sample diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index 4f2d10d..e377b4b 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -13,15 +13,15 @@ rule fastp: ), html=temp("results/{project}/trimmed/fastp/{sample}.html"), json="results/{project}/output/report/prerequisites/qc/{sample}.fastp.json", + log: + "logs/{project}/fastp/{sample}.log", + threads: 16 params: adapters=get_adapters, extra="--qualified_quality_phred {phred} --length_required {minlen}".format( phred=(config["quality-criteria"]["min-PHRED"]), minlen=(config["quality-criteria"]["min-length-reads"]), ), - log: - "logs/{project}/fastp/{sample}.log", - threads: 16 wrapper: "v7.1.0/bio/fastp" @@ -30,19 +30,19 @@ rule fastqc: input: rules.fastp.output.trimmed, output: - zip=temp( - expand( - "results/{{project}}/trimmed/fastp/{{sample}}.{read}_fastqc.zip", - read=["1", "2"], - ) + zip1=temp( + "results/{project}/trimmed/fastp/{sample}.1_fastqc.zip", + ), + zip2=temp( + "results/{project}/trimmed/fastp/{sample}.2_fastqc.zip", ), - threads: 4 - resources: - mem_mb=1024, log: "logs/{project}/fastqc/{sample}.log", conda: "../envs/fastqc.yaml" + threads: 4 + resources: + mem_mb=1024, shell: "fastqc --memory {resources.mem_mb} --threads {threads} " "--format fastq --quiet {input} > {log} 2>&1" @@ -68,11 +68,11 @@ rule multiqc: labels={"sample": "all samples"}, ), "results/{project}/output/report/all/multiqc_{project}.zip", + log: + "logs/{project}/multiqc.log", params: extra=("--title 'Results for data from {project} project'"), use_input_files_only=True, - log: - "logs/{project}/multiqc.log", wrapper: "v8.1.1/bio/multiqc" @@ -92,13 +92,13 @@ rule qc_summary: output: csv="results/{project}/output/report/all/quality_summary.csv", vis_csv=temp("results/{project}/output/report/all/quality_summary_visual.csv"), - params: - other_host=config["host-filtering"]["do-host-filtering"], - hostname=config["host-filtering"]["host-name"], log: "logs/{project}/report/qc_summary.log", conda: "../envs/python.yaml" + params: + other_host=config["host-filtering"]["do-host-filtering"], + hostname=config["host-filtering"]["host-name"], script: "../scripts/quality_summary.py" @@ -117,16 +117,16 @@ rule qc_summary_report: }, ) ), + log: + "logs/{project}/report/qc_summary_rbt_csv.log", + conda: + "../envs/rbt.yaml" params: pin_until="sample", styles="resources/report/tables/", name="quality_summary", header="Quality summary based on fastp report and mapping to host genome(s)", pattern=config["tablular-config"], - log: - "logs/{project}/report/qc_summary_rbt_csv.log", - conda: - "../envs/rbt.yaml" shell: "rbt csv-report {input} --pin-until {params.pin_until} {output} && " "(sed -i '{params.pattern} {params.header}' " diff --git a/workflow/rules/report.smk b/workflow/rules/report.smk index e8c401e..421644f 100644 --- a/workflow/rules/report.smk +++ b/workflow/rules/report.smk @@ -26,7 +26,11 @@ rule snakemake_report: sample=get_samples(), ), expand( - "results/{{project}}/output/report/{sample}/{sample}_kaiju.out.html", + "results/{{project}}/output/report/{sample}/{sample}_reads_kaiju.out.html", + sample=get_samples(), + ), + expand( + "results/{{project}}/output/report/{sample}/{sample}_contigs_kaiju.out.html", sample=get_samples(), ), output: diff --git a/workflow/rules/resistance.smk b/workflow/rules/resistance.smk new file mode 100644 index 0000000..d049578 --- /dev/null +++ b/workflow/rules/resistance.smk @@ -0,0 +1,249 @@ +# Resistance analysis +if not config["card"]["use-local"]: + + rule download_CARD_data: + output: + json=get_card_db_file(), + log: + "logs/CARD_data_download.log", + conda: + "../envs/unix.yaml" + threads: 30 + params: + download=config["card"]["url"], + folder=lambda wildcards, output: Path(output.json).parent, + filename=lambda wildcards, output: Path(output.json).name, + shell: + "(cd {params.folder} && " + "wget {params.download} && " + "tar -xvf data ./{params.filename}) > {log} 2>&1" + + +rule CARD_load_DB_for_reads: + input: + get_card_db_file(), + output: + temp(touch("results/CARD_load_DB_for_reads.done")), + log: + "logs/CARD_load_DB_for_reads.log", + conda: + "../envs/card.yaml" + shell: + "rgi clean --local && " "rgi load --card_json {input} --local > {log} 2>&1" + + +rule CARD_annotation: + input: + json=get_card_db_file(), + load=rules.CARD_load_DB_for_reads.output, + output: + done=temp(touch("results/CARD_annotation.done")), + ann=get_card_annotation_file(), + log: + "logs/CARD_annotation.log", + conda: + "../envs/card.yaml" + params: + folder=lambda wildcards, input: Path(input.json).parent, + file=lambda wildcards, input: Path(input.json).name, + shell: + "(cd {params.folder}/ && " + "rgi card_annotation -i {params.file}) && " + "rgi load -i {input.json} --card_annotation {output.ann} --local > {log} 2>&1" + + +rule CARD_read_run: + input: + fa=get_filtered_gz_fastqs, + db=rules.CARD_annotation.output, + output: + txt="results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", + log: + "logs/{project}/ARGs/reads/{sample}.log", + conda: + "../envs/card.yaml" + threads: 64 + params: + folder=lambda wildcards, output: Path(output.txt).parent, + shell: + "rgi bwt -1 {input.fa[0]} -2 {input.fa[1]} " + "-o {params.folder}/{wildcards.sample} --local " + "--clean -n {threads} > {log} 2>&1" + + +""" +if is_uniCARD_fasta: + + rule uniCARD_makeDB: + input: + db=get_uniCARD_db(), + output: + dmdb=get_unicard_dmnd(), + log: + "logs/unicard_makeDB.log", + conda: + "../envs/diamond.yaml" + threads: 60 + params: + filename=lambda wildcards, input: Path(input.db).name, + filename_wo_ext=Path(get_uniCARD_db_wo_ext()).name, + folder=lambda wildcards, input: Path(input.db).parent, + shell: + "(cd {params.folder} && diamond makedb --in {params.filename} -d {params.filename_wo_ext}) > {log} 2>&1" +""" + + +rule uniCARD_run: + input: + faa=get_proteins, + dmdb=get_uniCARD_db(), + output: + tsv=temp( + "results/{project}/output/resistance/uniCARD/assembly/{sample}_out.tsv" + ), + log: + "logs/{project}/uniCARD/{sample}.log", + # ensure priority over potential contig classification + priority: 1 + conda: + "../envs/diamond.yaml" + threads: 60 + # to ensure only one of these commands are run at the same time + resources: + heavy=1, + params: + db_wo_ext=get_uniCARD_db_wo_ext(), + max_targets=15, + shell: + "diamond blastp -q {input.faa} -d {params.db_wo_ext} " + "-p {threads} --max-target-seqs {params.max_targets} " + "-o {output.tsv} > {log} 2>&1" + + +rule gz_uniCARD_tsv: + input: + unicard_in=rules.uniCARD_run.output.tsv, + output: + tsv=temp( + "results/{project}/output/resistance/uniCARD/assembly/{sample}_out.tsv.gz" + ), + log: + "logs/{project}/uniCARD/{sample}_gz.log", + conda: + "../envs/unix.yaml" + threads: 16 + shell: + "(gzip -k {input.unicard_in}) > {log} 2>&1" + + +rule filter_uniCARD: + input: + tsv=rules.gz_uniCARD_tsv.output.tsv, + json=get_CARD_hierarchy(), + output: + best="results/{project}/output/resistance/uniCARD/assembly/{sample}.csv", + log: + "logs/{project}/uniCARD/{sample}_filter.log", + conda: + "../envs/python.yaml" + threads: 20 + params: + # filter window for filtering UniCARD results, default: 1% + filter_window=0.01, + shell: + "python workflow/scripts/uniCARD_filtering.py --infile {input.tsv} " + "--card_hierarchy {input.json} --outfile {output.best} " + "--filter_window {params.filter_window} > {log} 2>&1" + + +""" +rule CARD_read_sample_summary: + input: + txt=rules.CARD_read_run.output.txt, + output: + csv="results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", + params: + case="reads", + log: + "logs/{project}/ARGs/reads/{sample}.log", + threads: 2 + conda: + "../envs/python.yaml" + script: + "../scripts/arg_summary_sample.py" + + + +# updates CARD database to use for contigs instead of reads +# read based classification must be finished before +rule CARD_load_DB: + input: + db=get_card_db_file(), + read_args=expand( + "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", + sample=get_samples(), + project=get_project(), + ), + output: + touch("results/CARD_load_DB.done"), + log: + "logs/CARD_load_DB.log", + conda: + "../envs/card.yaml" + shell: + "rgi clean --local && " + "rgi load --card_json {input.db} --local > {log} 2>&1" + + +rule CARD_assembly_run: + input: + faa=rules.gzip_proteins.output.faa, + db=rules.CARD_load_DB.output, + output: + txt="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.txt", + json="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.json", + params: + path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), + log: + "logs/{project}/ARGs/assembly/{sample}.log", + threads: 64 + conda: + "../envs/card.yaml" + shell: + "rgi main -i {input.faa} -o {params.path_wo_ext} " + "-t protein -a DIAMOND --local " + "-n {threads} --clean > {log} 2>&1" + + +use rule CARD_read_sample_summary as CARD_assembly_sample_summary with: + input: + txt=rules.CARD_assembly_run.output.txt, + json=rules.CARD_assembly_run.output.json, + output: + csv="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}_assembly_ARGs.csv", + params: + case="assembly", + log: + "logs/{project}/ARGs/assembly/{sample}.log", + + +use rule CARD_assembly_run as CARD_mag_run with: + input: + fa="results/{project}/output/fastas/{sample}/mags/{binID}.fa.gz", + db=rules.CARD_load_DB.output, + output: + txt="results/{project}/output/resistance/CARD/mags/{sample}/{binID}.txt", + params: + path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), + log: + "logs/{project}/ARGs/mags/{sample}/{binID}.log", + + +rule wrap_mag_ARGs: + input: + get_mag_ARGs, + output: + touch("results/{project}/output/ARGs/mags/{sample}/all_mags.done"), + log: + "logs/{project}/ARGs/mags/{sample}/all_mags.log", +""" diff --git a/workflow/scripts/contig_classification.py b/workflow/scripts/contig_classification.py new file mode 100644 index 0000000..6f3611e --- /dev/null +++ b/workflow/scripts/contig_classification.py @@ -0,0 +1,75 @@ +import pandas as pd +import re +import sys + +sys.stderr = open(snakemake.log[0], "w") + +kaiju_out = snakemake.input.out +kaiju_species = snakemake.input.species +kaiju_genus = snakemake.input.genus +headers = snakemake.input.headers + +classification = snakemake.output.csv + + +## Prepare classification and length information for all contigs classified to at least genus level ## + +df_kaiju_raw = pd.read_table(kaiju_out, names=["classified", "contig", "taxid"]) +kaiju_df = df_kaiju_raw[df_kaiju_raw["classified"] != "U"].copy() + +species_df = pd.read_csv(kaiju_species, sep="\t") +# remove rows for unclassified and virus +species_df = species_df.iloc[:-3].copy() + +genus_df = pd.read_csv(kaiju_genus, sep="\t") +# remove rows for unclassified and virus +genus_df = genus_df.iloc[:-3].copy() + +# make taxon name to taxon id dictionary for all occurring species and genera +species_map = dict(zip(species_df["taxon_id"], species_df["taxon_name"])) +genus_map = dict(zip(genus_df["taxon_id"], genus_df["taxon_name"])) + +# kaiju_df is the original DataFrame with only classified contigs +kaiju_df["species_name"] = kaiju_df["taxid"].map(species_map) +kaiju_df["genus_name"] = kaiju_df["taxid"].map(genus_map) + +# removing all not classified to species or genus level +kaiju_df_red = kaiju_df[ + ~(kaiju_df["species_name"].isna() & kaiju_df["genus_name"].isna()) +].copy() +# adding info about the level of classification +kaiju_df_red["level"] = ( + kaiju_df_red["species_name"].notna().map({True: "species", False: "genus"}) +) +kaiju_df_red["classification"] = kaiju_df_red["species_name"].combine_first( + kaiju_df_red["genus_name"] +) +kaiju_df_red = kaiju_df_red.drop( + columns=["species_name", "genus_name", "classified", "taxid"] +) + +# get contig length +records = [] + +with open(headers) as f: + for line in f: + # Example header: >k141_570066 flag=0 multi=1.0000 len=450 + contig_match = re.search(r"^>(\S+)", line) + length_match = re.search(r"len=(\d+)", line) + + if contig_match and length_match: + contig = contig_match.group(1) + length = int(length_match.group(1)) + records.append((contig, length)) + +# Convert to DataFrame +length_df = pd.DataFrame(records, columns=["contig", "contig_len"]) + +# merge contig classification with length information +merged_df = kaiju_df_red.merge(length_df, on="contig", how="left") +cols = merged_df.columns.tolist() +cols.insert(1, cols.pop(cols.index("contig_len"))) +merged_df = merged_df[cols] + +# Write out contig classification file +merged_df.to_csv(classification, index=False) diff --git a/workflow/scripts/uniCARD_filtering.py b/workflow/scripts/uniCARD_filtering.py new file mode 100644 index 0000000..0e62acd --- /dev/null +++ b/workflow/scripts/uniCARD_filtering.py @@ -0,0 +1,171 @@ +import json +import pandas as pd +import argparse + +# identifier for CARD derived sequences +WT_PREFIX = "CARD_WT:" +CARD_PREFIX = "CARD:" + +# columns to reduce DIAMOND output +DIAMOND_COLS = [ + "query", + "target", + "%identity", + "length", + "#mismatches", + "#gaps", + "query_start", + "query_end", + "target_start", + "target_end", + "evalue", + "bitscore", +] + + +def write_filtered_ARGs(card_classes, unicard_in, outfile, filter_window): + + try: + # read in CARD drugs and classes json file + with open(card_classes) as f: + data = json.load(f) + + card_df = pd.DataFrame.from_dict(data, orient="index") + card_df.index.name = "ARO_ID" + card_df = card_df.reset_index() + + if unicard_in.endswith(".gz"): + df_unicard = pd.read_table( + unicard_in, names=DIAMOND_COLS, compression="gzip" + ) + else: + df_unicard = pd.read_table(unicard_in, names=DIAMOND_COLS) + + ## top 10 ## + # Step 1: Keep only 10 best hits per query + df_top10 = df_unicard[df_unicard.groupby("query").cumcount() < 10] + + # Step 2: Keep all rows for a query **only if at least one Target starts with CARD_PREFIX, + # so wild-type CARD entries are not kept (starting with WT_PREFIX) + df_top10 = df_top10[ + df_top10["query"].isin( + df_top10[df_top10["target"].str.startswith(CARD_PREFIX)]["query"] + ) + ] + + ## filter besthit CARD entries + + # filter for the first row per query + df_firsts = df_top10[df_top10.groupby("query").cumcount() < 1] + # 1. check if best hit is CARD wild-type -> remove from df + df_firsts_reduced = df_firsts[ + ~df_firsts["query"].isin( + df_firsts[df_firsts["target"].str.startswith(WT_PREFIX)]["query"] + ) + ] + # 2. from remaining save those, which are CARD hits + df_bests = df_firsts_reduced[ + df_firsts_reduced["query"].isin( + df_firsts_reduced[ + df_firsts_reduced["target"].str.startswith(CARD_PREFIX) + ]["query"] + ) + ] + + # remove first hit CARD queries from first row per query df + df_firsts_remaining = df_firsts_reduced[ + ~df_firsts_reduced["query"].isin(df_bests["query"]) + ] + + df_filtered = df_top10[df_top10["query"].isin(df_firsts_remaining["query"])] + + filtered_card_rows = [] + # go through all remaining best hits and compare to CARD hits + for idx, row in df_firsts_remaining.iterrows(): + # get all rows with same query to compare + compare_df = df_filtered[df_filtered["query"] == row["query"]] + compare_df = compare_df[ + compare_df["target"].str.startswith(CARD_PREFIX) + & ~compare_df["target"].str.contains("_snp_") + ] + # check for each row in this comparison df if the CARD entry has same values than first match + for idx_card, row_card in compare_df.iterrows(): + if ( + row["%identity"] == row_card["%identity"] + and row["bitscore"] == row_card["bitscore"] + ): + filtered_card_rows.append(row_card) + break + # additional filtering in default 1% window below best hit + elif (row["%identity"] * (1 - filter_window)) <= row_card[ + "%identity" + ] and (row["bitscore"] * (1 - filter_window)) <= row_card["bitscore"]: + filtered_card_rows.append(row_card) + break + + filtered_card_df = pd.DataFrame(filtered_card_rows) + + df_bests = pd.concat([df_bests, filtered_card_df]) + + # Adding contig column + df_bests.insert( + 1, "contig", df_top10["query"].map(lambda x: x.rsplit("_", 1)[0]) + ) + # extract ARO_ID, for variants it's eg 'CARD:3003394_snp_R123S' + df_bests.insert( + 3, + "ARO_ID", + df_top10["target"].map( + lambda x: x.split(CARD_PREFIX, 1)[-1].split("_", 1)[0] + ), + ) + + # Adding CARD drug info: + merged_df = pd.merge(df_bests, card_df, on="ARO_ID", how="left") + + # Write besthits to output file + merged_df.to_csv(outfile, index=False) + + except Exception as e: + print(f"Error filtering UniCARD results: {e}") + raise + + +def main(): + parser = argparse.ArgumentParser(description="Filter UniCARD results") + + parser.add_argument( + "--infile", + required=True, + help="DIAMOND output tsv file with UniCARD results to filter, may be gz", + ) + parser.add_argument( + "--card_hierarchy", required=True, help="Path to CARDs hierarchy json file" + ) + parser.add_argument( + "--outfile", + required=True, + help="Output csv file to store final filtered ARGs", + ) + parser.add_argument( + "--filter_window", + type=float, + default=0.01, + help="Percent window for filtering (default 0.01 (1%%))", + ) + + args = parser.parse_args() + + infile = args.infile + outfile = args.outfile + card_hierarchy = args.card_hierarchy + filter_window = args.filter_window + + print("Filtering UniCARD diamond results...") + write_filtered_ARGs(card_hierarchy, infile, outfile, filter_window) + + print(f"Filtering done! Results saved to: {outfile}") + + +if __name__ == "__main__": + main() From 4e26a60a47814a46456745a87092bc28c1073f53 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Wed, 1 Apr 2026 14:43:13 +0000 Subject: [PATCH 09/15] fix metabat processing bug --- workflow/Snakefile | 12 +++++------- workflow/scripts/metabat_postprocess.py | 25 +++++++++++++++---------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/workflow/Snakefile b/workflow/Snakefile index d9992ae..009e616 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -85,10 +85,8 @@ rule all: project=get_project(), sample=get_samples(), ), - - -"""expand( - "results/{project}/output/proteins/{sample}/{sample}_proteins.faa.gz", - project=get_project(), - sample=get_samples(), -),""" + expand( + "results/{project}/output/classification/assembly/{sample}/kaiju.out.gz", + project=get_project(), + sample=get_samples(), + ), diff --git a/workflow/scripts/metabat_postprocess.py b/workflow/scripts/metabat_postprocess.py index a38d269..eac3999 100644 --- a/workflow/scripts/metabat_postprocess.py +++ b/workflow/scripts/metabat_postprocess.py @@ -42,14 +42,19 @@ def append_string_to_file(filepath, append_string, mod_ext): for filename in os.listdir(directory) if filename.endswith(mod_ext) ] - - print("starting to concatenate files") - for modfile in modfiles: - # Concatenate all modified files into a single file - cmd = f"cat {modfile} >> {output_tsv}" + if len(modfiles) == 0: + print("no metabat bins") + # create an empty file + cmd = f"touch {output_tsv}" + os.system(cmd) + else: + print("starting to concatenate files") + for modfile in modfiles: + # Concatenate all modified files into a single file + cmd = f"cat {modfile} >> {output_tsv}" + os.system(cmd) + + print("all files concatenated") + # remove all modified files after concatenating + cmd = f"rm {directory}/*{mod_ext}" os.system(cmd) - - print("all files concatenated") - # remove all modified files after concatenating - cmd = f"rm {directory}/*{mod_ext}" - os.system(cmd) From ba00e07f3abdb5972f325d4d9d160492c481d56f Mon Sep 17 00:00:00 2001 From: josefawelling Date: Thu, 9 Apr 2026 14:49:37 +0000 Subject: [PATCH 10/15] add resistance abundance on contigs output, formatting --- README.md | 133 ++++++--------- workflow/Snakefile | 8 +- workflow/rules/analysis.smk | 167 ------------------- workflow/rules/classify.smk | 4 +- workflow/rules/common.smk | 4 +- workflow/rules/resistance.smk | 49 ++++-- workflow/scripts/assembly_summary.py | 57 ++++--- workflow/scripts/uniCARD_abundance.py | 52 ++++++ workflow/scripts/uniCARD_abundance_sample.py | 95 +++++++++++ 9 files changed, 275 insertions(+), 294 deletions(-) create mode 100644 workflow/scripts/uniCARD_abundance.py create mode 100644 workflow/scripts/uniCARD_abundance_sample.py diff --git a/README.md b/README.md index 730ba36..cd87352 100644 --- a/README.md +++ b/README.md @@ -4,90 +4,37 @@ [![GitHub actions status](https://github.com///workflows/Tests/badge.svg?branch=main)](https://github.com///actions?query=branch%3Amain+workflow%3ATests) -ResMAG is a state-of-the-art and user-friendly Snakemake workflow designed for the analysis of metagenomic data. It integrates multiple bioinformatics tools and algorithms to facilitate key steps in metagenome analysis, including bin refinement, metagenome-assembled genome (MAG) reconstruction, taxonomic classification of MAGs, and identification of antibiotic resistance genes.
    +ResMAG is a state-of-the-art and user-friendly Snakemake workflow designed for the analysis of metagenomic data. It integrates multiple bioinformatics tools and algorithms to facilitate key steps in metagenome analysis, including quality control, assembling, bin refinement, metagenome-assembled genome (MAG) reconstruction and taxonomic classification. ResMAG has a special focus on highly diverse samples, such as wastewater, and the identification of antibiotic resistance genes.
    -### Key Features +--- -**Binning Techniques**: Employ a collection of five state-of-the-art binning tools to partition metagenomic contigs into individual bins, allowing for comprehensive and accurate analysis.
    +## Key Features -**MAG Reconstruction**: Utilize cutting-edge algorithms to reconstruct high-quality metagenome-assembled genomes (MAGs) from sequencing data.
    +**Binning Techniques**: Employ a collection of two state-of-the-art binning tools to partition metagenomic contigs into individual bins, allowing for comprehensive and accurate analysis.
    -**Taxonomic Classification**: Apply advanced taxonomic classification methods to assign taxonomic labels to MAGs and identify the microbial community composition within the metagenomic samples.
    +**MAG Reconstruction**: Utilize cutting-edge algorithms to reconstruct high-quality metagenome-assembled genomes (MAGs), especially from highly diverse microbial communities, like wastewater.
    + +**Taxonomic Classification**: Apply advanced taxonomic classification methods to assign taxonomic labels to reads, contigs and MAGs and identify the microbial community composition within the metagenomic samples.
    **Antibiotic Resistance Gene Identification**: Perform in-depth analysis to detect and characterize antibiotic resistance genes within the metagenomic data, providing valuable insights into antimicrobial resistance profiles.
    **Performance Refinement**: Continuously optimize the pipeline by incorporating the latest advancements in metagenomics research, ensuring the highest accuracy and efficiency in metagenomic data analysis.
    -### Overview -```mermaid -%%{init: { - 'theme':'base', - 'themeVariables': { - 'secondaryColor': '#fff', - 'tertiaryColor': '#fff', - 'tertiaryBorderColor' : '#fff'} - }}%% - -flowchart TB; - - subgraph " " - direction TB - - %% Nodes - A[/short reads/] - B["QC
    fastp and fastQC"] - C["Filtering human reads
    minimap2"] - L["optional
    Filtering host reads
    minimap2"] - N["taxonomic classification
    Kaiju"] - D["Assembly
    MegaHIT"] - E["Binning
    MetaBAT and MetaCoAG"] - F["Bin refinement
    DAS Tool"] - G["Bin quality control
    CheckM2"] - M[/MAGs/] - H["Resistance analysis
    CARD"] - I["taxonomic classification
    GTDBTk"] - J[/MultiQC report/] - K[/Assembly summary/] - - %% input & output node design - classDef in_output fill:#fff,stroke:#cde498,stroke-width:4px - class A,J,K,M in_output - %% rule node design - classDef rule fill:#cde498,stroke:#000 - class B,C,D,E,F,G,H,I,L,N rule - - %% Node links - A --> B - B --> C - B --- J - C --> L - L --> D - C --> N - L --> N - C --> D - D --> E - D --- K - E --> F - F --> G - G --"filtering"--> M - M --- H - M --- I - - end +--- -``` +## Overview -## Usage +--- -### Preparations +## Preparations -#### Clone workflow -To prepare the workflow clone it to your desired working folder via git or your preferred IDE -``` -git clone https://github.com/IKIM-Essen/metagenomics_res.git -``` +Create a snakemake environment using [mamba](https://mamba.readthedocs.io/en/latest/) via: -#### Download GTDB + ```mamba create -c conda-forge -c bioconda -n snakemake snakemake snakemake-storage-plugin-fs``` + +For installation details, see the [instructions in the Snakemake documentation](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). + +### Download GTDB The GTDB needs to be downloaded and decompressed, it requires about 140 Gb. 1. Change to the directory where the GTDB should be stored @@ -101,17 +48,18 @@ The GTDB needs to be downloaded and decompressed, it requires about 140 Gb. tar xzf gtdbtk_data.tar.gz ``` 4. After successful step 3: the archive can be removed -5. Please specify the path to your decompressed GTDB in the config file (see [Configuring workflow](#configuring-workflow)) - +5. Please specify the path to your decompressed GTDB in the config file (see [Configuring workflow](#configuring-the-workflow)) -#### Install Snakemake -Create a snakemake environment using [mamba](https://mamba.readthedocs.io/en/latest/) via: +--- - ```mamba create -c conda-forge -c bioconda -n snakemake snakemake snakemake-storage-plugin-fs``` +## Usage -For installation details, see the [instructions in the Snakemake documentation](https://snakemake.readthedocs.io/en/stable/getting_started/installation.html). +Please obtain a copy of this workflow by [cloning](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) this repository. +``` +git clone https://github.com/IKIM-Essen/metagenomics_res.git +``` -#### Configuring workflow +### Configuring the workflow 1. Edit the `config/config.yaml` file: - Specify a project name (`project-name`) - Specify filtering options for human reads (`human-filtering`) @@ -139,43 +87,62 @@ Executing the workflow: using `$N` cores. It is recommended to use all available cores. -The usage of this workflow is described in the [Snakemake Workflow Catalog](https://snakemake.github.io/snakemake-workflow-catalog/?usage=%2F). +--- ## Output +--- + ## Tools A list of tools used in the pipeline: | Tool | Link | | --------- | ------------------------------------------------- | +| CARD | https://doi.org/10.1093/nar&gkac920 | | CheckM2 | https://doi.org/10.1038/s41592-023-01940-w | | CoverM | https://github.com/wwood/CoverM | | DAS Tool | https://doi.org/10.1038/s41564-018-0171-1 | +| DIAMOND | https://doi.org/10.1038/s41592-021-01101-x | | fastp | https://doi.org/10.1093/bioinformatics/bty560 | | FastQC | www.bioinformatics.babraham.ac.uk/projects/fastqc | +| geNomad | https://doi.org/10.1038/s41587-023-01953-y | +| GTDB-Tk | https://doi.org/10.1093/bioinformatics/btz848 | | Kaiju | https://doi.org/10.1038/ncomms11257 | | MEGAHIT | https://doi.org/10.1093/bioinformatics/btv033 | | MetaBAT | http://dx.doi.org/10.7717/peerj.1165 | | MetaCoAG | https://doi.org/10.1101/2021.09.10.459728 | | minimap2 | https://doi.org/10.1093/bioinformatics/bty191 | | MultiQC | www.doi.org/10.1093/bioinformatics/btw354 | -| pandas | pandas.pydata.org | +| pprodigal | https://github.com/sjaenick/pprodigal | +| Rust-Bio | https://doi.org/10.1093/bioinformatics/btv573 | | samtools | https://doi.org/10.1093/gigascience/giab008 | | Snakemake | www.doi.org/10.12688/f1000research.29032.1 | +| UniCARD | https://github.com/IKIM-Essen/uniCARD | + +--- + +## Contributions -## Contact Information +Pull requests and feature suggestions are very welcome! Feel free to fork and submit improvements. +For any other questions, or feedback, please contact the project maintainer Josefa Welling (@josefawelling) at josefa.welling@uk-essen.de. -For any questions, or feedback, please contact the project maintainer at josefa.welling@uk-essen.de. We appreciate your input and support in using and improving ResMAG. +We appreciate your input and support in using and improving ResMAG. + +--- ## Acknowledgements We would like to express our gratitude towards Katharina Block, Adrian Doerr, Miriam Balzer, Alexander Thomas, Johannes Köster, Ann-Kathrin Doerr and the IKIM who have contributed to the development and testing of ResMAG. Their valuable insights and feedback have been helpful throughout the creation of the workflow. +--- + ## Citation -A paper is on its way. If you use ResMAG in your work before the paper, then please consider citing this GitHub. +A paper is on its way. If you use ResMAG in your work, don't forget to give credits to the authors by citing the URL of this repository. + +--- ## License -ResMAG is released under the [BSD-2 Clause](https://www.open-xchange.com/hubfs/2_Clause_BSD_License.pdf?hsLang=en). Please review the license file for more details. \ No newline at end of file +ResMAG is licensed under the [BSD-2 Clause](https://www.open-xchange.com/hubfs/2_Clause_BSD_License.pdf?hsLang=en). \ No newline at end of file diff --git a/workflow/Snakefile b/workflow/Snakefile index 009e616..3d0b47e 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -75,16 +75,10 @@ rule all: sample=get_samples(), ), expand( - "results/{project}/output/resistance/uniCARD/assembly/{sample}.csv", + "results/{project}/output/resistance/uniCARD/{project}_drug_class_abundance_all_possible.csv", project=get_project(), - sample=get_samples(), ), #contig classification - expand( - "results/{project}/output/classification/assembly/{sample}/{sample}_classified_contigs.csv", - project=get_project(), - sample=get_samples(), - ), expand( "results/{project}/output/classification/assembly/{sample}/kaiju.out.gz", project=get_project(), diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 1427fbc..196df26 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -94,170 +94,3 @@ rule move_genomad_output: sum_folder=lambda wildcards, input: Path(input.plasmid_tsv).parent, shell: "scp {params.sum_folder}/*.tsv {params.sum_folder}/*.json {params.outdir}/ > {log} 2>&1" - - -""" -# Resistance analysis -if not config["card"]["use-local"]: - - rule download_CARD_data: - output: - json=get_card_db_file(), - log: - "logs/CARD_data_download.log", - conda: - "../envs/unix.yaml" - threads: 30 - params: - download=config["card"]["url"], - folder=lambda wildcards, output: Path(output.json).parent, - filename=lambda wildcards, output: Path(output.json).name, - shell: - "(cd {params.folder} && " - "wget {params.download} && " - "tar -xvf data ./{params.filename}) > {log} 2>&1" - - -rule CARD_load_DB_for_reads: - input: - get_card_db_file(), - output: - temp(touch("results/CARD_load_DB_for_reads.done")), - log: - "logs/CARD_load_DB_for_reads.log", - conda: - "../envs/card.yaml" - shell: - "rgi clean --local && " "rgi load --card_json {input} --local > {log} 2>&1" - - -rule CARD_annotation: - input: - json=get_card_db_file(), - load=rules.CARD_load_DB_for_reads.output, - output: - done=temp(touch("results/CARD_annotation.done")), - ann=get_card_annotation_file(), - log: - "logs/CARD_annotation.log", - conda: - "../envs/card.yaml" - params: - folder=lambda wildcards, input: Path(input.json).parent, - file=lambda wildcards, input: Path(input.json).name, - shell: - "(cd {params.folder}/ && " - "rgi card_annotation -i {params.file}) && " - "rgi load -i {input.json} --card_annotation {output.ann} --local > {log} 2>&1" - - -rule CARD_read_run: - input: - fa=get_filtered_gz_fastqs, - db=rules.CARD_annotation.output, - output: - txt="results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", - log: - "logs/{project}/ARGs/reads/{sample}.log", - conda: - "../envs/card.yaml" - threads: 64 - params: - folder=lambda wildcards, output: Path(output.txt).parent, - shell: - "rgi bwt -1 {input.fa[0]} -2 {input.fa[1]} " - "-o {params.folder}/{wildcards.sample} --local " - "--clean -n {threads} > {log} 2>&1" - - - -rule CARD_read_sample_summary: - input: - txt=rules.CARD_read_run.output.txt, - output: - csv="results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", - params: - case="reads", - log: - "logs/{project}/ARGs/reads/{sample}.log", - threads: 2 - conda: - "../envs/python.yaml" - script: - "../scripts/arg_summary_sample.py" - - - -# updates CARD database to use for contigs instead of reads -# read based classification must be finished before -rule CARD_load_DB: - input: - db=get_card_db_file(), - read_args=expand( - "results/{project}/output/resistance/CARD/reads/{sample}/{sample}.gene_mapping_data.txt", - sample=get_samples(), - project=get_project(), - ), - output: - touch("results/CARD_load_DB.done"), - log: - "logs/CARD_load_DB.log", - conda: - "../envs/card.yaml" - shell: - "rgi clean --local && " - "rgi load --card_json {input.db} --local > {log} 2>&1" - - -rule CARD_assembly_run: - input: - faa=rules.gzip_proteins.output.faa, - db=rules.CARD_load_DB.output, - output: - txt="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.txt", - json="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}.json", - params: - path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), - log: - "logs/{project}/ARGs/assembly/{sample}.log", - threads: 64 - conda: - "../envs/card.yaml" - shell: - "rgi main -i {input.faa} -o {params.path_wo_ext} " - "-t protein -a DIAMOND --local " - "-n {threads} --clean > {log} 2>&1" - - -use rule CARD_read_sample_summary as CARD_assembly_sample_summary with: - input: - txt=rules.CARD_assembly_run.output.txt, - json=rules.CARD_assembly_run.output.json, - output: - csv="results/{project}/output/resistance/CARD/assembly/{sample}/{sample}_assembly_ARGs.csv", - params: - case="assembly", - log: - "logs/{project}/ARGs/assembly/{sample}.log", - - -use rule CARD_assembly_run as CARD_mag_run with: - input: - fa="results/{project}/output/fastas/{sample}/mags/{binID}.fa.gz", - db=rules.CARD_load_DB.output, - output: - txt="results/{project}/output/resistance/CARD/mags/{sample}/{binID}.txt", - params: - path_wo_ext=lambda wildcards, output: Path(output.txt).with_suffix(""), - log: - "logs/{project}/ARGs/mags/{sample}/{binID}.log", - - -rule wrap_mag_ARGs: - input: - get_mag_ARGs, - output: - touch("results/{project}/output/ARGs/mags/{sample}/all_mags.done"), - log: - "logs/{project}/ARGs/mags/{sample}/all_mags.log", -""" diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index c90b279..5a40f69 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -108,7 +108,7 @@ rule gzip_kaiju_contigs: conda: "../envs/unix.yaml" shell: - "gzip {input.report} > {log} 2>&1" + "gzip -k {input.report} > {log} 2>&1" use rule kaiju_table_reads as kaiju_table_contigs with: @@ -124,8 +124,8 @@ use rule kaiju_table_reads as kaiju_table_contigs with: use rule kaiju2krona_reads as kaiju2krona_contigs with: input: - db_files=get_kaiju_files(), report=rules.run_kaiju_contigs.output.report, + db_files=get_kaiju_files(), output: krona=temp( "results/{project}/output/classification/assembly/{sample}/kaiju.out.krona" diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index cd38d08..3a0a358 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -148,8 +148,10 @@ def get_card_db_file(): def get_card_annotation_file(): version = config["card"]["version"] + filename = "".join(["card_database_", version, ".fasta"]) folder = str(Path(get_card_db_file()).parent) - path = "".join([folder, "/card_database_", version, ".fasta"]) + + path = "/".join([folder, filename]) return path diff --git a/workflow/rules/resistance.smk b/workflow/rules/resistance.smk index d049578..b6a89e5 100644 --- a/workflow/rules/resistance.smk +++ b/workflow/rules/resistance.smk @@ -141,7 +141,7 @@ rule filter_uniCARD: tsv=rules.gz_uniCARD_tsv.output.tsv, json=get_CARD_hierarchy(), output: - best="results/{project}/output/resistance/uniCARD/assembly/{sample}.csv", + csv="results/{project}/output/resistance/uniCARD/assembly/{sample}.csv", log: "logs/{project}/uniCARD/{sample}_filter.log", conda: @@ -152,28 +152,55 @@ rule filter_uniCARD: filter_window=0.01, shell: "python workflow/scripts/uniCARD_filtering.py --infile {input.tsv} " - "--card_hierarchy {input.json} --outfile {output.best} " + "--card_hierarchy {input.json} --outfile {output.csv} " "--filter_window {params.filter_window} > {log} 2>&1" -""" -rule CARD_read_sample_summary: +rule resistance_abundance_per_sample: input: - txt=rules.CARD_read_run.output.txt, + csv=rules.filter_uniCARD.output.csv, + json=get_CARD_hierarchy(), output: - csv="results/{project}/output/resistance/CARD/reads/{sample}/{sample}_read_ARGs.csv", - params: - case="reads", + abd_class="results/{project}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_drug_class_abundance.csv", + abd_ARGs="results/{project}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_ARGs_abundance.csv", log: - "logs/{project}/ARGs/reads/{sample}.log", - threads: 2 + "logs/{project}/uniCARD/{sample}_abundance.log", conda: "../envs/python.yaml" + threads: 5 script: - "../scripts/arg_summary_sample.py" + "../scripts/uniCARD_abundance_sample.py" + +rule resistance_abundance_all: + input: + contig_files=expand( + "results/{{project}}/output/classification/assembly/{sample}/{sample}_classified_contigs.csv", + sample=get_samples(), + ), + abd_class_files=expand( + "results/{{project}}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_drug_class_abundance.csv", + sample=get_samples(), + ), + abd_ARGs_files=expand( + "results/{{project}}/output/resistance/uniCARD/per_sample_resistance_abundance/{sample}_ARGs_abundance.csv", + sample=get_samples(), + ), + output: + abd_class_possible="results/{project}/output/resistance/uniCARD/{project}_drug_class_abundance_all_possible.csv", + abd_class_present="results/{project}/output/resistance/uniCARD/{project}_drug_class_abundance_all_present.csv", + abd_ARGs_possible="results/{project}/output/resistance/uniCARD/{project}_ARGs_abundance_all_possible.csv", + abd_ARGs_present="results/{project}/output/resistance/uniCARD/{project}_ARGs_abundance_all_present.csv", + log: + "logs/{project}/uniCARD/all_abundances.log", + conda: + "../envs/python.yaml" + threads: 5 + script: + "../scripts/uniCARD_abundance.py" +""" # updates CARD database to use for contigs instead of reads # read based classification must be finished before rule CARD_load_DB: diff --git a/workflow/scripts/assembly_summary.py b/workflow/scripts/assembly_summary.py index f0c5efd..535e0a3 100644 --- a/workflow/scripts/assembly_summary.py +++ b/workflow/scripts/assembly_summary.py @@ -10,16 +10,18 @@ csv_mags = snakemake.input.csv_mags csv_bins = snakemake.input.csv_bins -df=pd.read_csv(qc_csv,index_col="sample") -samples=df.index.to_list() +df = pd.read_csv(qc_csv, index_col="sample") +samples = df.index.to_list() summary_dict = {} for sample in samples: summary_sample_dict = {} - summary_sample_dict["#reads_after_filtering"] = df.loc[sample,"#reads_after_filtering"] + summary_sample_dict["#reads_after_filtering"] = df.loc[ + sample, "#reads_after_filtering" + ] for asbl_log in asbl_logs: - if asbl_log.rfind(f'{sample}_megahit') >= 0: + if asbl_log.rfind(f"{sample}_megahit") >= 0: with open(asbl_log, "r") as a_log: for line in a_log: if line.find("total") >= 0: @@ -29,38 +31,44 @@ sep = sep.split() if "contigs" in sep: colname = "#contigs" - value=int(sep[0]) + value = int(sep[0]) elif "total" in sep: colname = "assembly_size(bp)" - value=int(sep[1]) + value = int(sep[1]) elif "max" in sep: colname = "longest_contig(bp)" - value=int(sep[1]) + value = int(sep[1]) elif "avg" in sep: colname = "average_contig_size(bp)" - value=int(sep[1]) + value = int(sep[1]) elif "N50" in sep: colname = "N50(bp)" - value=int(sep[1]) - + value = int(sep[1]) + summary_sample_dict[colname] = value for txt in txts: - if txt.rfind(f'{sample}_reads_mapped') >= 0: + if txt.rfind(f"{sample}_reads_mapped") >= 0: with open(txt) as t: - mapped=int(t.readline()) + mapped = int(t.readline()) summary_sample_dict["#assembled_reads"] = mapped - summary_sample_dict["%assembled_reads"] = round(((mapped/summary_sample_dict["#reads_after_filtering"]) * 100),2) + summary_sample_dict["%assembled_reads"] = round( + ((mapped / summary_sample_dict["#reads_after_filtering"]) * 100), 2 + ) break - + for bin_file in csv_bins: if bin_file.rfind(sample) >= 0: - bin_df=pd.read_csv(bin_file) + bin_df = pd.read_csv(bin_file) summary_sample_dict["#bins"] = len(bin_df) - mag_file=[file for file in csv_mags if os.path.basename(os.path.dirname(file)) == sample][0] - mag_df=pd.read_csv(mag_file) + mag_file = [ + file + for file in csv_mags + if os.path.basename(os.path.dirname(file)) == sample + ][0] + mag_df = pd.read_csv(mag_file) summary_sample_dict["#MAGs"] = len(mag_df) summary_dict[sample] = summary_sample_dict @@ -74,19 +82,22 @@ summary_df.to_csv(snakemake.output.csv) + ## edit decimal separator for visual of table def format_int_with_commas(x): - #Formats an integer with commas as thousand separators + # Formats an integer with commas as thousand separators return f"{x:,}" + def format_float_with_percent(x): - #Formats a float with percent sign + # Formats a float with percent sign return f"{x:.2f}%" -cols=[col for col in summary_df.columns.to_list() if col.find("%")>=0] -summary_df[cols]=summary_df[cols].map(format_float_with_percent) -cols=[col for col in summary_df.columns.to_list() if col.find("%")<0] -summary_df[cols]=summary_df[cols].map(format_int_with_commas) +cols = [col for col in summary_df.columns.to_list() if col.find("%") >= 0] +summary_df[cols] = summary_df[cols].map(format_float_with_percent) + +cols = [col for col in summary_df.columns.to_list() if col.find("%") < 0] +summary_df[cols] = summary_df[cols].map(format_int_with_commas) summary_df.to_csv(snakemake.output.vis_csv) diff --git a/workflow/scripts/uniCARD_abundance.py b/workflow/scripts/uniCARD_abundance.py new file mode 100644 index 0000000..3235bf4 --- /dev/null +++ b/workflow/scripts/uniCARD_abundance.py @@ -0,0 +1,52 @@ +import os +import pandas as pd +import sys + +sys.stderr = open(snakemake.log[0], "w") + + +contig_files = snakemake.input.contig_files +abd_class_files = snakemake.input.abd_class_files +abd_ARGs_files = snakemake.input.abd_ARGs_files + + +## Combining contig and ARG information for all samples +classes_all_samples = [] +names_all_samples = [] +for contig_file in contig_files: + sample_id = os.path.basename(os.path.dirname(contig_file)) + + classification_df = pd.read_csv(contig_file) + abd_class_file = [file for file in abd_class_files if file.find(sample_id) >= 0][0] + abd_classes_df = pd.read_csv(abd_class_file) + abd_ARGs_file = [file for file in abd_ARGs_files if file.find(sample_id) >= 0][0] + abd_names_df = pd.read_csv(abd_ARGs_file) + + merged_classes_df = classification_df.merge( + abd_classes_df, on="contig", how="inner" + ) + merged_classes_df.insert(0, "sampleID", sample_id) + + classes_all_samples.append(merged_classes_df) + + merged_names_df = classification_df.merge(abd_names_df, on="contig", how="inner") + merged_names_df.insert(0, "sampleID", sample_id) + + names_all_samples.append(merged_names_df) + + +final_classes_df = pd.concat(classes_all_samples, ignore_index=True) +final_classes_df.to_csv(snakemake.output.abd_class_possible, index=False) + +# in antibiotic classes per sample file all antibiotic classes in CARD are included +# All that do not appear in at least one sample are removed: +final_classes_df_red = final_classes_df.loc[:, (final_classes_df != 0).any(axis=0)] +final_classes_df_red.to_csv(snakemake.output.abd_class_present, index=False) + +final_names_df = pd.concat(names_all_samples, ignore_index=True) +final_names_df.to_csv(snakemake.output.abd_ARGs_possible, index=False) + +# in ARG per sample file all ARGs in CARD are included +# All that do not appear in at least one sample are removed: +final_names_df_red = final_names_df.loc[:, (final_names_df != 0).any(axis=0)] +final_names_df_red.to_csv(snakemake.output.abd_ARGs_present, index=False) diff --git a/workflow/scripts/uniCARD_abundance_sample.py b/workflow/scripts/uniCARD_abundance_sample.py new file mode 100644 index 0000000..1f5553a --- /dev/null +++ b/workflow/scripts/uniCARD_abundance_sample.py @@ -0,0 +1,95 @@ +import pandas as pd +import json +import sys + +sys.stderr = open(snakemake.log[0], "w") + +card_drugs = snakemake.input.json +args_file = snakemake.input.csv +# abundance outfiles +abd_class = snakemake.output.abd_class +abd_ARGs = snakemake.output.abd_ARGs + +## Write drug classes and ARO names abundance: +with open(card_drugs) as f: + data = json.load(f) + +# Convert to DataFrame +df_card = pd.DataFrame.from_dict(data, orient="index") + +# Split semicolon-separated values into lists +df_card["classes"] = df_card["classes"].str.split("; ") +df_card["antibiotics"] = df_card["antibiotics"].str.split("; ") + +# Flatten and get unique values +all_classes = pd.Series( + [item for sublist in df_card["classes"].dropna() for item in sublist] +).unique() +all_classes.sort() +all_antibiotics = pd.Series( + [item for sublist in df_card["antibiotics"].dropna() for item in sublist if item] +).unique() +all_aros = df_card["ARO_name"].unique() +all_aros.sort() + + +##Create classes and ARG names abundance files +df_args = pd.read_csv(args_file) +df_args_for_classes = df_args.copy() +# Step 1: Preprocess the classes column into lists +df_args_for_classes["class_list"] = ( + df_args_for_classes["classes"] + .fillna("") + .apply(lambda x: [c.strip() for c in x.split(";")]) +) + +# Step 2: Create one-hot columns per class +for c in all_classes: + df_args_for_classes[c] = df_args_for_classes["class_list"].apply( + lambda lst: int(c in lst) + ) + +# Step 3: Group by contig +agg_dict = {"ARO_ID": "count"} +agg_dict.update( + {c: "sum" for c in all_classes} # was max before +) # "sum" counts how many times the class appears + + +class_abundance_df = df_args_for_classes.groupby("contig").agg(agg_dict).reset_index() +# Write out abundance file +class_abundance_df = class_abundance_df.rename(columns={"ARO_ID": "#ARGs"}) +class_abundance_df.to_csv(abd_class, index=False) + +df_args["present"] = 1 + +# Step 2: Pivot to get binary presence/absence for all ARO_name per contig +aro_matrix = df_args.pivot_table( + index="contig", + columns="ARO_name", + values="present", + aggfunc="sum", # <-- changed from "max" to "sum" + fill_value=0, +) + +# Step 3: Count ARGs (genes) per contig +gene_counts = df_args.groupby("contig").size().rename("#ARGs") + +# Step 4: Join gene counts and presence matrix +names_abundance_df = pd.concat([gene_counts, aro_matrix], axis=1).reset_index() + +# Find missing AROs +missing_aros = [aro for aro in all_aros if aro not in names_abundance_df.columns] + +# Create a DataFrame with 0s for missing AROs (same index as result) +missing_df = pd.DataFrame(0, index=names_abundance_df.index, columns=missing_aros) + +# Concatenate in one go — avoids fragmentation +names_abundance_df = pd.concat([names_abundance_df, missing_df], axis=1) + +# Optional: reorder columns — num_genes first, then sorted AROs +aro_columns = sorted([a for a in all_aros if a in names_abundance_df.columns]) +names_abundance_df = names_abundance_df[["contig", "#ARGs"] + aro_columns] + +# Write out abundance file +names_abundance_df.to_csv(abd_ARGs, index=False) From 79e315d6443bb663ba4d5147f2cbb68cee278c33 Mon Sep 17 00:00:00 2001 From: Josefa Welling <82578997+josefawelling@users.noreply.github.com> Date: Thu, 16 Apr 2026 09:11:50 +0200 Subject: [PATCH 11/15] format classify.smk --- workflow/rules/classify.smk | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 5a40f69..b021589 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -16,13 +16,16 @@ if not config["kaiju"]["use-local"]: "wget -c {params.download} -O - | " "tar -zxv -C {params.db_folder}) > {log} 2>&1" + # classification of reads rule run_kaiju_reads: input: db_files=get_kaiju_files(), fastqs=get_filtered_gz_fastqs, output: - report=temp("results/{project}/output/classification/reads/{sample}/kaiju.out"), + report=temp( + "results/{project}/output/classification/reads/{sample}/kaiju.out" + ), params: evalue=0.00001, threads: 64 @@ -85,7 +88,9 @@ rule run_kaiju_contigs: db_files=get_kaiju_files(), assembly=get_gz_assembly, output: - report=temp("results/{project}/output/classification/assembly/{sample}/kaiju.out"), + report=temp( + "results/{project}/output/classification/assembly/{sample}/kaiju.out" + ), threads: 60 log: "logs/{project}/contig_classification/{sample}_kaiju_run.log", @@ -130,7 +135,8 @@ use rule kaiju2krona_reads as kaiju2krona_contigs with: krona=temp( "results/{project}/output/classification/assembly/{sample}/kaiju.out.krona" ), - html=report("results/{project}/output/report/{sample}/{sample}_contigs_kaiju.out.html", + html=report( + "results/{project}/output/report/{sample}/{sample}_contigs_kaiju.out.html", htmlindex="index.html", category="5. Taxonomic classification", subcategory="5.3 Contig classification", From 5e9b1ee8990c156296e773f1e6836beb486091b9 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Thu, 16 Apr 2026 07:58:55 +0000 Subject: [PATCH 12/15] formatting --- .github/workflows/main.yml | 51 ++++++++++++++++++------------------- workflow/rules/classify.smk | 8 ++++-- 2 files changed, 31 insertions(+), 28 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 0dee238..3be8ff8 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -2,16 +2,15 @@ name: Tests on: push: - branches: [ main ] + branches: [main] pull_request: - branches: [ main ] - + branches: [main] jobs: Formatting: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v6 - name: Formatting uses: github/super-linter@v5 env: @@ -23,33 +22,33 @@ jobs: Linting: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 - - name: Lint workflow - uses: snakemake/snakemake-github-action@v1 - with: - directory: . - snakefile: workflow/Snakefile - args: "--lint" + - uses: actions/checkout@v6 + - name: Lint workflow + uses: snakemake/snakemake-github-action@v2 + with: + directory: . + snakefile: workflow/Snakefile + args: "--lint" # Testing: - # runs-on: ubuntu-latest - # needs: - # - Linting - # - Formatting +# runs-on: ubuntu-latest +# needs: +# - Linting +# - Formatting # steps: - # - uses: actions/checkout@v3 +# - uses: actions/checkout@v3 # - name: Test workflow - # uses: snakemake/snakemake-github-action@v1 - # with: - # directory: .tests - # snakefile: workflow/Snakefile - # args: "--use-conda --show-failed-logs --cores 3 --conda-cleanup-pkgs cache --all-temp" +# uses: snakemake/snakemake-github-action@v1 +# with: +# directory: .tests +# snakefile: workflow/Snakefile +# args: "--use-conda --show-failed-logs --cores 3 --conda-cleanup-pkgs cache --all-temp" # - name: Test report - # uses: snakemake/snakemake-github-action@v1 - # with: - # directory: .test - # snakefile: workflow/Snakefile - # args: "--report report.zip" +# uses: snakemake/snakemake-github-action@v1 +# with: +# directory: .test +# snakefile: workflow/Snakefile +# args: "--report report.zip" diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index 5a40f69..7e752ad 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -16,6 +16,7 @@ if not config["kaiju"]["use-local"]: "wget -c {params.download} -O - | " "tar -zxv -C {params.db_folder}) > {log} 2>&1" + # classification of reads rule run_kaiju_reads: input: @@ -85,7 +86,9 @@ rule run_kaiju_contigs: db_files=get_kaiju_files(), assembly=get_gz_assembly, output: - report=temp("results/{project}/output/classification/assembly/{sample}/kaiju.out"), + report=temp( + "results/{project}/output/classification/assembly/{sample}/kaiju.out" + ), threads: 60 log: "logs/{project}/contig_classification/{sample}_kaiju_run.log", @@ -130,7 +133,8 @@ use rule kaiju2krona_reads as kaiju2krona_contigs with: krona=temp( "results/{project}/output/classification/assembly/{sample}/kaiju.out.krona" ), - html=report("results/{project}/output/report/{sample}/{sample}_contigs_kaiju.out.html", + html=report( + "results/{project}/output/report/{sample}/{sample}_contigs_kaiju.out.html", htmlindex="index.html", category="5. Taxonomic classification", subcategory="5.3 Contig classification", From 28b46f9ecd2eb76c48552f4f6c526f25a13dd017 Mon Sep 17 00:00:00 2001 From: Josefa Welling <82578997+josefawelling@users.noreply.github.com> Date: Thu, 16 Apr 2026 10:06:53 +0200 Subject: [PATCH 13/15] format classify.smk --- workflow/rules/classify.smk | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index b021589..7e752ad 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -23,9 +23,7 @@ rule run_kaiju_reads: db_files=get_kaiju_files(), fastqs=get_filtered_gz_fastqs, output: - report=temp( - "results/{project}/output/classification/reads/{sample}/kaiju.out" - ), + report=temp("results/{project}/output/classification/reads/{sample}/kaiju.out"), params: evalue=0.00001, threads: 64 From 41e4d163a788ea779eb7937542d6a4426f81070f Mon Sep 17 00:00:00 2001 From: josefawelling Date: Thu, 16 Apr 2026 08:53:09 +0000 Subject: [PATCH 14/15] format --- config/config.yaml | 24 +++++++++--------- workflow/rules/analysis.smk | 2 +- workflow/rules/bin_qc.smk | 32 ++++++++++++------------ workflow/rules/classify.smk | 6 ++--- workflow/rules/common.smk | 28 ++++----------------- workflow/rules/host_filtering.smk | 22 ++++++++--------- workflow/rules/resistance.smk | 41 ++++++++++++++----------------- 7 files changed, 66 insertions(+), 89 deletions(-) diff --git a/config/config.yaml b/config/config.yaml index bf7dfc1..61b3a1a 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -15,11 +15,11 @@ adapter-seqs: "-a CTGTCTCTTATACACATCT -g AGATGTGTATAAGAGACAG" # reads are always filtered for human contamination human-filtering: # if you have a human reference genome on your device, - # you can set use-local to True and specify the path to the human reference genome - use-local: True + # you can set use-shared to True and specify the path to the human reference genome + use-shared: True # path to a locally stored human reference genome local-path: /groups/ds/databases_refGenomes/refGenomes/latest/human/GRCh38_latest_genomic.fna.gz - # if use-local = False the reference genome is downloaded via the following url + # if use-shared = False the reference genome is downloaded via the following url download-path: https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.fna.gz # if filtering for non-human host is required, @@ -44,33 +44,33 @@ MAG-criteria: gtdb: /groups/ds/databases_refGenomes/databases/latest/gtdb/release226/ kaiju: - use-local: True + use-shared: True # path of kaiju DB .fmi file - # if use-local False: DB is downloaded and saved to this path + # if use-shared False: DB is downloaded and saved to this path fmi-file: /groups/ds/databases_refGenomes/databases/latest/kaiju/refseq_nr/unpacked_2024-08-13/kaiju_db_refseq_nr.fmi download: https://kaiju-idx.s3.eu-central-1.amazonaws.com/2024/kaiju_db_refseq_nr_2024-08-13.tgz checkm2: - use-local: True - #if use-local False a folder named CheckM2_database will be created under the given path + use-shared: True + #if use-shared False a folder named CheckM2_database will be created under the given path db-folder: /groups/ds/databases_refGenomes/databases/latest/checkM2/ # folder of genomad DB genomad: - use-local: True + use-shared: True # a folder named genomad_db is expected under the given path - # if use-local False this folder will be created + # if use-shared False this folder will be created db-folder: /groups/ds/databases_refGenomes/databases/latest/ card: - use-local: True + use-shared: True db-file: /groups/ds/databases_refGenomes/databases/latest/card/card.json version: v4.0.1 url: https://card.mcmaster.ca/latest/data uniCARD: - #path to the .dmnd file of UniCARD database - db-file: /groups/ds/databases_refGenomes/databases/latest/uniCARD/uniCARD.dmnd + #path to the fasta.gz file of UniCARD database + db-file: /groups/ds/databases_refGenomes/databases/latest/uniCARD/uniCARD.fasta.gz # path to CARDs hierarchy file hierarchy-json: /groups/ds/databases_refGenomes/databases/latest/uniCARD/CARD_hierarchy_v4.0.0.json diff --git a/workflow/rules/analysis.smk b/workflow/rules/analysis.smk index 196df26..6a5262e 100644 --- a/workflow/rules/analysis.smk +++ b/workflow/rules/analysis.smk @@ -37,7 +37,7 @@ rule gzip_proteins: # Plasmid analysis -if not config["genomad"]["use-local"]: +if not config["genomad"]["use-shared"]: rule load_genomad_DB: output: diff --git a/workflow/rules/bin_qc.smk b/workflow/rules/bin_qc.smk index c8f1b79..519b926 100644 --- a/workflow/rules/bin_qc.smk +++ b/workflow/rules/bin_qc.smk @@ -1,15 +1,15 @@ ## bin QC -if not config["checkm2"]["use-local"]: +if not config["checkm2"]["use-shared"]: rule checkm2_DB_download: output: dbfile=get_checkm2_db(), - params: - direct=lambda wildcards, output: Path(output.dbfile).parent.parent, log: "logs/checkm2_DB_download.log", conda: "../envs/checkm2.yaml" + params: + direct=lambda wildcards, output: Path(output.dbfile).parent.parent, shell: "checkm2 database --download --path {params.direct} > {log} 2>&1" @@ -20,13 +20,13 @@ rule checkm2_run: dbfile=get_checkm2_db(), output: stats="results/{project}/output/report/{sample}/checkm2_quality_report.tsv", - params: - outdir="results/{project}/qc/checkm2/{sample}/", log: "logs/{project}/checkm2/{sample}.log", - threads: 24 conda: "../envs/checkm2.yaml" + threads: 24 + params: + outdir="results/{project}/qc/checkm2/{sample}/", shell: "(checkm2 predict -x fa.gz --threads {threads} --force " "--input {input.bins}/ --output-directory {params.outdir}/ " @@ -43,14 +43,14 @@ rule bin_summary_sample: csv_bins="results/{project}/output/report/{sample}/{sample}_bin_summary.csv", csv_tax="results/{project}/output/report/{sample}/{sample}_bin_taxonomy.csv", csv_mags="results/{project}/output/report/{sample}/{sample}_mags_summary.csv", - params: - max_cont=config["MAG-criteria"]["max-contamination"], - min_comp=config["MAG-criteria"]["min-completeness"], log: "logs/{project}/bin_summary/{sample}.log", - threads: 4 conda: "../envs/python.yaml" + threads: 4 + params: + max_cont=config["MAG-criteria"]["max-contamination"], + min_comp=config["MAG-criteria"]["min-completeness"], script: "../scripts/bin_summary_sample.py" @@ -68,14 +68,14 @@ use rule qc_summary_report as bin_sample_report with: labels={"sample": "{sample}"}, ) ), + log: + "logs/{project}/report/{sample}/bin_rbt_csv.log", params: pin_until="bin", styles="resources/report/tables/", name="{sample}_bin_summary", header="Bin summary for sample {sample}", pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/bin_rbt_csv.log", use rule qc_summary_report as taxonomy_report with: @@ -91,14 +91,14 @@ use rule qc_summary_report as taxonomy_report with: labels={"sample": "{sample}"}, ) ), + log: + "logs/{project}/report/{sample}/taxonomy_rbt_csv.log", params: pin_until="bin", styles="resources/report/tables/", name="{sample}_taxonomy_summary", header="Taxonomy summary for sample {sample}", pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/taxonomy_rbt_csv.log", use rule qc_summary_report as mag_report with: @@ -114,11 +114,11 @@ use rule qc_summary_report as mag_report with: labels={"sample": "{sample}"}, ) ), + log: + "logs/{project}/report/{sample}/mag_rbt_csv.log", params: pin_until="MAG", styles="resources/report/tables/", name="{sample}_MAG_summary", header="MAG summary for sample {sample}", pattern=config["tablular-config"], - log: - "logs/{project}/report/{sample}/mag_rbt_csv.log", diff --git a/workflow/rules/classify.smk b/workflow/rules/classify.smk index b021589..c785336 100644 --- a/workflow/rules/classify.smk +++ b/workflow/rules/classify.smk @@ -1,5 +1,5 @@ # if there is no local database version to use, it is downloaded -if not config["kaiju"]["use-local"]: +if not config["kaiju"]["use-shared"]: rule download_kaiju: output: @@ -23,9 +23,7 @@ rule run_kaiju_reads: db_files=get_kaiju_files(), fastqs=get_filtered_gz_fastqs, output: - report=temp( - "results/{project}/output/classification/reads/{sample}/kaiju.out" - ), + report=temp("results/{project}/output/classification/reads/{sample}/kaiju.out"), params: evalue=0.00001, threads: 64 diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 3a0a358..7169886 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -98,7 +98,7 @@ def get_human_ref_download(): def get_human_ref(): - if config["human-filtering"]["use-local"]: + if config["human-filtering"]["use-shared"]: path = config["human-filtering"]["local-path"] return path else: @@ -164,33 +164,15 @@ def get_CARD_hierarchy(): def get_uniCARD_db_wo_ext(): - db_wo_ext = str(Path(get_uniCARD_db()).with_suffix("")) - return db_wo_ext - - -""" -def is_uniCARD_fasta(): - filename = str(Path(get_uniCARD_db()).name) - ext = "".join(Path(filename).suffixes) - if ext in ["fasta", "fasta.gz"]: - return True - else: - return False - - -def get_uniCARD_db_wo_ext(): - db_wo_ext = str(Path(get_uniCARD_db()).with_suffix("")) - if is_uniCARD_fasta(): - db_wo_ext = db_wo_ext.with_suffix("") + db_wo_ext = str(Path(get_uniCARD_db()).with_suffix("").with_suffix("")) return db_wo_ext def get_unicard_dmnd(): - if is_uniCARD_fasta(): - return "".join([get_uniCARD_db_wo_ext(), ".dmnd"]) - else: - return get_uniCARD_db() + return "".join([get_uniCARD_db_wo_ext(), ".dmnd"]) + +""" def get_mag_fa(wildcards): folder = "results/{}/output/fastas/{}/mags/".format( diff --git a/workflow/rules/host_filtering.smk b/workflow/rules/host_filtering.smk index 258e1fd..e896199 100644 --- a/workflow/rules/host_filtering.smk +++ b/workflow/rules/host_filtering.smk @@ -1,20 +1,20 @@ from pathlib import Path # Download human reference genome if not local file is given -if not config["human-filtering"]["use-local"]: +if not config["human-filtering"]["use-shared"]: rule download_human_ref: output: fasta=get_human_ref(), - params: - download=get_human_ref_download(), - folder=lambda wildcards, output: Path(output.fasta).parent, log: "logs/human_ref_download.log", group: "refGenome_depended" conda: "../envs/unix.yaml" + params: + download=get_human_ref_download(), + folder=lambda wildcards, output: Path(output.fasta).parent, shell: "(mkdir -p {params.folder} && " "cd {params.folder} && " @@ -27,11 +27,11 @@ rule map_to_human: ref=get_human_ref(), output: bam=temp("results/{project}/human_filtering/alignments/{sample}.bam"), - threads: 64 log: "logs/{project}/human_filtering/map_to_human_{sample}.log", conda: "../envs/minimap2.yaml" + threads: 64 shell: "(minimap2 -a -xsr -t {threads} {input.ref} {input.fastqs} | " "samtools view -bh | " @@ -43,11 +43,11 @@ rule index_human_alignment: rules.map_to_human.output.bam, output: bai=temp("results/{project}/human_filtering/alignments/{sample}.bam.bai"), - threads: 20 log: "logs/{project}/human_filtering/index_human_alignment_{sample}.log", conda: "../envs/minimap2.yaml" + threads: 20 shell: "samtools index {input} > {log} 2>&1" @@ -63,11 +63,11 @@ rule filter_human: read=["R1", "R2"], ) ), - threads: 64 log: "results/{project}/output/report/prerequisites/qc/{sample}_filter_human.log", conda: "../envs/minimap2.yaml" + threads: 64 shell: "(samtools fastq --threads {threads} -F 3584 -f 77 " "-o {output.filtered[0]} {input.bam} && " @@ -82,10 +82,10 @@ rule gzip_filtered_reads: "results/{project}/output/filtered_reads/{sample}_{read}.fastq.gz", log: "logs/{project}/human_filtering/gzip_{sample}_{read}.log", - threads: 20 priority: 1 conda: "../envs/unix.yaml" + threads: 20 shell: "gzip -k {input} > {log} 2>&1" @@ -100,18 +100,18 @@ if config["host-filtering"]["do-host-filtering"]: ref=get_host_ref(), output: bam=temp("results/{project}/host_filtering/alignments/{sample}.bam"), - threads: 20 log: "logs/{project}/host_filtering/map_to_host_{sample}.log", + threads: 20 use rule index_human_alignment as index_host_alignment with: input: rules.map_to_host.output.bam, output: bai=temp("results/{project}/host_filtering/alignments/{sample}.bam.bai"), - threads: 3 log: "logs/{project}/host_filtering/index_host_alignment_{sample}.log", + threads: 3 rule filter_host: input: @@ -124,11 +124,11 @@ if config["host-filtering"]["do-host-filtering"]: read=["R1", "R2"], ) ), - threads: 64 log: "results/{project}/output/report/prerequisites/qc/{sample}_filter_host.log", conda: "../envs/minimap2.yaml" + threads: 64 shell: "(samtools fastq -F 3584 -f 77 {input.bam} | " "gzip -c > {output.filtered[0]} && " diff --git a/workflow/rules/resistance.smk b/workflow/rules/resistance.smk index b6a89e5..95a7d14 100644 --- a/workflow/rules/resistance.smk +++ b/workflow/rules/resistance.smk @@ -1,5 +1,5 @@ # Resistance analysis -if not config["card"]["use-local"]: +if not config["card"]["use-shared"]: rule download_CARD_data: output: @@ -71,32 +71,28 @@ rule CARD_read_run: "--clean -n {threads} > {log} 2>&1" -""" -if is_uniCARD_fasta: - - rule uniCARD_makeDB: - input: - db=get_uniCARD_db(), - output: - dmdb=get_unicard_dmnd(), - log: - "logs/unicard_makeDB.log", - conda: - "../envs/diamond.yaml" - threads: 60 - params: - filename=lambda wildcards, input: Path(input.db).name, - filename_wo_ext=Path(get_uniCARD_db_wo_ext()).name, - folder=lambda wildcards, input: Path(input.db).parent, - shell: - "(cd {params.folder} && diamond makedb --in {params.filename} -d {params.filename_wo_ext}) > {log} 2>&1" -""" +rule uniCARD_makeDB: + input: + db=get_uniCARD_db(), + output: + dmdb=get_unicard_dmnd(), + log: + "logs/unicard_makeDB.log", + conda: + "../envs/diamond.yaml" + threads: 60 + params: + filename=lambda wildcards, input: Path(input.db).name, + filename_wo_ext=Path(get_uniCARD_db_wo_ext()).name, + folder=lambda wildcards, input: Path(input.db).parent, + shell: + "(cd {params.folder} && diamond makedb --in {params.filename} -d {params.filename_wo_ext}) > {log} 2>&1" rule uniCARD_run: input: faa=get_proteins, - dmdb=get_uniCARD_db(), + dmdb=get_unicard_dmnd(), output: tsv=temp( "results/{project}/output/resistance/uniCARD/assembly/{sample}_out.tsv" @@ -113,6 +109,7 @@ rule uniCARD_run: heavy=1, params: db_wo_ext=get_uniCARD_db_wo_ext(), + #db_wo_ext=lambda w, input: os.path.splitext(input.dmdb)[0], max_targets=15, shell: "diamond blastp -q {input.faa} -d {params.db_wo_ext} " From 447717a64078d209dc5572dc316d96c86b878cc1 Mon Sep 17 00:00:00 2001 From: josefawelling Date: Thu, 16 Apr 2026 12:32:38 +0000 Subject: [PATCH 15/15] update readme --- README.md | 7 ++++++- workflow/rules/resistance.smk | 3 +-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index cd87352..f0c7b89 100644 --- a/README.md +++ b/README.md @@ -50,6 +50,10 @@ The GTDB needs to be downloaded and decompressed, it requires about 140 Gb. 4. After successful step 3: the archive can be removed 5. Please specify the path to your decompressed GTDB in the config file (see [Configuring workflow](#configuring-the-workflow)) +### Create UniCARD database +Please follow the steps outlined in the [UniCARD repository](https://github.com/IKIM-Essen/uniCARD) to create a UniCARD database for resistance gene annotation. +And specify the path to your database in the config file (see [Configuring workflow](#configuring-the-workflow)) + --- ## Usage @@ -65,7 +69,8 @@ git clone https://github.com/IKIM-Essen/metagenomics_res.git - Specify filtering options for human reads (`human-filtering`) - Specify host filtering options, if you have a non-human host (`host-filtering`) - Specify options for different databases: - - GTDB database needs to be downloaded before (see [Download GTDB](#Download-GTDB)) + - GTDB database needs to be downloaded before (see [Download GTDB](#download-gtdb)) + - UniCARD database needs to be create before (see [Create UniCARD database](#create-unicard-database)) - other databases (kaiju, CheckM2, CARD, genomad) can be given as a local path or downloaded when running the pipeline 2. Provide sample information in the `config/pep/samples.csv` file while keeping the header and the format as: diff --git a/workflow/rules/resistance.smk b/workflow/rules/resistance.smk index 95a7d14..08dfcef 100644 --- a/workflow/rules/resistance.smk +++ b/workflow/rules/resistance.smk @@ -108,8 +108,7 @@ rule uniCARD_run: resources: heavy=1, params: - db_wo_ext=get_uniCARD_db_wo_ext(), - #db_wo_ext=lambda w, input: os.path.splitext(input.dmdb)[0], + db_wo_ext=lambda w, input: os.path.splitext(input.dmdb)[0], max_targets=15, shell: "diamond blastp -q {input.faa} -d {params.db_wo_ext} "