Step 1: Load the example methylation dataset

The example data are distributed with myoRDNA in the inst/extdata directory.

We use system.file() to retrieve the location of the dataset within the installed package. This approach ensures that the example works regardless of where the package is installed.

library(myoRDNA)
library(methylKit)
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#> 
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#> 
#>     data
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#>     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#>     rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#>     tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo

file_path <- system.file(
  "extdata",
  "filtered_meth_raw_chr1_chrRDNAm_small.rds",
  package = "myoRDNA",
  mustWork = TRUE
)

loaded_data <- readRDS(file_path)

The RDS file contains previously processed methylation data from chromosome 1 and the rDNA chromosome (chrRDNAm). The underlying data retain information about methylation coverage, methylated read counts, unmethylated read counts and genomic positions.

As the data have already been processed and saved, users do not need to repeat the original sequencing alignment and methylation extraction procedures to follow this tutorial.

Step 2: Define experimental conditions

The dataset contains six OLD samples and four YOUNG samples. We define a binary treatment vector to distinguish these experimental groups:

• 0: OLD
• 1: YOUNG

The order of the treatment assignments must correspond exactly to the order of the samples in loaded_data.

treatment <- c(
  0, 0, 0, 0, 0, 0,
  1, 1, 1, 1
)

num_samples <- length(treatment)

This is a necessary component required for interacting with methylKit and this can be simply carried over for use in myoRDNA. The number of samples needs to be a numeric value and should be either the number of samples in the whole dataset or across two different conditions. We also define a named list containing the sequencing accession numbers associated with each condition.

samples <- list(
  'OLD' = c(
    "SRR28431752", "SRR28431753", "SRR28431754", "SRR28431755", 
    "SRR28431756", "SRR28431757"
  ),
  'YOUNG' = c(
    "SRR28431748", "SRR28431749", "SRR28431750", "SRR28431751"
  )
)

This is an important structure and needs to be of the same format as shown above for tracking samples when comparing rDNA metrics (such as copy number and aggregate methylation proportion) in the later steps. Please ensure that the sample ids match the sample ids for the corresponding data in methylKit, otherwise output will be NULL.

Step 3: Construct a methylRawList object

This step is only necessary for following along and using the rds data from the myoRDNA, if the user is using these features directly interfacing with the methylKit object than this isn’t needed.

We construct a methylRawList object by combining the individual sample objects with their treatment assignments.

stopifnot(length(loaded_data) == length(treatment))

loaded_data <- methods::new(
  "methylRawList",
  loaded_data,
  treatment = treatment
)

A methylRawList stores the methylation measurements from multiple samples alongside their experimental group assignments. This structure allows methylation data to be processed consistently across samples and enables subsequent comparisons between experimental conditions.

Step 4: Extract rDNA methylation data

The example dataset includes methylation measurements from chromosome 1 and the rDNA chromosome.

We use get_rDNA_data() to extract only the measurements associated with the rDNA chromosome.

meth_rDNA <- get_rDNA_data(
  loaded_data,
  num_samples,
  rDNA_chr_name = "chrRDNAm"
)

The rDNA_chr_name argument specifies the chromosome identifier used in the adapted reference genome.

In this example, the mouse rDNA chromosome is named chrRDNAm. This identifier must match the chromosome name in the original methylation data.

The resulting meth_rDNA object contains the rDNA-associated methylation measurements for all ten samples, preserving the information required for downstream analysis. If the name of the chromosome associated with rDNA is not known (in this case chrRDNAm) the user can make use of the identify unique chromosomes.

unique_names <- get_unique_chr_names(
    loaded_data,
    num_samples
)

This returns an array with the unique chromosome names that are included, additionally it performs a check to verify that all samples have consistent chromosome names.

Step 5: Adjust rDNA genomic coordinates

Depending on how the adapted reference genome was constructed, rDNA genomic coordinates may require adjustment to correspond to the reference sequence used for downstream analysis.

The shifting() function applies a common coordinate transformation to the selected rDNA data.

First, we specify the coordinate cutoff:

cutoff <- 35000

We then apply the coordinate shift:

meth_rDNA_shifted <- shifting(
  meth_rdna_raw = meth_rDNA,
  cutoff = cutoff,
  treatment = treatment
)

