Licensing

The methylInheritance package and the underlying methylInheritance code are distributed under the Artistic license 2.0. You are free to use and redistribute this software.

Citing

If you use this package for a publication, we would ask you to cite the following:

Pascal Belleau, Astrid Deschênes, Marie-Pier Scott-Boyer, Romain Lambrot, Mathieu Dalvai, Sarah Kimmins, Janice Bailey, Arnaud Droit; Inferring and modeling inheritance of differentially methylated changes across multiple generations, Nucleic Acids Research, Volume 46, Issue 14, 21 August 2018, Pages e85. DOI: https://doi.org/10.1093/nar/gky362

Introduction

DNA methylation plays an important role in the biology of tissue development and diseases. High-throughput sequencing techniques enable genome-wide detection of differentially methylated elements (DME), commonly sites (DMS) or regions (DMR). The analysis of treatment effects on DNA methylation, from one generation to the next (inter-generational) and across generations that were not exposed to the initial environment (trans-generational) represent complex designs. Due to software design, the detection of DME is usually made on each generation separately. However, the common DME between generations due to randomness is not negligible when the number of DME detected in each generation is high. To judge the effect on DME that is inherited from a treatment in previous generation, the observed number of conserved DME must be compared to the randomly expected number.

We present a permutation analysis, based on Monte Carlo sampling, aim to infer a relation between the number of conserved DME from one generation to the next to the inheritance effect of treatment and to dismiss stochastic effect. It is used as a robust alternative to inference based on parametric assumptions.

The methylInheritance package can perform a permutation analysis on both differentially methylated sites (DMS) and differentially methylated tiles (DMT) using the methylKit package.

Loading methylInheritance package

As with any R package, the methylInheritance package should first be loaded with the following command:

library(methylInheritance)

Description of the permutation analysis

The permutation analysis is a statistical significance test in which the distribution of the test statistic under the null hypothesis is obtained by calculating the values of the test statistic under rearrangement of the labels on the observed data points. The rearrangement of the labels is done through repeated random sampling (Legendre and Legendre 1998, 142–57).

Null Hypothesis: The number of conserved DME correspond to a number that can be obtained through a randomness analysis.

Alternative Hypothesis: The number of conserved DME do not correspond to a number that can be obtained through a randomness analysis.

A typical methylInheritance analysis consists of the following steps:

1. Process to a differentially methylation analysis on each generation separately using real dataset with the methylKit package.
2. Calculate the number of conserved differentially methylated elements between two consecutive generations (F1 and F2, F2 and F3, etc..). The number of conserved differentially methylated elements is also calculated for three or more consecutive generations, always starting with the first generation (F1 and F2 and F3, F1 and F2 and F3 and F4, etc..). Those results are considered the reference values.
3. Fix a threshold (conventionally 0.05) that is used as a cutoff between the null and alternative hypothesis.
4. Process to a differential methylation analysis on each shuffled dataset. Each generation is analysed separately using the methylKit package.
5. Calculate the significant level for each consecutive subset of generations. The significant level is defined as the percentage of results equal or higher than the reference values. The reference values are added to the analysis so that it becomes impossible for the test to conclude that no value is as extreme as, or more extreme than the reference values.

All those steps have been encoded in the methylInheritance package.

Case study

## Load dataset containing information over three generations
data(demoForTransgenerationalAnalysis)

## The length of the dataset corresponds to the number of generation
## The generations are stored in order (first entry = first generation, 
## second entry = second generation, etc..)
length(demoForTransgenerationalAnalysis)
## [1] 3


## All samples related to one generation are contained in a methylRawList 
## object.
## The methylRawList object contains two Slots:
## 1- treatment: a numeric vector denoting controls and cases.
## 2- .Data: a list of methylRaw objects. Each object stores the raw 
##           mehylation data of one sample.


## A section of the methylRaw object containing the information of the 
## first sample from the second generation 
head(demoForTransgenerationalAnalysis[[2]][[1]]) 
##   chr start  end strand coverage numCs numTs
## 1   S  1000 1000      +       88    87     1
## 2   S  1038 1038      +       77    76     1
## 3   S  1061 1061      +       84    84     0
## 4   S  1098 1098      +      100    46    54
## 5   S  1758 1758      +       93    93     0
## 6   S  1801 1801      +       77    76     1

