interpret_peaky.RdExtracts the marginal posterior probabilities of contact and other relevant statistics from the RJMCMC results and merges these with bait-associated information.
interpret_peaky(bait, peaks, omega_power, log_file = NA)
| 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. |
|---|---|
| peaks | The models fitted by peaky. |
| omega_power | The same value as used when running peaky, i.e. the expected decay of adjusted read counts around a truly interacting prey. See details. |
A data table containing bait-associated information, posterior probabilities of contact and other statistics.
The steepness of the function to be fitted to putative peaks is determined by \(\omega\) according to \(\beta \exp{- \abs{\omega * d}}\), where \(\beta\) represents peak height and \(d\) the distance from the center of the peak in bp.
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) relevant_bait = BTS[baitID==618421] omega_power=4.7 PKS = peaky(relevant_bait, omega_power, iterations=1e6) interpret_peaky(relevant_bait, PKS, omega_power) }