diff --git a/CHANGELOG.md b/CHANGELOG.md index 316d446..9ea3f9b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/sc2ts/inference.py b/sc2ts/inference.py index 9b8a522..a543388 100644 --- a/sc2ts/inference.py +++ b/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 diff --git a/tests/test_inference.py b/tests/test_inference.py index 4bd70cf..5d838bb 100644 --- a/tests/test_inference.py +++ b/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,