Package {Rduckhts}

DuckDB HTS File Reader Extension for R

Description

The Rduckhts package provides an interface to the DuckDB HTS (High Throughput Sequencing) file reader extension from within R. It enables reading common bioinformatics file formats such as VCF/BCF, SAM/BAM/CRAM, FASTA, FASTQ, GFF, GTF, and tabix-indexed files directly from R using SQL queries via DuckDB.

Author(s)

DuckHTS Contributors

References

https://github.com/RGenomicsETL/duckhts

See Also

Useful links:

• https://github.com/RGenomicsETL/duckhts

• https://rgenomicsetl.r-universe.dev/Rduckhts

• Report bugs at https://github.com/RGenomicsETL/duckhts/issues

Detect Complex Types in DuckDB Table

Description

Identifies columns in a DuckDB table that contain complex types (ARRAY or MAP) that will be returned as R lists.

Usage

detect_complex_types(con, table_name)

Arguments

Value

A data frame with columns that have complex types, showing column_name, column_type, and a description of R type.

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
bcf_path <- system.file("extdata", "vcf_file.bcf", package = "Rduckhts")
rduckhts_bcf(con, "variants", bcf_path, overwrite = TRUE)
complex_cols <- detect_complex_types(con, "variants")
print(complex_cols)
dbDisconnect(con, shutdown = TRUE)

DuckDB to R Type Mappings

Description

The mapping covers the most common data types used in HTS file processing:

• BIGINT <-> double (not integer due to 64-bit overflow protection)

• DOUBLE <-> numeric/double

• VARCHAR <-> character/string

• BOOLEAN <-> logical

• ARRAY types (e.g., VARCHAR[], BIGINT[]) <-> list

• MAP types (e.g., MAP(VARCHAR, VARCHAR)) <-> data.frame

Important notes:

• 64-bit integers (BIGINT, UBIGINT) become double to prevent overflow

• DATE/TIME values return as Unix epoch numbers (double)

• MAP types become data frames with 'key' and 'value' columns

• ARRAY types become vectors (which are lists in R terminology)

Usage

duckdb_type_mappings()

Details

Returns a named list mapping between DuckDB and R data types. This is useful for understanding type conversions when reading HTS files or when specifying column types in tabix functions.

Value

A named list with two elements:

duckdb_to_r

Named character vector mapping DuckDB types to R types

r_to_duckdb

Named character vector mapping R types to DuckDB types

Examples

mappings <- duckdb_type_mappings()
mappings$duckdb_to_r["BIGINT"]
mappings$r_to_duckdb["integer"]

Bootstrap the duckhts extension sources into the R package

Description

Copies extension source files from the parent duckhts repository into inst/duckhts_extension/ so the R package becomes self-contained. Run this before R CMD build to prepare the source tarball.

Usage

duckhts_bootstrap(repo_root = NULL)

Arguments

Value

Invisibly returns the destination directory.

Retired manual DuckHTS extension builder

Description

The former manual builder duplicated the package configure path and could mutate an installed package tree. Build Rduckhts from a source tarball with R CMD build and install that tarball instead.

Usage

duckhts_build(build_dir = NULL, make = NULL, force = FALSE, verbose = TRUE)

Arguments

Value

This function always raises an error.

Load the duckhts extension into a DuckDB connection

Description

Compatibility wrapper around rduckhts_connect() and rduckhts_load(). Prefer the rduckhts_* functions in new code.

Usage

duckhts_load(con = NULL, extension_path = NULL)

Arguments

Value

The DuckDB connection (invisibly).

Extract Array Elements Safely

Description

Helper function to safely extract elements from DuckDB arrays (returned as R lists) with proper error handling.

Usage

extract_array_element(array_col, index = NULL, default = NA)

Arguments

Value

The array element at the specified index, or full array if index is NULL

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
bcf_path <- system.file("extdata", "vcf_file.bcf", package = "Rduckhts")
rduckhts_bcf(con, "variants", bcf_path, overwrite = TRUE)
data <- dbGetQuery(con, "SELECT ALT FROM variants LIMIT 5")
first_alt <- extract_array_element(data$ALT, 1)
all_alts <- extract_array_element(data$ALT)
dbDisconnect(con, shutdown = TRUE)

Extract MAP Keys and Values

Description

Helper function to work with DuckDB MAP data (returned as data frames). Can extract keys, values, or search for specific key-value pairs.

Usage

extract_map_data(map_col, operation = "keys", default = NA)

Arguments

Value

Extracted data based on the operation

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
gff_path <- system.file("extdata", "gff_file.gff.gz", package = "Rduckhts")
rduckhts_gff(con, "annotations", gff_path, attributes_map = TRUE, overwrite = TRUE)
data <- dbGetQuery(con, "SELECT attributes FROM annotations LIMIT 5")
keys <- extract_map_data(data$attributes, "keys")
name_values <- extract_map_data(data$attributes, "Name")
dbDisconnect(con, shutdown = TRUE)

Normalize R Data Types to DuckDB Types for Tabix

Description

Normalizes R data type names to their corresponding DuckDB types for use with tabix readers. This function handles common R type name variations and maps them to appropriate DuckDB column types.

Usage

normalize_tabix_types(types)

Arguments

Details

The function performs the following normalizations:

• Integer types (integer, int, int32, int64) -> BIGINT

• Numeric types (numeric, double, float) -> DOUBLE

• Character types (character, string, chr) -> VARCHAR

• Logical types (logical, bool, boolean) -> BOOLEAN

• Other types -> Converted to uppercase as-is

If an empty vector is provided, it returns the empty vector unchanged.

Value

A character vector of normalized DuckDB type names suitable for tabix columns.

See Also

rduckhts_tabix for using normalized types with tabix readers, duckdb_type_mappings for the complete type mapping table.

Examples

normalize_tabix_types(c("integer", "character", "numeric"))
normalize_tabix_types(c("int", "string", "float"))

Project ancestry directly from an indexed BAM or CRAM

Description

The panel must be a committed, versioned canonical panel intersected with the reference products. It has Somalier panel columns assembly, site_index, region, position, allele_a, and allele_b. Each count is extracted by the DuckHTS panel reader on the caller's connection. The count relation and frequency view are discarded on return; no private connection is created.

Usage

rduckhts_ancestry_bam(
  con,
  source_path,
  sample_id,
  reference_path,
  panel_table,
  reference_table,
  loadings_table,
  correction_table = NULL,
  frequency_method = c("allele_fraction", "called_genotype"),
  min_depth = 7,
  min_cor = 0.4,
  ...
)

Arguments

Value

Proportions and matching audit with frequency_method and min_depth.

Estimate ancestry from a materialized 'read_geno' relation

Description

Converts each diploid, biallelic call to ALT dosage / 2. A haploid, partially missing, or multiallelic call has NULL dosage and is counted in the missing-site audit. Original sample indices identify the calls unless a 'read_bcf_samples' relation is provided for original-header sample names. That relation must give every selected sample index one distinct non-null name.

Usage

rduckhts_ancestry_geno(
  con,
  geno_table,
  reference_table,
  loadings_table,
  correction_table,
  samples_table = NULL,
  sum_to_one = TRUE,
  min_cor = 0.4,
  non_reference_only
)

Arguments

Value

A row per sample and group with the matching audit.

Create a deterministic ancestry site panel for BAM/CRAM counts

Description

Sites are biallelic unambiguous SNVs present in every reference group and every PC loading. Positions are positive one-based whole numbers and all group alleles must be present. With 'candidate_table', only matching candidate loci and alleles are retained. Otherwise one site is selected per chromosome/spaced genomic window. Results are sorted by chromosome and position and assigned dense zero-based panel ordinals. Keep the returned panel SHA-256 and the reference release with the materialized panel as its versioned identity.

Usage

rduckhts_ancestry_panel(
  con,
  reference_table,
  loadings_table,
  table_name,
  assembly,
  candidate_table = NULL,
  spacing_bp = 5000,
  max_sites = 17000,
  overwrite = FALSE
)

Arguments

Value

A one-row data frame with panel SHA-256 and selected site count.

Estimate projected ancestry proportions from frequency or dosage relations

Description

Input columns are sample_id, chromosome, position, allele_a, allele_b and frequency (or dosage), describing allele_b. Positions are positive one-based whole numbers. Reference data may be a keyed wide relation with consecutive PC1..PCn columns and one frequency column per group, or a long relation with group_id and frequency plus a keyed pc/loading relation. Corrections contain pc and coefficient. Groups (1..30), PCs (1..64), frequencies and loadings must be complete and finite at each contributing reference locus; group frequencies must lie in [0, 1]. Duplicate loci matching the input are rejected. Long references require complete keyed rows during pivoting. Inputs and references must share an assembly, uppercase biallelic SNV alleles and reference orientation. Chromosome names are compared with 'duckhts_contig_key()': one leading 'chr' is removed in any letter case, 'M' and 'MT' become 'MT', and 'X' and 'Y' are uppercased. Any other name must match byte for byte, so '01' does not match '1'. A reference, wide or long, whose chromosome column has an integer type compares that key as an integer, so 'chr01' joins '1'. Input duplicates at a sample/locus, palindromic alleles and missing frequencies are dropped; other alleles match directly, reversed, strand complemented or both. Reversed alleles use 1-frequency. Missing genotypes are not imputed. Audit counts include dropped physical input rows. The solver and correlation gates use full precision; returned proportions are rounded to seven decimals. Aligned scratch Parquet in 'tempdir()' is removed on return. Correlation failures return NULL proportions and a status.

Usage

rduckhts_ancestry_proportions(
  con,
  input_table,
  reference_table,
  loadings_table,
  correction_table = NULL,
  input_kind = c("frequency", "dosage"),
  sum_to_one = TRUE,
  min_cor = 0.4,
  table_name = NULL,
  overwrite = FALSE,
  group_ids = NULL
)

Arguments

Value

A row per sample and reference group, with status and matching audit. Data-frame results carry 'aligned_bytes', the size of temporary aligned Parquet written during the call.

Create SAM/BAM/CRAM Table

Description

Creates a DuckDB table from SAM, BAM, or CRAM files using the DuckHTS extension.

Usage

