inst/tinytest/test_bcf.R

In Rduckhts: 'DuckDB' High Throughput Sequencing File Formats Reader Extension

library(tinytest)
library(DBI)

test_bcf_string_format_lists <- function() {
  con <- rduckhts_connect()
  on.exit(
    {
      dbDisconnect(con, shutdown = TRUE)
    },
    add = TRUE
  )

  bcf_path <- system.file(
    "extdata",
    "format_string_list.vcf",
    package = "Rduckhts"
  )
  fixed_path <- system.file(
    "extdata",
    "fixed_count_arrays.vcf",
    package = "Rduckhts"
  )
  mapping_path <- system.file(
    "extdata",
    "mapping_number_families.vcf",
    package = "Rduckhts"
  )
  vep_path <- system.file("extdata", "test_vep.vcf", package = "Rduckhts")
  expect_true(file.exists(bcf_path))
  expect_true(file.exists(fixed_path))
  expect_true(file.exists(mapping_path))
  expect_true(file.exists(vep_path))

  expect_silent(rduckhts_bcf(
    con,
    "bcf_string_lists",
    bcf_path,
    decompression_threads = 0,
    overwrite = TRUE
  ))
  expect_silent(rduckhts_bcf(
    con,
    "bcf_fixed_counts",
    fixed_path,
    overwrite = TRUE
  ))
  expect_error(
    rduckhts_bcf(
      con,
      "bcf_bad_threads",
      bcf_path,
      decompression_threads = -1,
      overwrite = TRUE
    ),
    pattern = "decompression_threads must be a single whole number"
  )
  expect_error(
    rduckhts_bcf(
      con,
      "bcf_bad_csq_types",
      vep_path,
      additional_csq_column_types = "DISTANCE BogusType",
      overwrite = TRUE
    ),
    pattern = "additional_csq_column_types must use bcftools-style"
  )

  type_row <- DBI::dbGetQuery(
    con,
    "SELECT typeof(FORMAT_LAA_SAMPLE1) AS t FROM bcf_string_lists LIMIT 1"
  )
  expect_equal(type_row$t[1], "VARCHAR[]")

  first_row <- DBI::dbGetQuery(
    con,
    paste(
      "SELECT list_extract(FORMAT_LAA_SAMPLE1, 1) AS laa1,",
      "list_extract(FORMAT_LAA_SAMPLE1, 2) AS laa2",
      "FROM bcf_string_lists WHERE POS = 73"
    )
  )
  expect_equal(first_row$laa1[1], "1")
  expect_equal(first_row$laa2[1], "2")

  missing_row <- DBI::dbGetQuery(
    con,
    paste(
      "SELECT FORMAT_LAA_SAMPLE1 IS NULL AS is_null,",
      "FORMAT_LAA_SAMPLE1::VARCHAR AS laa",
      "FROM bcf_string_lists WHERE POS = 263"
    )
  )
  expect_true(isTRUE(missing_row$is_null[1]))
  expect_true(is.na(missing_row$laa[1]))

  fixed_type <- DBI::dbGetQuery(
    con,
    "SELECT typeof(INFO_SB) AS info_t, typeof(FORMAT_HQ_S1) AS fmt_t FROM bcf_fixed_counts LIMIT 1"
  )
  expect_equal(fixed_type$info_t[1], "INTEGER[]")
  expect_equal(fixed_type$fmt_t[1], "INTEGER[]")

  fixed_row <- DBI::dbGetQuery(
    con,
    paste(
      "SELECT INFO_SB::VARCHAR AS info_sb,",
      "FORMAT_HQ_S1::VARCHAR AS format_hq",
      "FROM bcf_fixed_counts WHERE POS = 100"
    )
  )
  expect_equal(fixed_row$info_sb[1], "[10, 20, 30, 40]")
  expect_equal(fixed_row$format_hq[1], "[5, 9]")

  quoted_mapping <- DBI::dbQuoteString(con, mapping_path)

  number_dot_row <- DBI::dbGetQuery(
    con,
    sprintf(
      paste(
        "SELECT typeof(INFO_VSTR) AS info_vstr_type, INFO_VSTR::VARCHAR AS info_vstr,",
        "typeof(INFO_VINT) AS info_vint_type, INFO_VINT::VARCHAR AS info_vint,",
        "typeof(FORMAT_LAA_S1) AS fmt_laa_type, FORMAT_LAA_S1::VARCHAR AS fmt_laa",
        "FROM read_bcf(%s) WHERE POS = 200"
      ),
      quoted_mapping
    )
  )
  expect_equal(number_dot_row$info_vstr_type[1], "VARCHAR[]")
  expect_equal(number_dot_row$info_vstr[1], "[red, blue, green]")
  expect_equal(number_dot_row$info_vint_type[1], "INTEGER[]")
  expect_equal(number_dot_row$info_vint[1], "[1, 2, 3]")
  expect_equal(number_dot_row$fmt_laa_type[1], "VARCHAR[]")
  expect_equal(number_dot_row$fmt_laa[1], "[1, 2]")

}

