Preserve VCF fields defined in header that don't appear otherwise

Author: tomwhite

bio2zarr/icf.py

@@ -79,6 +79,14 @@ def full_name(self):
             return self.name
         return f"{self.category}/{self.name}"
 
+    @property
+    def max_number(self):
+        if self.vcf_number in ("R", "A", "G", "."):
+            return self.summary.max_number
+        else:
+            # use declared number if larger than max found
+            return max(self.summary.max_number, int(self.vcf_number))
+
     def smallest_dtype(self):
         """
         Returns the smallest dtype suitable for this field based

bio2zarr/vcz.py

@@ -136,12 +136,14 @@ def from_field(
             dimensions.append("samples")
         if array_name is None:
             array_name = prefix + vcf_field.name
+
         # TODO make an option to add in the empty extra dimension
-        if (
-            vcf_field.summary.max_number > 0 and vcf_field.vcf_number in ("R", "A", "G")
-        ) or (vcf_field.summary.max_number > 1 or vcf_field.full_name == "FORMAT/LAA"):
-            shape.append(vcf_field.summary.max_number)
-            chunks.append(vcf_field.summary.max_number)
+        max_number = vcf_field.summary.max_number
+        if (max_number > 0 and vcf_field.vcf_number in ("R", "A", "G")) or (
+            max_number > 1 or vcf_field.full_name == "FORMAT/LAA"
+        ):
+            shape.append(max_number)
+            chunks.append(max_number)
             # TODO we should really be checking this to see if the named dimensions
             # are actually correct.
             if vcf_field.vcf_number == "R":

tests/data/vcf/chr22.vcf.gz

@@ -563,6 +563,21 @@ def test_gts(self, ds):
         assert "call_genotype" not in ds
 
 
+class TestChr22Example:
+    data_path = "tests/data/vcf/chr22.vcf.gz"
+
+    @pytest.fixture(scope="class")
+    def ds(self, tmp_path_factory):
+        out = tmp_path_factory.mktemp("data") / "example.vcf.zarr"
+        icf.convert([self.data_path], out, worker_processes=0)
+        return sg.load_dataset(out)
+
+    def test_call_SB(self, ds):
+        # fixes https://github.com/sgkit-dev/bio2zarr/issues/355
+        assert ds.call_SB.dims == ("variants", "samples", "FORMAT_SB_dim")
+        assert ds.call_SB.shape == (100, 100, 4)
+
+
 class Test1000G2020Example:
     data_path = "tests/data/vcf/1kg_2020_chrM.vcf.gz"