fix subsetting variants and samples at the same time for bcf

Zilong-Li · Nov 27, 2024 · b703214 · b703214
1 parent 9522938
commit b703214
Show file tree

Hide file tree

Showing 8 changed files with 32 additions and 12 deletions.
diff --git a/news.org b/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=

diff --git a/test/Makefile b/test/Makefile
@@ -27,6 +27,6 @@ clean:
 	rm -f *.o *.bin
 
 test:
-	./bcf-reader.bin --success
-	./bcf-writer.bin --success
-	./bcf-variant.bin --success
+	./bcf-reader.bin --success 2> /dev/null
+	./bcf-writer.bin --success 2> /dev/null
+	./bcf-variant.bin --success 2> /dev/null
diff --git a/test/bcf-reader.cpp b/test/bcf-reader.cpp
@@ -158,7 +158,7 @@ TEST_CASE("throw error if the region is not valid", "[bcf-reader]")
     BcfReader * br = nullptr;
     try
     {
-        br = new BcfReader("test.vcf.gz", "XXXX");
+        br = new BcfReader("test-region.vcf.gz", "XXXX");
         delete br;
     }
     catch(exception & e)
@@ -170,25 +170,41 @@ TEST_CASE("throw error if the region is not valid", "[bcf-reader]")
 TEST_CASE("can check the status of a region", "[bcf-reader]")
 {
     int status;
-    BcfReader br("test.vcf.gz");
+    BcfReader br("test-region.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]")
 {
-    // BcfReader br("test.vcf.gz", "chr21:5030089-5030090");
     int nVariants = -1;
-    BcfReader br("test.vcf.gz");
+    BcfReader br("test-region.vcf.gz");
     nVariants = br.getVariantsCount("chr21:5030089-5030090");
     REQUIRE(nVariants == 0);
     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++;
+    }
+}
diff --git a/test/test.vcf.gz → test/test-region.vcf.gz b/test/test.vcf.gz → test/test-region.vcf.gz
diff --git a/test/test.vcf.gz.tbi → test/test-region.vcf.gz.tbi b/test/test.vcf.gz.tbi → test/test-region.vcf.gz.tbi
diff --git a/test/test.bcf b/test/test.bcf
diff --git a/test/test.bcf.csi b/test/test.bcf.csi
diff --git a/vcfpp.h b/vcfpp.h
@@ -1555,7 +1555,7 @@ class BcfReader
     std::shared_ptr<hts_itr_t> itr; // hts iterator
     kstring_t s = {0, 0, NULL}; // kstring
     std::string fname;
-    bool isBcf; // if the input file is bcf or vcf;
+    bool isBcf = false; // if the input file is bcf or vcf;
 
   public:
     /// a BcfHeader object
@@ -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);
             }