Search Shortcut cmd + k | ctrl + k
duckhts

DuckDB extension for reading HTS file formats via htslib

Maintainer(s): sounkou-bioinfo

Installing and Loading

INSTALL duckhts FROM community;
LOAD duckhts;

Example

-- Load the extension
LOAD duckhts;

-- Read a VCF/BCF file (tidy FORMAT columns)
SELECT CHROM, POS, REF, ALT, SAMPLE_ID
FROM read_bcf('test/data/formatcols.vcf.gz', tidy_format := true)
LIMIT 5;

-- Read a BAM/SAM file
SELECT QNAME, RNAME, POS, READ_GROUP_ID, SAMPLE_ID
FROM read_bam('test/data/rg.sam.gz')
LIMIT 5;

About duckhts

DuckHTS provides DuckDB table functions for common high-throughput sequencing (HTS) formats using htslib. Query VCF, BCF, BAM, CRAM, FASTA, FASTQ, GTF, GFF, and tabix-indexed files directly in SQL.

The extension also includes sequence utility UDFs, SAM flag predicate helpers, and HTS metadata helpers.

Functions included in this extension:

Diagnostics

  • duckhts_htslib_version(): 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.
  • duckhts_htslib_features(): Return the htslib runtime feature bitfield reported by hts_features(). Use duckhts_htslib_feature_string() for the corresponding build description.
  • duckhts_htslib_feature_string(): 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.
  • duckhts_simd_backend(): 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.
  • duckhts_simd_requested_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.
  • duckhts_simd_backend_compiled(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.
  • duckhts_simd_backend_cpu_supported(backend): 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.
  • duckhts_simd_backend_available(backend): 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.
  • duckhts_simd_info(): Return one row per known concrete DuckHTS SIMD backend with extension-owned selectable, compiled, CPU-supported, available, selected, requested, and dispatch-mode diagnostics. Availability is the intersection of compiled and CPU/runtime-supported. selectable reports whether the backend has a selectable implementation path; explicit selection still requires available = TRUE. selected is TRUE when the current dispatch table uses that backend for at least one logical kernel. auto is a selection request and is not a concrete backend row.
  • duckhts_simd_kernel_info(): 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.
  • duckhts_simd_set_backend(backend): 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.
  • duckhts_duckdb_type_supported(type_name): 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.
  • duckhts_duckdb_supports_variant(): Return whether the currently open DuckDB runtime advertises the VARIANT logical type. Use this to gate optional SQL that depends on DuckDB VARIANT support.
  • duckhts_duckdb_supports_geometry(): Return whether the currently open DuckDB runtime advertises the GEOMETRY logical type. Use this to gate optional SQL that depends on DuckDB GEOMETRY support.

Variant Annotation

  • duckvep_ensembl_regions(core_schema, reference_chunks_table, assembly, species_id := 1): Match a tiled reference sequence to one Ensembl core assembly and assign dense DuckVEP sequence-region ordinals. core_schema must contain the Ensembl seq_region and coord_system tables. reference_chunks_table must contain chrom, zero-based start, half-open end, and seq columns, normally materialized from fasta_nuc(…, include_seq := TRUE). Chunks must be contiguous from zero, their sequence lengths must match their intervals, and every FASTA contig must match exactly one same-length Ensembl region. Regions absent from the supplied FASTA are deliberately excluded. The supplied reference therefore defines the exact modeled path set: a primary-assembly FASTA produces a primary-path model, while alternate haplotypes or patches require their named sequence paths in the reference chunks. This macro does not apply Ensembl assembly_exception projections, merge X/Y PAR paths, or mark circular sequence regions; ordinary mitochondrial coordinates can be modeled, but an interval that wraps from the end of a circular sequence to its start is not accepted.
  • duckvep_ensembl_transcripts(core_schema, reference_chunks_table, assembly, species_id := 1): Build the validated DuckVEP transcript relation from Ensembl core tables and matching tiled FASTA sequence. This supplied builder compiles the VEP-116 Ensembl core selection; it does not build VEP –refseq or –merged transcript sets. Before dense ordinals are assigned, transcripts must be current and have a non-empty stable ID, and artifact-biotype or readthrough_tra transcripts are excluded. The result contains the exact resident-loader columns, stable/source identifiers, biotypes, versioned RefSeq accessions from MANE Select and MANE Plus Clinical attributes when present, a nested ranked-exon projection, mature miRNA cDNA attributes projected into genomic exon segments, and supported Translation SeqEdits. MANE accessions remain cold DuckDB columns while only selection flags enter the resident C model. GENCODE Basic and Primary membership is retained in transcript_flags for late filtering; this macro does not prefilter the model to either set. Prepared sequence contains the phase-adjusted CDS plus the complete transcript-oriented pre-CDS and post-CDS spliced sequence; the legacy post_cds_bases field remains as the first three post-CDS bases. Codon-table IDs come from the Ensembl seq_region_attrib relation and default to table 1 exactly as VEP 116 does; every BioPerl/VEP-supported NCBI table ID is accepted and invalid or conflicting source attributes reject the import. Single-residue initial_met, _selenocysteine, amino_acid_sub, and _stop_codon_rt edits are retained as a sparse reference-peptide overlay. Transcript-level sequence corrections, including inserted or deleted transcript bases relative to the genome, and other Translation SeqEdit shapes fail closed by withholding sequence with an explicit reason. Transcript rows remain specific to the exact selected sequence path, including separate X/Y PAR and alternate-haplotype stable IDs. The macro does not collapse alt_allele equivalence groups, use IS_PAR metadata to merge consequence rows, or lift transcripts across a circular sequence origin.
  • duckvep_ensembl_regulation_features(funcgen_schema, regions_table): Build the VEP-116-compatible interval relation for Ensembl RegulatoryFeature and MotifFeature rows on the sequence regions selected for a DuckVEP model. funcgen_schema must contain regulatory_feature, feature_type, and motif_feature; regions_table is the output of duckvep_ensembl_regions(…). Matching VEP's database and cache source, epigenetically_modified_region (EMAR) RegulatoryFeature rows are excluded before dense ordinals are assigned. The result maps source sequence-region IDs to model ordinals, validates one-based inclusive coordinates and source feature types, and preserves stable/source IDs, feature metadata, binding-matrix IDs, regulatory-build IDs, and motif scores as cold relation columns. feature_kind is the compact resident code: 1 regulatory region or 2 transcription-factor binding site. Pass the five hot columns regulation_feature_index, seq_region, feature_start, feature_end, and feature_kind to duckvep_model_load(interval_feature_query := …); the remaining columns stay relational for late projection.
  • duckvep_model_receipt(regions_table, transcripts_table, source_name, source_version, assembly, source_manifest_sha256, reference_sha256, transcript_filter, regulation_features_table := NULL): Create one deterministic receipt for prepared DuckVEP region and transcript relations plus an optional prepared regulation-feature relation. The original eight positional arguments remain valid. When regulation_features_table is supplied by name, the receipt also validates its dense ordinals, region agreement, interval geometry, and feature kinds; counts regulatory regions and motif features; and includes their five resident columns in the model hash. It returns the declared provenance as source_name, source_version, assembly, source_manifest_sha256, reference_sha256, and transcript_filter; counts transcript, sequence, mature-miRNA, and peptide-edit content; and hashes every semantic resident-model field in stable ordinal order. source_manifest_sha256 should identify a canonical manifest containing every exact core and funcgen input used. The declared transcript_filter or source manifest should also record whether the reference contains only primary paths or includes alternate/patch paths and whether assembly_exception preprocessing was applied. The receipt contains no clock time, so rebuilding identical inputs yields the same model hash.
  • duckvep_model_load(name, sequence_region_query, transcript_query, exon_query, mature_mirna_query := NULL, peptide_edit_query := NULL, interval_feature_query := NULL, reference_fasta := NULL, transcript_coverage_complete := FALSE): Load one immutable consequence model into the current DuckDB database instance under a caller-chosen name. The resident kernel is source-neutral: callers may prepare any transcript catalog that satisfies this model contract. Three required query strings read committed, non-temporary DuckDB relations: sorted UINTEGER sequence-region ordinals; dense transcript rows with genomic span, strand, gene ordinal, flags, optional CDS span/sequence/codon table; and transcript-ordered exon rows with genomic and cDNA spans plus phase. Each exon cDNA span must be contiguous in transcript order and the same length as its genomic exon span. Transcript models that insert or delete bases relative to the genome therefore need a future explicit mapping and are not accepted by this interface. The optional mature_mirna_query reads transcript_index UINTEGER, mature_mirna_start UBIGINT, and mature_mirna_end UBIGINT ordered by transcript and start. The optional peptide_edit_query reads transcript_index UINTEGER, protein_position UINTEGER, and alternate_amino_acid VARCHAR ordered uniquely by transcript and protein position. The optional interval_feature_query reads regulation_feature_index UINTEGER, seq_region UINTEGER, feature_start UINTEGER, feature_end UINTEGER, and feature_kind UTINYINT ordered by region, start, and index; kind 1 is RegulatoryFeature and kind 2 is MotifFeature. Loading keeps those hot columns in a compact SoA and builds a separate cgranges seed index, while identifiers and funcgen metadata remain ordinary DuckDB columns. The transcript projection accepts 11 base columns, 12 columns with the legacy post_cds_bases BLOB of up to three bases, or 13 columns ending in complete transcript-oriented pre_cds_sequence and post_cds_sequence BLOBs. Complete flanks enable VEP start/stop predicates for length-changing edits that cross a CDS start or end. reference_fasta enables reference-validated, VEP-style 3-prime HGVS shifting. When supplied, the sequence-region query must return seq_region UINTEGER, sequence_length UBIGINT, and seq_region_name VARCHAR; loading requires an existing FASTA index and verifies every declared name and exact length without creating or modifying the index. The named model pins open read descriptors for the validated FASTA, .fai, and optional .gzi. Linux workers reopen those descriptors through /proc/self/fd, Windows retains the resolved source under deny-write sharing, and other POSIX workers use independent resolved-source handles with identity checks rather than shared /dev/fd seek state. Each annotation worker owns its mutable faidx handle, reuses contained sorted reference requests, and rejects detectable replacement or in-place source mutation around a cache-miss fetch. Sequence-region ordinals are exact model paths: loading performs no contig-synonym, PAR, patch, haplotype, assembly_exception, or circular-origin projection. Callers must apply any declared path mapping before loading or constructing events. Explicit NULL for any optional query or reference parameter is equivalent to omission. The default is a partial transcript model: a variant with no loaded transcript is unresolved, not intergenic. Setting transcript_coverage_complete := TRUE requires sequence_length UBIGINT and permits supported intergenic results after coordinate bounds checks. Loading validates and narrows the model once and returns one TRUE row. Several named models may coexist.
  • duckvep_model_drop(name): Remove a named resident DuckVEP consequence model and release its transcript and regulation-feature interval indexes, sequences, and cached worker state. Returns FALSE when the name is absent or the model is in use by an annotation vector.
  • duckvep_allele_geometry(position, reference, alternate): Interpret one literal biallelic allele with the same lossless geometry authority used by DuckVEP. position is positive and one-based; REF and ALT are non-empty A/C/G/T/N strings of at most 65,535 bases and must differ after ASCII case folding. kind_code is stable: 0 SNV, 1 insertion, 2 deletion, 3 length-changing replacement, and 4 multi-base substitution. anchor_side_code is 0 for no retained anchor, 1 for a left anchor, and 2 for a right anchor. All coordinate fields are non-null zero-based half-open intervals except insertion_boundary0, which is non-null only for an interbase insertion. Difference offsets and lengths are non-null indices into the uploaded REF and ALT. The result separates the uploaded REF span, VEP-116 VariationFeature span, and fully minimized semantic reference edit. A pure insertion has equal feature/edit start and end plus insertion_boundary0. This function does not split multiallelic records or left-align alleles: run duckhts_bcftools_norm(…, split_multiallelic := TRUE) first when canonical exact keys or cross-source joins require normalization. Use feature_start0/feature_end0 for VEP consequence feature overlap, edit_start0/edit_end0 for affected-reference supplementary providers, and an explicit insertion-flank policy when interbase is TRUE; the uploaded raw REF span is not a biological affected interval.
  • duckvep_annotate(events_table, model_name, hgvs := FALSE, upstream_distance := 5000, downstream_distance := 5000, rich := FALSE): Annotate a narrow canonical event relation against one resident DuckVEP model through one fixed-schema, pipe-friendly surface. events_table names a coordinate-ordered relation with event_index UBIGINT, seq_region UINTEGER, positive one-based position, reference VARCHAR, alternate VARCHAR, end_position UBIGINT, structural_type VARCHAR, copy_change VARCHAR, mate_seq_region UINTEGER, and mate_position UBIGINT. One row represents one ALT allele. Keep genotype, phase-set, confidence-interval, raw-ALT, orientation, and other provenance in the source relation; join selected annotations back by event_index instead of multiplying wide payloads by transcript count. Literal reference/alternate rows are small variants; a typed end/span or supported symbolic ALT is structural; a complete mate coordinate pair is a breakend. Contradictory or incomplete geometry fails instead of being guessed, including disagreement between supported symbolic ALT and explicit structural_type values. For VCF span SVs, prepare position as POS + 1 after removing the left anchor; insertion position is the interbase site after POS. DEL implies LOSS and DUP/TDUP imply GAIN when copy_change is NULL. Filtering the globally coordinate-ordered event stream preserves the order required by each native transcript and regulation-feature sweep; the macro does not hide a sort. The event seq_region is an exact ordinal in the loaded model; annotation does not strip contig prefixes or project PAR, patch, alternate-haplotype, or circular-origin coordinates. X/Y PAR and alternate-path transcript rows therefore remain path-specific exactly as prepared. Each event is annotated independently: genotype and phase-set columns are not combined into haplotypes. The stable compact contract includes SO/region masks, IMPACT/status/reason codes, one-based cDNA/CDS/protein positions, amino-acid bytes, NMD prediction/escape codes, and overlap-object ordinals. protein_position is the affected one-based protein residue/codon ordinal used by codon-level evidence such as PS1/PM5. rich := TRUE adds Consequence, IMPACT, region, amino-acid, NMD, overlap-object, duckvep_status, and duckvep_reason text from the same native pass while retaining every compact field; FALSE or explicit NULL leaves those presentation columns NULL. hgvs := TRUE fuses independent-event HGVSc/HGVSn/HGVSp into the small-variant lane without repeating candidate discovery; it composes with rich := TRUE in one pass; structural and breakend HGVS remain NULL. DuckHTS audit status/reason fields distinguish an unresolved computation from a genuinely absent consequence and are not Ensembl CSQ fields. Zero disables either directional transcript window. Decode consequence_mask only after filtering by joining duckvep_so_terms(); join transcript_index, gene_index, or regulation_feature_index to the prepared model relations for stable identifiers, biotype, MANE/GENCODE/CCDS flags, and other presentation metadata. This typed rich projection does not claim the complete serialized VEP CSQ vocabulary: exon/intron ordinals, codon strings, distance text, CANONICAL/APPRIS/TSL, existing variation, and supplementary predictors require their declared model or annotation relations.
  • duckvep_so_terms(): Return the generated Sequence Ontology metadata behind DuckVEP's stable consequence-mask bits: bit index, single-bit mask, term, impact code and label, VEP severity rank, and evaluator tier. Filter and aggregate compact annotation rows numerically first, then expand only selected mask bits with (annotation.consequence_mask & term.consequence_mask) <> 0. The relation is generated from the same pinned VEP-116 class-model metadata as the native consequence engine; it is not a second hand-maintained term inventory.

Readers

  • 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'): Read VCF and BCF variant data with typed INFO, FORMAT, typed CSQ/ANN/BCSQ subfields, optional tidy sample output, optional bcftools-style CSQ type overrides, explicit scan_mode control ('auto' or full-file 'sequential'), optional htslib decompression worker threads via decompression_threads (default 0 for single-threaded reads), and decode_error_policy ('null', 'warn', or 'error') for corrupt BCF header-vs-payload decode mismatches. Comma-separated indexed regions use htslib's native multi-region iterator and emit a record once when requested regions overlap.
  • read_bcf_v2(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, samples_file := NULL, exclude_samples := NULL, exclude_samples_file := NULL, force_samples := FALSE, include_info := TRUE, include_format := TRUE, include_vep := TRUE, info_fields := NULL, format_fields := NULL, vep_fields := NULL): Experimental clone of read_bcf(…) for non-destructive benchmarking of current-schema scan optimizations. It preserves the read_bcf wide/tidy output schema by default, streams full-file scans by default instead of auto-iterating indexed contigs, accepts the same scan_mode and decode_error_policy ('null', 'warn', or 'error') contract as read_bcf, and shares read_bcf's native deduplicating multi-region iterator. It sets htslib VCF text max_unpack from DuckDB projection pushdown so unprojected INFO/FORMAT text can be skipped earlier, keeps INFO/FORMAT/GT decode buffers and parsed VEP/CSQ records in local scan state across output chunks, repeats VEP-derived values correctly across tidy sample rows, supports v2-only htslib sample pushdown via comma-separated samples/exclude_samples strings or samples_file/exclude_samples_file lists, can bind only requested INFO/FORMAT/VEP fields through include_info/include_format/include_vep and comma-separated info_fields/format_fields/vep_fields controls, rejects format_fields when no samples are selected, treats tidy_format := true with include_format := false as one row per variant because no SAMPLE_ID/FORMAT columns are bound, and counts VCF text rows without parsing when no columns are projected.
  • read_bcf_appender(path, target_table, region := NULL, tidy_format := TRUE, overwrite := TRUE, include_file_offset := FALSE, region_threads := 1, decompression_threads := 0, decode_error_policy := 'null'): Experimental side-effecting benchmark helper that reads BCF/VCF with htslib and appends a narrow scalar tidy table through DuckDB's chunk appender API. The target table is created in the default schema with columns CHROM, POS, REF, ALT, SAMPLE_ID, and FORMAT_GT; set include_file_offset := TRUE to add FILE_OFFSET as a file-order token. Comma-separated regions use htslib's native deduplicating multi-region iterators. region_threads controls at most 16 worker handles claiming primary-contig jobs iteratively: every contig keeps its complete requested interval set in one native iterator, so repeated, overlapping, or record-spanning intervals produce the same target schema and row multiset for every thread count. Header contig count does not create worker threads or open handles, and a request limited to one contig uses one iterator. Parallel row arrival is unordered; use ORDER BY FILE_OFFSET when file order matters. decompression_threads configures htslib workers per open region-worker handle, so total helper threads can scale as region workers times decompression workers. The operation runs inside an internal transaction and surfaces malformed-record errors as DuckDB errors; decode_error_policy controls corrupt BCF FORMAT/GT header-vs-payload mismatches with 'null', 'warn', or 'error'. target_table must be an unqualified SQL identifier. This is intended to compare appender-based materialization against read_bcf(…) query pipelines, not as a full replacement for the typed read_bcf surface.
  • 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): Read SAM, BAM, and CRAM alignments with optional typed SAMtags and auxiliary tag maps. Use sequence_encoding := 'nt16' to return SEQ as UTINYINT[], quality_representation := 'phred' to return QUAL as UTINYINT[], and cigar_representation := 'binary' to return packed BAM CIGAR operations as UINTEGER[] instead of SAM text. scan_mode := 'sequential' forces a full-file streaming scan instead of index-backed count/parallel paths and is incompatible with region. decompression_threads controls per-file htslib worker threads and defaults to 2; use 0 to disable worker threads.
  • read_fasta(path, region := NULL, index_path := NULL, gzi_path := NULL, sequence_encoding := 'string', scan_mode := 'auto'): Read FASTA records or indexed FASTA regions as sequence rows. Use sequence_encoding := 'nt16' to return SEQUENCE as UTINYINT[] (htslib nt16 4-bit codes) instead of VARCHAR. For bgzipped FASTA, gzi_path may point to an explicit .gzi sidecar when it is not colocated with the FASTA. scan_mode := 'sequential' forces full-file streaming/counting instead of index-backed count paths and is incompatible with region.
  • read_bed(path, region := NULL, index_path := NULL, scan_mode := 'auto'): Read BED3-BED12 interval files with canonical typed columns and optional tabix-backed region filtering. scan_mode := 'sequential' forces full-file streaming/counting instead of index-backed count paths and is incompatible with region.
  • fasta_nuc(path, bed_path := NULL, bin_width := NULL, region := NULL, index_path := NULL, gzi_path := NULL, bed_index_path := NULL, include_seq := FALSE): Compute bedtools nuc-style nucleotide composition for supplied BED intervals or generated fixed-width bins over a FASTA reference. For bgzipped FASTA, gzi_path may point to an explicit .gzi sidecar when it is not colocated with the FASTA.
  • read_fastq(path, interleaved := FALSE, mate_path := NULL, sequence_encoding := 'string', quality_representation := 'string', input_quality_encoding := 'phred33', scan_mode := 'auto'): 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.
  • read_bigwig(path, region := NULL, blocks_per_iteration := 64): Read stored BigWig signal intervals as CHROM, START0, END0, and VALUE. Coordinates are zero-based half-open. region uses htslib's one-based inclusive syntax; comma-separated requests are merged per contig and emit each stored interval once. Full scans claim nonempty contigs across DuckDB workers, and multi-region scans claim merged ranges. Each worker owns its mutable file handle and iterator. blocks_per_iteration controls libBigWig iterator batching, not the number of DuckDB workers. Local, native remote, and browser wasm reads all use DuckHTS's htslib hFILE transport rather than a second libcurl stack.
  • read_gff(path, header_names := NULL, header := FALSE, column_types := NULL, auto_detect := FALSE, attributes_map := FALSE, attributes_list := FALSE, attributes_pairs := FALSE, strict := FALSE, region := NULL, index_path := NULL, scan_mode := 'auto'): Read GFF annotations with optional raw scalar and richer list/pair parsed attribute columns, strict GFF3 structural validation, and indexed region filtering. Comma-separated indexed regions use htslib's native multi-region iterator and emit a row once when requested regions overlap. scan_mode := 'sequential' forces full-file streaming/counting instead of index-backed count paths and is incompatible with region.
  • read_gtf(path, header_names := NULL, header := FALSE, column_types := NULL, auto_detect := FALSE, attributes_map := FALSE, attributes_list := FALSE, attributes_pairs := FALSE, region := NULL, index_path := NULL, scan_mode := 'auto'): Read GTF annotations with optional raw scalar and richer list/pair parsed attribute columns and indexed region filtering. Comma-separated indexed regions use htslib's native multi-region iterator and emit a row once when requested regions overlap. scan_mode := 'sequential' forces full-file streaming/counting instead of index-backed count paths and is incompatible with region.
  • read_tabix(path, header_names := NULL, header := FALSE, column_types := NULL, auto_detect := FALSE, region := NULL, index_path := NULL, scan_mode := 'auto'): Read generic tabix-indexed text data with optional header handling and type inference. Comma-separated indexed regions use htslib's native multi-region iterator and emit a row once when requested regions overlap. scan_mode := 'sequential' forces full-file streaming/counting instead of index-backed count paths and is incompatible with region.
  • fasta_index(path, index_path := NULL): Build a FASTA index (.fai) and return a single row with columns success (BOOLEAN) and index_path (VARCHAR).
  • hts_union_query(reader, pattern, params := ''): 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.
  • hts_region_union_query(reader, path, regions, params := ''): Generate a UNION ALL BY NAME query string that reads one HTS file through separate per-region table-function scans. The result adds filename, duckhts_region_shard_id, and duckhts_region_shard columns so benchmark queries can compare region-sharded execution and explicit ordered writes. Pass regions as a DuckDB list of region strings; params is appended to each reader call and should not include region. The generated UNION ALL does not deduplicate records; adjacent or overlapping shards can duplicate BAM alignments or BCF/VCF records that span shard boundaries, so add explicit shard-local filters or downstream deduplication when exact once-only semantics are required.