test_bcf_projection_sql <- function() {
  con <- rduckhts_connect()
  on.exit(
    {
      dbDisconnect(con, shutdown = TRUE)
    },
    add = TRUE
  )
  bcf_path <- system.file("extdata", "formatcols.vcf.gz", package = "Rduckhts")
  expect_true(file.exists(bcf_path))
  quoted_path <- DBI::dbQuoteString(con, bcf_path)

  out <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "SELECT count(*) AS n, min(POS) AS min_pos, max(POS) AS max_pos, ",
        "min(FORMAT_S) AS fmt_s ",
        "FROM read_bcf(%s, tidy_format := true)"
      ),
      quoted_path
    )
  )
  expect_equal(out$n[1], 3)
  expect_equal(out$min_pos[1], 100)
  expect_equal(out$max_pos[1], 100)
  expect_equal(out$fmt_s[1], "a")

  expect_equal(dbGetQuery(con, paste(
    "SELECT count(*) AS n FROM duckdb_functions()",
    "WHERE function_name = 'read_bcf_v2'"
  ))$n[[1]], 0)

  vep_tidy_path <- system.file(
    "extdata",
    "test_vep_tidy.vcf",
    package = "Rduckhts"
  )
  expect_true(file.exists(vep_tidy_path))
  quoted_vep_tidy <- DBI::dbQuoteString(con, vep_tidy_path)
  vep_tidy <- DBI::dbGetQuery(
    con,
    sprintf(
      paste(
        "SELECT SAMPLE_ID, VEP_SYMBOL::VARCHAR AS symbol,",
        "VEP_DISTANCE::VARCHAR AS distance, FORMAT_GT",
        "FROM read_bcf(%s, tidy_format := true)",
        "ORDER BY SAMPLE_ID"
      ),
      quoted_vep_tidy
    )
  )
  expect_equal(vep_tidy$SAMPLE_ID, c("S1", "S2"))
  expect_equal(vep_tidy$symbol, c("[GENE1]", "[GENE1]"))
  expect_equal(vep_tidy$distance, c("[12]", "[12]"))
  expect_equal(vep_tidy$FORMAT_GT, c("0/1", "1/1"))

  indexed_bcf_path <- system.file("extdata", "vcf_file.bcf", package = "Rduckhts")
  quoted_indexed_bcf <- DBI::dbQuoteString(con, indexed_bcf_path)

  gt_edge_cases <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "SELECT string_agg(FORMAT_GT, ',' ORDER BY POS, SAMPLE_ID) AS gt ",
        "FROM read_bcf(%s, tidy_format := true) ",
        "WHERE FORMAT_GT IN ('2', '0/300', '240/260')"
      ),
      quoted_indexed_bcf
    )
  )
  expect_equal(gt_edge_cases$gt[1], "2,0/300,240/260")

  ploidy_edge_path <- system.file(
    "extdata",
    "genotype_ploidy_edge_cases.vcf",
    package = "Rduckhts"
  )
  expect_true(file.exists(ploidy_edge_path))
  ploidy_edges <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "SELECT string_agg(SAMPLE_ID || '=' || FORMAT_GT, ',' ORDER BY SAMPLE_ID) AS gt ",
        "FROM read_bcf(%s, tidy_format := true)"
      ),
      DBI::dbQuoteString(con, ploidy_edge_path)
    )
  )
  expect_equal(
    ploidy_edges$gt[1],
    "BIG=0/10,DIP=0|1,HAP=1,TET=0/1/2/3,TRI=0/1/2"
  )

  phased_fields_path <- system.file(
    "extdata",
    "phased_genotype_fields.vcf",
    package = "Rduckhts"
  )
  expect_true(file.exists(phased_fields_path))
  quoted_phased <- DBI::dbQuoteString(con, phased_fields_path)
  phased_gt <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "SELECT string_agg(SAMPLE_ID || '=' || FORMAT_GT || ':PS=' || coalesce(CAST(FORMAT_PS AS VARCHAR), 'NA'), ',' ORDER BY POS, SAMPLE_ID) AS gt ",
        "FROM read_bcf(%s, tidy_format := true)"
      ),
      quoted_phased
    )
  )
  expect_equal(
    phased_gt$gt[1],
    paste(
      "S1=0|1:PS=10,S2=1|2:PS=10,S1=1|0:PS=20,S2=0/1:PS=NA",
      "S1=1:PS=30,S2=0|1|1:PS=30,S1=0|1|1|2:PS=40,S2=2|2:PS=40",
      sep = ","
    )
  )
  phased_ploidy_lengths <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "SELECT string_agg(SAMPLE_ID || '=' || FORMAT_GT || ':PL=' || length(FORMAT_PL) || ':GP=' || length(FORMAT_GP) || ':DS=' || length(FORMAT_DS), ',' ORDER BY POS, SAMPLE_ID) AS lengths ",
        "FROM read_bcf(%s, tidy_format := true) ",
        "WHERE POS >= 30"
      ),
      quoted_phased
    )
  )
  expect_equal(
    phased_ploidy_lengths$lengths[1],
    "S1=1:PL=2:GP=2:DS=1,S2=0|1|1:PL=4:GP=4:DS=1,S1=0|1|1|2:PL=15:GP=15:DS=2,S2=2|2:PL=6:GP=6:DS=2"
  )

  tidy_chunk_path <- system.file(
    "extdata",
    "tidy_chunk_boundary.vcf",
    package = "Rduckhts"
  )
  expect_true(file.exists(tidy_chunk_path))
  tidy_chunk <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "WITH t AS (",
        "SELECT SAMPLE_ID, FORMAT_GT ",
        "FROM read_bcf(%s, tidy_format := true)",
        ") ",
        "SELECT ",
        "(SELECT count(*) FROM t) AS n, ",
        "(SELECT count(DISTINCT SAMPLE_ID) FROM t) AS samples, ",
        "(SELECT string_agg(SAMPLE_ID || '=' || FORMAT_GT, ',' ORDER BY SAMPLE_ID) ",
        " FROM t WHERE SAMPLE_ID IN ('S0001', 'S2048', 'S2049', 'S2053')) AS edge_gt"
      ),
      DBI::dbQuoteString(con, tidy_chunk_path)
    )
  )
  expect_equal(tidy_chunk$n[1], 2053)
  expect_equal(tidy_chunk$samples[1], 2053)
  expect_equal(tidy_chunk$edge_gt[1], "S0001=0/1,S2048=0/1,S2049=0/1,S2053=0/1")

  ensembl_path <- system.file(
    "extdata",
    "ensembl_release_consequences.vcf",
    package = "Rduckhts"
  )
  expect_true(file.exists(ensembl_path))
  ensembl_counts <- DBI::dbGetQuery(
    con,
    sprintf(
      paste0(
        "SELECT count(*) AS records, ",
        "sum(list_count(VEP_Consequence)) AS csq_entries, ",
        "sum(list_count(INFO_VE)) AS ve_entries ",
        "FROM read_bcf(%s)"
      ),
      DBI::dbQuoteString(con, ensembl_path)
    )
  )
  expect_equal(ensembl_counts$records[1], 2)
  expect_equal(ensembl_counts$csq_entries[1], 6)
  expect_equal(ensembl_counts$ve_entries[1], 9)

}

