# methcp: User’s Guide

Boying Gong

#### 2020-04-12

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

Load packages.

    library(bsseq)
    library(MethCP)

# 2 Two-group comparison

In this section, we use the CpG methylation data from an Arabidopsis dataset available from GEO with accession number 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.

## 2.1 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.

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

    bs_object

    ## An object of type 'BSseq' with
    ##   82103 methylation loci
    ##   6 samples
    ## has not been smoothed
    ## All assays are in-memory

Header of the raw file for one of the samples.

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

    ##    Chr Pos M Cov
    ## 1 chr1 109 3   3
    ## 2 chr1 110 4   4
    ## 3 chr1 115 3   3
    ## 4 chr1 116 4   4
    ## 5 chr1 161 0   1
    ## 6 chr1 162 1   1

## 2.2 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.

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

    ## two groups detected:
    ##  will calculate methylation difference as the difference of
    ## treatment (group: 1) - control (group: 0)
    ## two groups detected:
    ##  will calculate methylation difference as the difference of
    ## treatment (group: 1) - control (group: 0)

    obj_methylKit

    ## MethCP object with 2 chromosomes, 81936 methylation loci
    ## test: methylKit
    ## group1: control1 control2 control3
    ## group2: treatment1 treatment2 treatment3
    ## has not been segmented

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.

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

    ## Filtering out NAs and infinite values.
    ## Cytosine counts before filtering: 5.
    ## Cytosine counts after filtering: 3.

    obj

    ## MethCP object with 1 chromosomes, 3 methylation loci
    ## test: myTest
    ## group1: notApplicable
    ## group2: notApplicable
    ## has not been segmented

## 2.3 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.

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

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

    obj_methylKit

    ## MethCP object with 2 chromosomes, 81936 methylation loci
    ## test: methylKit
    ## group1: control1 control2 control3
    ## group2: treatment1 treatment2 treatment3
    ## has been segmented

## 2.4 Significant regions

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

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

    region_methylKit <- getSigRegion(obj_methylKit)
    head(region_methylKit)

    ##      seqnames start   end nC.valid nC mean.diff mean.cov  region.pval
    ## 1.2      chr1 26853 27083       10 10   -0.3126   6.8500 6.306067e-14
    ## 1.5      chr1 30691 30768       11 11   -0.1776   8.3788 4.028997e-05
    ## 1.13     chr1 47817 47966       16 16   -0.2058   7.1354 4.314515e-08
    ## 1.21     chr1 65562 65735       10 10   -0.5742   5.2167 0.000000e+00
    ## 1.36     chr1 92672 92891       10 10    0.4067   8.9667 0.000000e+00
    ## 1.38     chr1 94712 95021       12 12    0.1304   7.9444 1.295822e-10

# 3 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 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.

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

    head(meta)

    ##   Condition TimeName Time Replicate                  SampleName
    ## 1     Col-0       DS    0         1   GSM2481288_mC_Col-0_DS_r1
    ## 2     Col-0       DS    0         2   GSM2481289_mC_Col-0_DS_r2
    ## 3     Col-0     0DAI    1         1 GSM2481290_mC_Col-0_0DAI_r1
    ## 4     Col-0     0DAI    1         2 GSM2481291_mC_Col-0_0DAI_r2
    ## 5     Col-0     1DAI    2         1 GSM2481292_mC_Col-0_1DAI_r1
    ## 6     Col-0     1DAI    2         2 GSM2481293_mC_Col-0_1DAI_r2

Read the counts data.

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

    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.

    obj <- calcLociStatTimeCourse(bs_object, meta)

    obj

    ## MethCP object with 1 chromosomes, 26001 methylation loci
    ## test: TimeCourse
    ## group1: GSM2481288_mC_Col-0_DS_r1 GSM2481289_mC_Col-0_DS_r2 GSM2481290_mC_Col-0_0DAI_r1 GSM2481291_mC_Col-0_0DAI_r2 GSM2481292_mC_Col-0_1DAI_r1 GSM2481293_mC_Col-0_1DAI_r2 GSM2481294_mC_Col-0_2DAI_r1 GSM2481295_mC_Col-0_2DAI_r2 GSM2481296_mC_Col-0_3DAI_r1 GSM2481297_mC_Col-0_3DAI_r2 GSM2481298_mC_Col-0_4DAI_r1 GSM2481299_mC_Col-0_4DAI_r2
    ## group2: GSM2481300_mC_rdd_DS_r1 GSM2481301_mC_rdd_DS_r2 GSM2481302_mC_rdd_0DAI_r1 GSM2481303_mC_rdd_0DAI_r2 GSM2481304_mC_rdd_1DAI_r1 GSM2481305_mC_rdd_1DAI_r2 GSM2481306_mC_rdd_2DAI_r1 GSM2481307_mC_rdd_2DAI_r2 GSM2481308_mC_rdd_3DAI_r1 GSM2481309_mC_rdd_3DAI_r2 GSM2481310_mC_rdd_4DAI_r1 GSM2481311_mC_rdd_4DAI_r2
    ## has not been segmented

