---
title: "methcp: User’s Guide"
author: "Boying Gong"
date: "`r Sys.Date()`"
vignette: >
    %\VignetteIndexEntry{methcp: User’s Guide}
    %\VignetteEngine{knitr::rmarkdown}
    \usepackage[utf8]{inputenc}
output: 
    BiocStyle::html_document:
        toc_float: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = FALSE)
```

# Introduction

`methcp` is a differentially methylated region 
(DMR) detecting method for whole-genome bisulfite sequencing (WGBS) 
data. It is applicable for a wide range of experimental designs.
In this document, we provide examples for two-group comparisons and 
time-course analysis.  

`methcp` identifies DMRs based on change point detection, which 
naturally segments the genome and provides region-level 
differential analysis. We direct the interested reader to our paper
[here](https://link.springer.com/chapter/10.1007/978-3-030-17083-7_5)

Load packages. 

```{r, message=FALSE}
library(bsseq)
library(MethCP)
```

# Two-group comparison

In this section, we use the CpG methylation data from an Arabidopsis dataset 
available from GEO with accession number
[GSM954584](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM954584). We 
take a subset of chromosome 1 and 2 from each of the samples and perform 
differential analysis between the wild-type plants and the H2A.Z mutant 
plants, which we refer to as `treatment` and `control` in the rest of 
the document. 

## Read data

We use a well-developed Bioconductor package `bsseq` to load the store the raw 
data. Below is an example of how to read raw counts using `bsseq`. 

We provide a helper function `createBsseqObject` to create a bsseq object
when the data for each sample is stored in a separate text file. 

For more operations regarding the `bsseq` object, or to create a `bsseq` 
object customized to your file format, please refer to their 
[User's Guide](http://bioconductor.org/packages/release/bioc/html/bsseq.html). 

```{r 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')
```

The `bsseq` object.

```{r showBSobject}
bs_object
```

Header of the raw file for one of the samples.

```{r}
dt <- read.table(
            raw_files[1], stringsAsFactors = FALSE, header = TRUE)
head(dt)
```

## Calculate the statistics

We calculate the per-cytosine statistics using two different test `DSS` and 
`methylKit` using the function `calcLociStat`. This function returns a 
`MethCP` object that is used for the future segmentation step. We allow 
parallelized computing when there are multiple chromosomes in the dataset. 
Note that the number of cores used as indicated in the parameter `mc.cores` 
should not exceeds the number of chromosomes in the dataset.

```{r 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", mc.cores = 1)
obj_methylKit <- calcLociStat(
    bs_object, group1, group2, test = "methylKit", mc.cores = 1)
```

```{r obj_DSS}
obj_DSS
```

```{r obj_methylKit}
obj_methylKit
```

In cases the user wants to use their pre-calculated test statistics for 
experiments other than two-group comparison and time course data, we use the 
calculated statistics and create a `MethCP` object.

```{r 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"
)
```

```{r}
obj
```

## Segmentation

`segmentMethCP` performs segmentation on a `MethCP` object. We allow 
parallelized computing when there are multiple chromosomes in the dataset. 
Different from `calcLociStat` function in the previous section, we do not put 
any constraint on the number of cores used. Please see the documentation for 
adjusting the parameters used in the segmentation.

```{r segmentation}
obj_DSS <- segmentMethCP(
    obj_DSS, bs_object, region.test = "weighted-coverage", 
    mc.cores = 1)

obj_methylKit <- segmentMethCP(
    obj_methylKit, bs_object, region.test = "fisher", mc.cores = 1)
```

```{r}
obj_DSS
```

```{r}
obj_methylKit
```

## Significant regions

Use function `getSigRegion` on a `MethCP` object to get the list of DMRs.

```{r}
region_DSS <- getSigRegion(obj_DSS)
head(region_DSS)
```

```{r}
region_methylKit <- getSigRegion(obj_methylKit)
head(region_methylKit)
```

# A time-course example

MethCP is flexible for a wide variety of experimental designs. We apply MethCP 
on an Arabidopsis thaliana seed germination dataset available from the GEO 
with accession number 
[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE94712](GSE94712) 
The data is generated from two replicates of dry seed and germinating seeds of 
wild-type plants (Col-0) and ros1 dml2 dml3 (rdd) triple demethylase mutant 
plants at 0-4 days after imbibition for 4 days (DAI). For cytosine-based 
statistics, we fit linear models on the methylation ratios and test the 
differences between the time coefficient of condition Col-0 and condition rdd.

Read the meta data.

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

head(meta)
```

Read the counts data.

```{r}
# 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)
```

Apply coverage filter to make sure each loci has total coverage (summed across 
samples) more than 3 for each condition.

```{r}
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, ]
```

Calculate the statistics. A dataframe of the meta data will be passed to 
function `calcLociStatTimeCourse`. Note that there must be columns named 
`Condition`, `Time` and `SampleName` in the dataframe.

```{r}
obj <- calcLociStatTimeCourse(bs_object, meta)
```

```{r}
obj
```

Segmentation.

```{r}
obj <- segmentMethCP(obj, bs_object, region.test = "stouffer", mc.cores = 1)
```

Get the DMRs.

```{r}
regions <- getSigRegion(obj)
```

```{r}
head(regions)
```


# Session info

```{r}
sessionInfo() 
```