Converters

  • 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'): Build a DuckDB COPY statement that converts read_bcf(…) output to Parquet with DuckHTS key-value metadata, preserved or corrected VCF header text, user metadata from metadata := map(…), selected columns, filters, and partition columns. Optional metadata_json_file is caller-managed and requires DuckDB's json extension to be available when the builder is called; use metadata maps for CRAN/offline-safe workflows. Clients execute the returned SQL string; the R wrapper is a thin DBI helper that does this execution.
  • 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'): Build a DuckDB COPY statement that converts read_bam(…) output to Parquet with DuckHTS key-value metadata, preserved or corrected SAM header text, user metadata from metadata := map(…), selected columns, filters, and partition columns. Optional metadata_json_file is caller-managed and requires DuckDB's json extension to be available when the builder is called; use metadata maps for CRAN/offline-safe workflows. Clients execute the returned SQL string; the R wrapper is a thin DBI helper that does this execution.
  • 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'): Build a DuckDB COPY statement that converts read_gff(…) output to Parquet with DuckHTS key-value metadata, preserved or corrected GFF/tabix header text, user metadata from metadata := map(…), selected columns, filters, and partition columns. Optional metadata_json_file is caller-managed and requires DuckDB's json extension to be available when the builder is called; use metadata maps for CRAN/offline-safe workflows. Clients execute the returned SQL string; the R wrapper is a thin DBI helper that does this execution.
  • 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'): Build a DuckDB COPY statement that converts read_tabix(…) output to Parquet with DuckHTS key-value metadata, preserved or corrected tabix header text, user metadata from metadata := map(…), selected columns, filters, and partition columns. Optional metadata_json_file is caller-managed and requires DuckDB's json extension to be available when the builder is called; use metadata maps for CRAN/offline-safe workflows. Clients execute the returned SQL string; the R wrapper is a thin DBI helper that does this execution.