## The treatment vector for each generation
## The number of treatments and controls is the same in each generation
## However, it could also be different.
## Beware that getTreatment() is a function from the methylKit package.
getTreatment(demoForTransgenerationalAnalysis[[1]])
##  [1] 0 0 0 0 0 0 1 1 1 1 1 1
getTreatment(demoForTransgenerationalAnalysis[[2]])
##  [1] 0 0 0 0 0 0 1 1 1 1 1 1
getTreatment(demoForTransgenerationalAnalysis[[3]])
##  [1] 0 0 0 0 0 0 1 1 1 1 1 1

## The observation analysis is only done on differentially methylated sites
runObservation(methylKitData = demoForTransgenerationalAnalysis, 
                        type = "sites",     # Only sites
                        outputDir = "demo_01",   # RDS result files are saved 
                                                 # in the directory
                        nbrCores = 1,       # Number of cores used 
                        minReads = 10,      # Minimum read coverage
                        minMethDiff = 10,   # Minimum difference in methylation 
                                            # to be considered DMS
                        qvalue = 0.01,
                        vSeed = 2101)       # Ensure reproducible results
## [1] 0

## The results can be retrived using loadAllRDSResults() method
observedResults <- loadAllRDSResults(
                    analysisResultsDir = "demo_01/",  # Directory containing
                                                      # the observation
                                                      # results
                    permutationResultsDir = NULL, 
                    doingSites = TRUE, 
                    doingTiles = FALSE)

observedResults
## Permutation Analysis
## 
## Number of Generations:  3 
## Number of Permutations:  0 
## 
## Observation Results: 
##        SOURCE ELEMENT ANALYSIS   TYPE RESULT
## 1 OBSERVATION   SITES       i2 HYPER1     16
## 2 OBSERVATION   SITES       i2 HYPER2     16
## 3 OBSERVATION   SITES       i2  HYPO1      2
## 4 OBSERVATION   SITES       i2  HYPO2      4
## 5 OBSERVATION   SITES     iAll  HYPER     14
## 6 OBSERVATION   SITES     iAll   HYPO      2

## The permutation analysis is only done on differentially methylated sites
runPermutation(methylKitData = demoForTransgenerationalAnalysis, # multi-generational dataset
                        type = "sites",     # Only sites
                        outputDir = "demo_02",   # RDS permutation files are 
                                                 # saved in the directory
                        runObservationAnalysis = FALSE,
                        nbrCores = 1,           # Number of cores used
                        nbrPermutations = 2,    # Should be much higher for a
                                                # real analysis
                        minReads = 10,          # Minimum read coverage
                        minMethDiff = 10,   # Minimum difference in methylation
                                            # to be considered DMS
                        qvalue = 0.01,
                        vSeed = 2101)         # Ensure reproducible results
## [1] 0

## The results can be retrived using loadAllRDSResults() method
permutationResults <- loadAllRDSResults(
                    analysisResultsDir = NULL, 
                    permutationResultsDir = "demo_02",   # Directory containing
                                                    # the permutation
                                                    # results
                    doingSites = TRUE, 
                    doingTiles = FALSE)

permutationResults
## Permutation Analysis
## 
## Number of Generations:  0 
## Number of Permutations:  2 
## 
## Observation Results: 
##  No observation result.

## Merge observation and permutation results
allResults <- mergePermutationAndObservation(permutationResults = 
                                                    permutationResults,
                                    observationResults = observedResults)
allResults
## Permutation Analysis
## 
## Number of Generations:  3 
## Number of Permutations:  2 
## 
## Observation Results: 
##        SOURCE ELEMENT ANALYSIS   TYPE RESULT
## 1 OBSERVATION   SITES       i2 HYPER1     16
## 2 OBSERVATION   SITES       i2 HYPER2     16
## 3 OBSERVATION   SITES       i2  HYPO1      2
## 4 OBSERVATION   SITES       i2  HYPO2      4
## 5 OBSERVATION   SITES     iAll  HYPER     14
## 6 OBSERVATION   SITES     iAll   HYPO      2

## Conserved differentially methylated sites between F1 and F2 generations
F1_and_F2_results <- extractInfo(allResults = allResults, type = "sites", 
                                    inter = "i2", position = 1)

head(F1_and_F2_results)
##    TYPE RESULT      SOURCE
## 1  HYPO      2 OBSERVATION
## 2 HYPER     16 OBSERVATION
## 3  HYPO     13 PERMUTATION
## 4 HYPER      2 PERMUTATION
## 5  HYPO      0 PERMUTATION
## 6 HYPER      1 PERMUTATION

## Differentially conserved sites between F1 and F2 generations
F1_and_F2 <- extractInfo(allResults = demoResults, type = "sites", 
                            inter = "i2", position = 1)
