Description Usage Arguments Details Value Examples
Generalized Linear Model for group comparison of information
divergence variables yielded by function
estimateDivergence
from MethylIT R package output.
Basically, this a wrapping function to perform the fitting of generalized
linear models with glm
from 'stats' package to any
variable of interest given in GRanges objects of
estimateDivergence
output.
1 2 3 4 5 6  divTest(GR, control.names, treatment.names, glm.family = Gamma(link =
"log"), var.weights = FALSE, weights = NULL, varFilter = 0,
meanFilter = 0, FilterLog2FC = TRUE, Minlog2FC = 1,
divPerBp = 0.001, minInd = 2, pAdjustMethod = NULL, scaling = 1L,
pvalCutOff = 0.05, saveAll = FALSE, num.cores = 1, tasks = 0L,
verbose = TRUE, ...)

GR 
GRanges objects including control and treatment samples containing an information divergence of methylation levels. The names for each column must coincide with the names given for parameters: 'control.names' and 'treatment.names'. 
control.names 
Names/IDs of the control samples, which must be included in the variable GR in a metacolumn. 
treatment.names 
Names/IDs of the treatment samples, which must be included in the variable GR in a metacolumn. 
glm.family, link 
Parameter to be passed to function

var.weights 
Logical (default: FALSE). Whether to use group variances as weights. 
weights 
An optional list of two numeric vectors of ‘prior weights’ to be used in the fitting process. One vector of weights for the control and one for the treatment. Each vector with length equal to length(GR) (default: NULL). NonNULL weights can be used to indicate that different observations have different dispersions (with the values in weights being inversely proportional to the dispersions). 
varFilter 
Numeric (default: 0). GLM will be performed only for those rows (ranges denoting genomic regions) where the group variance is greater the number specified by varFilter. 
meanFilter 
Numeric (default: 0). GLM will be performed only for those rows (ranges denoting genomic regions) where the absolute difference of group means is greater the number specified by meanFilter. 
FilterLog2FC 
if TRUE, the results are filtered using the minimun absolute value of log2FoldChanges observed to accept that a gene in the treatment is differentially expressed in respect to the control. 
Minlog2FC 
minimum logarithm base 2 of fold changes 
divPerBp 
At least for one group the mean divergence per bp must be equal to or greater than 'divPerBp' (default divPerBp = 0.001). 
minInd 
Integer (Default: 2). At least one group must have 'minInd' individuals with a divergence value greater than zero. 
pAdjustMethod 
Method used to adjust the results; default: "NULL"
(see 
scaling 
integer (default 1). Scaling factor estimate the signal density as: scaling * "DIMPCountPerBp". For example, if scaling = 1000, then signal density denotes the number of DIMPs in 1000 bp. 
pvalCutOff 
cutoff used then a pvalue adjustment is performed 
saveAll 
if TRUE all the temporal results that passed filters 'varFilter' and are 'meanFilter' returned. If FALSE, only the comparisons that passed filters 'varFilter', 'meanFilter', and pvalue < pvalCutOff or adj.pvalue < pvalCutOff (if pAdjustMethod is not NULL) are returned. 
num.cores 
The number of cores to use, i.e. at most how many child
processes will be run simultaneously (see

tasks 
integer(1). The number of tasks per job. Value must be a scalar
integer >= 0L. In this documentation a job is defined as a single call to
a function, such as bplapply, bpmapply etc. A task is the division of the
X argument into chunks. When tasks == 0 (default), X is divided as evenly
as possible over the number of workers (see

verbose 
if TRUE, prints the function log to stdout 
... 
Additional parameters passed to 
The default parameter setting glm.family = Gamma(link = "log") is
thought to perform the group comparison of the sums of absolute
differences of methylation levels (total variation distance (TVD) at
genebody DIMPs on DMGs). The sums of Hellinger divergence (HD, at
genebody DIMPs on DMGs) can be tested with this setting as well. Both
TVD and HD follow asymptotic Chisquare distribution and, consequently,
so do the sum of TVD and the sum of HD. The Chisquare distribution is
a particular case of Gamma distribution:
f(xa,s) = 1/(s^a Gamma(a)) x^(a1) e^(x/s)
Chisquare density is derived after replacing a = n/2 and s = 2:
f(xn) = 1/(2^(n/2) Gamma(n/2)) x^(n/21) e^(x/2)
The original GRanges object with the columns "beta", "log2FC", "pvalue", "adj.pval" (if pAdjustMethod requested), "CT.divPerBp" and "TT.divPerBp" (divergence per base pairs), and "divPerBpVariation added.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52  ## Gene annotation
genes < GRanges(seqnames = "1",
ranges = IRanges(start = c(3631, 6788, 11649),
end = c(5899, 9130, 13714)),
strand = c("+", "", ""))
mcols(genes) < data.frame(gene_id = c("AT1G01010", "AT1G01020",
"AT1G01030"))
# === The number of cytosine sites to generate ===
sites = 11001
# == Set a seed for pseudorandom number generation ===
set.seed(123)
alpha.ct < 0.09
alpha.tt < 0.2
# === Simulate samples ===
ref = simulateCounts(num.samples = 2, sites = sites, alpha = alpha.ct,
beta = 0.5, size = 50, theta = 4.5, sample.ids = "R1")
# Control group
ctrl = simulateCounts(num.samples = 2, sites = sites, alpha = alpha.ct,
beta = 0.5, size = 50, theta = 4.5,
sample.ids = c("C1", "C2"))
# Treatment group
treat = simulateCounts(num.samples = 2, sites = sites, alpha = alpha.tt,
beta = 0.5, size = 50, theta = 4.5,
sample.ids = c("T1", "T2"))
# === Estime Divergences ===
HD = estimateDivergence(ref = ref$R1, indiv = c(ctrl, treat),
Bayesian = TRUE, num.cores = 1L, percentile = 1,
verbose = FALSE)
nlms < nonlinearFitDist(HD, column = 4, verbose = FALSE)
## Next, the potential signal can be estimated
PS < getPotentialDIMP(LR = HD, nlms = nlms, div.col = 4, alpha = 0.05)
## The cutpoint estimation used to discriminate the signal from the noise
cutpoints < estimateCutPoint(PS, control.names = c("C1", "C2"),
treatment.names = c("T1", "T2"),
div.col = 4, verbose = FALSE)
## DIMPs are selected using the cupoints
DIMPs < selectDIMP(PS, div.col = 9, cutpoint = min(cutpoints$cutpoint))
## Finally DIMPs statistics genes
tv_DIMPs = getGRegionsStat(GR = DIMPs, grfeatures = genes, stat = "sum",
absolute = TRUE, column = 7L)
GR_tv_DIMP = uniqueGRanges(tv_DIMPs, type = "equal", chromosomes = "1")
colnames(mcols(GR_tv_DIMP)) < c("C1", "C2", "T1", "T2")
res < divTest(GR=GR_tv_DIMP, control.names = c("C1", "C2"),
treatment.names = c("T1", "T2"))

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.