test_bcf_materialization <- function() {
  con <- rduckhts_connect()
  on.exit(dbDisconnect(con, shutdown = TRUE))
  path <- system.file("extdata", "vcf_file.bcf", package = "Rduckhts")
  quoted <- dbQuoteString(con, path)
  expect_equal(dbGetQuery(con, paste(
    "SELECT count(*) AS n FROM duckdb_functions()",
    "WHERE function_name IN ('read_bcf_v2','read_bcf_appender')"))$n[[1]], 0)
  dbExecute(con, sprintf(paste(
    "CREATE TABLE bcf_materialized AS SELECT * FROM read_bcf(%s,",
    "region := '1:3062915-3062915,1:3062918-3062918,1:3062915-3062915,2:3199812-3199812',",
    "tidy_format := true)"), quoted))
  actual <- dbGetQuery(con, "SELECT * FROM bcf_materialized ORDER BY CHROM,POS,SAMPLE_ID")
  expected <- dbGetQuery(con, sprintf(paste(
    "SELECT * FROM read_bcf(%s, scan_mode := 'sequential', tidy_format := true)",
    "WHERE (CHROM='1' AND POS IN (3062915,3062918)) OR (CHROM='2' AND POS=3199812)",
    "ORDER BY CHROM,POS,SAMPLE_ID"), quoted))
  expect_equal(nrow(actual), 6L)
  expect_equal(actual, expected)
  dbExecute(con, sprintf(paste(
    "CREATE OR REPLACE TABLE bcf_materialized AS SELECT * FROM read_bcf(%s,",
    "region := 'absent_contig:1-1000', tidy_format := true)"), quoted))
  expect_equal(dbGetQuery(con, "SELECT count(*) AS n FROM bcf_materialized")$n[[1]], 0)
}

