fix branch lengths in conversion scripts

ktmeaton · Apr 27, 2021 · c309726 · c309726
1 parent 8fba8a6
commit c309726
Show file tree

Hide file tree

Showing 4 changed files with 33 additions and 18 deletions.
diff --git a/workflow/rules/phylogeny.smk b/workflow/rules/phylogeny.smk
@@ -158,6 +158,6 @@ rule beast:
         cut -f 1,21,22 {input.tsv} > {output.latlon};
         cp {input.aln} {output.aln};
         cp {input.divtree} {output.divtree};
-        cp {input.timetree} {output.timetree};
+        python3 {scripts_dir}/nexus2newick.py {input.timetree} {output.timetree};
         cp {input.constant_sites} {output.constant_sites};
         """
diff --git a/workflow/scripts/newick2nexus.py b/workflow/scripts/newick2nexus.py
@@ -6,4 +6,20 @@
 tree_output = sys.argv[2]
 
 tree = Phylo.read(tree_input, "newick")
-Phylo.write(tree, tree_output, "nexus")
+
+# Get minimum branch length for output precision
+min_branch_length = 1e10
+
+for c in tree.find_clades():
+    # Remove comments
+    c.comment = None
+    if c.branch_length and c.branch_length < min_branch_length:
+        min_branch_length = c.branch_length
+
+# This is a terrible way to get sig digits
+sig_dig = 5
+# If we're using scientific notation, it's really small
+if "e-" in str(min_branch_length):
+    sig_dig = int(str(min_branch_length).split("e-")[-1]) + 2
+
+Phylo.write(tree, tree_output, "nexus", format_branch_length="%1.{}f".format(sig_dig))
diff --git a/workflow/scripts/nexus2beast.py b/workflow/scripts/nexus2beast.py
diff --git a/workflow/scripts/nexus2newick.py b/workflow/scripts/nexus2newick.py
@@ -7,8 +7,21 @@
 
 tree = Phylo.read(tree_input, "nexus")
 
-# Remove comments
+# Get minimum branch length for output precision
+min_branch_length = 1e10
+
 for c in tree.find_clades():
+    # Remove comments
     c.comment = None
+    if c.branch_length and c.branch_length < min_branch_length:
+        min_branch_length = c.branch_length
+
+# This is a terrible way to get sig digits
+sig_dig = 5
+# If we're using scientific notation, it's really small
+if "e-" in str(min_branch_length):
+    sig_dig = int(str(min_branch_length).split("e-")[-1]) + 2
 
-Phylo.write(tree, tree_output, "newick")
+Phylo.write(
+    tree, tree_output, "newick", format_branch_length="%1.{}f".format(sig_dig),
+)