Initial first tests on real data

Author: jeromekelleher

sc2ts/inference.py

@@ -563,10 +563,17 @@ def extend(
     show_progress=False,
     max_submission_delay=None,
     max_daily_samples=None,
-    num_threads=None,
+    num_threads=0,
     precision=None,
     rng=None,
 ):
+    if num_mismatches is None:
+        num_mismatches = 3
+    if max_hmm_cost is None:
+        max_hmm_cost = 5
+    if min_group_size is None:
+        min_group_size = 10
+
     check_base_ts(base_ts)
     logger.info(
         f"Extend {date}; ts:nodes={base_ts.num_nodes};samples={base_ts.num_samples};"

tests/test_inference.py

@@ -485,106 +485,53 @@ def test_high_recomb_mutation(self):
         self.check_double_mirror(ts)
 
 
-# # @pytest.fixture
-# def small_sd_fixture():
-#     reference = core.get_reference_sequence()
-#     print(reference)
-#     fasta = {"REF": reference}
-#     rows = [{"strain": "REF"}]
-#     sd = convert.convert_alignments(rows, fasta)
-
-#     return sd
-
-# ref = core.get_reference_sequence()
-# with tsinfer.SampleData(sequence_length=len(ref)) as sd:
-#     sd.add_individual(
-#         metadata={
-#             "strain": "A",
-#             "date": "2019-12-30",
-#             "date_submitted": "2020-01-02",
-#         }
-#     )
-#     sd.add_individual(
-#         metadata={
-#             "strain": "B",
-#             "date": "2020-01-01",
-#             "date_submitted": "2020-02-02",
-#         }
-#     )
-#     sd.add_individual(
-#         metadata={
-#             "strain": "C",
-#             "date": "2020-01-01",
-#             "date_submitted": "2020-02-02",
-#         }
-#     )
-#     sd.add_individual(
-#         metadata={
-#             "strain": "D",
-#             "date": "2020-01-02",
-#             "date_submitted": "2022-02-02",
-#         }
-#     )
-#     sd.add_individual(
-#         metadata={
-#             "strain": "E",
-#             "date": "2020-01-06",
-#             "date_submitted": "2020-02-02",
-#         }
-#     )
-# for
-# return sd
-
-# class TestInitialTables:
-#     def test_site_schema(self):
-#         sd = small_sd_fixture()
-#         pass
-
-
-@pytest.mark.skip()
-class TestInference:
-    def test_small_sd_times(self, small_sd_fixture):
-        ts = sc2ts.infer(small_sd_fixture)
-        inference.validate(small_sd_fixture, ts)
-        # Day 0 is Jan 6, and ultimate ancestor is one day older than the
-        # real root (reference)
-        np.testing.assert_array_equal(ts.nodes_time, [9, 8, 7, 5, 5, 4, 0])
-
-    def test_small_sd_submission_delay(self, small_sd_fixture):
-        ts = sc2ts.infer(small_sd_fixture, max_submission_delay=100)
-        strains = [ts.node(u).metadata["strain"] for u in ts.samples()]
-        # Strain D should be filtered out.
-        assert strains == ["A", "B", "C", "E"]
-        with pytest.raises(ValueError):
-            inference.validate(small_sd_fixture, ts)
-        inference.validate(small_sd_fixture, ts, max_submission_delay=100)
-
-    def test_daily_prefix(self, tmp_path, sd_fixture):
-        prefix = str(tmp_path) + "/x"
-        ts = sc2ts.infer(sd_fixture, daily_prefix=prefix)
-        paths = sorted(list(tmp_path.glob("x*")))
-        dailies = [tskit.load(x) for x in paths]
-        assert len(dailies) > 0
-        ts.tables.assert_equals(dailies[-1].tables)
-
-    @pytest.mark.parametrize("num_mismatches", [1, 2, 4, 1000])
-    def test_integrity(self, sd_fixture, num_mismatches):
-        ts = sc2ts.infer(sd_fixture, num_mismatches=num_mismatches)
-        assert ts.sequence_length == 29904
-        inference.validate(sd_fixture, ts)
-
-
-@pytest.mark.skip()
-class TestSubsetInferenceDefaults:
-    def test_metadata(self, ts_fixture):
-        for node in ts_fixture.nodes():
-            if node.flags == 0:
-                assert node.metadata == {}
-            elif node.flags == tskit.NODE_IS_SAMPLE:
-                assert "strain" in node.metadata
-            else:
-                assert node.flags == tsinfer.NODE_IS_IDENTICAL_SAMPLE_ANCESTOR
-                assert node.metadata == {}
-
-    def test_integrity(self, sd_fixture, ts_fixture):
-        inference.validate(sd_fixture, ts_fixture)
+class TestRealData:
+    def test_first_day(self, tmp_path, alignments_store, metadata_db):
+        ts = sc2ts.extend(
+            alignment_store=alignments_store,
+            metadata_db=metadata_db,
+            base_ts=sc2ts.initial_ts(),
+            date="2020-01-19",
+            match_db=sc2ts.MatchDb.initialise(tmp_path / "match.db"),
+        )
+        # 25.00┊ 0 ┊
+        #      ┊ ┃ ┊
+        # 24.00┊ 1 ┊
+        #      ┊ ┃ ┊
+        # 0.00 ┊ 2 ┊
+        #     0 29904
+        assert ts.num_trees == 1
+        assert ts.num_nodes == 3
+        assert ts.num_samples == 2
+        assert ts.num_mutations == 3
+        assert list(ts.nodes_time) == [25, 24, 0]
+        assert ts.metadata["sc2ts"]["date"] == "2020-01-19"
+        assert ts.metadata["sc2ts"]["samples_strain"] == [
+            "Wuhan/Hu-1/2019",
+            "SRR11772659",
+        ]
+        assert list(ts.samples()) == [1, 2]
+        assert ts.node(1).metadata["strain"] == "Wuhan/Hu-1/2019"
+        assert ts.node(2).metadata["strain"] == "SRR11772659"
+        assert list(ts.mutations_node) == [2, 2, 2]
+        assert list(ts.mutations_time) == [0, 0, 0]
+        assert list(ts.mutations_site) == [8632, 17816, 27786]
+        sc2ts_md = ts.node(2).metadata["sc2ts"]
+        assert len(sc2ts_md["mutations"]) == 3
+        for mut_md, mut in zip(sc2ts_md["mutations"], ts.mutations()):
+            assert mut_md["derived_state"] == mut.derived_state
+            assert mut_md["site_id"] == mut.site
+            assert mut_md["site_position"] == ts.sites_position[mut.site]
+            assert mut_md["inherited_state"] == ts.site(mut.site).ancestral_state
+        assert sc2ts_md["path"] == [{"left": 0, "parent": 1, "right": 29904}]
+        assert sc2ts_md["qc"] == {
+            "num_masked_sites": 133,
+            "original_base_composition": {
+                "A": 8894,
+                "C": 5472,
+                "G": 5850,
+                "N": 121,
+                "T": 9566,
+            },
+            "original_md5": "e96feaa72c4f4baba73c2e147ede7502",
+        }