## ----setup, include = FALSE---------------------------------------------------
startTime <- Sys.time()
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    crop = NULL,
    message = FALSE,
    warning = FALSE
)

library(dnaEPICO)

can_run_preprocessing <- requireNamespace("minfiData", quietly = TRUE) &&
    requireNamespace("IlluminaHumanMethylation450kmanifest", quietly = TRUE) &&
    requireNamespace(
        "IlluminaHumanMethylation450kanno.ilmn12.hg19",
        quietly = TRUE
    )

can_run_sva <- can_run_preprocessing
can_run_glm <- requireNamespace(
    "IlluminaHumanMethylation450kanno.ilmn12.hg19",
    quietly = TRUE
)
can_run_lme <- can_run_glm && requireNamespace("lmerTest", quietly = TRUE)

## ----install-dnaepico, eval = FALSE-------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("paulYRP/dnaEPICO")
# 
# BiocManager::valid()

## ----load-dnaepico, eval = FALSE, message = FALSE-----------------------------
# library("dnaEPICO")

## ----makefile-example---------------------------------------------------------
makefilePath <- extractMake(destDir = tempdir(), overwrite = TRUE)
basename(makefilePath)

## ----preprocessing-pipeline-inputs, eval = can_run_preprocessing--------------
preprocessing_inputs <- dnaEPICO:::exampleMinfiIdatInputsDnaEpico(n = 6L)

names(preprocessing_inputs)
preprocessing_inputs$tempDir
basename(preprocessing_inputs$phenoFile)
head(
    basename(list.files(
        preprocessing_inputs$idatFolder,
        full.names = TRUE
    )),
    4
)
basename(preprocessing_inputs$probeExclusionPath)

## ----preprocessing-pipeline, eval = can_run_preprocessing---------------------
preprocessPipelineResult <- preprocessingMinfiEwasWater(
    phenoFile = preprocessing_inputs$phenoFile,
    idatFolder = preprocessing_inputs$idatFolder,
    outputLogs = file.path(preprocessing_inputs$tempDir, "logs"),
    nSamples = 6,
    SampleID = "Sample_Name",
    arrayType = preprocessing_inputs$arrayType,
    annotationVersion = preprocessing_inputs$annotationVersion,
    scriptLabel = "preprocessingMinfiEwasWater",
    baseDataFolder = file.path(preprocessing_inputs$tempDir, "rData"),
    figureBaseDir = file.path(preprocessing_inputs$tempDir, "figures"),
    detPThreshold = 1,
    normMethods = "quantile",
    sexColumn = "Sex",
    removeSexMismatch = FALSE,
    pvalThreshold = 1,
    chrToRemove = "",
    snpsToRemove = "SBE",
    mafThreshold = 1,
    probeExclusionPath = preprocessing_inputs$probeExclusionPath,
    plotGroupVar = "Sex",
    lcRef = "saliva",
    phenoOrder = "Sample_Name;Sex;Basename;Sentrix_ID;Sentrix_Position",
    lcPhenoDir = file.path(
        preprocessing_inputs$tempDir,
        "data",
        "preprocessingMinfiEwasWater"
    ),
    display = FALSE,
    verbose = FALSE,
    logs = TRUE,
    saveOutputs = TRUE
)

preprocess_paths <- c(
    phenoLC = file.path(
        preprocessing_inputs$tempDir,
        "data",
        "preprocessingMinfiEwasWater",
        "phenoLC.csv"
    ),
    rgset = file.path(
        preprocessing_inputs$tempDir,
        "rData",
        "preprocessingMinfiEwasWater",
        "objects",
        "RGSet.RData"
    ),
    beta = file.path(
        preprocessing_inputs$tempDir,
        "rData",
        "preprocessingMinfiEwasWater",
        "metrics",
        "beta_NomFilt_MSetF_Flt_Rxy_Ds_Rc.RData"
    ),
    qc = file.path(
        preprocessing_inputs$tempDir,
        "figures",
        "preprocessingMinfiEwasWater",
        "qc",
        "quality_control(MSet).tiff"
    )
)