## Differentially conserved sites between F2 and F3 generations
F2_and_F3 <- extractInfo(allResults = demoResults, type = "sites", 
                            inter = "i2", position = 2)
## Differentially conserved sites between F1 and F2 and F3 generations
F2_and_F3 <- extractInfo(allResults = demoResults, type = "sites", 
                            inter = "iAll", position = 1)

## Show graph and significant level for differentially conserved sites 
## between F1 and F2 
output <- plotGraph(F1_and_F2)

Possibility to restart a permutation analysis

When a large number of permutations is processed, the time needed to process them all may be long (especially when the number of available CPU is limited). Furthermore, some permutations can fail due to parallelization problems.

The methylInheritance package offers the possibility to restart an analysis and run only missing permutation results. To take advantage of this option, the outputDir parameter must not be NULL so that permutation results are saved in RDS files. When the restartCalculation is set to TRUE, the method will load the permutation results present in RDS files (when available) and only rerun permutations that don’t have an associated RDS file.

## The permutation analysis is only done on differentially methylated tiles
## The "output" directory must be specified
## The "vSeed" must be specified to ensure reproducible results
## The "restartCalculation" is not important the first time the analysis is run
permutationResult <- runPermutation(
                        methylKitData = demoForTransgenerationalAnalysis, # multi-generational dataset
                        type = "tiles",     # Only tiles
                        outputDir = "test_restart",   # RDS files are created
                        runObservationAnalysis = TRUE,
                        nbrCores = 1,           # Number of cores used
                        nbrPermutations = 2,    # Should be much higher for a
                                                # real analysis
                        vSeed = 212201,     # Ensure reproducible results
                        restartCalculation = FALSE)

## Assume that the process was stopped before it has done all the permutations

## The process can be restarted
## All parameters must be identical to the first analysis except "restartCalculation"
## The "restartCalculation" must be set to TRUE
permutationResult <- runPermutation(
                        methylKitData = demoForTransgenerationalAnalysis, # multi-generational dataset
                        type = "tiles",     # Only tiles
                        outputDir = "test_restart",   # RDS files are created
                        runObservationAnalysis = TRUE,
                        nbrCores = 1,           # Number of cores used
                        nbrPermutations = 2,    # Should be much higher for a
                                                # real analysis
                        vSeed = 212201,     # Ensure reproducible results
                        restartCalculation = TRUE)         

Format multigenerational dataset into an input

The permutation analysis needs a list of methylRawList objects as input. A methylRawList is a list of methylRaw objects. The methylRawList and methylRaw objects are defined in the methylKit package.

To create a methylRawList, all samples (cases and controls) from the same generation must be first separately transformed into a methylRaw object. The S4 methylRaw class extends data.frame class and has been created to store raw methylation data. The raw methylation is essentially percent methylation values and read coverage values per base or region.

Excluding the data.frame section, the slots present in the methylRaw class are:

• sample.id: a string, the sample identification
• assembly: a string, the genomic assembly
• context: a string, the methylation context, as an exemple, CpG, CpH, etc…
• resolution: a string, the resolution of methylation information, mainly ‘base’ or ‘region’

## The list of methylRaw objects for all controls and cases related to F1
f1_list <- list()
f1_list[[1]] <- new("methylRaw", 
                    data.frame(chr = c("chr21", "chr21"), 
                    start = c(9764513, 9764542), 
                    end = c(9764513, 9764542), strand = c("+", "-"), 
                    coverage = c(100, 15), numCs = c(88, 2), 
                    numTs = c(100, 15) - c(88, 2)), 
                    sample.id = "F1_control_01", assembly = "hg19", 
                    context = "CpG", resolution = 'base')
f1_list[[2]] <- new("methylRaw", 
                    data.frame(chr = c("chr21", "chr21"), 
                    start = c(9764513, 9764522), 
                    end = c(9764513, 9764522), strand = c("-", "-"), 
                    coverage = c(38, 21), numCs = c(12, 2), 
                    numTs = c(38, 21) - c(12, 2)), 
                    sample.id = "F1_case_02", assembly = "hg19", 
                    context = "CpG", resolution = 'base')

## The list of methylRaw objects for all controls and cases related to F2
f2_list <- list()
f2_list[[1]] <- new("methylRaw", 
                    data.frame(chr = c("chr21", "chr21"), 
                    start = c(9764514, 9764522), 
                    end = c(9764514, 9764522), strand = c("+", "+"), 
                    coverage = c(40, 30), numCs = c(0, 2), 
                    numTs = c(40, 30) - c(0, 2)), 
                    sample.id = "F2_control_01", assembly = "hg19", 
                    context = "CpG", resolution = 'base')
