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
## 1 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGAAGG_report.txt 1145278 974438
## 2 2017-08-02_HL2F5BBXX_2_AAGAGGCA_AGGATG_report.txt 763055 633756
## 3 2017-08-02_HL2F5BBXX_2_AAGAGGCA_ATCAAG_report.txt 847927 717059
## unmapped
## 1 170840
## 2 129299
## 3 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 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so
##
## 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
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] BSgenome.Hsapiens.NCBI.GRCh38_1.3.1000
## [2] BSgenome.Mmusculus.UCSC.mm10_1.4.0
## [3] BSgenome_1.52.0
## [4] rtracklayer_1.44.0
## [5] Biostrings_2.52.0
## [6] XVector_0.24.0
## [7] GenomicRanges_1.36.0
## [8] GenomeInfoDb_1.20.0
## [9] IRanges_2.18.0
## [10] S4Vectors_0.22.0
## [11] BiocGenerics_0.30.0
## [12] scmeth_1.4.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-6 matrixStats_0.54.0
## [3] bit64_0.9-7 bsseq_1.20.0
## [5] progress_1.2.0 httr_1.4.0
## [7] tools_3.6.0 R6_2.4.0
## [9] DT_0.5 HDF5Array_1.12.0
## [11] DBI_1.0.0 lazyeval_0.2.2
## [13] colorspace_1.4-1 permute_0.9-5
## [15] tidyselect_0.2.5 prettyunits_1.0.2
## [17] bit_1.1-14 curl_3.3
## [19] compiler_3.6.0 Biobase_2.44.0
## [21] DelayedArray_0.10.0 labeling_0.3
## [23] scales_1.0.0 readr_1.3.1
## [25] rappdirs_0.3.1 stringr_1.4.0
## [27] digest_0.6.18 Rsamtools_2.0.0
## [29] R.utils_2.8.0 rmarkdown_1.12
## [31] pkgconfig_2.0.2 htmltools_0.3.6
## [33] dbplyr_1.4.0 limma_3.40.0
## [35] regioneR_1.16.0 htmlwidgets_1.3
## [37] rlang_0.3.4 RSQLite_2.1.1
## [39] shiny_1.3.2 DelayedMatrixStats_1.6.0
## [41] jsonlite_1.6 crosstalk_1.0.0
## [43] BiocParallel_1.18.0 gtools_3.8.1
## [45] R.oo_1.22.0 dplyr_0.8.0.1
## [47] RCurl_1.95-4.12 magrittr_1.5
## [49] GenomeInfoDbData_1.2.1 Matrix_1.2-17
## [51] Rcpp_1.0.1 munsell_0.5.0
## [53] Rhdf5lib_1.6.0 R.methodsS3_1.7.1
## [55] stringi_1.4.3 yaml_2.2.0
## [57] SummarizedExperiment_1.14.0 zlibbioc_1.30.0
## [59] rhdf5_2.28.0 plyr_1.8.4
## [61] BiocFileCache_1.8.0 AnnotationHub_2.16.0
## [63] grid_3.6.0 blob_1.1.1
## [65] promises_1.0.1 crayon_1.3.4
## [67] lattice_0.20-38 GenomicFeatures_1.36.0
## [69] hms_0.4.2 locfit_1.5-9.1
## [71] knitr_1.22 pillar_1.3.1
## [73] reshape2_1.4.3 biomaRt_2.40.0
## [75] XML_3.98-1.19 glue_1.3.1
## [77] evaluate_0.13 annotatr_1.10.0
## [79] data.table_1.12.2 BiocManager_1.30.4
## [81] httpuv_1.5.1 gtable_0.3.0
## [83] purrr_0.3.2 assertthat_0.2.1
## [85] ggplot2_3.1.1 xfun_0.6
## [87] mime_0.6 xtable_1.8-4
## [89] later_0.8.0 tibble_2.1.1
## [91] GenomicAlignments_1.20.0 AnnotationDbi_1.46.0
## [93] memoise_1.1.0 interactiveDisplayBase_1.22.0