Coverage

  • read_pileup(path, region := NULL, index_path := NULL, min_mapq := 0, flag_mask := 1796): 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.
  • bam_bin_counts(path, bin_width, chrom := NULL, reference := NULL, index_path := NULL, mapq := 0, require_flags := 0, exclude_flags := 0, rmdup := 'none', stats := NULL): 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.
  • 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): 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.
  • 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): Write native mosdepth-compatible coverage outputs for indexed BAM or CRAM input. Produces mosdepth-style summary, global distribution, per-base BED.gz + CSI, optional window/BED region outputs, optional quantized BED.gz + CSI, and optional threshold counts for by; fast_mode defaults to FALSE to match upstream mosdepth, default mode performs CIGAR-aware coverage with mate-overlap correction, fragment_mode switches coverage to full-fragment insert spans for proper pairs, use_median switches by outputs from mean to median, read_groups filters by comma-separated RG IDs, min_frag_len and max_frag_len filter on absolute template length, fasta is required for CRAM when htslib needs a reference, precision_digits controls decimal places in the text outputs, and processing_threads enables parallel contig processing (0 = sequential, >0 = number of worker threads).

Intervals

  • duckhts_cgranges_create(name): Create an empty session-scoped cgranges registry entry that can be populated with intervals and finalized for overlap queries.
  • duckhts_cgranges_add(name, chrom, start, end[, label]): Append an interval to a session-scoped cgranges registry entry before finalization. Labels may be BIGINT-like, DOUBLE, VARCHAR, or BOOLEAN.
  • duckhts_cgranges_index(name): Finalize a populated cgranges registry entry and build its immutable overlap index for subsequent queries.
  • duckhts_cgranges_destroy(name): Destroy a session-scoped cgranges registry entry and release its indexed interval storage when it is not in active use.
  • duckhts_cgranges_from_query(name, query, chrom_col, start_col, end_col[, label_col]): Execute a SQL query on an extension-owned DuckDB connection, append its interval rows into a session-scoped cgranges registry entry, and leave the populated index ready for explicit finalization with duckhts_cgranges_index(…).
  • duckhts_cgranges_from_table(name, table_name, chrom_col, start_col, end_col[, label_col]): Reserved convenience constructor for bulk cgranges population from a table name. The current implementation is intentionally deferred and directs callers to duckhts_cgranges_from_query(…).
  • duckhts_cgranges_has_overlap(name, chrom, start, end[, mode]): 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.
  • duckhts_cgranges_count_overlaps(name, chrom, start, end[, mode]): 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.
  • duckhts_cgranges_overlaps_list(name, chrom, start, end[, mode]): 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.
  • duckhts_cgranges_overlaps(name, chrom, start, end, mode := 'overlap', query_row_id := NULL): 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.
  • duckhts_cgranges_overlaps_bulk(name, query, chrom_col, start_col, end_col, mode := 'overlap', query_row_id_col := NULL): Run a SQL query that yields overlap probes, stream those rows through a finalized session-scoped cgranges registry entry, and return one row per matching indexed interval. The probe query runs on the extension-owned helper connection, so it must reference regular tables/views rather than connection-local temp tables. When query_row_id_col is omitted, query_row_id defaults to the 1-based probe row ordinal.
  • regionkey(chrom, start, end, strand := 0): 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.
  • regionkey_hex(rk): Render a RegionKey as its lowercase 16-character hexadecimal string representation.
  • parse_regionkey_hex(hex): Parse a 16-character hexadecimal RegionKey string back into its UBIGINT code. Invalid or non-hex strings return NULL.
  • encode_regionkey(chrom_code, start, end, strand_code): Encode the raw upstream RegionKey fields directly: chromosome code, 0-based start, 0-based end, and strand code (0 = unknown, 1 = +, 2 = -).
  • extract_regionkey_chrom(rk): Extract the raw upstream RegionKey chromosome code.
  • extract_regionkey_startpos(rk): Extract the raw upstream RegionKey 0-based start position.
  • extract_regionkey_endpos(rk): Extract the raw upstream RegionKey 0-based end position.
  • extract_regionkey_strand(rk): Extract the raw upstream RegionKey strand code (0 = unknown, 1 = +, 2 = -).
  • decode_regionkey(rk): Decode a RegionKey into its raw upstream numeric fields: chrom_code, start, end, and strand_code.
  • reverse_regionkey(rk): Decode a RegionKey into a STRUCT with chrom, chrom_code, start, end, strand, and strand_code.
  • extend_regionkey(rk, size): Extend a RegionKey interval by a fixed number of bases on both sides, clamping to the official 28-bit RegionKey position range.
  • are_overlapping_regions(chrom_a, start_a, end_a, chrom_b, start_b, end_b): Return TRUE when two explicit 0-based half-open intervals overlap on the same canonical chromosome.
  • are_overlapping_region_regionkey(chrom, start, end, rk): Return TRUE when a 0-based half-open interval overlaps the supplied RegionKey interval.
  • are_overlapping_regionkeys(rka, rkb): Return TRUE when two RegionKeys overlap.

