knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
We will analyze a subset of one of the sample dataset for illustration purposes.
library(parallel) library(npbin) library(data.table) minimum_coverage <- 5 # minimum total coverage allowed n_cores <- detectCores() # the number of cores to be used, can ONLY be 1 if run on Windows. dt <- ctcf colnames(dt)
dt.ct <- data.table(dt)[m >= minimum_coverage, ]
for illustration purpose, keep only the 2000 points, remove this line will end up with analyzing the whole data set. It could be slow if only one core is used.
dt.ct <- dt.ct[1:2000, ] dt.ct[, p_hat:=xm / m] n <- nrow(dt.ct)
NPBin
n_breaks <- 11 # number of breaks spline_order <- 4 # order of splines breaks <- seq(0, 1, length.out = n_breaks) pi_init <- initialize_weights( dt.ct, n_breaks = n_breaks, spline_order = spline_order, plot = TRUE ) # initialized the weights using the histogram of p_hat
estimate the overall model
overall_model_estimate <- emBinBspl( dt.ct[, xm], dt.ct[, m], breaks = breaks, k = spline_order, pi.init = pi_init, ncores = n_cores, err.max = 1e-3, iter.max = 200 ) names(overall_model_estimate)
estimate the null model
null_model_estimate <- estNull( dt.ct[, xm], dt.ct[, m], overall_model_estimate, init = NULL, iter.max = 200, ncores = n_cores, ub = rep(log(1e4), 2), err.max = 1e-4 ) names(null_model_estimate)
null_model_estimate[["coef.null"]]
Calculate the overdispersion/correlation parameter
1 / sum( null_model_estimate[["coef.null"]][["shape1"]], null_model_estimate[["coef.null"]][["shape2"]], 1 )
Plot the estimated null distribution
plot_estimated_null( dt.ct[, p_hat], shape1_shape2 = c( null_model_estimate[["coef.null"]][["shape1"]], null_model_estimate[["coef.null"]][["shape2"]] ) )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.