f2_list[[2]] <- new("methylRaw", 
                    data.frame(chr = c("chr21", "chr21"), 
                    start = c(9764513, 9764533), 
                    end = c(9764513, 9764533), strand = c("+", "-"), 
                    coverage = c(33, 23), numCs = c(12, 1), 
                    numTs = c(33, 23) - c(12, 1)), 
                    sample.id = "F2_case_01", assembly = "hg19", 
                    context = "CpG", resolution = 'base')

## The list to use as input for methylInheritance 
final_list <- list()

## The methylRawList for F1 - the first generation is on the first position
final_list[[1]] <- new("methylRawList", f1_list, treatment = c(0,1))
## The methylRawList for F2 - the second generation is on the second position
final_list[[2]] <- new("methylRawList", f2_list, treatment = c(0,1))

## A list of methylRawList ready for methylInheritance
final_list
## [[1]]
## methylRawList object with 2 methylRaw objects
## 
## methylRaw object with 2 rows
## --------------
##     chr   start     end strand coverage numCs numTs
## 1 chr21 9764513 9764513      +      100    88    12
## 2 chr21 9764542 9764542      -       15     2    13
## --------------
## sample.id: F1_control_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## methylRaw object with 2 rows
## --------------
##     chr   start     end strand coverage numCs numTs
## 1 chr21 9764513 9764513      -       38    12    26
## 2 chr21 9764522 9764522      -       21     2    19
## --------------
## sample.id: F1_case_02 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## treatment: 0 1 
## 
## [[2]]
## methylRawList object with 2 methylRaw objects
## 
## methylRaw object with 2 rows
## --------------
##     chr   start     end strand coverage numCs numTs
## 1 chr21 9764514 9764514      +       40     0    40
## 2 chr21 9764522 9764522      +       30     2    28
## --------------
## sample.id: F2_control_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## methylRaw object with 2 rows
## --------------
##     chr   start     end strand coverage numCs numTs
## 1 chr21 9764513 9764513      +       33    12    21
## 2 chr21 9764533 9764533      -       23     1    22
## --------------
## sample.id: F2_case_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## treatment: 0 1

Another approach is to transform the files that contain the raw methylation information into a format that can be read by the methylKit methRead function. The methRead function implements methods that enable the creation of methylRawList objects.

Here is one valid file format among many (tab separated):

chrBase     chr     base    strand  coverage    freqC   freqT
1.176367    1       176367  R       29          100.00  0.00
1.176392    1       176392  R       58          100.00  0.00
1.176422    1       176422  R       29          3.45    96.55
1.176552    1       176552  R       58          96.55   3.45

library(methylKit)

## The methylRawList for F1
generation_01 <- methRead(location = list("demo/F1_control_01.txt", "demo/F1_case_01.txt"), 
                    sample.id = list("F1_control_01", "F1_case_01"), 
                    assembly = "hg19", treatment = c(0, 1), context = "CpG")

## The methylRawList for F2
generation_02 <- methRead(location = list("demo/F2_control_01.txt", "demo/F2_case_01.txt"), 
                    sample.id = list("F2_control_01", "F2_case_01"), 
                    assembly = "hg19", treatment = c(0, 1), context = "CpG")

## A list of methylRawList ready for methylInheritance
final_list <- list(generation_01, generation_02)
final_list
## [[1]]
## methylRawList object with 2 methylRaw objects
## 
## methylRaw object with 2 rows
## --------------
##   chr   start     end strand coverage numCs numTs
## 1  21 9764513 9764513      +      100    88    12
## 2  21 9764542 9764542      -       15     2    13
## --------------
## sample.id: F1_control_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## methylRaw object with 2 rows
## --------------
##   chr   start     end strand coverage numCs numTs
## 1  21 9764513 9764513      -       38    12    26
## 2  21 9764522 9764522      -       21     2    19
## --------------
## sample.id: F1_case_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## treatment: 0 1 
## 
## [[2]]
## methylRawList object with 2 methylRaw objects
## 
## methylRaw object with 2 rows
## --------------
##   chr   start     end strand coverage numCs numTs
## 1  21 9764514 9764514      +       40     0    40
## 2  21 9764522 9764522      +       30     2    28
## --------------
## sample.id: F2_control_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## methylRaw object with 2 rows
## --------------
##   chr   start     end strand coverage numCs numTs
## 1  21 9764513 9764513      +       33    12    21
## 2  21 9764533 9764533      -       23     1    22
## --------------
## sample.id: F2_case_01 
## assembly: hg19 
## context: CpG 
## resolution: base 
## 
## treatment: 0 1

