Compute methylation levels

This function computes the methylation levels. This is an auxiliary function that can be applied to explore the data.

meth_levels(
  GR,
  x,
  columns = c(mC1 = 1, uC1 = 2, mC2 = NULL, uC2 = NULL),
  Bayesian = FALSE,
  init.pars = NULL,
  via.optim = TRUE,
  min.coverage = 4,
  tv = FALSE,
  bay.tv = FALSE,
  filter = FALSE,
  preserve.dt = FALSE,
  loss.fun = c("linear", "huber", "smooth", "cauchy", "arctg"),
  num.cores = 1,
  tasks = 0L,
  verbose = TRUE,
  ...
)

# S4 method for ANY,data.frame
meth_levels(
  GR,
  x,
  columns = c(mC1 = 1, uC1 = 2, mC2 = NULL, uC2 = NULL),
  Bayesian = FALSE,
  init.pars = NULL,
  via.optim = TRUE,
  min.coverage = 4,
  tv = FALSE,
  bay.tv = FALSE,
  filter = FALSE,
  preserve.dt = FALSE,
  loss.fun = c("linear", "huber", "smooth", "cauchy", "arctg"),
  num.cores = 1,
  tasks = 0L,
  verbose = TRUE,
  ...
)

# S4 method for GRanges,ANY
meth_levels(
  GR,
  x,
  columns = c(mC1 = 1, uC1 = 2, mC2 = NULL, uC2 = NULL),
  Bayesian = FALSE,
  init.pars = NULL,
  via.optim = TRUE,
  min.coverage = 4,
  tv = FALSE,
  bay.tv = FALSE,
  filter = FALSE,
  preserve.dt = FALSE,
  loss.fun = c("linear", "huber", "smooth", "cauchy", "arctg"),
  num.cores = 1,
  tasks = 0L,
  verbose = TRUE,
  ...
)

# S4 method for list,ANY
meth_levels(
  GR,
  x = NULL,
  columns = c(mC1 = 1, uC1 = 2, mC2 = NULL, uC2 = NULL),
  Bayesian = FALSE,
  init.pars = NULL,
  via.optim = TRUE,
  min.coverage = 4,
  tv = FALSE,
  bay.tv = FALSE,
  filter = FALSE,
  preserve.dt = FALSE,
  loss.fun = c("linear", "huber", "smooth", "cauchy", "arctg"),
  num.cores = detectCores() - 1,
  tasks = 0L,
  verbose = TRUE,
  ...
)

Arguments

GR, x

A GRanges-class object (GR) or 'data.frame' (x) with a matrix of counts in the meta-columns (methylated mC and unmethylated uC cytosines) or a list of GRanges-class objects.

columns

Vector of one or two integer numbers denoting the indexes of the columns where the methylated and unmethylated read counts are found. Unless specified in the parameter 'columns', the methylation counts must be given in the first four columns: 'mC1' and 'uC1' methylated and unmethylated counts for control sample, and 'mC2' and 'uC2' methylated and unmethylated counts for treatment sample, respectively.

Bayesian

logical(1). Whether to perform the estimations based on posterior estimations of methylation levels.

init.pars

initial parameter values. Defaults is NULL and an initial guess is estimated using optim function. If the initial guessing fails initial parameter values are to alpha = 1 & beta = 1, which imply the parsimony pseudo-counts greater than zero.

via.optim

Optional. Only used if Bayesian = TRUE Whether to estimate beta distribution parameters via optim or nls.lm. If any of this approaches fail then parameters used init.pars will be returned.

min.coverage

An integer or an integer vector of length 2. Cytosine sites where the coverage in both samples, 'x' and 'y', are less than min.coverage' are discarded. The cytosine site is preserved, however, if the coverage is greater than 'min.coverage' in at least one sample. If 'min.coverage' is an integer vector, then the corresponding min coverage is applied to each sample.

tv

logical(1). Whether to compute the total variation distance at each cytosine site. That is, the difference of methylation levels.

bay.tv

logical(1). Whether to compute the total variation distance at each cytosine site based on Bayesian estimation of methylation levels.

filter

logical(1). Optional. If TRUE, then only cytosine sites with \(coverages > min.coverage\) are including in the computation.

preserve.dt

logical(1). Option of whether to preserve all the metadata from the original 'data.frame' or GRanges-class object.

loss.fun

Described in estimateBetaDist.

num.cores, tasks

Parameters for parallel computation using package BiocParallel-package: the number of cores to use, i.e. at most how many child processes will be run simultaneously (see bplapply and the number of tasks per job (only for Linux OS).

verbose

if TRUE, prints the function log to stdout

...

Optional parameter values for: maxiter, ftol, ptol, and gradtol from nlsLM and nlm functions.

Author

Robersy Sanchez (https://genomaths.com)

Examples

## The read count data are created
num.samples <- 250
s <- 1:num.samples
gr <- data.frame(chr = 'chr1', start = s, end = s,
                strand = sample(c("+", "-"), num.samples, replace = TRUE),
                mCc = rnbinom(size = num.samples, mu = 4, n = 500),
                uCc = rnbinom(size = num.samples, mu = 4, n = 500),
                mCt = rnbinom(size = num.samples, mu = 4, n = 500),
                uCt = rnbinom(size = num.samples, mu = 4, n = 500))

gr <- makeGRangesFromDataFrame(gr, keep.extra.columns = TRUE)

gr <- meth_levels(GR = gr,
                columns = c(mC1 = 1, uC1 = 2,
                            mC2 = 3, uC2 = 4),
                preserve.dt = TRUE,
                Bayesian = TRUE, tv = TRUE, bay.tv = TRUE,
                num.cores = 1)
#> 
#> *** Estimating betaBinomial-posteriors... 
#> 
#> *** Estimating betaBinomial-posteriors...