The function identifies the maximum genomic start position across the input samples and calculates the shift as

\[ s = x_{\max} - c, \] where \(x_{\max}\) denotes the maximum genomic start position and \(c\) is the user-defined cutoff. The transformed genomic coordinates are

\[ x_{\mathrm{shifted}} = x_{\mathrm{original}} - s. \] The same transformation is applied to the start and end positions of every retained methylation record. Records falling outside the permitted coordinate range are removed, and the transformed data are returned as a methylRawList containing the original treatment assignments.

The cutoff of 35,000 is specific to the adapted reference used in this example. Users analysing data generated with a different reference should select an appropriate cutoff based on their reference sequence and coordinate convention. The shifting is specific to the human and mouse genome, not enough is known about rDNA contig in the rat for this shifting.

The need for shifting is obvious when plotting the numCs against genomic position. The activity will start just after genomic position at the point \(s\). The shifting is vital for regional analysis of the 18S for differential methylation. However, functions surrounding the rDNA copy number will produce an interpretative result without shifting.

Absolute rDNA copy number

The absolute rDNA copy-number estimator uses the ratio of rDNA sequencing coverage to reference genomic coverage.

For a diploid genome, the estimator is defined as

\[ \boxed{ CN_{\mathrm{absolute}} = 2 \times \frac{\overline{C}_{\mathrm{rDNA}}} {\overline{C}_{\mathrm{genome}}} } \]

where:

• \(\overline{C}_{\mathrm{rDNA}}\) is the coverage estimate for the rDNA reference sequence.
• \(\overline{C}_{\mathrm{genome}}\) is the corresponding genomic reference coverage.
• The factor of two accounts for the two copies of each autosome in a diploid genome.

For a haploid genome, the estimator becomes

\[ CN_{\mathrm{absolute}} = \frac{\overline{C}_{\mathrm{rDNA}}} {\overline{C}_{\mathrm{genome}}}. \]

More generally, for a genome with reference ploidy \(P\),

\[ CN_{\mathrm{absolute}} = P \times \frac{\overline{C}_{\mathrm{rDNA}}} {\overline{C}_{\mathrm{genome}}}. \]

The current implementation supports diploid (\(P=2\)) and non-diploid (\(P=1\)) calculations.

absolute_CN <- absolute_rDNA_CN(
  Meth_object = loaded_data,
  num_samples = num_samples,
  autosomes = TRUE,
  rDNA_chr = "chrRDNAm",
  diploid = TRUE
)

absolute_CN
#>      Sample_ID Genome_Depth rDNA_Depth Absolute_rDNA_CN
#> 1  SRR28431752       3.3421   27.04982         16.18732
#> 2  SRR28431753       2.8245   28.40987         20.11675
#> 3  SRR28431754       5.8205   68.65034         23.58916
#> 4  SRR28431755       4.4878   57.31855         25.54416
#> 5  SRR28431756       3.6623   26.87927         14.67890
#> 6  SRR28431757       3.8956   38.12884         19.57534
#> 7  SRR28431748       3.7810   38.35561         20.28860
#> 8  SRR28431749       3.3540   56.07577         33.43815
#> 9  SRR28431750       4.0604   62.75377         30.91014
#> 10 SRR28431751       3.9073   65.39894         33.47526

Relative rDNA copy number

Relative rDNA copy number provides a normalised measure of rDNA abundance that can be used to compare samples processed using compatible sequencing protocols. Unlike absolute copy-number estimation, relative estimation is supported for both WGBS and RRBS data in myoRDNA.

The relative estimator is calculated using the ratio of summed rDNA coverage to summed genomic reference coverage:

\[ \boxed{ CN_{\mathrm{relative}} = \frac{\displaystyle\sum_{i=1}^{N_r} C_{\mathrm{rDNA},i}} {\displaystyle\sum_{j=1}^{N_g} C_{\mathrm{genome},j}} } \]

where:

• \(C_{\mathrm{rDNA},i}\) denotes the coverage of the \(i\)-th retained rDNA position.
• \(C_{\mathrm{genome},j}\) denotes the coverage of the \(j\)-th retained genomic reference position.
• \(N_r\) and \(N_g\) are the numbers of retained rDNA and genomic reference positions, respectively.

