From e545bc250046a9785c4a337212eca251a58979b8 Mon Sep 17 00:00:00 2001 From: Rina Ahmed-Begrich Date: Wed, 9 Sep 2026 15:15:37 +0200 Subject: [PATCH 1/4] update ntsynt tool versions --- workflow/envs/ntsynt.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/workflow/envs/ntsynt.yml b/workflow/envs/ntsynt.yml index 7c3e9db..410766f 100644 --- a/workflow/envs/ntsynt.yml +++ b/workflow/envs/ntsynt.yml @@ -4,5 +4,5 @@ channels: - bioconda - nodefaults dependencies: - - ntsynt=1.0.5 - - ntsynt-viz=1.0.4 + - ntsynt=1.0.9 + - ntsynt-viz=1.1.1 From 9fae1ec17c4830dec31314a17601e3724a93bd8a Mon Sep 17 00:00:00 2001 From: Rina Ahmed-Begrich Date: Wed, 9 Sep 2026 15:18:23 +0200 Subject: [PATCH 2/4] feat: add pairwise whole-genome-comparison and variant calling --- .test/config/config.yml | 11 +++ config/config.yml | 11 +++ config/schemas/config.schema.yml | 50 +++++++++++++ workflow/envs/vcfutils.yml | 13 ++++ workflow/rules/common.smk | 11 +++ workflow/rules/qc.smk | 114 ++++++++++++++++++++++++++++++ workflow/scripts/create_dotplot.R | 46 ++++++++++++ 7 files changed, 256 insertions(+) create mode 100644 workflow/envs/vcfutils.yml create mode 100644 workflow/scripts/create_dotplot.R diff --git a/.test/config/config.yml b/.test/config/config.yml index 88685e9..b3fdaf7 100644 --- a/.test/config/config.yml +++ b/.test/config/config.yml @@ -50,3 +50,14 @@ synteny: extra: "" viz_scale: "1e6" viz_extra: "--normalize" + +reference_comparison: + skip: False + minimap2: + extra: "-c --cs" + sorting: "coordinate" + sort_extra: "--no-PG" + paftools: + extra: "" + bcftools: + extra: "" diff --git a/config/config.yml b/config/config.yml index eb4a738..42c448c 100644 --- a/config/config.yml +++ b/config/config.yml @@ -50,3 +50,14 @@ synteny: extra: "" viz_scale: "1e6" viz_extra: "--normalize" + +reference_comparison: + skip: False + minimap2: + extra: "-c --cs" + sorting: "coordinate" + sort_extra: "--no-PG" + paftools: + extra: "" + bcftools: + extra: "" \ No newline at end of file diff --git a/config/schemas/config.schema.yml b/config/schemas/config.schema.yml index 9014c40..d62199c 100644 --- a/config/schemas/config.schema.yml +++ b/config/schemas/config.schema.yml @@ -200,6 +200,55 @@ properties: - extra - viz_extra - viz_scale + reference_comparison: + type: object + properties: + skip: + type: boolean + description: Whether to skip whole-genome reference comparison + default: false + minimap2: + type: object + properties: + extra: + type: string + description: Extra command-line arguments for minimap2 analysis + default: "-c --cs" + sorting: + type: string + description: Sorting method for minimap2 output (e.g., "coordinate" or "queryname") + default: "coordinate" + sort_extra: + type: string + description: Extra command-line arguments for sorting minimap2 output + default: "--no-PG" + required: + - extra + - sorting + - sort_extra + paftools: + type: object + properties: + extra: + type: string + description: Extra command-line arguments for paftools analysis + default: "" + required: + - extra + bcftools: + type: object + properties: + extra: + type: string + description: Extra command-line arguments for bcftools analysis + default: "" + required: + - extra + required: + - skip + - minimap2 + - paftools + - bcftools required: - samplesheet - tool @@ -213,3 +262,4 @@ required: - checkm - rgi - synteny + - reference_comparison diff --git a/workflow/envs/vcfutils.yml b/workflow/envs/vcfutils.yml new file mode 100644 index 0000000..ce51ca1 --- /dev/null +++ b/workflow/envs/vcfutils.yml @@ -0,0 +1,13 @@ +name: vcfutils +channels: + - conda-forge + - bioconda + - nodefaults +dependencies: + - minimap2=2.31 + - bcftools=1.24 + - perl-vcftools-vcf=0.1.15 + - snippy=4.0.2 + - htslib=1.24 + - r-pafr=0.0.2 + - r-readr=2.2.0 diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 1ff0a7e..006b5c4 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -96,6 +96,17 @@ def get_final_input(wildcards): inputs += expand( "results/qc/genome_synteny/ntSynt-viz_ribbon-plot.pdf", ) + if ( + config["reference"]["fasta"] != "" + and not config["reference_comparison"]["skip"] + ): + inputs += expand( + "results/qc/reference_comparison/all_merged_aln.vcf.gz", + ) + inputs += expand( + "results/qc/reference_comparison/{sample}_dotplot.pdf", + sample=samples.index, + ) return inputs diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index 3d2128a..57b086c 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -305,3 +305,117 @@ rule viz_synteny: echo "Clean intermediate files." >>{log} rm -f ./ntSynt-viz.*.tsv ./ntSynt-viz_* ./ntSynt.*.tsv """ + + +rule minimap2_paf: + input: + target=config["reference"]["fasta"], + query=get_fasta, + output: + "results/qc/reference_comparison/{sample}_aln.paf", + log: + "results/qc/reference_comparison/logs/{sample}_aln.log", + threads: max(workflow.cores * 0.25, 1) + params: + extra=config["reference_comparison"]["minimap2"]["extra"], + sorting=config["reference_comparison"]["minimap2"]["sorting"], + sort_extra=config["reference_comparison"]["minimap2"]["sort_extra"], + message: + """--- Running assembly-to-assembly mapping using minimap2 ---""" + wrapper: + "v9.7.0/bio/minimap2/aligner" + + +rule paftools_vcf: + input: + target=config["reference"]["fasta"], + query=rules.minimap2_paf.output, + output: + tab="results/qc/reference_comparison/{sample}_aln.tab", + vcf="results/qc/reference_comparison/{sample}_aln.vcf", + vcf_gz="results/qc/reference_comparison/{sample}_aln.vcf.gz", + index="results/qc/reference_comparison/{sample}_aln.vcf.gz.tbi", + log: + "results/qc/reference_comparison/logs/{sample}_vcf.log", + conda: + "../envs/vcfutils.yml" + threads: max(workflow.cores * 0.25, 1) + params: + extra=config["reference_comparison"]["paftools"]["extra"], + message: + """--- Running assembly-to-assembly mapping using minimap2 ---""" + shell: + """ + sort -k6,6 -k8,8n {input.query} \ + | paftools.js call -f {input.target} - \ + {params.extra} 2>{log} \ + | bcftools reheader \ + --threads {threads} \ + -n {wildcards.sample} \ + | vcf-annotate --fill-type \ + >{output.vcf} 2>>{log} + snippy-vcf_to_tab --ref {input.target} \ + --vcf {output.vcf} \ + >{output.tab} 2>>{log} + bgzip -c {output.vcf} >{output.vcf_gz} 2>>{log} + tabix {output.vcf_gz} + """ + + +rule merge_vcfs: + input: + vcfs=expand( + "results/qc/reference_comparison/{sample}_aln.vcf.gz", sample=samples.index + ), + output: + merged_vcf_gz="results/qc/reference_comparison/all_merged_aln.vcf.gz", + index="results/qc/reference_comparison/all_merged_aln.vcf.gz.tbi", + log: + "results/qc/reference_comparison/logs/merge_vcfs.log", + conda: + "../envs/vcfutils.yml" + threads: max(workflow.cores * 0.5, 1) + params: + extra=config["reference_comparison"]["bcftools"]["extra"], + num_input=len(samples.index), + message: + """--- Merging individual VCF files into a single multi-sample VCF ---""" + shell: + """ + if [ {params.num_input} -eq 1 ]; then + echo "Only one VCF file present. Copying to merged output." >{log} + cp {input.vcfs} {output.merged_vcf_gz} + else + echo "Merging {params.num_input} VCF files using bcftools..." >{log} + bcftools merge \ + {input.vcfs} \ + {params.extra} \ + --output {output.merged_vcf_gz} \ + >>{log} 2>&1 + fi + tabix {output.merged_vcf_gz} + """ + + +rule dotplot: + input: + query=rules.minimap2_paf.output, + output: + pdf="results/qc/reference_comparison/{sample}_dotplot.pdf", + png="results/qc/reference_comparison/{sample}_dotplot.png", + log: + path="results/qc/reference_comparison/logs/{sample}_dotplot.log", + conda: + "../envs/vcfutils.yml" + threads: 1 + params: + query=lambda wc: samples.loc[samples["sample"] == wc.sample, "strain"].values[0], + target=( + config["reference"]["name"] + if config["reference"]["name"] + else os.basename(config["reference"]["fasta"]) + ), + message: + """--- Generating dotplot for assembly-to-assembly comparison ---""" + script: + "../scripts/create_dotplot.R" diff --git a/workflow/scripts/create_dotplot.R b/workflow/scripts/create_dotplot.R new file mode 100644 index 0000000..a2e60ec --- /dev/null +++ b/workflow/scripts/create_dotplot.R @@ -0,0 +1,46 @@ +library(pafr, quietly=TRUE) +library(readr, quietly=TRUE) + +paf_file <- snakemake@input[["query"]] +output_file <- snakemake@output[["pdf"]] +output_png <- snakemake@output[["png"]] +target_name <- snakemake@params[["target"]] +query_name <- snakemake@params[["query"]] + +messages <- c() + +alignments <- read_paf(paf_file) + +# Create plot without printing to avoid Rplots.pdf +p <- dotplot(alignments, xlab = query_name, ylab = target_name) + + theme_bw() + + theme( + axis.text.x = element_text(size = 12, face="bold"), + axis.text.y = element_text(size = 12, face="bold"), + axis.title.x = element_text(size = 14, face = "bold"), + axis.title.y = element_text(size = 14, face = "bold") + ) + +ggsave( + p, + filename = output_file, + width = 8, + height = 8 +) + +ggsave( + p, + filename = output_png, + width = 8, + height = 8 +) + +messages <- append( + messages, + "Generated dotplot for assembly-to-assembly comparison." +) + +readr::write_lines( + file = snakemake@log[["path"]], + x = paste0("DOTPLOT: ", messages) +) \ No newline at end of file From d3d18824675ae314de8ff13658cdb0eac923962a Mon Sep 17 00:00:00 2001 From: Rina Ahmed-Begrich Date: Thu, 10 Sep 2026 17:06:59 +0200 Subject: [PATCH 3/4] feat: add variant annotation with snpEff. restructuring output folders. --- .test/config/config.yml | 5 + config/config.yml | 7 +- config/schemas/config.schema.yml | 29 ++- resources/images/dag.svg | 210 +++++++++++------ workflow/Snakefile | 1 + workflow/envs/vcfutils.yml | 2 + workflow/rules/common.smk | 31 ++- workflow/rules/genome_comp.smk | 388 +++++++++++++++++++++++++++++++ workflow/rules/qc.smk | 317 +------------------------ 9 files changed, 592 insertions(+), 398 deletions(-) create mode 100644 workflow/rules/genome_comp.smk diff --git a/.test/config/config.yml b/.test/config/config.yml index b3fdaf7..ed8bdfb 100644 --- a/.test/config/config.yml +++ b/.test/config/config.yml @@ -61,3 +61,8 @@ reference_comparison: extra: "" bcftools: extra: "" + snpeff: + build: + extra: "-noCheckCds -noCheckProtein -noLog -q" + annotate: + extra: "-nodownload -ud 0 -noLog" diff --git a/config/config.yml b/config/config.yml index 42c448c..0177bf1 100644 --- a/config/config.yml +++ b/config/config.yml @@ -60,4 +60,9 @@ reference_comparison: paftools: extra: "" bcftools: - extra: "" \ No newline at end of file + extra: "" + snpeff: + build: + extra: "-noCheckCds -noCheckProtein -noLog -q" + annotate: + extra: "-nodownload -ud 0 -noLog" diff --git a/config/schemas/config.schema.yml b/config/schemas/config.schema.yml index d62199c..5304a51 100644 --- a/config/schemas/config.schema.yml +++ b/config/schemas/config.schema.yml @@ -244,11 +244,30 @@ properties: default: "" required: - extra - required: - - skip - - minimap2 - - paftools - - bcftools + snpeff: + type: object + properties: + build: + type: object + properties: + extra: + type: string + description: Extra command-line arguments for snpEff build + default: "-noCheckCds -noCheckProtein -noLog -q" + required: + - extra + annotate: + type: object + properties: + extra: + type: string + description: Extra command-line arguments for snpEff annotate + default: "-nodownload -ud 0 -noLog" + required: + - extra + required: + - build + - annotate required: - samplesheet - tool diff --git a/resources/images/dag.svg b/resources/images/dag.svg index 686e596..596fe72 100644 --- a/resources/images/dag.svg +++ b/resources/images/dag.svg @@ -4,49 +4,49 @@ - - + + snakemake_dag - + 0 - -all + +all 1 - + annotate_prokka - + 1->0 - - + + 2 - + annotate_pgap - + 2->0 - - + + 3 - + prepare_yaml_files - + 3->2 @@ -54,11 +54,11 @@ 4 - + get_pgap_fasta - + 4->3 @@ -66,23 +66,23 @@ 5 - + annotate_bakta - + 5->0 - - + + 6 - + get_bakta_db - + 6->5 @@ -90,86 +90,158 @@ 7 - + quast - + 7->0 - - + + 8 - -checkm + +rgi_detection - + 8->0 - - + + 9 - -get_ceckm_db + +viz_synteny - - -9->8 - - + + +9->0 + + 10 - -rgi_detection + +synteny_detection - - -10->0 - - + + +10->9 + + 11 - -viz_synteny + +prepare_names - - -11->0 - - + + +11->9 + + 12 - -synteny_detection + +merge_vcfs - - -12->11 - - + + +12->0 + + 13 - -prepare_names - - - -13->11 - - + +paftools_vcf + + + +13->12 + + + + + +16 + +snpeff + + + +13->16 + + + + + +14 + +minimap2_paf + + + +14->13 + + + + + +18 + +dotplot + + + +14->18 + + + + + +15 + +snpeff_vcf_to_tab + + + +15->0 + + + + + +16->15 + + + + + +17 + +snpeff_build_db + + + +17->16 + + + + + +18->0 + + diff --git a/workflow/Snakefile b/workflow/Snakefile index 03a3b50..a22994f 100644 --- a/workflow/Snakefile +++ b/workflow/Snakefile @@ -25,6 +25,7 @@ configfile: "config/config.yml" include: "rules/common.smk" include: "rules/annotate.smk" include: "rules/qc.smk" +include: "rules/genome_comp.smk" # ----------------------------------------------------- diff --git a/workflow/envs/vcfutils.yml b/workflow/envs/vcfutils.yml index ce51ca1..165e8ad 100644 --- a/workflow/envs/vcfutils.yml +++ b/workflow/envs/vcfutils.yml @@ -8,6 +8,8 @@ dependencies: - bcftools=1.24 - perl-vcftools-vcf=0.1.15 - snippy=4.0.2 + - snpeff=5.4.0c-0 + - snakemake-wrapper-utils=0.9.0 - htslib=1.24 - r-pafr=0.0.2 - r-readr=2.2.0 diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk index 006b5c4..c0f6786 100644 --- a/workflow/rules/common.smk +++ b/workflow/rules/common.smk @@ -44,7 +44,7 @@ def get_fasta_ntsynt(wildcards): def get_panaroo_gff(wildcards): return expand( - "results/qc/panaroo/{tool}/prepare/{sample}.gff", + "results/panaroo/{tool}/prepare/{sample}.gff", tool=wildcards.tool, sample=samples.index, ) @@ -52,7 +52,7 @@ def get_panaroo_gff(wildcards): def get_panaroo_fasta(wildcards): return expand( - "results/qc/panaroo/{tool}/prepare/{sample}.fna", + "results/panaroo/{tool}/prepare/{sample}.fna", tool=wildcards.tool, sample=samples.index, ) @@ -72,12 +72,12 @@ def get_final_input(wildcards): ) if len(samples.index) > 1 and not config["panaroo"]["skip"]: inputs += expand( - "results/qc/panaroo/{tool}/summary_statistics.txt", + "results/panaroo/{tool}/summary_statistics.txt", tool=config["tool"], ) if len(samples.index) > 1 and not config["fastani"]["skip"]: inputs += expand( - "results/qc/fastani/summary.txt", + "results/fastani/summary.txt", ) if not config["checkm"]["skip"]: inputs += expand( @@ -85,7 +85,7 @@ def get_final_input(wildcards): ) if not config["rgi"]["skip"]: inputs += expand( - "results/qc/rgi/{sample}/result.{ext}", + "results/rgi/{sample}/result.{ext}", sample=samples.index, ext=["txt", "json"], ) @@ -94,17 +94,22 @@ def get_final_input(wildcards): len(samples.index) == 1 and config["reference"]["fasta"] != "" ): inputs += expand( - "results/qc/genome_synteny/ntSynt-viz_ribbon-plot.pdf", + "results/genome_synteny/ntSynt-viz_ribbon-plot.pdf", ) if ( config["reference"]["fasta"] != "" and not config["reference_comparison"]["skip"] ): inputs += expand( - "results/qc/reference_comparison/all_merged_aln.vcf.gz", + "results/reference_comparison/all_merged_aln.vcf.gz", ) + if config["reference"]["gff"] != "": + inputs += expand( + "results/reference_comparison/annotated_vcf/{sample}_annotated.tab", + sample=samples.index, + ) inputs += expand( - "results/qc/reference_comparison/{sample}_dotplot.pdf", + "results/reference_comparison/{sample}_dotplot.pdf", sample=samples.index, ) return inputs @@ -132,3 +137,13 @@ def format_bakta_locustag(raw): f"\nlocustag '{raw}' converted to '{cleaned}' to meet BAKTA requirements (between 3 and 12 alphanumeric uppercase characters, start with a letter)\n" ) return cleaned + + +def get_chromosome(): + """Get the chromosome name from the reference fasta file.""" + if config["reference"]["fasta"]: + with open(config["reference"]["fasta"], "r") as f: + for line in f: + if line.startswith(">"): + return line[1:].strip() + return None diff --git a/workflow/rules/genome_comp.smk b/workflow/rules/genome_comp.smk new file mode 100644 index 0000000..29096ee --- /dev/null +++ b/workflow/rules/genome_comp.smk @@ -0,0 +1,388 @@ +rule prepare_panaroo: + input: + fasta="results/annotation/{tool}/{sample}/{sample}.fna", + gff="results/annotation/{tool}/{sample}/{sample}.gff", + output: + fasta="results/panaroo/{tool}/prepare/{sample}.fna", + gff="results/panaroo/{tool}/prepare/{sample}.gff", + log: + "results/panaroo/{tool}/prepare/{sample}.log", + conda: + "../envs/panaroo.yml" + params: + remove_source=config["panaroo"]["remove_source"], + remove_feature=config["panaroo"]["remove_feature"], + message: + """--- Prepare input files for pan-genome alignment ---""" + shell: + """ + echo 'Preparing annotation for Panaroo:' >{log} + echo ' - formatting seqnames in FASTA files' >>{log} + awk '{{ sub(/>.*\\|/, ">"); sub(/[[:space:]].*$/, ""); print }}' \ + {input.fasta} >{output.fasta} 2>>{log} + echo ' - removing sequences and selected features in GFF files' >>{log} + awk ' /^##FASTA/ {{exit}} $2 !~ /{params.remove_source}/ && $3 !~ /{params.remove_feature}/ {{print}}' \ + {input.gff} >{output.gff} 2>>{log} + """ + + +rule panaroo: + input: + gff=get_panaroo_gff, + fasta=get_panaroo_fasta, + output: + stats="results/panaroo/{tool}/summary_statistics.txt", + log: + "results/panaroo/{tool}/panaroo.log", + conda: + "../envs/panaroo.yml" + threads: max(workflow.cores * 0.5, 1) + params: + outdir=lambda wc, output: os.path.dirname(output.stats), + extra=config["panaroo"]["extra"], + message: + """--- Running PANAROO to create pangenome from all annotations ---""" + shell: + """ + printf '%s\n' {input.gff} \ + | paste -d ' ' - <(printf '%s\n' {input.fasta}) \ + >{params.outdir}/input_files.txt + panaroo \ + -i {params.outdir}/input_files.txt \ + -o {params.outdir} \ + -t {threads} \ + {params.extra} \ + >{log} 2>&1 + """ + + +rule fastani: + input: + fasta=get_all_fasta, + output: + txt="results/fastani/summary.txt", + log: + "results/fastani/fastani.log", + conda: + "../envs/fastani.yml" + threads: max(workflow.cores * 0.5, 1) + params: + outdir=lambda wc, output: os.path.dirname(output.txt), + ref_fasta=( + [config["reference"]["fasta"]] if config["reference"]["fasta"] else [] + ), + extra=config["fastani"]["extra"], + message: + """--- Running FastANI to compare genome similarity (all vs all) ---""" + shell: + """ + printf '%s\n' {input.fasta} >{params.outdir}/input_files.txt + printf '%s\n' {params.ref_fasta} >>{params.outdir}/input_files.txt + fastANI \ + --ql {params.outdir}/input_files.txt \ + --rl {params.outdir}/input_files.txt \ + --output {output.txt} \ + --threads {threads} \ + {params.extra} \ + >{log} 2>&1 + """ + + +rule synteny_detection: + input: + fastas=get_fasta_ntsynt, + output: + tsv="results/genome_synteny/ntSynt.synteny_blocks.tsv", + fai=directory("results/genome_synteny/fai"), + log: + "results/genome_synteny/logs/ntsynt.log", + conda: + "../envs/ntsynt.yml" + threads: workflow.cores + params: + outdir=lambda wc, output: os.path.dirname(output.tsv), + divergence=config["synteny"]["divergence"], + extra=config["synteny"]["extra"], + message: + """--- Running ntSynt for multi-genome macrosynteny synteny detection ---""" + shell: + """ + ntSynt {input.fastas} \ + -d {params.divergence} \ + -t {threads} \ + --force \ + --prefix ntSynt \ + {params.extra} \ + >{log} 2>&1 + echo "Synteny detection completed. Moving results to output directory." >>{log} + mv -f ./ntSynt.* {params.outdir}/ + echo "Create fai output directory." >>{log} + mkdir -p {output.fai} + mv -f ./*.fai {output.fai}/ + echo "Remove intermediate files." >>{log} + rm -f ./*.fai ./*.tsv ./*.bf ./*.dot + """ + + +rule prepare_names: + output: + "results/genome_synteny/ntsynt-viz_name_conversion.tsv", + log: + "results/genome_synteny/logs/prepare_ntsynt-viz_names.log", + conda: + "../envs/base.yml" + threads: 1 + params: + sample_sheet=config["samplesheet"], + ref=( + (config["reference"]["fasta"], config["reference"]["name"]) + if config["reference"]["fasta"] + else [] + ), + message: + """--- Preparing name mapping file for ntSynt visualization ---""" + script: + "../scripts/prepare_names.py" + + +rule viz_synteny: + input: + blocks=rules.synteny_detection.output.tsv, + fai=rules.synteny_detection.output.fai, + names=rules.prepare_names.output, + output: + pdf="results/genome_synteny/ntSynt-viz_ribbon-plot.pdf", + log: + "results/genome_synteny/logs/ntsynt-viz.log", + conda: + "../envs/ntsynt.yml" + threads: 1 + params: + outdir=lambda wc, output: os.path.dirname(output.pdf), + fais=lambda wc, input: " ".join(glob.glob(os.path.join(input.fai, "*.fai"))), + scale=config["synteny"]["viz_scale"], + ref_fasta=( + [ + "--target-genome", + config["reference"]["name"].replace("_", "-").replace(" ", "_"), + ] + if config["reference"]["fasta"] and config["reference"]["name"] + else ( + ["--target-genome", config["reference"]["fasta"]] + if config["reference"]["fasta"] + else [] + ) + ), + extra=config["synteny"]["viz_extra"], + message: + """--- Running ntSynt-viz to generate multi-genome ribbon plots ---""" + shell: + """ + ntsynt_viz.py \ + --blocks {input.blocks} \ + --fais {params.fais} \ + --name_conversion {input.names} \ + {params.ref_fasta} \ + --scale {params.scale} \ + --format pdf \ + --prefix ntSynt-viz \ + {params.extra} \ + >{log} 2>&1 + echo "Synteny-viz completed. Moving results to output directory." >>{log} + mv -f ./ntSynt.* {params.outdir}/ + mv -f ./ntSynt-viz.* {params.outdir}/ + mv -f ./ntSynt-viz_* {params.outdir}/ + echo "Clean intermediate files." >>{log} + rm -f ./ntSynt-viz.*.tsv ./ntSynt-viz_* ./ntSynt.*.tsv + """ + + +rule minimap2_paf: + input: + target=config["reference"]["fasta"], + query=get_fasta, + output: + "results/reference_comparison/{sample}_aln.paf", + log: + "results/reference_comparison/logs/{sample}_aln.log", + threads: max(workflow.cores * 0.25, 1) + params: + extra=config["reference_comparison"]["minimap2"]["extra"], + sorting=config["reference_comparison"]["minimap2"]["sorting"], + sort_extra=config["reference_comparison"]["minimap2"]["sort_extra"], + message: + """--- Running assembly-to-assembly mapping using minimap2 ---""" + wrapper: + "v9.7.0/bio/minimap2/aligner" + + +rule paftools_vcf: + input: + target=config["reference"]["fasta"], + query=rules.minimap2_paf.output, + output: + vcf="results/reference_comparison/{sample}_aln.vcf", + vcf_gz="results/reference_comparison/{sample}_aln.vcf.gz", + index="results/reference_comparison/{sample}_aln.vcf.gz.tbi", + log: + "results/reference_comparison/logs/{sample}_vcf.log", + conda: + "../envs/vcfutils.yml" + threads: max(workflow.cores * 0.25, 1) + params: + extra=config["reference_comparison"]["paftools"]["extra"], + message: + """--- Running assembly-to-assembly mapping using minimap2 ---""" + shell: + """ + sort -k6,6 -k8,8n {input.query} \ + | paftools.js call -f {input.target} - \ + {params.extra} 2>{log} \ + | bcftools reheader \ + --threads {threads} \ + -n {wildcards.sample} \ + | vcf-annotate --fill-type \ + >{output.vcf} 2>>{log} + bgzip -c {output.vcf} >{output.vcf_gz} 2>>{log} + tabix {output.vcf_gz} + """ + + +rule merge_vcfs: + input: + vcfs=expand( + "results/reference_comparison/{sample}_aln.vcf.gz", sample=samples.index + ), + output: + merged_vcf_gz="results/reference_comparison/all_merged_aln.vcf.gz", + index="results/reference_comparison/all_merged_aln.vcf.gz.tbi", + log: + "results/reference_comparison/logs/merge_vcfs.log", + conda: + "../envs/vcfutils.yml" + threads: max(workflow.cores * 0.5, 1) + params: + extra=config["reference_comparison"]["bcftools"]["extra"], + num_input=len(samples.index), + message: + """--- Merging individual VCF files into a single multi-sample VCF ---""" + shell: + """ + if [ {params.num_input} -eq 1 ]; then + echo "Only one VCF file present. Copying to merged output." >{log} + cp {input.vcfs} {output.merged_vcf_gz} + else + echo "Merging {params.num_input} VCF files using bcftools..." >{log} + bcftools merge \ + {input.vcfs} \ + {params.extra} \ + --output {output.merged_vcf_gz} \ + >>{log} 2>&1 + fi + tabix -p vcf {output.merged_vcf_gz} + """ + + +rule snpeff_build_db: + input: + gff=config["reference"]["gff"], + fasta=config["reference"]["fasta"], + output: + db=directory("results/reference_comparison/annotated_vcf/snpeff/custom_db/ref"), + log: + "results/reference_comparison/annotated_vcf/snpeff/logs/build_db.log", + conda: + "../envs/vcfutils.yml" + params: + genome_name=get_chromosome(), + extra=config["reference_comparison"]["snpeff"]["build"]["extra"], + message: + """--- Build custom SnpEff database ---""" + shell: + """ + mkdir -p {output.db} + config="$(realpath $(whereis snpEff | awk '{{print $2}}')).config" + today="$(date +%Y-%m-%d)" + cp {input.fasta} {output.db}/sequences.fa + cp {input.gff} {output.db}/genes.gff + cp ${{config}} {output.db}/../snpeff.config + #echo -e "\n# automatic entry by workflow: snakemake-assembly-postprocessing" >>{output.db}/../snpeff.config + echo -e "ref.genome : snakemake-assembly-postprocessing reference" >>{output.db}/../snpeff.config + echo -e "\tref.chromosome: {params.genome_name}\n" >>{output.db}/../snpeff.config + echo -e "\tref.{params.genome_name}.codonTable: Bacterial_and_Plant_Plastid\n" >>{output.db}/../snpeff.config + echo -e "\tref.retrieval_date : ${{today}}\n" >>{output.db}/../snpeff.config + snpEff build {params.extra} -c {output.db}/../snpeff.config -gff3 -dataDir $(realpath {output.db}/../) ref >{log} 2>&1 + """ + + +rule snpeff: + input: + calls=rules.paftools_vcf.output.vcf_gz, + db=rules.snpeff_build_db.output.db, + output: + calls="results/reference_comparison/annotated_vcf/{sample}_annotated.vcf", + stats="results/reference_comparison/annotated_vcf/snpeff/{sample}_snpeff.html", + csvstats="results/reference_comparison/annotated_vcf/snpeff/{sample}_snpeff.csv", + log: + "results/reference_comparison/annotated_vcf/snpeff/logs/{sample}_snpeff.log", + resources: + java_opts="-XX:ParallelGCThreads=2", + mem_mb=4096, + params: + extra=f"-c results/reference_comparison/annotated_vcf/snpeff/custom_db/snpeff.config {config["reference_comparison"]["snpeff"]["annotate"]["extra"]}", + message: + "annotate variants using snpEff" + wrapper: + "v9.15.0/bio/snpeff/annotate" + + +rule snpeff_vcf_to_tab: + input: + vcf=rules.snpeff.output.calls, + gff=config["reference"]["gff"], + target=config["reference"]["fasta"], + output: + tab="results/reference_comparison/annotated_vcf/{sample}_annotated.tab", + vcf_gz="results/reference_comparison/annotated_vcf/{sample}_annotated.vcf.gz", + index="results/reference_comparison/annotated_vcf/{sample}_annotated.vcf.gz.tbi", + log: + "results/reference_comparison/annotated_vcf/snpeff/logs/{sample}_snpeff_tab.log", + conda: + "../envs/vcfutils.yml" + threads: 1 + message: + """--- Converting SnpEff annotated VCF to tabular format ---""" + shell: + """ + snippy-vcf_to_tab \ + --ref {input.target} \ + --gff {input.gff} \ + --vcf {input.vcf} \ + >{output.tab} 2>{log} + bgzip -c {input.vcf} >{output.vcf_gz} 2>>{log} + tabix -p vcf {output.vcf_gz} + """ + + +rule dotplot: + input: + query=rules.minimap2_paf.output, + output: + pdf="results/reference_comparison/{sample}_dotplot.pdf", + png="results/reference_comparison/{sample}_dotplot.png", + log: + path="results/reference_comparison/logs/{sample}_dotplot.log", + conda: + "../envs/vcfutils.yml" + threads: 1 + params: + query=lambda wc: samples.loc[samples["sample"] == wc.sample, "strain"].values[0], + target=( + config["reference"]["name"] + if config["reference"]["name"] + else os.path.splitext(os.path.basename(config["reference"]["fasta"]))[0] + ), + message: + """--- Generating dotplot for assembly-to-assembly comparison ---""" + script: + "../scripts/create_dotplot.R" diff --git a/workflow/rules/qc.smk b/workflow/rules/qc.smk index 57b086c..a644252 100644 --- a/workflow/rules/qc.smk +++ b/workflow/rules/qc.smk @@ -91,103 +91,13 @@ rule checkm: """ -rule fastani: - input: - fasta=get_all_fasta, - output: - txt="results/qc/fastani/summary.txt", - log: - "results/qc/fastani/fastani.log", - conda: - "../envs/fastani.yml" - threads: max(workflow.cores * 0.5, 1) - params: - outdir=lambda wc, output: os.path.dirname(output.txt), - ref_fasta=( - [config["reference"]["fasta"]] if config["reference"]["fasta"] else [] - ), - extra=config["fastani"]["extra"], - message: - """--- Running FastANI to compare genome similarity (all vs all) ---""" - shell: - """ - printf '%s\n' {input.fasta} >{params.outdir}/input_files.txt - printf '%s\n' {params.ref_fasta} >>{params.outdir}/input_files.txt - fastANI \ - --ql {params.outdir}/input_files.txt \ - --rl {params.outdir}/input_files.txt \ - --output {output.txt} \ - --threads {threads} \ - {params.extra} \ - >{log} 2>&1 - """ - - -rule prepare_panaroo: - input: - fasta="results/annotation/{tool}/{sample}/{sample}.fna", - gff="results/annotation/{tool}/{sample}/{sample}.gff", - output: - fasta="results/qc/panaroo/{tool}/prepare/{sample}.fna", - gff="results/qc/panaroo/{tool}/prepare/{sample}.gff", - log: - "results/qc/panaroo/{tool}/prepare/{sample}.log", - conda: - "../envs/panaroo.yml" - params: - remove_source=config["panaroo"]["remove_source"], - remove_feature=config["panaroo"]["remove_feature"], - message: - """--- Prepare input files for pan-genome alignment ---""" - shell: - """ - echo 'Preparing annotation for Panaroo:' >{log} - echo ' - formatting seqnames in FASTA files' >>{log} - awk '{{ sub(/>.*\\|/, ">"); sub(/[[:space:]].*$/, ""); print }}' \ - {input.fasta} >{output.fasta} 2>>{log} - echo ' - removing sequences and selected features in GFF files' >>{log} - awk ' /^##FASTA/ {{exit}} $2 !~ /{params.remove_source}/ && $3 !~ /{params.remove_feature}/ {{print}}' \ - {input.gff} >{output.gff} 2>>{log} - """ - - -rule panaroo: - input: - gff=get_panaroo_gff, - fasta=get_panaroo_fasta, - output: - stats="results/qc/panaroo/{tool}/summary_statistics.txt", - log: - "results/qc/panaroo/{tool}/panaroo.log", - conda: - "../envs/panaroo.yml" - threads: max(workflow.cores * 0.5, 1) - params: - outdir=lambda wc, output: os.path.dirname(output.stats), - extra=config["panaroo"]["extra"], - message: - """--- Running PANAROO to create pangenome from all annotations ---""" - shell: - """ - printf '%s\n' {input.gff} \ - | paste -d ' ' - <(printf '%s\n' {input.fasta}) \ - >{params.outdir}/input_files.txt - panaroo \ - -i {params.outdir}/input_files.txt \ - -o {params.outdir} \ - -t {threads} \ - {params.extra} \ - >{log} 2>&1 - """ - - rule rgi_detection: input: fasta=get_fasta, output: - multiext("results/qc/rgi/{sample}/result", ".txt", ".json"), + multiext("results/rgi/{sample}/result", ".txt", ".json"), log: - "results/qc/rgi/{sample}/result.log", + "results/rgi/{sample}/result.log", threads: max(workflow.cores * 0.25, 1) params: input_type="contig", @@ -196,226 +106,3 @@ rule rgi_detection: """--- Running RGI to detect antibiotic resistance genes ---""" wrapper: "https://raw.githubusercontent.com/MPUSP/mpusp-snakemake-wrappers/refs/heads/main/rgi" - - -rule synteny_detection: - input: - fastas=get_fasta_ntsynt, - output: - tsv="results/qc/genome_synteny/ntSynt.synteny_blocks.tsv", - fai=directory("results/qc/genome_synteny/fai"), - log: - "results/qc/genome_synteny/logs/ntsynt.log", - conda: - "../envs/ntsynt.yml" - threads: workflow.cores - params: - outdir=lambda wc, output: os.path.dirname(output.tsv), - divergence=config["synteny"]["divergence"], - extra=config["synteny"]["extra"], - message: - """--- Running ntSynt for multi-genome macrosynteny synteny detection ---""" - shell: - """ - ntSynt {input.fastas} \ - -d {params.divergence} \ - -t {threads} \ - --force \ - --prefix ntSynt \ - {params.extra} \ - >{log} 2>&1 - echo "Synteny detection completed. Moving results to output directory." >>{log} - mv -f ./ntSynt.* {params.outdir}/ - echo "Create fai output directory." >>{log} - mkdir -p {output.fai} - mv -f ./*.fai {output.fai}/ - echo "Remove intermediate files." >>{log} - rm -f ./*.fai ./*.tsv ./*.bf ./*.dot - """ - - -rule prepare_names: - output: - "results/qc/genome_synteny/ntsynt-viz_name_conversion.tsv", - log: - "results/qc/genome_synteny/logs/prepare_ntsynt-viz_names.log", - conda: - "../envs/base.yml" - threads: 1 - params: - sample_sheet=config["samplesheet"], - ref=( - (config["reference"]["fasta"], config["reference"]["name"]) - if config["reference"]["fasta"] - else [] - ), - message: - """--- Preparing name mapping file for ntSynt visualization ---""" - script: - "../scripts/prepare_names.py" - - -rule viz_synteny: - input: - blocks=rules.synteny_detection.output.tsv, - fai=rules.synteny_detection.output.fai, - names=rules.prepare_names.output, - output: - pdf="results/qc/genome_synteny/ntSynt-viz_ribbon-plot.pdf", - log: - "results/qc/genome_synteny/logs/ntsynt-viz.log", - conda: - "../envs/ntsynt.yml" - threads: 1 - params: - outdir=lambda wc, output: os.path.dirname(output.pdf), - fais=lambda wc, input: " ".join(glob.glob(os.path.join(input.fai, "*.fai"))), - scale=config["synteny"]["viz_scale"], - ref_fasta=( - [ - "--target-genome", - config["reference"]["name"].replace("_", "-").replace(" ", "_"), - ] - if config["reference"]["fasta"] and config["reference"]["name"] - else ( - ["--target-genome", config["reference"]["fasta"]] - if config["reference"]["fasta"] - else [] - ) - ), - extra=config["synteny"]["viz_extra"], - message: - """--- Running ntSynt-viz to generate multi-genome ribbon plots ---""" - shell: - """ - ntsynt_viz.py \ - --blocks {input.blocks} \ - --fais {params.fais} \ - --name_conversion {input.names} \ - {params.ref_fasta} \ - --scale {params.scale} \ - --format pdf \ - --prefix ntSynt-viz \ - {params.extra} \ - >{log} 2>&1 - echo "Synteny-viz completed. Moving results to output directory." >>{log} - mv -f ./ntSynt.* {params.outdir}/ - mv -f ./ntSynt-viz.* {params.outdir}/ - mv -f ./ntSynt-viz_* {params.outdir}/ - echo "Clean intermediate files." >>{log} - rm -f ./ntSynt-viz.*.tsv ./ntSynt-viz_* ./ntSynt.*.tsv - """ - - -rule minimap2_paf: - input: - target=config["reference"]["fasta"], - query=get_fasta, - output: - "results/qc/reference_comparison/{sample}_aln.paf", - log: - "results/qc/reference_comparison/logs/{sample}_aln.log", - threads: max(workflow.cores * 0.25, 1) - params: - extra=config["reference_comparison"]["minimap2"]["extra"], - sorting=config["reference_comparison"]["minimap2"]["sorting"], - sort_extra=config["reference_comparison"]["minimap2"]["sort_extra"], - message: - """--- Running assembly-to-assembly mapping using minimap2 ---""" - wrapper: - "v9.7.0/bio/minimap2/aligner" - - -rule paftools_vcf: - input: - target=config["reference"]["fasta"], - query=rules.minimap2_paf.output, - output: - tab="results/qc/reference_comparison/{sample}_aln.tab", - vcf="results/qc/reference_comparison/{sample}_aln.vcf", - vcf_gz="results/qc/reference_comparison/{sample}_aln.vcf.gz", - index="results/qc/reference_comparison/{sample}_aln.vcf.gz.tbi", - log: - "results/qc/reference_comparison/logs/{sample}_vcf.log", - conda: - "../envs/vcfutils.yml" - threads: max(workflow.cores * 0.25, 1) - params: - extra=config["reference_comparison"]["paftools"]["extra"], - message: - """--- Running assembly-to-assembly mapping using minimap2 ---""" - shell: - """ - sort -k6,6 -k8,8n {input.query} \ - | paftools.js call -f {input.target} - \ - {params.extra} 2>{log} \ - | bcftools reheader \ - --threads {threads} \ - -n {wildcards.sample} \ - | vcf-annotate --fill-type \ - >{output.vcf} 2>>{log} - snippy-vcf_to_tab --ref {input.target} \ - --vcf {output.vcf} \ - >{output.tab} 2>>{log} - bgzip -c {output.vcf} >{output.vcf_gz} 2>>{log} - tabix {output.vcf_gz} - """ - - -rule merge_vcfs: - input: - vcfs=expand( - "results/qc/reference_comparison/{sample}_aln.vcf.gz", sample=samples.index - ), - output: - merged_vcf_gz="results/qc/reference_comparison/all_merged_aln.vcf.gz", - index="results/qc/reference_comparison/all_merged_aln.vcf.gz.tbi", - log: - "results/qc/reference_comparison/logs/merge_vcfs.log", - conda: - "../envs/vcfutils.yml" - threads: max(workflow.cores * 0.5, 1) - params: - extra=config["reference_comparison"]["bcftools"]["extra"], - num_input=len(samples.index), - message: - """--- Merging individual VCF files into a single multi-sample VCF ---""" - shell: - """ - if [ {params.num_input} -eq 1 ]; then - echo "Only one VCF file present. Copying to merged output." >{log} - cp {input.vcfs} {output.merged_vcf_gz} - else - echo "Merging {params.num_input} VCF files using bcftools..." >{log} - bcftools merge \ - {input.vcfs} \ - {params.extra} \ - --output {output.merged_vcf_gz} \ - >>{log} 2>&1 - fi - tabix {output.merged_vcf_gz} - """ - - -rule dotplot: - input: - query=rules.minimap2_paf.output, - output: - pdf="results/qc/reference_comparison/{sample}_dotplot.pdf", - png="results/qc/reference_comparison/{sample}_dotplot.png", - log: - path="results/qc/reference_comparison/logs/{sample}_dotplot.log", - conda: - "../envs/vcfutils.yml" - threads: 1 - params: - query=lambda wc: samples.loc[samples["sample"] == wc.sample, "strain"].values[0], - target=( - config["reference"]["name"] - if config["reference"]["name"] - else os.basename(config["reference"]["fasta"]) - ), - message: - """--- Generating dotplot for assembly-to-assembly comparison ---""" - script: - "../scripts/create_dotplot.R" From e90e8f0470dfcc8ccc0461c5c93f233045262a4c Mon Sep 17 00:00:00 2001 From: Rina Ahmed-Begrich Date: Fri, 11 Sep 2026 14:34:50 +0200 Subject: [PATCH 4/4] chore: update readme. --- README.md | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 63afa8a..2cc3a43 100644 --- a/README.md +++ b/README.md @@ -30,15 +30,18 @@ _Workflow overview:_ ## Workflow overview -1. Parse `samples.csv` table containing the samples's meta data (`python`) +1. Parse `samples.csv` table containing the samples's meta data (`python`). 2. Annotate assemblies using one of the following tools: - 1. NCBI's Prokaryotic Genome Annotation Pipeline ([PGAP](https://github.com/ncbi/pgap)). Note: needs to be installed manually - 2. [prokka](https://github.com/tseemann/prokka), a fast and light-weight prokaryotic annotation tool - 3. [bakta](https://github.com/oschwengers/bakta), a fast, alignment-free annotation tool. Note: Bakta will automatically download its companion database from zenodo (light: 1.5 GB, full: 40 GB) -3. Predict antimicrobial resistance (AMR) genes using [RGI](https://github.com/arpcard/rgi) -4. Create a QC report for the assemblies using [Quast](https://github.com/ablab/quast) -5. Create a pangenome analysis (orthologs/homologs) using [Panaroo](https://gthlab.au/panaroo/) + 1. NCBI's Prokaryotic Genome Annotation Pipeline ([PGAP](https://github.com/ncbi/pgap)). Note: needs to be installed manually. + 2. [prokka](https://github.com/tseemann/prokka), a fast and light-weight prokaryotic annotation tool. + 3. [bakta](https://github.com/oschwengers/bakta), a fast, alignment-free annotation tool. Note: Bakta will automatically download its companion database from zenodo (light: 1.5 GB, full: 40 GB). +3. Predict antimicrobial resistance (AMR) genes using [RGI](https://github.com/arpcard/rgi). +4. Create a QC report for the assemblies using [Quast](https://github.com/ablab/quast). +5. Create a pangenome analysis (orthologs/homologs) using [Panaroo](https://gthlab.au/panaroo/). 6. Compute pairwise average nucleotide identity (ANI) between the assemblies using [FastANI](https://github.com/ParBLiSS/FastANI) and plot a phylogenetic tree based on the ANI distances. +7. Estimate genome completeness and contamination with [checkM2](https://github.com/chklovski/CheckM2). +8. Detect and visualize multi-genome synteny with [ntSynt](https://github.com/BirolLab/ntSynt) and [ntSynt-viz](https://github.com/BirolLab/ntSynt-viz). +9. Compute pairwise whole-genome reference comparison and annotate variant effects with [minimap2](https://github.com/lh3/minimap2), [vcftools](https://github.com/vcftools/vcftools), [bcftools](https://github.com/samtools/bcftools) and [snpEff](https://github.com/pcingola/snpeff). ## Installation @@ -119,6 +122,6 @@ snakemake --cores 2 --sdm conda apptainer --directory .test > Köster J., Mölder F., Jablonski K. P., Letcher B., Hall M. B., Tomkins-Tinch C. H., Sochat V., Forster J., Lee S., Twardziok S. O., Kanitz A., Wilm A., Holtgrewe M., Rahmann S., & Nahnsen S. _Sustainable data analysis with Snakemake_. F1000Research, 10:33, 10, 33, **2021**. https://doi.org/10.12688/f1000research.29032.2. -> Coombe L, Kazemi P, Wong J, Birol I, Warren RL. _ntSynt: multi-genome synteny detection using minimizer graph mappings_. BMC Biology. 23:367, **2025**. https://doi.org/10.1186/s12915-025-02455-w +> Coombe L, Kazemi P, Wong J, Birol I, Warren RL. _ntSynt: multi-genome synteny detection using minimizer graph mappings_. BMC Biology., 23:367, **2025**. https://doi.org/10.1186/s12915-025-02455-w -> Coombe L, Warren RL, Birol I. _ntSynt-viz: Visualizing synteny patterns across multiple genomes_. bioRxiv 2025.01.15.633221. https://doi.org/10.1101/2025.01.15.633221 \ No newline at end of file +> Coombe L, Warren RL, Birol I. _ntSynt-viz: Visualizing synteny patterns across multiple genomes_. J. Evol. Biol., **2026**. https://doi.org/10.1093/jeb/voag079 \ No newline at end of file