Identify bug in A1/A2 allele handling

Author: jeromekelleher

bio2zarr/plink.py

@@ -12,9 +12,9 @@
 
 @dataclasses.dataclass
 class PlinkPaths:
-    bed_path: pathlib.Path
-    bim_path: pathlib.Path
-    fam_path: pathlib.Path
+    bed_path: str
+    bim_path: str
+    fam_path: str
 
 
 class PlinkFormat(vcz.Source):
@@ -25,13 +25,12 @@ def __init__(self, prefix):
         # TODO we will need support multiple chromosomes here to join
         # plinks into on big zarr. So, these will require multiple
         # bed and bim files, but should share a .fam
-        self.prefix = pathlib.Path(prefix)
+        self.prefix = str(prefix)
         paths = PlinkPaths(
-            self.prefix.with_suffix(".bed"),
-            self.prefix.with_suffix(".bim"),
-            self.prefix.with_suffix(".fam"),
+            self.prefix + ".bed",
+            self.prefix + ".bim",
+            self.prefix + ".fam",
         )
-
         self.bed = bed_reader.open_bed(
             paths.bed_path,
             bim_location=paths.bim_path,
@@ -73,8 +72,8 @@ def iter_id(self, start, stop):
         yield from self.bed.sid[start:stop]
 
     def iter_alleles_and_genotypes(self, start, stop, shape, num_alleles):
-        ref_field = self.bed.allele_1
-        alt_field = self.bed.allele_2
+        alt_field = self.bed.allele_1
+        ref_field = self.bed.allele_2
         bed_chunk = self.bed.read(slice(start, stop), dtype=np.int8).T
         gt = np.zeros(shape, dtype=np.int8)
         phased = np.zeros(shape[:-1], dtype=bool)
@@ -224,7 +223,7 @@ def validate(bed_path, zarr_path):
     root = zarr.open(store=zarr_path, mode="r")
     call_genotype = root["call_genotype"][:]
 
-    bed = bed_reader.open_bed(bed_path, count_A1=False, num_threads=1)
+    bed = bed_reader.open_bed(bed_path + ".bed", count_A1=False, num_threads=1)
 
     assert call_genotype.shape[0] == bed.sid_count
     assert call_genotype.shape[1] == bed.iid_count

tests/test_plink.py

@@ -31,8 +31,17 @@ def bed_path(self, tmp_path_factory):
             ],
             dtype=np.int8,
         )
-        bed_reader.to_bed(path, dosages.T)
-        return path
+        m = 7
+        d = {
+            "chromosome": ["chr1"] * m,
+            "sid": [f"id{j}" for j in range(m)],
+            "bp_position": range(1, m + 1),
+            "allele_1": ["A"] * m,
+            "allele_2": ["T"] * m,
+            "iid": [f"s{j}" for j in range(3)],
+        }
+        bed_reader.to_bed(path, dosages.T, properties=d)
+        return path.with_suffix("")
 
     @pytest.fixture(scope="class")
     def ds(self, tmp_path_factory, bed_path):
@@ -57,6 +66,22 @@ def test_genotypes(self, ds):
             ],
         )
 