test_bcf_record_cache <- function() {
  con <- rduckhts_connect()
  on.exit(dbDisconnect(con, shutdown = TRUE))
  path <- system.file("extdata", "bcf_cache_lifecycle.vcf", package = "Rduckhts")
  expect_true(file.exists(path))
  out <- dbGetQuery(con, sprintf(paste(
    "SELECT POS, SAMPLE_ID, VEP_SYMBOL::VARCHAR AS symbol, INFO_DP,",
    "FORMAT_GT, FORMAT_AD::VARCHAR AS ad, FORMAT_ST::VARCHAR AS st",
    "FROM read_bcf(%s, tidy_format := true) ORDER BY POS, SAMPLE_ID"
  ), dbQuoteString(con, path)))
  expect_equal(as.integer(out$POS), rep(c(10L, 20L, 30L), each = 2L))
  expect_equal(out$SAMPLE_ID, rep(c("S1", "S2"), 3L))
  expect_equal(out$symbol, rep(c("[GENE1]", NA, "[GENE2, GENE3]"), each = 2L))
  expect_equal(out$INFO_DP, rep(c(7L, NA_integer_, 1234L), each = 2L))
  expect_equal(out$FORMAT_GT, c("0/1", "1/1", "./.", "0", "0|1|2", "2/2"))
  expect_equal(out$ad, c("[6, 1]", "[0, 7]", NA, NA, "[4, 5, 6]", "[0, 0, 9]"))
  expect_equal(out$st, c("[x, y]", "[z]", NA, NA, "[longer, value, three]", "[last]"))

  path <- system.file("extdata", "tidy_chunk_boundary.vcf", package = "Rduckhts")
  out <- dbGetQuery(con, sprintf(paste(
    "SELECT count(*) AS n, count(*) FILTER (WHERE VEP_SYMBOL = ['GENE1']",
    "AND VEP_DISTANCE = [12] AND FORMAT_GT = '0/1') AS correct",
    "FROM read_bcf(%s, tidy_format := true)"
  ), dbQuoteString(con, path)))
  expect_equal(out$n[[1]], 2053)
  expect_equal(out$correct[[1]], 2053)
}

