dnaEPICO 0.99.38
This vignette describes the file-based dnaEPICO workflow. Functions return
structured objects and write files when saveOutputs = TRUE. Use this workflow
for reproducible directories, Make execution, and resumable analyses on local
or HPC systems.
The dnaEPICO package uses methods from several external
packages. Because no single manuscript describes all components, the guidance
below explains how to cite dnaEPICO depending on the functions you use.
This citation guidance is adapted from the vignettes and user guides of minfi, ENmix, and wateRmelon.
make f3 MODEL=model1, please cite Aryee et al. (2014); Fortin, Triche Jr., and Hansen (2017); Xu, Niu, and Taylor (2021);
Pidsley et al. (2013); Maksimovic, Gordon, and Oshlack (2012); Fortin et al. (2014); Triche et al. (2013); Touleimat and Tost (2012); Murat et al. (2020).make f4 MODEL=model1, please cite Aryee et al. (2014); Fortin, Triche Jr., and Hansen (2017); Xu, Niu, and Taylor (2021);
Pidsley et al. (2013); Maksimovic, Gordon, and Oshlack (2012); Fortin et al. (2014); Triche et al. (2013); Touleimat and Tost (2012); Murat et al. (2020); Marschner (2011).make f3lme MODEL=model1, please cite Aryee et al. (2014); Fortin, Triche Jr., and Hansen (2017);
Xu, Niu, and Taylor (2021); Pidsley et al. (2013); Maksimovic, Gordon, and Oshlack (2012); Fortin et al. (2014); Triche et al. (2013); Touleimat and Tost (2012); Murat et al. (2020); Bates et al. (2015).make all MODEL=model1, please cite Aryee et al. (2014); Fortin, Triche Jr., and Hansen (2017); Xu, Niu, and Taylor (2021);
Pidsley et al. (2013); Maksimovic, Gordon, and Oshlack (2012); Fortin et al. (2014); Triche et al. (2013); Touleimat and Tost (2012); Murat et al. (2020); Marschner (2011); Bates et al. (2015).dnaEPICO is built on core Bioconductor infrastructure for high-dimensional genomic data, with a focus on Illumina DNA methylation arrays. This vignette assumes familiarity with a general DNA methylation workflow. For an introduction, see this tutorial: https://paulyrp.github.io/2025-cpgpneurogenomics-workshop/tutorial.html. It covers the main concepts and analysis steps.
Preprocessing and quality control are performed using established Bioconductor
tools, including minfi, ENmix, and
wateRmelon. Downstream statistical modelling uses CpG-wise
generalised linear models and longitudinal mixed-effects models fitted with
lmerTest/lme4 or nlme. The nlme engine can include AR1 or CAR1
within-participant residual correlation. Users are expected to have basic
familiarity with R, command-line execution, and Illumina IDAT file structures.
For an introduction to Bioconductor, see the installation guide.
The lists of excluded probes depend on the Illumina methylation-array platform.
For Illumina HumanMethylationEPIC v2.0, use the cross-reactive probe-exclusion file from Peters et al. (2024).
For Illumina MethylationEPIC, also known as the 850k array, use
the probe-exclusion resources from
Pidsley et al. (2016).
The supporting files commonly used together are
13059_2016_1066_MOESM1_ESM.csv,
13059_2016_1066_MOESM4_ESM.csv,
13059_2016_1066_MOESM5_ESM.csv, and
13059_2016_1066_MOESM6_ESM.csv.
For Illumina HumanMethylation450k, use the cross-reactive and polymorphic probe resources from Chen et al. (2013).
Multiple probe-exclusion files can be supplied as a semicolon-separated
value in probeExclusionPath. dnaEPICO reads probe IDs from each file,
uses probeExclusionIdColumn when supplied, or otherwise auto-detects
common probe-ID columns such as ProbeID, TargetID, IlmnID, and
Name. The unique union of all probe IDs is then used to filter the
normalised object.
For EPICv2, setting useEpicV2Manifest = TRUE also retrieves the
expanded Peters et al. manifest from AnnotationHub resource
AH116484. Probes
flagged in selected manifest columns are added to the same exclusion set.
By default, probes flagged by CH_WGBS_evidence, CH_BLAT, or
MissingPos are removed, while MismatchPos is retained unless explicitly
enabled.
Install the development version of dnaEPICO and its dependencies with BiocManager:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("paulYRP/dnaEPICO")
BiocManager::valid()library("dnaEPICO")The package includes a template Makefile that can be copied into a project directory and adapted to local paths and model definitions.
makefilePath <- extractMake(destDir = tempdir(), overwrite = TRUE)
basename(makefilePath)
#> [1] "Makefile"After extraction, the file should be edited so that the project-specific paths, model labels, and cluster settings match the target analysis environment.
Arguments:
destDir = tempdir() writes the template to a temporary directory so the
example does not modify the working tree.overwrite = TRUE allows the template to be replaced if it already exists in
that temporary location.The file-based route usually starts with preprocessingMinfiEwasWater(). This
step writes the filtered RGSet, normalized metric matrices, quality-control
figures, and the phenoLC.csv file that will be consumed by svaEnmix() and
preprocessingPheno().
The next chunk shows how to recreate the temporary files used in this example. It writes a small phenotype table to a temporary directory, copies a small set of IDAT files, and records the probe-exclusion reference used during probe filtering.
preprocessing_inputs <- dnaEPICO:::exampleMinfiIdatInputsDnaEpico(n = 6L)
names(preprocessing_inputs)
#> [1] "tempDir" "idatFolder" "phenoFile"
#> [4] "targets" "arrayType" "annotationVersion"
#> [7] "probeExclusionPath"
preprocessing_inputs$tempDir
#> [1] "/tmp/RtmpBnguYi/dnaEPICO-idat-example-1259cbcf22fbf"
basename(preprocessing_inputs$phenoFile)
#> [1] "pheno.csv"
head(
basename(list.files(
preprocessing_inputs$idatFolder,
full.names = TRUE
)),
4
)
#> [1] "5723646052_R02C02_Grn.idat" "5723646052_R02C02_Red.idat"
#> [3] "5723646052_R04C01_Grn.idat" "5723646052_R04C01_Red.idat"
basename(preprocessing_inputs$probeExclusionPath)
#> [1] "12864_2024_10027_MOESM8_ESM.csv"preprocessing_inputs is a small list returned by the internal example helper.
names(preprocessing_inputs) shows its main elements: tempDir is the
temporary working directory, idatFolder contains the copied example IDAT
files, phenoFile is the phenotype table used by the function, targets is the
same phenotype information already loaded into memory, arrayType and
annotationVersion describe the array platform, and probeExclusionPath
points to the probe-exclusion reference used during probe filtering.
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
)
#> Plotting STAINING .jpg
#> Plotting EXTENSION .jpg
#> Plotting HYBRIDIZATION .jpg
#> Plotting TARGET_REMOVAL .jpg
#> Plotting BISULFITE_CONVERSION_I .jpg
#> Plotting BISULFITE_CONVERSION_II .jpg
#> Plotting SPECIFICITY_I .jpg
#> Plotting SPECIFICITY_II .jpg
#> Plotting NON-POLYMORPHIC .jpg
#> Plotting NEGATIVE .jpg
#> Plotting NORM_A .jpg
#> Plotting NORM_C .jpg
#> Plotting NORM_G .jpg
#> Plotting NORM_T .jpg
#> Plotting NORM_ACGT .jpg
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)
#> [1] "dnaEPICO_preprocessingMinfiEwasWater"
basename(preprocess_paths)
#> [1] "phenoLC.csv"
#> [2] "RGSet.RData"
#> [3] "beta_NomFilt_MSetF_Flt_Rxy_Ds_Rc.RData"
#> [4] "quality_control(MSet).tiff"
file.exists(preprocess_paths)
#> [1] TRUE TRUE TRUE TRUEThe returned object still has class dnaEPICO_preprocessingMinfiEwasWater, but
the printed vectors in this chunk focus on the written outputs. In
preprocess_paths, phenoLC is the phenotype table with estimated cell
proportions, rgset is the filtered RGSet, beta is the saved beta matrix
used later by phenotype preparation, and qc is an example quality-control
figure. basename(preprocess_paths) shows the file names, while
file.exists(preprocess_paths) confirms that each one was written successfully.
Arguments:
phenoFile and idatFolder define the phenotype table and the IDAT folder
used to start the pipeline.SampleID = "Sample_Name" is the key used to align phenotype rows with the
methylation objects.arrayType and annotationVersion define the array manifest and annotation
used by minfi.baseDataFolder, figureBaseDir, and lcPhenoDir define where serialized
objects, figures, and phenoLC.csv are written.detPThreshold, normMethods, pvalThreshold, chrToRemove,
snpsToRemove, mafThreshold, and probeExclusionPath define the main
preprocessing and filtering choices.sexColumn identifies reported sex for comparison with methylation-predicted
sex. removeSexMismatch = FALSE retains and reports confirmed mismatches;
TRUE removes samples only when both values are known and disagree. Missing
or unknown sex values are retained.plotGroupVar, lcRef, and phenoOrder control QC grouping, cell-type
estimation, and the structure of the exported phenotype table.saveOutputs = TRUE is what turns this function into the first file-producing
step of the pipeline.The next example shows svaEnmix() in file-writing mode. The function still
returns a structured result object, but the savedFiles element now records the
main output files written to disk.
In the file-based route, svaEnmix() consumes the phenoLC.csv and
RGSet.RData files written by preprocessingMinfiEwasWater(). The next chunk
shows how those temporary file paths are reconstructed locally.
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))
#> [1] "phenoLC.csv" "RGSet.RData"
file.exists(c(sva_targets_file, sva_rgset_file))
#> [1] TRUE TRUEsva_targets_file and sva_rgset_file point to the Step 1 outputs:
phenoLC.csv, containing the phenotype table and cell-composition estimates,
and RGSet.RData, containing the filtered methylation object.
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
)
#> 3 surrogate variables explain 91.17398 % of
#> data variation
class(svaPipelineResult$savedFiles)
#> [1] "dnaEPICO_svaEnmix_paths"
names(svaPipelineResult$savedFiles)
#> [1] "svaRData" "svaCSV" "phenoWithSva" "dataDir" "rDir"
basename(unlist(svaPipelineResult$savedFiles, use.names = FALSE))
#> [1] "svaMatrix.RData" "svaMatrix.csv" "phenoLC.csv" "svaEnmix"
#> [5] "svaEnmix"The returned savedFiles object records the main paths written by this step.
names(svaPipelineResult$savedFiles) identifies the output groups: svaRData
is the serialized surrogate-variable matrix, svaCSV is the CSV version of
that matrix, phenoWithSva is the same phenoLC.csv after the SVA columns are
appended, and dataDir / rDir are the parent directories used for the saved
outputs. The basename(...) call shows the file names generated from those
paths. The phenotype is validated, written to a temporary file, read back, and
then replaced with rollback protection. Sample identifiers, row order,
cell-composition columns, and all other existing columns must remain unchanged.
Existing PC names are rejected before replacement.
Arguments:
phenoFile and rgsetData identify the phenotype table and saved RGSet
produced by the earlier preprocessing step.SampleID, SentrixIDColumn, and SentrixPositionColumn specify how
samples and array positions are identified in the phenotype table.ctrlSvaPercVar = 0.90 keeps enough control-derived surrogate variables to
explain 90% of the control-probe variance.ctrlSvaFlag = 1 enables the ENmix control-based SVA workflow.dataBaseDir and rBaseDir define where file outputs are written when
saveOutputs = TRUE.display = FALSE and verbose = FALSE keep the console output quiet, while
logs = TRUE writes the log files used by the report.preprocessingPheno() is the step that usually prepares the largest number of
pipeline-ready files because it writes timepoint-specific tables, combined
longitudinal tables, and the Clock Foundation export inputs.
The next chunk shows how to recreate the temporary phenotype and matrix files used in this example. The helper writes a phenotype table plus aligned beta, M-value, and copy-number objects to a temporary directory.
pheno_inputs <- dnaEPICO:::examplePreprocessingPhenoStateDnaEpico()
names(pheno_inputs)
#> [1] "tempDir" "pheno" "phenoPath"
#> [4] "betaPath" "mPath" "cnPath"
#> [7] "metricsData" "timepointData" "combinedData"
#> [10] "clockFoundation" "preprocessingData"
pheno_inputs$tempDir
#> [1] "/tmp/RtmpBnguYi/dnaEPICO-preprocessingPheno-example-1259cb51eabc6d"
basename(c(
pheno_inputs$phenoPath,
pheno_inputs$betaPath,
pheno_inputs$mPath,
pheno_inputs$cnPath
))
#> [1] "phenoLC.csv" "beta.RData" "m.RData" "cn.RData"pheno_inputs is a list that bundles the temporary files and the aligned
in-memory objects used in this stage. names(pheno_inputs) shows that it
contains the temporary directory, the phenotype table, the saved beta, M-value,
and copy-number files, and precomputed helper objects such as metricsData,
timepointData, combinedData, and clockFoundation.
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)
#> [1] "timepoints" "combinedPheno"
#> [3] "combinedPhenoMethylation" "methylationScale"
#> [5] "methylationObjectPrefix" "betaCSV"
#> [7] "betaZIP" "phenoCF"
#> [9] "combinedPhenoB"
names(phenoPipelineResult$savedFiles$timepoints)
#> [1] "1" "2"
basename(unlist(phenoPipelineResult$savedFiles$timepoints[["1"]]))
#> [1] "phenoT1.csv" "betaT1.RData" "mT1.RData" "cnT1.RData"
#> [5] "phenoBT1.RData" "phenoBT1.RData"
basename(unlist(phenoPipelineResult$savedFiles[c(
"combinedPheno",
"combinedPhenoMethylation",
"betaCSV",
"phenoCF"
)]))
#> [1] "phenoT1T2.csv" "phenoBT1T2.RData" "beta.csv" "phenoCF.csv"These outputs are the files most commonly consumed by downstream modeling
functions. names(phenoPipelineResult$savedFiles) shows the high-level output
groups returned by the function.
names(phenoPipelineResult$savedFiles$timepoints)
lists the available timepoint-specific subsets. The basename(...) calls then
show examples of the concrete files written for one timepoint and for the
combined outputs. In this structure, combinedPhenoMethylation is the merged
phenotype-plus-methylation object that becomes the direct input to the GLM and
LME functions. It defaults to beta values, but can be written as M-values or
copy-number values when methylationScale is changed. The betaCSV, betaZIP,
and phenoCF outputs remain beta-based Clock Foundation exports.
Arguments:
phenoFile, betaPath, mPath, and cnPath point to the phenotype table
and the aligned methylation matrices from the previous workflow stages.SampleID = "Sample_Name" defines the sample key used during alignment.timeVar = "Timepoint", timepoints = "1,2", and
combineTimepoints = "1,2" determine the timepoint-specific outputs and the
combined longitudinal object.methylationScale = "beta" selects the merged modeling scale. Use "m" for
M-values or "cn" for copy-number values. Clock Foundation exports continue
to use beta values.outputPheno, outputRData, outputRDataMerge, and outputDir define the
directories used for phenotype exports, metric objects, merged objects, and
Clock Foundation inputs.saveOutputs = TRUE is what makes this step useful in a file-based pipeline,
because later modeling steps consume these written files.The modeling function also separates in-memory results from file outputs. The next example runs a small GLM analysis and prints the names of the files written to disk.
The next chunk shows how to recreate the temporary merged phenotype-plus-beta input file used by the GLM example.
glm_inputs <- dnaEPICO:::exampleMethylationGLMStateDnaEpico()
names(glm_inputs)
#> [1] "tempDir" "inputPath" "preparedData" "modelResults"
#> [5] "modelSummaries" "annotationData"
glm_inputs$tempDir
#> [1] "/tmp/RtmpBnguYi/dnaEPICO-glm-example-1259cb49f03d73"
basename(glm_inputs$inputPath)
#> [1] "phenoBT1.RData"
file.exists(glm_inputs$inputPath)
#> [1] TRUEglm_inputs is a list returned by the GLM example helper. names(glm_inputs)
shows that it contains the temporary directory, the saved merged
phenotype-plus-beta input file (inputPath), and the corresponding in-memory
objects produced during helper construction: preparedData, modelResults, and
modelSummaries.
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)
#> [1] "dnaEPICO_methylationGLM_paths"
names(glmPipelineResult$savedFiles)
#> [1] "modelFiles" "summaryFiles"
#> [3] "summaryTxtFiles" "significantCpGFiles"
#> [5] "annotatedGLM" "annotatedGLMText"
#> [7] "annotatedGLMReportMetadata" "annotatedGLMDictionary"
#> [9] "annotatedGLMMetadata" "vennDSheets"The returned savedFiles object records the groups of outputs written by the
GLM step. names(glmPipelineResult$savedFiles) includes compact phenotype
summary files, significant-CpG exports, summary text files, and the annotated
results workbook. Inspect this object to identify the written modeling outputs
and their groups.
Arguments:
inputPheno points to the merged phenotype-plus-beta object generated by
preprocessingPheno().phenotypes lists the outcome variables that will be modeled one at a time.covariates lists the adjustment variables included in every model.factorVars identifies the categorical variables that should be treated as
factors before model fitting.scaleVars identifies numeric fixed-effect variables to centre and divide by
their sample standard deviations before fitting.<Phenotype>_Model.Message; model attempts without
a returned p-value are counted in workbook metadata and are not added to the
annotated result table.cpgLimit = 2 keeps the example fast by fitting only two CpG models.nCores = 1 keeps the example deterministic and lightweight. In production,
nCores is a maximum; the effective worker count is also limited by the
workload, detected CPUs, and available memory.annotationPackage and annotationCols control which array annotation is
merged into the final summary tables.gencodeHub = TRUE retrieves the package-managed gene resource from
AnnotationHub and adds direct gene-body and nearest strand-specific TSS
annotations. This option requires a GRCh38 array annotation; the GENCODE
release is read from the resource metadata and retained in the output
column names, workbook metadata, and dictionary.outputRData, outputPlots, significantCpGDir, summaryTxtDir, and
annotatedGLMOut define where compact phenotype summaries, plots, text
summaries, and the
annotated GLM workbook are written.resumeFromSummary = TRUE reuses a complete phenotype summary only when it
was generated from the same input file and model configuration. A phenotype
without a complete summary restarts from its first CpG.reportAssetsDir writes the compressed report table directly into the report
project, leaving the annotated workbook directory free of report sidecars.saveOutputs = TRUE enables the file-based pipeline behaviour expected by
Makefile and HPC use.The longitudinal modeling step follows the same pattern as the GLM function, but the saved outputs now correspond to mixed-effects models, longitudinal summary tables, and interaction-specific result files.
The next chunk shows how to recreate the temporary combined longitudinal phenotype-plus-beta input file used by the LME example.
lme_inputs <- dnaEPICO:::exampleMethylationLMEStateDnaEpico()
names(lme_inputs)
#> [1] "tempDir" "inputPath" "preparedData" "modelResults"
#> [5] "modelSummaries" "annotationData"
lme_inputs$tempDir
#> [1] "/tmp/RtmpBnguYi/dnaEPICO-lme-example-1259cb76ae68b"
basename(lme_inputs$inputPath)
#> [1] "phenoBT1T2.RData"
file.exists(lme_inputs$inputPath)
#> [1] TRUElme_inputs plays the same role for the longitudinal example. The helper
returns a list with the temporary directory, the saved combined longitudinal
input file (inputPath), and the in-memory objects used to build that example:
preparedData, modelResults, and modelSummaries.
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)
#> [1] "dnaEPICO_methylationLME_paths"
names(lmePipelineResult$savedFiles)
#> [1] "modelFiles" "summaryFiles"
#> [3] "summaryTxtFiles" "significantInteractionFiles"
#> [5] "annotatedLME" "annotatedLMEText"
#> [7] "annotatedLMEReportMetadata" "annotatedLMEDictionary"
#> [9] "annotatedLMEMetadata" "vennDSheets"The returned savedFiles object records the groups of outputs written by the
LME step. names(lmePipelineResult$savedFiles) includes compact phenotype
summary files, summary text files, significant interaction exports, and the
annotated longitudinal results table. Inspect this object to identify the
written longitudinal modeling outputs and their groups.
Arguments:
inputPheno points to the combined longitudinal phenotype-plus-beta object.personVar defines the participant identifier used for the random effect.SampleID identifies the sample-level column used to derive personVar when
that participant column is absent. Automatic derivation removes a terminal
A or B; other identifier formats require an explicit participant column.timeVar defines the repeated-measures time variable used for longitudinal
summaries and preprocessing checks.phenotypes, covariates, and optional interactionTerm define the
fixed-effect portion of the mixed model.omnibusTest = TRUE adds one joint fixed-effect F test per CpG. With an
interaction, all estimable phenotype-by-interaction coefficients are tested
together; without an interaction, the phenotype main effect is tested.omnibusDdf selects "Satterthwaite" or "Kenward-Roger" denominator
degrees of freedom. Kenward-Roger testing requires pbkrtest and is more
computationally intensive for an EWAS.<Phenotype>_Model.Message; model attempts without a returned coefficient
or omnibus p-value are counted in workbook metadata and are not added to
the annotated result table.cpgLimit = 2 keeps the example short by fitting only two CpG models.nCores = 1 keeps the example lightweight. In production, nCores is a
maximum; the effective worker count is also limited by the workload,
detected CPUs, and available memory.saveSignificantInteractions, significantInteractionDir, and
significantInteractionPval control the export of interaction-specific
results.summaryTxtDir defines where plain-text model summaries are written.annotationPackage, annotationCols, and annotatedLMEOut define the
annotation resources and output location for the final longitudinal summary
workbook.gencodeHub = TRUE applies the same release-aware AnnotationHub annotation
to lmerTest/lme4 and nlme results when the array coordinates use GRCh38.reportAssetsDir writes the compressed longitudinal report table directly
into the report project rather than beside the workbook.resumeFromSummary = TRUE reuses a complete phenotype summary only when it
was generated from the same input file and model configuration. A phenotype
without a complete summary restarts from its first CpG.saveOutputs = TRUE enables the file-based longitudinal workflow expected by
Makefile and HPC use.For a full-cohort longitudinal interaction analysis, the following structure fits a participant-level random intercept and returns one joint interaction p-value per CpG in addition to the component coefficient p-values:
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"
)This model uses the fixed effects Timepoint * Profession together with
(1 | Participant). A significant omnibus result indicates that at least one
timepoint contrast differs by profession; the component coefficients are used
to identify its direction and location.
After the file-based examples above have been run, dnamReport() can assemble
the phenotype table, ENmix control figures, quality-control figures, batch
effect figures, metrics figures, model annotation tables, and logs into one
website report. The example below uses the concrete output paths created by the
previous chunks.
The preprocessing and SVA examples share preprocessing_inputs$tempDir, while
the GLM and LME examples write model outputs into their own temporary
directories. Because dnamReport() accepts one path per tab, those outputs can
be passed directly.
When running this section from an R session, dnamReport() still needs the
Quarto command line interface because the function renders the .qmd files into
the final website. Quarto is not required to install or load dnaEPICO, or to
run its preprocessing and statistical-modeling functions. Check whether R can
find Quarto with:
Sys.which("quarto")
system2("quarto", "--version")If Sys.which("quarto") returns an empty string, install the Quarto command
line interface from the
Quarto get-started page before
rendering the report. The R package quarto provides R helper functions for an
existing Quarto installation; it does not replace the Quarto command line
interface used by dnamReport().
If Quarto is installed but not on PATH, set the executable path before calling
dnamReport():
Sys.setenv(QUARTO_BIN = "/full/path/to/quarto")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$outputFileFor each CpG, methylationGLM() fits a regression model of the form:
\[ \text{Methylation}_{CpG_i} = \beta_0 + \beta_1 \cdot \text{Phenotype} + \sum_{k = 1}^{K} \beta_{k + 1} \cdot \text{Covariate}_k + \varepsilon \]
where \(\text{Methylation}_{CpG_i}\) is the methylation value at CpG \(i\), \(\beta_0\) is the intercept, the phenotype term is the variable of interest, and the remaining fixed effects represent the listed covariates.
Setting omnibusTest = TRUE adds one joint Wald F test per CpG using
car::linearHypothesis(). Without an interaction, the complete phenotype term
is tested. With an interaction, the complete phenotype-by-interaction term is
tested. One-degree-of-freedom terms, including a numeric phenotype entered as
one linear term, are also tested and reproduce the corresponding coefficient
p-value. Omnibus p-values are adjusted across valid CpGs within each phenotype
and tested term using padjmethod.
If no polygenic risk score is included, a model can be read schematically as:
\[ \begin{aligned} \text{Methylation}_{CpG_i} =\;& \beta_0 + \beta_1 \cdot \text{Phenotype} + \beta_2 \cdot \text{Age} + \beta_3 \cdot \text{Sex} + \beta_4 \cdot \text{Ethnicity} \\ &+ \varepsilon \end{aligned} \]
When prsMap is supplied, a phenotype-specific PRS term is appended only for
the matching phenotype. For example, the mapping
"Pheno1:PRS_1,Pheno2:PRS_2" means that:
Pheno1 include PRS_1Pheno2 include PRS_2For a phenotype mapped to a PRS, the model becomes:
\[ \text{Methylation}_{CpG_i} = \beta_0 + \beta_1 \cdot \text{Phenotype} + \sum_{k = 1}^{K} \beta_{k + 1} \cdot \text{Covariate}_k + \beta_{\text{PRS}} \cdot \text{PRS} + \varepsilon \]
If interactionTerm is supplied, the fixed effects include its main effect and
its interaction with the phenotype.
The longitudinal function methylationLME() follows the same file-based
pattern, but uses a mixed-effects model with a participant-level random
intercept. In schematic form:
\[ \begin{aligned} \text{Methylation}_{CpG_i} =\;& \beta_0 + \beta_1 \cdot \text{Phenotype} + \beta_2 \cdot \text{Interaction} + \beta_3 \cdot (\text{Phenotype} \times \text{Interaction}) + \sum_{k = 1}^{K} \beta_{k + 3} \cdot \text{Covariate}_k + \beta_{\text{PRS}} \cdot \text{PRS} + b_{\text{person}} + \varepsilon \end{aligned} \]
where \(b_{\text{person}}\) is the subject-specific random intercept. As in the
GLM workflow, the PRS term is included only when the current phenotype is
matched in prsMap.
Common arguments in this longitudinal stage are:
inputPheno for the combined longitudinal phenotype-plus-beta objectpersonVar for the participant identifier used in the random effectSampleID for the configured sample identifier used only when personVar
must be derivedtimeVar for repeated-measures summaries and preprocessing checksphenotypes, covariates, and optional interactionTerm for the
fixed-effect
structureprsMap for phenotype-specific PRS termscpgLimit and nCores for computational controlannotationPackage and annotationCols for annotated summariesRunning the exported workflow requires GNU Make 4.3 or later. When
extractMake() creates the Makefile, it records the absolute path to the
Rscript executable from the current R installation in the RSCRIPT Make
variable. Therefore, Rscript does not need to be on the system path. Report
targets also require Quarto. GNU Make and Quarto are not required to install or
load dnaEPICO; they support the exported workflow and report-rendering
features, respectively. The exported rules use relative paths and support
project directories containing spaces. A minimal layout is:
my_project/
Makefile
metadata/
pheno_model1.csv
data/
preprocessingMinfiEwasWater/
idats/
<IDAT files>The exported template should then be edited so that at least the model name, phenotype path, and directory variables match the project. A minimal example is:
MODEL = model1
PHENO_FILE = metadata/pheno_model1.csv
LOGS_DIR = logs
DATA_DIR = data
RDATA_DIR = rData
RESULTS_DIR = results
FIGURES_DIR = figuresA typical command sequence is:
make step1 MODEL=model1
make step2 MODEL=model1
make step3 MODEL=model1
make f3 MODEL=model1
make f4 MODEL=model1
make f3lme MODEL=model1
make all MODEL=model1
make status MODEL=model1
make clean MODEL=model1To use another R installation, override the recorded executable when invoking Make. This is useful after loading an R module on a high-performance computing system or after moving the Makefile to another computer:
make f4 MODEL=model1 RSCRIPT=/path/to/RscriptAlternatively, regenerate the Makefile with extractMake() from the R
installation that will run the workflow.
The targets run these workflow stages:
make step1 MODEL=model1 runs preprocessing only.make step2 MODEL=model1 runs the SVA step using the files written by Step 1.make step3 MODEL=model1 runs phenotype preparation after Step 2 appends
surrogate-variable columns to phenoLC.csv.make f3 MODEL=model1 runs Steps 1-3 and generates a preprocessing report.make f4 MODEL=model1 runs Steps 1-4 and generates a report with
preprocessing and GLM sections.make f3lme MODEL=model1 runs Steps 1-3 and LME, then generates a report with
preprocessing and LME sections.make all MODEL=model1 runs Steps 1-5 and generates the complete report.make status MODEL=model1 reports available outputs.make clean MODEL=model1 removes outputs for the selected model while
preserving shared raw inputs.If several models are declared in MODELS, the grouped multi-model targets can
be used instead:
make f3_models
make f4_models
make f3lme_models
make modelsThese targets dispatch every model listed in MODELS:
make f3_models runs the f3 route for each declared model.make f4_models runs the f4 route for each declared model.make f3lme_models runs the f3lme route for each declared model.make models runs the full all pipeline for each declared model.Model-level Venn output is disabled when its Makefile values are NULL. Set
VENND_GLM_PHENOTYPES or VENND_LME_PHENOTYPES to comma-separated phenotype
names to compare all coefficient p-value columns generated for those
phenotypes. A factor or interaction phenotype therefore contributes every
applicable coefficient column; the interaction variable is not repeated in the
selection. Optional VENND_GLM_LABELS and VENND_LME_LABELS values are applied
in that expanded coefficient order and retain their supplied letter case.
Omnibus p-values are selected independently with
VENND_GLM_OMNIBUS_PHENOTYPES or VENND_LME_OMNIBUS_PHENOTYPES, with optional
positional labels in the corresponding *_OMNIBUS_LABELS value. These options
require the relevant omnibus test to be enabled. Each requested analysis writes
nominal and genome-wide UCSC and release-labelled GENCODE figures in
figures/<model>/methylationGLM or figures/<model>/methylationLME. Nominal,
suggestive, and genome-wide membership tables are added to the annotated
workbook between metadata and the final dictionary sheet. Venn progress and
errors are retained in the corresponding model-analysis log. When four or more
sets are resolved, a ranked intersection figure accompanies each Venn diagram
to keep the overlap counts readable.
Manhattan plots and Venn membership use the same 9e-8 genome-wide
significance threshold. This threshold controls figure classifications and
membership tables; it does not alter model p-values or their FDR adjustment.
Set GENCODE_HUB = TRUE for a GRCh38 model when GENCODE output is required.
The R functions retrieve the internal AnnotationHub resource ID, validate its
assembly and gene-level schema, and derive the release label from the resource
metadata. Users do not provide a separate release value or local annotation
file. AnnotationHub uses its normal local cache after the first retrieval.
GENCODE_HUB = FALSE performs no GENCODE retrieval. Models using an hg19 array
annotation, including the 450K hg19 example, must keep this value FALSE.
Status and cleanup targets are available for routine use:
make status MODEL=model1
make clean MODEL=model1
make clean_modelsstatus reports available outputs, clean removes one model, and
clean_models removes every model listed in MODELS. Cleanup removes the
corresponding directories under data, rData, reports, logs, figures,
and results, and removes empty output roots.
Each step writes outputs consumed by downstream targets. Step 2 creates
data/<model>/svaEnmix/.sva_complete only after the updated phenoLC.csv has
been validated and installed. Downstream targets depend on both files, so Make
reruns SVA when Step 1 replaces phenoLC.csv and skips SVA when the completed
inputs are unchanged.
The exported Makefile contains a report rule whose output is
$(STEP6)/$(MODEL)/docs/index.html. With the default directory settings, this
becomes reports/<model>/docs/index.html. Users usually do not need to call
that file target directly: the grouped targets above depend on it, so Make
refreshes the report automatically after the required pipeline outputs are
available.
For example, the f3 route builds the report from the Data, Quality Control,
Batch Effect, Metrics, Report, and Logs inputs; Quality Control contains both
ENmix controls and methylation quality-control outputs. The f4 route adds only
the GLM page and GLM log, the f3lme route adds only the LME page and LME log,
and all includes both modeling sections. Each report contains the model
sections defined by its route.
The report uses a responsive three-panel dashboard. The dark left panel holds the persistent tab navigation and the dnaEPICORM logo, the centre panel uses the available page width for the selected table or figure, and the right panel shows section notes that update for the selected figure or workbook sheet. On smaller displays the navigation becomes a compact top menu and the notes move below the main content.
Figure-bearing tabs display one selected image at a time with previous/next, mouse-controlled zoom, scrolling, one card expansion control, and a download link outside the image region. Workbook tabs provide a sheet selector without loading every sheet into one visible panel. Package-generated plots adapt point density, labels, dimensions, and longitudinal traces to the number of samples or CpGs, while minfi, limma, and ENmix plots retain their established plotting methods.
The generated figures contain axis labels, legends, thresholds, and statistical data labels where required, but no embedded title or descriptive caption. The corresponding report card supplies the figure title, and explicit file names identify the analysis, variable, diagnostic, annotation source, and version as applicable. SVA technical-factor association figures are omitted when no association can be estimated, with the condition retained in the analysis log.
The report target reads the same types of outputs generated in the examples
above, organised under the Makefile project layout: the phenotype table in
data/, quality-control and metric figures in figures/, the detection
P-value object in rData/, model annotation tables in data/, and workflow
logs in logs/.
The generated site is opened from the report target path:
$(STEP6)/$(MODEL)/docs/index.htmlInternally, the Makefile calls dnamReport() with one argument per tab, using
the project paths defined by the Makefile variables. The report target requires
the Quarto command line interface in the same environment where make is run.
If Quarto is not already available, install the Quarto command line interface
from the Quarto get-started page, or
use the installation method provided by the local computing environment. Before
launching the report target, check:
quarto --versionIf Quarto is installed outside PATH, export its executable path before running
make:
export QUARTO_BIN=/full/path/to/quartoThe full Makefile exported by extractMake() is shown below. Adapt this
template
for the project before running the pipeline in the target analysis environment.
# ===============================================
# USER CONFIGURATION
# ===============================================
MAKEFLAGS += --output-sync
# Rscript from the R installation that exported this Makefile.
# Override this value when using another R installation.
RSCRIPT ?= /home/biocbuild/bbs-3.24-bioc/R/bin/Rscript
# ===============================================
# MODEL SELECTION
# ===============================================
# Single model
MODEL ?= model1
# Multiple models for parallel runs
MODELS = modelA modelB modelC
# ===============================================
# PER-MODEL OVERRIDES
# ===============================================
ifeq ($(MODEL), modelA)
PHENO_FILE = $(DATA_DIR)/$(STEP1)/phenoA.csv
else ifeq ($(MODEL), modelB)
PHENO_FILE = $(DATA_DIR)/$(STEP1)/phenoB.csv
else ifeq ($(MODEL), modelC)
PHENO_FILE = $(DATA_DIR)/$(STEP1)/phenoC.csv
else
PHENO_FILE = $(DATA_DIR)/$(STEP1)/pheno.csv
endif
# ===============================================
# DIRECTORIES
# ===============================================
LOGS_DIR = logs
DATA_DIR = data
RDATA_DIR = rData
RESULTS_DIR = results
FIGURES_DIR = figures
STEP1 = preprocessingMinfiEwasWater
STEP2 = svaEnmix
STEP3 = preprocessingPheno
STEP4 = methylationGLM
STEP5 = methylationLME
STEP6 = reports
METRICS_DIR = metrics
IDAT_DIR = idats
OBJ_DIR = objects
MERGE_DIR = mergeData
MODEL_DIR = models
CPG_DIR = cpgs
SUMMARY_DIR = summary
ENMIX_DIR = enmix
SVA_DIR = sva
QC_DIR = qc
# Optional R library path for local or HPC environments
R_DIR = NULL
# ===============================================
# GLOBAL PARAMETERS
# ===============================================
# STEP 1-5 PARAMETERS
SEED = 123
SEP_TYPE = NULL
SAMPLE_ID = Sample_Name
N_SAMPLES = 30
ARRAY_TYPE = IlluminaHumanMethylationEPICv2
ANNOTATION_VERSION = 20a1.hg38
IDAT_FORCE = FALSE
TIFF_WIDTH = 2000
TIFF_HEIGHT = 1000
TIFF_RES = 150
SEX_COLUMN = Gender
REMOVE_SEX_MISMATCH = FALSE
SENTRIX_ID_COLUMN = Sentrix_ID
SENTRIX_POSITION_COLUMN = Sentrix_Position
BASENAME_COLUMN = Basename
TIME_VAR = Timepoint
METHYLATION_SCALE = beta
PHENO_ORDER = $(SAMPLE_ID);$(TIME_VAR);$(SEX_COLUMN);PredSex;$(BASENAME_COLUMN);$(SENTRIX_ID_COLUMN);$(SENTRIX_POSITION_COLUMN)
# STEP 4-5 PARAMETERS
CPG_PREFIX = cg
CPG_LIMIT = NA
PRS_MAP = NULL
SUMMARY_PVAL = NA
N_CORES = 64
SAVE_TXT_SUMMARIES = TRUE
RESUME_FROM_SUMMARY = TRUE
CHUNK_SIZE = 10000
FDR_THRESHOLD = 0.05
PADJ_METHOD = fdr
ANNOTATION_PACKAGE = IlluminaHumanMethylationEPICv2anno.20a1.hg38
ANNOTATION_COLS = Name,chr,pos,UCSC_RefGene_Group,UCSC_RefGene_Name,Relation_to_Island,GencodeV41_Group
GENCODE_HUB = FALSE
# ===============================================
# STEP 1 PARAMETERS
# ===============================================
QC_CUTOFF = 10.5
DET_PTYPE = m+u
DET_PTHRESHOLD = 0.05
PVAL_THRESHOLD = 0.01
CHR_TO_REMOVE = chrX,chrY
SNPS_TO_REMOVE = SBE,CpG
PROBE_EXCLUSION_FILE = $(DATA_DIR)/$(STEP1)/12864_2024_10027_MOESM8_ESM.csv
PROBE_EXCLUSION_ID_COLUMN = NULL
USE_EPICV2_MANIFEST = FALSE
EPICV2_MANIFEST_CH_WGBS_EVIDENCE = TRUE
EPICV2_MANIFEST_CH_BLAT = TRUE
EPICV2_MANIFEST_MISSING_POS = TRUE
EPICV2_MANIFEST_MISMATCH_POS = FALSE
MAF_THRESHOLD = 0.1
PLOT_GROUP_VAR = TreatmentGroup
LC_REF = salivaEPIC
# ===============================================
# STEP 2 PARAMETERS
# ===============================================
CTRL_SVA_PERC_VAR = 0.90
CTRL_SVA_FLAG = 1
# ===============================================
# STEP 3 PARAMETERS
# ===============================================
TIMEPOINTS = 1,2
COMBINE_TIMEPOINTS = 1,2
# ===============================================
# STEP 4 PARAMETERS
# ===============================================
PHENOTYPES_GLM = TreatmentGroup
COVARIATES_GLM = Sex
FACTOR_VARS_GLM = Sex,TreatmentGroup
SCALE_VARS_GLM = NULL
GLM_LIBS = glm2
INTERACTION_GLM = NULL
GLM_OMNIBUS_TEST = FALSE
SUMMARY_RESIDUAL_SD = TRUE
SAVE_SIGNIFICANT_CPGS = TRUE
SIGNIFICANT_CPG_PVAL = 0.1
VENND_GLM_PHENOTYPES = NULL
VENND_GLM_LABELS = NULL
VENND_GLM_OMNIBUS_PHENOTYPES = NULL
VENND_GLM_OMNIBUS_LABELS = NULL
# ==============================================
# STEP 5 PARAMETERS
# ==============================================
PHENOTYPES_LME = TreatmentGroup
COVARIATES_LME = Sex
FACTOR_VARS_LME = Sex,TreatmentGroup
SCALE_VARS_LME = NULL
PERSON_VAR = person
LME_LIBS = lme4,lmerTest
LME_CORRELATION_STRUCTURE = none
LME_CORRELATION_VAR = NULL
INTERACTION_LME = NULL
LME_OMNIBUS_TEST = FALSE
LME_OMNIBUS_DDF = Satterthwaite
SAVE_SIGNIFICANT_INTERACTIONS = TRUE
SIGNIFICANT_INTERACTION_PVAL = 0.1
VENND_LME_PHENOTYPES = NULL
VENND_LME_LABELS = NULL
VENND_LME_OMNIBUS_PHENOTYPES = NULL
VENND_LME_OMNIBUS_LABELS = NULL
# ===============================================
# LOAD dnaEPICO PIPELINE RULES
# ===============================================
DNAPIPE_EMPTY :=
DNAPIPE_SPACE := $(DNAPIPE_EMPTY) $(DNAPIPE_EMPTY)
DNAPIPE_MK_RAW := $(shell "$(RSCRIPT)" -e "cat(system.file('extdata','make','Makefile.rules.pipeline',package='dnaEPICO'))")
DNAPIPE_MK := $(subst $(DNAPIPE_SPACE),\ ,$(DNAPIPE_MK_RAW))
include $(DNAPIPE_MK)
The exported Makefile is organised into sections so that project-specific choices can be edited without changing the pipeline rules themselves.
MAKEFLAGS += --output-sync keeps parallel Make output grouped by recipe,
which makes logs easier to read.RSCRIPT identifies the R installation used by the workflow. Its default is
recorded by extractMake() and can be overridden when Make is invoked.MODEL defines a single pipeline configuration to run.MODELS defines a set of configurations that can be launched in parallel by
the grouped Make targets.In single-model use, commands such as make f4 MODEL=model1 run one workflow.
In multi-model use, commands such as make f4_models evaluate the same template
for each name listed in MODELS.
PHENO_FILE is reassigned according to the selected MODEL.This block is useful when different cohorts, case definitions, or phenotype encodings should reuse the same analysis pipeline without duplicating the full Makefile.
LOGS_DIR, DATA_DIR, RDATA_DIR, RESULTS_DIR, and FIGURES_DIR define
the top-level folder structure.STEP1 to STEP6 define the canonical names of the major pipeline stages.METRICS_DIR, IDAT_DIR, OBJ_DIR, MERGE_DIR, MODEL_DIR, CPG_DIR,
SUMMARY_DIR, ENMIX_DIR, SVA_DIR, and QC_DIR define the subdirectory
names used by the rules file.R_DIR optionally points to a custom R library path on shared systems. Its
default value is NULL, which uses the default library path.These variables control where files are written and how paths are composed across all stages of the pipeline.
These parameters are shared across multiple workflow steps.
CPG_PREFIX: prefix used to identify CpG columns.CPG_LIMIT: optional cap on the number of CpGs processed; NA uses all.PRS_MAP: optional mapping from phenotype names to PRS variables; NULL
omits PRS terms from the models.SUMMARY_PVAL: optional p-value cutoff for summary tables; NA keeps all
summary rows.N_CORES: number of cores used within each R job.SAVE_TXT_SUMMARIES: whether plain-text summaries should be written.RESUME_FROM_SUMMARY: whether complete phenotype summaries generated from
the same input files and model configurations may be reused after a model
stage stops.CHUNK_SIZE: number of CpGs processed per chunk; NULL lets the R helpers
choose the chunk size automatically. Chunks bound in-memory computation;
they are not saved as restart checkpoints.FDR_THRESHOLD: false-discovery-rate threshold used for filtering.PADJ_METHOD: multiple-testing correction method.ANNOTATION_PACKAGE: array annotation package used for annotated summaries.ANNOTATION_COLS: annotation fields to merge into output tables.GENCODE_HUB: TRUE adds the package-managed release-aware GENCODE
AnnotationHub resource to GRCh38 GLM and LME workbooks; FALSE leaves the
array annotation unchanged.QC_CUTOFF: sample-quality threshold applied during preprocessing.DET_PTYPE: detection p-value method.DET_PTHRESHOLD: detection p-value cutoff.PVAL_THRESHOLD: probe-level p-value filter applied after preprocessing.CHR_TO_REMOVE: chromosomes to exclude from downstream matrices.SNPS_TO_REMOVE: SNP-related probe categories to remove.PROBE_EXCLUSION_FILE: semicolon-separated probe-exclusion reference files.PROBE_EXCLUSION_ID_COLUMN: optional column or semicolon-separated columns
containing probe IDs; NULL auto-detects each reference file.USE_EPICV2_MANIFEST: whether to also remove EPICv2 probes flagged in the
Peters et al. expanded manifest from AnnotationHub resource AH116484.EPICV2_MANIFEST_CH_WGBS_EVIDENCE, EPICV2_MANIFEST_CH_BLAT,
EPICV2_MANIFEST_MISSING_POS, and EPICV2_MANIFEST_MISMATCH_POS: manifest
flags included in the EPICv2 probe-exclusion set.MAF_THRESHOLD: minor-allele-frequency threshold used during SNP filtering.PLOT_GROUP_VAR: phenotype variable used for grouped QC plots.LC_REF: reference panel used for cell-type estimation.These values control the main preprocessing and filtering choices in
preprocessingMinfiEwasWater().
CTRL_SVA_PERC_VAR: target proportion of control-probe variance explained by
the retained surrogate variables.CTRL_SVA_FLAG: enables the control-based ENmix SVA workflow.These values determine how svaEnmix() estimates and retains surrogate
variables. The completed step updates the Step 1 phenoLC.csv in place and
records successful completion with .sva_complete; it does not create a second
phenotype file.
TIMEPOINTS: comma-separated timepoints exported separately.COMBINE_TIMEPOINTS: comma-separated timepoints combined into the
longitudinal object, or NULL to skip the combined object.These values determine how preprocessingPheno() splits and recombines the
phenotype-methylation data. The shared METHYLATION_SCALE value selects beta,
M-value, or copy-number inputs and is interpreted case-insensitively.
PHENOTYPES_GLM: comma-separated phenotype variables that will be modelled
one at a time by methylationGLM().COVARIATES_GLM: comma-separated adjustment variables included in each GLM.FACTOR_VARS_GLM: GLM variables treated as categorical.SCALE_VARS_GLM: numeric GLM fixed-effect variables centred and divided by
their sample standard deviations; NULL retains their original scales.GLM_LIBS: GLM implementation used by methylationGLM().INTERACTION_GLM: optional interaction term added to the fixed-effect part
of the GLM; NULL omits the interaction.GLM_OMNIBUS_TEST: whether one joint Wald F test is reported for the
phenotype or phenotype-by-interaction term at each CpG.SUMMARY_RESIDUAL_SD: whether residual standard deviations are reported in
summaries.SAVE_SIGNIFICANT_CPGS: whether tables of significant CpGs are exported.SIGNIFICANT_CPG_PVAL: p-value threshold used for those significant-CpG
exports.VENND_GLM_PHENOTYPES: phenotype names whose coefficient p-value columns are
compared in GLM Venn outputs; NULL disables coefficient Venn output.VENND_GLM_LABELS: optional labels applied in the expanded coefficient order.VENND_GLM_OMNIBUS_PHENOTYPES: phenotype names whose omnibus p-values are
compared; NULL disables omnibus Venn output.VENND_GLM_OMNIBUS_LABELS: optional labels for the selected omnibus sets.Model.Message column and aggregate availability counts
are written to metadata.PHENOTYPES_LME: comma-separated phenotype variables that will be modelled
one at a time by methylationLME().COVARIATES_LME: comma-separated adjustment variables included in each LME.FACTOR_VARS_LME: LME variables treated as categorical.SCALE_VARS_LME: numeric LME fixed-effect variables centred and divided by
their sample standard deviations; NULL retains their original scales.PERSON_VAR: participant identifier used as the random effect grouping
factor.LME_LIBS: mixed-effects libraries used by methylationLME().LME_CORRELATION_STRUCTURE: residual correlation structure for the nlme
LME backend; use none, AR1, or CAR1.LME_CORRELATION_VAR: variable used to order repeated observations within
participant for AR1 or CAR1; required for those structures.INTERACTION_LME: optional interaction term used in the longitudinal model;
NULL omits the interaction.LME_OMNIBUS_TEST: whether lmerTest/lme4 models calculate one joint test for
the selected phenotype main effect or phenotype-by-interaction term.LME_OMNIBUS_DDF: denominator degrees-of-freedom method for the joint test;
use Satterthwaite or Kenward-Roger.SAVE_SIGNIFICANT_INTERACTIONS: whether significant interaction tables are
exported.SIGNIFICANT_INTERACTION_PVAL: p-value threshold for those exported
interaction results.VENND_LME_PHENOTYPES: phenotype names whose coefficient p-value columns are
compared in LME Venn outputs; NULL disables coefficient Venn output.VENND_LME_LABELS: optional labels applied in the expanded coefficient order.VENND_LME_OMNIBUS_PHENOTYPES: phenotype names whose omnibus p-values are
compared; NULL disables omnibus Venn output.VENND_LME_OMNIBUS_LABELS: optional labels for the selected omnibus sets.Model.Message column and aggregate availability counts
are written to metadata.The template uses NA for optional numeric limits or filters where the meaning
is “use all” or “keep all”, and NULL when an optional separator, variable,
path, mapping, or model term should be absent. CHUNK_SIZE uses NULL for
automatic chunking because this disables a user-specified tuning value rather
than applying a filter.
DNAPIPE_MK resolves the packaged Makefile.rules.pipeline file.include $(DNAPIPE_MK) imports the actual pipeline rules after the user
configuration has been defined.This keeps the template compact while the package maintains the rules in one location.
Use the file-based workflow when you need:
The interactive workflow supports direct inspection of returned objects; the file-based workflow supports reproducible multi-step execution.
Date the vignette was generated.
#> [1] "2026-08-18 17:33:48 EDT"Wallclock time spent generating the vignette.
#> Time difference of 2.255 minsR session information.
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.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] parallel stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.1
#> [2] IlluminaHumanMethylation450kmanifest_0.4.0
#> [3] minfi_1.59.1
#> [4] bumphunter_1.55.1
#> [5] locfit_1.5-9.12
#> [6] iterators_1.0.14
#> [7] foreach_1.5.2
#> [8] Biostrings_2.81.6
#> [9] XVector_0.53.0
#> [10] SummarizedExperiment_1.43.0
#> [11] Biobase_2.73.2
#> [12] MatrixGenerics_1.25.0
#> [13] matrixStats_1.5.0
#> [14] GenomicRanges_1.65.1
#> [15] Seqinfo_1.3.0
#> [16] IRanges_2.47.2
#> [17] S4Vectors_0.51.6
#> [18] BiocGenerics_0.59.12
#> [19] generics_0.1.4
#> [20] dnaEPICO_0.99.38
#> [21] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.6.1 BiocIO_1.23.3
#> [3] bitops_1.1-0 filelock_1.0.3
#> [5] tibble_3.3.1 preprocessCore_1.75.0
#> [7] XML_3.99-0.23 lifecycle_1.0.5
#> [9] httr2_1.3.0 Rdpack_2.6.6
#> [11] doParallel_1.0.17 lattice_0.23-1
#> [13] MASS_7.3-66 base64_2.0.2
#> [15] scrime_1.3.7 magrittr_2.0.5
#> [17] openxlsx_4.2.8.1 limma_3.69.4
#> [19] sass_0.4.10 rmarkdown_2.31
#> [21] jquerylib_0.1.4 yaml_2.3.12
#> [23] otel_0.2.0 zip_3.0.2
#> [25] doRNG_1.8.6.3 askpass_1.2.1
#> [27] minqa_1.2.8 DBI_1.3.0
#> [29] RColorBrewer_1.1-3 abind_1.4-8
#> [31] quadprog_1.5-8 purrr_1.2.2
#> [33] RCurl_1.98-1.19 rappdirs_0.3.4
#> [35] ggrepel_0.9.8 irlba_2.3.7
#> [37] rentrez_1.2.4 genefilter_1.95.0
#> [39] annotate_1.91.0 DelayedMatrixStats_1.35.0
#> [41] codetools_0.2-20 DelayedArray_0.39.5
#> [43] xml2_1.6.0 tidyselect_1.2.1
#> [45] glm2_1.2.1 farver_2.1.2
#> [47] lme4_2.0-6 beanplot_1.3.1
#> [49] BiocFileCache_3.3.0 dynamicTreeCut_1.63-1
#> [51] illuminaio_0.55.0 GenomicAlignments_1.49.1
#> [53] jsonlite_2.0.0 multtest_2.69.0
#> [55] survival_3.8-9 tools_4.6.1
#> [57] Rcpp_1.1.2 glue_1.8.1
#> [59] SparseArray_1.13.2 BiocBaseUtils_1.15.1
#> [61] xfun_0.60 dplyr_1.2.1
#> [63] HDF5Array_1.41.2 withr_3.0.3
#> [65] numDeriv_2016.8-1.1 BiocManager_1.30.27
#> [67] fastmap_1.2.0 boot_1.3-32
#> [69] rhdf5filters_1.25.4 openssl_2.4.2
#> [71] caTools_1.18.4 digest_0.6.39
#> [73] R6_2.6.1 RPMM_1.25
#> [75] gtools_3.9.5 dichromat_2.0-1
#> [77] RSQLite_3.53.3 cigarillo_1.3.1
#> [79] h5mread_1.5.0 minfiData_0.59.0
#> [81] tidyr_1.3.2 data.table_1.18.4
#> [83] rtracklayer_1.73.0 httr_1.4.8
#> [85] S4Arrays_1.13.0 pkgconfig_2.0.3
#> [87] gtable_0.3.6 blob_1.3.0
#> [89] S7_0.2.2 siggenes_1.87.0
#> [91] impute_1.87.0 htmltools_0.5.9
#> [93] bookdown_0.47 geneplotter_1.91.0
#> [95] scales_1.4.0 png_0.1-9
#> [97] reformulas_0.4.4 knitr_1.51
#> [99] tzdb_0.5.0 rjson_0.2.23
#> [101] nloptr_2.2.1 nlme_3.1-170
#> [103] curl_7.1.0 cachem_1.1.0
#> [105] rhdf5_2.57.10 KernSmooth_2.23-27
#> [107] BiocVersion_3.24.0 AnnotationDbi_1.75.2
#> [109] restfulr_0.0.17 GEOquery_2.81.28
#> [111] pillar_1.11.1 grid_4.6.1
#> [113] reshape_0.8.10 vctrs_0.7.3
#> [115] gplots_3.3.0 ENmix_1.49.3
#> [117] dbplyr_2.6.0 xtable_1.8-8
#> [119] cluster_2.1.8.3 evaluate_1.0.5
#> [121] readr_2.2.0 GenomicFeatures_1.65.0
#> [123] cli_3.6.6 compiler_4.6.1
#> [125] Rsamtools_2.29.0 rlang_1.3.0
#> [127] crayon_1.5.3 rngtools_1.5.2
#> [129] labeling_0.4.3 nor1mix_1.3-3
#> [131] mclust_6.1.3 plyr_1.8.9
#> [133] stringi_1.8.9 BiocParallel_1.47.0
#> [135] lmerTest_3.2-1 Matrix_1.7-6
#> [137] ExperimentHub_3.3.2 hms_1.1.4
#> [139] sparseMatrixStats_1.25.0 bit64_4.8.2
#> [141] ggplot2_4.0.3 Rhdf5lib_2.1.0
#> [143] KEGGREST_1.53.6 statmod_1.5.2
#> [145] AnnotationHub_4.3.2 rbibutils_2.4.1
#> [147] memoise_2.0.1 bslib_0.12.0
#> [149] bit_4.6.0Use the Bioconductor support site for
package questions. Add the dnaEPICO tag, review
previous posts, and include a
small reproducible example with session information.
Aryee, Martin J, Andrew E Jaffe, Hector Corrada Bravo, Christine Ladd-Acosta, Andrew P Feinberg, Kasper D Hansen, and Rafael A Irizarry. 2014. “Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays.” Bioinformatics 30 (10): 1363–9. https://doi.org/10.1093/bioinformatics/btu049.
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.
Fortin, Jean-Philippe, Aurélie Labbe, Mathieu Lemire, Brent W Zanke, Thomas J Hudson, Elana J Frtig, Celia MT Greenwood, and Kasper D Hansen. 2014. “Functional normalization of 450k methylation array data improves replication in large cancer studies.” Genome Biology 15 (11): 503. https://doi.org/10.1186/s13059-014-0503-2.
Fortin, Jean-Philippe, Timothy Triche Jr., and Kasper D Hansen. 2017. “Preprocessing, Normalization and Integration of the Illumina HumanMethylationEPIC Array with Minfi.” Bioinformatics 33 (4): 558–60. https://doi.org/10.1093/bioinformatics/btw691.
Maksimovic, Jovana, Lavinia Gordon, and Alicia Oshlack. 2012. “SWAN: Subset quantile Within-Array Normalization for Illumina Infinium HumanMethylation450 BeadChips.” Genome Biology 13 (6): R44. https://doi.org/10.1186/gb-2012-13-6-r44.
Marschner, Ian. 2011. “Glm2: Fitting Generalized Linear Models with Convergence Problems.” The R Journal 3 (2): 12–15. https://doi.org/10.32614/RJ-2011-012.
Murat, Kubra, Björn Grüning, Pawel W Poterlowicz, Gareth Westgate, Desmond J Tobin, and Krzysztof Poterlowicz. 2020. “Ewastools: Infinium Human Methylation BeadChip pipeline for population epigenetics integrated into Galaxy.” GigaScience 9 (5): giaa049. https://doi.org/10.1093/gigascience/giaa049.
Pidsley, Ruth, Ching Ching Y Wong, Matteo Volta, Katie Lunnon, Jonathan Mill, and Leonard C Schalkwyk. 2013. “A data-driven approach to preprocessing Illumina 450K methylation array data.” BMC Genomics 14: 293. https://doi.org/10.1186/1471-2164-14-293.
Touleimat, Nizar, and Jörg Tost. 2012. “Complete Pipeline for Infinium() Human Methylation 450K BeadChip Data Processing Using Subset Quantile Normalization for Accurate DNA Methylation Estimation.” Epigenomics 4 (3): 325–41. https://doi.org/10.2217/epi.12.21.
Triche, Timothy J, Daniel J Weisenberger, David Van Den Berg, Peter W Laird, and Kimberly D Siegmund. 2013. “Low-level processing of Illumina Infinium DNA Methylation BeadArrays.” Nucleic Acids Research 41 (7): e90. https://doi.org/10.1093/nar/gkt090.
Xu, Zongli, Li Niu, and Jack A Taylor. 2021. “The ENmix DNA methylation analysis pipeline for Illumina BeadChip and comparisons with seven other preprocessing pipelines.” Clinical Epigenetics 13 (1): 216. https://doi.org/10.1186/s13148-021-01207-1.