inst/tinytest/test_intervals.R

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

library(tinytest)
library(DBI)

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

  fasta_path <- system.file("extdata", "ce.fa", package = "Rduckhts")
  bed_path <- system.file("extdata", "targets.bed", package = "Rduckhts")
  fasta_index_path <- tempfile("duckhts_ce_", fileext = ".fai")

  expect_true(file.exists(fasta_path))
  expect_true(file.exists(bed_path))
  expect_true(rduckhts_fasta_index(
    con,
    fasta_path,
    index_path = fasta_index_path
  )$success[1])

  expect_silent(rduckhts_bed(con, "targets", bed_path, overwrite = TRUE))
  expect_equal(
    DBI::dbGetQuery(con, "SELECT count(*) AS n FROM targets")$n[1],
    4
  )

  typed_row <- DBI::dbGetQuery(
    con,
    paste(
      "SELECT chrom, start, \"end\", name, score, strand, thick_start, block_count",
      "FROM targets ORDER BY chrom, start LIMIT 1"
    )
  )
  expect_equal(typed_row$chrom[1], "CHROMOSOME_I")
  expect_equal(typed_row$start[1], 0)
  expect_equal(typed_row$end[1], 10)
  expect_equal(typed_row$name[1], "target1")
  expect_equal(typed_row$score[1], "100")
  expect_equal(typed_row$strand[1], "+")
  expect_equal(typed_row$thick_start[1], 0)
  expect_equal(typed_row$block_count[1], 2)

  nuc <- rduckhts_fasta_nuc(
    con,
    fasta_path,
    bed_path = bed_path,
    index_path = fasta_index_path
  )
  first_nuc <- nuc[nuc$chrom == "CHROMOSOME_I" & nuc$start == 0, , drop = FALSE]
  expect_equal(nrow(first_nuc), 1)
  expect_equal(first_nuc$pct_at[1], 0.4)
  expect_equal(first_nuc$pct_gc[1], 0.6)
  expect_equal(first_nuc$num_a[1], 2)
  expect_equal(first_nuc$num_c[1], 4)
  expect_equal(first_nuc$num_g[1], 2)
  expect_equal(first_nuc$num_t[1], 2)
  expect_equal(first_nuc$seq_len[1], 10)

  bins <- rduckhts_fasta_nuc(
    con,
    fasta_path,
    bin_width = 10,
    region = "CHROMOSOME_I:1-20",
    index_path = fasta_index_path
  )
  expect_equal(nrow(bins), 2)
  expect_equal(sum(bins$seq_len), 20)

  seq_nuc <- rduckhts_fasta_nuc(
    con,
    fasta_path,
    bed_path = bed_path,
    index_path = fasta_index_path,
    include_seq = TRUE
  )
  first_seq <- seq_nuc[
    seq_nuc$chrom == "CHROMOSOME_I" & seq_nuc$start == 0,
    "seq",
    drop = TRUE
  ]
  expect_equal(first_seq[1], "GCCTAAGCCT")

  # Deliberately malformed FAI metadata must fail, including after a valid BED row.
  bad_fasta <- system.file("extdata", "faidx_invalid.fa", package = "Rduckhts")
  bad_bed <- system.file("extdata", "faidx_invalid.bed", package = "Rduckhts")
  expect_error(rduckhts_fasta_nuc(con, bad_fasta, bin_width = 1,
    region = "offset_wrap:5-5", include_seq = TRUE), "fasta_nuc: failed to fetch")
  expect_error(rduckhts_fasta_nuc(con, bad_fasta, bin_width = 20,
    region = "truncated:1-20"), "fasta_nuc: failed to fetch")
  expect_error(rduckhts_fasta_nuc(con, bad_fasta, bed_path = bad_bed,
    include_seq = TRUE), "fasta_nuc: failed to fetch")
  recovered <- rduckhts_fasta_nuc(con, bad_fasta, bin_width = 4,
    region = "ok:1-4", include_seq = TRUE)
  expect_equal(recovered$seq, "AcGT")
  expect_equal(recovered$seq_len, 4)

  n_fasta_path <- system.file("extdata", "nuc_with_n.fa", package = "Rduckhts")
  n_bed_path <- system.file("extdata", "nuc_with_n.bed", package = "Rduckhts")
  n_fasta_index_path <- tempfile("duckhts_nuc_n_", fileext = ".fai")
  expect_true(file.exists(n_fasta_path))
  expect_true(file.exists(n_bed_path))
  expect_true(rduckhts_fasta_index(
    con,
    n_fasta_path,
    index_path = n_fasta_index_path
  )$success[1])

  n_nuc <- rduckhts_fasta_nuc(
    con,
    n_fasta_path,
    bed_path = n_bed_path,
    index_path = n_fasta_index_path
  )
  expect_equal(n_nuc$pct_at[1], 2 / 6)
  expect_equal(n_nuc$pct_gc[1], 4 / 6)
  expect_equal(n_nuc$num_a[1], 1)
  expect_equal(n_nuc$num_c[1], 2)
  expect_equal(n_nuc$num_g[1], 2)
  expect_equal(n_nuc$num_t[1], 1)
  expect_equal(n_nuc$num_n[1], 4)
  expect_equal(n_nuc$seq_len[1], 10)

  gz_path <- tempfile("targets", fileext = ".bed.gz")
  tbi_path <- paste0(gz_path, ".tbi")
  expect_true(rduckhts_bgzip(
    con,
    bed_path,
    output_path = gz_path,
    threads = 1,
    keep = TRUE,
    overwrite = TRUE
  )$success[1])
  expect_true(rduckhts_tabix_index(
    con,
    gz_path,
    preset = "bed",
    index_path = tbi_path,
    threads = 1
  )$success[1])
  expect_silent(rduckhts_bed(
    con,
    "targets_idx",
    gz_path,
    region = "CHROMOSOME_I:1-20",
    index_path = tbi_path,
    overwrite = TRUE
  ))
  expect_equal(
    DBI::dbGetQuery(con, "SELECT count(*) AS n FROM targets_idx")$n[1],
    2
  )

  bed_count_only <- DBI::dbGetQuery(
    con,
    sprintf(
      "SELECT COUNT(*) AS n FROM read_bed('%s', index_path := '%s')",
      gz_path,
      tbi_path
    )
  )
  expect_equal(bed_count_only$n[1], 4)
}

