inst/doc/updating-methylSig-code.R

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

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

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

## ----eval = FALSE-------------------------------------------------------------
#  meth = methylSigReadData(
#      fileList = files,
#      pData = pData,
#      assembly = 'hg19',
#      destranded = TRUE,
#      maxCount = 500,
#      minCount = 10,
#      filterSNPs = TRUE,
#      num.cores = 1,
#      fileType = 'cytosineReport')

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

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

## ----eval = FALSE-------------------------------------------------------------
#  # For genomic windows, tiles = NULL
#  windowed_meth = methylSigTile(meth, tiles = NULL, win.size = 10000)
#  
#  # For pre-defined tiles, tiles should be a GRanges object.

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

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

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

## ----eval = FALSE-------------------------------------------------------------
#  result = methylSigCalc(
#      meth = meth,
#      comparison = 'DR_vs_DS',
#      dispersion = 'both',
#      local.info = FALSE,
#      local.winsize = 200,
#      min.per.group = c(3,3),
#      weightFunc = methylSig_weightFunc,
#      T.approx = TRUE,
#      num.cores = 1)

## ----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_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)

## ----eval = FALSE-------------------------------------------------------------
#  contrast = matrix(c(0,1), ncol = 1)
#  result_dss = methylSigDSS(
#      meth = meth,
#      design = design1,
#      formula = '~ group',
#      contrast = contrast,
#      group.term = 'group',
#      min.per.group=c(3,3))

## ----filter_by_group_coverage2, eval = FALSE----------------------------------
#  # IF NOT DONE PREVIOUSLY
#  # 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_dss_fit_simple------------------------------------------------------
# Test the simplest model with an intercept and Type
diff_fit_simple = diff_dss_fit(
    bs = bs,
    design = bsseq::pData(bs),
    formula = as.formula('~ Type'))

## ----diff_dss_test_simple-----------------------------------------------------
# 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)

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

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