Relative Enrichment

The first metric that we use is a relative enrichment metric, known as relH. This was first introduced in the MEDIPs package, via the MEDIPS.CpGenrich() function, and is a metric to express the proportion of CG dinucleotides, compared to the reference genome. The relative enrichment relH is calculated in the following way: * For each mapped DNA fragment, retrieve the end-to-end DNA spanned on the reference sequence. This is used instead of the reads themselves, as a fragment will often be longer than the read(s) sequenced at its ends, and may therefore contain CGs not covered by the reads. * In these sequences, count how many CGs are present in total, and divide by the total number of base pairs covered. * Divide by the same ratio, calculated for the entire genome.

This is therefore a measure of how many more CGs are present in the fragments covered within the bam file, compared to the reference genome.

Within mesa, this metric is calculated through makeQset, as part of the information added to the libraries table on the qseaSet. For the MBD-Seq data that we have studied, we apply a minimum relH of 2.5, with 3 as a higher cutoff. This value is typically around 4-5 in many of our studies.

We can use the addLibraryInformation function to add the information in the libraries table directly onto the sampleTable:

exampleTumourNormal %>% 
    addLibraryInformation() %>% 
    getSampleTable() %>% 
    dplyr::select(sample_name, relH)
#>          sample_name     relH
#> Colon1_N    Colon1_N 5.567633
#> Colon1_T    Colon1_T 4.921130
#> Colon2_N    Colon2_N 5.399262
#> Colon2_T    Colon2_T 5.686246
#> Lung1_N      Lung1_N 4.463002
#> Lung1_T      Lung1_T 4.406991
#> Lung2_N      Lung2_N 4.602440
#> Lung2_T      Lung2_T 4.753918
#> Lung3_N      Lung3_N 5.111422
#> Lung3_T      Lung3_T 5.078306

When a non-enriched (“Input”) sample is used to calculate copy number variation, an equivalent calculation is also performed on these non-enriched fragments. This is also stored in the libraries table, in columns starting with input:

exampleTumourNormal %>% 
    addLibraryInformation() %>% 
    getSampleTable() %>% 
    dplyr::select(sample_name, input_relH)
#>          sample_name input_relH
#> Colon1_N    Colon1_N   1.506917
#> Colon1_T    Colon1_T   1.446767
#> Colon2_N    Colon2_N   1.608130
#> Colon2_T    Colon2_T   1.474401
#> Lung1_N      Lung1_N   1.437550
#> Lung1_T      Lung1_T   1.481526
#> Lung2_N      Lung2_N   1.418541
#> Lung2_T      Lung2_T   1.400016
#> Lung3_N      Lung3_N   1.644814
#> Lung3_T      Lung3_T   1.623274

We typically see values of the input_relH between 1.3-1.6, suggesting that there are more CGs seen in the sequenced fragments than would be expected from the genomic sequence. This is presumably due to the evolutionary suppression of CG dinucleotides outside of specific contexts, meaning they are less common in those areas of the genome that are repetitive and hence less easy to map. This is likely to be dependent on the particular sequencing methods and the sample type under consideration. We therefore need to see an increase between the enriched and non-enriched data:

exampleTumourNormal %>% 
    addLibraryInformation() %>% 
    getSampleTable() %>% 
    dplyr::select(sample_name, relH_MeCap = relH, relH_Input = input_relH) %>% 
    pivot_longer(matches("relH"), 
                names_prefix = "relH", 
                names_to = "enrichment", 
                values_to = "relH") %>%
    ggplot(aes(x = sample_name, y = relH, col = enrichment)) +
    geom_point() +
    theme_bw()

There is a second relative enrichment metric defined by the MEDIPs package, known as GoGe, that is also calculated globally in a similar way but considers the ratio between the number of CGs present in the fragments compared to the number of Cs and Gs independently. This is calculated and stored in the libraries information alongside the relH metric.

Fragments without pattern

The libraries tab also calculates how many fragments do not contain a CG at all. This can be measure of background noise.

Hyper-stably methylated regions