rduckhts_bam(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  reference = NULL,
  standard_tags = FALSE,
  auxiliary_tags = FALSE,
  sequence_encoding = NULL,
  quality_representation = NULL,
  cigar_representation = NULL,
  scan_mode = NULL,
  decompression_threads = 2,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
bam_path <- system.file("extdata", "range.bam", package = "Rduckhts")
rduckhts_bam(con, "reads", bam_path, overwrite = TRUE)
dbGetQuery(con, "SELECT COUNT(*) FROM reads WHERE FLAG & 4 = 0")
dbDisconnect(con, shutdown = TRUE)

Native BAM/CRAM BED Regional Coverage Summary

Description

Computes samtools coverage-like regional summaries for BAM or CRAM input over a BED target set, with DuckHTS-specific pre/post-filter and strand-aware post-filter outputs.

Usage

rduckhts_bam_bed_coverage(
  con,
  path,
  bed_path,
  reference = NULL,
  index_path = NULL,
  bed_index_path = NULL,
  mapq = 0,
  min_baseq = 0,
  min_read_len = 0,
  require_flags = 0,
  exclude_flags = 1796,
  min_depth = 1,
  max_depth = 1e+06,
  decompression_threads = 0,
  fragment_mode = FALSE,
  strand_outputs = TRUE,
  processing_threads = 0
)

Arguments

Value

A data frame with one row per BED interval and pre/post regional summaries

Native Fixed-Width BAM/CRAM Bin Counts

Description

Count read starts into fixed-width genomic bins with optional duplicate handling and optional per-bin GC and MAPQ summary statistics.

Usage

rduckhts_bam_bin_counts(
  con,
  path,
  bin_width,
  chrom = NULL,
  include_unmapped = FALSE,
  reference = NULL,
  index_path = NULL,
  mapq = 0,
  require_flags = 0,
  exclude_flags = 0,
  rmdup = "none",
  stats = NULL
)

Arguments

Value

A data frame with one row per fixed-width bin across the selected contig span, including zero-count bins, plus total, forward, reverse, and optional GC/MAPQ summary columns

Convert SAM/BAM/CRAM reader output to Parquet with DuckHTS metadata

Description

Thin DBI wrapper around extension macro 'duckhts_bam_convert_parquet_sql(...)'.

Usage

rduckhts_bam_convert_parquet(
  con,
  path,
  output,
  columns = NULL,
  region = NULL,
  index_path = NULL,
  reference = NULL,
  standard_tags = FALSE,
  auxiliary_tags = FALSE,
  sequence_encoding = NULL,
  quality_representation = NULL,
  cigar_representation = NULL,
  decompression_threads = 2,
  where = NULL,
  compression = "zstd",
  row_group_size = 100000L,
  partition_by = NULL,
  include_metadata = TRUE,
  header_text = NULL,
  metadata = NULL,
  metadata_json_file = NULL,
  write_format_version = "1",
  overwrite = FALSE
)

Arguments

Value

Invisibly returns 'output'.

Build BAM or CRAM Index

Description

Builds a BAM or CRAM index using the DuckHTS extension.

Usage

rduckhts_bam_index(con, path, index_path = NULL, min_shift = 0, threads = 4)

Arguments

Value

A data frame with 'success', 'index_path', and 'index_format'

BAM/CRAM Insert-Size Counts

Description

Reduces eligible records into sparse orientation and absolute-TLEN counts using Fulcrum Genomics riker 0.4.1 semantics (MIT). Only non-first segments count; QC-fail records are excluded. The unsuffixed wrapper uses the native reducer; the '_sql' wrapper uses the SQL macro. The SQL arm projects binary CIGAR and flag/coordinate columns only; CRAM decoding also requests sequence and tags internally for long-CIGAR restoration. The native arm avoids these fields. For weighted metrics, materialize the counts and call ‘duckhts_insert_size_metrics(’counts_table')' in SQL.

Usage

rduckhts_bam_insert_size_counts(
  con,
  path,
  reference = NULL,
  include_duplicates = FALSE,
  threads = 1,
  decompression_threads = 0,
  scan_mode = c("auto", "indexed", "sequential"),
  index_path = NULL,
  reference_index_path = NULL
)

rduckhts_bam_insert_size_counts_sql(
  con,
  path,
  reference = NULL,
  include_duplicates = FALSE,
  threads = 1,
  decompression_threads = 0,
  scan_mode = c("auto", "sequential"),
  index_path = NULL
)

Arguments

Value

Data frame with pair_orientation, insert_size and read_pairs. Rows are unordered; absent bins are omitted.

Count Agreement Between Overlapping Mates

Description

Compares read 1 and read 2 of one fragment at each reference position where both mates have an aligned base, and counts the pairs of base qualities and bases. A real allele sits on the molecule and both mates show it; a sequencing error usually shows in one mate only. So the disagreement rate is about twice the instrument error rate, with no reference, mask or contamination model. Errors made before sequencing, such as PCR or damage errors, show in both mates and are not seen.

Usage

rduckhts_bam_mate_overlap_counts(
  con,
  path,
  region = NULL,
  index_path = NULL,
  reference = NULL,
  reference_index_path = NULL,
  min_mapq = 20,
  require_flags = 0,
  exclude_flags = 3844,
  max_pending = 1e+06
)

Arguments

Details

Alignments must be in coordinate order. The first mate of a pair is held until its mate arrives; more than 'max_pending' mates held at once is an error.

Value

A data frame with 'read1_quality' and 'read2_quality' ('NA' when the read stores no qualities), 'read1_base' and 'read2_base' on the forward strand of the reference, and 'positions'. Read 1 is the mate with flag 0x40. Rows with different bases are the disagreements.

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
bam <- system.file("extdata", "empty-tids.bam", package = "Rduckhts")
counts <- rduckhts_bam_mate_overlap_counts(con, bam)
sum(counts$positions[counts$read1_base != counts$read2_base]) / sum(counts$positions)
dbDisconnect(con, shutdown = TRUE)

Count Aligned Bases Against the Reference

Description

Counts the aligned bases of a SAM, BAM or CRAM file against the reference by mate, cycle, base quality and substitution. With a mask of the known variants of the sample or its population, the mismatches that remain are mostly errors of the read, the library or the alignment. A variant that the mask does not hold is still counted, so their rate is an upper bound on the error rate.

Usage

rduckhts_bam_mismatch_counts(
  con,
  path,
  reference,
  region = NULL,
  mask = NULL,
  index_path = NULL,
  reference_index_path = NULL,
  mask_index_path = NULL,
  min_mapq = 20,
  require_flags = 0,
  exclude_flags = 3844,
  indel_flank = 5
)

Arguments

Details

Contig names of the alignment header, the reference and the mask are compared byte for byte. A contig that the reference or the mask does not know is an error.

Value

A data frame with 'mate' (1 or 2 for a paired read, 0 otherwise), 'cycle', 'base_quality' ('NA' when the read stores no qualities), 'reference_base' and 'read_base' on the strand of the read, and 'bases'. Rows with equal bases are the matches.

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
extdata <- function(name) system.file("extdata", name, package = "Rduckhts")
counts <- rduckhts_bam_mismatch_counts(
  con, extdata("bam_mismatch.bam"), extdata("bam_mismatch.fa"),
  mask = extdata("bam_mismatch.mask.bcf")
)
sum(counts$bases[counts$reference_base != counts$read_base]) / sum(counts$bases)
dbDisconnect(con, shutdown = TRUE)

Read multiple BAM/SAM files into a DuckDB table

Description

Read and combine multiple BAM/SAM files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_bam_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  reference = NULL,
  standard_tags = FALSE,
  auxiliary_tags = FALSE,
  sequence_encoding = NULL,
  quality_representation = NULL,
  cigar_representation = NULL,
  scan_mode = NULL,
  decompression_threads = 2,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

Create VCF/BCF Table

Description

Creates a DuckDB table from a VCF or BCF file using the DuckHTS extension. This follows the RBCFTools pattern of creating a table that can be queried.

Usage

rduckhts_bcf(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  tidy_format = FALSE,
  additional_csq_column_types = NULL,
  scan_mode = NULL,
  decompression_threads = 0,
  decode_error_policy = "null",
  overwrite = FALSE,
  samples = NULL,
  regions_var = NULL,
  include_file_offset = FALSE
)

Arguments

Value

Invisible TRUE on success

Examples

library(DBI)
library(duckdb)

con <- rduckhts_connect()
bcf_path <- system.file("extdata", "vcf_file.bcf", package = "Rduckhts")
rduckhts_bcf(con, "variants", bcf_path, overwrite = TRUE)
dbGetQuery(con, "SELECT * FROM variants LIMIT 2")
dbDisconnect(con, shutdown = TRUE)

Convert VCF/BCF reader output to Parquet with DuckHTS metadata

Description

'rduckhts_bcf_convert_parquet()' is a thin DBI wrapper around the extension macro 'duckhts_bcf_convert_parquet_sql(...)'. The extension macro builds the 'COPY ... TO ... (FORMAT PARQUET, KV_METADATA ...)' statement, including DuckHTS write-format metadata, optional user metadata maps, selected columns, optional SQL filter text, 'tidy_format', and the VCF header under the 'vcf_header' key. The wrapper executes the generated statement.

Usage

rduckhts_bcf_convert_parquet(
  con,
  path,
  output,
  columns = NULL,
  region = NULL,
  index_path = NULL,
  tidy_format = FALSE,
  additional_csq_column_types = NULL,
  decompression_threads = 0,
  where = NULL,
  compression = "zstd",
  row_group_size = 100000L,
  partition_by = NULL,
  include_metadata = TRUE,
  header_text = NULL,
  metadata = NULL,
  metadata_json_file = NULL,
  write_format_version = "1",
  overwrite = FALSE
)

Arguments

Value

Invisibly returns 'output'.

Build VCF or BCF Index

Description

Builds a TBI or CSI index for a VCF/BCF file using the DuckHTS extension.

Usage

rduckhts_bcf_index(con, path, index_path = NULL, min_shift = NULL, threads = 4)

Arguments

Value

A data frame with 'success', 'index_path', and 'index_format'

Read multiple VCF/BCF files into a DuckDB table

Description

Read and combine multiple VCF/BCF files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_bcf_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  tidy_format = FALSE,
  additional_csq_column_types = NULL,
  scan_mode = NULL,
  decompression_threads = 0,
  .params = NULL,
  overwrite = FALSE,
  samples = NULL
)

Arguments

Value

Invisible TRUE on success.

Read the VCF/BCF Sample Catalog

Description

Return one row per selected sample with its original-header zero-based 'sample_index' and 'sample_name'. Join this relation to 'read_geno()' calls from the same unchanged file; names are not duplicated into every call.

Usage

rduckhts_bcf_samples(con, path, samples = NULL)

Arguments

Value

A data frame with 'sample_index' and 'sample_name' columns.

Normalize Variant Alleles with bcftools-style Semantics

Description

Applies the DuckHTS 'duckhts_bcftools_norm(...)' table macro to rows from a SQL query or table expression. Input rows must expose chromosome, 1-based position, reference allele, and alternate allele columns. Alternate alleles may be supplied either as a comma-delimited 'VARCHAR' or as a 'VARCHAR[]' list, matching the common DuckDB representations used by plain tables and 'read_bcf(...)'.

Usage

rduckhts_bcftools_norm(
  con,
  query,
  fasta_ref,
  chrom_col = "chrom",
  pos_col = "pos",
  ref_col = "ref",
  alt_col = "alt",
  split_multiallelic = FALSE,
  end_pos_col = NULL,
  svlen_col = NULL,
  fasta_index_path = NULL,
  gzi_path = NULL
)

Arguments

Value

A data frame with the original columns plus 'pos_normed', 'end_pos_normed', 'ref_normed', 'alt_normed', 'normed', and 'norm_status'. In split mode the result additionally includes 'alt_index'.

Create BED Table

Description

Creates a DuckDB table from a BED file using the DuckHTS extension.

Usage

rduckhts_bed(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  scan_mode = NULL,
  overwrite = FALSE,
  error_policy = NULL
)

Arguments

Value

Invisible TRUE on success

Read multiple BED files into a DuckDB table

Description

Read and combine multiple BED files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_bed_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  scan_mode = NULL,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

BGZF Decompress a File

Description

Decompresses a BGZF file using the DuckHTS extension.

Usage

rduckhts_bgunzip(
  con,
  path,
  output_path = NULL,
  threads = 4,
  keep = TRUE,
  overwrite = FALSE
)

Arguments

Value

A data frame describing the created output file

BGZF Compress a File

Description

Compresses a plain file to BGZF using the DuckHTS extension.

Usage

rduckhts_bgzip(
  con,
  path,
  output_path = NULL,
  threads = 4,
  level = -1,
  keep = TRUE,
  overwrite = FALSE
)

Arguments

Value

A data frame describing the created BGZF file

Create a BigWig Signal Table

Description

Materializes stored zero-based, half-open BigWig intervals through the DuckHTS extension. Region filters use htslib's one-based inclusive syntax; a character vector is combined into one multi-region request and overlapping requests emit each stored interval once.

Usage

rduckhts_bigwig(
  con,
  table_name,
  path,
  region = NULL,
  blocks_per_iteration = 64L,
  overwrite = FALSE
)

Arguments

Value

Invisibly returns TRUE.

Create a DuckDB connection with bundled DuckHTS loaded

Description

Creates a DuckDB connection that explicitly permits loading the package-built DuckHTS extension, disables automatic installation and loading of unrelated known DuckDB extensions, and loads the bundled DuckHTS extension.

Usage

rduckhts_connect(
  dbdir = ":memory:",
  read_only = FALSE,
  bigint = "numeric",
  config = list(),
  extension_path = NULL
)

Arguments

Details

Current versions of the duckdb R package can disable extension loading on Linux builds that use a C++ standard library other than libstdc++. DuckHTS is compiled locally with the same package toolchain, so this helper explicitly enables extension loading for this package-owned connection. Package-owned connections use per-session DuckDB extension/secret storage by default when the installed duckdb version supports that setting.

DuckDB reuses a live database instance when the same file-backed dbdir is opened again, ignoring new driver and configuration settings. This helper rejects such reuse rather than silently weakening its connection policy. Use rduckhts_load() on the existing connection if it already permits unsigned, driver-level extension loading, or close the existing database instance before calling this helper.

The settings allow_unsigned_extensions, autoinstall_known_extensions, and autoload_known_extensions are controlled by this helper and override entries with those names in config. Other named DuckDB settings are passed through unchanged.

Value

A DuckDB connection with DuckHTS loaded. File-backed connections also have connection-local TEMP macros installed without modifying the database catalog.

Examples

con <- rduckhts_connect()
DBI::dbGetQuery(con, "SELECT duckhts_htslib_version() AS version")
DBI::dbDisconnect(con, shutdown = TRUE)

Fit Read Error and Contamination from Allele Read Counts

Description

Fits, for each sample, the read error, the share of reads from a second genome, the allele balance at heterozygous sites, two spreads, an artefact weight and an excess of homozygosity by maximum likelihood from allele read counts at sites of known population frequency. The input is the relation of [rduckhts_roh_counts()]; the fit is meant to supply that decode's 'seq_error' and 'contamination'.

Usage

rduckhts_count_error_fit(
  con,
  counts_table,
  block_bases = 1e+07,
  freq_bins = 64,
  max_depth = 1000,
  min_block_sites = 1000,
  table_name = NULL,
  overwrite = FALSE,
  max_cells = 5e+06,
  max_cell_bytes = 2^30
)

Arguments

Details

The likelihood depends on a site only through its two counts and its frequency, so the fit runs on a histogram of (block, frequency bin, depth, alt count) cells that DuckDB builds and the native code holds in a bounded buffer. Each block of 'block_bases' bases with at least 'min_block_sites' sites is refitted for contamination with the other parameters held, and 'contamination_sd' is the standard deviation of those block estimates divided by the square root of their number: an interval from the blocks, not from the likelihood, because sites are linked. Each such block is also refitted for the artefact weight alone, and 'artefact_weight_sd' is formed the same way; both are NA with fewer than three blocks. [rduckhts_roh_counts()] with 'fit_table' uses 'artefact_weight_sd' for its 'stable_fraction'. The same cells are also fitted with a parent or child as the second genome ('contamination_relative', 'log_likelihood_relative'), since the unrelated model gives a relative about half its share.

'homozygosity_excess' is a nuisance parameter: it also absorbs a difference between the sample's population and the frequencies given, so it is not inbreeding. A second genome of another population biases 'contamination' upward under the receiver's frequencies.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given) with one row per sample: 'sample', 'sites', 'reads', 'mean_depth', 'blocks', 'seq_error', 'contamination', 'contamination_sd', 'homozygosity_excess', 'allele_balance', 'spread_hom', 'spread_het', 'artefact_weight', 'artefact_weight_sd', 'log_likelihood', 'contamination_relative', 'log_likelihood_relative', 'status' ('ok', 'no_sites' for a sample whose sites were all skipped, 'few_sites', 'few_blocks', 'at_bound', 'no_convergence' or 'no_convergence_relative') and 'method'.

