Contents

library(methodical)
library(TumourMethData)
library(BSgenome.Hsapiens.UCSC.hg19)

0.1 Introduction

Most functions from methodical take as input RangedSummarizedExperiment objects with methylation data. If there are many samples, there can be many millions or, in the case of WGBS data, even billions of data points within a DNA methylation dataset. It can be unfeasible to load all this data into memory at once. This problem can be overcome by using DelayedArrays backed by HDF5 files, enabling data to be read into memory only as needed. Methodical provides a suite of functions for working with such DNA methylation RangedSummarizedExperiments, including functions to extract methylation values for sites overlapping genomic regions of interest GRanges, to liftover the methylation sites from one genome build to another and to mask methylation sites overlapping certain genomic regions, e.g.  repeats.

Here we demonstrate this different functionality using a dataset downloaded from TumourMethData.

# Download RangedSummarizedExperiment with methylation data for prostate metastases from TumourMethData
mcrpc_wgbs_hg38_chr11 = TumourMethData::download_meth_dataset(dataset = "mcrpc_wgbs_hg38_chr11")
#> [1] "A HDF5 SummarizedExperiment is already present in /tmp/RtmppevoFB/mcrpc_wgbs_hg38_chr11 and is being returned"

0.2 Extracting methylation data from methylation RangedSummarizedExperiments

We’ll first demonstrate how to extract methylation data for individual CpG sites in mcrpc_wgbs_hg38_chr11 overlapping a supplied GRanges object with extractGRangesMethSiteValues, using the gene body as an example.

# Create a GRanges with the hg38 genomic coordinates for the GSTP1, including 
# 2 KB upstream of its designated start in Ensembl
gstp1_start_site_region <- GRanges("chr11:67581742-67586656:+")

# Extract methylation values for CpG sites overlapping GSTP1 gene body
gstp1_cpg_methylation <- extractGRangesMethSiteValues(
  meth_rse = mcrpc_wgbs_hg38_chr11, genomic_regions = gstp1_start_site_region)

# View the first few rows and columns of the result. 
# extractGRangesMethSiteValues returns a row for each methylation site and a 
# separate column for each sample where row names give the coordinates of the 
# methylation sites in character format. 
gstp1_cpg_methylation[1:6, 1:6]
#>                  DTB_003   DTB_005    DTB_008    DTB_009    DTB_011   DTB_018
#> chr11:67581759 0.1764706 0.2682927 0.23076923 0.10714286 0.09677419 0.7045455
#> chr11:67581799 0.1944444 0.1470588 0.12500000 0.10344828 0.04166667 0.4871795
#> chr11:67581838 0.3157895 0.1538462 0.15789474 0.13333333 0.10000000 0.6923077
#> chr11:67581849 0.1388889 0.0000000 0.07692308 0.06896552 0.16129032 0.4000000
#> chr11:67581866 0.3488372 0.2000000 0.17073171 0.08823529 0.29032258 0.7575758
#> chr11:67581875 0.1590909 0.0000000 0.15384615 0.08823529 0.25000000 0.3142857