test_bcf_index_snapshot <- function() {
  sys.source(system.file("tinytest", "bcf_snapshot_query.R", package = "Rduckhts",
                         mustWork = TRUE), envir = environment())
  con <- rduckhts_connect()
  on.exit(dbDisconnect(con, shutdown = TRUE))
  directory <- tempfile("duckhts-bcf-snapshot-")
  dir.create(directory)
  on.exit(unlink(directory, recursive = TRUE), add = TRUE)
  bcf_scan_expected_tables(con)
  fixture <- function(suffix) system.file("extdata", paste0("bcf_scan_contigs.", suffix),
                                         package = "Rduckhts", mustWork = TRUE)
  quote <- function(value) as.character(dbQuoteString(con, value))
  # Both formats have two occupied contigs. Exercise auto discovery and explicit
  # paths, full/region/sequential scans, complete wide/tidy rows and index counts.
  for (kind in c("bcf", "tbi", "csi")) for (explicit in c(FALSE, TRUE)) {
    suffix <- if (kind == "bcf") "bcf" else "full.vcf.gz"
    index_suffix <- if (kind == "bcf") "csi" else kind
    path <- file.path(directory, paste0("input.", suffix))
    expect_true(file.copy(fixture(suffix), path, overwrite = TRUE))
    index <- if (explicit) file.path(directory, "private.index") else paste0(path, ".", index_suffix)
    hidden <- paste0(index, ".hidden")
    min_shift <- if (kind == "tbi") 0L else 14L
    build <- sprintf("SELECT * FROM bcf_index(%s,index_path:=%s,min_shift:=%d,threads:=1)",
                     quote(path), quote(index), min_shift)
    for (workers in c(1L, 4L)) for (tidy in c(FALSE, TRUE)) for (mode in c("auto", "region", "sequential")) {
      dbExecute(con, sprintf("SET threads=%d", workers))
      dbGetQuery(con, build)
      source <- sprintf("read_bcf(%s,tidy_format:=%s,decompression_threads:=0%s%s)",
        quote(path), tolower(tidy), if (explicit) paste0(",index_path:=", quote(index)) else "",
        switch(mode, auto = "", region = ",region:='chr3:20-40,chr3:20-20'",
               sequential = ",scan_mode:='sequential'"))
      # Compare every typed column in a deterministic row order.
      projection <- "SELECT POS, * EXCLUDE (POS) FROM "
      sql <- paste0(projection, source, " ORDER BY ALL")
      expected <- dbGetQuery(con, paste0(projection, if (tidy) "expected_tidy" else "expected_sites",
        if (mode == "region") " WHERE CHROM='chr3'" else "", " ORDER BY ALL"))
      expect_equal(dbGetQuery(con, sql), expected)
      # Rewrites in bind callbacks occur after read_bcf binds and before it scans.
      for (replacement in c("empty", "shifted", "identical")) {
        replacement_path <- if (replacement == "identical") path else
          fixture(paste0(replacement, ".", if (kind == "bcf") "bcf" else "vcf.gz"))
        for (count in c(FALSE, TRUE)) {
          bound_sql <- .bcf_snapshot_query(con, source, path, index,
                                           replacement_path, min_shift, count)
          observed <- dbGetQuery(con, bound_sql)
          if (count) expect_equal(observed$n[[1]], nrow(expected))
          else expect_equal(observed, expected)
        }
      }
      # Before reader binding, the shifted BCF offset hits an invalid BGZF header;
      # shifted Tabix offsets are beyond this local VCF's EOF (bcftools agrees).
      # Sequential mode ignores either index and must return the complete file.
      shifted_path <- fixture(paste0("shifted.", if (kind == "bcf") "bcf" else "vcf.gz"))
      control_sql <- .bcf_snapshot_query(con, source, path, index, shifted_path,
                                         min_shift, replace_before_bind = TRUE)
      if (mode == "sequential") {
        expect_equal(dbGetQuery(con, control_sql), expected)
      } else if (kind == "bcf") {
        expect_error(dbGetQuery(con, control_sql),
                     pattern = "read_bcf: failed to read or parse BCF/VCF record")
      } else {
        expect_equal(dbGetQuery(con, control_sql), expected[FALSE, ])
      }
      # Completed results also remain valid after later file changes. This does
      # not test scan overlap: dbSendQuery can finish the scan before returning.
      fetch_after_index_change <- function(query, damage) {
        dbGetQuery(con, build)
        if (damage == "identical") {
          index_bytes <- readBin(index, "raw", n = file.info(index)$size)
        }
        result <- dbSendQuery(con, query)
        on.exit({
          dbClearResult(result)
          if (file.exists(hidden)) unlink(hidden)
        })
        if (damage == "missing") {
          expect_true(file.rename(index, hidden))
        } else if (damage == "corrupt") {
          writeLines("broken index", index)
        } else if (damage == "identical") {
          writeBin(index_bytes, index)
        } else {
          replacement <- paste0(damage, ".", if (kind == "bcf") "bcf.csi" else paste0("vcf.gz.index.", kind))
          expect_true(file.copy(fixture(replacement), index, overwrite = TRUE))
        }
        dbFetch(result)
      }
      for (damage in c("missing", "corrupt", "empty", "shifted", "identical")) {
        expect_equal(fetch_after_index_change(sql, damage), expected)
        count_sql <- paste("SELECT count(*) AS n FROM", source)
        expect_equal(fetch_after_index_change(count_sql, damage)$n[[1]], nrow(expected))
        expect_equal(dbGetQuery(con, "SELECT 4242 AS n")$n, 4242L)
        if (damage %in% c("missing", "corrupt")) {
          if (mode == "region") expect_error(dbGetQuery(con, sql), pattern = "region query requires an index")
          else expect_equal(dbGetQuery(con, sql), expected) # New bind: streaming fallback.
        }
      }
    }
    # A replacement data+index pair supplied before bind defines a new plan.
    expect_true(file.copy(fixture(paste0("shifted.", if (kind == "bcf") "bcf" else "vcf.gz")),
                          path, overwrite = TRUE))
    dbGetQuery(con, build)
    expect_equal(dbGetQuery(con, sql), expected)
    auto_sql <- sub(",scan_mode:='sequential'", "", sql, fixed = TRUE)
    expect_equal(dbGetQuery(con, auto_sql), expected)
    region_sql <- sub(",scan_mode:='sequential'", ",region:='chr3:20-40,chr3:20-20'", sql, fixed = TRUE)
    expected_region <- expected[expected$CHROM == "chr3", ]
    rownames(expected_region) <- NULL
    expect_equal(dbGetQuery(con, region_sql), expected_region)
    unlink(c(path, index))
  }
}

