scmeth Vignette
Contents
Introduction
Installation
Input files
Usage
4.1 Report
4.2 Functions
1. Introduction
Though a small chemical change in the genome, DNA methylation has significant impact in several diseases, developmental processes and other biological changes. Hence methylation data should be analyzed carefully to gain biological insights. scmeth package offers a few functions to assess the quality of the methylation data.
This bioconductor package contains functions to perform quality control and preprocessing analysis for single-cell methylation data. scmeth is especially customized to be used with the output from the FireCloud implementation of methylation pipeline. In addition to individual functions, report function in the package provides all inclusive report using most of the functions. If users prefer they can just use the report function to gain summary of their data.
2. Installation and package loading
scmeth is available in bioconductor and can be downloaded using the following commands
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("scmeth")Load the package
3. Input files
Main input for most of the function is a bsseq object. In the FireCloud implementation it is stored as hdf5 file which can be read via loadHDF5SummarizedExperiment function in HDF5Array package. Code chunk below shows how it can be loaded.
4. Usage
4.1 Report
A comprehensive quality control report can be generated in the package via report function. report function takes the bs object, the directory where the report should be saved, organism that this data is obtained from and the genomic build. Following is an example usage of the report function.
report function also takes two optional arguments: subSample and offset. With subsample parameter users can choose how many CpGs to consider in the analysis out of all the CpG information available. Our analysis have shown that conducting analysis with one million CpGs is suffice to capture the quality of the samples. However when more CpGs are added they reflect the full data precisely. Offset parameter is to avoid any telomere region when subsetting from the beginning of the data. Hence Offset parameter would avoid first n number of CpGs.
Command above generated an html report named qcReport.html. It will be stored in the indicated directory.
4.2 Functions
This section will elaborate on some of the main functions and show the usage of these functions based on a sample data set that comes along with the package.
scmeth package contains several functions to assess different metrics and success of the sequencing process.
coverage
One main metric is the CpG coverage. Coverage of the CpG can be assessed in different ways. Very basic one is to check how many CpG were observed in each sample. coverage function can be used to get this information.
Loading the data
directory <- system.file("extdata", "bismark_data", package='scmeth')
bsObject <- HDF5Array::loadHDF5SummarizedExperiment(directory)## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGAAGG_report.txt
## 756
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGGATG_report.txt
## 424
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_ATCAAG_report.txt
## 333readmetrics
Read information is important to assess whether sequencing and alignment succeeded. readmetrics function outputs a visualization showing number of reads seen in each samples and of those reads what proportion of them were mapped to the reference genome.
## sample total mapped unmapped
## 1 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGAAGG_report.txt 1145278 974438 170840
## 2 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGGATG_report.txt 763055 633756 129299
## 3 2017-08-02_HL2F5BBXX_2_AAGAGGCA_ATCAAG_report.txt 847927 717059 130868repmask
CpG Islands are characterized by their high GC content, high level of observed to expected ratio of CpGs and length over 500 bp. However some repeat regions in the genome also fit the same criteria although they are not bona fide CpG Island. Therefore it is important to see how many CpGs are observed in the non repeat regions of the genome. repMask functions provide information on the CpG coverage in non repeat regions of the genome. In order to build the repeat mask regions of the genome repmask function will require the organism and the genome build information.
library(BSgenome.Mmusculus.UCSC.mm10)
load(system.file("extdata", 'bsObject.rda', package='scmeth'))
repMask(bs, Mmusculus, "mm10")## coveredCpgs
## sc-RRBS_zyg_01_chr1 12208
## sc-RRBS_zyg_02_chr1 3056
## sc-RRBS_zyg_03_chr1 6666Coverage by Chromosome
There are several other ways the number of CpGs captured can be visualized. One of the way is to observe how the CpGs are distributed across different chromosomes. *chromosomeCoverage outputs CpG coverage by individual chromosomes.(Since the example data only contains information in chromosome 1 only the CpGs covered in chromosome 1 is shown.)
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGAAGG_report.txt
## 1 756
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGGATG_report.txt
## 1 424
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_ATCAAG_report.txt
## 1 333featureCoverage
Another way to observe the distribution of CpGs is to classify them by the genomic features they belong. Some of the features are very specific to the CpG dense regions such as CpG Islands, CpG Shores, CpG Shelves etc. Others are general genomic features such as introns, exons, promoters etc. This information can be obtained by featureCoverage function. In addition to the bs object this function requires the genomic features of interest and the genome build. Each element in the table represents the fraction of CpGs seen in particular cell in specific region compared to all the CpGs seen in that region.
cpgDensity
CpGs are not distributed across the genome uniformly. Most of the genome contains very low percentage of CpGs except for the CpG dense regions, i.e. CpG islands. Bisulfite sequencing targets all the CpGs across the genome, however reduced representation bisulfite sequencing (RRBS) target CpG dense CpG islands. Therefore CpG density plot will be a great diagnostic to see whether the protocol succeeded. In order to calculate the CpG density a window length should be specified. By default cpgDensity function chooses 1kB regions. Therefore CpG density plot can be used to check whether the protocol specifically targeted CpG dense or CpG sparse regions or whether CpGs were obtained uniformly across the regions.
downsample
In addition to the CpG coverage, methylation data can be assessed via down sampling analysis, methylation bias plot and methylation distribution. Down sampling analysis is a technique to assess whether the sequencing process achieved the saturation level in terms of CpG capture. In order to perform down sampling analysis CpGs that are covered at least will be sampled via binomial random sampling with given probability. At each probability level the number of CpGs captured is assessed. If the number of CpG captured attains a plateau then the sequencing was successful. downsample function provides a matrix of CpG coverage for each sample at various down sampling rates. The report renders this information into a plot. Downsampling rate ranges from 0.01 to 1, however users can change the downsampling rates.
mbiasPlot
Methylation bias plot shows the methylation along the reads. In a high quality samples methylation across the read would be more or less a horizontal line. However there could be fluctuations in the beginning or the end of the read due to the quality of the bases. Single cell sequencing samples also can show jagged trend in the methylation bias plot due to low read count. Methylation bias can be assessed via mbiasPlot function. This function takes the mbias file generated from FireCloud pipeline and generates the methylation bias plot.
methylationBiasFile <- '2017-04-21_HG23KBCXY_2_AGGCAGAA_TATCTC_pe.M-bias.txt'
mbiasList <- mbiasplot(mbiasFiles=system.file("extdata", methylationBiasFile,
package='scmeth'))
mbiasDf <- do.call(rbind.data.frame, mbiasList)
meanTable <- stats::aggregate(methylation ~ position + read, data=mbiasDf, FUN=mean)
sdTable <- stats::aggregate(methylation ~ position + read, data=mbiasDf, FUN=sd)
seTable <- stats::aggregate(methylation ~ position + read, data=mbiasDf, FUN=function(x){sd(x)/sqrt(length(x))})
sum_mt<-data.frame('position'=meanTable$position,'read'=meanTable$read,
'meth'=meanTable$methylation, 'sdMeth'=sdTable$methylation,
'seMeth'=seTable$methylation)
sum_mt$upperCI <- sum_mt$meth + (1.96*sum_mt$seMeth)
sum_mt$lowerCI <- sum_mt$meth - (1.96*sum_mt$seMeth)
sum_mt$read_rep <- paste(sum_mt$read, sum_mt$position, sep="_")
g <- ggplot2::ggplot(sum_mt)
g <- g + ggplot2::geom_line(ggplot2::aes_string(x='position', y='meth',
colour='read'))
g <- g + ggplot2::geom_ribbon(ggplot2::aes_string(ymin = 'lowerCI',
ymax = 'upperCI', x='position', fill = 'read'),
alpha=0.4)
g <- g + ggplot2::ylim(0,100) + ggplot2::ggtitle('Mbias Plot')
g <- g + ggplot2::ylab('methylation')
g
methylationDist
methylationDist function provides the methylation distribution of the samples. In this visualization methylation is divided into quantiles and ordered according to the cells with the lowest methylation to highest methylation. In single cell analysis almost all CpGs will be in the highest quantile or the lowest quantile. This visualization provides information on whether there are cells with intermediate methylation. Ideally , in single cell methylation most methylation should be either 1 or 0. If there are large number of intermediate methylation this indicates there might be some error in sequencing.
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGAAGG_report.txt
## 0.7841647
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGGATG_report.txt
## 0.6940448
## 2017-08-02_HL2F5BBXX_2_AAGAGGCA_ATCAAG_report.txt
## 0.8130466bsConversionPlot
Another important metric in methylation analysis is the bisulfite conversion rate. Bisulfite conversion rate indicates out of all the Cytosines in the non CpG context what fraction of them were methylated. Ideally this number should be 1 or 100% indicating none of the non CpG context cytosines are methylated. However in real data this will not be the case, yet bisulfite conversion rate below 95% indicates some problem with sample preparation. bsConversionPlot function generates a plot showing this metric for each sample.
## sample bsc
## 1 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGAAGG_report.txt 0.9940704
## 2 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGGATG_report.txt 0.9951570
## 3 2017-08-02_HL2F5BBXX_2_AAGAGGCA_ATCAAG_report.txt 0.9932266## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] BSgenome.Hsapiens.NCBI.GRCh38_1.3.1000
## [2] BSgenome.Mmusculus.UCSC.mm10_1.4.3
## [3] BSgenome_1.75.0
## [4] rtracklayer_1.67.0
## [5] BiocIO_1.17.1
## [6] Biostrings_2.75.3
## [7] XVector_0.47.0
## [8] GenomicRanges_1.59.1
## [9] GenomeInfoDb_1.43.2
## [10] IRanges_2.41.2
## [11] S4Vectors_0.45.2
## [12] BiocGenerics_0.53.3
## [13] generics_0.1.3
## [14] scmeth_1.27.0
## [15] rmarkdown_2.29
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 bitops_1.0-9
## [3] bsseq_1.43.1 permute_0.9-7
## [5] rlang_1.1.4 magrittr_2.0.3
## [7] matrixStats_1.4.1 compiler_4.4.2
## [9] RSQLite_2.3.9 DelayedMatrixStats_1.29.0
## [11] GenomicFeatures_1.59.1 png_0.1-8
## [13] vctrs_0.6.5 reshape2_1.4.4
## [15] stringr_1.5.1 pkgconfig_2.0.3
## [17] crayon_1.5.3 fastmap_1.2.0
## [19] dbplyr_2.5.0 labeling_0.4.3
## [21] Rsamtools_2.23.1 tzdb_0.4.0
## [23] UCSC.utils_1.3.0 purrr_1.0.2
## [25] bit_4.5.0.1 xfun_0.49
## [27] beachmat_2.23.5 zlibbioc_1.52.0
## [29] cachem_1.1.0 jsonlite_1.8.9
## [31] blob_1.2.4 rhdf5filters_1.19.0
## [33] DelayedArray_0.33.3 Rhdf5lib_1.29.0
## [35] BiocParallel_1.41.0 parallel_4.4.2
## [37] R6_2.5.1 bslib_0.8.0
## [39] stringi_1.8.4 limma_3.63.2
## [41] jquerylib_0.1.4 Rcpp_1.0.13-1
## [43] SummarizedExperiment_1.37.0 knitr_1.49
## [45] R.utils_2.12.3 readr_2.1.5
## [47] Matrix_1.7-1 tidyselect_1.2.1
## [49] abind_1.4-8 yaml_2.3.10
## [51] codetools_0.2-20 curl_6.0.1
## [53] lattice_0.22-6 tibble_3.2.1
## [55] regioneR_1.39.0 plyr_1.8.9
## [57] withr_3.0.2 Biobase_2.67.0
## [59] KEGGREST_1.47.0 evaluate_1.0.1
## [61] BiocFileCache_2.15.0 pillar_1.10.0
## [63] BiocManager_1.30.25 filelock_1.0.3
## [65] MatrixGenerics_1.19.0 DT_0.33
## [67] vroom_1.6.5 RCurl_1.98-1.16
## [69] BiocVersion_3.21.1 hms_1.1.3
## [71] ggplot2_3.5.1 sparseMatrixStats_1.19.0
## [73] munsell_0.5.1 scales_1.3.0
## [75] gtools_3.9.5 glue_1.8.0
## [77] maketools_1.3.1 tools_4.4.2
## [79] AnnotationHub_3.15.0 data.table_1.16.4
## [81] sys_3.4.3 locfit_1.5-9.10
## [83] GenomicAlignments_1.43.0 buildtools_1.0.0
## [85] XML_3.99-0.17 rhdf5_2.51.1
## [87] grid_4.4.2 crosstalk_1.2.1
## [89] AnnotationDbi_1.69.0 colorspace_2.1-1
## [91] GenomeInfoDbData_1.2.13 HDF5Array_1.35.2
## [93] restfulr_0.0.15 annotatr_1.33.0
## [95] cli_3.6.3 rappdirs_0.3.3
## [97] S4Arrays_1.7.1 dplyr_1.1.4
## [99] gtable_0.3.6 R.methodsS3_1.8.2
## [101] sass_0.4.9 digest_0.6.37
## [103] SparseArray_1.7.2 farver_2.1.2
## [105] rjson_0.2.23 htmlwidgets_1.6.4
## [107] R.oo_1.27.0 memoise_2.0.1
## [109] htmltools_0.5.8.1 lifecycle_1.0.4
## [111] httr_1.4.7 mime_0.12
## [113] statmod_1.5.0 bit64_4.5.2