Next we’ll show how to summarize the methylation of CpGs over regions of interest with summarizeRegionMethylation, using CpG islands as an example. summarizeRegionMethylation processes the the supplied genomic regions in chunks so that the methylation data for CpG sites overlapping the regions of interest is not all loaded into memory at once. The parameter max_sites_per_chunk controls the approximate number of CpG sites maximally read into memory at once and defaults to floor(62500000/ncol(meth_rse). Several chunks can be processed in parallel using BiocParallel via the BPPARAM argument which takes a BiocParallelParam object. The number of workers indicated by BiocParallelParam determines the number of chunks that will be processed in parallel. Some experimentation may be needed to find the optimal choices for max_sites_per_chunk and the number of workers in terms of speed and memory usage.

# Load CpG islands annotation for hg38 
cpg_island_annotation <- annotatr::build_annotations(genome = "hg38", annotation = "hg38_cpgs")
#> Building CpG islands...
#> Building CpG shores...
#> Building CpG shelves...
#> Building inter-CpG-islands...
names(cpg_island_annotation) <- cpg_island_annotation$id

# Filter for annotation for chr11
cpg_island_annotation = cpg_island_annotation[seqnames(cpg_island_annotation) == "chr11"]

# Convert into a GRangesList with separate GRanges for islands, shores, shelves and inter island regions
cpg_island_annotation <- GRangesList(split(cpg_island_annotation, cpg_island_annotation$type))

# Create a BPPARAM class
BPPARAM = BiocParallel::bpparam()

# Summarize methylation levels for CpG islands
cpg_island_methylation <- summarizeRegionMethylation(
  meth_rse = mcrpc_wgbs_hg38_chr11, genomic_regions = cpg_island_annotation$hg38_cpg_islands, 
  BPPARAM = BPPARAM, summary_function = colMeans)
#> There are some seqlevels from genomic_regions missing from meth_rse
#> Summarizing genomic region methylation

# Print a few rows for the first few samples of the result
cpg_island_methylation[1000:1006, 1:6]
#>                  DTB_003     DTB_005     DTB_008     DTB_009     DTB_011
#> island:14913 0.124852221 0.009667637 0.008330832 0.014607967 0.039786356
#> island:14914 0.009844099 0.029795906 0.006258614 0.072477345 0.052693357
#> island:14915 0.001953125 0.004184060 0.000000000 0.000000000 0.003801843
#> island:14916 0.006799628 0.002174912 0.006262350 0.003554305 0.007567667
#> island:14917 0.901040137 0.834174923 0.819908342 0.865086671 0.806623071
#> island:14918 0.940876775 0.940707064 0.950919937 0.959274590 0.907710156
#> island:14919 0.026562413 0.024696883 0.019078006 0.017051253 0.034368160
#>                  DTB_018
#> island:14913 0.012036163
#> island:14914 0.050911768
#> island:14915 0.002777778
#> island:14916 0.004661724
#> island:14917 0.838317611
#> island:14918 0.930573772
#> island:14919 0.022032495

0.3 Plotting methylation values in genomic regions and annotating genomic regions in plots

We can plot methylation values extracted from a genomic region for a single sample using the plotMethylationValues() function. We’ll demonstrate this using the values we extracted in the region surrounding the 5’ end of GSTP1 for the DTB_003 prostate metastasis sample. We indicate the sample we want to plot with the sample_name parameter.

# Plot the methylation values along the GSTP1 gene body for one prostate metastasis sample.
gstp1_methylation_plot = plotMethylationValues(gstp1_cpg_methylation, sample_name = "DTB_003")
print(gstp1_methylation_plot)

Additionally, we can also annotate our plots using the annotatePlot() function. It uses a GRangeList provided with the annotation_grl parameter to create an annotation plot showing the regions in the GRangesList which overlap the genomic region displayed in the main plot. Each of the GRanges objects making up the GRangesList is given a different colour in the annotation plot and the names of these component GRanges are indicated. We can control the colours used with the grl_colours parameter if we don’t want to use the default colours.

If we provide a GRanges object with the location of a transcription start site to the reference_tss, parameter, we can show the distance of methylation sites upstream and downstream of this.

By default, the main plot and annotation plot are combined into a single plot and returned. The annotation_plot_proportion parameter sets the proportion of the total plot height dedicated to the annotation plot. We can instead return the annotation plot by itself by setting the annotation_plot_only parameter to TRUE.

We’ll annotate the location of CpG islands, CpG shores, CpG shores and inter CpG island regions for gstp1_methylation_plot using the cpg_island_annotation GRangesList we created.

# Annotate gstp1_methylation_plot with cpg_island_annotation
annotatePlot(meth_site_plot = gstp1_methylation_plot, 
  annotation_grl = cpg_island_annotation, annotation_plot_proportion = 0.3, 
  grl_colours = c("#DEEBF7", "#9ECAE1", "#4292C6", "#08519C"))
#> Warning in get_plot_component(plot, "guide-box"): Multiple components found;
#> returning the first one. To return all, use `return_all = TRUE`.

# Create same plot, except showing the distance to the GSTP1 start site on the x-axis
annotatePlot(meth_site_plot = gstp1_methylation_plot, 
  annotation_grl = cpg_island_annotation, 
  reference_tss = GRanges("chr11:67583742"),annotation_plot_proportion = 0.3, 
  grl_colours = c("#DEEBF7", "#9ECAE1", "#4292C6", "#08519C"))
#> Warning in get_plot_component(plot, "guide-box"): Multiple components found;
#> returning the first one. To return all, use `return_all = TRUE`.

# Return the annotation plot by itself
annotatePlot(meth_site_plot = gstp1_methylation_plot, 
  annotation_grl = cpg_island_annotation, annotation_plot_proportion = 0.3, 
  grl_colours = c("#DEEBF7", "#9ECAE1", "#4292C6", "#08519C"), annotation_plot_only = TRUE)

0.4 Masking regions in a methylation RangedSummarizedExperiments

We may want to mask certain regions in a methylation RangedSummarizedExperiment. With the maskRangesInRSE() function, we can mask regions across all samples (which could be useful for repetitive regions or polymorphic regions) or on a sample by sample basis (which could be appropriate for different regions that are known to be mutated in different tumour samples). All methylation values within the masked regions will be set to NA.

The mask_ranges argument takes either a GRanges or GRangesList with the regions that should be masked. If a GRanges object is provided, all methylation sites overlapping these regions will be masked across all samples. If a GRangesList is provided, the names of the component GRanges should match sample names in the RangedSummarizedExperiment and in each sample, the methylation sites overlapping the regions in its corresponding GRangesList entry will be masked.

We will demonstrate how to mask LINE repetitive regions across all samples.

# Download repetitive sequences from AnnotationHub and filter for LINE elements
repeat_annotation_hg38 <- AnnotationHub::AnnotationHub()[["AH99003"]]
#> loading from cache
line_elements_hg38 <- repeat_annotation_hg38[repeat_annotation_hg38$repClass == "SINE"]

# Mask LINE elements in mcrpc_wgbs_hg38_chr11
mcrpc_wgbs_hg38_chr11_lines_masked <- maskRangesInRSE(rse = mcrpc_wgbs_hg38_chr11, 
  mask_ranges = line_elements_hg38)

# Extract the methylation values for one of the LINE elements in the 
# unmasked and masked RSE
extractGRangesMethSiteValues(meth_rse = mcrpc_wgbs_hg38_chr11, 
  genomic_regions = line_elements_hg38[1000])[, 1:6]
#> [1] DTB_003 DTB_005 DTB_008 DTB_009 DTB_011 DTB_018
#> <0 rows> (or 0-length row.names)
extractGRangesMethSiteValues(meth_rse = mcrpc_wgbs_hg38_chr11_lines_masked, 
  genomic_regions = line_elements_hg38[1000])[, 1:6]
#> [1] DTB_003 DTB_005 DTB_008 DTB_009 DTB_011 DTB_018
#> <0 rows> (or 0-length row.names)

0.5 Lifting over methylation sites in a RangedSummarizedExperiment from one genome build to another

Sometimes we may want to work with a different genome build to that used to construct a methylation RangedSummarizedExperiment. We can easily liftover the genomic coordinates of the methylation sites using liftoverMethRSE() function. To do this, we need a liftover chain file for the appropriate source and target genome builds and which we provide to the chain argument.

All methylation sites which cannot be mapped to the target genome build and those which result in many-to-one mappings are removed. We also need to decide whether we want to remove methylation sites in the source genome build which map to multiple sites in the target genome build. We do this using the remove_one_to_many_mapping argument, which has a default value of TRUE. We can also remove any regions which do not map to desired regions in the target genome, for example CpG sites, by providing a GRanges object to the argument permitted_target_regions.

We will demonstrate how to liftover mcrpc_wgbs_hg38_chr11 to hg19.

# Create a DNAStringSet for chromosome11
chr11_dss = setNames(DNAStringSet(BSgenome.Hsapiens.UCSC.hg19[["chr11"]]), "chr11")

# Get CpG sites for hg19 for chromsome 11
hg19_cpgs <- methodical::extractMethSitesFromGenome(genome = chr11_dss)

# Download hg38 to hg19 liftover chain from AnnotationHub
hg38tohg19Chain <- AnnotationHub::AnnotationHub()[["AH14108"]]
#> loading from cache

# Liftover mcrpc_wgbs_hg38_chr11 to mcrpc_wgbs_hg19_chr11
mcrpc_wgbs_hg19_chr11 <- liftoverMethRSE(meth_rse = mcrpc_wgbs_hg38_chr11, chain = hg38tohg19Chain, 
  remove_one_to_many_mapping = TRUE, permitted_target_regions = hg19_cpgs)
#> 42343 non-mapping sites removed
#> 0 one-to-many mapping sites removed
#> 1773 many-to-one mapping sites removed
#> 46319 sites not overlapping permitted target regions removed

# Compare the dimensions of mcrpc_wgbs_hg38_chr11 and mcrpc_wgbs_hg19_chr11. 
# 1,423,050 methylation sites could not be lifted over from hg38 to hg19. 
dim(mcrpc_wgbs_hg38_chr11)
#> [1] 1333114     100
dim(mcrpc_wgbs_hg19_chr11)
#> [1] 1286553     100

# chr1:921635 should be lifted over to chr1:857015 so confirm that they have 
# the same methylation values in hg38 and hg19
rtracklayer::liftOver(GRanges("chr11:67581759"), hg38tohg19Chain)
#> GRangesList object of length 1:
#> [[1]]
#> GRanges object with 1 range and 0 metadata columns:
#>       seqnames    ranges strand
#>          <Rle> <IRanges>  <Rle>
#>   [1]    chr11  67349230      *
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths
extractGRangesMethSiteValues(mcrpc_wgbs_hg38_chr11, GRanges("chr11:67581759"))[, 1:8]
#>                  DTB_003   DTB_005   DTB_008   DTB_009    DTB_011   DTB_018
#> chr11:67581759 0.1764706 0.2682927 0.2307692 0.1071429 0.09677419 0.7045455
#>                  DTB_019   DTB_020
#> chr11:67581759 0.1538462 0.1612903
extractGRangesMethSiteValues(mcrpc_wgbs_hg19_chr11, GRanges("chr11:67349230"))[, 1:8]
#>                  DTB_003   DTB_005   DTB_008   DTB_009    DTB_011   DTB_018
#> chr11:67349230 0.1764706 0.2682927 0.2307692 0.1071429 0.09677419 0.7045455
#>                  DTB_019   DTB_020
#> chr11:67349230 0.1538462 0.1612903

0.6 SessionInfo

sessionInfo()
#> R version 4.4.0 beta (2024-04-15 r86425)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.72.0                  
#>  [3] BiocIO_1.14.0                     Biostrings_2.72.0                
#>  [5] XVector_0.44.0                    HDF5Array_1.32.0                 
#>  [7] DelayedArray_0.30.0               SparseArray_1.4.0                
#>  [9] S4Arrays_1.4.0                    abind_1.4-5                      
#> [11] Matrix_1.7-0                      rtracklayer_1.64.0               
#> [13] AnnotationHub_3.12.0              BiocFileCache_2.12.0             
#> [15] dbplyr_2.5.0                      rhdf5_2.48.0                     
#> [17] TumourMethData_1.1.0              DESeq2_1.44.0                    
#> [19] methodical_1.0.0                  SummarizedExperiment_1.34.0      
#> [21] Biobase_2.64.0                    MatrixGenerics_1.16.0            
#> [23] matrixStats_1.3.0                 ggplot2_3.5.1                    
#> [25] GenomicRanges_1.56.0              GenomeInfoDb_1.40.0              
#> [27] IRanges_2.38.0                    S4Vectors_0.42.0                 
#> [29] BiocGenerics_0.50.0               BiocStyle_2.32.0                 
#> 
#> loaded via a namespace (and not attached):
#>  [1] jsonlite_1.8.8           magrittr_2.0.3           magick_2.8.3            
#>  [4] GenomicFeatures_1.56.0   farver_2.1.1             rmarkdown_2.26          
#>  [7] zlibbioc_1.50.0          vctrs_0.6.5              memoise_2.0.1           
#> [10] Rsamtools_2.20.0         RCurl_1.98-1.14          RcppRoll_0.3.0          
#> [13] tinytex_0.50             htmltools_0.5.8.1        curl_5.2.1              
#> [16] Rhdf5lib_1.26.0          sass_0.4.9               bslib_0.7.0             
#> [19] plyr_1.8.9               cachem_1.0.8             GenomicAlignments_1.40.0
#> [22] mime_0.12                lifecycle_1.0.4          iterators_1.0.14        
#> [25] pkgconfig_2.0.3          R6_2.5.1                 fastmap_1.1.1           
#> [28] GenomeInfoDbData_1.2.12  digest_0.6.35            colorspace_2.1-0        
#> [31] AnnotationDbi_1.66.0     ExperimentHub_2.12.0     regioneR_1.36.0         
#> [34] RSQLite_2.3.6            filelock_1.0.3           labeling_0.4.3          
#> [37] fansi_1.0.6              httr_1.4.7               compiler_4.4.0          
#> [40] bit64_4.0.5              withr_3.0.0              BiocParallel_1.38.0     
#> [43] DBI_1.2.2                R.utils_2.12.3           rappdirs_0.3.3          
#> [46] rjson_0.2.21             tools_4.4.0              R.oo_1.26.0             
#> [49] glue_1.7.0               restfulr_0.0.15          rhdf5filters_1.16.0     
#> [52] grid_4.4.0               reshape2_1.4.4           generics_0.1.3          
#> [55] gtable_0.3.5             tzdb_0.4.0               R.methodsS3_1.8.2       
#> [58] data.table_1.15.4        hms_1.1.3                utf8_1.2.4              
#> [61] BiocVersion_3.19.1       foreach_1.5.2            pillar_1.9.0            
#> [64] stringr_1.5.1            vroom_1.6.5              dplyr_1.1.4             
#> [67] lattice_0.22-6           bit_4.0.5                tidyselect_1.2.1        
#> [70] locfit_1.5-9.9           knitr_1.46               bookdown_0.39           
#> [73] xfun_0.43                annotatr_1.30.0          stringi_1.8.3           
#> [76] UCSC.utils_1.0.0         yaml_2.3.8               evaluate_0.23           
#> [79] codetools_0.2-20         tibble_3.2.1             BiocManager_1.30.22     
#> [82] cli_3.6.2                munsell_0.5.1            jquerylib_0.1.4         
#> [85] Rcpp_1.0.12              png_0.1-8                XML_3.99-0.16.1         
#> [88] parallel_4.4.0           readr_2.1.5              blob_1.2.4              
#> [91] bitops_1.0-7             scales_1.3.0             purrr_1.0.2             
#> [94] crayon_1.5.2             rlang_1.1.3              cowplot_1.1.3           
#> [97] KEGGREST_1.44.0