class(preprocessPipelineResult)
basename(preprocess_paths)
file.exists(preprocess_paths)

## ----sva-pipeline-inputs, eval = can_run_sva, message = FALSE, warning = FALSE----
sva_targets_file <- file.path(
    preprocessing_inputs$tempDir,
    "data",
    "preprocessingMinfiEwasWater",
    "phenoLC.csv"
)
sva_rgset_file <- file.path(
    preprocessing_inputs$tempDir,
    "rData",
    "preprocessingMinfiEwasWater",
    "objects",
    "RGSet.RData"
)

basename(c(sva_targets_file, sva_rgset_file))
file.exists(c(sva_targets_file, sva_rgset_file))

## ----sva-pipeline, eval = can_run_sva, message = FALSE, warning = FALSE-------
svaPipelineResult <- svaEnmix(
    phenoFile = sva_targets_file,
    rgsetData = sva_rgset_file,
    outputLogs = file.path(preprocessing_inputs$tempDir, "logs"),
    SampleID = "Sample_Name",
    arrayType = "IlluminaHumanMethylation450k",
    annotationVersion = "ilmn12.hg19",
    SentrixIDColumn = "Sentrix_ID",
    SentrixPositionColumn = "Sentrix_Position",
    ctrlSvaPercVar = 0.90,
    ctrlSvaFlag = 1,
    scriptLabel = "svaEnmix",
    dataBaseDir = file.path(preprocessing_inputs$tempDir, "data"),
    rBaseDir = file.path(preprocessing_inputs$tempDir, "rData"),
    figureBaseDir = file.path(preprocessing_inputs$tempDir, "figures"),
    display = FALSE,
    verbose = FALSE,
    logs = TRUE,
    saveOutputs = TRUE
)

class(svaPipelineResult$savedFiles)
names(svaPipelineResult$savedFiles)
basename(unlist(svaPipelineResult$savedFiles, use.names = FALSE))

## ----preprocessingPheno-pipeline-inputs, message = FALSE, warning = FALSE-----
pheno_inputs <- dnaEPICO:::examplePreprocessingPhenoStateDnaEpico()

names(pheno_inputs)
pheno_inputs$tempDir
basename(c(
    pheno_inputs$phenoPath,
    pheno_inputs$betaPath,
    pheno_inputs$mPath,
    pheno_inputs$cnPath
))

## ----preprocessingPheno-pipeline, message = FALSE, warning = FALSE------------
phenoPipelineResult <- preprocessingPheno(
    phenoFile = pheno_inputs$phenoPath,
    betaPath = pheno_inputs$betaPath,
    mPath = pheno_inputs$mPath,
    cnPath = pheno_inputs$cnPath,
    SampleID = "Sample_Name",
    timeVar = "Timepoint",
    timepoints = "1,2",
    combineTimepoints = "1,2",
    outputPheno = file.path(pheno_inputs$tempDir, "data", "preprocessingPheno"),
    outputRData = file.path(
        pheno_inputs$tempDir,
        "rData",
        "preprocessingPheno",
        "metrics"
    ),
    outputRDataMerge = file.path(
        pheno_inputs$tempDir,
        "rData",
        "preprocessingPheno",
        "mergeData"
    ),
    sexColumn = "Sex",
    outputLogs = file.path(pheno_inputs$tempDir, "logs"),
    outputDir = file.path(pheno_inputs$tempDir, "clockFoundation"),
    verbose = FALSE,
    logs = TRUE,
    saveOutputs = TRUE
)

names(phenoPipelineResult$savedFiles)
names(phenoPipelineResult$savedFiles$timepoints)
basename(unlist(phenoPipelineResult$savedFiles$timepoints[["1"]]))
basename(unlist(phenoPipelineResult$savedFiles[c(
    "combinedPheno",
    "combinedPhenoMethylation",
    "betaCSV",
    "phenoCF"
)]))