Quality Control

  • duckhts_fastq_qc(sequence, quality [, max_cycles]): 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.

Metadata

  • detect_quality_encoding(path, max_records := 10000): Inspect a FASTQ file's observed quality ASCII range and report compatible legacy encodings with a heuristic guessed encoding.
  • duckhts_samtools_idxstats(path, output := NULL, index_path := NULL, threads := 0, overwrite := FALSE): 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.
  • read_hts_header(path, format := NULL, mode := NULL): 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.
  • read_hts_index(path, format := NULL, index_path := NULL): Inspect high-level HTS index metadata such as sequence names and mapped counts.
  • read_hts_index_spans(path, format := NULL, index_path := NULL): Expand index metadata into span and chunk rows suitable for low-level index inspection.
  • read_hts_index_raw(path, format := NULL, index_path := NULL): Return the raw on-disk HTS index blob together with basic identifying metadata.

Compression

  • bgzip(path, output_path := NULL, threads := 4, level := -1, keep := TRUE, overwrite := FALSE): Compress a plain file to BGZF and return the created output path and byte counts.
  • bgunzip(path, output_path := NULL, threads := 4, keep := TRUE, overwrite := FALSE): Decompress a BGZF-compressed file and return the created output path and byte counts.

Indexing

  • bam_index(path, index_path := NULL, min_shift := 0, threads := 4): Build a BAM or CRAM index and report the written index path and format.
  • bcf_index(path, index_path := NULL, min_shift := NULL, threads := 4): Build a TBI or CSI index for a VCF or BCF file and report the written index path and format.
  • 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): Build a tabix index for a BGZF-compressed text file using a preset or explicit coordinate columns.

