Marginal reconstructions sometimes use joint ancestral state reconstruction

[jameshadfield]

While exploring how marginal reconstruction worked in the various pieces in TreeTime I noticed that branch_length_mode='marginal' would sometimes perform joint ancestral reconstructions. I presume this is a bug!

1. TreeTime.run uses the 'marginal_sequences' key:

treetime/treetime/treetime.py

Lines 210 to 212 in df6d2c0

	# determine how to reconstruct and sample sequences
	seq_kwargs = {
	'marginal_sequences': sequence_marginal or (self.branch_length_mode == 'marginal'),

2. The first time we call TreeAnc.infer_ancestral_sequences we remap this to the marginal argument that function expects:

treetime/treetime/treetime.py

Lines 237 to 239 in df6d2c0

	self.infer_ancestral_sequences(
	infer_gtr=infer_gtr, marginal=seq_kwargs['marginal_sequences'], **seq_kwargs
	)

3. Later on we just pass seq_kwargs to infer_ancestral_sequences, and thus the reconstruction runs with the default marginal=False parameters, which means "joint" in all the testing I've done

treetime/treetime/treetime.py

Lines 342 to 349 in df6d2c0

	if need_new_time_tree:
	self.make_time_tree(**tt_kwargs)
	if self.aln:
	ndiff = self.infer_ancestral_sequences('ml', **seq_kwargs)
	else: # no refinements, just iterate
	if self.aln:
	ndiff = self.infer_ancestral_sequences('ml', **seq_kwargs)
	self.make_time_tree(**tt_kwargs)

Suggested fix: chance the call signature of (3) to look more like (2)

---

Some quick and dirty logging which helped me track this down:

diff --git a/treetime/clock_tree.py b/treetime/clock_tree.py
index 8e51266c..ceef2899 100644
--- a/treetime/treeanc.py
+++ b/treetime/treeanc.py
@@ -546,6 +546,9 @@ class TreeAnc(object):
         self.logger(
             'TreeAnc.infer_ancestral_sequences with method: %s, %s' % (method, 'marginal' if marginal else 'joint'), 1
         )
+        
+        self.logger(f"TreeAnc.infer_ancestral_sequences. Marginal arg: {marginal}. Method = {method}", 1)
+
 
         if method.lower() in ['ml', 'probabilistic']:
             if marginal:
diff --git a/treetime/treetime.py b/treetime/treetime.py
index ea72c2b2..7fc74be2 100644
--- a/treetime/treetime.py
+++ b/treetime/treetime.py
@@ -339,6 +339,8 @@ class TreeTime(ClockTree):
                 assign_gamma(self.tree)
                 need_new_time_tree = True
 
+            self.logger(f"TreeTime::run marginal={seq_kwargs.get('marginal', '(key missing)')}  marginal_sequences={seq_kwargs.get('marginal_sequences', '(key missing)')}", 1)
+
             if need_new_time_tree:
                 self.make_time_tree(**tt_kwargs)
                 if self.aln:

[ivan-aksamentov]

Fixed in #602, released in 0.12.1

[jameshadfield]

Thanks for the quick fix ⭐