## ----glm-pipeline-inputs, eval = can_run_glm, message = FALSE, warning = FALSE----
glm_inputs <- dnaEPICO:::exampleMethylationGLMStateDnaEpico()

names(glm_inputs)
glm_inputs$tempDir
basename(glm_inputs$inputPath)
file.exists(glm_inputs$inputPath)

## ----glm-pipeline, eval = can_run_glm, message = FALSE, warning = FALSE-------
glmPipelineResult <- methylationGLM(
    inputPheno = glm_inputs$inputPath,
    phenotypes = "status",
    covariates = "sex,age",
    factorVars = "status,sex",
    scaleVars = "age",
    cpgLimit = 2,
    nCores = 1,
    outputLogs = file.path(glm_inputs$tempDir, "logs"),
    outputRData = file.path(glm_inputs$tempDir, "rData", "models"),
    outputPlots = file.path(glm_inputs$tempDir, "figures"),
    significantCpGDir = file.path(glm_inputs$tempDir, "significant"),
    summaryTxtDir = file.path(glm_inputs$tempDir, "summaries"),
    summaryPval = 1,
    significantCpGPval = 1,
    annotationPackage = "IlluminaHumanMethylation450kanno.ilmn12.hg19",
    annotationCols = "Name,chr,pos",
    annotatedGLMOut = file.path(glm_inputs$tempDir, "annotated"),
    reportAssetsDir = file.path(
        preprocessing_inputs$tempDir,
        "reports", "example", "assets", "results", "glm_results"
    ),
    display = FALSE,
    verbose = FALSE,
    logs = TRUE,
    saveOutputs = TRUE
)

class(glmPipelineResult$savedFiles)
names(glmPipelineResult$savedFiles)

## ----lme-pipeline-inputs, eval = can_run_lme, message = FALSE, warning = FALSE----
lme_inputs <- dnaEPICO:::exampleMethylationLMEStateDnaEpico()

names(lme_inputs)
lme_inputs$tempDir
basename(lme_inputs$inputPath)
file.exists(lme_inputs$inputPath)

## ----lme-pipeline, eval = can_run_lme, message = FALSE, warning = FALSE-------
lmePipelineResult <- methylationLME(
    inputPheno = lme_inputs$inputPath,
    outputLogs = file.path(lme_inputs$tempDir, "logs"),
    outputRData = file.path(lme_inputs$tempDir, "rData", "models"),
    outputPlots = file.path(lme_inputs$tempDir, "figures"),
    personVar = "person",
    SampleID = "Sample_Name",
    timeVar = "Timepoint",
    phenotypes = "score",
    covariates = "sex",
    factorVars = "sex",
    scaleVars = "score",
    cpgLimit = 2,
    nCores = 1,
    summaryPval = 1,
    saveSignificantInteractions = TRUE,
    significantInteractionDir = file.path(
        lme_inputs$tempDir,
        "results",
        "cpgs",
        "methylationLME"
    ),
    significantInteractionPval = 1,
    saveTxtSummaries = TRUE,
    summaryTxtDir = file.path(
        lme_inputs$tempDir,
        "results",
        "summary",
        "methylationLME"
    ),
    annotationPackage = "IlluminaHumanMethylation450kanno.ilmn12.hg19",
    annotationCols = "Name,chr,pos",
    annotatedLMEOut = file.path(lme_inputs$tempDir, "annotated"),
    reportAssetsDir = file.path(
        preprocessing_inputs$tempDir,
        "reports", "example", "assets", "results", "lme_results"
    ),
    display = FALSE,
    verbose = FALSE,
    logs = TRUE,
    saveOutputs = TRUE
)

class(lmePipelineResult$savedFiles)
names(lmePipelineResult$savedFiles)

