Seeding by tree

.gitignore

@@ -3,8 +3,14 @@
 *.tsz
 *.nwk
 *.nex
+# Generated on demand alongside a FASTA
+*.fai
 _version.py
 
+# Editor swap files
+*.swp
+*.swo
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

CHANGELOG.md

@@ -4,10 +4,27 @@
 
 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 `seed_groups`, replacing the previous `include_samples` option. It is a
+  list of seed *groups*, each a list of strain IDs that are inserted into the
+  ARG together as a single local tree. A tree is inferred over the group's
+  haplotypes, the haplotype of the group's inferred ancestor is matched against
+  the ARG with recombination disallowed, and the whole group is then attached at
+  that single placement. Seed samples are usually far diverged from the current
+  ARG, so matching each one individually gave spurious recombinations or matches
+  to derived samples with large numbers of reversions; the inferred ancestor is
+  much closer to the contemporaneous ARG. A group is inserted on the minimum
+  date over its members, and members with later dates keep their real dates via
+  negative ("in the future") node times. The previous `include_samples` formats
+  (bare strain IDs, or `(strain, match_date)` tuples) are no longer supported.
+
+  Note that the `hmm_match` metadata of a seed node records the match that placed
+  its *group*, not a per-sample match, so its mutation count is not the number of
+  mutations between that sample and its parent.
+
+- Fix retrospective matching on a day with no new samples. The table of samples
+  already in the ARG was only rebuilt when there were samples to match, so on
+  such a day the retrospective query ran against the previous day's table and
+  could add a second copy of a sample that had just been inserted.
 
 - 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

docs/example_config.toml

@@ -82,9 +82,19 @@ num_threads=-1
 # limit.
 memory_limit=32
 
-# A list of sample IDs (strings) for unconditional inclusion (e.g., to 
-# help seed major saltation events).
-include_samples=[]
+# A list of "seed" groups for unconditional inclusion (e.g., to help seed
+# major saltation events). Each entry is a list of sample IDs that are
+# inserted into the ARG together as a single local tree: a tree is inferred
+# over the group's haplotypes, the haplotype of the group's inferred ancestor
+# is matched against the ARG without recombination, and the whole group is
+# attached at that single placement.
+#
+# A group is inserted on the *minimum* date over its members. Members with
+# later dates keep their real dates, and so are given negative ("in the
+# future") node times until their real date is reached.
+#
+#   seed_groups=[["BA.1_strain_a", "BA.1_strain_b"], ["BA.2_strain_c"]]
+seed_groups=[]
 
 # Override specific parameter values over a time period.
 [[override]]

sc2ts/dataset.py

@@ -646,9 +646,11 @@ def tmp_dataset(path, alignments, date="2020-01-01", contig_id=None):
     # Minimal hacky thing for testing. Should refactor into something more useful.
     if contig_id is None:
         contig_id = core.REFERENCE_STRAIN
+    if isinstance(date, str):
+        date = [date] * len(alignments)
     sequence_length = len(next(iter(alignments.values()))) + 1
     Dataset.new(path, sequence_length=sequence_length, contig_id=contig_id)
     Dataset.append_alignments(path, alignments)
-    df = pd.DataFrame({"strain": alignments.keys(), "date": [date] * len(alignments)})
+    df = pd.DataFrame({"strain": alignments.keys(), "date": date})
     Dataset.add_metadata(path, df.set_index("strain"))
     return Dataset(path)