calcLociStat: Calculate the per-cytosine statistics for simple two-group...

In MethCP: Differential methylation anlsysis for bisulfite sequencing data

Description

calcLociStat calculates per-cytosine based statistics between two population groups.

Usage

calcLociStat(
    bs.object, group1, group2, test = c("DSS", "methylKit"),
    BPPARAM = bpparam())

Arguments

bs.object	a BSseq object from the bsseq package.
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().

Details

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.

Value

a MethCP object that is not segmented.

Examples

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

MethCP documentation built on Nov. 8, 2020, 8:24 p.m.