rduckhts_roh_counts: Find Runs of Homozygosity from Read Counts

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

Find Runs of Homozygosity from Read Counts

Description

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

Usage

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

Arguments

con	A DuckDB connection with DuckHTS loaded.
counts_table	Name of a table or view of read counts (see Details).
seq_error	Per-read probability of showing the other allele, in '(0, 0.5)'.
contamination	Fraction of reads from a contaminating individual, in '[0, 1)'.
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]'.
rec_rate	Optional constant recombination rate per base pair ('-M').
table_name	Optional output table. 'NULL' returns a data frame.
overwrite	Whether an existing output table may be replaced.
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.
fit_table	Optional table or view with one row per 'sample', containing non-NULL 'seq_error', 'contamination', 'allele_balance' and 'artefact_weight', and an 'artefact_weight_sd' column, such as the output of [rduckhts_count_error_fit()]. These replace the scalar arguments for each matching 'sample_id'. A sample of 'counts_table' with no fit row, several fit rows, a NULL value or a value out of range is an error; fit rows of other samples are not checked. A NULL 'artefact_weight_sd' counts as 0. Filter fit status in the caller's SQL; the spreads are not used.
allele_balance	Share of a heterozygote's reads that show the counted allele, in '[0.35, 0.65]'. The default 0.5 assumes no mapping bias.
artefact_weight	Chance that a site's counts are uniform on '0..depth' whatever its genotype, in '[0, 0.05]'. The default 0 gives such sites no term. Values outside either range are errors, not clipped.

Details

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

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

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

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

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

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

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

Value

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

Examples

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

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