test_interval_readers()

test_bed_error_policy <- function() {
  con <- rduckhts_connect()
  on.exit(dbDisconnect(con, shutdown = TRUE))
  path <- system.file("extdata", "bed_error_policy.bed", package = "Rduckhts")
  expect_true(file.exists(path))
  expect_error(rduckhts_bed(con, "bed_default", path), "fewer than 3")
  expect_error(rduckhts_bed(con, "bed_invalid", path, error_policy = "other"),
    "error_policy must be")
  expect_silent(rduckhts_bed(con, "bed_skip", path, error_policy = "skip"))
  expect_equal(DBI::dbGetQuery(con, "SELECT count(*) AS n FROM bed_skip")$n, 5)
  expect_silent(rduckhts_bed(con, "bed_report", path,
    error_policy = "report", scan_mode = "sequential"))
  errors <- DBI::dbGetQuery(con,
    "SELECT line_number, raw_line FROM bed_report WHERE error IS NOT NULL ORDER BY line_number")
  # These are physical lines 5 and 10 of the bundled text fixture.
  expect_equal(errors$line_number, c(5, 10))
  expect_equal(errors$raw_line, c("chr1\t7", "chr2"))
  expect_equal(DBI::dbGetQuery(con, "SELECT count(*) AS n FROM bed_report")$n, 7)
  # An invalid policy is rejected before overwrite = TRUE drops the existing table.
  expect_error(rduckhts_bed(con, "bed_skip", path, overwrite = TRUE,
    error_policy = "other"), "error_policy must be")
  expect_true(DBI::dbExistsTable(con, "bed_skip"))
  # A combination the SQL binder rejects (report with a region) also leaves it intact.
  n_before <- DBI::dbGetQuery(con, "SELECT count(*) AS n FROM bed_skip")$n
  expect_error(rduckhts_bed(con, "bed_skip", path, overwrite = TRUE,
    region = "chr1:1-10", error_policy = "report"), "physical line numbers")
  expect_equal(DBI::dbGetQuery(con, "SELECT count(*) AS n FROM bed_skip")$n, n_before)
  expect_error(rduckhts_bed(con, "bed_region", path,
    region = "chr1:1-10", error_policy = "report"), "physical line numbers")
}

test_bed_error_policy()