Examples

con <- rduckhts_connect()
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE counts AS SELECT 'S1' AS sample_id, '1' AS chrom,",
  "i * 2500 AS pos, 0.1 + 0.8 * ((i * 37) % 100) / 100.0 AS af,",
  "30 - ((i * 7) % 3) * 15 AS ref_count, ((i * 7) % 3) * 15 AS alt_count",
  "FROM range(1, 4001) AS t(i)"
))
rduckhts_count_error_fit(con, "counts")[, c("sample", "seq_error", "contamination", "status")]
DBI::dbDisconnect(con, shutdown = TRUE)

Estimate a Conditional Second-Genome Share from Low-Pass Counts

Description

Uses first-read and discordant-pair moments with fixed external calibration. This estimator is separate from rduckhts_count_error_fit(). It is suitable for low-pass, cell-free DNA or low-input DNA counts only under the supplied population frequencies and calibration; the share is not an unconditional biological fraction. Matched-frequency unrelated and IBD1 hypotheses can have identical moments and are reported as unresolved.

Usage

rduckhts_count_mixture_lowpass(
  con,
  counts_table,
  seq_error,
  artefact_weight,
  homozygosity_excess,
  block_bases = 1e+07,
  max_depth = 8,
  min_blocks = 30,
  max_blocks = 512,
  max_candidates = 16,
  max_state_bytes = 2^28,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Value

One row per sample, candidate and relationship model. IBD1 rows require exact candidate_af == sample_af at every input site, including zero-read and excluded sites. Columns include the bounded share (0 to 0.5), block-jackknife standard error and score-inverted limits, singles share, pair residual, relationship diagnostic, QC read/pair totals, calibration, score rank and status. The discordance ratio divides discordant pairs by SUM(choose(depth, 2) * sample_af * (1 - sample_af)). When table_name is supplied, returns invisible TRUE.

Correct Bin Read Counts for GC Content

Description

Gives each bin of each sample an expected read count from the median read rate of the sample's autosomal training bins in the same GC stratum, so that 'ratio = reads / expected' is about 1 for two copies, 0.5 for one and 1.5 for three. The output feeds [rduckhts_depth_region_fraction()] and [rduckhts_depth_sex_fraction()]. Normalization uses eligible autosomal bins with positive reads; empty bins retain zero ratios when their bin and stratum rules pass.

Usage

rduckhts_depth_bins(
  con,
  bins_table,
  min_mappability = 0.9,
  gc_step = 0.01,
  min_stratum_bins = 10,
  outlier_quantile = 0.01,
  max_ratio = 4,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Details

The GC stratum of a bin is 'floor(gc / gc_step)'. Training bins are autosomal bins that pass the bin rules, hold at least one read and have a read rate inside the 'outlier_quantile' and '1 - outlier_quantile' quantiles of the sample. A bin is kept unless it lacks a count or a GC value, has a mappability below 'min_mappability' (when the relation has a 'mappability' column), falls in a stratum with fewer than 'min_stratum_bins' training bins, or has a ratio above 'max_ratio'; 'reason' names the rule. Only ‘max_ratio' looks at a bin’s count, so a bin with no reads stays kept. Autosomal means a contig whose 'duckhts_contig_key()' is not X, Y or MT. Restrict contigs and remove problem bins in SQL first.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given) with one row per input row: 'sample_id', 'chrom', 'start', 'end', 'reads', 'gc', 'mappability', 'gc_stratum', 'stratum_bins', 'expected', 'ratio', 'log2_ratio', 'reason', 'kept' and 'method'.

See Also

[rduckhts_bam_bin_counts()], [rduckhts_fasta_nuc()]

Examples

con <- rduckhts_connect()
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE bins AS SELECT 'S1' AS sample_id, CAST(c AS VARCHAR) AS chrom,",
  "b * 1000000 AS start, (b + 1) * 1000000 AS \"end\",",
  "0.35 + 0.25 * ((b * 37 + c * 11) % 100) / 100.0 AS gc,",
  "CAST(600 + 4 * ((b * 13 + c * 7) % 11) AS BIGINT) AS reads",
  "FROM range(1, 23) t(c), range(0, 100) r(b)"
))
head(rduckhts_depth_bins(con, "bins")[, c("chrom", "start", "reads", "ratio", "kept")])
DBI::dbDisconnect(con, shutdown = TRUE)

Estimate the Share of Cells with a One-Copy Change of a Region

Description

For each sample and named region, compares the depth of the region's kept bins with the depth of the sample's other kept autosomal bins: 'ratio = (C_r / E_r) / (C_b / E_b)' with 'C' and 'E' the sums of reads and expected counts, and 'shift = 2 * (ratio - 1)', the share of cells with a one-copy gain (positive) or loss (negative) relative to two copies when the baseline is unchanged. A one-copy change in every cell gives 1 or -1.

Usage

rduckhts_depth_region_fraction(
  con,
  bins_table,
  regions_table,
  reference_table = NULL,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Details

Without 'reference_table', the interval treats bin counts as independent and quasi-Poisson, with the dispersion of the baseline bins, and is 'shift +- 1.959964 shift_se'. With 'reference_table', corrected bins of reference samples, the reported ‘shift' is the sample’s shift minus the mean shift of the references (leaving out the tested sample if it is a reference), and the interval uses Student's t with 'reference_samples - 1' degrees of freedom. Other changes in the sample move the baseline.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given) with one row per sample and region: 'sample_id', 'region', 'bins', 'baseline_bins', 'ratio', 'shift', 'shift_se', 'shift_low', 'shift_high', 'dispersion', 'reference_shift', 'reference_samples', 'status' ('ok', 'no_bins', 'no_baseline' or 'few_reference_samples') and 'method'.