Variants

  • variantkey(chrom, pos, ref, alt): 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.
  • variantkey_hex(vk): Render a VariantKey as its lowercase 16-character hexadecimal string representation.
  • parse_variantkey_hex(hex): Parse a 16-character hexadecimal VariantKey string back into its UBIGINT code. Invalid or non-hex strings return NULL.
  • encode_variantkey(chrom_code, pos0, refalt_code): Encode the raw upstream VariantKey fields directly: chromosome code, 0-based position, and 31-bit REF+ALT code.
  • extract_variantkey_chrom(vk): Extract the raw upstream VariantKey chromosome code.
  • extract_variantkey_pos(vk): Extract the raw upstream VariantKey 0-based position field.
  • extract_variantkey_refalt(vk): Extract the raw upstream 31-bit VariantKey REF+ALT code.
  • decode_variantkey(vk): Decode a VariantKey into its raw upstream numeric fields: chrom_code, pos0, and refalt_code.
  • reverse_variantkey(vk): 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.
  • variantkey_range(chrom, pos_min, pos_max): 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.
  • duckhts_contig_key(contig): 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.
  • 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): 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).
  • duckdb_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): DuckDB-specific wrapper over bcftools_liftover that takes either a table name or a derived-table expression plus column-name strings for chrom/pos/ref/alt and returns the lifted table. The no_left_align parameter mirrors –no-left-align in bcftools +liftover.
  • bcftools_norm_row(chrom, pos, ref, alt, fasta_ref, end_pos := NULL, svlen := NULL, fasta_index_path := NULL, gzi_path := NULL): Normalize one variant row with bcftools/vt-style left-alignment semantics against a FASTA reference. The alt argument may be either a comma-delimited VARCHAR or a VARCHAR[] list. The returned STRUCT contains pos_normed, end_pos_normed, ref_normed, alt_normed (always VARCHAR[]), normed (TRUE/FALSE/NULL), and norm_status. Symbolic rows can use end_pos, and symbolic rows can use svlen. gVCF /<*> reference-block alleles pass through with GVCFReferenceBlock, and mixed real-plus-gVCF-symbolic alleles normalize the real alleles while preserving symbolic alleles and caller-supplied reference-block END in site-preserving output.
  • 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): DuckDB table macro wrapper over bcftools_norm_row that normalizes variants from a table or derived-table expression while preserving the original columns. The input ALT column may be either VARCHAR or VARCHAR[]. The result appends pos_normed, end_pos_normed, ref_normed, alt_normed, normed, and norm_status; with split_multiallelic := TRUE, multiallelic sites are split before normalization and alt_normed becomes VARCHAR plus alt_index. This is a vt/vcfnorm-style row/table transform, not a full-record VCF/BCF rewrite: genotype, PL/GP/DS, and PS fields are preserved as caller columns unless a separate full-record writer/remapper layer is used.
  • 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): Compute polygenic scores from one genotype BCF/VCF and one or more summary-statistics files with bcftools +score-compatible GT/DS/HDS/AP/GP/AS dosage semantics, sample subsetting, and region/target/FILTER-string controls. The second argument accepts a scalar path or a DuckDB LIST/array of paths; TSV/SSF summaries produce one PRS column per file in a single genotype scan, while GWAS-VCF summaries still produce one PRS column per FORMAT sample. Use summaries_list_file with a NULL second argument to read paths from a file or directory; list-file entries are interpreted as written, matching upstream bcftools +score --summaries behavior, while directory inputs scan supported regular summary files in lexicographic order and ignore index sidecars. Use log_path to write per-PRS loaded/matched/allele-mismatch/duplicate-marker audit counts.
  • 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): 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.
  • duckdb_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): DuckDB macro wrapper over bcftools_munge_row that maps source columns (via preset or explicit map) and returns normalized GWAS-VCF-style rows with lean outputs and 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 duckdb_munge_metal.
  • duckdb_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): Extended munge macro with METAL meta-analysis output columns. Same as duckdb_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). The R wrapper 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.

