Contents

The DMCFB package is a pipeline to identify differentially methyalted cytosines (DMC) in bisulfite sequencing data using Bayesian functional regression models. In what follows we provides some guidlines on how to read your data and analyze them.

1 Reading data

1.1 Reading bisulfite data (using files)

The R-method readBismark() is used to read bisulfite data files that are created by Bismark. Each file must include six columns, with no header, that represent

  • Chromosome
  • Start position in the chromosome
  • End position in the chromosome
  • Methylation level (m/(m+u))
  • Number of methylated reads (m)
  • Number of un-methylated reads (u)

and each row is a cytonsine (or a small region) in DNA.

The function readBismark(<files' paths>, <files' names>) has two inputs: ‘the paths of the files’ and ‘the names of the files’. Using this function an object of class BSDMC is created. Extra information about data such as Age, Gender, Group, etc, must be assigned to the object using DataFrame function. As an example, we have provided three files in the package that can be read as follows:

library(DMCFB)
#> Loading required package: SummarizedExperiment
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: parallel
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:parallel':
#> 
#>     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
#>     clusterExport, clusterMap, parApply, parCapply, parLapply,
#>     parLapplyLB, parRapply, parSapply, parSapplyLB
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, intersect,
#>     is.unsorted, lapply, mapply, match, mget, order, paste, pmax,
#>     pmax.int, pmin, pmin.int, rank, rbind, rownames, sapply,
#>     setdiff, sort, table, tapply, union, unique, unsplit, which,
#>     which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:base':
#> 
#>     expand.grid
#> Loading required package: IRanges
#> 
#> Attaching package: 'IRanges'
#> The following object is masked from 'package:grDevices':
#> 
#>     windows
#> Loading required package: GenomeInfoDb
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> Loading required package: DelayedArray
#> Loading required package: matrixStats
#> 
#> Attaching package: 'matrixStats'
#> The following objects are masked from 'package:Biobase':
#> 
#>     anyMissing, rowMedians
#> Loading required package: BiocParallel
#> 
#> Attaching package: 'DelayedArray'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
#> The following objects are masked from 'package:base':
#> 
#>     aperm, apply, rowsum
#> DMCFB package, Version 1.0.0, Released 2019-10-29
#> DMCFB is a pipeline for identifying differentially methylated 
#>     cytosines using a Bayesian functional regression model in bisulfite 
#>     sequencing data. By using a functional regression data model, it tries to 
#>     capture position-specific, group-specific and other covariates-specific 
#>     methylation patterns as well as spatial correlation patterns and unknown 
#>     underlying models of methylation data. It is robust and flexible with 
#>     respect to the true underlying models and inclusion of any covariates, and 
#>     the missing values are imputed using spatial correlation between positions 
#>     and samples. A Bayesian approach is adopted for estimation and inference in 
#>     the proposed method.
#> BugReports: https://github.com/shokoohi/DMCFB/issues
#> 
#> Attaching package: 'DMCFB'
#> The following object is masked from 'package:Biobase':
#> 
#>     combine
#> The following object is masked from 'package:BiocGenerics':
#> 
#>     combine
fn <- list.files(system.file("extdata",package = "DMCFB"))
fn.f <- list.files(system.file("extdata",package="DMCFB"), full.names=TRUE)
OBJ <- readBismark(fn.f, fn)
#> Processing sample blk.BCU1568_BC_BS_1 ... 
#> Processing sample blk.BCU173_TC_BS_1 ... 
#> Processing sample blk.BCU551_Mono_BS_1 ... 
#> Building BSDMC object.
cdOBJ <- DataFrame(Cell = factor(c("BC", "TC","Mono"),
levels = c("BC", "TC", "Mono")), row.names = c("BCU1568","BCU173","BCU551"))
colData(OBJ) <- cdOBJ
OBJ
#> class: BSDMC 
#> dim: 25668 3 
#> metadata(0):
#> assays(3): methReads totalReads methLevels
#> rownames(25668): 1 2 ... 25667 25668
#> rowData names(0):
#> colnames(3): BCU1568 BCU173 BCU551
#> colData names(1): Cell

1.2 Reading bisulfite data (using matrices)

Alternatively, one can use two integer matrices and a DataFrame to create BSDMC object using cBSDMC() function. One matrix includes the read-depth and the other one includes methylation reads. The columns of these matrices represent samples and the rows represent cytosine positions.

Additional information about the genomic positions and covariates must be stored in a DataFrame and then assign to the object.

The following exampel shows the details.

library(DMCFB)
set.seed(1980)
nr <- 1000
nc <- 8
metht <- matrix(as.integer(runif(nr * nc, 0, 100)), nr)
methc <- matrix(rbinom(n=nr*nc,c(metht),prob = runif(nr*nc)),nr,nc)
methl <- methc/metht
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])
OBJ2 <- cBSDMC(rowRanges=r1,methReads=methc,totalReads=metht,
  methLevels=methl,colData=cd1)
OBJ2
#> class: BSDMC 
#> dim: 1000 8 
#> metadata(0):
#> assays(3): methReads totalReads methLevels
#> rownames(1000): 1 2 ... 999 1000
#> rowData names(0):
#> colnames(8): A B ... G H
#> colData names(1): Group

2 Identifying DMCs

To identify DMCs, one need to use the function findDMCFB() function. The function