The other metric that we typically use for quality control is to consider windows that are known to be commonly methylated across a range of tissues and diseases. For instance, for the human genome, [Edgar et al (2014)]((https://pubmed.ncbi.nlm.nih.gov/25493099/) identified probes from the Illumina 450k methylation array that were consistently methylated in over 1000 samples. We therefore take the windows covered by these probes, and calculate beta values for them, followed by determining how many of these windows have a calculated beta value above a threshold (e.g. 0.8). As noted above, beta values will not be calculated if a sample does not have sufficient coverage and/or enrichment, meaning this metric incorporates both the coverage depth as well as how well the enrichment has worked. For our samples within the National Biomarker Centre, we currently require at least 40% of the hyperstable windows with a CpG_density above 5 to have a beta value above 0.8. Samples with lower enrichment scores (as described by the relH score) require a greater coverage depth to reach this threshold.

For hg38, these hyper-stable methylated probes are incorporated in a data object hg38UltraStableProbes, and a function addHyperStableFraction is provided to calculate this and add as a column, hyperStableEdgar, to the sampleTable for filtering. Note that this will not work with the provided example data, as there is no overlap with the windows in the example dataset (i.e. filterByOverlaps(exampleTumourNormal,hg38UltraStableProbes) has no windows).

Summary

We find that these three metrics, valid_fragments, relH and hyperStableEdgar are the most important considerations. For 300bp windows on the human genome, with T7-MBD-seq, we require a relH of at least 2.5 and hyperStableEdgar value of at least 0.4 (40%). We may also apply a valid_fragments cutoff of 3000000, although most samples which pass the other two metrics will achieve this anyway.

Session info

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.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.32.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=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: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] stringr_1.6.0  tidyr_1.3.2    ggplot2_4.0.3  mesa_0.99.7    dplyr_1.2.1   
#> [6] qsea_1.39.0    rmarkdown_2.32
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6             rjson_0.2.23             xfun_0.61               
#>  [4] bslib_0.12.0             Biobase_2.73.2           lattice_0.23-1          
#>  [7] vctrs_0.7.3              tools_4.6.1              bitops_1.1-0            
#> [10] generics_0.1.4           stats4_4.6.1             curl_8.0.0              
#> [13] parallel_4.6.1           tibble_3.3.1             pkgconfig_2.0.3         
#> [16] BiocBaseUtils_1.15.1     Matrix_1.7-6             data.table_1.18.6.1     
#> [19] BSgenome_1.81.1          RColorBrewer_1.1-3       S7_0.2.2                
#> [22] S4Vectors_0.51.10        cigarillo_1.3.1          lifecycle_1.0.5         
#> [25] farver_2.1.2             compiler_4.6.1           Rsamtools_2.29.0        
#> [28] textshaping_1.0.5        Biostrings_2.81.9        statmod_1.5.2           
#> [31] Seqinfo_1.3.2            codetools_0.2-20         GenomeInfoDb_1.49.1     
#> [34] htmltools_0.5.9          sys_3.4.3                buildtools_1.0.0        
#> [37] sass_0.4.10              RCurl_1.98-1.20          yaml_2.3.12             
#> [40] pillar_1.11.1            crayon_1.5.3             jquerylib_0.1.4         
#> [43] BiocParallel_1.47.0      cachem_1.1.0             limma_3.99.0            
#> [46] gtools_3.9.5             tidyselect_1.2.1         digest_0.6.39           
#> [49] stringi_1.8.9            purrr_1.2.2              restfulr_0.0.17         
#> [52] labeling_0.4.3           maketools_1.3.2          fastmap_1.2.0           
#> [55] grid_4.6.1               cli_3.6.6                magrittr_2.0.5          
#> [58] utf8_1.2.6               XML_3.99-0.25            withr_3.0.3             
#> [61] UCSC.utils_1.9.0         scales_1.4.0             XVector_0.53.0          
#> [64] httr_1.4.9               matrixStats_1.5.0        otel_0.2.0              
#> [67] ragg_1.5.2               zoo_1.9-1                memoise_2.0.1           
#> [70] evaluate_1.0.5           knitr_1.52               GenomicRanges_1.65.4    
#> [73] IRanges_2.47.5           BiocIO_1.23.3            rtracklayer_1.73.0      
#> [76] rlang_1.3.0              glue_1.8.1               HMMcopy_1.55.1          
#> [79] BiocGenerics_0.59.12     jsonlite_2.0.0           R6_2.6.1                
#> [82] GenomicAlignments_1.49.2 systemfonts_1.3.2