rduckhts_count_error_fit: Fit Read Error and Contamination from Allele Read Counts

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

Fit Read Error and Contamination from Allele Read Counts

Description

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

Usage

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

Arguments

con	A DuckDB connection with DuckHTS loaded.
counts_table	Name of a table or view with 'sample_id', 'chrom', 'pos', 'ref_count', 'alt_count' and 'af', as for [rduckhts_roh_counts()].
block_bases	Bases per block for the block refits, at least 1.
freq_bins	Frequency bins of the histogram, 1 to 1024.
max_depth	Sites with more reads are left out, 1 to 65535.
min_block_sites	Sites a block needs to be refitted, at least 1.
table_name	Optional output table. 'NULL' returns a data frame.
overwrite	Whether an existing output table may be replaced.
max_cells	Most histogram cells of one sample; an error above it.
max_cell_bytes	Most bytes the histogram buffers of all fits in the process hold at once; an error above it.

Details

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

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

Value

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

Examples

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

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