Allow empty ancestral alleles

tests/test_sgkit.py

@@ -615,9 +615,36 @@ def test_empty_alleles_not_at_end(self, tmp_path):
         )
         sgkit.save_dataset(ds, path)
         samples = tsinfer.SgkitSampleData(path)
-        with pytest.raises(ValueError, match="Empty alleles must be at the end"):
+        with pytest.raises(ValueError, match="empty alleles must be at the end"):
             tsinfer.infer(samples)
 
+    def test_empty_ancestral_alleles(self, tmp_path):
+        path = tmp_path / "data.zarr"
+        ds = sgkit.simulate_genotype_call_dataset(n_variant=3, n_sample=3, n_ploidy=1)
+        ds["variant_allele"] = (
+            ds["variant_allele"].dims,
+            np.array(
+                [["", "A", "C", ""], ["A", "C", "", ""], ["A", "C", "", ""]], dtype="S1"
+            ),
+        )
+        ds["variant_ancestral_allele"] = (
+            ["variants"],
+            np.array(["", "A", ""], dtype="S1"),
+        )
+        sgkit.save_dataset(ds, path)
+        samples = tsinfer.SgkitSampleData(path)
+        for v in samples.variants(recode_ancestral=True):
+            if v.site.id == 0:
+                assert v.site.ancestral_state in (b"", "")
+                assert len(v.alleles) == 3
+                assert v.alleles[0] in (b"", "")
+            elif v.site.id == 1:
+                assert v.site.ancestral_state in (b"A", "A")
+                assert len(v.alleles) == 2
+            elif v.site.id == 3:
+                assert v.site.ancestral_state is None
+                assert len(v.alleles) == 3
+
 
 class TestSgkitMatchSamplesToDisk:
     @pytest.mark.skipif(sys.platform == "win32", reason="No cyvcf2 on windows")

tsinfer/formats.py

@@ -2584,24 +2584,34 @@ def variants(self, sites=None, recode_ancestral=None):
             genos = genos.reshape(self.num_samples)
             aa = site.ancestral_allele
             alleles = site.alleles
-            if aa != MISSING_DATA and aa > 0 and recode_ancestral:
-                # Need to recode this site
-                alleles = site.reorder_alleles()
-                # re-map the genotypes
-                geno_map = np.arange(len(alleles) - MISSING_DATA, dtype=genos.dtype)
-                geno_map[MISSING_DATA] = MISSING_DATA
-                geno_map[aa] = 0
-                geno_map[0:aa] += 1
-                genos = geno_map[genos]
+            aa_empty = False
+            if aa != MISSING_DATA:
+                if alleles[aa] == b"" or alleles[aa] == "":
+                    aa_empty = True
+                if aa > 0 and recode_ancestral:
+                    # Need to recode this site
+                    alleles = site.reorder_alleles()
+                    # re-map the genotypes
+                    geno_map = np.arange(len(alleles) - MISSING_DATA, dtype=genos.dtype)
+                    geno_map[MISSING_DATA] = MISSING_DATA
+                    geno_map[aa] = 0
+                    geno_map[0:aa] += 1
+                    genos = geno_map[genos]
             # Filter out empty alleles, as sgkit pads with them so that all sites have
             # the same number of alleles. This is only safe if the empty
             # alleles are at the end of the list, so check this.
             non_empty_alleles = []
             empty_seen = False
-            for allele in alleles:
+            for i, allele in enumerate(alleles):
                 if allele != b"" and allele != "":
                     if empty_seen:
-                        raise ValueError("Empty alleles must be at the end")
+                        raise ValueError(
+                            f"Site {site.id} (pos {site.position}): empty alleles "
+                            f"must be at the end, but alleles are {alleles}"
+                        )
+                    non_empty_alleles.append(allele)
+                elif i == 0 and aa_empty:
+                    # Single empty allele allowed if it is the starting ancestral allele
                     non_empty_alleles.append(allele)
                 else:
                     empty_seen = True