# Working with RangedSummarizedExperiments Big DNA Methylation Data

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

## 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")
    #> see ?TumourMethData and browseVignettes('TumourMethData') for documentation
    #> loading from cache

## 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

    # 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

## 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 `plotRegionValues()` 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 = plotRegionValues(gstp1_cpg_methylation, sample_name = "DTB_003")
    print(gstp1_methylation_plot)
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

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"))
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

    # 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"))
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

    # 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)
    #> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

## 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"]]
    #> downloading 1 resources
    #> retrieving 1 resource
    #> 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)

## 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"]]
    #> downloading 1 resources
    #> retrieving 1 resource
    #> 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

## SessionInfo

    sessionInfo()
    #> R version 4.5.2 (2025-10-31)
    #> Platform: x86_64-pc-linux-gnu
    #> Running under: Ubuntu 24.04.3 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.UCSC.hg19_1.4.3 BSgenome_1.78.0
    #>  [3] BiocIO_1.20.0                     Biostrings_2.78.0
    #>  [5] XVector_0.50.0                    rtracklayer_1.70.1
    #>  [7] AnnotationHub_4.0.0               BiocFileCache_3.0.0
    #>  [9] dbplyr_2.5.1                      rhdf5_2.54.1
    #> [11] TumourMethData_1.8.0              DESeq2_1.50.2
    #> [13] methodical_1.6.0                  SummarizedExperiment_1.40.0
    #> [15] Biobase_2.70.0                    MatrixGenerics_1.22.0
    #> [17] matrixStats_1.5.0                 ggplot2_4.0.1
    #> [19] GenomicRanges_1.62.1              Seqinfo_1.0.0
    #> [21] IRanges_2.44.0                    S4Vectors_0.48.0
    #> [23] BiocGenerics_0.56.0               generics_0.1.4
    #> [25] BiocStyle_2.38.0
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] DBI_1.2.3                bitops_1.0-9             httr2_1.2.2
    #>   [4] rlang_1.1.7              magrittr_2.0.4           otel_0.2.0
    #>   [7] compiler_4.5.2           RSQLite_2.4.5            mgcv_1.9-4
    #>  [10] GenomicFeatures_1.62.0   reshape2_1.4.5           png_0.1-8
    #>  [13] vctrs_0.7.1              stringr_1.6.0            pkgconfig_2.0.3
    #>  [16] crayon_1.5.3             fastmap_1.2.0            labeling_0.4.3
    #>  [19] Rsamtools_2.26.0         rmarkdown_2.30           tzdb_0.5.0
    #>  [22] UCSC.utils_1.6.1         purrr_1.2.1              bit_4.6.0
    #>  [25] xfun_0.56                cachem_1.1.0             cigarillo_1.0.0
    #>  [28] GenomeInfoDb_1.46.2      jsonlite_2.0.0           blob_1.3.0
    #>  [31] rhdf5filters_1.22.0      DelayedArray_0.36.0      Rhdf5lib_1.32.0
    #>  [34] BiocParallel_1.44.0      parallel_4.5.2           R6_2.6.1
    #>  [37] stringi_1.8.7            bslib_0.9.0              RColorBrewer_1.1-3
    #>  [40] jquerylib_0.1.4          iterators_1.0.14         Rcpp_1.1.1
    #>  [43] knitr_1.51               R.utils_2.13.0           readr_2.1.6
    #>  [46] splines_4.5.2            Matrix_1.7-4             tidyselect_1.2.1
    #>  [49] abind_1.4-8              yaml_2.3.12              codetools_0.2-20
    #>  [52] curl_7.0.0               plyr_1.8.9               regioneR_1.42.0
    #>  [55] lattice_0.22-7           tibble_3.3.1             withr_3.0.2
    #>  [58] KEGGREST_1.50.0          S7_0.2.1                 evaluate_1.0.5
    #>  [61] ExperimentHub_3.0.0      pillar_1.11.1            BiocManager_1.30.27
    #>  [64] filelock_1.0.3           foreach_1.5.2            vroom_1.6.7
    #>  [67] RCurl_1.98-1.17          hms_1.1.4                BiocVersion_3.22.0
    #>  [70] scales_1.4.0             glue_1.8.0               maketools_1.3.2
    #>  [73] tools_4.5.2              sys_3.4.3                data.table_1.18.0
    #>  [76] locfit_1.5-9.12          GenomicAlignments_1.46.0 buildtools_1.0.0
    #>  [79] XML_3.99-0.20            cowplot_1.2.0            grid_4.5.2
    #>  [82] AnnotationDbi_1.72.0     nlme_3.1-168             HDF5Array_1.38.0
    #>  [85] restfulr_0.0.16          annotatr_1.36.0          cli_3.6.5
    #>  [88] rappdirs_0.3.4           S4Arrays_1.10.1          dplyr_1.1.4
    #>  [91] gtable_0.3.6             R.methodsS3_1.8.2        sass_0.4.10
    #>  [94] digest_0.6.39            SparseArray_1.10.8       rjson_0.2.23
    #>  [97] farver_2.1.2             memoise_2.0.1            htmltools_0.5.9
    #> [100] R.oo_1.27.1              lifecycle_1.0.5          h5mread_1.2.1
    #> [103] httr_1.4.7               bit64_4.6.0-1