inst/doc/using-methylSig.R

In methylSig: MethylSig: Differential Methylation Testing for WGBS and RRBS Data

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(methylSig)

## ----install, eval=FALSE------------------------------------------------------
#  devtools::install_github('sartorlab/methylSig')

## ----read---------------------------------------------------------------------
files = c(
    system.file('extdata', 'bis_cov1.cov', package='methylSig'),
    system.file('extdata', 'bis_cov2.cov', package='methylSig')
)

bsseq_stranded = bsseq::read.bismark(
    files = files,
    colData = data.frame(row.names = c('test1','test2')),
    rmZeroCov = FALSE,
    strandCollapse = FALSE
)

## ----bsseq_access-------------------------------------------------------------
# pData
bsseq::pData(bsseq_stranded)

# GRanges
GenomicRanges::granges(bsseq_stranded)

# Coverage matrix
bsseq::getCoverage(bsseq_stranded, type = 'Cov')

# Methylation matrix
bsseq::getCoverage(bsseq_stranded, type = 'M')

## ----filter_by_coverage-------------------------------------------------------
# Load data for use in the rest of the vignette
data(BS.cancer.ex, package = 'bsseqData')
bs = BS.cancer.ex[1:10000]

bs = filter_loci_by_coverage(bs, min_count = 5, max_count = 500)

## ----filter_by_location-------------------------------------------------------
# Show locations of bs
GenomicRanges::granges(bs)

# Construct GRanges object
remove_gr = GenomicRanges::GRanges(
    seqnames = c('chr21', 'chr21', 'chr21'),
    ranges = IRanges::IRanges(
        start = c(9411552, 9411784, 9412099),
        end = c(9411552, 9411784, 9412099)
    )
)

bs = filter_loci_by_location(bs = bs, gr = remove_gr)

# Show removal
GenomicRanges::granges(bs)

## ----tile_by_windows----------------------------------------------------------
windowed_bs = tile_by_windows(bs = bs, win_size = 10000)

# Show tiling
GenomicRanges::granges(windowed_bs)

## ----tile_by_regions----------------------------------------------------------
# Collapsed promoters on chr21 and chr22
data(promoters_gr, package = 'methylSig')

promoters_bs = tile_by_regions(bs = bs, gr = promoters_gr)

## ----filter_by_group_coverage-------------------------------------------------
# Look a the phenotype data for bs
bsseq::pData(bs)

# Require at least two samples from cancer and two samples from normal
bs = filter_loci_by_group_coverage(
    bs = bs,
    group_column = 'Type',
    c('cancer' = 2, 'normal' = 2))

## ----diff_binomial------------------------------------------------------------
# Test cancer versus normal
diff_gr = diff_binomial(
    bs = bs,
    group_column = 'Type',
    comparison_groups = c('case' = 'cancer', 'control' = 'normal'))

diff_gr

## ----diff_methylsig-----------------------------------------------------------
# Test cancer versus normal with dispersion from both groups
diff_gr = diff_methylsig(
    bs = bs,
    group_column = 'Type',
    comparison_groups = c('case' = 'cancer', 'control' = 'normal'),
    disp_groups = c('case' = TRUE, 'control' = TRUE),
    local_window_size = 0,
    t_approx = TRUE,
    n_cores = 1)

diff_gr

## ----add_numerical_covariate--------------------------------------------------
bsseq::pData(bs)$num_covariate = c(84, 96, 93, 10, 18, 9)

## ----diff_dss_fit_simple------------------------------------------------------
diff_fit_simple = diff_dss_fit(
    bs = bs,
    design = bsseq::pData(bs),
    formula = as.formula('~ Type'))

## ----diff_dss_fit_paired------------------------------------------------------
# Paired-test
diff_fit_paired = diff_dss_fit(
    bs = bs,
    design = bsseq::pData(bs),
    formula = '~ Type + Pair')

## ----diff_dss_fit_num---------------------------------------------------------
# Numerical covariate test
diff_fit_num = diff_dss_fit(
    bs = bs,
    design = bsseq::pData(bs),
    formula = '~ num_covariate')

## ----diff_dss_fit_model-------------------------------------------------------
diff_fit_simple$X

diff_fit_paired$X

diff_fit_num$X

## ----contrast-----------------------------------------------------------------
# Test the simplest model for cancer vs normal
# Note, 2 rows corresponds to 2 columns in diff_fit_simple$X
simple_contrast = matrix(c(0,1), ncol = 1)

# Test the paired model for cancer vs normal
# Note, 4 rows corresponds to 4 columns in diff_fit_paired$X
paired_contrast = matrix(c(0,1,0,0), ncol = 1)

# Test the numerical covariate
num_contrast = matrix(c(0,1), ncol = 1)

## ----diff_dss_test_simple-----------------------------------------------------
diff_simple_gr = diff_dss_test(
    bs = bs,
    diff_fit = diff_fit_simple,
    contrast = simple_contrast,
    methylation_group_column = 'Type',
    methylation_groups = c('case' = 'cancer', 'control' = 'normal'))

diff_simple_gr

## ----diff_dss_test_paired-----------------------------------------------------
diff_paired_gr = diff_dss_test(
    bs = bs,
    diff_fit = diff_fit_paired,
    contrast = paired_contrast,
    methylation_group_column = 'Type',
    methylation_groups = c('case' = 'cancer', 'control' = 'normal'))

diff_paired_gr

## ----diff_dss_test_num--------------------------------------------------------
diff_num_gr = diff_dss_test(
    bs = bs,
    diff_fit = diff_fit_num,
    contrast = num_contrast,
    methylation_group_column = 'num_covariate')

diff_num_gr

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