Sequence UDFs

  • seq_revcomp(sequence): 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.
  • seq_canonical(sequence): 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.
  • seq_hash_2bit(sequence): 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.
  • seq_encode_4bit(sequence): Encode an IUPAC DNA sequence as a list of 4-bit base codes, preserving ambiguity symbols including N.
  • seq_decode_4bit(codes): Decode a list of 4-bit IUPAC DNA base codes back into a sequence string.
  • seq_gc_content(sequence): 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.
  • seq_kmers(sequence, k, canonical := FALSE): Expand a sequence into positional k-mers with optional canonicalization.

SAM Flag UDFs

  • sam_flag_bits(flag): 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.
  • sam_flag_has(flag, mask): Test whether any bits from the provided SAM flag mask are set in a flag value.
  • is_forward_aligned(flag): 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.
  • is_paired(flag): Test whether the SAM flag indicates that the template has multiple segments in sequencing (0x1).
  • is_proper_pair(flag): Test whether the SAM flag indicates that each segment is properly aligned according to the aligner (0x2).
  • is_unmapped(flag): Test whether the read itself is unmapped according to the SAM flag.
  • is_next_segment_unmapped(flag): Test whether the next segment in the template is flagged as unmapped (0x8).
  • is_reverse_complemented(flag): Test whether SEQ is stored reverse complemented (0x10); for mapped reads this corresponds to reverse-strand alignment.
  • is_next_segment_reverse_complemented(flag): Test whether SEQ of the next segment in the template is stored reverse complemented (0x20).
  • is_first_segment(flag): Test whether the read is marked as the first segment in the template.
  • is_last_segment(flag): Test whether the read is marked as the last segment in the template.
  • is_secondary(flag): Test whether the alignment is marked as secondary.
  • is_qc_fail(flag): Test whether the read failed vendor or pipeline quality checks.
  • is_duplicate(flag): Test whether the alignment is flagged as a duplicate.
  • is_supplementary(flag): Test whether the alignment is marked as supplementary.