+    def test_variant_position(self, ds):
+        nt.assert_array_equal(ds["variant_position"], np.arange(1, 8))
+
+    def test_variant_contig(self, ds):
+        nt.assert_array_equal(ds["variant_contig"], np.zeros(7))
+        nt.assert_array_equal(ds["contig_id"], ["chr1"])
+
+    def test_variant_id(self, ds):
+        nt.assert_array_equal(ds["variant_id"], [f"id{j}" for j in range(7)])
+
+    def test_variant_allele(self, ds):
+        nt.assert_array_equal(ds["variant_allele"], [["T", "A"] for _ in range(7)])
+
+    def test_sample_id(self, ds):
+        nt.assert_array_equal(ds["sample_id"], ["s0", "s1", "s2"])
+
     def test_missing_dependency(self):
         with mock.patch(
             "importlib.import_module",
@@ -75,25 +100,22 @@ def test_missing_dependency(self):
 
 class TestEqualSgkit:
     def test_simulated_example(self, tmp_path):
-        data_path = "tests/data/plink/"
-        bed_path = data_path + "plink_sim_10s_100v_10pmiss.bed"
-        fam_path = data_path + "plink_sim_10s_100v_10pmiss.fam"
-        bim_path = data_path + "plink_sim_10s_100v_10pmiss.bim"
+        prefix = "tests/data/plink/plink_sim_10s_100v_10pmiss"
         sg_ds = sgkit.io.plink.read_plink(
-            bed_path=bed_path, fam_path=fam_path, bim_path=bim_path
+            bed_path=prefix + ".bed", fam_path=prefix + ".fam", bim_path=prefix + ".bim"
         )
         out = tmp_path / "example.plink.zarr"
-        plink.convert(prefix=data_path + "/plink_sim_10s_100v_10pmiss", out=out)
+        plink.convert(prefix=prefix, out=out)
         ds = sg.load_dataset(out)
         nt.assert_array_equal(ds.call_genotype.values, sg_ds.call_genotype.values)
         nt.assert_array_equal(
             ds.call_genotype_mask.values, sg_ds.call_genotype_mask.values
         )
         # sgkit doesn't have phased
         nt.assert_array_equal(ds.variant_position.values, sg_ds.variant_position.values)
-        nt.assert_array_equal(
-            ds.variant_allele.values, sg_ds.variant_allele.values.astype("U")
-        )
+        # nt.assert_array_equal(
+        #     ds.variant_allele.values, sg_ds.variant_allele.values.astype("U")
+        # )
         nt.assert_array_equal(ds.variant_contig.values, sg_ds.variant_contig.values)
         nt.assert_array_equal(ds.variant_id.values, sg_ds.variant_id.values)
         # print(sg_ds.variant_id.values)
@@ -137,7 +159,7 @@ class TestExample:
 
     @pytest.fixture(scope="class")
     def ds(self, tmp_path_factory):
-        path = "tests/data/plink/example.bed"
+        path = "tests/data/plink/example"
         out = tmp_path_factory.mktemp("data") / "example.plink.zarr"
         plink.convert(path, out)
         return sg.load_dataset(out)
@@ -211,7 +233,7 @@ def test_sgkit(self, ds):
 class TestSimulatedExample:
     @pytest.fixture(scope="class")
     def ds(self, tmp_path_factory):
-        path = "tests/data/plink/plink_sim_10s_100v_10pmiss.bed"
+        path = "tests/data/plink/plink_sim_10s_100v_10pmiss"
         out = tmp_path_factory.mktemp("data") / "example.plink.zarr"
         plink.convert(path, out)
         return sg.load_dataset(out)
@@ -287,7 +309,7 @@ def test_contig_arrays(self, ds):
     def test_chunk_size(
         self, ds, tmp_path, variants_chunk_size, samples_chunk_size, worker_processes
     ):
-        path = "tests/data/plink/plink_sim_10s_100v_10pmiss.bed"
+        path = "tests/data/plink/plink_sim_10s_100v_10pmiss"
         out = tmp_path / "example.zarr"
         plink.convert(
             path,
@@ -318,7 +340,7 @@ def test_chunk_size(
 def test_by_validating(
     tmp_path, variants_chunk_size, samples_chunk_size, worker_processes
 ):
-    path = "tests/data/plink/plink_sim_10s_100v_10pmiss.bed"
+    path = "tests/data/plink/plink_sim_10s_100v_10pmiss"
     out = tmp_path / "example.zarr"
     plink.convert(
         path,
@@ -372,7 +394,7 @@ def multi_contig_bed_path(self, tmp_path_factory):
                 # Format: fam_id, ind_id, pat_id, mat_id, sex, phenotype
                 f.write(f"fam{i} ind{i} 0 0 0 -9\n")
 
-        return path
+        return path.with_suffix("")
 
     @pytest.fixture(scope="class")
     def ds(self, tmp_path_factory, multi_contig_bed_path):