calcLociStatTimeCourse: Calculate the per-cytosine statistics for time-course data.

In boyinggong/methcp: Differential methylation anlsysis for bisulfite sequencing data

Description

For each cytosine, calcLociStatTimeCourse fits a linear model on the arcsin-tranformed methylation ratios, and test the differences of the slope between the treatment and the control group.

Usage

calcLociStatTimeCourse(
    bs.object, meta, force.slope = FALSE,
    BPPARAM = bpparam())

Arguments

bs.object	a BSseq object from the bsseq package.
meta	a data.frame. See details.
force.slope	if TRUE, we force the slope in the linear model to be the same between two conditions. Otherwise, the slopes are fitted separately but not tested.
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

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. meta must contain columns named Condition, Time and SampleName in the dataframe. They are used to fit the linear model.

Value

A MethCP object that is not segmented.

Examples

library(bsseq)
# Simulate a small dataset with 2000 cyotsine and 10 samples,
# 5 in the treatment group and 5 in the control group. The
# methylation ratio are generated using Binomial distribution
# with probability 0.3, 0.4, 0.5, 0.6 and 0.7 for 5 time points.
nC <- 2000
nsamples <- 5
sim_cov <- rnbinom(10*nC, 5, 0.5) + 5
sim_cov <- matrix(sim_cov, ncol = 10)
time_point <- rep(1:nsamples, 2)
ratios <- time_point/10 + 0.2
sim_M <- sapply(1:(2*nsamples), function(i){
    sapply(sim_cov[, i], function(j) rbinom(1, j, ratios[i]))
})
sim_M <- matrix(sim_M, ncol = 2*nsamples)
# methylation ratios in the DMRs in the treatment group are
# generated using Binomial(0.3)
DMRs <- c(600:622, 1089:1103, 1698:1750)
sim_M[DMRs, 1:5] <- sapply(
    sim_cov[DMRs, 1:5], function(x) rbinom(1, x, 0.3))
# sample names
sample_names <- c(paste0("treatment", 1:nsamples),
paste0("control", 1:nsamples))
colnames(sim_cov) <- sample_names
colnames(sim_M) <- sample_names

# create a bs.object
bs_object_ts <- BSseq(gr = GRanges(
    seqnames = "Chr01", IRanges(
        start = (1:nC)*10, width = 1)),
    Cov = sim_cov, M = sim_M, sampleNames = sample_names)
DMRs_pos_ts <- DMRs*10
meta <- data.frame(
    Condition = rep(
        c("treatment", "control"),
        each = nsamples),
    SampleName = sample_names,
    Time = time_point)
obj_ts <- calcLociStatTimeCourse(bs_object_ts, meta)
obj_ts

boyinggong/methcp documentation built on April 25, 2021, 9 a.m.