Examples

con <- rduckhts_connect()
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE bins AS SELECT 'S1' AS sample_id, CAST(c AS VARCHAR) AS chrom,",
  "b * 1000000 AS start, (b + 1) * 1000000 AS \"end\",",
  "0.35 + 0.25 * ((b * 37 + c * 11) % 100) / 100.0 AS gc,",
  "CAST((600 + 4 * ((b * 13 + c * 7) % 11)) *",
  "CASE WHEN c = 8 THEN 1.1 ELSE 1 END AS BIGINT) AS reads",
  "FROM range(1, 23) t(c), range(0, 100) r(b)"
))
rduckhts_depth_bins(con, "bins", table_name = "corrected")
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE chromosomes AS SELECT CAST(c AS VARCHAR) AS region,",
  "CAST(c AS VARCHAR) AS chrom, 0 AS start, 100000000 AS \"end\" FROM range(1, 23) t(c)"
))
fractions <- rduckhts_depth_region_fraction(con, "corrected", "chromosomes")
fractions[fractions$region == "8", c("region", "shift", "shift_low", "shift_high")]
DBI::dbDisconnect(con, shutdown = TRUE)

Estimate the Share of a Second Genome of the Other Sex from X and Y Depth

Description

Treats each sample as a major genome and a second genome of the other sex with a DNA share 'f' below 0.5. 'x_copies' and 'y_copies' are twice the depth ratio of the kept X or Y bins to the kept autosomal bins. For an XX major genome 'x_copies = 2 - f' and 'y_copies = f'; for an XY one 'x_copies = 1 + f' and 'y_copies = 1 - f'. The major genome is XX when 'x_copies >= 1.5'. 'fraction_x' and 'fraction_y' are the two estimates, unclipped; 'fraction' combines them by inverse variance and 'agreement_z' compares them, since an X or Y copy change of the major genome moves them differently.

Usage

rduckhts_depth_sex_fraction(
  con,
  bins_table,
  reference_table = NULL,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Details

Without 'reference_table', one X copy is 1, two are 2, one Y copy is 1 and a genome without Y has no Y depth. With 'reference_table', its XY samples set the depth of one X and one Y copy, and its XX samples, when present, set two X copies and the Y background; the tested sample is left out of its own reference. Remove the pseudo-autosomal regions of X, and Y bins that are not single-copy, in SQL first.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given) with one row per sample: 'sample_id', 'autosomal_bins', 'x_bins', 'y_bins', 'dispersion', 'x_copies', 'x_copies_se', 'y_copies', 'y_copies_se', 'major_sex', 'fraction_x', 'fraction_x_se', 'fraction_y', 'fraction_y_se', 'fraction', 'fraction_se', 'fraction_low', 'fraction_high', 'agreement_z', 'reference_samples', 'status' ('ok', 'no_autosomal_bins', 'no_xy_reference', 'bad_reference', 'no_sex_bins', 'no_x_bins', 'no_y_bins' or 'x_y_disagree') and 'method'.

Examples

con <- rduckhts_connect()
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE bins AS SELECT 'S1' AS sample_id, chrom,",
  "b * 1000000 AS start, (b + 1) * 1000000 AS \"end\",",
  "0.35 + 0.25 * ((b * 37 + length(chrom) * 11) % 100) / 100.0 AS gc,",
  "CAST((600 + 4 * ((b * 13) % 11)) * CASE chrom WHEN 'X' THEN 0.95 WHEN 'Y' THEN 0.05",
  "ELSE 1 END AS BIGINT) AS reads",
  "FROM (SELECT CAST(c AS VARCHAR) AS chrom FROM range(1, 23) t(c)",
  "UNION ALL SELECT 'X' UNION ALL SELECT 'Y') g, range(0, 60) r(b)"
))
rduckhts_depth_bins(con, "bins", table_name = "corrected")
rduckhts_depth_sex_fraction(con, "corrected")[, c("major_sex", "fraction", "status")]
DBI::dbDisconnect(con, shutdown = TRUE)

Decode Tumour Fraction, Ploidy and Copy-Number Segments from Corrected Bins

Description

Fits the copy-number state model of ichorCNA (Adalsteinsson, Ha, Freeman et al. 2017; GavinHaLab/ichorCNA tag v0.6.1) to each sample's kept corrected bins: a hidden Markov model over copy numbers 1 to 'max_cn' with a Student's t emission on the log ratio, whose state means follow from the normal fraction 'n' and tumour ploidy 'phi', fitted by expectation maximisation with the upstream priors from every start of 'ploidies' x 'normal_fractions'. Each start gives one solution, ranked by its objective (data log-likelihood plus log prior). [rduckhts_depth_state_solutions()] returns the solutions; [rduckhts_depth_state_segments()] returns their Viterbi segments.

Usage

rduckhts_depth_state_solutions(
  con,
  bins_table,
  normal_fractions = c(0.5, 0.6, 0.7, 0.8, 0.9),
  ploidies = c(2, 3),
  max_cn = 7,
  txn_e = 0.9999999,
  txn_strength = 1e+07,
  max_iter = 20,
  min_segment_bins = 50,
  alt_fraction_threshold = 0.05,
  close_margin = 10,
  max_bins = 1e+06,
  table_name = NULL,
  overwrite = FALSE
)

rduckhts_depth_state_segments(
  con,
  bins_table,
  normal_fractions = c(0.5, 0.6, 0.7, 0.8, 0.9),
  ploidies = c(2, 3),
  max_cn = 7,
  txn_e = 0.9999999,
  txn_strength = 1e+07,
  max_iter = 20,
  min_segment_bins = 50,
  alt_fraction_threshold = 0.05,
  close_margin = 10,
  max_bins = 1e+06,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Details

The model learns from autosomal bins (contigs whose 'duckhts_contig_key()' is not X, Y or MT) and decodes every kept bin. A kept bin without a ratio takes every state with likelihood 1, as upstream carries it. Upstream's guard sets the fraction to zero ('status' 'guard_zero') when the longest altered training segment has at most 'min_segment_bins' bins and at most 'alt_fraction_threshold' of the training bins are altered; a solution ranked below the best with an objective within 'close_margin' of it is 'close_to_best', the rest are 'ok'. A sample without copy-number change therefore reports no fraction, by the guard or by a fit whose normal fraction is near 1; in that fit the state means lie within the noise of each other and whole chromosomes can take any state, so its copy numbers are not evidence. Fraction and ploidy trade against each other, so all solutions are returned; a ploidy-3 explanation of a ploidy-2 genome usually ranks second. A genome with one altered region does not identify the fraction either: a one-copy gain in a share 'm' of cells and a two-copy gain in 'm / 2' give the same ratio, and the prior decides, as in upstream. The port and its kept upstream quirks are documented in the DuckHTS catalog entry and 'design/depth_fraction_lane.md'.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given). Solutions have one row per sample and start: 'sample_id', 'solution' (rank, 1 is best), 'start_normal', 'start_ploidy', 'normal_fraction', 'tumour_fraction', 'ploidy', 'log_likelihood', 'altered_fraction', 'iterations', 'status', 'segments', 'bins', 'upstream' and 'method'. Segments have one row per sample, solution and segment: the solution columns through 'status', then 'chrom', 'start', 'end', 'bins', 'state' (the copy number of the Viterbi state), ‘copy_number' (after upstream’s integer correction of high-level amplifications), 'median_log2_ratio', 'upstream' and 'method'.

Examples

con <- rduckhts_connect()
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE bins AS SELECT 'S1' AS sample_id, CAST(c AS VARCHAR) AS chrom,",
  "b * 1000000 AS start, (b + 1) * 1000000 AS \"end\", true AS kept,",
  "log2(CASE WHEN c = 8 THEN 1.15 ELSE 1 END) +",
  "0.05 * (((b * 37 + c * 11) % 100) / 50.0 - 1) AS log2_ratio",
  "FROM range(1, 23) t(c), range(0, CASE WHEN c = 8 THEN 60 ELSE 20 END) r(b)"
))
# Chromosome 8 has 60 bins so the gain passes the 50-bin guard. One start
# keeps the example short; the default grid has ten.
best <- rduckhts_depth_state_solutions(con, "bins", normal_fractions = 0.7, ploidies = 2)
best[, c("tumour_fraction", "ploidy", "status")]
segments <- rduckhts_depth_state_segments(con, "bins", normal_fractions = 0.7, ploidies = 2)
segments[segments$chrom == "8", c("chrom", "copy_number", "bins")]
DBI::dbDisconnect(con, shutdown = TRUE)

Detect FASTQ Quality Encoding

Description

Inspects a FASTQ file's observed quality ASCII range and reports compatible legacy encodings with a heuristic guessed encoding.

Usage

rduckhts_detect_quality_encoding(con, path, max_records = 10000)

Arguments

Value

A data frame with the detected quality encoding summary

Create FASTA Table

Description

Creates a DuckDB table from FASTA files using the DuckHTS extension.

Usage

