inst/doc/methcp.R

In MethCP: Differential methylation anlsysis for bisulfite sequencing data

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)

## ---- message=FALSE-----------------------------------------------------------
library(bsseq)
library(MethCP)

## ----readData-----------------------------------------------------------------
# The dataset is consist of 6 samples. 3 samples are H2A.Z mutant 
# plants, and 3 samples are controls.
sample_names <- c(
    paste0("control", seq_len(3)), 
    paste0("treatment", seq_len(3))
)

# Get the vector of file path and names 
raw_files <- system.file(
    "extdata", paste0(sample_names, ".txt"), package = "MethCP")

# load the data
bs_object <- createBsseqObject(
    files = raw_files, sample_names = sample_names, 
    chr_col = 'Chr', pos_col = 'Pos', m_col = "M", cov_col = 'Cov')

## ----showBSobject-------------------------------------------------------------
bs_object

## -----------------------------------------------------------------------------
dt <- read.table(
            raw_files[1], stringsAsFactors = FALSE, header = TRUE)
head(dt)

## ----calcStat-----------------------------------------------------------------
# the sample names of the two groups to compare. They should be subsets of the 
# sample names provided when creating the `bsseq` objects.
group1 <- paste0("control", seq_len(3))
group2 <- paste0("treatment", seq_len(3))

# Below we calculate the per-cytosine statistics using two different 
# test `DSS` and `methylKit`. The users may pick one of the two for their
#  application.
# obj_DSS <- calcLociStat(bs_object, group1, group2, test = "DSS")
obj_methylKit <- calcLociStat(
    bs_object, group1, group2, test = "methylKit")

## ----obj_methylKit------------------------------------------------------------
obj_methylKit

## ----createmethcp-------------------------------------------------------------
data <- data.frame(
    chr = rep("Chr01", 5),
    pos = c(2, 5, 9, 10, 18),
    effect.size = c(1,-1, NA, 9, Inf),
    pvals = c(0, 0.1, 0.9, NA, 0.02))
obj <- MethCPFromStat(
    data, test.name="myTest",
    pvals.field = "pvals",
    effect.size.field="effect.size",
    seqnames.field="chr",
    pos.field="pos"
)

## -----------------------------------------------------------------------------
obj

## ----segmentation-------------------------------------------------------------
# obj_DSS <- segmentMethCP(
#     obj_DSS, bs_object, region.test = "weighted-coverage")

obj_methylKit <- segmentMethCP(
    obj_methylKit, bs_object, region.test = "fisher")

## -----------------------------------------------------------------------------
obj_methylKit

## -----------------------------------------------------------------------------
# region_DSS <- getSigRegion(obj_DSS)
# head(region_DSS)

## -----------------------------------------------------------------------------
region_methylKit <- getSigRegion(obj_methylKit)
head(region_methylKit)

## -----------------------------------------------------------------------------
meta_file <- system.file(
    "extdata", "meta_data.txt", package = "MethCP")
meta <- read.table(meta_file, sep = "\t", header = TRUE)

head(meta)

## -----------------------------------------------------------------------------
# Get the vector of file path and names 
raw_files <- system.file(
    "extdata", paste0(meta$SampleName, ".tsv"), package = "MethCP")

# read files
bs_object <- createBsseqObject(
    files = raw_files, sample_names = meta$SampleName, 
    chr_col = 1, pos_col = 2, m_col = 4, cov_col = 5, header = TRUE)

## -----------------------------------------------------------------------------
groups <- split(seq_len(nrow(meta)), meta$Condition)
coverages <- as.data.frame(getCoverage(bs_object, type = "Cov"))
filter <- rowSums(coverages[, meta$SampleName[groups[[1]]]] != 0) >= 3 &
    rowSums(coverages[, meta$SampleName[groups[[2]]]] != 0) >= 3
bs_object <- bs_object[filter, ]

## -----------------------------------------------------------------------------
obj <- calcLociStatTimeCourse(bs_object, meta)

## -----------------------------------------------------------------------------
obj

## -----------------------------------------------------------------------------
obj <- segmentMethCP(obj, bs_object, region.test = "stouffer")

## -----------------------------------------------------------------------------
regions <- getSigRegion(obj)

## -----------------------------------------------------------------------------
head(regions)

## -----------------------------------------------------------------------------
sessionInfo()