bcf_scan_expected_tables <- function(con) {
  dbExecute(con, "CREATE TABLE expected_sites(CHROM VARCHAR, POS BIGINT, ID VARCHAR, REF VARCHAR,
    ALT VARCHAR[], QUAL DOUBLE, FILTER VARCHAR[], INFO_DP INTEGER,
    FORMAT_GT_S1 VARCHAR, FORMAT_GT_S2 VARCHAR)")
  dbExecute(con, "INSERT INTO expected_sites VALUES
    ('chr1',10,'known','A',['C'],60,['PASS'],7,'0/1','1/1'),
    ('chr3',20,'duplicate','G',['T'],50,['PASS'],8,'./1','0/0'),
    ('chr3',20,'duplicate','G',['T'],50,['PASS'],8,'./1','0/0'),
    ('chr3',20,'allele','G',['A'],30,['PASS'],9,'1/1','./.'),
    ('chr3',40,'last','T',['G'],NULL,['PASS'],NULL,'1','0|1')")
  dbExecute(con, "CREATE VIEW expected_tidy AS
    SELECT CHROM,POS,ID,REF,ALT,QUAL,FILTER,INFO_DP,SAMPLE_ID,FORMAT_GT
    FROM expected_sites, LATERAL (VALUES ('S1',FORMAT_GT_S1),('S2',FORMAT_GT_S2))
    AS samples(SAMPLE_ID,FORMAT_GT)")
}

test_bcf_contig_dictionary <- function() {
  con <- rduckhts_connect()
  on.exit(dbDisconnect(con, shutdown = TRUE))
  bcf_scan_expected_tables(con)
  for (threads in c(1L, 4L)) for (tidy in c(FALSE, TRUE)) {
    dbExecute(con, sprintf("SET threads=%d", threads))
    expected_table <- if (tidy) "expected_tidy" else "expected_sites"
    expected <- dbGetQuery(con, paste("SELECT * FROM", expected_table, "ORDER BY ALL"))
    expected_region <- dbGetQuery(con, paste("SELECT * FROM", expected_table,
                                            "WHERE CHROM='chr3' ORDER BY ALL"))
    projection <- if (tidy) "CHROM,POS,ID,SAMPLE_ID,FORMAT_GT" else "CHROM,POS,ID,FORMAT_GT_S1"
    expected_projection <- dbGetQuery(con, paste("SELECT", projection, "FROM", expected_table, "ORDER BY ALL"))
    for (header in c("full", "partial", "none")) {
      path <- system.file("extdata", paste0("bcf_scan_contigs.", header, ".vcf.gz"), package = "Rduckhts")
      expect_true(nzchar(path))
      expect_false(file.exists(paste0(path, ".tbi")) || file.exists(paste0(path, ".csi")))
      source <- sprintf("read_bcf(%s, tidy_format:=%s, decompression_threads:=0",
                        dbQuoteString(con, path), tolower(tidy))
      # No-index auto fallback must have the same complete rows.
      expect_equal(dbGetQuery(con, paste0("SELECT * FROM ", source, ") ORDER BY ALL")), expected)
      for (kind in c("tbi", "csi")) {
        index <- paste0(path, ".index.", kind)
        expect_true(file.exists(index))
        indexed <- paste0(source, ", index_path:=", dbQuoteString(con, index))
        for (mode in c("auto", "sequential")) {
          sql <- paste0(indexed, ", scan_mode:=", dbQuoteString(con, mode), ")")
          expect_equal(dbGetQuery(con, paste("SELECT * FROM", sql, "ORDER BY ALL")), expected)
          expect_equal(dbGetQuery(con, paste("SELECT", projection, "FROM", sql, "ORDER BY ALL")),
                       expected_projection)
        }
        region <- paste0(indexed, ", region:='chr3:20-40,chr3:20-20')")
        expect_equal(dbGetQuery(con, paste("SELECT * FROM", region, "ORDER BY ALL")), expected_region)
      }
    }
    bcf <- system.file("extdata", "bcf_scan_contigs.bcf", package = "Rduckhts")
    expect_equal(dbGetQuery(con, sprintf("SELECT * FROM read_bcf(%s,tidy_format:=%s) ORDER BY ALL",
                                         dbQuoteString(con, bcf), tolower(tidy))), expected)
    empty <- system.file("extdata", "bcf_scan_contigs.empty.vcf.gz", package = "Rduckhts")
    for (kind in c("tbi", "csi")) {
      result <- dbGetQuery(con, sprintf("SELECT * FROM read_bcf(%s,index_path:=%s,tidy_format:=%s)",
        dbQuoteString(con, empty), dbQuoteString(con, paste0(empty, ".index.", kind)), tolower(tidy)))
      expect_equal(result, expected[FALSE, ])
    }
  }
  path <- system.file("extdata", "bcf_scan_contigs.partial.vcf.gz", package = "Rduckhts")
  expect_silent(rduckhts_bcf(con, "partial_vcf", path,
                            index_path = paste0(path, ".index.tbi"), overwrite = TRUE))
  expect_equal(dbGetQuery(con, "SELECT * FROM partial_vcf ORDER BY ALL"),
               dbGetQuery(con, "SELECT * FROM expected_sites ORDER BY ALL"))
}

test_bcf_literal_contigs <- function() {
  con <- rduckhts_connect()
  on.exit(dbDisconnect(con, shutdown = TRUE))
  dbExecute(con, "CREATE TABLE expected(CHROM VARCHAR, POS BIGINT, ID VARCHAR, REF VARCHAR,
    ALT VARCHAR[], QUAL DOUBLE, FILTER VARCHAR[], FORMAT_GT_S1 VARCHAR)")
  dbExecute(con, "INSERT INTO expected VALUES
    ('chr1',150,'ordinary','A',['C'],60,['PASS'],'0/1'),
    ('chr1:100-200',10,'literal','G',['T'],50,['PASS'],'1|1'),
    ('chr1:100-200',10,'literal','G',['T'],50,['PASS'],'1|1')")
  dbExecute(con, "CREATE VIEW expected_tidy AS SELECT * EXCLUDE(FORMAT_GT_S1),
    'S1' AS SAMPLE_ID, FORMAT_GT_S1 AS FORMAT_GT FROM expected")
  for (threads in c(1L, 4L)) for (tidy in c(FALSE, TRUE)) {
    dbExecute(con, sprintf("SET threads=%d", threads))
    expected <- dbGetQuery(con, paste("SELECT * FROM", if (tidy) "expected_tidy" else "expected",
                                     "ORDER BY ALL"))
    selected <- expected[expected$CHROM == "chr1:100-200", ]
    rownames(selected) <- NULL
    for (suffix in c("bcf", "vcf.gz")) {
      path <- system.file("extdata", paste0("bcf_literal_contigs.", suffix),
                          package = "Rduckhts", mustWork = TRUE)
      source <- sprintf("read_bcf(%s,tidy_format:=%s,decompression_threads:=0",
                        dbQuoteString(con, path), tolower(tidy))
      for (mode in c("auto", "sequential")) {
        reader <- paste0(source, ",scan_mode:=", dbQuoteString(con, mode), ")")
        expect_equal(dbGetQuery(con, paste("SELECT * FROM", reader, "ORDER BY ALL")), expected)
        expect_equal(dbGetQuery(con, paste("SELECT count(*) AS n FROM", reader))$n[[1]], 3)
        expect_equal(dbGetQuery(con, paste("SELECT CHROM,POS FROM", reader, "ORDER BY ALL")),
                     expected[c("CHROM", "POS")])
      }
      reader <- paste0(source, ",region:='{chr1:100-200}:1-20,{chr1:100-200}:10-10')")
      expect_equal(dbGetQuery(con, paste("SELECT * FROM", reader, "ORDER BY ALL")), selected)
      expect_error(dbGetQuery(con, paste0("SELECT * FROM ", source, ",region:='chr1:100-200')")),
                   pattern = "region list: invalid item")
    }
  }
}

test_bcf_literal_contigs()
test_bcf_contig_dictionary()
test_bcf_index_snapshot()
test_bcf_record_cache()
test_bcf_string_format_lists()
test_bcf_projection_sql()
test_bcf_materialization()