CIGAR Utils

  • cigar_has_soft_clip(cigar): 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.
  • cigar_has_hard_clip(cigar): 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.
  • cigar_left_soft_clip(cigar): 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.
  • cigar_right_soft_clip(cigar): 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.
  • cigar_query_length(cigar): 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.
  • cigar_aligned_query_length(cigar): 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.
  • cigar_reference_length(cigar): 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.
  • cigar_has_op(cigar, op): 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.

Operational notes:

  • Reader table functions accept scan_mode := 'sequential' to force full-file streaming/counting instead of index-backed count/parallel paths where applicable; region queries remain index-backed and are incompatible with sequential mode.
  • Paired FASTQ is supported via mate_path or interleaved := true.
  • CRAM reads are supported with an explicit reference file.
  • GTF/GFF attributes can be returned as a parsed MAP with attributes_map := true via read_gff or read_gtf.
  • Optional SAM tag columns and an auxiliary tag map are available through standard_tags and auxiliary_tags.
  • Tabix readers support header, header_names, type inference with auto_detect, and explicit column_types.
  • MSVC builds (windows_amd64/windows_arm64) are not supported; use MinGW/RTools on Windows.

Added Functions

function_name function_type description comment examples
_duckvep_annotate_breakend_compact scalar NULL NULL  
_duckvep_annotate_breakend_rich scalar NULL NULL  
_duckvep_annotate_small_compact scalar NULL NULL  
_duckvep_annotate_small_hgvs scalar NULL NULL  
_duckvep_annotate_small_rich scalar NULL NULL  
_duckvep_annotate_small_rich_hgvs scalar NULL NULL  
_duckvep_annotate_structural_compact scalar NULL NULL  
_duckvep_annotate_structural_rich scalar NULL NULL  
are_overlapping_region_regionkey scalar NULL NULL  
are_overlapping_regionkeys scalar NULL NULL  
are_overlapping_regions scalar NULL NULL  
bam_bin_counts table NULL NULL  
bam_index table NULL NULL  
bcf_index table NULL NULL  
bcftools_liftover scalar NULL NULL  
bcftools_munge_row scalar NULL NULL  
bcftools_norm_row scalar NULL NULL  
bcftools_score table NULL NULL  
bgunzip table NULL NULL  
bgzip table NULL NULL  
cigar_aligned_query_length scalar NULL NULL  
cigar_has_hard_clip scalar NULL NULL  
cigar_has_op scalar NULL NULL  
cigar_has_soft_clip scalar NULL NULL  
cigar_left_soft_clip scalar NULL NULL  
cigar_query_length scalar NULL NULL  
cigar_reference_length scalar NULL NULL  
cigar_right_soft_clip scalar NULL NULL  
decode_regionkey scalar NULL NULL  
decode_variantkey scalar NULL NULL  
detect_quality_encoding table NULL NULL  
duckdb_liftover table_macro NULL NULL  
duckdb_munge table_macro NULL NULL  
duckdb_munge_metal table_macro NULL NULL  
duckdb_munge_preset_map macro NULL NULL  
duckdb_munge_resolved_map macro NULL NULL  
duckhts_alt_to_list scalar NULL NULL  
duckhts_bam_bed_coverage table NULL NULL  
duckhts_bam_convert_parquet_sql macro NULL NULL  
duckhts_bcf_convert_parquet_sql macro NULL NULL  
duckhts_bcftools_norm table_macro NULL NULL  
duckhts_cgranges_add scalar NULL NULL  
duckhts_cgranges_count_overlaps scalar NULL NULL  
duckhts_cgranges_create scalar NULL NULL  
duckhts_cgranges_destroy scalar NULL NULL  
duckhts_cgranges_from_query scalar NULL NULL  
duckhts_cgranges_from_table scalar NULL NULL  
duckhts_cgranges_has_overlap scalar NULL NULL  
duckhts_cgranges_index scalar NULL NULL  
duckhts_cgranges_overlaps table NULL NULL  
duckhts_cgranges_overlaps_bulk table NULL NULL  
duckhts_cgranges_overlaps_list scalar NULL NULL  
duckhts_contig_key scalar NULL NULL  
duckhts_duckdb_supports_geometry macro NULL NULL  
duckhts_duckdb_supports_variant macro NULL NULL  
duckhts_duckdb_type_supported macro NULL NULL  
duckhts_fastq_qc aggregate NULL NULL  
duckhts_gff_convert_parquet_sql macro NULL NULL  
duckhts_htslib_feature_string scalar NULL NULL  
duckhts_htslib_features scalar NULL NULL  
duckhts_htslib_version scalar NULL NULL  
duckhts_json_path_key macro NULL NULL  
duckhts_mosdepth table NULL NULL  
duckhts_parquet_copy_sql macro NULL NULL  
duckhts_parquet_ident_list macro NULL NULL  
duckhts_parquet_metadata_args macro NULL NULL  
duckhts_quote_ident macro NULL NULL  
duckhts_quote_string macro NULL NULL  
duckhts_raw_header_text macro NULL NULL  
duckhts_samtools_idxstats table NULL NULL  
duckhts_simd_backend scalar NULL NULL  
duckhts_simd_backend_available scalar NULL NULL  
duckhts_simd_backend_compiled scalar NULL NULL  
duckhts_simd_backend_cpu_supported scalar NULL NULL  
duckhts_simd_info table NULL NULL  
duckhts_simd_kernel_info table NULL NULL  
duckhts_simd_requested_backend scalar NULL NULL  
duckhts_simd_set_backend table NULL NULL  
duckhts_tabix_convert_parquet_sql macro NULL NULL  
duckvep_allele_geometry scalar NULL NULL  
duckvep_annotate table_macro NULL NULL  
duckvep_ensembl_regions table_macro NULL NULL  
duckvep_ensembl_regulation_features table_macro NULL NULL  
duckvep_ensembl_transcripts table_macro NULL NULL  
duckvep_model_drop scalar NULL NULL  
duckvep_model_load table NULL NULL  
duckvep_model_receipt table_macro NULL NULL  
duckvep_so_terms table NULL NULL  
encode_regionkey scalar NULL NULL  
encode_variantkey scalar NULL NULL  
extend_regionkey scalar NULL NULL  
extract_regionkey_chrom scalar NULL NULL  
extract_regionkey_endpos scalar NULL NULL  
extract_regionkey_startpos scalar NULL NULL  
extract_regionkey_strand scalar NULL NULL  
extract_variantkey_chrom scalar NULL NULL  
extract_variantkey_pos scalar NULL NULL  
extract_variantkey_refalt scalar NULL NULL  
fasta_index table NULL NULL  
fasta_nuc table NULL NULL  
hts_region_union_query macro NULL NULL  
hts_union_query macro NULL NULL  
is_duplicate scalar NULL NULL  
is_first_segment scalar NULL NULL  
is_forward_aligned scalar NULL NULL  
is_last_segment scalar NULL NULL  
is_next_segment_reverse_complemented scalar NULL NULL  
is_next_segment_unmapped scalar NULL NULL  
is_paired scalar NULL NULL  
is_proper_pair scalar NULL NULL  
is_qc_fail scalar NULL NULL  
is_reverse_complemented scalar NULL NULL  
is_secondary scalar NULL NULL  
is_supplementary scalar NULL NULL  
is_unmapped scalar NULL NULL  
parse_regionkey_hex scalar NULL NULL  
parse_variantkey_hex scalar NULL NULL  
read_bam table NULL NULL  
read_bcf table NULL NULL  
read_bcf_appender table NULL NULL  
read_bcf_v2 table NULL NULL  
read_bed table NULL NULL  
read_bigwig table NULL NULL  
read_fasta table NULL NULL  
read_fastq table NULL NULL  
read_gff table NULL NULL  
read_gtf table NULL NULL  
read_hts_header table NULL NULL  
read_hts_index table NULL NULL  
read_hts_index_raw table_macro NULL NULL  
read_hts_index_spans table NULL NULL  
read_pileup table NULL NULL  
read_tabix table NULL NULL  
regionkey scalar NULL NULL  
regionkey_hex scalar NULL NULL  
reverse_regionkey scalar NULL NULL  
reverse_variantkey scalar NULL NULL  
sam_flag_bits scalar NULL NULL  
sam_flag_has scalar NULL NULL  
seq_canonical scalar NULL NULL  
seq_decode_4bit scalar NULL NULL  
seq_encode_4bit scalar NULL NULL  
seq_gc_content scalar NULL NULL  
seq_hash_2bit scalar NULL NULL  
seq_kmers table NULL NULL  
seq_revcomp scalar NULL NULL  
tabix_index table NULL NULL  
variantkey scalar NULL NULL  
variantkey_hex scalar NULL NULL  
variantkey_range scalar NULL NULL  

Overloaded Functions

This extension does not add any function overloads.

Added Types

This extension does not add any types.

Added Settings

This extension does not add any settings.