Match seed samples without recombination

CHANGELOG.md

@@ -4,6 +4,11 @@
 
 In development
 
+- Seed samples listed in `include_samples` are now matched separately from the
+  daily samples, with recombination effectively disallowed. These samples are
+  usually far diverged from the current ARG, which previously led to spurious
+  recombinations.
+
 - Add basic support for non-SARS-CoV-2 genomes via an optional reference FASTA.
   Supply `--reference` to `import-alignments` and a `reference_fasta` key in the
   inference config; both default to the built-in SARS-CoV-2 reference, so

sc2ts/inference.py

@@ -33,6 +33,10 @@
 MISSING = -1
 DELETION = core.IUPAC_ALLELES.index("-")
 
+# Value of num_mismatches which effectively disallows recombination during
+# matching, since solve_num_mismatches drives rho down to its 1e-200 floor.
+NO_RECOMBINATION_NUM_MISMATCHES = 1000
+
 
 def get_provenance_dict(command, args, start_time):
     document = {
@@ -675,21 +679,37 @@ def _extend(
             logger.info(f"Subset from {len(samples)} to {max_daily_samples}")
             samples = rng.sample(samples, max_daily_samples)
 
-    samples = samples + unconditional_include_samples
-    samples.sort(key=lambda s: s.strain)
-
     ts = increment_time(date, base_ts)
-    if len(samples) > 0:
-        match_samples(
-            date,
-            samples,
-            base_ts=base_ts,
-            num_mismatches=num_mismatches,
-            deletions_as_missing=deletions_as_missing,
-            show_progress=show_progress,
-            num_threads=num_threads,
-            memory_limit=memory_limit,
-        )
+    if len(samples) + len(unconditional_include_samples) > 0:
+        if len(samples) > 0:
+            match_samples(
+                date,
+                samples,
+                base_ts=base_ts,
+                num_mismatches=num_mismatches,
+                deletions_as_missing=deletions_as_missing,
+                show_progress=show_progress,
+                num_threads=num_threads,
+                memory_limit=memory_limit,
+            )
+        if len(unconditional_include_samples) > 0:
+            # Seed samples are usually far diverged from the current ARG, and
+            # matching them with the standard num_mismatches gives spurious
+            # recombinations. Match them without recombination instead.
+            match_samples(
+                date,
+                unconditional_include_samples,
+                base_ts=base_ts,
+                num_mismatches=NO_RECOMBINATION_NUM_MISMATCHES,
+                deletions_as_missing=deletions_as_missing,
+                show_progress=show_progress,
+                num_threads=num_threads,
+                memory_limit=memory_limit,
+            )
+
+        samples = samples + unconditional_include_samples
+        samples.sort(key=lambda s: s.strain)
+
         characterise_match_mutations(base_ts, samples)
         characterise_recombinants(base_ts, samples)
 
@@ -2075,7 +2095,7 @@ def rematch_recombinant(base_ts, recomb_ts, node_id, num_mismatches):
     match_tsinfer(
         samples=[sample],
         ts=base_ts,
-        num_mismatches=1000,
+        num_mismatches=NO_RECOMBINATION_NUM_MISMATCHES,
         mismatch_threshold=2 * original_cost,
     )
     result.no_recomb_match = sample.hmm_match

tests/test_inference.py

@@ -577,6 +577,13 @@ def test_2020_02_02_include_samples(
         u = ts.samples()[ts.metadata["sc2ts"]["samples_strain"].index("SRR11597115")]
         assert ts.nodes_flags[u] & tskit.NODE_IS_SAMPLE > 0
         assert ts.nodes_flags[u] & sc2ts.NODE_IS_UNCONDITIONALLY_INCLUDED > 0
+        # Seed samples are matched without recombination, so the node must
+        # have a single parent edge spanning the full sequence.
+        assert ts.nodes_flags[u] & sc2ts.NODE_IS_RECOMBINANT == 0
+        edges = [e for e in ts.edges() if e.child == u]
+        assert len(edges) == 1
+        assert edges[0].left == 0
+        assert edges[0].right == ts.sequence_length
 
     def test_2020_02_02_mutation_overlap(
         self,