diff --git a/inst/extdata/test.bcf b/inst/extdata/test.bcf
new file mode 100644
index 0000000..3a30747
Binary files /dev/null and b/inst/extdata/test.bcf differ
diff --git a/inst/extdata/test.bcf.csi b/inst/extdata/test.bcf.csi
new file mode 100644
index 0000000..d115bd0
Binary files /dev/null and b/inst/extdata/test.bcf.csi differ
diff --git a/tests/testthat/test-vcf-reader.R b/tests/testthat/test-vcf-reader.R
index 5e4c0d4..5ed05a9 100644
--- a/tests/testthat/test-vcf-reader.R
+++ b/tests/testthat/test-vcf-reader.R
@@ -1,5 +1,6 @@
 library(testthat)
 
+bcffile <- system.file("extdata", "test.bcf", package="vcfppR")
 vcffile <- system.file("extdata", "raw.gt.vcf.gz", package="vcfppR")
 svfile <- system.file("extdata", "sv.vcf.gz", package="vcfppR")
 
@@ -268,4 +269,14 @@ test_that("can set region back and forth", {
   
 })
 
-
+test_that("bcfreader: subset samples and region for BCF", {
+  br <- vcfreader$new(bcffile, "1:10000-13000", "I1,I10")
+  expect_true(br$variant())
+  expect_identical(br$genotypes(F), c(1L, 1L, 0L, 1L))
+  expect_true(br$variant())
+  expect_identical(br$genotypes(F), c(0L, 0L, 0L, 0L))
+  expect_true(br$variant())
+  expect_identical(br$genotypes(F), c(1L, 0L, 0L, 0L))
+  expect_true(br$variant())
+  expect_identical(br$genotypes(F), c(1L, 1L, 0L, 1L))
+})
diff --git a/tests/testthat/test-vcf-table.R b/tests/testthat/test-vcf-table.R
index dc71f86..aacfbb8 100644
--- a/tests/testthat/test-vcf-table.R
+++ b/tests/testthat/test-vcf-table.R
@@ -1,5 +1,6 @@
 library(testthat)
 
+bcffile <- system.file("extdata", "test.bcf", package="vcfppR")
 vcffile <- system.file("extdata", "raw.gt.vcf.gz", package="vcfppR")
 svfile <- system.file("extdata", "sv.vcf.gz", package="vcfppR")
 
@@ -76,3 +77,9 @@ test_that("extract EV (String, number=1) from sv file", {
 test_that("query a non-existing region", {
   expect_error(res <- vcftable(vcffile, region = "chr14:0-1000"))
 })
+
+test_that("subset samples and region for BCF", {
+  res <- vcftable(bcffile, region =  "1:10000-13000", samples =  "I1,I30", vartype = "snps", format = "DS")
+  m <- matrix(c(2.000, 2.000,  0.195, 0.000, 1.000, 0.000, 2.000, 0.036), byrow = TRUE, ncol = 2)
+  expect_equal(res[[10]], m, tolerance = 1e6)
+})