Modeling
4 May 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] 995se = se[se_ok,]NOTE: 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:
smry = DML(se, ~tissue + sex)
smry## 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"head(test_result)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.
library(ggplot2)
ggplot(test_result) + geom_point(aes(Est_sexMale, -log10(Pval_sexMale)))
Likewise, we can ask whether DNA methylation might be different between fat and colon. We can do
ggplot(test_result) + geom_point(aes(Est_tissueFat, -log10(Pval_tissueFat)))
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:
smry2 = DML(se, ~ age + sex)
test_result2 = summaryExtractTest(smry2) %>% arrange(Est_age)Let’s verify the CpGs positively associated with age.
test_result2 %>% dplyr::select(Probe_ID, Est_age, Pval_age) %>% taildf = 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:
dmContrasts(smry) # pick a contrast from below## [1] "X.Intercept." "tissueCecum" "tissueEsophagus" "tissueFat"
## [5] "tissueHeart" "sexMale"merged = DMR(se, smry, "sexMale", platform="MM285") # merge CpGs to regions## Merging correlated CpGs ... Done.
## Generated 509 segments.
## Combine p-values ...
## - 25 significant segments.
## - 1 significant segments (after BH).
## Done.merged %>% dplyr::filter(Seg_Pval_adj < 0.01)See Supplemental Vignette for track-view visualization of the data.
Session Info
sessionInfo()## R version 4.3.0 (2023-04-21)
## Platform: aarch64-apple-darwin22.3.0 (64-bit)
## Running under: macOS Ventura 13.2.1
##
## Matrix products: default
## BLAS: /Users/zhouw3/.Renv/versions/4.3.0/lib/R/lib/libRblas.dylib
## LAPACK: /Users/zhouw3/.Renv/versions/4.3.0/lib/R/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggplot2_3.4.2 tidyr_1.3.0
## [3] dplyr_1.1.2 SummarizedExperiment_1.30.1
## [5] Biobase_2.60.0 GenomicRanges_1.52.0
## [7] GenomeInfoDb_1.36.0 IRanges_2.34.0
## [9] S4Vectors_0.38.1 MatrixGenerics_1.12.0
## [11] matrixStats_0.63.0 sesame_1.18.4
## [13] sesameData_1.18.0 ExperimentHub_2.8.0
## [15] AnnotationHub_3.8.0 BiocFileCache_2.8.0
## [17] dbplyr_2.3.2 BiocGenerics_0.46.0
## [19] rmarkdown_2.21 R6_2.5.1
##
## 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] compiler_4.3.0 RSQLite_2.3.1
## [7] mgcv_1.8-42 png_0.1-8
## [9] vctrs_0.6.2 reshape2_1.4.4
## [11] stringr_1.5.0 pkgconfig_2.0.3
## [13] crayon_1.5.2 fastmap_1.1.1
## [15] XVector_0.40.0 ellipsis_0.3.2
## [17] labeling_0.4.2 utf8_1.2.3
## [19] promises_1.2.0.1 tzdb_0.3.0
## [21] preprocessCore_1.62.1 purrr_1.0.1
## [23] bit_4.0.5 xfun_0.39
## [25] zlibbioc_1.46.0 cachem_1.0.8
## [27] jsonlite_1.8.4 blob_1.2.4
## [29] highr_0.10 later_1.3.1
## [31] DelayedArray_0.26.1 BiocParallel_1.34.0
## [33] interactiveDisplayBase_1.38.0 parallel_4.3.0
## [35] bslib_0.4.2 stringi_1.7.12
## [37] RColorBrewer_1.1-3 jquerylib_0.1.4
## [39] Rcpp_1.0.10 knitr_1.42
## [41] wheatmap_0.2.0 readr_2.1.4
## [43] splines_4.3.0 httpuv_1.6.9
## [45] Matrix_1.5-4 tidyselect_1.2.0
## [47] yaml_2.3.7 codetools_0.2-19
## [49] curl_5.0.0 lattice_0.21-8
## [51] tibble_3.2.1 plyr_1.8.8
## [53] shiny_1.7.4 withr_2.5.0
## [55] KEGGREST_1.40.0 evaluate_0.20
## [57] Biostrings_2.68.0 pillar_1.9.0
## [59] BiocManager_1.30.20 filelock_1.0.2
## [61] generics_0.1.3 RCurl_1.98-1.12
## [63] BiocVersion_3.17.1 hms_1.1.3
## [65] munsell_0.5.0 scales_1.2.1
## [67] BiocStyle_2.27.2 xtable_1.8-4
## [69] glue_1.6.2 tools_4.3.0
## [71] grid_4.3.0 AnnotationDbi_1.62.1
## [73] colorspace_2.1-0 nlme_3.1-162
## [75] GenomeInfoDbData_1.2.10 cli_3.6.1
## [77] rappdirs_0.3.3 fansi_1.0.4
## [79] S4Arrays_1.0.1 gtable_0.3.3
## [81] sass_0.4.6 digest_0.6.31
## [83] farver_2.1.1 memoise_2.0.1
## [85] htmltools_0.5.5 lifecycle_1.0.3
## [87] httr_1.4.5 mime_0.12
## [89] bit64_4.0.5