rduckhts_fasta(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  gzi_path = NULL,
  sequence_encoding = NULL,
  scan_mode = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success

Build FASTA Index

Description

Builds a FASTA index (.fai) using the DuckHTS extension.

Usage

rduckhts_fasta_index(con, path, index_path = NULL)

Arguments

Value

A data frame with columns 'success' and 'index_path'

Read multiple FASTA files into a DuckDB table

Description

Read and combine multiple FASTA files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_fasta_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  gzi_path = NULL,
  sequence_encoding = NULL,
  scan_mode = NULL,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

Compute FASTA Interval Nucleotide Composition

Description

Computes bedtools nuc-style nucleotide composition over either a BED file or generated fixed-width bins.

Usage

rduckhts_fasta_nuc(
  con,
  path,
  bed_path = NULL,
  bin_width = NULL,
  region = NULL,
  index_path = NULL,
  gzi_path = NULL,
  bed_index_path = NULL,
  include_seq = FALSE
)

Arguments

Value

A data frame with interval composition statistics

Create FASTQ Table

Description

Creates a DuckDB table from FASTQ files using the DuckHTS extension.

Usage

rduckhts_fastq(
  con,
  table_name,
  path,
  mate_path = NULL,
  interleaved = FALSE,
  sequence_encoding = NULL,
  quality_representation = NULL,
  input_quality_encoding = NULL,
  scan_mode = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success

Read multiple FASTQ files into a DuckDB table

Description

Read and combine multiple FASTQ files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_fastq_multi(
  con,
  table_name,
  files,
  mate_path = NULL,
  interleaved = FALSE,
  sequence_encoding = NULL,
  quality_representation = NULL,
  input_quality_encoding = NULL,
  scan_mode = NULL,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

List DuckHTS Extension Functions

Description

Returns the package-bundled function catalog generated from the top-level functions.yaml manifest in the duckhts repository.

Usage

rduckhts_functions(category = NULL, kind = NULL)

Arguments

Value

A data frame describing the extension functions, including the DuckDB function name, kind, category, signature, return type, optional R helper wrapper, short description, and example SQL.

Examples

catalog <- rduckhts_functions()
subset(catalog, category == "Sequence UDFs", select = c("name", "description"))
subset(rduckhts_functions(kind = "table"), select = c("name", "r_wrapper"))

Create GenBank Feature Table

Description

Creates a DuckDB table from a GenBank flat file using the DuckHTS extension. Features are emitted in read_gff's column shape, so a GenBank record substitutes for a GFF without a schema change. Locations built with join() or order() give one row per segment in biological order, complement(...) sets strand "-", and the GFF3 phase of each CDS segment is carried from /codon_start across segments. Parent links a feature to the gene sharing its /locus_tag wherever that gene appears in the record, repeated qualifiers become one key with comma-joined values, and valueless qualifiers read true. Records stream one at a time, and a record without a terminating // or with a malformed location is an error naming the feature and line.

Usage

rduckhts_genbank(
  con,
  table_name = NULL,
  path,
  attributes_map = FALSE,
  overwrite = FALSE,
  attributes = NULL
)

Arguments

Value

Invisible TRUE on success

Write GenBank Sequence as FASTA

Description

Writes the ORIGIN sequence of each record in a GenBank flat file to a FASTA file using the DuckHTS extension. Records are written under the same name rduckhts_genbank reports as seqname, so feature coordinates land on the contig of that name; the DEFINITION follows the name without its trailing period, as in NCBI's FASTA export. The FASTA is written to a temporary file beside output_path and renamed into place only after the input has been read to a clean end, so a failure never leaves a partial output and an existing file is never lost. Records without an ORIGIN block are skipped; zero written records is an error.

Usage

rduckhts_genbank_to_fasta(
  con,
  path,
  output_path = NULL,
  line_width = 70,
  overwrite = FALSE
)

Arguments

Value

A data frame with columns success, output_path and records_written

Read Record-Major Genotypes

Description

Read typed GT/PS calls without repeating variant text per sample. Each record has a zero-based scan-local 'record_index' and a list of calls with original-header sample indices, nullable allele indices, per-slot phase bits, and nullable scalar phase sets. An absent GT has NULL allele/phase lists; a missing allele still occupies a slot. Phase bits follow HTSlib decoding, including its leading-slot convention for VCF versions before 4.4. PS cardinality excludes vector-end padding retained after sample selection.

Usage

rduckhts_geno(
  con,
  table_name = NULL,
  path,
  region = NULL,
  index_path = NULL,
  samples = NULL,
  non_reference_only = FALSE,
  scan_mode = "auto",
  decompression_threads = 0,
  decode_error_policy = "null",
  overwrite = FALSE,
  format_fields = NULL,
  raw_gt = FALSE,
  include_filter = FALSE,
  regions_var = NULL
)

Arguments

Details

Full scans preserve the input stream when assigning ordinals. Indexed regions use HTSlib's union order and start a new ordinal at zero. Use SQL 'ORDER BY record_index' when order matters; the ordinal is not a persistent file locator. Empty sample selection and sparse calls preserve every selected variant row.

Value

A data frame when 'table_name' is 'NULL', otherwise invisible 'TRUE'.

See Also

[rduckhts_bcf_samples()], [rduckhts_geno_sites()]

Examples

con <- rduckhts_connect()
path <- system.file("extdata", "geno_calls.bcf", package = "Rduckhts")
rduckhts_geno(con, "calls", path, non_reference_only = TRUE)
DBI::dbGetQuery(con, "SELECT record_index, len(calls) AS n FROM calls ORDER BY record_index")
DBI::dbDisconnect(con, shutdown = TRUE)

Fetch Exact Alleles at Requested Sites

Description

Fetch full-site genotypes for requested (chromosome, position, REF, ALT) alleles from an indexed VCF/BCF. Positions become 0-based half-open intervals for one native indexed scan through 'read_geno(regions := ...)', then exact alleles are restored with an equality join. Everything runs on 'con': the request snapshot and the plan variable are session-local, are removed on exit even after an error, and no private connection is opened.

Usage

rduckhts_geno_sites(
  con,
  path,
  sites,
  source_build,
  table_name = NULL,
  samples = NULL,
  format_fields = NULL,
  index_path = NULL,
  decode_error_policy = "error",
  overwrite = FALSE
)

Arguments

Details

Coordinates in 'sites' are 1-based VCF positions. Requested alleles are compared verbatim: nothing is trimmed, left-aligned, strand-flipped, split or lifted over, and 'source_build' is only a caller-asserted label that every request must carry (contig names cannot prove assembly identity). The 'match_status' labels are position-level diagnostics, not claims about normalization equivalence.

The result holds every request column followed by:

'record_index'

Scan-local record ordinal (not a persistent locator).

'full_alt'

The record's complete 'ALT' list in file order.

'alt_index'

1-based ordinal of the matching ALT slot.

'calls'

The full-site 'calls' list of [rduckhts_geno()], including 'format_fields' members such as 'GP', 'DS' or 'HS'.

'match_status'

‘"matched"': chrom, pos, REF and ALT equal a record’s ALT slot (one row per matching record and ALT slot, so physical duplicate records and duplicated ALT slots each return a row). '"allele_not_at_site"': a record at chrom and pos has an equal REF but none of its ALTs equals 'alt'. '"ref_mismatch"': records exist at chrom and pos but none has an equal REF. '"absent"': no record starts at chrom and pos (an unknown contig is absent).

Unmatched requests have one row with 'NULL' record and call fields; no reference genotype is fabricated. Reference and missing calls at a matched site are returned because 'non_reference_only' is fixed to 'FALSE'. Rows carry no implicit order: use 'ORDER BY request_id, record_index, alt_index'.

Value

A data frame, or invisibly 'table_name' when it is given.

See Also

[rduckhts_geno()]

Examples

con <- rduckhts_connect()
path <- system.file("extdata", "geno_sites.bcf", package = "Rduckhts")
DBI::dbWriteTable(con, "requests", data.frame(
  request_id = 1:2, build = "GRCh38", chrom = "chr1", pos = c(100L, 100L),
  ref = "A", alt = c("G", "T")))
rduckhts_geno_sites(con, path, "requests", "GRCh38", format_fields = c("GP", "DS"))
DBI::dbDisconnect(con, shutdown = TRUE)

Create GFF3 Table

Description

Creates a DuckDB table from GFF3 files using the DuckHTS extension.

Usage

rduckhts_gff(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  scan_mode = NULL,
  attributes_map = FALSE,
  attributes_list = FALSE,
  attributes_pairs = FALSE,
  strict = FALSE,
  overwrite = FALSE,
  attributes = NULL
)

Arguments

Value

Invisible TRUE on success

Convert GFF3 reader output to Parquet with DuckHTS metadata

Description

Thin DBI wrapper around extension macro 'duckhts_gff_convert_parquet_sql(...)'.

Usage

rduckhts_gff_convert_parquet(
  con,
  path,
  output,
  columns = NULL,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  attributes_map = FALSE,
  attributes_list = FALSE,
  attributes_pairs = FALSE,
  strict = FALSE,
  where = NULL,
  compression = "zstd",
  row_group_size = 100000L,
  partition_by = NULL,
  include_metadata = TRUE,
  header_text = NULL,
  metadata = NULL,
  metadata_json_file = NULL,
  write_format_version = "1",
  overwrite = FALSE
)

Arguments

Value

Invisibly returns 'output'.

Read multiple GFF files into a DuckDB table

Description

Read and combine multiple GFF3 files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_gff_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  scan_mode = NULL,
  attributes_map = FALSE,
  attributes_list = FALSE,
  attributes_pairs = FALSE,
  strict = FALSE,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

Create GTF Table

Description

Creates a DuckDB table from GTF files using the DuckHTS extension.

Usage

rduckhts_gtf(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  scan_mode = NULL,
  attributes_map = FALSE,
  attributes_list = FALSE,
  attributes_pairs = FALSE,
  overwrite = FALSE,
  attributes = NULL
)

Arguments

Value

Invisible TRUE on success

Read multiple GTF files into a DuckDB table

Description

Read and combine multiple GTF files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_gtf_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  scan_mode = NULL,
  attributes_map = FALSE,
  attributes_list = FALSE,
  attributes_pairs = FALSE,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

Read HTS Header Metadata

Description

Reads file header records from HTS-supported formats using the DuckHTS extension.

Usage

rduckhts_hts_header(con, path, format = NULL, mode = NULL)

Arguments

Value

A data frame with parsed header metadata.

Read HTS Index Metadata

Description

Reads index metadata from HTS-supported index files via DuckHTS.

Usage

rduckhts_hts_index(con, path, format = NULL, index_path = NULL)

Arguments

Value

A data frame with index metadata.

Read Raw HTS Index Blob

Description

Returns raw index metadata blob data for a file index.

Usage

rduckhts_hts_index_raw(con, path, format = NULL, index_path = NULL)

Arguments

Value

A data frame with raw index blob metadata.

Read HTS Index Spans

Description

Returns index span-oriented metadata for planning range workloads.

Usage

rduckhts_hts_index_spans(con, path, format = NULL, index_path = NULL)

Arguments

Value

A data frame with span-oriented index metadata.

Get the Installed htslib Linking Contract

Description

Resolves headers, an exact shared or static library, linker flags, enabled features, and build identity from the installed Rduckhts package. With validation enabled, the receipt and public headers are also compared with the htslib version reported by the loaded DuckHTS extension.

Usage

rduckhts_htslib_config(link = NULL, validate = TRUE)

Arguments

Details

The shared contract is currently available on native Unix builds. MinGW and browser-wasm builds expose the static contract unless a shared htslib was explicitly built. Static consumers must review 'static_license_note'.

Value

An object of class 'rduckhts_htslib_config'. Its 'cppflags' and 'ldflags' elements can be consumed by a downstream package configure script; 'duckdb_platform' records the extension footer platform and the remaining fields form the versioned build receipt.

Examples

## Not run: 
config <- rduckhts_htslib_config()
config$cppflags
config$ldflags
config$features

## End(Not run)

Inspect the Loaded htslib Build

Description

Returns the version and build features reported by the htslib library that is actually loaded with DuckHTS.

Usage

rduckhts_htslib_info(con)

Arguments

Value

A one-row data frame with 'version', 'feature_bits', and 'feature_string' columns.

Examples

## Not run: 
con <- rduckhts_connect()
rduckhts_htslib_info(con)
DBI::dbDisconnect(con, shutdown = TRUE)

## End(Not run)

Return the Loaded htslib Version

Description

Return the Loaded htslib Version

Usage

rduckhts_htslib_version(con)

Arguments

Value

The runtime htslib semantic version as a character scalar.

Install DuckHTS macros on a connection

Description

Executes the ordered definitions exported by the loaded DuckHTS extension as connection-local TEMP macros. Use this on DBI or pool connections that did not come from rduckhts_connect(). The operation is idempotent and shadows same-named persistent macros without deleting them. Inspect old persistent macros in duckdb_functions() using the file's database_name; drop them explicitly if no longer wanted.

Usage

rduckhts_install_macros(con)

Arguments

Details

This function does not begin or commit a transaction. Inside a caller transaction its DDL participates in that transaction: rollback removes newly installed TEMP macros, and the caller can install them again after rollback. Call it separately for each connection, including pool connections.

Value

Invisibly returns the definitions SHA-256 digest.

Lift Over Variant Coordinates Against a Query

Description

Applies the DuckHTS 'duckhts_liftover(...)' table macro to rows from a SQL query or table expression with chromosome and position columns, plus optional reference and alternate alleles.

Usage

rduckhts_liftover(
  con,
  query,
  chain_path,
  dst_fasta_ref,
  chrom_col = "chrom",
  pos_col = "pos",
  ref_col = NULL,
  alt_col = NULL,
  src_fasta_ref = NULL,
  max_snp_gap = 1,
  max_indel_inc = 250,
  lift_mt = FALSE,
  end_pos_col = NULL,
  no_left_align = FALSE
)

Arguments

Value

A data frame with source columns, lifted coordinates/alleles, and warnings.

Load DuckHTS Extension

Description

Loads the DuckHTS extension into an existing DuckDB connection. This must be called before using HTS reader functions on a connection not created by rduckhts_connect().

Usage

rduckhts_load(con, extension_path = NULL)

Arguments

Details

DuckDB 1.5.0 or newer is required. SQL registration errors retain the underlying DuckDB diagnostic in the initialization error.

The connection must permit unsigned extension loading. With current versions of the duckdb R package, its driver must also permit extension loading. Prefer rduckhts_connect() when Rduckhts owns the connection. File-backed connections receive TEMP macros on this connection after LOAD. Use rduckhts_install_macros() on additional connections.

Value

TRUE if the extension was loaded successfully.

Examples

con <- rduckhts_connect()
DBI::dbDisconnect(con, shutdown = TRUE)

Native mosdepth-Compatible Coverage Outputs

Description

Writes native mosdepth-compatible coverage outputs for indexed BAM or CRAM input.

Usage

rduckhts_mosdepth(
  con,
  prefix,
  path,
  chrom = NULL,
  by = NULL,
  fasta = NULL,
  read_groups = NULL,
  no_per_base = FALSE,
  threads = 2,
  processing_threads = 2,
  flag = 1796,
  include_flag = 0,
  fast_mode = FALSE,
  fragment_mode = FALSE,
  use_median = FALSE,
  mapq = 0,
  min_frag_len = -1,
  max_frag_len = -1,
  precision_digits = 2,
  quantize = NULL,
  thresholds = NULL,
  index_path = NULL,
  overwrite = FALSE
)

Arguments

Value

A data frame describing the written output paths

Munge Summary Statistics Rows

Description

Applies the DuckHTS 'duckhts_munge(...)' table macro to rows from a SQL query or table expression, using either an upstream-style preset, a named column map, or a two-column mapping file. When no mapping mode is provided, the bundled 'colheaders.tsv' alias file is used by default.

Usage

rduckhts_munge(
  con,
  query,
  fasta_ref = NULL,
  preset = NULL,
  column_map = NULL,
  column_map_file = NULL,
  iffy_tag = "IFFY",
  mismatch_tag = "REF_MISMATCH",
  ns = NULL,
  nc = NULL,
  ne = NULL
)

Arguments

Value

A data frame with normalized GWAS-VCF-style variant/effect columns.

Create BAM Pileup Table

Description

Creates a DuckDB table from a region-scoped BAM pileup using the DuckHTS extension. This is a compact base/quality pileup view backed by htslib's pileup engine; it is not samtools mpileup text parity.

Usage

rduckhts_pileup(
  con,
  table_name,
  path,
  region,
  index_path = NULL,
  min_mapq = 0,
  flag_mask = 1796,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success

Find Runs of Homozygosity

Description

Detect runs of homozygosity in a VCF or BCF with the two-state hidden Markov model of 'bcftools roh': an autozygous state and a Hardy-Weinberg state, decoded with the Viterbi algorithm. Segments, marker counts and start and end positions equal 'bcftools roh'; quality is the mean forward-backward phred score, equal to the one decimal 'bcftools' prints. The model is a port of 'vcfroh.c' and 'HMM.c' (MIT, Genome Research Ltd), run by the native 'duckhts_roh_segments()' kernel over one list per sample and chromosome.

Usage

rduckhts_roh(
  con,
  path,
  af_tag = NULL,
  af_table = NULL,
  genetic_map = NULL,
  hw_to_az = 6.7e-08,
  az_to_hw = 5e-09,
  gt_error = NULL,
  rec_rate = NULL,
  samples = NULL,
  table_name = NULL,
  overwrite = FALSE,
  reference_table = NULL,
  proportions_table = NULL,
  af_clamp = 0.001,
  max_sites = 2e+07,
  max_site_bytes = 2^32,
  site_rule = "matching",
  af_order = NULL
)

Arguments

Details

Allele frequencies come from an INFO tag ('af_tag', like '–AF-tag', which must be declared 'Type=Float,Number=A') or from a relation ('af_table', like '–AF-file') with columns 'chrom', 'pos', 'ref', 'alt' and 'af', matched to each record by chromosome, position, REF and its ALT alleles joined by commas. Sites without a usable frequency (absent, missing or exactly 0) are skipped, as are records with more than one ALT or none. With 'af_table', 'site_rule = "bcftools"' chooses the record of a repeated position as 'bcftools roh –AF-file' does (see the argument).

Alternatively, 'reference_table' and 'proportions_table' select the ancestry-tuned path. The reference is long-format with columns 'chromosome', 'position', 'allele_a', 'allele_b', 'group_id' and 'frequency' (frequency of 'allele_b'); proportions have 'sample_id', 'group_id' and 'proportion'. Their group sets must match exactly, and each VCF sample must have a proportions row for every group. Per-site AF is the sum of each group's frequency weighted by the sample's proportion, with the proportions divided by their sum (which must be above 0 and at most 1, so 'sum_to_one = FALSE' results and rounded proportions are accepted). REF='allele_a', ALT='allele_b' uses that AF; reversed alleles use '1 - AF'. Reference alleles must be on the forward strand of the VCF's assembly, as in a FASTA-anchored panel; no strand flip is attempted, so palindromic (A/T, C/G) sites are oriented by REF like any other, and a record whose alleles match neither order is not used and does not claim a repeated position. Contig names are matched with 'duckhts_contig_key()' once per distinct name (one leading 'chr' removed, M/MT written as MT, X and Y uppercased), so 'chr1' and '1' match, as do 'chrX' and 'X'; accessions, patches and numeric sex chromosomes are not mapped. 'af_clamp' limits nonzero-clamp frequencies to '[af_clamp, 1-af_clamp]'; the default keeps zero population frequencies from being treated as impossible, and zero disables clamping. Only called sites are frequency-weighted. Supply exactly one frequency source: 'af_tag', 'af_table', or both ancestry relations.

Genotype evidence is FORMAT/PL unless 'gt_error' is given, which uses the diploid GT calls with that phred error ('-G') and needs no PL. With no recombination map, transitions are the per-base-pair probabilities 'hw_to_az' and 'az_to_hw' compounded over the physical distance between sites, exactly as 'bcftools roh' does without '-m' or '-M'. 'rec_rate' is a constant rate per base pair ('-M'). 'genetic_map' names a relation with columns 'chrom', 'pos' and 'cm' (cumulative centimorgans, the IMPUTE2 format), interpolated as '-m' does; chromosomes it does not cover are skipped.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given) with one row per run: 'sample', 'chrom', 'start' and 'end' (one-based, inclusive, the first and last marker), 'length' ('end - start + 1'), 'n_markers' and 'quality'. Rows are unordered.

Examples

con <- rduckhts_connect()
path <- system.file("extdata", "roh_fixture.vcf.gz", package = "Rduckhts")
roh <- rduckhts_roh(con, path, af_tag = "AF")
roh[order(roh$sample, roh$chrom, roh$start), ]
DBI::dbDisconnect(con, shutdown = TRUE)

Find Runs of Homozygosity from Read Counts

Description

Detect runs of homozygosity from allele read counts, for example counts at panel sites from a BAM or CRAM, with the two-state model and transitions of [rduckhts_roh()]. Genotype likelihoods come from a binomial read model with a per-read sequencing error and an optional contamination fraction instead of FORMAT/PL or GT. This emission is a DuckHTS extension; 'bcftools roh' has no read-count mode.

Usage

rduckhts_roh_counts(
  con,
  counts_table,
  seq_error = 0.001,
  contamination = 0,
  genetic_map = NULL,
  hw_to_az = 6.7e-08,
  az_to_hw = 5e-09,
  rec_rate = NULL,
  table_name = NULL,
  overwrite = FALSE,
  max_sites = 2e+07,
  max_site_bytes = 2^32,
  fit_table = NULL,
  allele_balance = 0.5,
  artefact_weight = 0
)

Arguments

Details

'counts_table' names a table or view with columns 'sample_id', 'chrom', 'pos' (one-based), 'ref_count', 'alt_count' and 'af', one row per sample and site. 'alt_count' counts reads that show the allele whose population frequency is 'af'; 'ref_count' counts reads that show the other allele. The frequency can come from any source joined in SQL beforehand: an INFO tag, a frequency table or ancestry-weighted frequencies. Somalier site counts map as 'ref_count = a', 'alt_count = b', with 'af' the frequency of 'allele_b'.

Each read shows the counted allele with probability '(1 - contamination) * q + contamination * c', where 'q' is 'seq_error', 'allele_balance' or '1 - seq_error' for zero, one or two copies, and 'c' is the mean of 'q' over the genotype of a contaminating individual of the same population (Hardy-Weinberg at ‘af'). With chance 'artefact_weight' a site’s counted reads are instead uniform on '0..depth', whatever its genotype. The four numbers mean what [rduckhts_count_error_fit()] reports under these names. Without the contamination term, contaminant reads at homozygous sites look like heterozygous evidence, which shortens runs and lowers FROH.

Sites with no reads, a NULL count, or an 'af' that is NULL, NaN or 0 are skipped; a repeated position keeps one row. Negative counts are an error. The model assumes biallelic sites, independent reads (no duplicate or strand modelling) and one error rate for all reads. Reference-biased mapping makes the other allele of a heterozygote rarer; 'allele_balance' takes one value per sample, while the bias is stronger where flanking heterozygosity is high, so such heterozygotes can still look homozygous and lengthen runs in divergent regions. The contamination fraction is supplied by the caller, for example from [rduckhts_count_error_fit()]. A contaminant from a different population than 'af' describes is approximated by 'af'.

'seq_error' is an error of one read, and the reads of a site multiply. With 'artefact_weight = 0', a site with balanced reads excludes both homozygous genotypes whatever 'seq_error' is: the decode then corresponds to the genotype decode with no tolerated genotype error, and reports fewer bases in runs than 'rduckhts_roh(gt_error = 30)' on the same sample. With 'artefact_weight = w > 0', one site at depth 'd' changes the odds between the two states by a factor of at most '1 + (d + 1) * (1 - w) / w', as 'gt_error' caps one called genotype at '10^(gt_error / 10)'.

The runs depend on 'artefact_weight' more than on the other fitted numbers. With 'fit_table', each sample is therefore decoded three times from one collection of its sites: at its 'artefact_weight' 'w', and at 'w - artefact_weight_sd' and 'w + artefact_weight_sd', each kept in '[0, 0.05]'. Each run of the decode at 'w' gets 'stable_fraction': the share of its bases that are also in a run of both other decodes. A run with 'stable_fraction' 1 does not depend on where 'w' lies within one standard error. Filter on it to keep the runs, or the run bases, that the interval of the fit does not move; a low value marks a run whose call hinges on the artefact weight.

'gained_bases' guards the other direction: the bases that the decode at 'w + artefact_weight_sd' calls next to the run but outside every run at 'w'. A run at 'w + sd' meets a run at 'w' when they share a base or abut. A base of such a run before the first run at 'w' it meets counts for that run; any other counts for the run at 'w' before it, so a gap that the higher weight closes between two runs counts once, for the run before the gap. A run at 'w + sd' that meets no run at 'w' has no row, so its bases are not reported. A high 'gained_bases' marks a run that a higher weight would extend or join to its neighbour.

A NULL 'artefact_weight_sd' (fewer than three blocks in the fit) counts as 0 and gives ‘stable_fraction' 1 and 'gained_bases' 0, so check the fit’s 'blocks' too. The decode then takes about three times as long.

Value

A data frame (or invisible 'TRUE' when 'table_name' is given) with the columns of [rduckhts_roh()]: one row per run, unordered. With 'fit_table' it also has 'stable_fraction' and 'gained_bases' (see Details).

Examples

con <- rduckhts_connect()
DBI::dbExecute(con, paste(
  "CREATE TEMP TABLE counts AS SELECT 'S1' AS sample_id, '1' AS chrom,",
  "i * 25000 AS pos, 0.5 AS af,",
  "CASE WHEN i BETWEEN 40 AND 120 OR i % 2 = 0 THEN 30 ELSE 15 END AS ref_count,",
  "CASE WHEN i BETWEEN 40 AND 120 OR i % 2 = 0 THEN 0 ELSE 15 END AS alt_count",
  "FROM range(1, 161) AS t(i)"
))
rduckhts_roh_counts(con, "counts")
DBI::dbDisconnect(con, shutdown = TRUE)

samtools idxstats-Compatible Alignment Summary

Description

Writes samtools idxstats-compatible alignment summary output for BAM, CRAM, or SAM input.

Usage

rduckhts_samtools_idxstats(
  con,
  path,
  output = NULL,
  index_path = NULL,
  threads = 0,
  overwrite = FALSE
)

Arguments

Value

A data frame with 'success', 'path', 'output_path', 'used_index_fast_path', and 'error_message'

Compute Polygenic Scores

Description

Calls the DuckHTS 'bcftools_score(...)' table function to compute sample-level polygenic scores from one genotype VCF/BCF file and one or more summary-statistics files.

Usage

rduckhts_score(
  con,
  bcf_path,
  summary_path = NULL,
  use = NULL,
  columns = "PLINK",
  columns_file = NULL,
  q_score_thr = NULL,
  summaries_list_file = NULL,
  log_path = NULL,
  use_variant_id = FALSE,
  counts = FALSE,
  samples = NULL,
  force_samples = FALSE,
  regions = NULL,
  regions_file = NULL,
  regions_overlap = 1,
  targets = NULL,
  targets_file = NULL,
  targets_overlap = 0,
  apply_filters = NULL,
  include = NULL,
  exclude = NULL
)

Arguments

Value

A data frame with one row per sample and score/count columns.

DuckHTS SIMD backend diagnostics

Description

Inspect or explicitly select the SIMD dispatch policy used by bundled DuckHTS byte-oriented helper kernels such as seq_gc_content(...).

Usage

rduckhts_simd_backend(con)

rduckhts_simd_requested_backend(con)

rduckhts_simd_backend_compiled(con, backend)

rduckhts_simd_backend_cpu_supported(con, backend)

rduckhts_simd_backend_available(con, backend)

rduckhts_simd_info(con)

rduckhts_simd_kernel_info(con)

rduckhts_simd_set_backend(con, backend = "auto")

Arguments

Value

rduckhts_simd_backend(), rduckhts_simd_requested_backend(), and rduckhts_simd_set_backend() return a character scalar. rduckhts_simd_backend_compiled(), rduckhts_simd_backend_cpu_supported(), and rduckhts_simd_backend_available() return logical scalars. rduckhts_simd_info() returns the extension-owned backend inventory table with one row per known backend; availability means compiled and CPU/runtime supported; SQL-level backend changes are process-wide and use the one-row duckhts_simd_set_backend(...) table function; the selectable column reports whether the backend has a selectable implementation path. Explicit selection still requires available = TRUE. rduckhts_simd_kernel_info() returns one row per logical SIMD kernel and is the authoritative diagnostic for mixed per-kernel auto-dispatch.

Examples

## Not run: 
con <- rduckhts_connect()
rduckhts_simd_info(con)
rduckhts_simd_kernel_info(con)
rduckhts_simd_backend_available(con, "scalar")
rduckhts_simd_set_backend(con, "scalar")
rduckhts_simd_set_backend(con, "auto")
DBI::dbDisconnect(con, shutdown = TRUE)

## End(Not run)

Extract Panel-Aligned Counts from BAM or CRAM

Description

Count observed A, B, and other query bases at every site in a typed panel. The panel is prepared once, then each scan worker owns one indexed multi-region scan and independent BAM/CRAM, index, reference, pileup, and overlap state. 'worker_count' bounds the number of DuckDB-scheduled panel shards; the connection's thread setting bounds how many can run concurrently. 'decompression_threads' separately controls htslib decompression workers per source handle. Valid uncovered sites are measured zero depth; reference or alignment-header mismatches remain rows with NULL counts and a named status. The native function reads the panel from a local Parquet file on a private in-memory instance. A panel table or view is therefore written to a scratch Parquet file for the call, which also makes caller-local temporary relations and uncommitted rows visible.

Usage

rduckhts_somalier_bam_counts(
  con,
  source_path,
  sample_id,
  reference_path,
  panel_table = NULL,
  panel_parquet = NULL,
  index_path = NULL,
  reference_index_path = NULL,
  min_mapq = 1,
  min_baseq = 0,
  require_flags = 0,
  exclude_flags = 1796,
  overlap_policy = c("hileup_v0.1.0", "none"),
  decompression_threads = 0,
  worker_count = 1,
  max_depth = 1e+05,
  max_overlap_qnames = 1e+05,
  max_sites = 1e+06,
  max_region_bytes = 67108864,
  remote_block_bytes = 1048576,
  remote_cache_bytes = 67108864,
  reference_cache_bytes = 67108864,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Value

A data frame ordered by 'site_index' if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Estimate Per-Sample Contamination with CHARR

Description

Apply the Somalier-derived CHARR estimator to measured A/B/other count evidence aligned to an ordered panel and population-B allele frequencies. ‘frequency_table' contains the panel’s six identity columns plus 'population_b_af'; it must cover every panel site exactly once. Evidence and frequency identities, coordinates, and A/B orientation are checked against the panel before their digests are derived. A result with no usable homozygous-like evidence has status 'no_evidence' and a NULL estimate.

Usage

rduckhts_somalier_charr(
  con,
  evidence_table = NULL,
  evidence_parquet = NULL,
  panel_table = NULL,
  panel_parquet = NULL,
  frequency_table = NULL,
  frequency_parquet = NULL,
  table_name = NULL,
  sample_ids = NULL,
  min_depth = 15,
  max_depth = 1e+06,
  hom_minor_rate = 0.12,
  hom_tail_alpha = 0.002,
  max_threshold_work = 1.6e+07,
  max_sites = 1e+06,
  overwrite = FALSE
)

Arguments

Details

Each of the evidence, panel, and frequency inputs is supplied as exactly one named table/view or ordinary Parquet path. Optional sample selection is exact: every requested ID must occur. Results contain the panel and frequency digests, usable-site denominators, numerical status, and all filter settings.

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Select a Somalier Panel from Population Variants

Description

Select typed population variants on the caller's DuckDB connection. A table or view may be populated by 'read_bcf()'; Parquet and VCF/BCF paths are also supported. Input columns follow 'read_bcf()' ('CHROM', 'POS', 'REF', 'ALT', 'FILTER', and declared 'INFO_*' fields). Interval relations have 'chrom', 'start', 'end' in zero-based half-open BED coordinates, as returned by 'read_bed()'. A gnotate exclusion relation has 'chrom', 'pos', 'ref', 'alt' in one-based VCF coordinates.

Usage

rduckhts_somalier_find_sites(
  con,
  source_table = NULL,
  source_parquet = NULL,
  source_vcf = NULL,
  assembly,
  min_af = 0.15,
  min_an = 115000,
  af_field = "AF",
  an_field = "AN",
  snp_dist = 10000,
  target_af = 0.48,
  include_table = NULL,
  exclude_table = NULL,
  gnotate_exclude_table = NULL,
  compatibility_mode = c("somalier_v0.3.4", "generic"),
  max_autosomal = 65535,
  max_x = 10001,
  max_y = 5001,
  table_name = NULL,
  overwrite = FALSE,
  sex_spacing = c("upstream", "enforced"),
  tie_order = c("input", "lexical")
)

Arguments

Details

'somalier_v0.3.4' names the pinned human-contig rules. Upstream skips autosomal REF=C and restricts X to the inclusive zero-based range 2781479..154931044. By default X/Y selections do not enter the spacing state, as in upstream v0.3.4. 'sex_spacing = "enforced"' applies the X/Y distances to every selection. Upstream's stable AF-score sort retains input order on ties; this selector uses scan order by default, or lexical position and allele order with 'tie_order = "lexical"'. For a table or view, scan order is only reliable when its query explicitly orders the records. Only uppercase A/C/G/T biallelic SNVs enter the canonical panel. 'generic' does not interpret contig names as human sex chromosomes. The result can include X/Y, which 'duckhts_somalier_import_sites()' keeps as sex-chromosome sites for 'rduckhts_somalier_sex()'. Gnotate zip input is not supported: use an allele-specific exclusion relation. 'assembly' labels the input coordinates; this function does not run the liftover steps in Somalier's 'scripts/find_sites.sh'.

Value

A canonical panel relation with selection score, source AN and spacing-distance diagnostics.

Import an Already Selected Somalier Sites VCF or BCF

Description

Convert an existing Somalier-compatible sites file into the canonical typed panel and population-frequency relation used by DuckHTS extraction, relatedness, contamination, and sex functions. REF and ALT are oriented into lexical A/B order and alternate-allele frequency is flipped with the alleles, so 'population_b_af' always describes 'allele_b'. Records on the exact Somalier v0.3.4 X/Y aliases are kept as sex-chromosome sites, numbered after every autosomal site; relatedness and contamination functions ignore them and 'rduckhts_somalier_sex()' reads them. No pseudo-autosomal filtering happens on import or extraction, as in Somalier: select X sites outside the PAR when choosing the sites file. At least one autosomal record is required.

Usage

rduckhts_somalier_import_sites(
  con,
  path,
  assembly,
  max_sites = 1e+06,
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Details

This function does not select sites from a population VCF. Somalier's 'find-sites' algorithm has separate AF/AN, QC, interval-exclusion, and spacing semantics and is not implied by importing its output.

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Estimate Directional Contamination Against Matched Anchors

Description

Estimate one contamination fraction for each explicitly ordered receiver/anchor pair. The pair relation must contain distinct nonempty 'receiver_id' and 'anchor_id' columns. Every listed sample must have one count tuple for every ordered panel site. Both samples and population-B frequencies use the panel's A/B orientation; the anchor supplies the expected uncontaminated receiver genotype and is not interpreted as the contaminating donor. Reversing a pair is therefore a different analysis.

Usage

rduckhts_somalier_matched_contamination(
  con,
  evidence_table = NULL,
  evidence_parquet = NULL,
  panel_table = NULL,
  panel_parquet = NULL,
  frequency_table = NULL,
  frequency_parquet = NULL,
  pairs_table = NULL,
  pairs_parquet = NULL,
  table_name = NULL,
  min_depth = 15,
  max_depth = 1e+06,
  hom_minor_rate = 0.05,
  hom_tail_alpha = 0.001,
  error_rate = 0.002,
  min_probability = 1e-10,
  min_prior_frequency = 1e-06,
  alpha_min = 0,
  alpha_max = 1,
  grid_step = 0.01,
  refine_tolerance = 1e-10,
  max_evaluations = 4096,
  max_threshold_work = 1.6e+07,
  max_sites = 1e+06,
  overwrite = FALSE
)

Arguments

Details

Sample identifiers and assembly strings are limited to 1,024 bytes, including when prepared profiles are persisted and supplied directly.

The returned score omits alpha-independent binomial coefficients. It compares alpha candidates for the same observed receiver/anchor pair and is not an absolute likelihood comparable between pairs. No usable evidence yields status 'no_evidence' with NULL alpha and relative log-likelihood.

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Compare Somalier-Derived Sample Sketches

Description

Compute the named relatedness and concordance statistics for either every distinct pair in one sketch relation or the ordered pairs in 'pairs_table'. The sketch relation must contain one non-NULL 'sketch' struct per sample. Selected pairs are an ordinary relation with 'sample_a' and 'sample_b' columns. Missing or duplicate sample and pair identities error instead of silently dropping or multiplying requested comparisons. The native kernel checks assembly, ordered-panel digest, classification settings, mask shape, and mask contents. For all pairs it checks the contents of each sketch once and then compares the checked sketches; for selected pairs it checks both sketches of each comparison. No SQL row-order guarantee is implied.

Usage

rduckhts_somalier_relatedness(
  con,
  sketches_table = NULL,
  sketches_parquet = NULL,
  pairs_table = NULL,
  table_name = NULL,
  max_sites = 1e+06,
  max_pairs = 5e+07,
  overwrite = FALSE
)

Arguments

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Infer X/Y Dosage Sex Evidence

Description

Report, per sample, the X and Y depth relative to autosomal depth, the heterozygous and homozygous-alternate counts of usable X sites, and a call of '"XX"', '"XY"' or '"ambiguous"' following Somalier v0.3.4 ('relate.nim'). The counts relation comes from 'rduckhts_somalier_vcf_counts()' or 'rduckhts_somalier_bam_counts()' over a panel that includes X/Y sites (see 'rduckhts_somalier_import_sites()'); every selected sample needs one row for every panel site.

Usage

rduckhts_somalier_sex(
  con,
  counts_table = NULL,
  counts_parquet = NULL,
  panel_table = NULL,
  panel_parquet = NULL,
  table_name = NULL,
  sample_ids = NULL,
  min_depth = 7,
  min_sex_site_depth = 7,
  min_het_balance = 0.3,
  hom_balance_cutoff = 0.01,
  min_usable_x_sites = 11,
  xy_max_het_ratio = 0.05,
  xx_min_het_ratio = 0.4,
  y_signal_min = 0.4,
  y_gate = c("sample", "cohort"),
  overwrite = FALSE
)

Arguments

Details

Depth ratios are '2 * mean A+B depth over usable sites / autosomal mean depth', with the autosomal mean taken over sites passing the relatedness allele-balance gate. The call is '"XY"' when at least 'min_usable_x_sites' X sites are usable and their het/hom-alt ratio is below 'xy_max_het_ratio', '"XX"' when it exceeds 'xx_min_het_ratio', and otherwise '"ambiguous"'; an XX dosage with a present Y signal is also '"ambiguous"'. Metrics are 'NA' without usable sites or autosomal depth, with a 'status' naming why. The call is evidence for review, not a diagnosis or determination of sex.

Somalier applies its Y check only when the cohort has Y depth, which makes a sample's call depend on the rest of the batch. 'y_gate = "sample"' (the default) applies it whenever the sample has a Y signal; '"cohort"' reproduces Somalier's gate, computed among the selected samples, and 'cohort_has_y' reports it.

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Prepare Somalier-Derived Sample Sketches

Description

Build packed, panel-verified relatedness sketches from measured A/B/other count evidence. The panel must contain 'assembly', zero-based 'site_index', 'region', one-based 'position', and uppercase single-base 'allele_a' and 'allele_b'. Evidence must contain the same site identity columns plus 'sample_id' and nullable count columns 'a', 'b', and 'other'. All three counts are NULL for unavailable evidence; three measured zeros are not unavailable. The native SQL preparation checks every evidence site's geometry and A/B orientation against the ordered panel before computing its digest. Panel alleles must be distinct uppercase single-base A/C/G/T with lexical A < B. Rows on the exact Somalier v0.3.4 X/Y aliases are sex-chromosome sites that must follow every autosomal site in 'site_index'; they and their evidence are ignored here, so results equal those on the autosomal-only panel. Other contig aliases cannot be classified biologically from the region string. The three-state calculation assumes diploid sites; count evidence alone does not prove sample ploidy.

Usage

rduckhts_somalier_sketches(
  con,
  evidence_table = NULL,
  evidence_parquet = NULL,
  panel_table = NULL,
  panel_parquet = NULL,
  table_name = NULL,
  sample_ids = NULL,
  min_depth = 7,
  min_het_balance = 0.3,
  hom_balance_cutoff = 0.01,
  max_sites = 1e+06,
  overwrite = FALSE
)

Arguments

Details

Supply each source as either a table/view name or an ordinary Parquet path. Parquet inputs are exposed through query-scoped temporary views; no private sketch format or user-supplied panel digest is involved. The result has one 'sketch' struct per selected sample. Its packed words can be persisted with DuckDB's usual Parquet 'COPY' statement.

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Extract Panel-Aligned Counts from VCF or BCF

Description

Produce the complete sample-by-panel count relation consumed by the Somalier-derived relatedness, contamination and sex functions. 'FORMAT/AD' must declare 'Number=R,Type=Integer'; A and B slots are matched by exact REF/ALT identity, and 'other' sums only the remaining declared-allele slots. Missing sites and unavailable AD remain rows with three NULL counts, distinct from measured zero depth. The panel can be any typed table/view or ordinary Parquet file with the canonical six panel identity columns.

Usage

rduckhts_somalier_vcf_counts(
  con,
  path,
  panel_table = NULL,
  panel_parquet = NULL,
  samples = NULL,
  filter_policy = c("pass_or_unapplied", "include_all", "error"),
  table_name = NULL,
  overwrite = FALSE
)

Arguments

Value

A data frame if 'table_name' is 'NULL'; otherwise invisible 'TRUE'.

Create Tabix-Indexed File Table

Description

Creates a DuckDB table from any tabix-indexed file using the DuckHTS extension.

Usage

rduckhts_tabix(
  con,
  table_name,
  path,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  scan_mode = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success

Convert generic tabix reader output to Parquet with DuckHTS metadata

Description

Thin DBI wrapper around extension macro 'duckhts_tabix_convert_parquet_sql(...)'.

Usage

rduckhts_tabix_convert_parquet(
  con,
  path,
  output,
  columns = NULL,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  where = NULL,
  compression = "zstd",
  row_group_size = 100000L,
  partition_by = NULL,
  include_metadata = TRUE,
  header_text = NULL,
  metadata = NULL,
  metadata_json_file = NULL,
  write_format_version = "1",
  overwrite = FALSE
)

Arguments

Value

Invisibly returns 'output'.

Build Tabix Index

Description

Builds a tabix index for a BGZF-compressed text file using the DuckHTS extension.

Usage

rduckhts_tabix_index(
  con,
  path,
  preset = "vcf",
  index_path = NULL,
  min_shift = 0,
  threads = 4,
  seq_col = NULL,
  start_col = NULL,
  end_col = NULL,
  comment_char = NULL,
  skip_lines = NULL
)

Arguments

Value

A data frame with 'success', 'index_path', and 'index_format'

Read multiple tabix-indexed files into a DuckDB table

Description

Read and combine multiple tabix-indexed files via UNION ALL BY NAME, materialising the result as a DuckDB table. Each row includes a filename column identifying its source file.

Usage

rduckhts_tabix_multi(
  con,
  table_name,
  files,
  region = NULL,
  index_path = NULL,
  header = NULL,
  header_names = NULL,
  auto_detect = NULL,
  column_types = NULL,
  scan_mode = NULL,
  .params = NULL,
  overwrite = FALSE
)

Arguments

Value

Invisible TRUE on success.

Setup HTSlib Environment

Description

Sets the 'HTS_PATH' environment variable to point to the bundled htslib plugins directory. This enables remote file access via libcurl plugins (e.g., s3://, gs://, http://) when plugins are available.

Usage

setup_hts_env(plugins_dir = NULL)

Arguments

Details

Call this before the process opens its first HTS file. htslib discovers dynamic plugins on first file access and does not rescan 'HTS_PATH' later.

Value

Invisibly returns the previous value of 'HTS_PATH' (or 'NA' if unset).

Examples

## Not run: 
setup_hts_env()

plugins_path <- tempfile("hts_plugins_")
dir.create(plugins_path)
setup_hts_env(plugins_dir = plugins_path)
unlink(plugins_path, recursive = TRUE)

## End(Not run)