diff --git a/news.org b/news.org index 50ebf90..91aa4e3 100644 --- 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 index a1e0108..0216b80 100644 --- 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 index 11fe014..2dc4e4f 100644 --- 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 ds; + int n = 1; + while(br.getNextVariant(var)) + { + var.getFORMAT("DS", ds); + if(n == 1) REQUIRE(ds == vector{2, 2}); + if(n == 2) REQUIRE(ds == vector{0.195, 0.000}); + if(n == 3) REQUIRE(ds == vector{1.000, 0.000}); + if(n == 4) REQUIRE(ds == vector{2.000, 0.036}); + n++; + } +} diff --git a/test/test.vcf.gz b/test/test-region.vcf.gz similarity index 100% rename from test/test.vcf.gz rename to test/test-region.vcf.gz diff --git a/test/test.vcf.gz.tbi b/test/test-region.vcf.gz.tbi similarity index 100% rename from test/test.vcf.gz.tbi rename to test/test-region.vcf.gz.tbi diff --git a/test/test.bcf b/test/test.bcf new file mode 100644 index 0000000..3a30747 Binary files /dev/null and b/test/test.bcf differ diff --git a/test/test.bcf.csi b/test/test.bcf.csi new file mode 100644 index 0000000..d115bd0 Binary files /dev/null and b/test/test.bcf.csi differ diff --git a/vcfpp.h b/vcfpp.h index 0755cbe..354d657 100644 --- a/vcfpp.h +++ b/vcfpp.h @@ -1555,7 +1555,7 @@ class BcfReader std::shared_ptr 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); }