In anthony-aylward/npbin:

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 <- dnase
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"]]
  )
)

anthony-aylward/npbin documentation built on Aug. 22, 2019, 8:08 a.m.