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

Extension Function Reference

Generated from functions.yaml.

Export the ordered DuckHTS macro definitions for connection-local installation.

Signature:

duckhts_macro_definitions()

Returns:

TABLE

LOAD creates macros in a writable in-memory default database. For file-backed or read-only default databases LOAD creates no macros: execute the exported TEMP statements in install_order on each connection that needs them (R: rduckhts_install_macros). The decision concerns the database instance's default database; a connection that USEs another catalog before or after LOAD installs the TEMP statements too. TEMP macros shadow, but do not delete, persistent macros in older files. The definitions_sha256 value identifies the entire ordered definition set.

install_order UINTEGER, name VARCHAR, public BOOLEAN (listed in functions.yaml), sql VARCHAR (CREATE OR REPLACE TEMP MACRO statement), definitions_sha256 VARCHAR.

SELECT install_order, name, sql FROM duckhts_macro_definitions() ORDER BY install_order;

Return the runtime version reported by the htslib library loaded with DuckHTS. Rduckhts uses this value to reject a downstream linking receipt whose source/header version does not match the loaded library.

Signature:

duckhts_htslib_version()

Returns:

VARCHAR

SELECT duckhts_htslib_version();

Return the htslib runtime feature bitfield reported by hts_features(). Use duckhts_htslib_feature_string() for the corresponding build description.

Signature:

duckhts_htslib_features()

Returns:

UINTEGER

SELECT duckhts_htslib_features();

Return htslib's runtime build-feature description, including configured transports, compression libraries, compiler, and build flags. DuckHTS snapshots it once while loading the extension so parallel SQL calls read immutable text.

Signature:

duckhts_htslib_feature_string()

Returns:

VARCHAR

SELECT duckhts_htslib_feature_string();

Return the current DuckHTS SIMD dispatch label. For explicit scalar or concrete backend requests this is the requested policy; for auto it is the single selected backend when all logical kernels resolve to the same backend, or mixed when per-kernel auto-dispatch resolves to multiple backends. Use duckhts_simd_kernel_info() for per-kernel details.

Signature:

duckhts_simd_backend()

Returns:

VARCHAR

SELECT duckhts_simd_backend();

Return the current explicit SIMD backend request, usually auto unless SELECT backend FROM duckhts_simd_set_backend('auto'|'scalar'|backend) was called. The selected per-kernel backend may differ under auto-dispatch across x86, ARM, wasm, and scalar-only builds.

Signature:

duckhts_simd_requested_backend()

Returns:

VARCHAR

SELECT duckhts_simd_requested_backend();

Return whether a concrete DuckHTS SIMD backend was compiled into this build. This is independent of whether the current CPU/runtime supports executing that backend; for example avx512 can be compiled but not CPU-supported on the running host.

Signature:

duckhts_simd_backend_compiled(backend)

Returns:

BOOLEAN

SELECT duckhts_simd_backend_compiled('scalar');

SELECT duckhts_simd_backend_compiled('avx512');

Return whether the current CPU/runtime supports a concrete DuckHTS SIMD backend, independent of whether DuckHTS compiled an implementation for it. Availability is the intersection of compiled and CPU-supported.

Signature:

duckhts_simd_backend_cpu_supported(backend)

Returns:

BOOLEAN

SELECT duckhts_simd_backend_cpu_supported('avx2');

SELECT duckhts_simd_backend_cpu_supported('avx512');

Return whether a concrete SIMD backend is usable in the current process. Availability means the backend is compiled into DuckHTS and supported by the current CPU/runtime. auto is a selection request rather than a concrete backend and is not reported as available here.

Signature:

duckhts_simd_backend_available(backend)

Returns:

BOOLEAN

SELECT duckhts_simd_backend_available('scalar');

SELECT duckhts_simd_backend_available('avx2');

SELECT duckhts_simd_backend_available('avx512');

Report compiled, runtime-supported and selected status for each concrete DuckHTS SIMD backend.

Signature:

duckhts_simd_info()

Returns:

table

Rows include selectable, compiled, CPU-supported, available, selected, requested and dispatch-mode fields. available requires both compiled and CPU/runtime-supported. Explicit selection requires available=TRUE and a selectable implementation path. selected means at least one logical kernel uses that backend. auto is a request, not a concrete backend row.

SELECT * FROM duckhts_simd_info();

SELECT backend FROM duckhts_simd_info() WHERE available;

Return one row per logical DuckHTS SIMD kernel showing the concrete backend selected by the current immutable dispatch table, the selected capability, the requested backend policy, whether scalar was used as a per-kernel fallback, and the dispatch mode. This is the authoritative diagnostic for mixed auto-dispatch when different kernels resolve to different backends.

Signature:

duckhts_simd_kernel_info()

Returns:

table

SELECT * FROM duckhts_simd_kernel_info();

SELECT kernel, selected_backend FROM duckhts_simd_kernel_info();

Explicitly select the DuckHTS SIMD dispatch policy for this process using a one-row table-function call and return the current dispatch label in a backend column. Use auto for per-kernel runtime dispatch or scalar for a portable baseline; unavailable platform-specific requests such as avx512 on non-AVX-512 CPUs raise an error instead of silently falling back.

Signature:

duckhts_simd_set_backend(backend)

Returns:

table(backend VARCHAR)

SELECT backend FROM duckhts_simd_set_backend('scalar');

SELECT backend FROM duckhts_simd_set_backend('auto');

Return whether the currently open DuckDB runtime advertises a logical type with the given name through duckdb_types(). This is a catalog-level runtime probe for feature gating SQL/macros across DuckDB versions.

Signature:

duckhts_duckdb_type_supported(candidate_type_name)

Returns:

BOOLEAN

SELECT duckhts_duckdb_type_supported('VARIANT');

SELECT duckhts_duckdb_type_supported('GEOMETRY');

Return whether the currently open DuckDB runtime advertises the VARIANT logical type. Use this to gate optional SQL that depends on DuckDB VARIANT support.

Signature:

duckhts_duckdb_supports_variant()

Returns:

BOOLEAN

SELECT duckhts_duckdb_supports_variant();

Return whether the currently open DuckDB runtime advertises the GEOMETRY logical type. Use this to gate optional SQL that depends on DuckDB GEOMETRY support.

Signature:

duckhts_duckdb_supports_geometry()

Returns:

BOOLEAN

SELECT duckhts_duckdb_supports_geometry();

Read VCF/BCF with header-typed INFO/FORMAT, typed CSQ/ANN/BCSQ annotations, sample selection, optional tidy sample rows and optional record file positions.

Signature:

read_bcf(path, region := NULL, index_path := NULL, tidy_format := FALSE, additional_csq_column_types := NULL, scan_mode := 'auto', decompression_threads := 0, decode_error_policy := 'null', samples := NULL, regions := NULL, include_file_offset := FALSE)

Returns:

table

INFO/FORMAT Type and Number determine SQL types and scalar/list shape, including nonstandard declarations for named tags; no schema repair is inferred. additional_csq_column_types overrides typed annotation subfields. Missing numeric/string list positions are NULL; numeric vector-end padding terminates lists. Absent fields are NULL; explicit numeric missing lists are [NULL], without inferred allele/genotype cardinality. Tidy rows repeat complete record annotation across output chunks.

FILTER is VARCHAR[] in header order. PASS is [PASS], an unapplied dot is NULL, and named failures retain their identifiers.

decode_error_policy is null, warn or error for header/payload type mismatches and oversized numeric scalars. Multiple scalar elements before vector-end padding count as malformed, including missing elements. FORMAT null/warn withholds that tag for every selected sample on the record.

HTSlib selectors: NULL/'-' keeps all, empty string keeps none, comma-separated names include, leading '^' excludes. Unknown names error. Selection retains header order; read_bcf_samples() supplies original-header indices/names.

Comma-separated indexed regions use a union iterator, emitting each physical record once across overlaps. NULL/empty region means no filter; empty list items and malformed known-contig intervals error. Unknown contigs follow HTSlib's skip policy.

Typed regions: regions := STRUCT(chrom VARCHAR, start BIGINT, "end" BIGINT)[] selects 0-based half-open intervals (VCF POS p is [p-1, p)) with one native multi-region iterator built directly from header/index contig ids, not region strings. NULL or omitted regions keeps the ordinary behaviour; a non-NULL empty list selects zero records, including COUNT(), and never falls back to a scan (list() over no rows is NULL, so wrap it as coalesce(list(...), []::STRUCT(chrom VARCHAR, start BIGINT, "end" BIGINT)[])). region and non-NULL regions are mutually exclusive; regions requires a usable index (.tbi/.csi) and is incompatible with scan_mode := 'sequential'. chrom is a literal contig name, so names such as 'HLA-A01:01' or 'chr1:2' need no braces; a contig absent from the header/index yields no records for that interval, like an unknown contig in region. NULL fields, start < 0 and end <= start error. Intervals are sorted and overlapping or adjacent ones merged per contig without broadening; a physical record is returned once however many intervals overlap it, and physical duplicates stay separate. Records that merely overlap an interval (a deletion spanning it) are returned: filter POS in SQL for exact positions. Caps: 1,000,000 intervals and 128 MiB of native payload, checked before storage; they do not bound the caller's own list aggregation. Prepare the plan on the caller connection with SET VARIABLE p = (SELECT coalesce(list({'chrom': chrom, 'start': pos - 1, 'end': pos}), []::STRUCT(chrom VARCHAR, start BIGINT, "end" BIGINT)[]) FROM (SELECT DISTINCT chrom, pos FROM requests)) and pass regions := getvariable('p'). An unknown variable name is NULL to DuckDB (a full scan), so rduckhts_bcf(regions_var =) and rduckhts_geno(regions_var =) check that the variable exists.

scan_mode='sequential' streams without loading an index. auto streams if no usable index was available at bind; region queries require one. Index-only counts do not validate data contents. Automatic/region plans retain the parsed bind-time index, including remote/explicit index_path: subsequent index removal, replacement or corruption does not change that bound plan. SQL prepared statements can rebind at execution and then resolve the current index; PREPARE does not pin an index across executions. Supply a matching data/index pair whenever a query binds and keep data/header contents unchanged during its execution. This is not a data snapshot or validation of an unrelated index.

include_file_offset := TRUE adds FILE_OFFSET UBIGINT after FILTER: the BGZF virtual position immediately after each record, not its start, as in read_bam. Sequential, indexed, parallel, region and typed-region scans of one unchanged file give a record the same value, and it increases strictly in file order, so ORDER BY FILE_OFFSET restores file order, also among records at one position. Tidy rows of a record share it. BGZF input includes BCF written by bcftools -Ou, which is level-0 BGZF. Uncompressed VCF, plain gzip VCF and a BCF stream without BGZF blocks return NULL. Values are not comparable across files.

decompression_threads controls per-file HTSlib workers, default 0; it is separate from DuckDB scan parallelism.

SELECT CHROM, POS, REF, ALT FROM read_bcf('vcf_file.bcf') LIMIT 5;

SELECT CHROM, POS, ID, FILE_OFFSET FROM read_bcf('cohort.vcf.gz', include_file_offset := true) ORDER BY FILE_OFFSET;

Read one row per VCF/BCF record with typed arbitrary-ploidy GT/PS calls, selected FORMAT fields and optional original VCF genotype text.

Signature:

read_geno(path, region := NULL, index_path := NULL, samples := NULL, non_reference_only := FALSE, scan_mode := 'auto', decompression_threads := 0, decode_error_policy := 'null', format_fields := NULL, raw_gt := FALSE, include_filter := FALSE, regions := NULL)

Returns:

table

Return record_index UBIGINT, CHROM, POS, ID, REF, ALT VARCHAR[] and calls STRUCT(sample_index UINTEGER, alleles INTEGER[], phase_before BOOLEAN[], phase_set BIGINT)[]. Sample indices are zero-based original-header positions; join read_bcf_samples() from the same unchanged file for names. Absent GT has NULL allele/phase lists; missing allele slots remain NULL entries. Phase flags follow HTSlib, including its pre-VCF-4.4 leading-slot convention. Absent PS is NULL; PS must declare Number=1,Type=Integer.

samples uses read_bcf's selectors. non_reference_only removes calls without a called ALT, never records. Zero selected samples yields empty calls on every row. Calls decode only when projected; no normalization, depth or phase inference is performed.

Full scans stream in physical input order regardless of index availability, with one scan worker and separate HTSlib decompression threads. record_index is zero-based and scan-local, not a persistent file locator. Indexed region unions restart at zero in HTSlib iterator order; overlapping regions visit a physical record once while preserving distinct duplicates. Use ORDER BY record_index when order matters. Sequential mode rejects regions and typed regions.

format_fields adds header-typed members in calls.format, such as AD/DP/GQ; NULL/[] preserves default schema. Header Type/Number determine scalar/list shape. Missing elements retain ordinals, vector-end padding terminates lists, absent fields are NULL and explicit numeric missing lists are [NULL], without inferred A/R/G padding. Missing GT does not hide selected values unless the call is filtered. Unknown, empty, NULL or case-insensitively duplicate selections error. Lookup uses exact header spelling: GT/PS cannot be reselected; declared lowercase gt/ps are distinct tags.

include_filter := TRUE adds the same physical record's FILTER as VARCHAR[]. PASS is [PASS], an unapplied dot is NULL, and named failing filters retain header order. FALSE preserves the default output schema.

Header/type/cardinality mismatches follow decode_error_policy. Scalar cardinality excludes vector-end padding, including after sample selection, but counts missing elements. Multiple scalar elements error or make that FORMAT tag NULL for every selected sample under null/warn; list cardinality is not inferred. Physical read/allocation failures always error.

raw_gt := TRUE appends calls.raw_gt VARCHAR after any format struct, preserving exact original VCF separators, leading phase markers and allele spelling. Absent GT is NULL; literal '.' stays text. FALSE preserves the default schema. BCF rejects TRUE even when calls is unprojected. Text follows original-header sample selection and scan ordinals, including region unions; it does not emulate a downstream parser.

regions := STRUCT(chrom VARCHAR, start BIGINT, "end" BIGINT)[] has the same contract as read_bcf: 0-based half-open, NULL vs empty, literal contig names, index required, 1,000,000 interval and 128 MiB caps; see read_bcf. Ordinals restart at zero in HTSlib iterator order.

Exact alleles: regions plus SQL restore requested (chrom, pos, REF, ALT) alleles without changing the reader. Keep requests unchanged, plan I/O with SELECT DISTINCT chrom, pos, read with non_reference_only := false and format_fields such as ['GP','DS','HS'] so hom-ref and missing calls are returned, expand ALT with its ordinal, and join on equality: WITH alleles AS (SELECT record_index, CHROM, POS, REF, ALT AS full_alt, calls, unnest(ALT) AS match_alt, generate_subscripts(ALT, 1) AS alt_index FROM read_geno('f.bcf', regions := getvariable('p'), non_reference_only := false)) SELECT s.*, a.record_index, a.full_alt, a.alt_index, a.calls FROM requests s LEFT JOIN alleles a ON s.chrom = a.CHROM AND s.pos = a.POS AND s.ref = a.REF AND s.alt = a.match_alt. Use equality on the expanded ALT: list_contains(a.ALT, s.alt) in the join condition plans a nested-loop join. Multi-ALT records keep every ALT in file order with the matched 1-based alt_index, and a request matching two ALT slots or two physical records returns a row each. GP, DS, HS and phase are full-site values: no allele splitting, GP remapping, dosage selection, trimming, strand flip or liftover, so a shifted indel representation is not matched. An ALT-less record disappears under unnest, so diagnosing unmatched requests needs a second step; rduckhts_geno_sites() composes the complete recipe on the caller connection, labelling each request matched, allele_not_at_site, ref_mismatch or absent (position-level diagnostics only, never a fabricated genotype). record_index is scan-local; order with ORDER BY request_id, record_index, alt_index. The reader does not deduplicate requests.

SELECT record_index, CHROM, POS, unnest(calls) AS call FROM read_geno('geno_calls.bcf', non_reference_only := TRUE) ORDER BY record_index;

SET VARIABLE p = [{'chrom': 'chr1', 'start': 99, 'end': 100}, {'chrom': 'HLA-A*01:01', 'start': 9, 'end': 10}]; SELECT record_index, CHROM, POS, REF, ALT FROM read_geno('geno_sites.bcf', regions := getvariable('p'), non_reference_only := FALSE, format_fields := ['GP', 'DS', 'HS']);

Read the typed VCF/BCF sample catalog as sample_index UINTEGER and sample_name VARCHAR without reading records. Indices are zero-based positions in the original header, remain stable under selection and join read_geno calls from the same unchanged file. NULL or '-' selects all; an empty string selects none; comma-separated names include samples; '^' excludes them. Names are validated by HTSlib, selected rows retain header order, and unknown names error.

Signature:

read_bcf_samples(path, samples := NULL)

Returns:

table

SELECT * FROM read_bcf_samples('geno_calls.bcf', samples := '^S1') ORDER BY sample_index;

Read SAM/BAM/CRAM alignments with optional typed SAM tags, auxiliary maps and packed sequence, quality or CIGAR output.

Signature:

read_bam(path, standard_tags := FALSE, auxiliary_tags := FALSE, region := NULL, index_path := NULL, reference := NULL, sequence_encoding := 'string', quality_representation := 'string', cigar_representation := 'string', scan_mode := 'auto', decompression_threads := 2)

Returns:

table

AUXILIARY_TAGS leaves out a tag whose key is not valid UTF-8, since a MAP key cannot be NULL; a value that is not valid UTF-8 is NULL.

sequence_encoding='nt16' returns SEQ UTINYINT[]; quality_representation='phred' returns QUAL UTINYINT[]; cigar_representation='binary' returns packed BAM CIGAR UINTEGER[] instead of SAM text.

FILE_OFFSET is the BGZF virtual position immediately after each compressed BAM record, not its start. SAM (including compressed SAM), CRAM and non-BGZF BAM return NULL. ORDER BY FILE_OFFSET orders one unchanged BGZF BAM; arrival order is not guaranteed and offsets are not comparable across files.

scan_mode='sequential' streams rather than using indexed count/parallel paths and rejects region. NULL/empty region means no filter; empty comma-separated items and malformed known-contig intervals error. Unknown contigs follow HTSlib's skip policy. An indexed full-file scan is split by reference across DuckDB threads and reports the index's row total to the planner as an estimate. Both need the optional per-reference statistics of a BAI or CSI: when a reference has alignments and no statistics, the scan is one sequential stream and no estimate is reported. A CRAM scan is split by reference and reports no estimate. In a region list, '.' is the whole file and '*' is the reads without coordinates. htslib locates both from the same optional statistics; for a BAM whose index lacks them, read_bam reads those two items sequentially, so the result does not depend on the index statistics.

decompression_threads controls per-file HTSlib workers, default 2; use 0 to disable. It does not set DuckDB processing parallelism.

A CRAM scan asks htslib only for the fields its projected columns need, and every projected column keeps its value. SEQ and QUAL are decoded only when SEQ, QUAL, CIGAR, MD, NM or AUXILIARY_TAGS is projected, because htslib makes MD and NM from the reference when the file does not store them. Tag data is decoded only when CIGAR, a tag column, AUXILIARY_TAGS, READ_GROUP_ID or SAMPLE_ID is projected; CIGAR needs both because a long CIGAR may be stored as a placeholder with its real value in a CG tag, which htslib restores only from tag data and SEQ. Without SEQ and QUAL, htslib does not load the reference for a slice that covers one reference.

With decompression_threads > 0, a CRAM scan that projects tag data keeps htslib's spare record buffers, which only grow: memory rises with the alignments read, toward (2 x decompression_threads + 1) slices x about 10,000 records x the largest record. A CRAM scan that decodes SEQ and QUAL with decompression_threads > 0 loads the whole reference sequence of each contig it reads. decompression_threads := 0 avoids both.

SELECT QNAME, FLAG, RNAME, POS FROM read_bam('range.bam') LIMIT 5;

Build COPY SQL for read_bcf() output with Parquet metadata, VCF header text and selected columns, filters or partitions.

Signature:

duckhts_bcf_convert_parquet_sql(path, output, columns := []::VARCHAR[], region := NULL, index_path := NULL, tidy_format := FALSE, additional_csq_column_types := NULL, decompression_threads := 0, where_sql := NULL, compression := 'zstd', row_group_size := 100000, partition_by := []::VARCHAR[], include_metadata := TRUE, header_text := NULL, metadata := map([]::VARCHAR[], []::VARCHAR[]), metadata_json_file := NULL, overwrite := FALSE, write_format_version := '1')

Returns:

VARCHAR

Include DuckHTS key/value metadata, preserve or correct header text, and add user metadata with metadata := map(...). metadata_json_file is caller-managed and requires DuckDB's json extension at builder invocation; use maps for offline/CRAN workflows.

The function returns SQL without executing it; the R wrapper executes it through DBI.

SELECT duckhts_bcf_convert_parquet_sql('cohort.vcf.gz', 'cohort.parquet', columns := ['CHROM','POS','REF','ALT'], metadata := map(['project'], ['cohort-a']));

Build COPY SQL for read_bam() output with Parquet metadata, SAM header text and selected columns, filters or partitions.

Signature:

duckhts_bam_convert_parquet_sql(path, output, columns := []::VARCHAR[], 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_sql := NULL, compression := 'zstd', row_group_size := 100000, partition_by := []::VARCHAR[], include_metadata := TRUE, header_text := NULL, metadata := map([]::VARCHAR[], []::VARCHAR[]), metadata_json_file := NULL, overwrite := FALSE, write_format_version := '1')

Returns:

VARCHAR

The function returns SQL without executing it; the R wrapper executes it through DBI.

SELECT duckhts_bam_convert_parquet_sql('sample.bam', 'sample.parquet', columns := ['QNAME','FLAG','RNAME','POS']);

Build COPY SQL for read_gff() output with Parquet metadata, GFF/tabix header text and selected columns, filters or partitions.

Signature:

