--- title: "Modeling" date: "`r BiocStyle::doc_date()`" package: sesame output: BiocStyle::html_document fig_width: 8 fig_height: 6 vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{3. Modeling} %\VignetteEncoding{UTF-8} --- # Linear Modeling Here I will try to address the analysis workflow of studying DNA methylation with one or more independent variable, be it tissue, cell type, sex, age, tumor/normal, treatment/control or a combination of these factors. The `DML` (short for Differential Methylation Locus) function models beta values (DNA methylation levels) on known co-variates (e.g., case vs control) using linear modeling. This general supervised learning framework identifies CpG loci whose differential methylation is associated with known co-variates. Let's first load needed packages ```{r message=FALSE} library(sesame) library(dplyr) library(tidyr) library(SummarizedExperiment) library(ggplot2) ``` In the following, we will use 10 mouse array samples as our example dataset. This dataset contains mouse samples from different tissues and individuals 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`: ```{r message=FALSE} se = sesameDataGet("MM285.10.tissues")[1:1000,] # use a random 1000 probes colData(se) ``` > **CRITICAL:** If your data is NA-masked (see the Basic Usage tutorial), it is critical to preclude CpGs with missing levels. For example, if a CpG does not contain any non-NA values for all male mice, one cannot assess sex-specific DNA methylation for this CpG. One needs to preclude such CpGs from the modeling of sex-specific differential methylation. This can be done using the `checkLevels` function. In the following, we will test for both sex and tissue as co-variates using the `checkLevels` function: ```{r} se_ok = (checkLevels(assay(se), colData(se)$sex) & checkLevels(assay(se), colData(se)$tissue)) sum(se_ok) se = 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]( https://stats.idre.ucla.edu/r/library/r-library-contrast-coding-systems- for-categorical-variables/)). For example, to use `Colon` as the reference tissue and `Female` as the reference sex, one can do the following ```{r} 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). All these results are returned in an object of class `DMLSummary`; ```{r} smry = DML(se, ~tissue + sex) smry ``` You can use `mc.cores` argument to parallelize the computing. ## Test Interpretation The `DMLSummary` object is a list of expanded `summary.lm` objects from the linear model (as would be returned by `summary(lm())`. The `summaryExtractTest` function abstracts a `DMLSummary` object into a data frame with rows corresponding to CpGs/loci and columns containing the slopes and p-values of each variable. ```{r} test_result = summaryExtractTest(smry) colnames(test_result) # the column names, show four groups of statistics 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 beta coefficient) 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_tissueFat` should be interpreted as the estimated beta value difference of Fat compared to the reference tissue (which is `Colon`, as set above. If reference is not set, the first level in alphabetic order is used as the reference level). There are is a special column named ``Est_`(Intercept)` ``. They correspond to the beta value of the reference level (in this case a Female Colon sample). * **Pval_** : The p-value of t-testing the slope (statistical significance of the difference). For example, `Pval_tissueFat` tests whether `Fat` is significantly different from `Colon` (the reference level) in DNA methylation. The ``Pval_`(Intercept)` `` tests whether the reference level is significantly different from zero. * **FPval_** : The p-value of the F-test contrasting the full model against a reduced model with the labeled contrasting held out. Note that "Pval_" and "FPval_" are equivalent when the contrast variable is a 2-level factor, i.e., 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). ## Goodness of fit One may want to ask a question like whether the methylation of the probe is tissue-specific rather than the methylation of the probe is specifically hypermethylated in fat. The former question ask about the contrast variable as a whole instead of a specific level of the contrast variable. To answer this question, we can use a F-test contasting the full model with a reduced model with the target contrast variable held out. This information is captured in the **FPval_** columns. For example, to get all CpGs that are methylated specific to sex, ```{r} test_result %>% dplyr::filter(FPval_sex < 0.05, Eff_sex > 0.1) %>% select(Pval_sexMale, Eff_sex) ``` Here we used 0.1 as the effect size threshold. This means DNA methylation difference under 0.1 (10%) is not considered a biologically meaningful difference. This is a valid assumption for homogenous cell composition as most cells would be either biallelically methylated, unmethylated or monoallelically methylated. We can define CpG methylation as sex-specific, tissue-specific or both ```{r} 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 ``` As you can see from the result, two probes are sex-specific and one probe 13 probes are tissue-specific. There is no overlap between probes whose methylation reading is differential across both contrasts. ## Pairwise Comparison From the test result, we can also ask whether the DNA methylation might be different between two sexes or two tissues. For example, `Est_sexMale` compares male from females. The following code creates a volcano plot. ```{r} 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 ```{r} ggplot(test_result) + geom_point(aes(Est_tissueFat, -log10(Pval_tissueFat))) ``` ## Continuous Variable # DMR Testing sex-specific differential methylation yields chrX-linked probes. ```{r} cf_list = summaryExtractCfList(smry) cf_list = DMR(se, cf_list$sexMale) topSegments(cf_list) %>% dplyr::filter(Seg.Pval.adj < 0.05) ``` ```{r message=FALSE, warning=FALSE, include=FALSE} library(sesame) library(dplyr) options(rmarkdown.html_vignette.check_title = FALSE) ``` # Track View To visualize all probes from a gene ```{r, message=FALSE, fig.width=6, fig.height=5} betas <- sesameDataGet('HM450.10.TCGA.PAAD.normal') visualizeGene('DNMT1', betas, platform='HM450') ``` To visualize probes from arbitrary region ```{r, message=FALSE, fig.width=6, fig.height=5} visualizeRegion( 'chr19',10260000,10380000, betas, platform='HM450', show.probeNames = FALSE) ``` To visualize by probe names ```{r, message=FALSE, fig.width=6} visualizeProbes(c("cg02382400", "cg03738669"), betas, platform='HM450') ```