DNA methylation studies have increased in number over the past decade thanks to the recent advances in next-generation sequencing (NGS) and microarray technology (MA), providing many data sets at high resolution, enabling researchers to understand methylation patterns and their regulatory roles in biological processes and diseases.
Notwithstanding that diverse methods and software have created ample opportunities for researchers to do quantitative analysis, they make it difficult for practitioners to choose the one that is suitable and efficient in analyzing DNA methylation data. Having examined most of differentially methylation identification tools for bisulfite sequencing (BS-Seq) data, we observed several drawbacks in the existing analytic tools. To address these issues we have developed a novel differentially methylated CpG site identification tool which is based on Hidden Markov models (HMM) called DMCHMM. This vignette provides some guidelines on how to use the package and analyze BS-Seq data.
Following topics will be discussed in this vignette:
Two different classes are defined by extending the SummarizedExperiment-class. The BSData-class is designed to hold BS-Seq data. Similarly, cBSData-method is defined to create a BSData object. This class includes two slots: the methReads, a matrix with columns representing samples and rows representing genomic positions (CpG sites) and elements of matrix representing methylation counts at each position in each sample; the totalReads, a matrix with similar columns and rows except the elements representing total number of reads.
For reading raw BS-Seq data we adopted The readBismark function from BiSeq package. The readBismark-method reads samples stored in different files with six columns of chromosome, start position, end position, methylation percentage, number of Cs and number of Ts.
Three data files are included in the DMCHMM package for illustration. The data can be imported using following code.
library(DMCHMM)
fn <- list.files(system.file("extdata",package = "DMCHMM"))
fn.f <- list.files(system.file("extdata",package="DMCHMM"), full.names=TRUE)
OBJ <- readBismark(fn.f, fn, mc.cores = 2)## Building BSData object.cdOBJ <- DataFrame(Cell = factor(c("BC", "TC","Mono"),
labels = c("BC", "TC", "Mono")), row.names = c("BCU1568","BCU173","BCU551"))
colData(OBJ) <- cdOBJ
OBJ## class: BSData
## dim: 25668 3
## metadata(0):
## assays(2): totalReads methReads
## rownames(25668): 1 2 ... 25667 25668
## rowData names(0):
## colnames(3): BCU1568 BCU173 BCU551
## colData names(1): CellThe above data set only include one sample for each cell type. We need more samples to be able to compare their methylations and find DMCs. For illustration we generate a sample of BS-Seq data as follows.
nr <- 150; nc <- 8
metht <- matrix(as.integer(runif(nr * nc, 0, 20)), nr)
methc <- matrix(rbinom(n=nr*nc,c(metht),prob = runif(nr*nc)),nr,nc)
r1 <- GRanges(rep("chr1", nr), IRanges(1:nr, width=1), strand="*")
names(r1) <- 1:nr
cd1 <- DataFrame(Group=rep(c("G1","G2"),each=nc/2),row.names=LETTERS[1:nc])
OBJ1 <- cBSData(rowRanges=r1,methReads=methc,totalReads=metht,colData=cd1)
OBJ1## class: BSData
## dim: 150 8
## metadata(0):
## assays(2): totalReads methReads
## rownames(150): 1 2 ... 149 150
## rowData names(0):
## colnames(8): A B ... G H
## colData names(1): GroupThere are two approaches to smoothed the data before testing for DMCs. Either EM or MCMC can be used to predict methylation levels utilizing HMM. The methHMEM-method which is developed to predict methylation levels. The output is a BSDMCs-class that can be either used to find DMCs or use MCMC algorithm to re-smooth the raw data. The process is as follows.
## class: BSDMCs
## dim: 150 8
## metadata(3): K Beta Pm
## assays(5): methReads totalReads methLevels methStates methVars
## rownames(150): 1 2 ... 149 150
## rowData names(0):
## colnames(8): A B ... G H
## colData names(1): GroupAlthough EM algorithm is a fast way to smooth the data but the results are not as good as the MCMC algorithm. The MCMC algorithm, however, is slow. In order to increase the speed, we first use methHMEM-method to find the HMM order for each sample and then we use methHMCMC-method to predict methylation levels. The procedure is as follows.
## class: BSDMCs
## dim: 150 8
## metadata(3): K Beta Pm
## assays(5): methReads totalReads methLevels methStates methVars
## rownames(150): 1 2 ... 149 150
## rowData names(0):
## colnames(8): A B ... G H
## colData names(1): GroupHaving smoothed the data using HMM, we run linear between predicted methylation levels and grouping covariate. In case other covariates exist, one can use the formula argument to specify a linear model. When there is no covariates no action is required. The following command identifies the DMCs. The results are stored in a BSDMCs-class and can be retrieved by calling metadata command.
## Warning in fdrtool(x, statistic = "pvalue", plot = FALSE, verbose = FALSE):
## There may be too few input test statistics for reliable FDR calculations!## DMCs pvalues qvalues DMCsGroupG1vsG2 methDirGroupG1vsG2
## 1 0 0.75597968 0.7017300 0 hypo
## 2 0 0.09361918 0.5513001 0 hypo
## 3 0 0.40777707 0.6108087 0 hyper
## 4 0 0.06379220 0.5513001 0 hypo
## 5 0 0.10955666 0.5513001 0 hyper
## 6 0 0.03267532 0.5513001 0 hypo## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] DMCHMM_1.16.0 fdrtool_1.2.16
## [3] BiocParallel_1.28.0 SummarizedExperiment_1.24.0
## [5] Biobase_2.54.0 GenomicRanges_1.46.0
## [7] GenomeInfoDb_1.30.0 IRanges_2.28.0
## [9] S4Vectors_0.32.0 BiocGenerics_0.40.0
## [11] MatrixGenerics_1.6.0 matrixStats_0.61.0
##
## loaded via a namespace (and not attached):
## [1] zoo_1.8-9 xfun_0.27 bslib_0.3.1
## [4] splines_4.1.1 lattice_0.20-45 htmltools_0.5.2
## [7] rtracklayer_1.54.0 yaml_2.2.1 survival_3.2-13
## [10] XML_3.99-0.8 rlang_0.4.12 jquerylib_0.1.4
## [13] calibrate_1.7.7 multcomp_1.4-17 GenomeInfoDbData_1.2.7
## [16] stringr_1.4.0 zlibbioc_1.40.0 Biostrings_2.62.0
## [19] mvtnorm_1.1-3 codetools_0.2-18 evaluate_0.14
## [22] restfulr_0.0.13 knitr_1.36 fastmap_1.1.0
## [25] parallel_4.1.1 TH.data_1.1-0 DelayedArray_0.20.0
## [28] jsonlite_1.7.2 XVector_0.34.0 Rsamtools_2.10.0
## [31] rjson_0.2.20 digest_0.6.28 stringi_1.7.5
## [34] BiocIO_1.4.0 grid_4.1.1 tools_4.1.1
## [37] bitops_1.0-7 sandwich_3.0-1 magrittr_2.0.1
## [40] sass_0.4.0 RCurl_1.98-1.5 crayon_1.4.1
## [43] MASS_7.3-54 Matrix_1.3-4 rmarkdown_2.11
## [46] R6_2.5.1 GenomicAlignments_1.30.0 compiler_4.1.1