## ----lme-omnibus-example, eval = FALSE----------------------------------------
# omnibusResult <- methylationLME(
#     inputPheno = "rData/model/preprocessingPheno/mergeData/phenoBT1T2T3.RData",
#     personVar = "Participant",
#     timeVar = "Timepoint",
#     phenotypes = "Timepoint",
#     covariates = "Age,Sex,PC1,PC2,PC3",
#     factorVars = "Timepoint,Profession,Sex",
#     interactionTerm = "Profession",
#     omnibusTest = TRUE,
#     omnibusDdf = "Satterthwaite"
# )

## ----quarto-check-r, eval = FALSE---------------------------------------------
# Sys.which("quarto")
# system2("quarto", "--version")

## ----quarto-path-r, eval = FALSE----------------------------------------------
# Sys.setenv(QUARTO_BIN = "/full/path/to/quarto")

## ----dnamReport-pipeline-example, eval = FALSE--------------------------------
# report_log_dir <- file.path(preprocessing_inputs$tempDir, "logs", "report")
# dir.create(report_log_dir, recursive = TRUE, showWarnings = FALSE)
# 
# report_logs <- c(
#     file.path(
#         preprocessing_inputs$tempDir,
#         "logs",
#         "log_preprocessingMinfiEwasWater.txt"
#     ),
#     file.path(pheno_inputs$tempDir, "logs", "log_preprocessingPheno.txt"),
#     file.path(preprocessing_inputs$tempDir, "logs", "log_svaEnmix.txt"),
#     file.path(glm_inputs$tempDir, "logs", "log_methylationGLM.txt"),
#     file.path(lme_inputs$tempDir, "logs", "log_methylationLME.txt")
# )
# 
# existing_logs <- report_logs[file.exists(report_logs)]
# if (length(existing_logs)) {
#     file.copy(existing_logs, report_log_dir, overwrite = TRUE)
# }
# 
# reportResult <- dnamReport(
#     outputDir = file.path(preprocessing_inputs$tempDir, "reports", "example"),
#     phenoTab = file.path(
#         preprocessing_inputs$tempDir,
#         "data",
#         "preprocessingMinfiEwasWater",
#         "phenoLC.csv"
#     ),
#     enmixTab = file.path(
#         preprocessing_inputs$tempDir,
#         "figures",
#         "preprocessingMinfiEwasWater",
#         "enmix"
#     ),
#     qcTab = file.path(
#         preprocessing_inputs$tempDir,
#         "figures",
#         "preprocessingMinfiEwasWater",
#         "qc"
#     ),
#     svaTab = file.path(preprocessing_inputs$tempDir, "figures", "svaEnmix"),
#     metricTab = file.path(
#         preprocessing_inputs$tempDir,
#         "figures",
#         "preprocessingMinfiEwasWater",
#         "metrics"
#     ),
#     glmTab = glmPipelineResult$savedFiles$annotatedGLM,
#     lmeTab = lmePipelineResult$savedFiles$annotatedLME,
#     logTab = report_log_dir,
#     detPPath = file.path(
#         preprocessing_inputs$tempDir,
#         "rData",
#         "preprocessingMinfiEwasWater",
#         "qc",
#         "detP_RGSet.RData"
#     ),
#     detPThreshold = 0.01,
#     verbose = FALSE,
#     logs = TRUE
# )
# 
# reportResult$outputFile

## ----makefile-template-full, echo = FALSE, results = "asis"-------------------
makefile_lines <- readLines(makefilePath, warn = FALSE)
cat("```make\n", paste(makefile_lines, collapse = "\n"), "\n```\n", sep = "")

## ----vignette-date, echo = FALSE----------------------------------------------
Sys.time()

## ----vignette-wallclock, echo = FALSE-----------------------------------------
totalTime <- diff(c(startTime, Sys.time()))
round(totalTime, digits = 3)

## ----vignette-session-info, echo = FALSE--------------------------------------
sessionInfo()