Modeling
24 October 2023
Differential Methylation
Here we demonstrate the analysis of DNA methylation dependent on one or more predictors. The predictors might be tissue, cell type, sex, age, tumor/normal, case/control or a combination of these factors. The DML (Differential Methylation Locus) function models β values (DNA methylation levels) using mixed linear models. This general supervised learning framework identifies CpG loci whose differential methylation is associated with known co-variates and can be used to perform epigenome-wide association studies (EWAS). Let’s first load the needed packages:
library(sesame)
library(SummarizedExperiment)
sesameDataCache() # required at new sesame installationIn the following, we will use an MM285 dataset of 10 mouse samples. This dataset contains mouse samples from different tissues and mice of different ages and sexes. The dataset is stored in a SummarizedExperiment object, which contains a data matrix combined with column-wise metadata accessible with colData:
se = sesameDataGet("MM285.10.SE.tissue")[1:1000,] # an arbitrary 1000 CpGs
cd = as.data.frame(colData(se)); rownames(cd) = NULL
cdCRITICAL: If your data contains NA, it is required that you exclude CpGs with missing levels. For example, one cannot assess sex-specific DNA methylation for a CpG that only has non-NA measurement on one sex. Exclusion of such CpGs for differential methylation modeling can be done using the checkLevels function. Here, we will check this for both sex and tissue:
se_ok = (checkLevels(assay(se), colData(se)$sex) &
checkLevels(assay(se), colData(se)$tissue))
sum(se_ok) # the number of CpGs that passes## [1] 995NOTE: If your model include discrete contrast variables like tissue and sex as in the current example, you should consider explicitly turning it into a factor variable with a reference level (we use treatment coding, see different coding systems).
For example, to use Colon as the reference tissue and Female as the reference sex, one can do the following
colData(se)$tissue <- relevel(factor(colData(se)$tissue), "Colon")
colData(se)$sex <- relevel(factor(colData(se)$sex), "Female")Then we will model DNA methylation variation treating tissue and sex as covariates. To do that we will call the DML function and specify the R formula ~tissue + sex. This function fits DNA methylation reading to a linear model and perform the corresponding slope test and goodness-of-fit test (F-test holding out each contrast variable). See also formula for how to specify lm/glm-like symbolic form for regression. All these results are returned in an object of class DMLSummary:
## DMLSummary Object with 995 Loci, 10 samples.## Contrasts: tissue, sexYou can use mc.cores argument to parallelize the computing.
Test Interpretation
The DMLSummary object is a list of slightly expanded summary.lm objects as is typically returned by summary(lm()). The summaryExtractTest function extracts some key test statistics from the DMLSummary object and store them in a data frame. Rows of the data frame correspond to CpGs/loci and columns contain the slopes and p-values of each variable.
test_result = summaryExtractTest(smry)
colnames(test_result) # the column names, show four groups of statistics## [1] "Probe_ID" "Est_X.Intercept." "Est_tissueCecum"
## [4] "Est_tissueEsophagus" "Est_tissueFat" "Est_tissueHeart"
## [7] "Est_sexMale" "Pval_X.Intercept." "Pval_tissueCecum"
## [10] "Pval_tissueEsophagus" "Pval_tissueFat" "Pval_tissueHeart"
## [13] "Pval_sexMale" "FPval_tissue" "FPval_sex"
## [16] "Eff_tissue" "Eff_sex"With the exception of the Intercept, there are four groups of columns, each starting with “Est_”, “Pval_”, “FPval_”, and “Eff_” as prefix. Here are what they represent:
Est_* : The slope estimate (aka the β coefficient, not to be confused with the DNA methylation β-value though) for continuous variable. DNA methylation difference of the current level with respect to the reference level for nominal contrast variables. Each suffix is concatenated from the contrast variable name (e.g., tissue, sex) and the level name if the contrast variable is discrete (e.g, Cecum, Esophagus, Fat). For example,
Est_tissueFatshould be interpreted as the estimated methylation level difference of Fat compared to the reference tissue (which isColon, as set above). If reference is not set, the first level in the alphabetic order is used as the reference level. There is a special column namedEst_`(Intercept)`. It corresponds to the base-level methylation of the reference (in this case a Female Colon sample).Pval_* : The unadjusted p-values of t-testing the slope. This represents the statistical significance of the methylation difference. For example,
Pval_tissueFattests whetherFatis significantly different fromColon(the reference level) in DNA methylation. ThePval_`(Intercept)`tests whether the reference level is significantly different from zero.FPval_* : The unadjusted p-value of the F-test contrasting the full model against a reduced model with the labeled contrast variable held out. Note that “Pval_” and “FPval_” are equivalent when the contrast variable is a 2-level factor, i.e., in the case of a pairwise comparison.
Eff_* : The effect size of each normial contrast variable. This is equivalent to the maximum slope subtracted by the minimum level including the reference level (0).
Multiple-testing adjustment can be done afterwards using R’s p.adjust function. It is integrated to the DMR function by default (see below).
Goodness of Fit
One may want to ask a question like
Is the CpG methylation tissue-specific?
rather than
Is the CpG more methylated in fat compared to liver?
The first question ask about the contrast variable as a whole while the second question concerns only a specific level in the contrast variable. To answer this question, we can use an F-test contasting the full model with a reduced model with the target contrast held out. By default, this statistics is already computed in the DML function. The test result is recorded in the FPval_ columns. For example, to get all CpGs that are methylated specific to sex,
library(dplyr)
library(tidyr)
test_result %>% dplyr::filter(FPval_sex < 0.05, Eff_sex > 0.1) %>%
select(FPval_sex, Eff_sex)Here we used 0.1 as the effect size threshold. This means DNA methylation difference under 0.1 (10%) is considered not biologically meaningful. This can be a valid assumption for homogenous cell composition as most cells would be either biallelically methylated, unmethylated or monoallelically methylated. But different threshold can be used in different analysis scenarios.
We can define CpG methylation as sex-specific, tissue-specific or both, by:
test_result %>%
mutate(sex_specific =
ifelse(FPval_sex < 0.05 & Eff_sex > 0.1, TRUE, FALSE)) %>%
mutate(tissue_specific =
ifelse(FPval_tissue < 0.05 & Eff_tissue > 0.1, TRUE, FALSE)) %>%
select(sex_specific, tissue_specific) %>% table## tissue_specific
## sex_specific FALSE TRUE
## FALSE 891 99
## TRUE 5 0As you can see from the result, some probes are sex-specific and others are tissue-specific. There is no overlap between probes whose methylation reading is differential along both contrasts.
Pairwise Comparison
From the test result, we can also ask whether the DNA methylation is different between two sexes or between two specific tissues. For example, Est_sexMale compares male from females. The following code creates a volcano plot.
Likewise, we can ask whether DNA methylation might be different between fat and colon. We can do
Continuous Predictors
The variable tested in the DML function can be continuous. Suppose we are interested in age besides sex. We will call the same function but with the following formula:
Let’s verify the CpGs positively associated with age.
df = data.frame(Age = colData(se)$age,
BetaValue = assay(se)[test_result2$Probe_ID[nrow(test_result2)],])
ggplot(df, aes(Age, BetaValue)) + geom_smooth(method="lm") + geom_point()## `geom_smooth()` using formula = 'y ~ x'
DMR
For a given contrast, one can merge neighboring CpGs that show consistent methylation variation into differentially methylated regions (DMRs). For example, we can merge sex-specific differential methylation identified above to chromosome X regions that show X-inactivation-related methylation difference. To do this, we need to pick a contrast:
## [1] "X.Intercept." "tissueCecum" "tissueEsophagus" "tissueFat"
## [5] "tissueHeart" "sexMale"## Merging correlated CpGs ... Done.
## Generated 509 segments.
## Combine p-values ...
## - 25 significant segments.
## - 1 significant segments (after BH).
## Done.See Supplemental Vignette for track-view visualization of the data.
Session Info
## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-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_GB 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: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_1.3.0 dplyr_1.1.3
## [3] ggplot2_3.4.4 tibble_3.2.1
## [5] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [7] MatrixGenerics_1.14.0 matrixStats_1.0.0
## [9] knitr_1.44 S4Vectors_0.40.0
## [11] IRanges_2.36.0 GenomicRanges_1.54.0
## [13] sesame_1.20.0 sesameData_1.19.0
## [15] ExperimentHub_2.10.0 AnnotationHub_3.10.0
## [17] BiocFileCache_2.10.0 dbplyr_2.3.4
## [19] BiocGenerics_0.48.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.1.3 bitops_1.0-7
## [3] rlang_1.1.1 magrittr_2.0.3
## [5] e1071_1.7-13 compiler_4.3.1
## [7] RSQLite_2.3.1 mgcv_1.9-0
## [9] png_0.1-8 vctrs_0.6.4
## [11] reshape2_1.4.4 stringr_1.5.0
## [13] pkgconfig_2.0.3 crayon_1.5.2
## [15] fastmap_1.1.1 XVector_0.42.0
## [17] ellipsis_0.3.2 fontawesome_0.5.2
## [19] labeling_0.4.3 utf8_1.2.4
## [21] promises_1.2.1 rmarkdown_2.25
## [23] tzdb_0.4.0 preprocessCore_1.64.0
## [25] purrr_1.0.2 bit_4.0.5
## [27] xfun_0.40 randomForest_4.7-1.1
## [29] zlibbioc_1.48.0 cachem_1.0.8
## [31] GenomeInfoDb_1.38.0 jsonlite_1.8.7
## [33] blob_1.2.4 later_1.3.1
## [35] DelayedArray_0.28.0 BiocParallel_1.36.0
## [37] interactiveDisplayBase_1.40.0 parallel_4.3.1
## [39] R6_2.5.1 bslib_0.5.1
## [41] stringi_1.7.12 RColorBrewer_1.1-3
## [43] jquerylib_0.1.4 Rcpp_1.0.11
## [45] wheatmap_0.2.0 readr_2.1.4
## [47] splines_4.3.1 httpuv_1.6.12
## [49] Matrix_1.6-1.1 tidyselect_1.2.0
## [51] abind_1.4-5 yaml_2.3.7
## [53] codetools_0.2-19 curl_5.1.0
## [55] lattice_0.22-5 plyr_1.8.9
## [57] shiny_1.7.5.1 withr_2.5.1
## [59] KEGGREST_1.42.0 evaluate_0.22
## [61] proxy_0.4-27 Biostrings_2.70.0
## [63] pillar_1.9.0 BiocManager_1.30.22
## [65] filelock_1.0.2 stats4_4.3.1
## [67] generics_0.1.3 RCurl_1.98-1.12
## [69] BiocVersion_3.18.0 hms_1.1.3
## [71] munsell_0.5.0 scales_1.2.1
## [73] BiocStyle_2.30.0 xtable_1.8-4
## [75] class_7.3-22 glue_1.6.2
## [77] tools_4.3.1 grid_4.3.1
## [79] AnnotationDbi_1.64.0 colorspace_2.1-0
## [81] nlme_3.1-163 GenomeInfoDbData_1.2.11
## [83] cli_3.6.1 rappdirs_0.3.3
## [85] fansi_1.0.5 S4Arrays_1.2.0
## [87] gtable_0.3.4 sass_0.4.7
## [89] digest_0.6.33 SparseArray_1.2.0
## [91] ggrepel_0.9.4 farver_2.1.1
## [93] memoise_2.0.1 htmltools_0.5.6.1
## [95] lifecycle_1.0.3 httr_1.4.7
## [97] mime_0.12 MASS_7.3-60
## [99] bit64_4.0.5