Obtains peak locations and heights based on adjusted readcounts

Fits additive models of peaks of varying strengths in various locations to the adjusted readcounts via RJMCMC, and stores these models on disk.

peaky(
  bait,
  omega_power,
  iterations = 1e+06,
  min_interactions = 20,
  log_file = NA
)

Arguments

bait	Data table containing the putative interactions of a bait, having the columns 'baitID', 'dist', and 'residual'. These report the bait ID, its distance to putative preys, and the adjusted readcounts for its interactions with them, respectively.
omega_power	Expected decay of adjusted read counts around a truly interacting prey. See details.
iterations	Number of models to parametrize. Greated numbers should lead to increased reproducibility.
log_file	Path to a log file.

Value

The output directory.

Details

The steepness of the function to be fitted to putative peaks is determined by \(\omega\) according to \(\beta\) exp(-\(|\omega * d|\)), where \(\beta\) represents peak height and \(d\) the distance from the center of the peak in bp.

Examples

base = system.file("extdata",package="peaky")
interactions_file = paste0(base,"/counts.tsv")
fragments_file = paste0(base,"/fragments.bed")
interactions = data.table::fread(interactions_file)
fragments = data.table::fread(fragments_file)
if (FALSE) {
BI = bin_interactions(interactions, fragments, bins=5)
models = by(BI$interactions, BI$interactions$dist.bin, model_bin, subsample_size=1000)
residuals = lapply(models, "[[", "residuals")
bins = split(BI$interactions, BI$interactions$dist.bin)

BTS = split_baits(bins, residuals)

peaky(BTS[baitID==618421], omega_power=-4.7, iterations=10e3)
}