Segmentation.

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

Get the DMRs.

    regions <- getSigRegion(obj)

    head(regions)

    ##     seqnames   start     end nC.valid  nC mean.diff mean.cov  region.pval
    ## 31         1  238367  239326       69  69   -0.1042   5.7144 0.0012861236
    ## 66         1  394175  395220       50  50   -0.2490   3.9292 0.0026789957
    ## 70         1  405983  408116       52  52   -0.1079   6.3197 0.0082183067
    ## 92         1  531037  532534       36  36   -0.1900   4.3507 0.0007397345
    ## 187        1  962256  966028      145 145   -0.1049   5.9529 0.0029463070
    ## 241        1 1258866 1259519       30  30   -0.1454   5.6153 0.0028433185

# 4 Session info

    sessionInfo()

    ## R version 4.0.0 alpha (2020-03-31 r78116)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.1252
    ## [3] LC_MONETARY=English_United States.1252
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.1252
    ##
    ## attached base packages:
    ## [1] stats4    parallel  stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] MethCP_1.1.6                bsseq_1.23.2
    ##  [3] SummarizedExperiment_1.17.5 DelayedArray_0.13.12
    ##  [5] matrixStats_0.56.0          Biobase_2.47.3
    ##  [7] GenomicRanges_1.39.3        GenomeInfoDb_1.23.16
    ##  [9] IRanges_2.21.8              S4Vectors_0.25.15
    ## [11] BiocGenerics_0.33.3         BiocStyle_2.15.6
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] splines_4.0.0            DelayedMatrixStats_1.9.1 R.utils_2.9.2
    ##  [4] gtools_3.8.2             assertthat_0.2.1         BiocManager_1.30.10
    ##  [7] BSgenome_1.55.4          GenomeInfoDbData_1.2.2   Rsamtools_2.3.7
    ## [10] yaml_2.2.1               pillar_1.4.3             numDeriv_2016.8-1.1
    ## [13] lattice_0.20-41          glue_1.4.0               limma_3.43.5
    ## [16] bbmle_1.0.23.1           digest_0.6.25            XVector_0.27.2
    ## [19] qvalue_2.19.0            colorspace_1.4-1         htmltools_0.4.0
    ## [22] Matrix_1.2-18            R.oo_1.23.0              plyr_1.8.6
    ## [25] pkgconfig_2.0.3          XML_3.99-0.3             emdbook_1.3.12
    ## [28] DSS_2.35.0               bookdown_0.18            zlibbioc_1.33.1
    ## [31] purrr_0.3.3              mvtnorm_1.1-0            snow_0.4-3
    ## [34] scales_1.1.0             HDF5Array_1.15.18        BiocParallel_1.21.2
    ## [37] tibble_3.0.0             ggplot2_3.3.0            ellipsis_0.3.0
    ## [40] cli_2.0.2                mclust_5.4.6             magrittr_1.5
    ## [43] crayon_1.3.4             evaluate_0.14            R.methodsS3_1.8.0
    ## [46] fansi_0.4.1              MASS_7.3-51.5            tools_4.0.0
    ## [49] data.table_1.12.8        lifecycle_0.2.0          stringr_1.4.0
    ## [52] Rhdf5lib_1.9.2           munsell_0.5.0            locfit_1.5-9.4
    ## [55] Biostrings_2.55.7        compiler_4.0.0           fastseg_1.33.0
    ## [58] rlang_0.4.5              rhdf5_2.31.10            grid_4.0.0
    ## [61] RCurl_1.98-1.1           methylKit_1.13.2         bitops_1.0-6
    ## [64] rmarkdown_2.1            DNAcopy_1.61.0           gtable_0.3.0
    ## [67] reshape2_1.4.4           R6_2.4.1                 GenomicAlignments_1.23.2
    ## [70] dplyr_0.8.5              knitr_1.28               rtracklayer_1.47.0
    ## [73] bdsmatrix_1.3-4          permute_0.9-5            stringi_1.4.6
    ## [76] Rcpp_1.0.4.6             vctrs_0.2.4              tidyselect_1.0.0
    ## [79] xfun_0.12                coda_0.19-3