library(DMCFB)
start.time <- Sys.time()
path0 <- "..//BCData/" # provide the path to the files
namelist.new <- list.files(path0,pattern="blk",full.names=F)
namelist.new.f <- list.files(path0,pattern="blk",full.names=T)
type <- NULL
for(i in seq_along(namelist.new)){
    type[i] <- unlist(strsplit(namelist.new[i], split=c('_'), fixed=TRUE))[2]
}
type
table(type)
indTC <- which(type=="TC")
indBC <- which(type=="BC")
indMono <- which(type=="Mono")
namelist.new <- namelist.new[c(indBC,indMono,indTC)]
namelist.new.f <- namelist.new.f[c(indBC,indMono,indTC)]
BLKDat <- readBismark(namelist.new.f, namelist.new)
colData1 <- DataFrame(Group = factor(
  c(rep("BC",length(indBC)), rep("Mono",length(indMono)), 
  rep("TC", length(indTC))), levels = c("BC", "Mono", "TC")), 
  row.names = colnames(BLKData))
colData(BLKDat) <- colData1
BLK.BC.Mono.TC <- sort(BLKDat)
DMC.obj = findDMCFB(object = BLKDat, bwa = 30, bwb = 30, nBurn = 300, nMC = 300,
  nThin = 1, alpha = 5e-5, pSize = 500, sfiles = FALSE)

3 Figures

To plot DMCs one can use the plotDMCFB() function to plot an BSDMC object that resulted from running findDMCFB() function. To illustrate use the following example:

library(DMCFB)
set.seed(1980)
nr <- 1000
nc <- 8
metht <- matrix(as.integer(runif(nr * nc, 0, 100)), nr)
methc <- matrix(rbinom(n=nr*nc,c(metht),prob = runif(nr*nc)),nr,nc)
methl <- methc/metht
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 <- cBSDMC(rowRanges=r1,methReads=methc,totalReads=metht,
  methLevels=methl,colData=cd1)
OBJ2 = findDMCFB(object = OBJ1, bwa = 30, bwb = 30, nBurn = 10, nMC = 10,
  nThin = 1, alpha = 0.05, pSize = 500, sfiles = FALSE)
#> ------------------------------------------------------------
#> Running Bayesian functional regression model ...
#> The priors's SD = 0.3027, estimated from data ...
#> Number of assigned cores: 1 ...
#> ------------------------------------------------------------
#> Fitted model:
#> logit(MethRead/ReadDepth) ~ F(Group)
#> ------------------------------------------------------------
#> Creating 2 batches of genomic positions ...
#> Running batch 1/2; chr1; 500 positions; Region [   1,  500]; Date 2019-10-30 01:07:08
#> Running batch 2/2; chr1; 500 positions; Region [ 501, 1000]; Date 2019-10-30 01:07:13
#> ------------------------------------------------------------
#> Combining 2 objects ...
#> Objects are combined.
#> ------------------------------------------------------------
#> Identifying DMCs ...
#> DMCs are identified.
#> ------------------------------------------------------------
#> Percentage of non-DMCs and DMCs:
#> Equal(%)   DMC(%) 
#>     23.8     76.2
#> ------------------------------------------------------------
#> Percentage of hyper-, hypo-, and equal-methylated positions:
#>        Equal(%) Hyper(%) Hypo(%)
#> G2vsG1        0        0       0
#> ------------------------------------------------------------
plotDMCFB(OBJ2, region = c(1,400), nSplit = 2)

4 Session info

sessionInfo()
#> R version 3.6.1 (2019-07-05)
#> 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] parallel  stats4    stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] DMCFB_1.0.0                 SummarizedExperiment_1.16.0
#>  [3] DelayedArray_0.12.0         BiocParallel_1.20.0        
#>  [5] matrixStats_0.55.0          Biobase_2.46.0             
#>  [7] GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
#>  [9] IRanges_2.20.0              S4Vectors_0.24.0           
#> [11] BiocGenerics_0.32.0         BiocStyle_2.14.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.2               lattice_0.20-38         
#>  [3] Rsamtools_2.2.0          Biostrings_2.54.0       
#>  [5] assertthat_0.2.1         digest_0.6.22           
#>  [7] foreach_1.4.7            R6_2.4.0                
#>  [9] evaluate_0.14            coda_0.19-3             
#> [11] httr_1.4.1               pillar_1.4.2            
#> [13] zlibbioc_1.32.0          rlang_0.4.1             
#> [15] data.table_1.12.6        minqa_1.2.4             
#> [17] nloptr_1.2.1             Matrix_1.2-17           
#> [19] rmarkdown_1.16           splines_3.6.1           
#> [21] lme4_1.1-21              stringr_1.4.0           
#> [23] fastDummies_1.5.0        RCurl_1.95-4.12         
#> [25] compiler_3.6.1           rtracklayer_1.46.0      
#> [27] xfun_0.10                pkgconfig_2.0.3         
#> [29] speedglm_0.3-2           htmltools_0.4.0         
#> [31] tidyselect_0.2.5         tibble_2.1.3            
#> [33] GenomeInfoDbData_1.2.2   bookdown_0.14           
#> [35] codetools_0.2-16         XML_3.98-1.20           
#> [37] crayon_1.3.4             dplyr_0.8.3             
#> [39] GenomicAlignments_1.22.0 MASS_7.3-51.4           
#> [41] bitops_1.0-6             grid_3.6.1              
#> [43] nlme_3.1-141             arm_1.10-1              
#> [45] magrittr_1.5             stringi_1.4.3           
#> [47] XVector_0.26.0           doParallel_1.0.15       
#> [49] benchmarkmeData_1.0.2    boot_1.3-23             
#> [51] iterators_1.0.12         tools_3.6.1             
#> [53] benchmarkme_1.0.2        glue_1.3.1              
#> [55] purrr_0.3.3              abind_1.4-5             
#> [57] yaml_2.2.0               BiocManager_1.30.9      
#> [59] knitr_1.25