The numerator and denominator therefore represent the total coverage across the rDNA and genomic reference positions used by the estimator. Unlike the absolute estimator, the relative estimator does not apply a diploidy correction in the current implementation.

relative_CN <- relative_rDNA_CN(
  Meth_object = loaded_data,
  num_samples = num_samples,
  autosomes = TRUE,
  rDNA_chr = "chrRDNAm",
  diploid = TRUE
)

relative_CN
#>      Sample_ID Total_Number_of_Reads Number_of_Reads_rDNA Relative_rDNA_CN
#> 1  SRR28431752                 33421                91753         2.745370
#> 2  SRR28431753                 28245                97247         3.442981
#> 3  SRR28431754                 58205               250917         4.310918
#> 4  SRR28431755                 44878               205487         4.578791
#> 5  SRR28431756                 36623                88836         2.425689
#> 6  SRR28431757                 38956               134061         3.441344
#> 7  SRR28431748                 37810               133746         3.537318
#> 8  SRR28431749                 33540               202041         6.023882
#> 9  SRR28431750                 40604               224533         5.529825
#> 10 SRR28431751                 39073               234586         6.003788

Considerations for copy-number estimation

Several factors should be considered when interpreting rDNA copy-number estimates. First, the repetitive nature of rDNA introduces potential alignment ambiguity. Reads derived from different rDNA repeat units will align to the same reference sequence.

Second, sufficient sequencing coverage is required to obtain stable estimates. In particular, absolute copy-number estimation depends on having a suitable genomic reference against which rDNA coverage can be normalised. This can become particularly difficult with PBAT which typically provides lower alignment rates.

Finally, the adapted reference genome and chromosome identifiers must be consistent with those used during sequencing alignment and methylation extraction. For reproducible analyses, users should document the reference genome, coverage-filtering thresholds, sequencing methodology and genomic regions used for normalisation.

Step 1: Partition rDNA into genomic regions

The sort_by_rDNA_region() function separates methylation measurements into their corresponding rDNA regions using predefined genomic coordinates.

The package supports the following reference assemblies:

For the example mouse WGBS dataset, we use the mm39 reference assembly.

region_res <- sort_by_rDNA_region(
  meth_rDNA,
  num_samples,
  species = "mm39",
  min_cov = 10
)
#> Species selected: mm39

The species argument specifies the reference assembly used to identify the rDNA regions. Users analysing rat or human data should specify rno7 or hg38, respectively.

The min_cov argument defines the minimum sequencing coverage required for a CpG position to be retained. In this example, CpG positions with fewer than ten sequencing reads are excluded.

The resulting region_res object organises methylation measurements by sample and rDNA region, allowing downstream calculations to be performed independently for each region. ### rDNA regions

Depending on the reference assembly, rDNA measurements can be assigned to the following annotated regions:

• 18S: Encodes the 18S ribosomal RNA component of the small ribosomal subunit.
• 5.8S: Encodes 5.8S ribosomal RNA.
• 28S: Encodes the 28S ribosomal RNA component of the large ribosomal subunit.
• 5’ ETS: The external transcribed spacer located upstream of the 18S region.
• ITS1: The internal transcribed spacer between the 18S and 5.8S regions.
• ITS2: The internal transcribed spacer between the 5.8S and 28S regions.
• 3’ ETS: The external transcribed spacer downstream of the 28S region.

These regions can be analysed individually to investigate spatial variation in rDNA methylation. It is important to note that the human and mouse have all regions in the list above. However, the rate is limited to the 18S, 28S, and 5.8S.

Step 2: Calculate aggregate methylation metrics

Following regional partitioning, the calculate_aggregate_meth() function can be used to summarise methylation measurements within each rDNA region.

For this example, we calculate aggregate methylation for the 18S region.

region <- "18S"

aggregate_18S <- calculate_aggregate_meth(
  region_res,
  region,
  conditions = samples,
  med = FALSE
)

The region argument specifies the rDNA region to analyse, while conditions associates samples with their experimental groups.

The med argument controls the choice of summary statistic. In the mean-based calculation shown here, it is set to FALSE.

