rduckhts_roh: Find Runs of Homozygosity

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

Find Runs of Homozygosity

Description

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

Usage

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

Arguments

con	A DuckDB connection with DuckHTS loaded.
path	One VCF/BCF path or URI.
af_tag	INFO tag holding the ALT allele frequency, for example '"AF"'.
af_table	Name of a table or view of site allele frequencies, instead of 'af_tag' or the ancestry inputs.
genetic_map	Optional name of a genetic-map table or view.
hw_to_az, az_to_hw	Per-base-pair transition probabilities Hardy-Weinberg to autozygous ('-a') and back ('-H'), each in '[0, 1]'.
gt_error	Optional phred error for GT-only emissions ('-G'), at least 0.
rec_rate	Optional constant recombination rate per base pair ('-M').
samples	Optional HTSlib sample selector: comma-separated inclusion, leading '^' exclusion, '"-"' for all, or '""' for none.
table_name	Optional output table. 'NULL' returns a data frame.
overwrite	Whether an existing output table may be replaced.
reference_table	Long-format reference with chromosome, position, allele_a, allele_b, group_id and frequency columns.
proportions_table	Long-format sample ancestry proportions with sample_id, group_id and proportion columns. Its group set must exactly match 'reference_table'.
af_clamp	Clamp ancestry-tuned frequencies to '[af_clamp, 1-af_clamp]'; zero disables clamping. Must be in '[0, 0.5)'.
max_sites	Most sites one sample and chromosome may hold, from 1 to 100,000,000. A larger group is an error.
max_site_bytes	Bound on the bytes the site buffers hold at once (16 bytes per site, 24 for read counts). 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. Exceeding it is an error: decode fewer samples per call or raise it.
site_rule	Which record of a position listed in 'af_table' gives its site. '"matching"' (the default) uses a record whose REF and ALT match a row of 'af_table' ('duckhts_roh_af_table'). '"bcftools"' follows 'bcftools roh –AF-file' ('duckhts_roh_af_table_bcftools'): in file order, a record of the other indel class than the row is passed over when the next record of the file has the same POS (on any contig); a record with no non-symbolic ALT, or several, claims nothing; and the first remaining record claims the position, which is skipped unless its alleles equal the row's (also when its first ALT is symbolic, which bcftools would decode through the second ALT). It needs one 'af_table' row per position, unless 'af_order' gives file order, and a sorted, BGZF-compressed VCF or BCF, and reads the file a second time without samples. Only 'af_table' takes '"bcftools"'.
af_order	Optional name of a column of 'af_table' that gives file order, for example a line number added by the caller, or 'FILE_OFFSET' when the frequencies come from 'read_bcf()'. With 'site_rule = "bcftools"' only. Several rows may then share a position: the row with the least value is used, as 'bcftools roh –AF-file' reads the first line of a position. A 'NULL' value, a value repeated within a position, or a name that is not a column of 'af_table' is an error.

Details

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

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

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

Value

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

Examples

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

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