calcLociStat {MethCP} | R Documentation |
calcLociStat
calculates per-cytosine based statistics
between two population groups.
calcLociStat( bs.object, group1, group2, test = c("DSS", "methylKit"), BPPARAM = bpparam())
bs.object |
a |
group1 |
a character vector containing the sample names of the treatment group. |
group2 |
a character vector containing the sample names of the control group. |
test |
a character string containing the names of the test to be performed per cytosine. |
BPPARAM |
An optional BiocParallelParam instance determining the parallel back-end to be used during evaluation, or a list of BiocParallelParam instances, to be applied in sequence for nested calls to BiocParallel functions. Default bpparam(). |
For each cytosine, calcLociStat
calculates a statistics
using either package DSS
or methylKit
to test the
differences between two groups, and returns a MethCP
object.
For customized per-cytosine statistics, please use the function
methcpFromStat
.
The input bs.object
is a BSseq
object from the bsseq
package which contains the raw data including coverges, methylated
counts and position infomation for every cytosine in the dataset.
a MethCP
object that is not segmented.
library(bsseq) library(GenomicRanges) library(IRanges) set.seed(0286374) # Similate a small dataset with 11 cyotsine and 6 samples, # 3 in the treatment group and 3 in the control group. The # methylation ratio are generated using Binomial distribution # with probability 0.3. nC <- 2000 sim_cov <- rnbinom(6*nC, 5, 0.5) + 5 sim_M <- vapply( sim_cov, function(x) rbinom(1, x, 0.3), FUN.VALUE = numeric(1)) sim_cov <- matrix(sim_cov, ncol = 6) sim_M <- matrix(sim_M, ncol = 6) # methylation ratios in the DMRs in the treatment group are # generated using Binomial(0.7) DMRs <- c(600:622, 1089:1103, 1698:1750) sim_M[DMRs, 1:3] <- vapply( sim_cov[DMRs, 1:3], function(x) rbinom(1, x, 0.7), FUN.VALUE = numeric(1)) # sample names sample_names <- c(paste0("treatment", 1:3), paste0("control", 1:3)) colnames(sim_cov) <- sample_names colnames(sim_M) <- sample_names # create a bs.object bs_object <- BSseq(gr = GRanges( seqnames = "Chr01", IRanges( start = (1:nC)*10, width = 1)), Cov = sim_cov, M = sim_M, sampleNames = sample_names) # methcp_obj1 <- calcLociStat( # bs_object, # group1 = paste0("treatment", 1:3), # group2 = paste0("control", 1:3), # test = "DSS") methcp_obj2 <- calcLociStat( bs_object, group1 = paste0("treatment", 1:3), group2 = paste0("control", 1:3), test = "methylKit")