More information about methRead function can be found in the documentation of the methylKit package.

Acknowledgment

We thank Marie-Pier Scott-Boyer for her advice on the vignette content.

Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.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.26.so;  LAPACK version 3.12.0
## 
## Random number generation:
##  RNG:     L'Ecuyer-CMRG 
##  Normal:  Inversion 
##  Sample:  Rejection 
##  
## 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] methylInheritance_1.31.0 methylKit_1.31.0         GenomicRanges_1.57.2    
## [4] GenomeInfoDb_1.41.2      IRanges_2.39.2           S4Vectors_0.43.2        
## [7] BiocGenerics_0.53.0      knitr_1.48               BiocStyle_2.35.0        
## 
## loaded via a namespace (and not attached):
##   [1] bitops_1.0-9                gridExtra_2.3              
##   [3] rlang_1.1.4                 magrittr_2.0.3             
##   [5] matrixStats_1.4.1           rebus.base_0.0-3           
##   [7] compiler_4.4.1              mgcv_1.9-1                 
##   [9] vctrs_0.6.5                 reshape2_1.4.4             
##  [11] stringr_1.5.1               pkgconfig_2.0.3            
##  [13] crayon_1.5.3                fastmap_1.2.0              
##  [15] XVector_0.45.0              labeling_0.4.3             
##  [17] utf8_1.2.4                  Rsamtools_2.21.2           
##  [19] rebus_0.1-3                 rmarkdown_2.28             
##  [21] UCSC.utils_1.1.0            xfun_0.48                  
##  [23] zlibbioc_1.51.2             cachem_1.1.0               
##  [25] jsonlite_1.8.9              highr_0.11                 
##  [27] DelayedArray_0.33.1         BiocParallel_1.41.0        
##  [29] parallel_4.4.1              R6_2.5.1                   
##  [31] bslib_0.8.0                 stringi_1.8.4              
##  [33] limma_3.61.12               rtracklayer_1.65.0         
##  [35] jquerylib_0.1.4             numDeriv_2016.8-1.1        
##  [37] Rcpp_1.0.13                 SummarizedExperiment_1.35.5
##  [39] R.utils_2.12.3              Matrix_1.7-1               
##  [41] splines_4.4.1               qvalue_2.37.0              
##  [43] abind_1.4-8                 yaml_2.3.10                
##  [45] codetools_0.2-20            curl_5.2.3                 
##  [47] lattice_0.22-6              tibble_3.2.1               
##  [49] rebus.datetimes_0.0-2       plyr_1.8.9                 
##  [51] Biobase_2.67.0              withr_3.0.2                
##  [53] coda_0.19-4.1               evaluate_1.0.1             
##  [55] rebus.numbers_0.0-1         mclust_6.1.1               
##  [57] Biostrings_2.75.0           pillar_1.9.0               
##  [59] BiocManager_1.30.25         MatrixGenerics_1.17.1      
##  [61] RCurl_1.98-1.16             emdbook_1.3.13             
##  [63] ggplot2_3.5.1               munsell_0.5.1              
##  [65] scales_1.3.0                gtools_3.9.5               
##  [67] glue_1.8.0                  maketools_1.3.1            
##  [69] tools_4.4.1                 BiocIO_1.17.0              
##  [71] sys_3.4.3                   data.table_1.16.2          
##  [73] GenomicAlignments_1.41.0    buildtools_1.0.0           
##  [75] mvtnorm_1.3-1               XML_3.99-0.17              
##  [77] grid_4.4.1                  bbmle_1.0.25.1             
##  [79] bdsmatrix_1.3-7             colorspace_2.1-1           
##  [81] nlme_3.1-166                GenomeInfoDbData_1.2.13    
##  [83] rebus.unicode_0.0-2         restfulr_0.0.15            
##  [85] cli_3.6.3                   fastseg_1.51.0             
##  [87] fansi_1.0.6                 S4Arrays_1.5.11            
##  [89] gtable_0.3.6                R.methodsS3_1.8.2          
##  [91] sass_0.4.9                  digest_0.6.37              
##  [93] SparseArray_1.5.45          farver_2.1.2               
##  [95] rjson_0.2.23                htmltools_0.5.8.1          
##  [97] R.oo_1.26.0                 lifecycle_1.0.4            
##  [99] httr_1.4.7                  statmod_1.5.0              
## [101] MASS_7.3-61