Better use of ambiguous sites in iqtree

Snakefile

@@ -141,7 +141,7 @@ rule ska_align:
     conda:
         config["poppipe_location"] + "/envs/ska.yml"
     shell:
-        "ska align -v --filter no-filter {input.skf} > {output.alignment} 2> {log}"
+        "ska align -v --filter no-const --no-gap-only-sites {input.skf} > {output.alignment} 2> {log}"
 
 # ska for mapping (needed for gubbins)
 rule ska_map:
@@ -206,7 +206,7 @@ rule gubbins:
     conda:
         config["poppipe_location"] + "/envs/gubbins.yml"
     threads:
-        4
+        2
     shell:
         "pushd output/strains/{wildcards.strain}/ && \
         run_gubbins.py map_variants.aln --prefix {params.prefix} \
@@ -231,8 +231,6 @@ rule bactdating:
         "logs/bactdating_{strain}.log"
     conda:
         config["poppipe_location"] + "/envs/transphylo.yml"
-    threads:
-        4
     shell:
         "Rscript --vanilla {params.script} output/strains/{wildcards.strain}/gubbins {input.metadata} {output.sorted} {output.rds} > {log}"
 
@@ -269,8 +267,6 @@ rule transphylo:
         "logs/transphylo_{strain}.log"
     conda:
         config["poppipe_location"] + "/envs/transphylo.yml"
-    threads:
-        4
     shell:
         "Rscript --vanilla {params.script} --rds {input.rds} --sorted {input.sorted} --output {output} \
         --gubbins {params.gubbins} --wshape {params.w_shape} \

config.yml

@@ -22,7 +22,7 @@ ska:
 iqtree:
   enabled: true
   mode: full
-  model: GTR+G # do not use ASC
+  model: 012310+G+ASC # +ASC recommended. May replace 012310 with GTR for single strand data
 
 # fastbaps options
 fastbaps:

scripts/run_iqtree.py

@@ -10,15 +10,21 @@ def midpoint_root(infile, outfile):
     t.set_outgroup(t.get_midpoint_outgroup())
     t.write(format=5, outfile=outfile) # format 5: internal and leaf branches + leaf names
 
+def iqtree_cmd(alignment, start_tree, threads, prefix, model, mode):
+    cmd = f"iqtree --quiet -st DNA -s {alignment} -t {start_tree} -T {threads} --prefix {prefix}"
+    if mode == "full":
+        cmd += f" -m {model}"
+    elif mode == "fast":
+        cmd += " --fast"
+    return cmd
+
 if snakemake.params.enabled:
-    iqtree_cmd = "iqtree --quiet -s " + snakemake.input.alignment + " -t " + snakemake.input.start_tree + \
-                 " -T " + str(snakemake.threads) + " --prefix " + snakemake.params.prefix
-    if snakemake.params.mode == "full":
-        iqtree_cmd += " -m " + snakemake.params.model
-    elif snakemake.params.mode == "fast":
-        iqtree_cmd += " --fast"
-
-    subprocess.run(iqtree_cmd, shell=True, check=True)
+    cmd = iqtree_cmd(snakemake.input.alignment, snakemake.input.start_tree, str(snakemake.threads), snakemake.params.prefix, snakemake.params.model, snakemake.params.mode)
+    ret_code = subprocess.run(cmd, shell=True)
+    # Automatically retry with varsites if ASC fails
+    if ret_code.returncode != 0:
+        cmd = iqtree_cmd(snakemake.params.prefix + ".varsites.phy", snakemake.input.start_tree, str(snakemake.threads), snakemake.params.prefix, snakemake.params.model, snakemake.params.mode)
+        subprocess.run(cmd, shell=True, check=True)
 else:
     copyfile(snakemake.input.start_tree, snakemake.output.unrooted)