The resulting object contains aggregate methylation measurements for the selected region, organised according to the supplied experimental conditions.

Mathematical formulation of methylation proportions

Let \(N_{C,i}\) and \(N_{T,i}\) denote the methylated and unmethylated read counts, respectively, at CpG position \(i\) within a selected rDNA region.

The sequencing coverage at that position is

\[ C_i = N_{C,i} + N_{T,i}. \]

The methylation proportion at each retained CpG position is defined as

\[ p_i = \frac{N_{C,i}} {N_{C,i}+N_{T,i}}. \]

Only CpG positions satisfying the specified minimum coverage threshold are included in the aggregation.

Statistical methodology

Step 1: Site-level differential methylation

Site-level differential methylation identifies individual CpG positions that exhibit differences in methylation between experimental conditions.

For each genomic position, myoRDNA fits a separate regression model using the methylated and unmethylated read counts across samples.

Here, we analyse the 18S region of the previously prepared mouse WGBS dataset.

site_results <- diff_meth(
  meth_data = meth_rDNA_shifted,
  samples = samples,
  species = "mm39",
  num_samples = num_samples,
  min_cov = 10,
  region = "18S",
  model_type = "binomial",
  p_adjust_method = "BH",
  site = TRUE
)
#> Species selected: mm39

Setting site = TRUE performs the differential methylation analysis independently at each retained CpG position.

The min_cov = 10 argument excludes measurements with fewer than ten sequencing reads.

To use either of the alternative regression models, change the model_type argument: binomial, quasi-binomial, and beta-binomial.

Beta-binomial regression can be computationally more demanding, particularly when thousands of CpG positions are analysed independently. As a result of the algorithm requiring the use of optimisation to estimate parameters in the model.

Step 2: Region-level differential methylation

Region-level differential methylation investigates differences across an entire annotated rDNA region, rather than fitting independent models at each CpG position.

For a selected region, the methylated and unmethylated read counts can be aggregated within each sample:

\[ M_{i,R} = \sum_{j \in R} N_{C,ij}, \]

\[ U_{i,R} = \sum_{j \in R} N_{T,ij}, \]

where \(R\) denotes the selected region. The aggregated counts are then used to estimate the association between methylation probability and experimental condition.

The following example performs region-level differential methylation for 18S.

region_results <- diff_meth(
  meth_data = meth_rDNA_shifted,
  samples = samples,
  species = "mm39",
  num_samples = num_samples,
  min_cov = 10,
  region = "18S",
  model_type = "binomial",
  p_adjust_method = "BH",
  site = FALSE
)
#> Species selected: mm39

Setting site = FALSE performs the analysis at the regional level. As with site-level analysis, the regression method can be changed using model_type.

Region-level aggregation provides a summary of methylation differences across the selected region. However, this approach can obscure opposing methylation changes at individual CpG positions and should therefore be considered complementary to site-level analysis.

Step 3: Interpreting differential methylation

For the binary treatment coding between CONTROL and CONDITION 1 used in this example, the estimated methylation probabilities are

\[ \widehat{p}_{\mathrm{CONTROL}} = \frac{1}{1+\exp(-\widehat{\beta}_0)} \] and

\[ \widehat{p}_{\mathrm{CONDITION \ 1}} = \frac{1}{ 1+\exp[-(\widehat{\beta}_0+\widehat{\beta}_1)] }. \] An odds ratio greater than one indicates higher fitted methylation odds in CONDITION 1 samples, whereas an odds ratio below one indicates lower fitted methylation odds.

The estimated methylation difference, expressed in percentage points, is

\[ \boxed{ \Delta_{\mathrm{meth}} = 100 \left( \widehat{p}_{\mathrm{CONDITION \ 1}} - \widehat{p}_{\mathrm{CONTROL}} \right) } \] A positive value therefore indicates higher estimated methylation in CONDITION 1 samples relative to OLD samples. Statistical evidence for a treatment effect is assessed through the corresponding regression coefficient and its associated hypothesis test.

Step 4: Multiple-testing correction

Site-level differential methylation involves testing multiple CpG positions, increasing the probability of obtaining false-positive results.