duckhts_gff_convert_parquet_sql(path, output, columns := []::VARCHAR[], region := NULL, index_path := NULL, header := NULL, header_names := []::VARCHAR[], auto_detect := NULL, column_types := []::VARCHAR[], attributes_map := FALSE, attributes_list := FALSE, attributes_pairs := FALSE, strict := FALSE, where_sql := NULL, compression := 'zstd', row_group_size := 100000, partition_by := []::VARCHAR[], include_metadata := TRUE, header_text := NULL, metadata := map([]::VARCHAR[], []::VARCHAR[]), metadata_json_file := NULL, overwrite := FALSE, write_format_version := '1')

Returns:

VARCHAR

The function returns SQL without executing it; the R wrapper executes it through DBI.

SELECT duckhts_gff_convert_parquet_sql('annotations.gff3.gz', 'annotations/', columns := ['seqname','feature','start','end'], partition_by := ['feature']);

Build COPY SQL for read_tabix() output with Parquet metadata, header text and selected columns, filters or partitions.

Signature:

duckhts_tabix_convert_parquet_sql(path, output, columns := []::VARCHAR[], region := NULL, index_path := NULL, header := NULL, header_names := []::VARCHAR[], auto_detect := NULL, column_types := []::VARCHAR[], where_sql := NULL, compression := 'zstd', row_group_size := 100000, partition_by := []::VARCHAR[], include_metadata := TRUE, header_text := NULL, metadata := map([]::VARCHAR[], []::VARCHAR[]), metadata_json_file := NULL, overwrite := FALSE, write_format_version := '1')

Returns:

VARCHAR

The function returns SQL without executing it; the R wrapper executes it through DBI.

SELECT duckhts_tabix_convert_parquet_sql('regions.tsv.gz', 'regions.parquet', header_names := ['chrom','pos','value'], auto_detect := TRUE);

Construct a region-scoped BAM pileup with one row per covered position, emitting chrom, 1-based position, depth, observed bases, and Phred+33 qualities after SAM flag and MAPQ filtering. This is a compact htslib pileup view, not samtools mpileup text parity.

Signature:

read_pileup(path, region := NULL, index_path := NULL, min_mapq := 0, flag_mask := 1796)

Returns:

table

SELECT chrom, pos, depth FROM read_pileup('range.bam', region := 'CHROMOSOME_I:1-200') LIMIT 5;

Read full FASTA records or indexed regions with text or packed sequence output.

Signature:

read_fasta(path, region := NULL, index_path := NULL, gzi_path := NULL, sequence_encoding := 'string', scan_mode := 'auto')

Returns:

table

sequence_encoding='nt16' returns SEQUENCE UTINYINT[] using HTSlib nt16 codes instead of VARCHAR. NAME is the literal indexed header name; HTSlib quoting permits comma/colon names. Repeated requests retain one row per interval.

For bgzipped FASTA, gzi_path selects a non-colocated .gzi sidecar. scan_mode='sequential' streams/counts without indexed count paths and rejects region. NULL/empty region means no filter; empty comma-separated items error.

SELECT NAME, length(SEQUENCE) FROM read_fasta('ce.fa');

Read BED3-BED12 interval files with canonical typed columns and optional tabix-backed region filtering.

Signature:

read_bed(path, region := NULL, index_path := NULL, scan_mode := 'auto', error_policy := 'error')

Returns:

table

scan_mode := 'sequential' forces full-file streaming/counting instead of index-backed count paths and is incompatible with region.

