feat: fix panaroo bug, added fastANI rule, harmonized multi-threading

Author: m-jahn

.test/config/config.yml

@@ -27,3 +27,7 @@ panaroo:
   remove_source: "cmsearch"
   remove_feature: "tRNA|rRNA|ncRNA|exon|sequence_feature"
   extra: "--clean-mode strict --remove-invalid-genes"
+
+fastani:
+  skip: False
+  extra: ""

config/README.md

@@ -9,6 +9,7 @@ A Snakemake workflow for the post-processing of microbial genome assemblies.
    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. Create a QC report for the assemblies using [Quast](https://github.com/ablab/quast)
 4. Create a pangenome analysis (orthologs/homologs) using [Panaroo](https://gthlab.au/panaroo/)
+5. 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.
 
 ## Running the workflow
 
@@ -22,4 +23,4 @@ The samplesheet table has the following layout:
 | EC2224 | "Streptococcus pyogenes" | SF370  | SPY       | assembly.fasta |
 | ...    | ...                      | ...    | ...       | ...            |
 
-**Note:** Pangenome analysis with `Panaroo` requires at least two samples.
+**Note:** Pangenome analysis with `Panaroo` and pairwise similarity analysis with `FastANI` requires at least two samples.

config/config.yml

@@ -27,3 +27,7 @@ panaroo:
   remove_source: "cmsearch"
   remove_feature: "tRNA|rRNA|ncRNA|exon|sequence_feature"
   extra: "--clean-mode strict --remove-invalid-genes"
+
+fastani:
+  skip: False
+  extra: ""

config/schemas/config.schema.yml

@@ -119,7 +119,17 @@ properties:
         type: string
         description: Extra command-line arguments for Panaroo
         default: "--clean-mode strict --remove-invalid-genes"
-
+  fastani:
+    type: object
+    properties:
+      skip:
+        type: boolean
+        description: Whether to skip FastANI analysis
+        default: false
+      extra:
+        type: string
+        description: Extra command-line arguments for FastANI
+        default: ""
 required:
   - samplesheet
   - tool
@@ -128,3 +138,4 @@ required:
   - bakta
   - quast
   - panaroo
+  - fastani

workflow/envs/fastani.yml

@@ -0,0 +1,7 @@
+name: panaroo
+channels:
+  - conda-forge
+  - bioconda
+  - nodefaults
+dependencies:
+  - fastani=1.34

workflow/envs/panaroo.yml

@@ -6,4 +6,5 @@ channels:
 dependencies:
   - numpy=1.26.4
   - scipy=1.11.4
-  - panaroo=1.5.2
+  - biopython=1.84
+  - panaroo=1.6.0

workflow/rules/annotate.smk

@@ -48,7 +48,6 @@ rule annotate_pgap:
         "results/annotation/pgap/logs/{sample}_pgap.log",
     conda:
         "../envs/base.yml"
-    threads: 1
     params:
         pgap=config["pgap"]["bin"],
         use_yaml_config=config["pgap"]["use_yaml_config"],
@@ -87,7 +86,7 @@ rule annotate_prokka:
         "results/annotation/prokka/logs/{sample}_prokka.log",
     conda:
         "../envs/prokka.yml"
-    threads: workflow.cores * 0.25
+    threads: max(workflow.cores * 0.5, 1)
     params:
         prefix=lambda wc: wc.sample,
         locustag=lambda wc: samples.loc[wc.sample]["id_prefix"],
@@ -127,7 +126,7 @@ rule get_bakta_db:
         "results/annotation/bakta/database/db.log",
     conda:
         "../envs/bakta.yml"
-    threads: workflow.cores * 0.25
+    threads: max(workflow.cores * 0.25, 1)
     params:
         download_db=config["bakta"]["download_db"],
         existing_db=config["bakta"]["existing_db"],
@@ -160,7 +159,7 @@ rule annotate_bakta:
         "results/annotation/bakta/logs/{sample}_bakta.log",
     conda:
         "../envs/bakta.yml"
-    threads: workflow.cores * 0.25
+    threads: max(workflow.cores * 0.25, 1)
     params:
         prefix=lambda wc: wc.sample,
         locustag=lambda wc: format_bakta_locustag(samples.loc[wc.sample]["id_prefix"]),

workflow/rules/common.smk

@@ -62,6 +62,11 @@ def get_final_input(wildcards):
             "results/qc/panaroo/{tool}/summary_statistics.txt",
             tool=config["tool"],
         )
+    if len(samples.index) > 1 and not config["fastani"]["skip"]:
+        inputs += expand(
+            "results/qc/fastani/{tool}/summary.txt",
+            tool=config["tool"],
+        )
     return inputs
 
 

workflow/rules/qc.smk

@@ -7,7 +7,7 @@ rule quast:
         "results/qc/quast/{tool}/quast.log",
     conda:
         "../envs/quast.yml"
-    threads: 4
+    threads: max(workflow.cores * 0.5, 1)
     params:
         outdir=lambda wc, output: os.path.dirname(output.report),
         ref_fasta=(
@@ -36,6 +36,40 @@ rule quast:
         """
 
 
+rule fastani:
+    input:
+        fasta=get_quast_fasta,
+    output:
+        txt="results/qc/fastani/{tool}/summary.txt",
+    log:
+        "results/qc/fastani/{tool}/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=(
+            " ".join(["-r", config["quast"]["reference_fasta"]])
+            if config["quast"]["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;
+        {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",
@@ -74,7 +108,7 @@ rule panaroo:
         "results/qc/panaroo/{tool}/panaroo.log",
     conda:
         "../envs/panaroo.yml"
-    threads: 4
+    threads: max(workflow.cores * 0.5, 1)
     params:
         outdir=lambda wc, output: os.path.dirname(output.stats),
         extra=config["panaroo"]["extra"],