To account for multiple comparisons, myoRDNA supports p-value adjustment using the Benjamini–Hochberg (BH) as a default procedure. However, if the user wishes to use a less or more conservative method this can be adjusted through the p_adjust_method to none for no adjustment or bonferroni or holm to name a few.

The BH method controls the false discovery rate (FDR), defined as the expected proportion of false discoveries among rejected hypotheses, under its applicable dependence assumptions.

The adjusted p-values are then returned to their original genomic positions.

For example, an FDR threshold of 0.05 can be used to identify CpG positions with BH-adjusted p-values below 0.05.

The scope of the correction depends on the set of tests included. When different rDNA regions are analysed separately, region-specific adjustments do not automatically control the FDR across all regions.

Beta-binomial optimisation and model diagnostics

Beta-binomial regression estimates the regression coefficients and dispersion parameter by maximising the likelihood of the observed methylation counts.

Because the beta-binomial likelihood can be difficult to optimise, particularly with limited biological replication or extreme methylation proportions, numerical convergence should be assessed before interpreting the fitted coefficients and p-values.

# Inspect fitted model
summary(fit)

# Check convergence code
fit$fit$convergence

# Check whether the Hessian is positive definite
fit$sdr$pdHess

# Run additional diagnostics
glmmTMB::diagnose(fit)

Session information

sessionInfo()
#> 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] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] methylKit_1.39.0     GenomicRanges_1.65.4 Seqinfo_1.3.2       
#> [4] IRanges_2.47.5       S4Vectors_0.51.10    BiocGenerics_0.59.12
#> [7] generics_0.1.4       myoRDNA_0.99.0       rmarkdown_2.32      
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1         dplyr_1.2.1              farver_2.1.2            
#>  [4] R.utils_2.13.0           Biostrings_2.81.9        S7_0.2.2                
#>  [7] bitops_1.1-0             fastmap_1.2.0            RCurl_1.98-1.20         
#> [10] GenomicAlignments_1.49.2 XML_3.99-0.25            digest_0.6.39           
#> [13] lifecycle_1.0.5          statmod_1.5.2            magrittr_2.0.5          
#> [16] compiler_4.6.1           rlang_1.3.0              sass_0.4.10             
#> [19] tools_4.6.1              yaml_2.3.12              data.table_1.18.6.1     
#> [22] rtracklayer_1.73.0       knitr_1.52               mclust_6.1.3            
#> [25] curl_8.0.0               plyr_1.8.9               RColorBrewer_1.1-3      
#> [28] BiocParallel_1.47.0      numDeriv_2016.8-1.1      sys_3.4.3               
#> [31] R.oo_1.27.1              grid_4.6.1               ggplot2_4.0.3           
#> [34] scales_1.4.0             gtools_3.9.5             MASS_7.3-66             
#> [37] bbmle_1.0.26             cli_3.6.6                mvtnorm_1.4-2           
#> [40] crayon_1.5.3             otel_0.2.0               httr_1.4.9              
#> [43] reshape2_1.4.5           rjson_0.2.23             BiocBaseUtils_1.15.1    
#> [46] bdsmatrix_1.3-7          qvalue_2.45.0            cachem_1.1.0            
#> [49] stringr_1.6.0            splines_4.6.1            parallel_4.6.1          
#> [52] XVector_0.53.0           restfulr_0.0.17          vctrs_0.7.3             
#> [55] Matrix_1.7-6             jsonlite_2.0.0           fastseg_1.59.0          
#> [58] maketools_1.3.2          limma_3.99.0             jquerylib_0.1.4         
#> [61] glue_1.8.1               emdbook_1.3.14           codetools_0.2-20        
#> [64] stringi_1.8.9            gtable_0.3.6             BiocIO_1.23.3           
#> [67] tibble_3.3.1             pillar_1.11.1            htmltools_0.5.9         
#> [70] R6_2.6.1                 evaluate_1.0.5           lattice_0.23-1          
#> [73] Biobase_2.73.2           R.methodsS3_1.8.2        Rsamtools_2.29.0        
#> [76] cigarillo_1.3.1          bslib_0.12.0             Rcpp_1.1.2              
#> [79] coda_0.19-4.1            xfun_0.61                buildtools_1.0.0        
#> [82] pkgconfig_2.0.3