read_bed is a lenient reader, not a BEDv1 validator (https://samtools.github.io/hts-specs/BEDv1.pdf): it also reads UCSC track files, skipping track and browser lines, and it only rejects data lines with fewer than three tab-delimited fields. Coordinates are not checked against BEDv1: a non-integer chromStart/chromEnd reads as NULL and chromEnd < chromStart is returned as is, so filter those in SQL when strict BED is required.

error_policy := 'error' aborts on data lines with fewer than three tab-delimited fields; 'skip' drops those lines; 'report' emits them with NULL normal columns and adds error (VARCHAR), line_number (BIGINT, 1-based physical line including headers/comments), and raw_line (VARCHAR, without newline). Good rows have NULL error and raw_line. Report requires a full-file scan, not a region query. Non-integer numeric fields become NULL, reversed intervals remain accepted, and blank/header/track/browser lines are ignored in all modes.

SELECT chrom, start, "end", name FROM read_bed('targets.bed') LIMIT 5;

Compute bedtools nuc-style nucleotide composition for supplied BED intervals or generated fixed-width bins over a FASTA reference. A failed reference fetch fails the query with the file and zero-based half-open interval; requested intervals are not silently omitted. For bgzipped FASTA, gzi_path may point to an explicit .gzi sidecar when it is not colocated with the FASTA.

Signature:

fasta_nuc(path, bed_path := NULL, bin_width := NULL, region := NULL, index_path := NULL, gzi_path := NULL, bed_index_path := NULL, include_seq := FALSE)

Returns:

table

SELECT chrom, start, "end", pct_gc FROM fasta_nuc('ce.fa', bin_width := 1000) LIMIT 5;

Create an empty session-scoped cgranges registry entry that can be populated with intervals and finalized for overlap queries.

Signature:

duckhts_cgranges_create(name)

Returns:

BOOLEAN

SELECT duckhts_cgranges_create('targets_idx');

Append an interval to a session-scoped cgranges registry entry before finalization. Labels may be BIGINT-like, DOUBLE, VARCHAR, or BOOLEAN.

Signature:

duckhts_cgranges_add(name, chrom, start, end[, label])

Returns:

BOOLEAN

SELECT duckhts_cgranges_add('targets_idx', 'chr1', 10, 20, 'exon1');

Finalize a populated cgranges registry entry and build its immutable overlap index for subsequent queries.

Signature:

duckhts_cgranges_index(name)

Returns:

BOOLEAN

SELECT duckhts_cgranges_index('targets_idx');

Destroy a session-scoped cgranges registry entry and release its indexed interval storage when it is not in active use.

Signature:

duckhts_cgranges_destroy(name)

Returns:

BOOLEAN

SELECT duckhts_cgranges_destroy('targets_idx');

Create, populate and finalize a session-scoped cgranges registry entry from the rows of a table or view, on the caller's connection.

Signature:

duckhts_cgranges_from_table(name, table_name, chrom_col, start_col, end_col)

Returns:

table(indexed BOOLEAN)

table_name is any relation visible to the calling connection, including TEMP tables, views and uncommitted rows; build it from a query with CREATE TEMP VIEW. chrom_col, start_col, end_col and label_col name its columns. Use duckhts_cgranges_from_table for unlabeled intervals or duckhts_cgranges_from_table_labelled with label_col. The entry is created, filled and finalized in one statement and returns one TRUE row, so duckhts_cgranges_index(...) is not needed afterwards.

The name must be new in the session. A failure while filling (NULL chrom, start or end, an unsupported label type, a coordinate outside int32) leaves the partly filled entry in place; remove it with duckhts_cgranges_destroy(...) before retrying.

Without label_col the label and interval_ordinal follow insertion order. A parallel scan of a large table does not fix that order, so pass label_col when a stable identity is needed.

duckhts_cgranges_from_table uses the base inputs. duckhts_cgranges_from_table_labelled requires label_col as an additional positional argument.

SELECT * FROM duckhts_cgranges_from_table('targets_idx', 'targets', 'chrom', 'start', 'end');

Build and index a labelled interval index from a table or view.

Signature:

duckhts_cgranges_from_table_labelled(name, table_name, chrom_col, start_col, end_col, label_col)

Returns:

table(indexed BOOLEAN)

Without label_col the label and interval_ordinal follow insertion order. A parallel scan of a large table does not fix that order, so pass label_col when a stable identity is needed.

duckhts_cgranges_from_table uses the base inputs. duckhts_cgranges_from_table_labelled requires label_col as an additional positional argument.

SELECT * FROM duckhts_cgranges_from_table_labelled('targets_idx', 'targets_table', 'chrom', 'start', 'end', 'label');

Vectorized scalar predicate for streaming provider rows through a finalized session-scoped cgranges index. Returns TRUE when the query interval overlaps at least one indexed interval, or when mode = 'contain' and it fully contains at least one indexed interval; NULL inputs return NULL.

Signature:

duckhts_cgranges_has_overlap(name, chrom, start, end[, mode])

Returns:

BOOLEAN

SELECT * FROM read_bed('queries.bed') WHERE duckhts_cgranges_has_overlap('targets_idx', chrom, start, "end");

Vectorized scalar overlap counter for streaming provider rows through a finalized session-scoped cgranges index. Returns the number of indexed intervals that overlap the query interval, or with mode = 'contain' the number fully contained by it; NULL inputs return NULL.

Signature:

duckhts_cgranges_count_overlaps(name, chrom, start, end[, mode])

Returns:

BIGINT

SELECT chrom, start, "end", duckhts_cgranges_count_overlaps('targets_idx', chrom, start, "end") AS n_targets FROM read_bed('queries.bed');

Vectorized scalar overlap expander for streaming provider rows through a finalized session-scoped cgranges index. Returns a LIST of hit STRUCTs that can be expanded with UNNEST, preserving provider columns while emitting one row per matching indexed interval. Because scalar return types are fixed, labels are returned as text with label_type describing the original cgranges label kind; NULL inputs return NULL.

Signature:

duckhts_cgranges_overlaps_list(name, chrom, start, end[, mode])

Returns:

STRUCT(interval_ordinal BIGINT, label VARCHAR, label_type VARCHAR, interval_chrom VARCHAR, interval_start INTEGER, interval_end INTEGER)[]

This is the bulk overlap path for any relation of probes. Expand the list in the SELECT list, as in the example: a select-list UNNEST streams the hits, whereas CROSS JOIN UNNEST over the same list plans a lateral join that was several times slower in measurement. Use row_number() OVER () or an existing key to identify probes.

SELECT q.chrom, q.start, q."end", unnest(duckhts_cgranges_overlaps_list('targets_idx', q.chrom, q.start, q."end")) AS hit FROM read_bed('queries.bed') AS q;

Query a finalized session-scoped cgranges registry entry and return one row per overlapping or containing indexed interval, preserving the original label type and interval coordinates.

Signature:

duckhts_cgranges_overlaps(name, chrom, start, end, mode := 'overlap', query_row_id := NULL)

Returns:

table

SELECT * FROM duckhts_cgranges_overlaps('targets_idx', 'chr1', 100, 150);

Read single-end, paired-end, or interleaved FASTQ files with optional legacy quality decoding. By default, FASTQ qualities are interpreted as modern Phred+33 input. Use sequence_encoding := 'nt16' to return SEQUENCE as UTINYINT[] and quality_representation := 'phred' to return QUALITY as UTINYINT[] instead of VARCHAR. input_quality_encoding accepts 'phred33', 'auto', 'phred64', or 'solexa64'. scan_mode := 'sequential' forces raw streaming/counting instead of index-backed count paths.

Signature:

read_fastq(path, interleaved := FALSE, mate_path := NULL, sequence_encoding := 'string', quality_representation := 'string', input_quality_encoding := 'phred33', scan_mode := 'auto')

Returns:

table

SELECT NAME, MATE FROM read_fastq('r1.fq', mate_path := 'r2.fq') LIMIT 5;

Read stored BigWig signal intervals as CHROM, START0, END0 and VALUE.

Signature:

read_bigwig(path, region := NULL, blocks_per_iteration := 64)

Returns:

table

Output intervals are zero-based half-open. region uses HTSlib one-based inclusive syntax; comma-separated requests merge per contig and emit each stored interval once.

Full scans parallelize over nonempty contigs and region scans over merged ranges, with worker-owned handles/iterators. blocks_per_iteration controls libBigWig batching, not DuckDB worker count. Local, native remote and browser wasm reads use HTSlib hFILE transport.

SELECT * FROM read_bigwig('scores.bw', region := 'chr1:100000-101000,chr2:200000-201000');

Aggregate canonical sequence and Phred+33 quality strings directly into exact read/base/Q20/Q30/Q40, nucleotide, quality-sum, and per-cycle sufficient statistics. The nested cycles list supports mean-quality, nucleotide-content, GC, and read-length curves without expanding one SQL row per base. Rows with any NULL input are ignored. Per-cycle state defaults to at most 1,048,576 cycles; pass a constant max_cycles per aggregate group to choose a larger explicit limit, up to 16,777,216.

Signature:

duckhts_fastq_qc(sequence, quality [, max_cycles])

Returns:

STRUCT

SELECT duckhts_fastq_qc(SEQUENCE, QUALITY) AS qc FROM read_fastq('reads.fastq.gz');

WITH q AS (SELECT duckhts_fastq_qc(SEQUENCE, QUALITY) AS qc FROM read_fastq('reads.fastq.gz')) SELECT cycle.* FROM q, UNNEST(qc.cycles) AS u(cycle);

Greedily select ranked Somalier candidate positions at a minimum genomic distance.

Signature:

duckhts_somalier_spacing(positions UBIGINT[], distance UBIGINT)

Returns:

BOOLEAN[]

positions must be ordered by AF rank within one chromosome. The returned mask has one element per input position; the first position wins each conflict, with equal positions conflicting. Positions are positive one-based integers; distance is positive and candidates are limited to one million per call.

Filter a typed read_bcf-style relation in SQL, group candidates by chromosome with list(CAST(POS AS UBIGINT) ORDER BY score, deterministic_tie_key), call this scalar on the ordered list, expand its Boolean mask, apply class-specific caps, and number the final rows by chromosome and position. rduckhts_somalier_find_sites() composes the complete recipe on the caller connection, including caller TEMP tables. No private connection or init-time macro is involved.

SELECT duckhts_somalier_spacing([100, 110, 109]::UBIGINT[], 10);

Import an already selected Somalier sites VCF/BCF as one canonical panel and population-frequency relation, autosomal and X/Y sites together.

Signature:

duckhts_somalier_import_sites(path, assembly_name, max_sites := 1000000)

Returns:

table

Each biallelic SNV is oriented into lexical allele_a/allele_b order. INFO/AF is the source ALT frequency and is transformed with REF/ALT, so population_b_af always describes allele_b. Dense zero-based site_index numbers every autosomal site first, in Somalier v0.3.4's lexical region then position order, and then the X/Y sites in the same order, so the autosomal ordinals are 0..n-1 whatever the X/Y region names sort as.

INFO/AF must declare Number=A,Type=Float. Retained records must be distinct uppercase A/C/G/T biallelic SNVs with one finite AF in [0,1]. Records on the exact Somalier v0.3.4 X/Y aliases (X, chrX, NC_000023.10, NC_000023.11, Y, chrY, NC_000024.9, NC_000024.10) are kept as sex-chromosome sites and at least one autosomal record is required to normalise their depth. No pseudo-autosomal filtering happens here or at count time, as in Somalier: choose X sites outside the PAR when selecting the sites file (rduckhts_somalier_find_sites does for human builds). Existing FILTER values are retained as provenance because this imports an already selected sites file; it does not apply selection policy.

The returned relation directly satisfies the canonical panel and population-frequency contracts and retains source REF, ALT, ALT frequency and FILTER. This is not Somalier find-sites: population AF/AN and QC filtering, interval exclusions, linkage spacing and target-frequency ranking are a separate panel-selection method. max_sites is an explicit input limit.

CREATE TABLE fingerprint_panel AS SELECT * FROM duckhts_somalier_import_sites('sites.vcf.gz', 'GRCh38');

Extract a complete panel-aligned A/B/other count relation from VCF/BCF FORMAT/AD.

Signature:

duckhts_somalier_vcf_counts(path, panel_table, samples := NULL, filter_policy := 'pass_or_unapplied')

Returns:

table

panel_table is the canonical typed panel relation with one assembly, dense zero-based site_index, region, one-based position and lexical uppercase allele_a/allele_b. Output has exactly one row per selected header sample and panel ordinal, including absent source records. samples uses HTSlib selection syntax and source_sample_index retains the original zero-based header ordinal.

FORMAT/AD must declare Number=R,Type=Integer. Exact REF/ALT identity selects A and B; remaining declared-allele AD slots are summed as other. Counts are never inferred from GT or DP. Missing records, alleles or AD produce three NULL counts and a named unavailable status; three zeros are measured evidence. Duplicate source records or allele identities, malformed cardinality and negative counts error.

Every panel site, autosomal or X/Y, is counted in the one pass. Symbolic alleles are unavailable in this scope. filter_policy='pass_or_unapplied' accepts PASS and dot but makes named failures unavailable; 'include_all' uses their AD; 'error' rejects them. X/Y sites ignore filter_policy and always use their AD, as Somalier extraction does. source_method, count_scope, status, FILTER, record/sample ordinals and matched allele slots preserve extraction provenance. Full scans are sequential and retain physical source-record identity; Parquet is optional downstream storage, not an input requirement.

CREATE TABLE allele_counts AS SELECT * FROM duckhts_somalier_vcf_counts('cohort.bcf', 'fingerprint_panel');

Extract complete panel-aligned A/B/other base counts from one indexed BAM or CRAM source.

Signature:

duckhts_somalier_bam_counts(source_path, panel_table, sample_id, reference_path, panel_parquet := NULL, index_path := NULL, reference_index_path := NULL, min_mapq := 1, min_baseq := 0, require_flags := 0, exclude_flags := 1796, overlap_policy := 'hileup_v0.1.0', decompression_threads := 0, worker_count := 1, max_depth := 100000, max_overlap_qnames := 100000, max_sites := 1000000, max_region_bytes := 67108864, remote_block_bytes := 1048576, remote_cache_bytes := 67108864, reference_cache_bytes := 67108864)

Returns:

table

panel_parquet supplies the canonical typed six-column panel; panel_table is retained for signature compatibility and must be NULL. The panel is read from a local Parquet file when the query is bound, on the calling thread: a private in-memory DuckDB instance is opened and closed inside bind to materialize and validate it once. A prepared statement keeps that panel and reads nothing again. worker_count then partitions the sites into at most 64 DuckDB scan jobs, each counting into its own copy of the sites. Each worker-local job owns one alignment handle, index, multi-region iterator, pileup, reference handle and bounded overlap workspace; no mutable htslib or faidx state is shared. DuckDB's connection thread setting limits concurrent jobs. decompression_threads separately controls htslib workers per alignment handle. Parallel output order is unspecified; use ORDER BY site_index when order matters.

Every panel site, autosomal or X/Y, is counted in the one scan. A and B are exact uppercase panel bases; other counts every other observed base. Valid uncovered sites emit measured 0/0/0. Missing reference contigs or positions, reference mismatches, alignment-header contig absence and positions beyond a declared alignment contig emit NULL counts with a named unavailable status. Deletions and reference skips are not observed bases. MAPQ, base quality, required/excluded flags and overlap suppression are explicit per-call settings.

The private reader keeps no connection into the caller's database, so a closed database is released; it also cannot see the caller's tables, so write the panel with COPY ... TO 'panel.parquet' first (rduckhts_somalier_bam_counts does this for a panel_table). panel_parquet must be a local ordinary Parquet file: remote URLs are not read. max_sites bounds the shared panel. max_depth, max_overlap_qnames and max_region_bytes bound each active scan job; concurrent workspace can therefore grow with worker_count. remote_block_bytes, remote_cache_bytes and reference_cache_bytes apply per worker-owned handle, and zero disables that cache override. Explicit non-colocated BAM/CRAM and FASTA indexes are honored; CRAM keeps the alignment file's original @SQ lengths and does not create a default FASTA-index sidecar. overlap_policy='hileup_v0.1.0' pins encounter-order mate suppression only, while 'none' counts every retained observation. The panel reader never runs on a DuckDB worker thread: Rtools-built DuckDB hosts destroy a thread_local on freed memory when a worker that closed a DuckDB instance exits (see Runtime compatibility in the README).

CREATE TABLE allele_counts AS SELECT * FROM duckhts_somalier_bam_counts('sample.bam', NULL, 'sample-1', 'reference.fa', panel_parquet := 'fingerprint_panel.parquet');

Solve a nearest-positive-definite constrained ancestry projection from aggregated PC products.

Signature:

duckhts_ancestry_proportions(x_pc_major, y_pc, group_count, sum_to_one := true)

Returns:

DOUBLE[]

x_pc_major is a flat DOUBLE[] in PC-major, group-minor order for X = P^T F_ref; y_pc is a DOUBLE[] with one value per PC for (P^T f) times the shrinkage correction. Groups are ordered identically within each PC. 1 to 30 groups and 1 to 64 PCs are supported. All values must be finite. sum_to_one=true enforces nonnegative q summing to one; false allows a sum at most one. The scalar returns full-precision group proportions in input order; the R wrapper rounds displayed proportions to seven decimals after correlation gating.

The Gram matrix X^T X follows Matrix::nearPD defaults: Dykstra alternating projections, relative eigen cutoff 1e-6, infinity-norm convergence tolerance 1e-7 and at most 100 iterations, then eigenvalue floor 1e-8 times the largest absolute eigenvalue with diagonal rescaling. A nonconverged repair errors. The native solver is independent of bigsnpr and quadprog code.

rduckhts_ancestry_proportions normalises long or keyed wide reference products to one row per locus with PC and group columns, checks long-reference completeness, unique input-matched reference loci and finite values at aligned sites, and matches sample_id, chromosome, positive one-based whole-number position and alleles from the caller's input-frequency relation. It reverses input frequency (1 - f) for reversed alleles, audits strand flips, and drops ambiguous, duplicate, missing and mismatched input sites. A bounded aligned Parquet relation holds the retained physical rows. Per-sample sums of loading * reference_frequency form X in PC-major, group-minor order; sums of loading * aligned_frequency times correction form y in PC order. The native scalar solves for q at full precision; cor(F_ref q, f) and per-group cor_each are evaluated before failed proportions are gated to NULL and returned proportions rounded. The scratch Parquet is removed after the call; no private connection or macro is created.

Input genotypes are diploid dosage / 2; NULL genotypes and frequencies are dropped, not imputed. bigsnpr imputes missing individual genotypes before calling snp_ancestry_summary, so parity comparisons must use the same retained sites. The default min_cor is 0.4 as in bigsnpr 1.12.21; callers may select a stricter gate for cohorts or low-depth count frequencies.

SELECT duckhts_ancestry_proportions([1.0, 0.0, 0.0, 1.0], [0.25, 0.75], 2, true);

Decode runs of homozygosity from one sample's sorted per-chromosome site lists, reproducing bcftools roh.

Signature:

duckhts_roh_segments(positions, af, pl, map_pos, map_cm, rec_rate, hw_to_az, az_to_hw)

Returns:

STRUCT(start BIGINT, "end" BIGINT, n_markers INTEGER, quality DOUBLE)[]

duckhts_roh_segments(positions, af, pl, map_pos, map_cm, rec_rate, hw_to_az, az_to_hw) takes phred-scaled likelihoods as INTEGER[][] (three values per site: RR, RA, AA). duckhts_roh_segments(positions, af, dosage, gt_error, map_pos, map_cm, rec_rate, hw_to_az, az_to_hw) takes INTEGER[] genotype dosages 0, 1 or 2 and a phred gt_error, as bcftools roh -G does. duckhts_roh_segments(positions, af, ref_count, alt_count, seq_error, contamination, map_pos, map_cm, rec_rate, hw_to_az, az_to_hw) takes INTEGER[] read counts of the other allele and of the allele af refers to, a per-read seq_error in (0, 0.5) and a contamination fraction in [0, 1); duckhts_roh_counts describes that read model, which is a DuckHTS extension with no bcftools counterpart. All arguments are positional; the duckhts_roh macro supplies the defaults.

positions is a BIGINT[] of one-based coordinates, sorted ascending and all in 1..2147483647; af is a DOUBLE[] of ALT allele frequencies and pl or dosage a list of the same length, one entry per site. Build them per sample and chromosome with list(... ORDER BY pos). Lists of unequal length, NULL positions, unsorted positions, an af outside [0, 1], a dosage outside 0..2 and a hw_to_az or az_to_hw outside [0, 1] are errors that name the problem. A NULL list returns NULL; an empty list returns an empty list. One list is decoded at a time, so memory follows the longest list.

As in bcftools roh, a site is skipped when its position repeats the previous entry, its af is NULL, NaN or exactly 0, its genotype evidence is NULL or has NULL elements, its PL does not have three values or has a negative value, all three PL values are equal, or the likelihoods sum to zero. PL values above 255 count as 255. Skipped sites are not markers and do not break a run.

A two-state hidden Markov model, autozygous (AZ) and Hardy-Weinberg (HW), with initial probabilities 1/2, emissions P(D|AZ) = P(D|RR)(1-f) + P(D|AA)f and P(D|HW) = P(D|RR)(1-f)^2 + 2P(D|RA)f(1-f) + P(D|AA)f^2, and transitions hw_to_az (default 6.7e-8) and az_to_hw (default 5e-9) per base pair compounded over the physical distance between sites. Segments are the Viterbi AZ runs; quality is the mean forward-backward phred score of the sites in the run. start and end are the first and last marker of the run, inclusive.

map_pos (BIGINT[], one-based, strictly ascending) and map_cm (DOUBLE[], cumulative centimorgans) are the nodes of a genetic map, the positions and third column of an IMPUTE2 map file. When both are given, the crossover terms between two sites are scaled by the map distance in Morgans, interpolated as bcftools roh -m does (the slope between the nodes that bracket the two sites, zero outside the map), and rec_rate, when given, scales it further. With no map and a positive rec_rate, the scale is the physical distance times rec_rate, as bcftools roh -M does. With neither, only the per-base-pair compounding applies. NULL map_pos and map_cm mean no map; giving one without the other, lists of unequal length or an empty map are errors.

The arithmetic follows the pinned bcftools vcfroh.c and HMM.c operation for operation, so positions and marker counts equal bcftools roh and qualities agree to the one decimal bcftools prints; differences could arise only from the floating-point contraction choices of a different compiler or CPU.

SELECT duckhts_roh_segments([100, 200, 300], [0.3, 0.3, 0.3], [[0, 30, 60], [0, 30, 60], [0, 30, 60]], NULL::BIGINT[], NULL::DOUBLE[], NULL::DOUBLE, 6.7e-8, 5e-9);

Find runs of homozygosity in a VCF/BCF with the bcftools roh model, using allele frequencies from an INFO tag.

Signature:

duckhts_roh(path, af_tag, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

path is a VCF/BCF with FORMAT/PL (default) or FORMAT/GT (gt_error set), read with read_bcf(tidy_format := true) on the caller's connection. af_tag names an INFO tag declared Type=Float, Number=A; its first value is the ALT allele frequency. Sites with an absent or missing tag, or a frequency of exactly 0, are skipped as bcftools roh skips them. Records with more than one ALT, no ALT, or a first ALT that is symbolic (<*>, ) are not analysed. samples uses read_bcf's selectors.

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 (mean forward-backward phred score), the RG lines of bcftools roh. Rows are not ordered; ORDER BY sample, chrom, start. Chromosomes are decoded separately and a run never spans two.

gt_error NULL uses PL. A numeric gt_error (phred, bcftools -G, for example 30) uses diploid GT calls 0/0, 0/1 and 1/1 and ignores PL; other or missing genotypes are skipped, and FORMAT/PL need not exist. gt_error is an error of the site: the called genotype has likelihood close to 1 and a genotype one allele away has 10^(-gt_error/10), so one call that contradicts a run counts against it by that factor and no more, whatever caused it. A run can therefore continue through an isolated heterozygous call; a larger gt_error makes that rarer.

hw_to_az and az_to_hw are the per-base-pair transition probabilities (bcftools -a and -H). rec_rate is a constant recombination rate per base pair (bcftools -M). duckhts_roh_map takes genetic_map, a table name, with columns chrom, pos and cm (cumulative centimorgans at one-based positions, for example an IMPUTE2 map), interpolated as bcftools -m does; chromosomes without map rows are skipped, as bcftools skips a chromosome with no map file.

DuckDB manages, and can spill, the scan and the joins. The sites of each sample and chromosome are held in native buffers of 16 bytes per site until that sample and chromosome is decoded, so memory follows samples times sites. max_sites caps one sample and chromosome (default 20,000,000; at most 100,000,000). max_site_bytes (default 4 GiB) bounds the buffers held at once: a decode does not grow a buffer when the bytes held by all ROH decodes in the process would pass its own max_site_bytes. Decodes that run at the same time with different values each apply their own, so the process total can pass the smaller value while the decode with the larger value allocates. Exceeding either limit is an error: decode fewer samples per query with samples, or raise the limit. When a query fails while the sites are being collected with 5 or more threads, DuckDB 1.5.x does not destroy the collecting states (https://github.com/duckdb/duckdb/issues/26871); their buffers stay charged against max_site_bytes until the process ends. On DuckDB 1.5.x, an out-of-memory error while these states, or the ordered genetic_map lists, are held can end the process instead of returning the error (https://github.com/duckdb/duckdb/issues/26441); raise memory_limit or decode fewer samples per query. Decoding one sample and chromosome also needs approximately 36.25 bytes per site of workspace in each thread: 4 bytes for position, 16 for emissions, 16 for forward/posterior probabilities and two traceback bits. For capacity n, these arrays use 36*n + ceil(n/4) + 16 bytes, excluding fixed tables and output segments. Model parameters and limits are checked even when the input has no record.

Duplicate positions on a chromosome keep the first record to arrive, as bcftools roh keeps the first record in the file. That is the file order when one scan thread reads both records; deduplicate upstream when it matters. With a frequency relation (duckhts_roh_af_table, duckhts_roh_ancestry), only the records that match a frequency row take part. bcftools' default AC/AN frequencies, --AF-dflt, --estimate-AF, --include/--exclude, --skip-indels, --ignore-homref, --buffer-size and --viterbi-training are not offered.

The chrom of genetic_map is compared byte for byte with the contig names of the sites. A chromosome that the map spells differently (1 and chr1) is skipped like a chromosome without map rows, without an error. Rename the map's contigs in a view first.

duckhts_roh uses the base inputs. duckhts_roh_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh('cohort.vcf.gz', 'AF') ORDER BY sample, chrom, start;

Find runs of homozygosity using an INFO allele frequency tag and a genetic map.

Signature:

duckhts_roh_map(path, af_tag, genetic_map, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

duckhts_roh uses the base inputs. duckhts_roh_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_map('cohort.vcf.gz', 'AF', 'genetic_map') ORDER BY sample, chrom, start;

Find runs of homozygosity in a VCF/BCF with the bcftools roh model, using allele frequencies from a caller relation.

Signature:

duckhts_roh_af_table(path, af_table, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

af_table is a table or view name with columns chrom, pos (one-based), ref, alt and af, like a bcftools --AF-file. A record takes the frequency of the row matching its chrom, pos, REF and its ALT alleles joined by commas; a record without a matching row is not used and does not claim a repeated position, and a record whose row has a NULL af is skipped. bcftools roh --AF-file chooses differently at a repeated position; duckhts_roh_af_table_bcftools follows it. Give one row per key; af must be in [0, 1], and exactly 0 skips the site. This is the entry point for per-sample frequencies, for example ancestry-tuned ones, by running it once per relation.

Everything else is as duckhts_roh: the same output, genotype evidence (PL, or GT with gt_error), transition parameters, genetic-map form, memory limits and limits.

The chrom of af_table is compared with the VCF's CHROM byte for byte. When the two spell contigs differently (1 and chr1), no record finds a frequency and the function returns no rows, without an error. Rename the contigs of af_table to the VCF's spelling in a view first; duckhts_contig_key() gives a common key.

duckhts_roh_af_table uses the base inputs. duckhts_roh_af_table_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_af_table('cohort.vcf.gz', 'site_frequencies') ORDER BY sample, chrom, start;

Find runs of homozygosity using matching allele frequencies and a genetic map.

Signature:

duckhts_roh_af_table_map(path, af_table, genetic_map, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

Everything else is as duckhts_roh: the same output, genotype evidence (PL, or GT with gt_error), transition parameters, genetic-map form, memory limits and limits.

duckhts_roh_af_table uses the base inputs. duckhts_roh_af_table_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_af_table_map('cohort.vcf.gz', 'site_frequencies', 'genetic_map') ORDER BY sample, chrom, start;

Find runs of homozygosity like duckhts_roh_af_table, choosing the record of a repeated listed position as bcftools roh --AF-file does.

Signature:

duckhts_roh_af_table_bcftools(path, af_table, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296, af_order := NULL)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

af_table is as for duckhts_roh_af_table. bcftools reads only the first line of a position, and a relation has no line order, so without af_order af_table must give at most one row per position: more is an error that names the position.

af_order names a column of af_table that gives the order of its lines, for example a line number the caller adds, or FILE_OFFSET when the frequencies come from read_bcf(). The name is matched without case. At each position the row with the least af_order value stands for the first line, and only that row is used. A NULL af_order value, a value repeated within a position (each error names the position), or a name that is not a column of af_table is an error.

In file order (read_bcf's FILE_OFFSET), as bcftools 1.23 and 1.24 do: a record whose indel class differs from the row's is passed over when the next record of the file has the same POS, on any contig and listed or not, since bcftools compares POS only; the row is an indel when any of its alleles is longer than one character, and a record is one when HTSlib's bcf_get_variant_types() sets VCF_INDEL (an MNP such as AC>GT and a complex ACG>AT are not; duckhts__htslib_is_indel(ref, alts) is the port). A record with no non-symbolic ALT, or several, claims nothing; the first remaining record claims the position, which is skipped unless that record's REF and ALT equal the row's. A claiming record whose first ALT is symbolic (<*>, ) is skipped too, where bcftools would decode its second allele.

path must be a BGZF-compressed VCF or BCF; other compression is an error. It must also be sorted, as bcftools roh requires, so that the records of a position and of a contig are contiguous in the file; that is not checked.

Besides the genotype scan of duckhts_roh_af_table, one more scan without samples reads CHROM, POS, REF, ALT and FILE_OFFSET; records whose POS af_table does not list are only counted per contig, so the extra state is one row per contig and per record at a listed POS. af_table is scanned three times rather than held; af_order adds its column to two of those scans and no other scan.

Everything else is as duckhts_roh: the same output, genotype evidence (PL, or GT with gt_error), transition parameters, genetic-map form and memory limits.

duckhts_roh_af_table_bcftools uses the base inputs. duckhts_roh_af_table_bcftools_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_af_table_bcftools('cohort.vcf.gz', 'site_frequencies') ORDER BY sample, chrom, start;

SELECT * FROM duckhts_roh_af_table_bcftools('cohort.vcf.gz', 'af_file_lines', af_order := 'line') ORDER BY sample, chrom, start;

Find runs of homozygosity using bcftools AF-file site selection and a genetic map.

Signature:

duckhts_roh_af_table_bcftools_map(path, af_table, genetic_map, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296, af_order := NULL)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

Everything else is as duckhts_roh: the same output, genotype evidence (PL, or GT with gt_error), transition parameters, genetic-map form and memory limits.

duckhts_roh_af_table_bcftools uses the base inputs. duckhts_roh_af_table_bcftools_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_af_table_bcftools_map('cohort.vcf.gz', 'af_file_lines', 'genetic_map', af_order := 'line') ORDER BY sample, chrom, start;

Find runs of homozygosity using each sample's ancestry-weighted allele frequencies.

Signature:

duckhts_roh_ancestry(path, reference_table, proportions_table, af_clamp := 1e-3, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

reference_table is a long relation with chromosome, position (one-based), allele_a, allele_b, group_id and frequency (frequency of allele_b). proportions_table has sample_id, group_id and proportion. Both relations must use exactly the same group IDs; every called VCF sample needs proportions. Contig names are matched with duckhts_contig_key(), resolved once per distinct name. It removes one leading chr case-insensitively, writes M/MT as MT and uppercases X and Y, so chr1 and 1 match, as do chrX and X. Accessions, patches and numeric sex chromosomes are not mapped. Each site must have one allele orientation across its group rows, and exactly one row per group.

For each called site and sample, AF is the sum over groups of proportion times group frequency. The reference is aggregated once per site and proportions once per sample, then joined to called sites. REF=allele_a and ALT=allele_b uses the weighted frequency; the reversed orientation uses 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. af_clamp in [0, 0.5) clamps AF to [af_clamp, 1-af_clamp]; zero disables clamping. The default avoids interpreting a population frequency of zero as impossible.

The output, genotype evidence, transitions, genetic-map form, memory limits and limits are as duckhts_roh_af_table.

duckhts_roh_ancestry uses the base inputs. duckhts_roh_ancestry_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_ancestry('cohort.vcf.gz', 'reference_long', 'sample_proportions', gt_error := 30) ORDER BY sample, chrom, start;

Find runs of homozygosity using ancestry-weighted frequencies and a genetic map.

Signature:

duckhts_roh_ancestry_map(path, reference_table, proportions_table, genetic_map, af_clamp := 1e-3, hw_to_az := 6.7e-8, az_to_hw := 5e-9, gt_error := NULL, rec_rate := NULL, samples := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

The output, genotype evidence, transitions, genetic-map form, memory limits and limits are as duckhts_roh_af_table.

duckhts_roh_ancestry uses the base inputs. duckhts_roh_ancestry_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_ancestry_map('cohort.vcf.gz', 'reference_long', 'sample_proportions', 'genetic_map', gt_error := 30) ORDER BY sample, chrom, start;

Find runs of homozygosity from allele read counts, with sequencing error, optional contamination, allele balance and a site artefact term.

Signature:

duckhts_roh_counts(counts_table, seq_error := 1e-3, contamination := 0.0, hw_to_az := 6.7e-8, az_to_hw := 5e-9, rec_rate := NULL, max_sites := 20000000, max_site_bytes := 4294967296, allele_balance := 0.5, artefact_weight := 0.0)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

counts_table is a table or view with sample_id, chrom, pos (one-based), ref_count, alt_count and af, one row per sample and site. alt_count counts reads showing the allele whose population frequency is af; ref_count counts reads showing the other allele. Any frequency source composes in SQL first: an INFO tag, a frequency table or ancestry-weighted frequencies. Somalier site counts map as ref_count = a, alt_count = b and af = frequency of allele_b.

seq_error, contamination, allele_balance and artefact_weight mean what duckhts_count_error_fit reports under these names. Use duckhts_roh_counts_fitted(counts_table, fit_table, ...) to supply all four per sample from a caller-named relation; in R, pass fit_table to rduckhts_roh_counts.

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). Reads are independent, so a site's alt count is binomial. With chance artefact_weight the alt count is instead uniform on 0 to depth, whatever the genotype: each genotype's likelihood is (1 - artefact_weight) * binomial + artefact_weight / (depth + 1), the form duckhts_count_error_fit uses. The resulting genotype likelihoods enter the bcftools roh model in place of PL. This emission is a DuckHTS extension; bcftools roh has no read-count mode.

allele_balance is the share of a heterozygote's reads that show the counted allele; values below 0.5 model reference-biased mapping. artefact_weight is the chance that a site's counts do not reflect its genotype (a mapping artefact, a paralogous read or a bad panel site). With artefact_weight w > 0, one site at depth d changes the odds between the two HMM states by a factor of at most 1 + (d + 1)(1 - w) / w (about 10,300 at d = 30 and w = 0.003), as gt_error caps one called genotype at 10^(gt_error / 10). The defaults, 0.5 and 0, give the decode without these terms, bit for bit.

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, so the read-count decode corresponds to the GT decode with no tolerated genotype error: it ends a run at a heterozygous site that duckhts_roh with gt_error := 30 would pass through, and reports fewer bases in runs. benchmarks/benchmark_roh_counts_validation.md measures this on three 1000 Genomes samples, at artefact_weight 0. A positive artefact_weight bounds each site's influence as gt_error does.

allele_balance must be in [0.35, 0.65] and artefact_weight in [0, 0.05]; a value outside is an error, never clipped, also when counts_table is empty. The fit measures allele_balance 0.4974 to 0.4984 and artefact_weight 0.18% to 0.30% of sites on three 30x chr20 samples (benchmarks/benchmark_count_error_fit.md); the ranges hold such values with wide margins. Beyond them a value points to copy number, a mixture or a broken panel, which one number per sample does not describe. The SQL tests check that allele_balance +-0.01 and artefact_weight up to 0.004 move no run end of the synthetic titration fixture by more than one site.

Sites with no reads, a NULL count, or an af that is NULL, NaN or 0 are skipped; a repeated position keeps one row, as in the other ROH macros. The relation need not be sorted. Negative counts, fractional counts, seq_error outside (0, 0.5) and contamination outside [0, 1) are errors.

Biallelic sites only; one error rate for all reads, with no base- or mapping-quality, duplicate or strand modelling, and binomial counts (the spreads of duckhts_count_error_fit are not used). allele_balance is one value per sample, but reference-biased mapping is stronger where flanking heterozygosity is high; such heterozygotes can still look homozygous and lengthen runs in divergent regions. The contaminant genotype is averaged per read rather than drawn once per site as in the fit; both give a read the same chance. A contaminant from a different population than af describes is approximated by af.

Output, transitions, rec_rate, the genetic-map form and the memory limits are as duckhts_roh, with 24 bytes per site.

duckhts_roh_counts uses the base inputs. duckhts_roh_counts_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_counts('site_counts', contamination := 0.02) ORDER BY sample, chrom, start;

Find runs of homozygosity from read counts using a genetic map.

Signature:

duckhts_roh_counts_map(counts_table, genetic_map, seq_error := 1e-3, contamination := 0.0, hw_to_az := 6.7e-8, az_to_hw := 5e-9, rec_rate := NULL, max_sites := 20000000, max_site_bytes := 4294967296, allele_balance := 0.5, artefact_weight := 0.0)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE)

Output, transitions, rec_rate, the genetic-map form and the memory limits are as duckhts_roh, with 24 bytes per site.

duckhts_roh_counts uses the base inputs. duckhts_roh_counts_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_counts_map('site_counts', 'genetic_map', contamination := 0.02) ORDER BY sample, chrom, start;

Find runs of homozygosity from allele read counts using per-sample fitted error, contamination, allele balance and artefact weight, with the share of each run that the artefact weight's interval does not move and the bases a higher weight would add next to it.

Signature:

duckhts_roh_counts_fitted(counts_table, fit_table, hw_to_az := 6.7e-8, az_to_hw := 5e-9, rec_rate := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE, stable_fraction DOUBLE, gained_bases BIGINT)

fit_table is a caller-named table or view with one row per sample and columns sample, seq_error, contamination, allele_balance, artefact_weight and artefact_weight_sd, such as the output of duckhts_count_error_fit; other columns are ignored. sample joins to counts_table.sample_id, with a NULL sample matching a NULL sample_id as in duckhts_roh_counts and duckhts_count_error_fit. For each sample of counts_table, a missing fit row, a NULL value of one of the four fitted numbers or more than one fit row is an explicit error naming the sample and problem, whatever DuckDB's scalar_subquery_error_on_multiple_rows setting. A relation without one of the five columns is a binder error. Fit rows of samples that counts_table does not hold are not checked, so empty counts give no rows; the other arguments are still checked. Each value must be in the range of the scalar argument of duckhts_roh_counts: seq_error in (0, 0.5), contamination in [0, 1), allele_balance in [0.35, 0.65] and artefact_weight in [0, 0.05], finite; a fitted value outside is an error, not clipped. artefact_weight_sd must be finite and at least 0; a NULL counts as 0. Fit status filtering belongs in the caller's SQL. The spreads are not used.

Counts, emissions, the run columns, transitions, genetic-map form and memory limits are identical to duckhts_roh_counts: sample through quality equal its output at each sample's four numbers. Only the four fitted numbers vary by sample. In R, pass fit_table to rduckhts_roh_counts; its scalar seq_error, contamination, allele_balance and artefact_weight arguments are unused.

The runs depend on artefact_weight more than on the other fitted numbers, and the fit gives an interval for it. Each sample is 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]. The rows are the runs of the decode at w; stable_fraction is the share of a run's bases that are also in a run of both other decodes, from 0 to 1. A run with stable_fraction 1 does not depend on where w lies within one standard error; a low value marks a run whose call hinges on the artefact weight. Filter on it, or weight run bases by it, to report the runs or bases the fit's interval does not move. stable_fraction sees only bases a run would lose; gained_bases guards the other direction (see Gained bases). A NULL artefact_weight_sd (a fit with fewer than three blocks) counts as 0, so every run gets stable_fraction 1 and gained_bases 0; check the fit's blocks or status as well. The two extra decodes share the site lists and the kernel workspace, so memory is that of duckhts_roh_counts plus, for the list a thread decodes, 24 bytes per run of the two extra decodes and 16 per run at w, and decode time is about three times as long.

gained_bases is the number of bases that the decode at w + artefact_weight_sd calls next to a run at w but outside every run at w, from the same decode as stable_fraction. A run at w + sd meets a run at w when the two share a base or abut (the end of one is one base before the start of the other). Each base of a run at w + sd that meets some run at w, and that lies outside every run at w, counts for exactly one run at w: a base before the first run at w that the run at w + sd meets counts for that run; any other base counts for the last run at w before it. So a gap that the higher weight closes between two runs counts once, for the run before the gap, and an extension past a run's end counts for that run. Bases of a run at w + sd that meets no run at w belong to no row and are not reported; a higher weight can call such new runs. A high gained_bases marks a run that a higher weight would extend or join to its neighbour.

duckhts_roh_counts_fitted uses the base inputs. duckhts_roh_counts_fitted_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_counts_fitted('site_counts', 'sample_fits') ORDER BY sample, chrom, start;

SELECT sample, sum(length) AS run_bases, sum(length * stable_fraction) AS stable_run_bases FROM duckhts_roh_counts_fitted('site_counts', 'sample_fits') GROUP BY sample;

Decode read-count ROH with per-sample fitted parameters and a genetic map.

Signature:

duckhts_roh_counts_fitted_map(counts_table, fit_table, genetic_map, hw_to_az := 6.7e-8, az_to_hw := 5e-9, rec_rate := NULL, max_sites := 20000000, max_site_bytes := 4294967296)

Returns:

table(sample VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, length BIGINT, n_markers INTEGER, quality DOUBLE, stable_fraction DOUBLE, gained_bases BIGINT)

duckhts_roh_counts_fitted uses the base inputs. duckhts_roh_counts_fitted_map requires genetic_map as an additional positional argument.

SELECT * FROM duckhts_roh_counts_fitted_map('site_counts', 'sample_fits', 'genetic_map') ORDER BY sample, chrom, start;

Fit read error, contamination and allele balance of each sample by maximum likelihood from allele read counts at sites of known population frequency.

Signature:

duckhts_count_error_fit(counts_table, block_bases := 10000000, freq_bins := 64, max_depth := 1000, min_block_sites := 1000, max_cells := 5000000, max_cell_bytes := 1073741824)

Returns:

table(sample VARCHAR, sites BIGINT, reads BIGINT, mean_depth DOUBLE, blocks INTEGER, seq_error DOUBLE, contamination DOUBLE, contamination_sd DOUBLE, homozygosity_excess DOUBLE, allele_balance DOUBLE, spread_hom DOUBLE, spread_het DOUBLE, artefact_weight DOUBLE, artefact_weight_sd DOUBLE, log_likelihood DOUBLE, contamination_relative DOUBLE, log_likelihood_relative DOUBLE, status VARCHAR, method VARCHAR)

counts_table is the relation duckhts_roh_counts takes: sample_id, chrom, pos (one-based), ref_count, alt_count and af, one row per sample and site, where alt_count counts reads showing the allele whose population frequency is af. Any frequency source composes in SQL first, for example ancestry-weighted frequencies. Sites with a NULL count, no reads, more than max_depth reads (the two counts are summed as BIGINT), or an af that is NULL, NaN, 0 or 1 are skipped; a sample whose sites are all skipped keeps its row. Negative or fractional counts are errors, as is a NULL, non-positive or fractional pos, or a fractional option. Block k of a chromosome holds the positions k * block_bases + 1 to (k + 1) * block_bases.

A site's sample genotype has 0, 1 or 2 copies of the counted allele with Hardy-Weinberg chances at af, moved toward homozygosity by homozygosity_excess; a second genome contributes a share contamination of the reads, with its own genotype drawn once per site (unrelated: Hardy-Weinberg; relative: one allele of the sample's genotype and one draw at af). A read shows the allele with chance seq_error, allele_balance or 1 - seq_error for 0, 1 or 2 copies, mixed over the two genomes. The alt count is beta-binomial with spread spread_hom at homozygous and spread_het at heterozygous sites (0 is the binomial). With chance artefact_weight a site's alt count is uniform on 0..depth whatever its genotype. All seven are fitted together by Nelder-Mead from two starts on the histogram of (block, frequency bin, depth, alt count) cells, so the likelihood depends on a site only through its counts and the mean frequency of its cell.

One row per sample: the fitted seven numbers, the log-likelihood, and contamination_relative with log_likelihood_relative from the same cells with a parent or child as the second genome; the sites and reads used and their mean depth. Each block of block_bases bases with at least min_block_sites sites is refitted for contamination with the other parameters held at the whole-sample fit, and contamination_sd is the standard deviation of those block estimates divided by the square root of their number (NULL with fewer than three blocks). Each such block is also refitted for artefact_weight alone, the other parameters held the same way, and artefact_weight_sd is formed from those estimates as contamination_sd is; duckhts_roh_counts_fitted uses it for stable_fraction. status is ok, no_sites (a sample whose sites were all skipped: sites 0, the numbers NULL), few_sites (under 1,000 sites), few_blocks (under three refitted blocks), at_bound (contamination, seq_error or artefact_weight at its upper bound), no_convergence or no_convergence_relative; method names the fit and its version.

seq_error is the chance that a read shows the other allele at a homozygous site. With a mask of known variants absent, it holds site artefacts that artefact_weight does not absorb. homozygosity_excess is a nuisance parameter: it also absorbs any difference between the sample's population and the frequencies given, so it is not inbreeding. contamination is the share of reads of a second genome of the frequencies' population; a second genome of another population biases it upward (benchmarks/benchmark_count_error_fit.md measures the size on a count-level mixture of two 1000 Genomes samples of different populations), and a relative is fitted at about half its share by the unrelated model, which contamination_relative reports under the other hypothesis. The intervals come from the blocks, not from the likelihood, because sites are linked.

Memory: DuckDB builds the histogram and can spill it; the native buffer of a sample holds 16 bytes per cell, at most max_cells cells. max_cell_bytes bounds the native histogram memory of the process at once, each limit with an explicit error: the buffers of the aggregate (16 bytes per cell, up to three times that while a buffer grows) and the working memory of the running fits (24 bytes per cell, 88 bytes per distinct (depth, alt count) pair and 24 bytes per depth up to the deepest). The BLOB that carries a sample's cells from the aggregate to its fit is DuckDB's memory, under its memory_limit. When a query fails while the histogram is being built with 5 or more threads, DuckDB 1.5.x does not destroy the states of the parallel aggregate (https://github.com/duckdb/duckdb/issues/26871); their buffers stay charged against max_cell_bytes until the process ends. On DuckDB 1.5.x, an out-of-memory error while the histogram states are held can end the process instead of returning the error (https://github.com/duckdb/duckdb/issues/26441); raise memory_limit or fit fewer samples per query. A sample has at most 65,536 blocks and 1,024 frequency bins. Depth is capped by max_depth, which bounds the cells; deeper sites are left out. One fit runs on one thread; the fits of a query's samples run in the threads that hold their rows, so more threads fit more samples at once (the report measures the gain at four threads).

SELECT sample, seq_error, contamination, contamination_sd, status FROM duckhts_count_error_fit('site_counts');

Estimate a conditional second-genome share from calibrated low-pass read and discordant-pair moments.

Signature:

duckhts_count_mixture_lowpass(counts_table, seq_error, artefact_weight, homozygosity_excess, block_bases := 10000000, max_depth := 8, min_blocks := 30, max_blocks := 512, max_candidates := 16, max_state_bytes := 268435456)

Returns:

table(sample VARCHAR, candidate VARCHAR, relationship_model VARCHAR, share DOUBLE, share_lower DOUBLE, share_upper DOUBLE, share_se DOUBLE, singles_share DOUBLE, singles_se DOUBLE, pair_residual DOUBLE, pair_residual_se DOUBLE, pair_denominator DOUBLE, discordance_ratio DOUBLE, seq_error DOUBLE, artefact_weight DOUBLE, homozygosity_excess DOUBLE, moment_score DOUBLE, max_block_leverage DOUBLE, sites BIGINT, usable_sites BIGINT, zero_depth_sites BIGINT, singleton_sites BIGINT, reads BIGINT, other_reads BIGINT, positive_pair_sites BIGINT, excluded_sites BIGINT, concordant_ref_pairs BIGINT, discordant_pairs BIGINT, concordant_alt_pairs BIGINT, unavailable_sites BIGINT, blocks BIGINT, pair_blocks BIGINT, score_rank BIGINT, jackknife_failed_refits BIGINT, at_bound BOOLEAN, matched_frequencies BOOLEAN, interval_hull BOOLEAN, status VARCHAR, relation_diagnostic VARCHAR, method VARCHAR)

counts_table contains sample_id, candidate, chrom, one-based whole pos, ref_count, alt_count, other_count, sample_af and candidate_af. The frequencies describe the allele counted by alt_count. Sample, candidate and chromosome identifiers are nonempty. One row per sample/candidate/chrom/pos is required; counts and sample_af must agree across candidates at a site. Counts are nonnegative whole BIGINT values. Frequencies are NULL or finite in [0,1]; a site missing any candidate row/frequency or any count is unavailable to every candidate of that sample. Zero-depth and singleton sites remain in QC. Depth is ref_count + alt_count; deeper-than-max_depth sites are excluded, not truncated. other_count contributes to QC only. Calibration is required and fixed: seq_error in [0,0.5), artefact_weight in [0,1), homozygosity_excess in [0,1]. Positions k * block_bases + 1 through (k + 1) * block_bases share a chromosome-local block.

One row per sample/candidate/relationship_model; unrelated is always present. An ibd1 row requires exact numerical candidate_af = sample_af at every site, including zero-depth, unavailable-count and excluded sites; NULL frequencies do not match. An empty input returns no rows; a sample with zero usable reads retains a no_sites row with NULL share. QC pair counts are choose(ref_count,2), ref_countalt_count and choose(alt_count,2), over usable, non-excluded sites only. pair_denominator is SUM(choose(depth,2)sample_af*(1-sample_af)); discordance_ratio uses that baseline-one denominator, not twice it, and is NULL if it is zero. blocks counts informative blocks (pair-bearing for matched frequencies, read-bearing otherwise); pair_blocks counts pair-bearing blocks. score_rank is the retained block covariance rank, using a relative eigenvalue cutoff of 1e-10. moment_score is the squared weighted moment discrepancy. max_block_leverage is the largest read or pair exposure fraction; jackknife_failed_refits counts unusable deletions. at_bound flags a share near 0 or 0.5. interval_hull flags limits spanning disconnected accepted share sets. method is duckhts_count_mixture_lowpass v1.

The share is constrained to [0,0.5]. First moments use the basis (1,f,g); discordant-pair moments use (1,f,g,ff,fg,gg). A two-pass generalized moment fit builds a rank-revealing weight matrix from centered block residuals at an exposure-normalized first fit. Genotypes are HWE with fixed sample homozygosity excess; unrelated draws a second independent HWE genotype at candidate_af. The matched IBD1 hypothesis draws one transmitted sample allele and one independent allele at sample_af. Sequencing error maps q to seq_error+(1-2seq_error)q. The uniform alt-count artefact has mean allele fraction 1/2 and discordant-pair probability 1/3 at every depth. singles_share is the unbounded calibrated read-moment regression diagnostic, with block-jackknife singles_se. pair_residual is observed minus unrelated-predicted discordant-pair probability at singles_share, with its deletion standard error. The whole conditional share is refitted after deleting each read-bearing block with calibration fixed; share_se is the delete-one-block jackknife standard error, not the spread of independently fitted block shares. Limits invert the whole-fit score, using jackknife variance to calibrate interior curvature and the block-covariance chi-square(1) score cutoff at a share bound. The 95 percent limits are clipped to [0,0.5]. Fewer than min_blocks informative blocks yields NULL limits and insufficient_blocks while retaining an identified point share; default 30, at least 2. Lowering the floor is for small fixtures and does not establish interval coverage. Failed refits, zero covariance rank or gross moment mismatch yield NULL limits. The gross-mismatch screen uses df + 2sqrt(dflog(1000)) + 2log(1000), df = score_rank-1, as a conservative chi-square tail threshold.

The linked phased chr20 campaign uses 1 Mb blocks and 50 replicates per stratum, not the genome-wide release target of 200. It does not establish reliable detection of 1-2 percent shares at 0.1x; many point-accuracy thresholds fail even at larger shares. Interval coverage alone is not evidence of adequate point precision. Matched-frequency IBD1 fits can be uninformative, and count thinning retains single-source calibration offsets and missing intervals. The autosomal campaign evaluates 90 strata with 200 replicates each, including both error calibrations; 75 strata are primary. At a 5 percent share, unrelated second genomes give a positive interval lower limit in 199-200 of 200 replicates at 0.1x, versus 160-165 of 200 for IBD1; at 0.4x all tested relationships give 200 of 200 positive limits. These call rates do not establish a calibrated detection limit: only 30 of 75 primary strata meet every declared criterion, and 5 of 15 primary null strata fail the false-positive upper-bound requirement. Frequency anchors include the simulated genomes. No universal share detection limit or genome-wide accuracy guarantee is established. See benchmarks/benchmark_count_mixture_lowpass.md and benchmarks/benchmark_count_mixture_lowpass_genome.md for all strata, denominators and failures.

The estimate is conditional on both populations' supplied frequencies and on external error, artefact and homozygosity calibration. It is not an unconditional biological fraction and is not calibrated for arbitrary low-input DNA or cell-free DNA. Matched-frequency unrelated and IBD1 models can have identical first and second moments: relation_diagnostic is unresolved rather than selecting a relationship by score. With population contrast and enough blocks, the block-standardized pair residual is compatible_unrelated within 1.96 standard errors, shared_allele_signal below that range, or model_mismatch above it; lack of information gives uninformative. A shared-allele signal is a diagnostic, not a relationship classification. Status is ok_conditional, no_sites, no_pairs, no_population_contrast, insufficient_blocks, nonidentified, uninformative or model_mismatch. Flat moments give a NULL share; no_population_contrast marks read-only data without frequency contrast, and nonidentified also marks covariance-rank failure. Interval coverage at the block floor is not guaranteed.

block_bases is positive whole, max_depth is 1..65535, min_blocks is 2..2147483647, max_blocks is 1..4096 (default 512), max_candidates is 1..64 per sample (default 16), and max_state_bytes is positive whole (default 256 MiB). Native state retains only fixed-size block moments, never sites or a depth histogram. Its aggregate buffers and fit copies use the shared native count-state byte budget; exceeding the active query limit returns an error. DuckDB owns exact binary per-block sums, SQL validation/joins/sorts/spill and the state BLOB, under its memory_limit, not max_state_bytes. Each non-negative moment accumulator uses 152 bytes independently of site count and rounds once to binary64; SQL validation and joins still require input-sized state. Blocks are sorted before fitting so input order and thread count do not change the result. One fit uses fixed-size numerical workspace and one thread; delete-one-block refits take quadratic work in the number of blocks. DuckDB 1.5.x may retain parallel aggregate states after an aggregate error (https://github.com/duckdb/duckdb/issues/26871); such retained buffers remain charged until process exit. This separate estimator does not change duckhts_count_error_fit.

SELECT sample, candidate, relationship_model, share, share_lower, share_upper, relation_diagnostic, status FROM duckhts_count_mixture_lowpass('lowpass_counts', 0.01, 0.02, 0.0);

Derive a stable SHA-256 identity for an ordered biallelic sample-fingerprinting panel.

Signature:

duckhts_somalier_panel_sha256(panel_table)

Returns:

VARCHAR

panel_table has assembly, dense zero-based site_index, region, positive one-based position, allele_a and allele_b. Assembly and region are each limited to 1,024 bytes before hashing. One nonempty assembly, unique physical region/position and uppercase single-base A/C/G/T alleles in lexical A < B order are required. The pinned Somalier v0.3.4 X/Y aliases are sex-chromosome sites; every other region is autosomal, since other aliases cannot be classified from a region string. Every autosomal site_index must precede every X/Y site_index. The digest commits to this domain (panel version 3: biallelic SNVs, autosomal then sex-chromosome sites) and every ordered site, independent of physical row order. Panels digested before this version have different identities.

SELECT duckhts_somalier_panel_sha256('fingerprint_panel');

Derive a stable identity for panel-aligned population-B allele frequencies.

Signature:

duckhts_somalier_frequency_sha256(frequency_table, panel_table)

Returns:

VARCHAR

frequency_table must cover every autosomal panel site exactly once (X/Y rows are ignored) with the same assembly, ordinal, coordinate and A/B orientation plus a finite population_b_af in [0,1]. The digest commits to the panel identity and every ordered frequency value.

SELECT duckhts_somalier_frequency_sha256('population_frequencies', 'fingerprint_panel');

Classify one measured A/B/other count tuple for Somalier-derived autosomal relatedness.

Signature:

duckhts_somalier_classify(a, b, other, min_depth, min_het_balance, hom_balance_cutoff)

Returns:

STRUCT(genotype TINYINT, middling BOOLEAN, unavailable BOOLEAN)

a, b and other are all measured or all NULL. Three zeros are measured zero-depth evidence; three NULL values are unavailable. Genotype is -1 unknown, 0 homozygous A, 1 heterozygous or 2 homozygous B.

Alleles must already use the panel's ordered A/B orientation. Relatedness classification preserves the pinned Somalier v0.3.4 10% other-read filter; contamination uses separate stricter eligibility.

SELECT duckhts_somalier_classify(20, 20, 0, 7, 0.3, 0.01);

Build one panel-verified packed relatedness sketch per sample from typed count evidence.

Signature:

duckhts_somalier_prepare_sketches(evidence_table, panel_table, min_depth, min_het_balance, hom_balance_cutoff, max_sites := 1000000)

Returns:

table(sketch STRUCT)

evidence_table contains sample_id, the six panel identity columns, and nullable a, b and other counts. Each sample must contain every panel ordinal exactly once. Panel and evidence rows on the X/Y aliases are ignored: this is an autosomal statistic, and results equal those on the autosomal-only panel. The reported panel identity is the whole panel's. Count channels already follow the panel's canonical lexical A/B order; changed coordinates or orientation error. Counts are never inferred from GT or DP. All-NULL tuples are unavailable, all-zero tuples are measured zero depth, and partial NULL tuples error.

The returned struct contains identities, classification settings, counters, a raw-count receipt, content integrity fields and three UBIGINT[] masks. The receipt binds every ordinal, availability state and A/B/other tuple even when changed counts retain the same genotype. It is an ordinary typed value suitable for Parquet, not a serialized native object. max_sites bounds each prepared sketch; sample_id and assembly are each limited to 1,024 bytes.

CREATE TABLE sample_sketches AS SELECT * FROM duckhts_somalier_prepare_sketches('allele_counts', 'fingerprint_panel', 7, 0.3, 0.01);

Verify persisted relatedness sketches against their retained raw count evidence.

Signature:

duckhts_somalier_verify_sketches(evidence_table, panel_table, sketches_table, max_sites := 1000000)

Returns:

BOOLEAN

Checks persisted classification settings before using them as rebuild parameters. With valid settings, it rebuilds every sample through the panel and count-validation path and compares the complete typed sketch. It returns false for invalid retained settings, changed A/B/other counts, changed availability, altered masks or receipts, and missing, extra or duplicate sample sketches. Invalid panel/evidence geometry and incomplete or duplicate site ordinals error.

evidence_table must contain exactly the samples represented by sketches_table. max_sites is a positive panel-site limit at most 100,000,000. Receipts detect accidental divergence between persisted evidence and sketches; they are not an authenticity mechanism.

SELECT duckhts_somalier_verify_sketches('allele_counts', 'fingerprint_panel', 'sample_sketches');

Compute fused Somalier-derived relatedness and concordance statistics for two prepared sketches.

Signature:

duckhts_somalier_relatedness(sketch_a, sketch_b, max_sites)

Returns:

STRUCT

The struct retains sample and panel identities, method version, status, jointly-called count, IBS0/IBS2, shared heterozygotes, heterozygote and homozygote denominators, middling/unavailable counters, relatedness and named concordance values. relatedness is 2(shared_hets - 2 IBS0)/max(1,het_ab). inferred_hom_concordance is matching_hom_count/max(1,min(callable_hom_count_a,callable_hom_count_b)). raw_hom_b_concordance is (shared_hom_b - 2 IBS0)/max(1,min(hom_b_count_a,hom_b_count_b)). adjusted_concordance applies pinned v0.3.4's middling and low-homozygote penalties and upper-range transform. Floating results are NULL when no site is jointly callable.

Both sketches must have identical panel digests, site counts and classification settings. The contents of both sketches are checked for every call, so comparing every pair of a cohort repeats that check for each pair; use duckhts_somalier_relatedness_all_pairs for all pairs. Mask words are borrowed directly; comparison allocates no pair-sized workspace. SQL chooses the requested pair relation and output ordering.

SELECT unnest(duckhts_somalier_relatedness(a.sketch, b.sketch, 1000000)) FROM selected_pairs p JOIN sample_sketches a ON a.sketch.sample_id = p.sample_a JOIN sample_sketches b ON b.sketch.sample_id = p.sample_b;

Compute Somalier-derived relatedness and concordance statistics for every unordered sample pair of one sketch relation.

Signature:

duckhts_somalier_relatedness_all_pairs(sketches_table, max_sites := 1000000, max_pairs := 50000000)

Returns:

table

sketches_table names a relation with one non-NULL sketch column per sample, as duckhts_somalier_prepare_sketches returns. Two sketches of one sample, a NULL sketch, a sketch with more sites than max_sites, or a sketch whose masks do not match its stored content digest is an error. All sketches must share assembly, panel digest, site count and classification settings; a mixed relation is an error. These checks run once over the relation, so they hold for a query that reads no pair column, for example count(*).

One row per unordered pair with sample_a < sample_b, with the fields of duckhts_somalier_relatedness as columns. The rows equal those of duckhts_somalier_relatedness over the same pairs. No row order is guaranteed.

Each sketch is content-checked once and the pairs are compared from the checked sketches, so the content checks grow with the number of samples and the comparisons with pairs times sites. The pair indexes are split into at most 65,536 blocks of consecutive indexes. DuckDB scans the blocks in chunks of up to 2,048, one chunk per task, so at most 32 of its threads compare pairs. The rows do not depend on the thread count.

The output has N(N-1)/2 rows for N samples. max_pairs must be a positive, finite whole number; a relation with more pairs than max_pairs is an error before any pair is compared. The default, 50,000,000, is 400 times the pairs of 500 samples. Memory model for N sketches of S sites, each with three masks of ceil(S/64) 64-bit words: the checked sketches form one list in sample_id order, which the query holds in a few shared copies. Each thread copies the list once for each chunk of blocks it scans and reads both sketches of every pair of that chunk in place from its copy. The workspace of a thread is therefore one copy of the list plus fixed 2,048-row vectors of pair indexes and results: it grows with N times S, not with the pairs, and the pair indexes are never collected. See benchmarks/benchmark_somalier_cohort_scaling.md for measured time and memory.

SELECT * FROM duckhts_somalier_relatedness_all_pairs('sample_sketches');

Verify a typed relatedness result against its two sealed sketches.

Signature:

duckhts_somalier_verify_relatedness(pair_result, sketch_a, sketch_b, max_sites)

Returns:

BOOLEAN

Returns false when a sketch digest, sample or panel identity, method/status, site denominator, integer statistic, floating statistic or status-dependent NULL value is inconsistent. The check recomputes the pair metrics from borrowed mask words without a per-pair allocation. It detects accidental corruption, not malicious rewriting of both sketches and their integrity fields.

max_sites is a positive panel-site limit at most 100,000,000. Invalid or NULL typed inputs return false.

SELECT duckhts_somalier_verify_relatedness(r.pair, a.sketch, b.sketch, 1000000) FROM pair_results r JOIN sample_sketches a ON a.sketch.sample_id = r.pair.sample_a JOIN sample_sketches b ON b.sketch.sample_id = r.pair.sample_b;

Estimate per-sample contamination with a bounded Somalier-derived CHARR reduction.

Signature:

duckhts_somalier_charr(evidence_table, panel_table, frequency_table, min_depth := 15, max_depth := 1000000, hom_minor_rate := 0.12, hom_tail_alpha := 0.002, max_threshold_work := 16000000, max_sites := 1000000)

Returns:

table(contamination STRUCT)

Count evidence and population_b_af must match the panel's full ordered site identity. Panel and evidence rows on the X/Y aliases are ignored: this is an autosomal statistic, and results equal those on the autosomal-only panel. The reported panel identity is the whole panel's. CHARR uses measured counts, its own homozygous-like binomial eligibility and the pinned 4% other-read filter; relatedness genotype masks are insufficient.

The struct retains sample, panel and frequency identities, method and numerical status, observed/unavailable/usable and homozygous-site denominators, estimate and every filter/limit. No usable evidence has status no_evidence and a NULL estimate, distinct from measured zero contamination.

A+B depth must not exceed 1,000,000. Distinct measured depths are certified once per call; max_threshold_work bounds their cumulative exact recurrence and continued-fraction steps and is at most 100,000,000. Exhaustion errors without publishing a partial result. sample_id, assembly and panel region are each limited to 1,024 bytes. Input order and parallel aggregate reduction order do not change the estimate.

Eligibility follows Somalier 0.3.4 CHARR except that DuckHTS computes high-depth binomial tails without upstream's numerical underflow. Very deep sites can therefore have different eligibility; a pinned upstream counterexample is retained in the conformance tests.

SELECT unnest(contamination) FROM duckhts_somalier_charr('allele_counts', 'fingerprint_panel', 'population_frequencies');

Estimate directional contamination for explicitly selected receiver/anchor sample pairs.

Signature:

duckhts_somalier_matched_contamination(evidence_table, panel_table, frequency_table, pairs_table, min_depth := 15, max_depth := 1000000, hom_minor_rate := 0.05, hom_tail_alpha := 0.001, error_rate := 0.002, min_probability := 1e-10, min_prior_frequency := 1e-6, alpha_min := 0, alpha_max := 1, grid_step := 0.01, refine_tolerance := 1e-10, max_evaluations := 4096, max_threshold_work := 16000000, max_sites := 1000000)

Returns:

table(contamination STRUCT)

pairs_table contains distinct receiver_id and anchor_id rows. Panel and evidence rows on the X/Y aliases are ignored: this is an autosomal statistic, and results equal those on the autosomal-only panel. The reported panel identity is the whole panel's. The anchor supplies the receiver's expected uncontaminated homozygous genotype; it is not assumed to identify the contaminating donor. Reversing a pair is a different fit.

The struct retains ordered sample, panel and frequency identities, method/status, observed/unavailable/usable denominators, alpha, evaluation count and every filter/search limit. No usable evidence returns NULL alpha. relative_log_likelihood omits alpha-independent binomial coefficients and is comparable only across alpha values for the same observations.

The query certifies each distinct measured depth once, then prepares one bounded site profile per distinct selected sample and one panel-aligned frequency profile before joining requested ordered pairs. The pair scalar borrows those DuckDB-owned lists and allocates no per-pair workspace. Panel and evidence cardinalities are checked before profile-list construction; max_sites is a per-call panel limit.

Panel assembly and region are each limited to 1,024 bytes. Persisted profiles separately limit sample_id and assembly to 1,024 bytes. A+B depth must not exceed max_depth. max_threshold_work bounds cumulative exact binomial certification steps and is at most 100,000,000; exhaustion errors without publishing profiles. max_sites is at most 100,000,000, and max_evaluations must fit the declared search workspace. The site profiles are ordered list aggregates, one per sample: on DuckDB 1.5.x, an out-of-memory error while they are built can end the process instead of returning the error (https://github.com/duckdb/duckdb/issues/26441), and a query that fails while they are built with 5 or more threads keeps their memory until the process ends (https://github.com/duckdb/duckdb/issues/26871); raise memory_limit or pass fewer samples per query.

DuckHTS searches a full 0.01 grid and refines a feasible local optimum rather than reproducing Somalier v0.3.4's fixed coarse/high-refinement sequence. The retained two-site witness fits alpha about 0.39759 versus upstream 0.440983. Both likelihoods and settings are retained in the differential test; results are not guaranteed bitwise identical to the CLI.

SELECT unnest(contamination) FROM duckhts_somalier_matched_contamination('allele_counts', 'fingerprint_panel', 'population_frequencies', 'receiver_anchor_pairs');

Report X/Y dosage evidence and a review-only XX, XY or ambiguous call per sample from panel counts.

Signature:

duckhts_somalier_sex(counts_table, panel_table, 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 := 'sample')

Returns:

table

counts_table is the sample-by-panel relation from duckhts_somalier_vcf_counts() or duckhts_somalier_bam_counts() over a panel that includes X/Y sites; each sample needs one row for every panel site with the panel's geometry and A/B orientation. Coverage is checked per sample from the row count plus the sum and XOR of a 64-bit hash of site_index against the panel's, because a per-sample DISTINCT set would hold samples x sites entries: a missing or repeated ordinal changes the count or the hashes and is rejected, while a different multiset of ordinals with the same count passes only if both hashes collide. Autosomal sites give the depth normaliser, X and Y sites the metrics. Without a usable X site the panel is still accepted and the status says so.

Follows Somalier v0.3.4 relate.nim. autosomal_depth_mean is the mean A+B depth of autosomal sites with other <= 10% of all reads and A+B >= min_depth. An X or Y site is usable when other <= 4%, A+B >= min_sex_site_depth and its B balance is homozygous (below hom_balance_cutoff, above one minus it, or A or B zero) or heterozygous (between min_het_balance and one minus it); x_hom_ref, x_het and x_hom_alt count usable X sites by genotype in the panel's A/B order (allele_a below allele_b), so hom_ref means A/A and hom_alt means B/B, not REF/REF and ALT/ALT of any sites VCF; Somalier 0.3.5 and later reorient its reported hom-ref and hom-alt counts through --sites, which this function does not do. x_depth_ratio and y_depth_ratio are 2 * mean A+B depth over usable sites / autosomal_depth_mean, so a diploid X is about 2 and a single X about 1. x_het_hom_alt_ratio is x_het / x_hom_alt (Infinity when there is no hom-alt site and some het site, NULL when both are zero), the statistic Somalier thresholds. Ratios are NULL without usable sites or autosomal depth; the integer counts stay exact zeros. Floating means can differ from Somalier's running mean in the last bits.

inferred_sex is XY when there are at least min_usable_x_sites usable X sites (Somalier: more than 10) and x_het_hom_alt_ratio < xy_max_het_ratio, XX when it is > xx_min_het_ratio, otherwise ambiguous. y_signal is present when y_depth_ratio > y_signal_min and absent otherwise, NULL without usable Y sites. An XX dosage with a present Y signal is reported as ambiguous with the review flag apparent_y_with_xx_dosage, and an XY dosage with an absent Y signal keeps XY with apparent_y_loss, as Somalier notes. status is ok, no_usable_x_sites, insufficient_x_sites or no_autosomal_depth.

Somalier applies its Y check only when the cohort has Y depth (more than five samples, or more than 10%, with usable Y sites), so a sample's call depends on who else is in the batch. The default y_gate := 'sample' applies the check whenever the sample itself has a Y signal; y_gate := 'cohort' reproduces Somalier's gate, and cohort_has_y reports it. Somalier also leaves the sex unset for samples that fail its autosomal quality gate; DuckHTS reports the X dosage regardless and reports no_autosomal_depth when the normaliser is missing.

The call is evidence of X/Y dosage for review, not a diagnosis or a legal determination of sex; aneuploidy, mosaicism, poor coverage and sample swaps all change it. Results also carry the whole-panel digest, panel denominators and every threshold used.

See benchmarks/benchmark_somalier_sex.md for measured scaling. One pass over the counts relation: state is one aggregate group of a few numbers per sample plus the panel; no product of samples and sites is retained beyond the counts relation itself.

SELECT sample_id, inferred_sex, status, x_depth_ratio, y_depth_ratio FROM duckhts_somalier_sex('allele_counts', 'fingerprint_panel');

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

Signature:

detect_quality_encoding(path, max_records := 10000)

Returns:

table

SELECT * FROM detect_quality_encoding('reads.fq.gz');

Read GFF annotations with optional parsed attributes, strict GFF3 validation and indexed region selection.

Signature:

read_gff(path, header_names := NULL, header := FALSE, column_types := NULL, auto_detect := FALSE, attributes_map := FALSE, attributes_list := FALSE, attributes_pairs := FALSE, attributes := []::VARCHAR[], strict := FALSE, region := NULL, index_path := NULL, scan_mode := 'auto')

Returns:

table

attributes_map leaves out an attribute whose key is not valid UTF-8, since a MAP key cannot be NULL; a value that is not valid UTF-8 is NULL.

A requested key absent from every feature yields NULL on every row, exactly as attributes_map[key] does. GFF declares no attribute vocabulary, so binding does not scan the file to validate key presence; that would require a full extra pass. Assert presence with an aggregate, for example: SELECT CASE WHEN count(ID) > 0 THEN true ELSE error('ID is all NULL') END FROM read_gff('annotations.gff', attributes := ['ID']);

attributes := ['Parent', 'ID'] appends VARCHAR columns named for the requested keys. Each value equals attributes_map[key], including NULL for absent keys; repeated keys use the first value and GFF3 percent encoding is retained. Key lookup is case-sensitive, as GFF3 tags are: attributes := ['id'] on a file that uses ID returns a column of NULLs. Keys must be nonempty and distinct from fixed and optional attribute column names; because column names ignore case, two keys that differ only in case are rejected. Only projected keys are parsed during scanning.

Comma-separated indexed regions emit each row once across overlaps. scan_mode='sequential' streams/counts instead of using indexed count paths and rejects region. NULL/empty region means no filter; empty comma-separated items and malformed known-contig intervals error. Unknown contigs follow HTSlib's skip policy.

SELECT seqname, feature, start, "end" FROM read_gff('gff_file.gff.gz') LIMIT 5;

Read GTF annotations with optional parsed attributes and indexed region selection.

Signature:

read_gtf(path, header_names := NULL, header := FALSE, column_types := NULL, auto_detect := FALSE, attributes_map := FALSE, attributes_list := FALSE, attributes_pairs := FALSE, attributes := []::VARCHAR[], region := NULL, index_path := NULL, scan_mode := 'auto')

Returns:

table

attributes_map leaves out an attribute whose key is not valid UTF-8, since a MAP key cannot be NULL; a value that is not valid UTF-8 is NULL.

A requested key absent from every feature yields NULL on every row, exactly as attributes_map[key] does. GTF declares no attribute vocabulary, so binding does not scan the file to validate key presence; that would require a full extra pass. Assert presence with an aggregate, for example: SELECT CASE WHEN count(gene_id) > 0 THEN true ELSE error('gene_id is all NULL') END FROM read_gtf('annotations.gtf', attributes := ['gene_id']);

attributes := ['gene_id', 'transcript_id'] appends VARCHAR columns named for the requested keys. Each value equals attributes_map[key], including NULL for absent keys; repeated keys use the first value. Key lookup is case-sensitive: attributes := ['Gene_ID'] on a file that uses gene_id returns a column of NULLs. Keys must be nonempty and distinct from fixed and optional attribute column names; because column names ignore case, two keys that differ only in case are rejected. Only projected keys are parsed during scanning.

SELECT seqname, feature, start, "end" FROM read_gtf('annotations.gtf.gz') LIMIT 5;

Read GenBank flat-file features in read_gff's column shape, with optional parsed qualifier MAP.

Signature:

read_genbank(path, attributes_map := FALSE, attributes := []::VARCHAR[])

Returns:

table

attributes_map leaves out a qualifier whose key is not valid UTF-8, since a MAP key cannot be NULL; a value that is not valid UTF-8 is NULL.

seqname is VERSION, else ACCESSION, else the LOCUS name; a segment on a remote accession (ACC.1:5..40) is reported under that accession. source is 'GenBank'. join()/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 (absent means 1). A span whose end precedes its start on a circular record wraps the origin into two segments, and a between site n^m is the zero-length site at n. The record-level source feature is dropped, and /translation is omitted as redundant with ORIGIN.

Synthesized GFF3 keys ID, Name and Parent accompany the original qualifiers. ID is gene- for genes and - otherwise, with n the feature's 0-based position in the file. Parent links a feature to the gene sharing its /locus_tag (or /gene) anywhere in the record. Repeated qualifiers become one key with comma-joined values, valueless qualifiers read 'true', and ; = & , % are percent-encoded. attributes_map := TRUE adds the same pairs as a MAP.

A requested key absent from every feature yields NULL on every row, exactly as attributes_map[key] does. GenBank declares no qualifier vocabulary, so binding does not scan the file to validate key presence; that would require a full extra pass. Assert presence with an aggregate, for example: SELECT CASE WHEN count(product) > 0 THEN true ELSE error('product is all NULL') END FROM read_genbank('phix174.gb', attributes := ['product']);

attributes := ['gene', 'product'] appends VARCHAR columns named for the requested keys, after attributes_map when that is requested. Each value equals attributes_map[key] byte for byte: repeated qualifiers are comma-joined, a valueless qualifier reads 'true', values stay percent-encoded, the synthesized ID, Name and Parent are addressable, and an absent key is NULL. Key lookup is case-sensitive, like qualifier keys themselves: attributes := ['ec_number'] on features that carry /EC_number returns a column of NULLs. Keys must be nonempty and distinct from the fixed columns and attributes_map when it is requested; because column names ignore case, two keys that differ only in case are rejected. Only projected keys are computed during scanning.

Records stream one at a time, so memory follows the largest record. A record without a terminating //, a FEATURES table with no sequence section, a malformed location, an unsupported location form (one-of, gap, bond, nested join/order, a.b), an origin-spanning span on a linear record, or a /codon_start outside 1..3 is an error naming the feature and line. Table layout and these rules follow BioPython's GenBank scanner, and test/scripts/genbank_oracle_test.py diffs the reader against it.

SELECT seqname, feature, start, "end" FROM read_genbank('phix174.gb') LIMIT 5;

SELECT feature, locus_tag, product FROM read_genbank('phix174.gb', attributes := ['locus_tag', 'product']) WHERE feature = 'CDS';

Write the ORIGIN sequence of each GenBank record as FASTA and return success, output_path and records_written.

Signature:

genbank_to_fasta(path, output_path := NULL, line_width := 70, overwrite := FALSE)

Returns:

table

Records are written under the same name read_genbank reports as seqname, so feature coordinates land on the contig of that name; the DEFINITION follows without its trailing period, as in NCBI's FASTA export. A segment on a remote accession is not written.

output_path defaults to path with .fa appended. 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; with overwrite := FALSE an existing output is an error. Records without an ORIGIN block are skipped; zero written records is an error. line_width must be at least 1.

SELECT * FROM genbank_to_fasta('phix174.gb', output_path := 'phix174.fa');

Read tabix-indexed text with optional header handling, inferred types and region selection.

Signature:

read_tabix(path, header_names := NULL, header := FALSE, column_types := NULL, auto_detect := FALSE, region := NULL, index_path := NULL, scan_mode := 'auto')

Returns:

table

SELECT * FROM read_tabix('meta_tabix.tsv.gz') LIMIT 5;

Build a FASTA index (.fai) and return a single row with columns success (BOOLEAN) and index_path (VARCHAR).

Signature:

fasta_index(path, index_path := NULL)

Returns:

table

SELECT * FROM fasta_index('ce.fa');

Compress a plain file to BGZF and return the created output path and byte counts.

Signature:

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

Returns:

table

Compression runs during query execution. Binding, PREPARE, and EXPLAIN without ANALYZE do not create output or remove input files. Each execution enforces overwrite; keep=FALSE removes the input only after successful completion.

SELECT * FROM bgzip('regions.bed');

Decompress a BGZF-compressed file and return the created output path and byte counts.

Signature:

bgunzip(path, output_path := NULL, threads := 4, keep := TRUE, overwrite := FALSE)

Returns:

table

Decompression runs during query execution. Binding, PREPARE, and EXPLAIN without ANALYZE do not create output or remove input files. Each execution enforces overwrite; keep=FALSE removes the input only after successful completion.

SELECT * FROM bgunzip('regions.bed.gz');

Build a BAM or CRAM index and report the written index path and format.

Signature:

bam_index(path, index_path := NULL, min_shift := 0, threads := 4)

Returns:

table

SELECT * FROM bam_index('range.bam');

Build a TBI or CSI index for a VCF or BCF file and report the written index path and format.

Signature:

bcf_index(path, index_path := NULL, min_shift := NULL, threads := 4)

Returns:

table

SELECT * FROM bcf_index('formatcols.vcf.gz');

Build a tabix index for a BGZF-compressed text file using a preset or explicit coordinate columns.

Signature:

tabix_index(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)

Returns:

table

SELECT * FROM tabix_index('gff_file.gff.gz', preset := 'gff');

Compose projected BAM/CRAM reads, flag predicates and packed-CIGAR lengths into riker 0.4.1 insert-size counts.

Signature:

duckhts_bam_insert_size_counts_sql(path, reference := NULL, include_duplicates := FALSE, decompression_threads := 0, scan_mode := 'auto', index_path := NULL)

Returns:

table(pair_orientation VARCHAR, insert_size BIGINT, read_pairs UBIGINT)

Fulcrum Genomics riker 0.4.1 (MIT), tag 40d5e3260c5b1b9d210c28f7c90a75387013b4ce, is the oracle. Its executable rules supersede the design note's draft rules; Picard is context only.

riker src/commands/isize.rs:197-218: paired, read and mate mapped, primary (not secondary or supplementary), not first-segment or QC-fail, nonzero TLEN, same reference, and no duplicate flag unless include_duplicates is true. No MAPQ, proper-pair or QNAME predicate; a second-segment flag is not required. Absolute signed-32-bit TLEN includes INT32_MIN as 2147483648.

riker src/sam/pair_orientation.rs:30-59: equal strands are TANDEM. For reverse reads compare mate alignment start with inclusive CIGAR reference end; for forward reads compare POS with POS + TLEN (no minus one). Strictly earlier forward end is FR; equality is RF. Missing CIGAR spans one reference base (noodles alignment_end, called at :49).

Unordered sparse rows; absent bins are omitted. read_bam projects only FLAG, RNAME, RNEXT, POS, PNEXT, TLEN and binary CIGAR. DuckDB's threads setting controls SQL processing; scan_mode is auto or sequential as in read_bam. decompression_threads controls htslib workers per handle. No SEQ/QUAL/tags/QNAME columns are emitted, but CRAM CIGAR projection currently requests SEQ/QUAL and AUX internally for long-CIGAR CG restoration (src/bam_reader.c bam_cram_required_fields).

duckhts_bam_insert_size_counts implements the same counts with native per-worker sparse histograms. benchmarks/benchmark_insert_size_counts.md records both arms against riker; native takes the unsuffixed name by measurement. SQL misses every declared 1x target on that workload. Both arms remain public.

CREATE TABLE insert_counts AS SELECT * FROM duckhts_bam_insert_size_counts_sql('sample.bam');

SELECT * FROM duckhts_insert_size_metrics('insert_counts', sample_name := 'S1');

Reduce BAM/CRAM records to sparse insert-size counts using a native riker 0.4.1 kernel.

Signature:

duckhts_bam_insert_size_counts(path, reference := NULL, include_duplicates := FALSE, threads := 1, decompression_threads := 0, scan_mode := 'auto', index_path := NULL, reference_index_path := NULL)

Returns:

table(pair_orientation VARCHAR, insert_size BIGINT, read_pairs UBIGINT)

Fulcrum Genomics riker 0.4.1 (MIT), tag 40d5e3260c5b1b9d210c28f7c90a75387013b4ce, is the executable oracle. Its rules supersede the design note's draft first-segment, QC-fail, POS + TLEN - 1 and mixed-case orientation rules. Picard InsertSizeMetricsCollector is context only, not the oracle.

riker src/commands/isize.rs:197-206: paired, read and mate mapped, primary (not secondary or supplementary), not first-segment or QC-fail, nonzero TLEN, and no duplicate flag unless include_duplicates is true. Neither MAPQ (including 0 and 255) nor proper-pair status is a predicate. Rejecting first segments does not require a second-segment flag. No QNAME matching is performed.

riker src/sam/pair_orientation.rs:30-59: require the same reference; equal strands are TANDEM. For a reverse read compare mate alignment start with inclusive CIGAR reference end; for a forward read compare alignment start with alignment start + TLEN (no minus one). Forward end strictly before reverse end is FR; equality is RF. Missing CIGAR spans one reference base, as in noodles alignment_end (riker src/sam/pair_orientation.rs:49). riker src/commands/isize.rs:211-218 counts unsigned absolute TLEN, including INT32_MIN as 2147483648.

Unordered sparse rows; absent bins are omitted. Indexed scans use per-worker file handles and histograms (threads 1..8), with one-million-base task windows and read-start ownership so overlapping reads count once. auto falls back to a single sequential worker if no index exists; indexed requires an index; sequential streams the file explicitly. SEQ, QUAL and tags are not decoded. decompression_threads is htslib workers per handle, independent of histogram workers. Each worker and the merged result are capped at one million distinct orientation/size bins; exceeding the cap is an error. Counts are checked for UBIGINT overflow. The reduction runs during table-function initialization; the stable C API offers no interrupt polling inside that scan. On wasm use threads=1 and decompression_threads=0.

Counts are over the whole file, not grouped by read group or sample. reference supplies the FASTA for CRAM. Sample labels and provenance belong in SQL; no riker text files or plots are written. SAM TLEN is a signed 32-bit value.

duckhts_bam_insert_size_counts_sql composes projected read_bam, flag predicates and packed-CIGAR lengths in SQL. benchmarks/benchmark_insert_size_counts.md records both arms against riker; native takes the unsuffixed name by measurement. Native BAM misses the declared 1x targets at 1, 4 and 8 threads; native CRAM meets all four. Both arms remain public and face the same targets.

CREATE TABLE insert_counts AS SELECT * FROM duckhts_bam_insert_size_counts('sample.bam');

SELECT * FROM duckhts_insert_size_metrics('insert_counts', sample_name := 'S1');

Compute riker 0.4.1 weighted insert-size metrics from a named sparse counts relation.

Signature:

duckhts_insert_size_metrics(counts_table, deviations := 10.0, min_frac := 0.05, sample_name := NULL)

Returns:

table(sample VARCHAR, pair_orientation VARCHAR, read_pairs HUGEINT, median_insert_size DOUBLE, median_absolute_deviation DOUBLE, mean_insert_size DOUBLE, standard_deviation DOUBLE, mode_insert_size BIGINT, min_insert_size BIGINT, max_insert_size BIGINT, histogram_cutoff BIGINT)

Fulcrum Genomics riker 0.4.1 (MIT), tag 40d5e326, is the oracle; Picard is context only. The executable rules supersede the design note's draft rules (see duckhts_bam_insert_size_counts).

counts_table names a table or view with pair_orientation (FR, RF, TANDEM), insert_size (whole 1..UINT32_MAX) and positive whole read_pairs. Duplicate orientation/size rows are summed. Empty input returns no rows. deviations and min_frac must be finite and nonnegative. min_frac greater than 1 yields no metric rows, as in riker; sample_name is caller supplied, not inferred from tags.

riker src/counter.rs:110-115 selects the largest mode on ties; :180-210 interpolates the two middle observations for even totals; :226-240 stores doubled absolute deviations as rounded, saturating u32 keys and halves their median. src/commands/isize.rs:245-277 retains orientations whose count is at least total * min_frac; original counts, median, MAD and min/max use every size. Only mean and population standard deviation trim to sizes <= median + deviations * MAD (src/counter.rs:154-175). No read-pair expansion is used; intermediates scale with distinct bins. Values are numeric; printf('%.2f', value) reproduces riker's two-decimal serialization (src/commands/isize.rs:328-340).

riker src/commands/isize.rs:236-290 uses the floor of the largest median + deviations * MAD among retained orientations. histogram_cutoff repeats that shared value, clamped to UINT32_MAX because larger cutoffs admit the same supported insert sizes. To reproduce the histogram, pivot all orientations' counts at sizes <= this cutoff, not just retained orientations, filling missing counts with zero. If no orientation passes min_frac, riker emits the entire raw histogram; use the counts relation without trimming. Empty input has neither metrics nor histogram rows.

SELECT * FROM duckhts_insert_size_metrics('insert_counts', deviations := 10, min_frac := 0.05, sample_name := 'S1');

Count BAM or CRAM read starts into fixed-width bins. Returns one row per bin across the selected contig span, including zero-count bins, with total, forward, and reverse counts; rmdup := 'streaming' applies the WisecondorX-style larp/larp2 consecutive-position deduplication, rmdup := 'flag' drops SAM duplicate-flagged reads, and stats := 'gc', 'mq', or 'gc,mq' adds per-bin pre/post-filter GC and MAPQ sufficient statistics, including reference GC when reference is provided.

Signature:

bam_bin_counts(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)

Returns:

table

include_unmapped := TRUE appends one synthetic row with chrom = '*' for no-coordinate records, even when no such records are present; start, end and bin_id are NULL. The default FALSE omits this row. SAM-flag filters, MAPQ filtering and flag-based duplicate removal apply to these records; streaming duplicate removal is not applied to them. Reference-GC fields are NULL for this row.

SELECT * FROM bam_bin_counts('fixture_mixed.cram', 5000, reference := 'fixture_ref.fa', rmdup := 'streaming', stats := 'gc,mq');

Correct read counts in fixed bins for GC content with a median read rate per GC stratum, giving each bin an expected count, a depth ratio and whether it is kept.

Signature:

duckhts_depth_bins(bins_table, min_mappability := 0.9, gc_step := 0.01, min_stratum_bins := 10, outlier_quantile := 0.01, max_ratio := 4.0)

Returns:

table(sample_id VARCHAR, chrom VARCHAR, start BIGINT, "end" BIGINT, reads BIGINT, gc DOUBLE, mappability DOUBLE, gc_stratum INTEGER, stratum_bins BIGINT, expected DOUBLE, ratio DOUBLE, log2_ratio DOUBLE, reason VARCHAR, kept BOOLEAN, method VARCHAR)

bins_table is a caller-named table or view with one row per sample and bin: sample_id, chrom, start and end (zero-based, half-open), reads (a whole number of reads, for example count_total of bam_bin_counts) and gc, the G+C share of the A, C, G and T reference bases of the bin from 0 to 1 (for example (num_c + num_g) / (num_a + num_c + num_g + num_t) of fasta_nuc). An optional mappability column from 0 to 1 is used when present. A NULL sample_id or chrom, start < 0 or end <= start, negative or fractional reads, and gc or mappability outside [0, 1] are errors. Restrict contigs (for example to the primary assembly) and remove known problem bins (gaps, centromeres, bins with many N bases, the pseudo-autosomal regions of X where that matters) in SQL first.

A bin's read rate is reads / (end - start). Its GC stratum is floor(gc / gc_step). The training bins of a sample are its autosomal bins that pass the bin rules, have at least one read and a read rate inside the outlier_quantile and 1 - outlier_quantile quantiles of those bins' rates. expected is the bin's length times the median rate of the training bins of its sample and stratum, scaled so that the median ratio of the sample's kept positive-read autosomal bins is 1. Empty bins do not set this scale and retain zero ratios when their bin and stratum rules pass. ratio = reads / expected is then about 1 for two copies, 0.5 for one and 1.5 for three; log2_ratio is NULL for a bin without reads. Autosomal means a contig whose duckhts_contig_key is not X, Y or MT. The strata are fitted on the sample's own autosomal bins, so a change that holds a large share of the training bins of some GC strata is partly absorbed by the correction: correct a whole genome rather than a few chromosomes (on five chromosomes, a planted chr18 gain was recovered at 0.6 to 0.7 of its size; benchmarks/benchmark_depth_fraction.md), or give the estimators expected counts from a reference panel instead.

kept is true when reason is NULL. reason names the first rule a bin fails: no_reads (NULL reads), no_gc, low_mappability (below min_mappability; only when mappability is given), sparse_stratum (fewer than min_stratum_bins training bins in the stratum), no_autosomal_bins (the sample has no autosomal bin to scale by) or high_ratio (ratio above max_ratio, typically collapsed repeats). Only high_ratio looks at a bin's count, so a bin with no reads, such as Y in a genome without Y, stays kept. mappability only filters; it is not a second correction.

gc_step must be in [1/2147483647, 1], so floor(gc / gc_step) fits the INTEGER gc_stratum output even at gc = 1. The median and quantile aggregates hold one value per autosomal bin of each sample; the joins build on samples x GC strata and on samples; the input is materialized once and the output has one row per input row, in no particular order. Bins are not normalised against a reference set; the estimators take one.

SELECT * FROM duckhts_depth_bins('bins') WHERE kept LIMIT 5;

CREATE TABLE bins AS SELECT 'S1' AS sample_id, b.chrom, b.start, b."end", b.count_total AS reads, (n.num_c + n.num_g) / nullif(n.num_a + n.num_c + n.num_g + n.num_t, 0) AS gc FROM bam_bin_counts('sample.bam', 1000000) b JOIN fasta_nuc('reference.fa', bin_width := 1000000) n USING (chrom, start) WHERE n.num_n < 0.05 * n.seq_len;

Estimate the share of cells with a one-copy gain or loss of each named region from corrected bin depth, with an interval from the bins or from a reference set.

Signature:

duckhts_depth_region_fraction(bins_table, regions_table)

Returns:

table(sample_id VARCHAR, region VARCHAR, bins BIGINT, baseline_bins BIGINT, ratio DOUBLE, shift DOUBLE, shift_se DOUBLE, shift_low DOUBLE, shift_high DOUBLE, dispersion DOUBLE, reference_shift DOUBLE, reference_samples BIGINT, status VARCHAR, method VARCHAR)

duckhts_depth_region_fraction gives the interval from the bins; duckhts_depth_region_fraction_ref takes reference_table. bins_table holds corrected bins, such as the output of duckhts_depth_bins: sample_id, chrom, start, end, reads, expected and kept. Only kept bins are used; a kept bin needs reads >= 0 and a finite expected > 0. regions_table has region, chrom, start and end (zero-based, half-open); several rows with one region name form one region, and its overlapping or abutting intervals are merged. A bin belongs to a region when its midpoint lies in one of the region's intervals. A NULL kept counts as not kept. Contig names are compared byte for byte.

With C and E the sums of reads and expected over kept bins, ratio = (C_r / E_r) / (C_b / E_b), where r is the region and b the baseline: the sample's kept autosomal bins outside the region. shift = 2 * (ratio - 1) is the share of cells with a one-copy gain (positive) or loss (negative) of the region relative to two copies, if the baseline is unchanged. A one-copy change in every cell gives 1 or -1; the X of an XY genome gives -1. Other changes in the sample move the baseline and so every shift: a 150 Mb gain in 30% of cells moves the other regions by about -0.014.

Bin counts are quasi-Poisson: the variance of a sum of counts is dispersion times that sum (at least one read), with dispersion the mean squared Pearson residual of the baseline bins, at least 1. shift_se follows by the delta method, treating bins as independent, and the interval is shift +- 1.959964 shift_se. Real libraries have correlated residual bias, so this interval is optimistic for large regions; a reference set measures that bias.

With reference_table, corrected bins of reference samples in the same layout, each reference sample's shift for the region is computed the same way, leaving out the tested sample when its sample_id is also in the reference. shift is the sample's shift minus the reference mean (reference_shift). Its variance is the sample's own bin variance, plus the variance of the reference shifts beyond their mean bin variance, plus the variance of the reference mean; the interval uses the 97.5% quantile of Student's t with reference_samples - 1 degrees of freedom. reference_samples counts the reference samples used.

One row per sample of bins_table and region. status is ok, no_bins (no kept bin in the region), no_baseline (fewer than two kept baseline bins or no baseline reads) or few_reference_samples (fewer than two; shift_se and the interval are NULL). method names the estimator and its version.

The narrow input columns (sample, contig, coordinates, reads, expected, kept) are materialized once, in DuckDB's buffer manager; the per-sample sums and the bin-region join read that copy. The join matches each kept bin to the merged intervals of its contig and keeps the pairs whose interval holds the bin's midpoint; the state per sample and region is a few sums, and the output has one row per sample and region. Depth does not see copy-neutral changes, and a constitutional copy-number variant looks like a change in every cell.

duckhts_depth_region_fraction uses the base inputs. duckhts_depth_region_fraction_ref requires reference_table as an additional positional argument.

SELECT sample_id, region, shift, shift_low, shift_high, status FROM duckhts_depth_region_fraction('corrected_bins', 'chromosomes');

Estimate regional depth shifts and reference-calibrated uncertainty.

Signature:

duckhts_depth_region_fraction_ref(bins_table, regions_table, reference_table)

Returns:

table(sample_id VARCHAR, region VARCHAR, bins BIGINT, baseline_bins BIGINT, ratio DOUBLE, shift DOUBLE, shift_se DOUBLE, shift_low DOUBLE, shift_high DOUBLE, dispersion DOUBLE, reference_shift DOUBLE, reference_samples BIGINT, status VARCHAR, method VARCHAR)

duckhts_depth_region_fraction uses the base inputs. duckhts_depth_region_fraction_ref requires reference_table as an additional positional argument.

SELECT sample_id, region, shift, shift_low, shift_high, reference_samples FROM duckhts_depth_region_fraction_ref('corrected_bins', 'chromosomes', 'reference_bins');

Estimate the share of a second genome of the other sex from the depth of X and Y relative to the autosomes.

Signature:

duckhts_depth_sex_fraction(bins_table)

Returns:

table(sample_id VARCHAR, autosomal_bins BIGINT, x_bins BIGINT, y_bins BIGINT, dispersion DOUBLE, x_copies DOUBLE, x_copies_se DOUBLE, y_copies DOUBLE, y_copies_se DOUBLE, major_sex VARCHAR, fraction_x DOUBLE, fraction_x_se DOUBLE, fraction_y DOUBLE, fraction_y_se DOUBLE, fraction DOUBLE, fraction_se DOUBLE, fraction_low DOUBLE, fraction_high DOUBLE, agreement_z DOUBLE, reference_samples BIGINT, status VARCHAR, method VARCHAR)

duckhts_depth_sex_fraction uses the ideal copy levels; duckhts_depth_sex_fraction_ref takes reference_table. bins_table holds corrected bins as for duckhts_depth_region_fraction. X, Y and MT are the contigs whose duckhts_contig_key is X, Y or MT; every other contig is autosomal. Remove the pseudo-autosomal regions of X, and any Y bins that are not single-copy in an XY genome, in SQL first: they are not one X or one Y copy. A NULL kept counts as not kept; a NULL chrom is an error.

x_copies and y_copies are twice the depth ratio of the kept X or Y bins to the kept autosomal bins (sums of reads over sums of expected), after the calibration below. The sample is a major genome and a second genome of the other sex with a DNA share f below 0.5. The major genome is XX when x_copies >= 1.5 (by y_copies < 0.5 without X 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. fraction_x and fraction_y are the two estimates, unclipped, so a value below 0 means no second genome is seen. fraction combines them by inverse variance; agreement_z = (fraction_x - fraction_y) / sqrt(fraction_x_se^2 + fraction_y_se^2). An X or Y copy change of the major genome, such as loss of Y in part of the cells or an X aneuploidy, moves the two estimates differently.

Without a reference, one X copy is 1, two are 2, one Y copy is 1 and the Y of a genome without Y is 0. With reference_table, corrected bins of reference samples with a sex column (XX or XY, one value per sample), the XY references set the depth of one X copy and one Y copy, and the XX references, when there are any, set two X copies and the Y background of a genome without Y; without XX references two X copies are twice one and the background is 0. The tested sample is left out of its own reference. Mismapped reads give an XX genome some Y depth, which without XX references biases fraction_y upward, and because its standard error comes from a few reads, also fraction.

Standard errors use the quasi-Poisson bin model of duckhts_depth_region_fraction with the dispersion of the autosomal bins; the calibration levels are treated as known. The interval is fraction +- 1.959964 fraction_se.

One row per sample of bins_table. status is ok, no_autosomal_bins, no_xy_reference (a reference without XY samples), bad_reference (levels that do not separate one copy from two or from none), no_sex_bins, no_x_bins or no_y_bins (the other chromosome gives the estimate), or x_y_disagree (absolute agreement_z above 3). A reference sex other than XX or XY is an error. method names the estimator and its version.

One pass over the input into a few sums per sample. The estimate is evidence of X and Y dosage, not a determination of sex.

duckhts_depth_sex_fraction uses the base inputs. duckhts_depth_sex_fraction_ref requires reference_table as an additional positional argument.

SELECT sample_id, major_sex, fraction, fraction_low, fraction_high, status FROM duckhts_depth_sex_fraction('corrected_bins');

Estimate sex-copy mixture fractions using XX and XY reference bins.

Signature:

duckhts_depth_sex_fraction_ref(bins_table, reference_table)

Returns:

table(sample_id VARCHAR, autosomal_bins BIGINT, x_bins BIGINT, y_bins BIGINT, dispersion DOUBLE, x_copies DOUBLE, x_copies_se DOUBLE, y_copies DOUBLE, y_copies_se DOUBLE, major_sex VARCHAR, fraction_x DOUBLE, fraction_x_se DOUBLE, fraction_y DOUBLE, fraction_y_se DOUBLE, fraction DOUBLE, fraction_se DOUBLE, fraction_low DOUBLE, fraction_high DOUBLE, agreement_z DOUBLE, reference_samples BIGINT, status VARCHAR, method VARCHAR)

One pass over the input into a few sums per sample. The estimate is evidence of X and Y dosage, not a determination of sex.

duckhts_depth_sex_fraction uses the base inputs. duckhts_depth_sex_fraction_ref requires reference_table as an additional positional argument.

SELECT sample_id, major_sex, fraction, fraction_low, fraction_high, status FROM duckhts_depth_sex_fraction_ref('corrected_bins', 'reference_bins_with_sex');

Fit the ichorCNA copy-number state model to corrected bins: tumour fraction, ploidy and rank of every start.

Signature:

duckhts_depth_state_solutions(bins_table, normal_fractions := [0.5, 0.6, 0.7, 0.8, 0.9], ploidies := [2.0, 3.0], max_cn := 7, txn_e := 0.9999999, txn_strength := 10000000.0, max_iter := 20, min_segment_bins := 50, alt_fraction_threshold := 0.05, close_margin := 10.0, max_bins := 1000000)

Returns:

table(sample_id VARCHAR, solution INTEGER, start_normal DOUBLE, start_ploidy DOUBLE, normal_fraction DOUBLE, tumour_fraction DOUBLE, ploidy DOUBLE, log_likelihood DOUBLE, altered_fraction DOUBLE, iterations INTEGER, status VARCHAR, segments BIGINT, bins BIGINT, upstream VARCHAR, method VARCHAR)

bins_table holds corrected bins with sample_id, chrom, start, end, log2_ratio and kept, as duckhts_depth_bins returns them. Only kept bins take part; a kept bin whose log2_ratio is NULL takes every state with likelihood 1, as upstream carries a bin without a ratio. Bins are ordered by chrom then start within a sample. The model learns from the autosomal bins (duckhts_contig_key not X, Y or MT) and decodes every bin. Restrict contigs and remove problem bins in SQL first.

A port of ichorCNA v0.6.1 (GavinHaLab/ichorCNA, Adalsteinsson, Ha, Freeman et al. 2017) for one sample with clonal states of copy number 1 to max_cn and the Student's t emission (nu 2.1) on the natural-log ratio, whose state means are log(((1 - n) c + 2 n) / ((1 - n) phi + 2 n)) for normal fraction n and ploidy phi. Expectation maximisation with the upstream priors (Beta(1, 1) on n, Gamma(4, 0.75) on phi, Gamma(2, beta_k) on each state precision, Dirichlet(200, 800 for two copies) on the initial distribution and Dirichlet(A * txn_strength) on the transitions) runs from each start of ploidies x normal_fractions for at most max_iter iterations until the objective changes by less than 1e-4 relative; a normal start of 0.95 runs only at ploidy 2, as upstream. The forward-backward and Viterbi passes follow HMMcopy 1.54.0. Upstream maximises each parameter with L-BFGS-B; this port climbs each one-dimensional objective from the same start, so values can differ in the last digits.

Solutions are ranked by their objective, the data log-likelihood plus log prior. Upstream's guard sets normal_fraction to 1 (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 lower-ranked solution whose objective is within close_margin of the best is close_to_best (a DuckHTS addition); the rest are ok. tumour_fraction is 1 - normal_fraction.

One row per sample and start with the solution columns, the number of its segments and bins, upstream ('ichorCNA v0.6.1') and method ('duckhts_depth_state v1'). altered_fraction is the share of training bins whose Viterbi state is not two copies. iterations counts the EM iterations whose parameters the solution reports.

One decode per sample and macro call, on one thread per sample. The model holds about max_cn * 28 + 40 bytes per bin and one 48-byte record per segment and start; max_bins bounds the bins of a sample and the start grid has at most 100 points. A sample without copy-number change gives no fraction: the guard sets it to zero, or the fit converges to a normal fraction near 1, where the state means lie within the noise of each other and whole chromosomes can take any state, as in upstream; the copy numbers of such a solution are not evidence, its fraction is. Fraction and ploidy trade against each other, so every start is returned, and a genome with one altered region does not identify the fraction: 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 between them, as in upstream; changes of more than one magnitude identify the fraction. Depth does not see copy-neutral changes, and a constitutional copy-number variant looks like a change in every cell. Parallelism across samples depends on how DuckDB spreads the per-sample rows over threads. Kept upstream quirks are listed in design/depth_fraction_lane.md.

test/sql/depth_state_model.test and the R tests recover a planted normal fraction of 0.7 within 0.05, ploidy 2 within 0.15, whole-chromosome and arm-level changes over at least 95% and 90% of their bins, and the guard on an unchanged genome. benchmarks/benchmark_depth_state_model.md runs the unmodified upstream code on the same corrected bins of 16 real low-pass composites: where both land on the same optimum they agree to four decimals with every bin in the same state; on flat objectives (no change, one small event) the two optimizers reach different local maxima; latency is linear in bins.

SELECT sample_id, tumour_fraction, ploidy, altered_fraction, status FROM duckhts_depth_state_solutions('corrected_bins') WHERE solution = 1;

Viterbi copy-number segments of every solution of the ichorCNA state model over corrected bins.

Signature:

duckhts_depth_state_segments(bins_table, normal_fractions := [0.5, 0.6, 0.7, 0.8, 0.9], ploidies := [2.0, 3.0], max_cn := 7, txn_e := 0.9999999, txn_strength := 10000000.0, max_iter := 20, min_segment_bins := 50, alt_fraction_threshold := 0.05, close_margin := 10.0, max_bins := 1000000)

Returns:

table(sample_id VARCHAR, solution INTEGER, start_normal DOUBLE, start_ploidy DOUBLE, normal_fraction DOUBLE, tumour_fraction DOUBLE, ploidy DOUBLE, log_likelihood DOUBLE, altered_fraction DOUBLE, iterations INTEGER, status VARCHAR, chrom VARCHAR, start BIGINT, end BIGINT, bins INTEGER, state INTEGER, copy_number INTEGER, median_log2_ratio DOUBLE, upstream VARCHAR, method VARCHAR)

The same relation and options as duckhts_depth_state_solutions; the two macros share one decode per call.

One row per sample, solution and segment: the solution columns of duckhts_depth_state_solutions, then chrom, start and end (the first bin's start and the last bin's end), bins, state (the copy number of the Viterbi state), copy_number (the state after upstream's integer correction, which re-reads segments at the highest state or with a median implying at least 1.2 max_cn copies from their median and the fitted purity and ploidy when the purity is at least 0.05), median_log2_ratio (the median of the segment's non-NULL ratios), upstream and method. A segment never crosses a contig.

As duckhts_depth_state_solutions. Segment order is unspecified; order by sample_id, solution, chrom and start when it matters.

SELECT chrom, start, "end", bins, copy_number, median_log2_ratio FROM duckhts_depth_state_segments('corrected_bins') WHERE sample_id = 'S1' AND solution = 1 AND copy_number <> 2 ORDER BY chrom, start;

Count aligned read bases against the reference by mate, cycle, base quality and substitution, optionally leaving out masked positions such as known variants.

Signature:

duckhts_bam_mismatch_counts(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)

Returns:

table

path is a SAM, BAM or CRAM file and reference is its indexed FASTA; fasta_index() builds the index. region is a comma-separated list of contigs or intervals and needs an alignment index; without it the whole file is read. Alignments are kept when MAPQ is at least min_mapq, every bit of require_flags is set and no bit of exclude_flags is set; the default 3844 leaves out unmapped, secondary, failed, duplicate and supplementary alignments.

One row for each combination that occurs: mate (1 or 2 for a paired read, 0 otherwise), cycle (the one-based position of the base in the read as sequenced, hard clips included), base_quality (NULL when the read stores no qualities), reference_base and read_base (both on the strand of the read, so a reverse alignment is complemented), and bases, the number of aligned bases. Rows with reference_base equal to read_base are the matches, so a rate is a ratio of two sums. Rows come in the order of mate, cycle, base_quality, reference_base, read_base.

Every aligned base (CIGAR M, = or X) of a kept alignment whose read base and reference base are both A, C, G or T. A read base written as = (the SAM sign for the reference base) is the reference base. A base within indel_flank read bases of an insertion, a deletion, a reference skip or a soft clip is left out, because bases next to a gap are often misaligned. A region selects alignments by overlap; all aligned bases of a selected alignment are counted, also those outside the region.

mask is an indexed BCF or a bgzip-compressed VCF with a tabix index, for example known variants. Every reference position of a mask record is left out. A record with an indel or a symbolic allele also leaves out indel_flank bases on each side. The span of a record is the one htslib gives it, the same for a BCF and a VCF mask: the longest of REF, INFO/END, INFO/SVLEN of a deletion, duplication or copy-number allele, and FORMAT/LEN of a reference block. A VCF mask line is parsed up to INFO unless its FORMAT names LEN, so the genotypes of a population VCF are not parsed; each line is still decompressed and read whole into one buffer. The mask removes the variants it holds and no others: a variant of the sample that the mask does not hold, such as a private, de novo or somatic variant or an allele of a second genome, is still counted as a mismatch. So with a mask of known variants the mismatches that remain are mostly errors of the read, the library or the alignment, and their rate is an upper bound on the error rate; it is the error rate only when the mask holds every variant site of the genomes in the sample.

The contig names of the alignment header, the reference and the mask are compared byte for byte. A contig with a kept alignment that the reference does not have is an error. So is a contig that the mask does not know, in its header or in its index: a silent miss would count known variants as errors.

One worker. Memory does not grow with the alignments read: 36.5 MB of counters, one reference window of 1 MiB with its mask, the alignment in hand and the longest mask line read, plus htslib's state for the file. For CRAM that includes the reference of the slice being decoded; a slice over a sparse region can span many megabases. A cycle above 1000 is counted as 1000 and a base quality above 93 as 93. The window moves along an alignment that reaches past it, so no alignment needs more reference in memory; an alignment with a reference skip of more than 1 MiB fetches a window on each side of the skip, and many such alignments are slow. Alignments are expected in coordinate order; an unsorted file is counted correctly but fetches a reference window for each alignment.

SELECT base_quality, sum(bases) FILTER (WHERE reference_base != read_base) / sum(bases) AS mismatch_rate FROM duckhts_bam_mismatch_counts('sample.bam', 'reference.fa', region := 'chr20:10000000-12000000', mask := 'known_variants.vcf.gz') GROUP BY base_quality ORDER BY base_quality;

Count where read 1 and read 2 of one fragment agree and disagree at the positions both align to, by base quality and base, as an estimate of sequencing error.

Signature:

duckhts_bam_mate_overlap_counts(path, region := NULL, index_path := NULL, reference := NULL, reference_index_path := NULL, min_mapq := 20, require_flags := 0, exclude_flags := 3844, max_pending := 1000000)

Returns:

table

path is a paired-end SAM, BAM or CRAM file in coordinate order; alignments out of order are an error. region is a comma-separated list of contigs or intervals and needs an alignment index; without it the whole file is read. reference is an indexed FASTA; fasta_index() builds the index. It is required for CRAM. Its only other use is to resolve a read base stored as '='; without a reference such a base is left out. Alignments are kept when MAPQ is at least min_mapq, every bit of require_flags is set and no bit of exclude_flags is set; the default 3844 leaves out unmapped, secondary, failed, duplicate and supplementary alignments.

One row for each combination that occurs: read1_quality and read2_quality (the base qualities, NULL when the read stores no qualities), read1_base and read2_base (both on the forward strand of the reference, so a real allele gives equal bases), and positions, the number of shared positions. Read 1 is the mate with flag 0x40 and read 2 the mate with flag 0x80, whichever comes first. Rows with read1_base different from read2_base are the disagreements, so a rate is a ratio of two sums. Rows come in the order of read1_quality, read2_quality, read1_base, read2_base, with a missing quality last.

Each reference position where both mates of a pair have an aligned base (CIGAR M, = or X), once per pair. Both mates are kept primary alignments on the same contig, the pair flag (0x1) is set and the mate-unmapped flag (0x8) is not, each names the other's start, and one is read 1 and the other read 2. Soft-clipped and inserted bases are not counted, nor are deleted or skipped (CIGAR N) positions, and a position where either base is not A, C, G or T is left out. 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, in the library, by PCR or by DNA damage, show in both mates and are not seen. Only fragments shorter than the two reads overlap, so the count leans toward short fragments and read ends. A region selects alignments by overlap, and a pair counts only when both of its mates are selected.

One worker. The counters are fixed at 1.2 MB, and a base quality above 93 is counted as 93. The first mate of a pair is held until its mate arrives, with its name, CIGAR, sequence and qualities: about 400 bytes plus 1.5 bytes per read base and 4 bytes per CIGAR operation. A held mate whose mate start the scan has passed without a match is dropped, and a change of contig drops every held mate. Mates are matched by read name. A record that has the name of a held mate but does not pair with it is left out, and the held mate stays held for its own mate. SAM gives one name to one template, so such a record occurs only in a malformed file. More than max_pending (1 to 100,000,000) mates held at once is an error that names the limit. Held entries are reused but kept until the scan ends, so held memory is at most max_pending times the size of the longest alignment held, and in practice follows the most mates held at once.

SELECT read1_quality, sum(positions) FILTER (WHERE read1_base != read2_base) / sum(positions) AS discordance FROM duckhts_bam_mate_overlap_counts('sample.bam', region := 'chr20:10000000-12000000') GROUP BY read1_quality ORDER BY read1_quality;

Compute samtools coverage-like regional summaries for BAM or CRAM input over a BED target set, returning one row per BED interval with DuckHTS-specific pre/post-filter read counts, covered bases, percentage covered, mean depth, mean baseQ, mean mapQ, and strand-specific post-filter summaries in read mode. Indexed BAM/CRAM input is required in the current implementation. decompression_threads controls htslib worker threads for BAM/CRAM decoding; use 0 to disable them.

Signature:

duckhts_bam_bed_coverage(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 := 1000000, decompression_threads := 0, fragment_mode := FALSE, strand_outputs := TRUE, processing_threads := 0)

Returns:

table

SELECT chrom, start, "end", numreads_post, covbases_post, coverage_post FROM duckhts_bam_bed_coverage('fixture_mixed.bam', 'fixture_mixed_regions.bed');

Write mosdepth-compatible coverage files from indexed BAM/CRAM.

Signature:

duckhts_mosdepth(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)

Returns:

table

Produce summary, global distribution, per-base BED.gz/CSI, optional window/BED-region results, quantized BED.gz/CSI and threshold counts for by. precision_digits sets text decimal places.

Default fast_mode=FALSE uses CIGAR-aware coverage and mate-overlap correction. fragment_mode counts full insert spans of proper pairs; use_median switches by output from mean to median. read_groups filters comma-separated RG IDs; min_frag_len/max_frag_len filter absolute template length. Supply fasta when CRAM requires reference.

processing_threads=0 is sequential; positive values select parallel contig-worker count.

SELECT * FROM duckhts_mosdepth('sample', 'range.cram', fasta := 'ce.fa', by := '1000', fragment_mode := TRUE, read_groups := '1', use_median := TRUE, min_frag_len := 50, max_frag_len := 500, quantize := ':1:4:', thresholds := '1,10,20', precision_digits := 4, overwrite := TRUE);

Write samtools idxstats-compatible TAB-delimited output for BAM, CRAM, or SAM input. Indexed BAM uses hts_idx_get_stat(...) for the fast path; CRAM, SAM, and unindexed BAM fall back to a full scan while preserving samtools-style contig rows plus the final * row.

Signature:

duckhts_samtools_idxstats(path, output := NULL, index_path := NULL, threads := 0, overwrite := FALSE)

Returns:

table

SELECT * FROM duckhts_samtools_idxstats('range.bam', output := 'range.idxstats.txt', overwrite := TRUE);

Inspect HTS headers in parsed, raw, or combined form across supported formats. Raw VCF/BCF mode includes the final #CHROM sample header line so the returned text is suitable for Parquet metadata and future VCF/BCF regeneration.

Signature:

read_hts_header(path, format := NULL, mode := NULL)

Returns:

table

key_values leaves out a pair whose key is not valid UTF-8, since a MAP key cannot be NULL; a value that is not valid UTF-8 is NULL.

SELECT record_type, id FROM read_hts_header('formatcols.vcf.gz') LIMIT 10;

Inspect high-level HTS index metadata such as sequence names and mapped counts.

Signature:

read_hts_index(path, format := NULL, index_path := NULL)

Returns:

table

SELECT seqname, index_type FROM read_hts_index('vcf_file.bcf');

Expand index metadata into span and chunk rows suitable for low-level index inspection.

Signature:

read_hts_index_spans(path, format := NULL, index_path := NULL)

Returns:

table

SELECT seqname, chunk_beg_vo, chunk_end_vo FROM read_hts_index_spans('vcf_file.bcf') LIMIT 5;

Return the raw on-disk HTS index blob together with basic identifying metadata.

Signature:

read_hts_index_raw(path, format := NULL, index_path := NULL)

Returns:

table

SELECT length(raw) FROM read_hts_index_raw('formatcols.vcf.gz');

Encode a normalized biallelic variant as an official VariantKey-compatible 64-bit unsigned integer. This DuckHTS wrapper accepts 1-based VCF/DuckHTS POS to match bcftools %VKX / +add-variantkey, internally converts to the upstream 0-based field, and preserves the official hashed nonreversible mode for large, ambiguous, and symbolic REF/ALT strings. Only CHROM, POS, REF, and ALT are encoded; END, SVLEN, mate breakend coordinates, and other SV metadata are not.

Signature:

variantkey(chrom, pos, ref, alt)

Returns:

UBIGINT

The chromosome part of a VariantKey is a code of the upstream format, not a contig key. A leading chr is removed in any letter case. An all-digit name becomes its number modulo 256, with no range check, so 257 shares the code of 1. X, Y, and M or MT, in any letter case, become 23, 24 and 25, so 25 and MT share a code. Any other name has no code, and neither has a number that comes to 0: variantkey() and regionkey() return NULL for it, so scaffolds and accessions have no key rather than a shared one. Use duckhts_contig_key() to join on contig names.

SELECT variantkey_hex(variantkey('1', 324684, 'C', 'G'));

Render a VariantKey as its lowercase 16-character hexadecimal string representation.

Signature:

variantkey_hex(vk)

Returns:

VARCHAR

SELECT variantkey_hex(variantkey('1', 324684, 'C', 'G'));

Parse a 16-character hexadecimal VariantKey string back into its UBIGINT code. Invalid or non-hex strings return NULL.

Signature:

parse_variantkey_hex(hex)

Returns:

UBIGINT

SELECT parse_variantkey_hex('08027a2588b00000');

Encode the raw upstream VariantKey fields directly: chromosome code, 0-based position, and 31-bit REF+ALT code.

Signature:

encode_variantkey(chrom_code, pos0, refalt_code)

Returns:

UBIGINT

SELECT variantkey_hex(encode_variantkey(1, 324683, 145752064));

Extract the raw upstream VariantKey chromosome code.

Signature:

extract_variantkey_chrom(vk)

Returns:

UTINYINT

SELECT extract_variantkey_chrom(parse_variantkey_hex('08027a2588b00000'));

Extract the raw upstream VariantKey 0-based position field.

Signature:

extract_variantkey_pos(vk)

Returns:

UINTEGER

SELECT extract_variantkey_pos(parse_variantkey_hex('08027a2588b00000'));

Extract the raw upstream 31-bit VariantKey REF+ALT code.

Signature:

extract_variantkey_refalt(vk)

Returns:

UINTEGER

SELECT extract_variantkey_refalt(parse_variantkey_hex('08027a2588b00000'));

Decode a VariantKey into its raw upstream numeric fields: chrom_code, pos0, and refalt_code.

Signature:

decode_variantkey(vk)

Returns:

STRUCT

SELECT (decode_variantkey(parse_variantkey_hex('08027a2588b00000'))).pos0;

Decode a VariantKey into a STRUCT with chrom, chrom_code, 1-based pos, upstream 0-based pos0, ref, alt, refalt_code, and reversible. For hashed nonreversible keys, reversible is FALSE and ref/alt are returned as NULL because DuckHTS v1 does not ship the optional NRVK lookup sidecar.

Signature:

reverse_variantkey(vk)

Returns:

STRUCT

SELECT (reverse_variantkey(parse_variantkey_hex('08027a2588b00000'))).ref;

Return the inclusive minimum and maximum VariantKey bounds for a chromosome plus 1-based VCF position range, suitable for numeric range filtering on precomputed VariantKeys.

Signature:

variantkey_range(chrom, pos_min, pos_max)

Returns:

STRUCT

SELECT variantkey_hex((variantkey_range('1', 100, 100)).min), variantkey_hex((variantkey_range('1', 100, 100)).max);

Encode a genomic interval as an official RegionKey-compatible 64-bit unsigned integer. Start and end use 0-based half-open interval semantics, matching BED-style coordinates; strand accepts -1, 0, or 1.

Signature:

regionkey(chrom, start, end, strand := 0)

Returns:

UBIGINT

SELECT regionkey_hex(regionkey('X', 1007, 1807, 1));

Render a RegionKey as its lowercase 16-character hexadecimal string representation.

Signature:

regionkey_hex(rk)

Returns:

VARCHAR

SELECT regionkey_hex(regionkey('X', 1007, 1807, 1));

Parse a 16-character hexadecimal RegionKey string back into its UBIGINT code. Invalid or non-hex strings return NULL.

Signature:

parse_regionkey_hex(hex)

Returns:

UBIGINT

SELECT parse_regionkey_hex('b80001f78000387a');

Encode the raw upstream RegionKey fields directly: chromosome code, 0-based start, 0-based end, and strand code (0 = unknown, 1 = +, 2 = -).

Signature:

encode_regionkey(chrom_code, start, end, strand_code)

Returns:

UBIGINT

SELECT regionkey_hex(encode_regionkey(23, 1007, 1807, 1));

Extract the raw upstream RegionKey chromosome code.

Signature:

extract_regionkey_chrom(rk)

Returns:

UTINYINT

SELECT extract_regionkey_chrom(parse_regionkey_hex('b80001f78000387a'));

Extract the raw upstream RegionKey 0-based start position.

Signature:

extract_regionkey_startpos(rk)

Returns:

UINTEGER

SELECT extract_regionkey_startpos(parse_regionkey_hex('b80001f78000387a'));

Extract the raw upstream RegionKey 0-based end position.

Signature:

extract_regionkey_endpos(rk)

Returns:

UINTEGER

SELECT extract_regionkey_endpos(parse_regionkey_hex('b80001f78000387a'));

Extract the raw upstream RegionKey strand code (0 = unknown, 1 = +, 2 = -).

Signature:

extract_regionkey_strand(rk)

Returns:

UTINYINT

SELECT extract_regionkey_strand(parse_regionkey_hex('b80001f78000387a'));

Decode a RegionKey into its raw upstream numeric fields: chrom_code, start, end, and strand_code.

Signature:

decode_regionkey(rk)

Returns:

STRUCT

SELECT (decode_regionkey(parse_regionkey_hex('b80001f78000387a'))).strand_code;

Decode a RegionKey into a STRUCT with chrom, chrom_code, start, end, strand, and strand_code.

Signature:

reverse_regionkey(rk)

Returns:

STRUCT

SELECT (reverse_regionkey(parse_regionkey_hex('b80001f78000387a'))).strand;

Extend a RegionKey interval by a fixed number of bases on both sides, clamping to the official 28-bit RegionKey position range.

Signature:

extend_regionkey(rk, size)

Returns:

UBIGINT

SELECT reverse_regionkey(extend_regionkey(regionkey('X', 10000, 20000, -1), 1000));

Return a conservative contig join key by removing one non-empty leading chr prefix case-insensitively and normalizing M/MT to MT. X and Y are uppercased; all other suffixes are preserved. This does not map numeric sex chromosomes, accessions, patches, or alternate loci.

Signature:

duckhts_contig_key(contig)

Returns:

VARCHAR

Readers and region arguments compare contig names byte for byte. This function is the explicit key for a join of two sources that spell contigs differently: join on duckhts_contig_key(a.chrom) = duckhts_contig_key(b.chrom), or rename one side. Two names of one source can share a key (1 and chr1), so a join on the key repeats rows when one side has both. duckhts_roh_ancestry and the R ancestry wrappers use this key. Functions with another rule state it under Contig names.

Leading zeros (01 is not 1), numeric sex and mitochondrial codes (23, 24, 25, 26), accessions such as NC_000001.11, patches, alternate loci and unplaced scaffolds. The letter case of other names is kept.

SELECT duckhts_contig_key('chr1'), duckhts_contig_key('chrM');

Return TRUE when two explicit 0-based half-open intervals overlap on the same canonical chromosome.

Signature:

are_overlapping_regions(chrom_a, start_a, end_a, chrom_b, start_b, end_b)

Returns:

BOOLEAN

SELECT are_overlapping_regions('1', 2, 4, '1', 3, 7);

Return TRUE when a 0-based half-open interval overlaps the supplied RegionKey interval.

Signature:

are_overlapping_region_regionkey(chrom, start, end, rk)

Returns:

BOOLEAN

SELECT are_overlapping_region_regionkey('X', 1008, 1800, parse_regionkey_hex('b80001f78000387a'));

Return TRUE when two RegionKeys overlap.

Signature:

are_overlapping_regionkeys(rka, rkb)

Returns:

BOOLEAN

SELECT are_overlapping_regionkeys(regionkey('X', 1007, 1807, 1), parse_regionkey_hex('b80001f78000387a'));

Row-oriented liftover kernel intended to mirror bcftools +liftover semantics as closely as possible while returning one STRUCT per input row with fields: src_chrom, src_pos, src_ref, src_alt, dest_chrom, dest_pos, dest_end, dest_ref, dest_alt, mapped, reverse_complemented, swap, reject_reason, and note. Set no_left_align := true to skip post-liftover left-alignment of lifted indels (mirrors --no-left-align in bcftools +liftover).

Signature:

bcftools_liftover(chrom, pos, ref, alt, chain_path, dst_fasta_ref, src_fasta_ref, max_snp_gap, max_indel_inc, lift_mt, end_pos, no_left_align)

Returns:

STRUCT

Input names and chain source names are compared after one rule: a leading chr is removed in any letter case; M, MT and 26 become MT; 23, 25, X, XY, XX, PAR1 and PAR2 become X; 24 and Y become Y. A lifted record carries the first name that the destination FASTA index has among: the chain's destination name; that name with chr added or removed; its form under the rule above, without and with chr; and for the mitochondrion MT, chrM and M. A FASTA that holds two of these names is not an error. A mitochondrial record that passes through without lifting takes the first of chrM, MT and M that the destination FASTA has.

SELECT (bcftools_liftover(chrom, pos, ref, alt, 'hg19ToHg38.over.chain.gz', 'hg38.fa', 'hg19.fa', 1, 250, false, NULL::BIGINT, false)).dest_pos FROM variants;

Lift variants from a table or derived-table expression with bcftools_liftover.

Signature:

duckhts_liftover(table_name, chrom_col, pos_col, ref_col := NULL, alt_col := NULL, chain_path := NULL, dst_fasta_ref := NULL, src_fasta_ref := NULL, max_snp_gap := 1, max_indel_inc := 250, lift_mt := false, end_pos_col := NULL, no_left_align := false)

Returns:

table

Column-name strings select chrom/pos/ref/alt. The no_left_align parameter mirrors --no-left-align in bcftools +liftover.

SELECT * FROM duckhts_liftover('variants', 'chrom', 'pos', ref_col := 'ref', alt_col := 'alt', chain_path := 'hg19ToHg38.over.chain.gz', dst_fasta_ref := 'hg38.fa');

SELECT * FROM duckhts_liftover('(SELECT chrom, pos, ref, alt FROM variants) AS v', 'chrom', 'pos', ref_col := 'ref', alt_col := 'alt', chain_path := 'hg19ToHg38.over.chain.gz', dst_fasta_ref := 'hg38.fa');

Normalize one variant against FASTA with bcftools/vt-style left alignment.

Signature:

bcftools_norm_row(chrom, pos, ref, alt, fasta_ref, end_pos := NULL, svlen := NULL, fasta_index_path := NULL, gzi_path := NULL)

Returns:

STRUCT

alt accepts comma-delimited VARCHAR or VARCHAR[]. Symbolic <DEL> may use end_pos; <DUP> may use svlen.

Return pos_normed, end_pos_normed, ref_normed, alt_normed (always VARCHAR[]), nullable normed and norm_status. gVCF <NON_REF>/<*> reference blocks pass through with GVCFReferenceBlock. Mixed real/gVCF-symbolic rows normalize real alleles while preserving symbolic alleles and supplied reference-block END.

The reference sequence of a record is looked up in the FASTA index under, in order: the record's name; the name without chr (any letter case), or with chr added when it has none; and for M, MT or chrM also MT, chrM and M. The first name that the index has is used. A FASTA that holds both 1 and chr1 is not an error.

SELECT nr.pos_normed, nr.ref_normed, nr.alt_normed FROM (SELECT bcftools_norm_row('chr1', 100, 'AC', 'A', 'ref.fa', NULL::BIGINT, NULL::BIGINT, NULL, NULL) AS nr);

Normalize variants from a table or derived-table expression while preserving input columns.

Signature:

duckhts_bcftools_norm(table_name, 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)

Returns:

table

ALT accepts VARCHAR or VARCHAR[]. Append pos_normed, end_pos_normed, ref_normed, alt_normed, normed and norm_status. split_multiallelic=TRUE splits sites before normalization; alt_normed becomes VARCHAR and alt_index is added.

This wraps bcftools_norm_row, not a full-record VCF/BCF rewrite. GT, PL/GP/DS and PS remain unchanged caller columns unless a separate writer/remapper updates them.

SELECT * FROM duckhts_bcftools_norm('variants', 'ref.fa');

SELECT * FROM duckhts_bcftools_norm('(SELECT CHROM AS chrom, POS AS pos, REF AS ref, ALT AS alt FROM read_bcf(''cohort.vcf.gz'')) AS v', 'ref.fa', split_multiallelic := TRUE);

Compute polygenic scores from genotype VCF/BCF and summary statistics using bcftools +score dosage semantics.

Signature:

bcftools_score(bcf_path, summary_path_or_list, 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)

Returns:

table

Support GT/DS/HDS/AP/GP/AS dosage, sample subsets and region/target/FILTER-string controls. The second argument accepts one path or a list. TSV/SSF inputs yield one PRS column per file in a single genotype scan; GWAS-VCF yields one per FORMAT sample.

With NULL second argument, summaries_list_file reads paths from a file or directory. List entries are interpreted as written; directories scan supported regular files lexicographically and omit index sidecars.

log_path writes per-PRS loaded/matched/allele-mismatch/duplicate-marker counts.

Chromosome names of the summary statistics are looked up in the genotype VCF header as the upstream plugin does, in this order: the exact name; the name without a lowercase chr prefix; chr plus a name of at most two characters; 23, 25, XY, XX, PAR1 and PAR2 as X or chrX; 24 as Y or chrY; 26, MT and chrM as MT or chrM. The lookup is case-sensitive, and M alone does not find MT. A marker whose chromosome is not found is skipped without an error.

SELECT * FROM bcftools_score('cohort.bcf', 'gwas.tsv.gz', columns := 'PLINK') LIMIT 5;

SELECT * FROM bcftools_score('cohort.bcf', ['score1.tsv.gz', 'score2.tsv.gz'], columns := 'GWAS-SSF') LIMIT 5;

SELECT * FROM bcftools_score('cohort.bcf', NULL, columns := 'GWAS-SSF', summaries_list_file := 'scores.list', log_path := 'score.log') LIMIT 5;

Normalize one summary-statistics row into GWAS-VCF-style fields (chrom/pos/ref/alt/effect metrics), resolving REF/ALT orientation against a FASTA reference and applying swap-aware sign/frequency/count transforms. The output flag alleles_swapped means REF/ALT orientation was swapped to match the FASTA reference.

Signature:

bcftools_munge_row(chrom, pos, a1, a2, id, p, z, or, beta, n, n_cas, n_con, info, frq, se, lp, ac, neff, neffdiv2, het_i2, het_p, het_lp, dire, fasta_ref, iffy_tag := 'IFFY', mismatch_tag := 'REF_MISMATCH', ns := NULL, nc := NULL, ne := NULL)

Returns:

STRUCT

SELECT (bcftools_munge_row('chr1', 12345, 'A', 'G', 'rs1', 0.01, NULL, NULL, 0.12, 1000, NULL, NULL, NULL, 0.3, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, 'ref.fa')).ref;

Map source columns and normalize summary statistics with bcftools_munge_row.

Signature:

duckhts_munge(table_name, preset := '', column_map := map([''], ['']), column_map_file := '', fasta_ref := NULL, iffy_tag := 'IFFY', mismatch_tag := 'REF_MISMATCH', ns := NULL, nc := NULL, ne := NULL)

Returns:

table

Select columns via preset or explicit map to return normalized GWAS-VCF-style rows with explicit alleles_swapped semantics.

Output columns: chrom, pos, id, ref, alt, alleles_swapped, filter, ns, ez, nc, es, se, lp, af, ac, ne (16 columns). For METAL meta-analysis output with SI/I2/CQ/ED columns, use duckhts_munge_metal.

SELECT * FROM duckhts_munge('gwas_table', preset := 'PLINK', fasta_ref := 'ref.fa');

SELECT * FROM duckhts_munge('(SELECT * FROM gwas_table) AS s', column_map := map(['CHR','BP','A1','A2','SNP'], ['chrom','pos','ea','nea','id']), fasta_ref := 'ref.fa');

Normalize summary statistics with additional METAL meta-analysis output columns.

Signature:

duckhts_munge_metal(table_name, preset := '', column_map := map([''], ['']), column_map_file := '', fasta_ref := NULL, iffy_tag := 'IFFY', mismatch_tag := 'REF_MISMATCH', ns := NULL, nc := NULL, ne := NULL)

Returns:

table

Same as duckhts_munge but additionally emits: si (imputation info, from INFO input), i2 (Cochran's I² heterogeneity, from HET_I2), cq (Cochran's Q -log10 p, from HET_LP or -log10(HET_P)), and ed (effect direction string, from DIRE; +/- flipped on allele swap).

rduckhts_munge() auto-dispatches to this macro when metal keys (INFO, HET_I2, HET_P, HET_LP, DIRE) are present in the resolved column map.

SELECT * FROM duckhts_munge_metal('metal_results', preset := 'METAL', fasta_ref := 'ref.fa');

SELECT * FROM duckhts_munge_metal('meta', column_map := map(['CHR','BP','A1','A2','SNP','HET_I2','DIRE'], ['chr','pos','ea','nea','snp','HetISq','Direction']), fasta_ref := 'ref.fa');

Generate a UNION ALL BY NAME query string that reads every file matching a glob pattern through the named reader function. The result includes a 'filename' column identifying the source file for each row. Assign to a variable with SET VARIABLE and execute via query(getvariable(...)). Optional params string is appended to each reader call. In R, use the typed rduckhts_*_multi() helpers instead, which accept file vectors with optional per-file parameters and create DuckDB tables directly.

Signature:

hts_union_query(reader, pattern, params := '')

Returns:

VARCHAR

SET VARIABLE q = hts_union_query('read_bam', 'samples/*.bam'); SELECT * FROM query(getvariable('q'));

SET VARIABLE q = hts_union_query('read_bcf', 'cohort/*.vcf.gz', 'tidy_format := true'); SELECT * FROM query(getvariable('q'));

Generate UNION ALL BY NAME SQL over separate per-region scans of one HTS file.

Signature:

hts_region_union_query(reader, path, regions, params := '')

Returns:

VARCHAR

regions is a list of region strings. params is appended to each reader call and must not include region. Output adds filename, duckhts_region_shard_id and duckhts_region_shard for shard provenance.

UNION ALL does not deduplicate. Adjacent or overlapping shards can repeat spanning BAM/VCF/BCF records; apply shard-local filters or explicit downstream deduplication when exactly-once output is required.

SET VARIABLE q = hts_region_union_query('read_bam', 'sample.bam', ['chr1:1-1000000','chr1:1000001-2000000']); SELECT * FROM query(getvariable('q'));

SET VARIABLE q = hts_region_union_query('read_bcf', 'cohort.vcf.gz', string_split('22:16000000-16999999,22:17000000-17999999', ','), 'tidy_format := true'); SELECT * FROM query(getvariable('q'));

Compute the reverse complement of a DNA sequence using A, C, G, T, and N bases. Overloaded: accepts either a VARCHAR text sequence (returns VARCHAR) or a UTINYINT[] of htslib nt16 codes as produced by read_bam(sequence_encoding := 'nt16') (returns UTINYINT[]); the nt16 overload is bit-identical to the text path after decoding, so BAM pipelines can reverse-complement without leaving the nt16 encoding.

Signature:

seq_revcomp(sequence)

Returns:

VARCHAR

SELECT seq_revcomp('ACGTN');

SELECT seq_revcomp(SEQ) FROM read_bam('reads.bam', sequence_encoding := 'nt16');

Return the lexicographically smaller of a sequence and its reverse complement. Overloaded: accepts either a VARCHAR text sequence (returns VARCHAR) or a UTINYINT[] of htslib nt16 codes as produced by read_bam(sequence_encoding := 'nt16') (returns UTINYINT[]); the nt16 overload compares by decoded base order and is bit-identical to the text path after decoding.

Signature:

seq_canonical(sequence)

Returns:

VARCHAR

SELECT seq_canonical('ACGTN');

SELECT seq_canonical(SEQ) FROM read_bam('reads.bam', sequence_encoding := 'nt16');

Encode a short DNA sequence as a 2-bit unsigned integer hash. Overloaded to also accept a UTINYINT[] of htslib nt16 codes (from read_bam(sequence_encoding := 'nt16')); non-ACGT codes yield NULL, bit-identical to the text path.

Signature:

seq_hash_2bit(sequence)

Returns:

UBIGINT

SELECT seq_hash_2bit('ACGT');

Encode an IUPAC DNA sequence as a list of 4-bit base codes, preserving ambiguity symbols including N.

Signature:

seq_encode_4bit(sequence)

Returns:

UTINYINT[]

SELECT seq_encode_4bit('ACGTRYSWKMBDHVN');

Decode a list of 4-bit IUPAC DNA base codes back into a sequence string.

Signature:

seq_decode_4bit(codes)

Returns:

VARCHAR

SELECT seq_decode_4bit(seq_encode_4bit('ACGTRYSWKMBDHVN'));

Compute GC fraction for a DNA sequence as a value between 0 and 1. Overloaded: accepts either a VARCHAR text sequence or a UTINYINT[] of htslib nt16 codes as produced by read_bam(sequence_encoding := 'nt16'); the nt16 overload classifies codes directly and is bit-identical to the text path, so BAM pipelines can compute GC without decoding sequences back to text.

Signature:

seq_gc_content(sequence)

Returns:

DOUBLE

SELECT seq_gc_content('ACGT');

SELECT seq_gc_content(SEQ) FROM read_bam('reads.bam', sequence_encoding := 'nt16');

Expand a sequence into positional k-mers with optional canonicalization.

Signature:

seq_kmers(sequence, k, canonical := FALSE)

Returns:

table

SELECT * FROM seq_kmers('ACGT', 2);

Decode a SAM flag into a struct of boolean bit fields using explicit SAM-oriented names such as is_paired, is_proper_pair, is_next_segment_unmapped, and is_supplementary.

Signature:

sam_flag_bits(flag)

Returns:

STRUCT

SELECT (sam_flag_bits(99)).is_proper_pair;

Test whether any bits from the provided SAM flag mask are set in a flag value.

Signature:

sam_flag_has(flag, mask)

Returns:

BOOLEAN

SELECT sam_flag_has(99, 2);

Test whether a mapped segment is aligned to the forward strand. Returns NULL for unmapped segments because SAM flag 0x10 does not define genomic strand when 0x4 is set.

Signature:

is_forward_aligned(flag)

Returns:

BOOLEAN

SELECT is_forward_aligned(0);

Test whether a CIGAR string contains any soft-clipped segment (S). Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_has_soft_clip(cigar[, strict])

Returns:

BOOLEAN

See cigar_query_length for full-input validation and missing-input behavior.

SELECT cigar_has_soft_clip('5S90M5S');

Test whether a CIGAR string contains any hard-clipped segment (H). Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_has_hard_clip(cigar[, strict])

Returns:

BOOLEAN

See cigar_query_length for full-input validation and missing-input behavior.

SELECT cigar_has_hard_clip('5H95M');

Return the left-end soft-clipped length from a CIGAR string, or zero if the alignment does not start with S. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_left_soft_clip(cigar[, strict])

Returns:

BIGINT

See cigar_query_length for full-input validation and missing-input behavior. The literal first op determines the left soft clip; a leading H is not skipped.

SELECT cigar_left_soft_clip('5S90M5S');

Return the right-end soft-clipped length from a CIGAR string, or zero if the alignment does not end with S. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_right_soft_clip(cigar[, strict])

Returns:

BIGINT

See cigar_query_length for full-input validation and missing-input behavior. The literal last op determines the right soft clip; a trailing H is not skipped.

SELECT cigar_right_soft_clip('5S90M5S');

Return the query-consuming length from a CIGAR string, counting M, I, S, =, and X. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_query_length(cigar[, strict])

Returns:

BIGINT

Text and packed CIGARs are validated in full. Supported ops are M, I, D, N, S, H, P, = and X, each with a positive length that fits BIGINT. Consumed query and reference spans must each fit BIGINT. Invalid ops, missing lengths, trailing digits, arithmetic overflow and NULL packed elements are invalid input. This validates operation syntax and numeric ranges, not biological ordering constraints.

The optional final positional BOOLEAN strict defaults to FALSE: invalid input returns NULL. TRUE uses the same grammar but raises a DuckDB error naming the function and failure. Where an operation can be identified, diagnostics give its 1-based packed-op index or the 1-based byte at which the text operation starts. No read identifier is inferred. A top-level SQL NULL argument, including strict, returns NULL. Empty text, '*' and an empty packed list return NULL in either policy.

SELECT cigar_query_length('5S90M5I');

SELECT cigar_query_length([84, 1440, 81]::UINTEGER[], TRUE);

Return the aligned query length from a CIGAR string, counting M, =, and X but excluding clips and insertions. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_aligned_query_length(cigar[, strict])

Returns:

BIGINT

See cigar_query_length for full-input validation and missing-input behavior.

SELECT cigar_aligned_query_length('5S90M5I');

Return the reference-consuming length from a CIGAR string, counting M, D, N, =, and X. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_reference_length(cigar[, strict])

Returns:

BIGINT

See cigar_query_length for full-input validation and missing-input behavior.

SELECT cigar_reference_length('90M5D');

Test whether a CIGAR string contains at least one instance of the requested operator. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_has_op(cigar, op[, strict])

Returns:

BOOLEAN

Uses the full-input validation of cigar_query_length, including the suffix after any matching op. The requested operator is one supported ASCII character, case-insensitive independently of the process locale.

The optional final positional BOOLEAN strict defaults to FALSE: an invalid operator or malformed CIGAR returns NULL. TRUE raises an error with the diagnostic conventions of cigar_query_length. A top-level SQL NULL argument, including strict, returns NULL. Empty text, '*' and an empty packed list return false for a valid requested operator in either policy.

SELECT cigar_has_op('5S90M5S', 'S');

Return the aligned blocks of a CIGAR as a struct of three parallel BIGINT lists: ref_start, query_start and width, one entry per M, = or X op in CIGAR order. Overloaded to also accept a UINTEGER[] binary CIGAR (as produced by read_bam(cigar_representation := 'binary')); the binary overload is bit-identical to the text path.

Signature:

cigar_aligned_blocks(cigar, pos[, strict])

Returns:

STRUCT

ref_start is pos plus the reference bases consumed before the block, so it carries whatever base pos uses; pass read_bam's 1-based POS for 1-based starts or 0 for offsets from the alignment start. query_start is the 0-based offset into the stored SEQ: soft clips count, hard clips do not. width is the op length. D and N advance the reference and split blocks, I advances the query and splits blocks, H and P consume nothing. Blocks are never merged, as in pysam get_blocks() and GenomicAlignments cigarRangesAlongReferenceSpace over M, = and X.

Uses the full-input CIGAR grammar of cigar_query_length. The half-open reference end, pos + cigar_reference_length(cigar), must also fit BIGINT. A negative pos is permitted. A valid CIGAR with no aligned op returns three empty lists.

The optional final positional BOOLEAN strict defaults to FALSE: malformed CIGAR or coordinate overflow returns NULL. TRUE raises an error with the diagnostic conventions of cigar_query_length. A top-level SQL NULL argument, including pos or strict, returns NULL. Empty text, '*' and an empty packed list return NULL in either policy.

SELECT (cigar_aligned_blocks('5S90M5S', 100)).ref_start;

SELECT UNNEST((b).ref_start) AS ref_start, UNNEST((b).width) AS width FROM (SELECT cigar_aligned_blocks(CIGAR, POS) AS b FROM read_bam('reads.bam', cigar_representation := 'binary'));

Test whether the SAM flag indicates that the template has multiple segments in sequencing (0x1).

Signature:

is_paired(flag)

Returns:

BOOLEAN

SELECT is_paired(99);

Test whether the SAM flag indicates that each segment is properly aligned according to the aligner (0x2).

Signature:

is_proper_pair(flag)

Returns:

BOOLEAN

SELECT is_proper_pair(99);

Test whether the read itself is unmapped according to the SAM flag.

Signature:

is_unmapped(flag)

Returns:

BOOLEAN

SELECT is_unmapped(4);

Test whether the next segment in the template is flagged as unmapped (0x8).

Signature:

is_next_segment_unmapped(flag)

Returns:

BOOLEAN

SELECT is_next_segment_unmapped(9);

Test whether SEQ is stored reverse complemented (0x10); for mapped reads this corresponds to reverse-strand alignment.

Signature:

is_reverse_complemented(flag)

Returns:

BOOLEAN

SELECT is_reverse_complemented(16);

Test whether SEQ of the next segment in the template is stored reverse complemented (0x20).

Signature:

is_next_segment_reverse_complemented(flag)

Returns:

BOOLEAN

SELECT is_next_segment_reverse_complemented(32);

Test whether the read is marked as the first segment in the template.

Signature:

is_first_segment(flag)

Returns:

BOOLEAN

SELECT is_first_segment(64);

Test whether the read is marked as the last segment in the template.

Signature:

is_last_segment(flag)

Returns:

BOOLEAN

SELECT is_last_segment(128);

Test whether the alignment is marked as secondary.

Signature:

is_secondary(flag)

Returns:

BOOLEAN

SELECT is_secondary(256);

Test whether the read failed vendor or pipeline quality checks.

Signature:

is_qc_fail(flag)

Returns:

BOOLEAN

SELECT is_qc_fail(512);

Test whether the alignment is flagged as a duplicate.

Signature:

is_duplicate(flag)

Returns:

BOOLEAN

SELECT is_duplicate(1024);

Test whether the alignment is marked as supplementary.

Signature:

is_supplementary(flag)

Returns:

BOOLEAN

SELECT is_supplementary(2048);

Any scripts or data that you put into this service are public.

Rduckhts documentation built on Oct. 25, 2026, 1:08 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

Rduckhts 'DuckDB' High Throughput Sequencing File Formats Reader Extension

inst/function_catalog/reference.md In Rduckhts: 'DuckDB' High Throughput Sequencing File Formats Reader Extension

Extension Function Reference

duckhts_macro_definitions

Installation

Columns

Examples

duckhts_htslib_version

Examples

duckhts_htslib_features

Examples

duckhts_htslib_feature_string

Examples

duckhts_simd_backend

Examples

duckhts_simd_requested_backend

Examples

duckhts_simd_backend_compiled

Examples

duckhts_simd_backend_cpu_supported

Examples

duckhts_simd_backend_available

Examples

duckhts_simd_info

Diagnostics

Examples

duckhts_simd_kernel_info

Examples

duckhts_simd_set_backend

Examples

duckhts_duckdb_type_supported

Examples

duckhts_duckdb_supports_variant

Examples

duckhts_duckdb_supports_geometry

Examples

read_bcf

Types

FILTER

Errors

Samples

Regions

Typed regions

Scanning

Offsets

Threads

Examples

read_geno

Output

Selection

Ordering

FORMAT

FILTER

Errors

Raw GT

Typed regions

Exact alleles

Examples

read_bcf_samples

Examples

read_bam

Map keys

Representation

Offsets

Scanning

Threads

CRAM fields

Limits

Examples

duckhts_bcf_convert_parquet_sql

Metadata

Execution

Examples

duckhts_bam_convert_parquet_sql

Metadata

Execution

Examples

duckhts_gff_convert_parquet_sql

Metadata

Execution

Rduckhts
'DuckDB' High Throughput Sequencing File Formats Reader Extension

inst/function_catalog/reference.md
In Rduckhts: 'DuckDB' High Throughput Sequencing File Formats Reader Extension