fix subsetting variants and samples at the same time for bcf

news.org

@@ -1,4 +1,7 @@
 #+title: News and Changes
+* v0.6.1
+- fix subsetting for bcf [[https://github.com/Zilong-Li/vcfppR/issues/11][#11]]
+
 * v0.6.0
 - add =BcfReader::getStatus=
 - fix =BcfReader::getVariantsCount=

test/bcf-reader.cpp

@@ -172,14 +172,14 @@ TEST_CASE("can check the status of a region", "[bcf-reader]")
     int status;
     BcfReader br("test.vcf.gz");
     status = br.getStatus("chr21:5030089-5030090");
-    REQUIRE(status == 0);  // valid but empty
+    REQUIRE(status == 0); // valid but empty
     status = br.getStatus("chr21:5030089-");
-    REQUIRE(status == 1);  // valid and not empty
+    REQUIRE(status == 1); // valid and not empty
     status = br.getStatus("chr22:5030089");
-    REQUIRE(status == -2);  // not valid or not found in the VCF
+    REQUIRE(status == -2); // not valid or not found in the VCF
     BcfReader br2("test-no-index.vcf.gz");
     status = br2.getStatus("chr21");
-    REQUIRE(status == -1);  // no index file found
+    REQUIRE(status == -1); // no index file found
 }
 
 TEST_CASE("can count the number of variants in a valid region", "[bcf-reader]")
@@ -192,3 +192,20 @@ TEST_CASE("can count the number of variants in a valid region", "[bcf-reader]")
     nVariants = br.getVariantsCount("chr21:5030089-");
     REQUIRE(nVariants == 13);
 }
+
+TEST_CASE("subset both region and samples for bcf", "[bcf-reader]")
+{
+    BcfReader br("test.bcf", "1:10000-13000", "I1,I30");
+    BcfRecord var(br.header);
+    vector<float> ds;
+    int n = 1;
+    while(br.getNextVariant(var))
+    {
+        var.getFORMAT("DS", ds);
+        if(n == 1) REQUIRE(ds == vector<float>{2, 2});
+        if(n == 2) REQUIRE(ds == vector<float>{0.195, 0.000});
+        if(n == 3) REQUIRE(ds == vector<float>{1.000, 0.000});
+        if(n == 4) REQUIRE(ds == vector<float>{2.000, 0.036});
+        n++;
+    }
+}

vcfpp.h

@@ -1748,6 +1748,7 @@ class BcfReader
             if(isBcf)
             {
                 ret = bcf_itr_next(fp.get(), itr.get(), r.line.get());
+                bcf_subset_format(r.header->hdr, r.line.get()); // has to be called explicitly for bcf
                 bcf_unpack(r.line.get(), BCF_UN_ALL);
                 return (ret >= 0);
             }