## ----model0, include=FALSE----------------------------------------------------
options(rmarkdown.html_vignette.check_title = FALSE)

## ----model1, message=FALSE, warning=FALSE-------------------------------------
library(sesame)
library(SummarizedExperiment)
sesameDataCache() # required at new sesame installation

## ----model2, message=FALSE----------------------------------------------------
se = sesameDataGet("MM285.10.SE.tissue")[1:1000,] # an arbitrary 1000 CpGs
cd = as.data.frame(colData(se)); rownames(cd) = NULL
cd

## ----model3-------------------------------------------------------------------
se_ok = (checkLevels(assay(se), colData(se)$sex) &
    checkLevels(assay(se), colData(se)$tissue))
sum(se_ok)                      # the number of CpGs that passes
se = se[se_ok,]

## ----model4-------------------------------------------------------------------
colData(se)$tissue <- relevel(factor(colData(se)$tissue), "Colon")
colData(se)$sex <- relevel(factor(colData(se)$sex), "Female")

## ----model5-------------------------------------------------------------------
smry = DML(se, ~tissue + sex)
smry

## ----model6-------------------------------------------------------------------
test_result = summaryExtractTest(smry)
colnames(test_result) # the column names, show four groups of statistics
head(test_result)

## ----model7, message = FALSE--------------------------------------------------
library(dplyr)
library(tidyr)
test_result %>% dplyr::filter(FPval_sex < 0.05, Eff_sex > 0.1) %>%
    select(FPval_sex, Eff_sex)

## ----model8-------------------------------------------------------------------
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

## ----model9-------------------------------------------------------------------
library(ggplot2)
ggplot(test_result) + geom_point(aes(Est_sexMale, -log10(Pval_sexMale)))

## ----model10------------------------------------------------------------------
ggplot(test_result) + geom_point(aes(Est_tissueFat, -log10(Pval_tissueFat)))

## ----model11------------------------------------------------------------------
smry2 = DML(se, ~ age + sex)
test_result2 = summaryExtractTest(smry2) %>% arrange(Est_age)

## ----model12------------------------------------------------------------------
test_result2 %>% dplyr::select(Probe_ID, Est_age, Pval_age) %>% tail
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()

## ----model13, eval=TRUE-------------------------------------------------------
dmContrasts(smry)                       # pick a contrast from below
merged = DMR(se, smry, "sexMale", platform="MM285") # merge CpGs to regions
merged %>% dplyr::filter(Seg_Pval_adj < 0.01)

## -----------------------------------------------------------------------------
sessionInfo()