Skip to content

Getting started

Jump to bottom
Paul Ruiz edited this page Aug 21, 2026 · 1 revision

dnaEPICO provides a structured workflow for preprocessing Illumina DNA methylation arrays, fitting CpG-wise statistical models, and generating an interactive report. It supports EPICv2, EPIC, and 450K arrays.

The workflow can be run locally or on a High-Performance Computing (HPC) system using GNU Make.

Step 1: Install the required software

Before installing dnaEPICO, confirm that the following software is available:

  • R version 4.4 or later
  • GNU Make
  • Quarto command-line interface
  • Git, when installing from GitHub

Check the installations from a terminal:

R --version
make --version
quarto --version
git --version

Quarto is required for generating the web report. Installation instructions are available from the Quarto website.

R and Bioconductor packages should be managed with BiocManager, as recommended in the Bioconductor installation guide.

Step 2: Install dnaEPICO

Bioconductor version

Install the released version with:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("dnaEPICO")

GitHub development version

Install the current development version directly from the dnaEPICO GitHub repository:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install(
    "paulYRP/dnaEPICO",
    ask = FALSE,
    update = FALSE
)

The two installation commands are alternatives; it is not necessary to run both.

Confirm that the package can be loaded:

library(dnaEPICO)

packageVersion("dnaEPICO")

dnaEPICO installs its required R packages automatically. Array annotation packages and cell-composition reference packages may need to be installed separately according to the selected array and reference panel.

For an EPICv2 analysis using the 20a1.hg38 annotation:

BiocManager::install(c(
    "IlluminaHumanMethylationEPICv2manifest",
    "IlluminaHumanMethylationEPICv2anno.20a1.hg38"
))

Step 3: Prepare the project structure

Create a project directory with the following initial structure:

my_project/
├── Makefile
└── data/
    └── preprocessingMinfiEwasWater/
        ├── idats/
        │   ├── 2034567890_R01C01_Grn.idat
        │   └── 2034567890_R01C01_Red.idat
        ├── pheno.csv
        └── 12864_2024_10027_MOESM8_ESM.csv

Only the input files and Makefile need to be prepared manually. dnaEPICO creates the model-specific data, rData, results, figures, logs, and reports directories during the analysis.

Phenotype file

The column names in pheno.csv must agree with the names configured in the Makefile.

A typical longitudinal phenotype file could contain:

Sample_Name Basename Sentrix_ID Sentrix_Position Sex Timepoint Participant TreatmentGroup Age
Sample01 2034567890_R01C01 2034567890 R01C01 F 1 P001 Control 35
Sample02 2034567890_R02C01 2034567890 R02C01 F 2 P001 Control 35
Sample03 2034567890_R03C01 2034567890 R03C01 M 1 P002 Treatment 42

The main columns are:

  • Sample_Name: unique identifier for each array sample.
  • Basename: IDAT basename without _Red.idat or _Grn.idat.
  • Sentrix_ID: Illumina chip or slide identifier.
  • Sentrix_Position: position of the sample on the chip.
  • Sex: reported sample sex used for sex prediction and mismatch checking.
  • Timepoint: visit or collection time.
  • Participant: participant identifier used by longitudinal models.
  • TreatmentGroup: example phenotype to be tested.
  • Age: example numeric covariate.

Additional phenotype, technical, cell-composition, and adjustment variables can be included.

Each sample must have one red and one green IDAT file:

2034567890_R01C01_Red.idat
2034567890_R01C01_Grn.idat

Probe-exclusion file

The EPICv2 off-target probe reference used in the example Makefile can be downloaded from the associated Figshare dataset.

dnaEPICO can also use multiple probe-exclusion files or the expanded EPICv2 manifest options provided in the Makefile.

Step 4: Create and configure the Makefile

From R, move to the project directory and export the current Makefile template:

library(dnaEPICO)

extractMake(
    destDir = "/path/to/my_project",
    overwrite = FALSE
)

If a Makefile already exists and should be replaced:

extractMake(
    destDir = "/path/to/my_project",
    overwrite = TRUE
)

Edit the exported Makefile before running the pipeline.

Basic project settings

MODEL ?= model1

PHENO_FILE = $(DATA_DIR)/$(STEP1)/pheno.csv

SAMPLE_ID = Sample_Name
SEX_COLUMN = Sex
TIME_VAR = Timepoint
PERSON_VAR = Participant

MODEL is a label used to separate the outputs from different analyses.

Array settings

For EPICv2:

ARRAY_TYPE = IlluminaHumanMethylationEPICv2
ANNOTATION_VERSION = 20a1.hg38
ANNOTATION_PACKAGE = IlluminaHumanMethylationEPICv2anno.20a1.hg38

These settings must be changed when analysing EPIC or 450K arrays.

Number of samples

Use all samples for the full analysis:

N_SAMPLES = NA

A numeric value, such as 30, limits preprocessing to the first 30 phenotype rows and is intended only for testing.

Sex-mismatch handling

REMOVE_SEX_MISMATCH = FALSE
  • FALSE reports confirmed mismatches without removing the samples.
  • TRUE removes samples when the reported and methylation-predicted sex are both known and disagree.

Samples with missing or unknown sex cannot be classified as mismatches and are retained.

Methylation scale

METHYLATION_SCALE = beta

Valid Makefile values are:

  • beta: beta values
  • m: M-values
  • cn: copy-number values

Use lowercase values in the Makefile.

Timepoints

TIMEPOINTS = 1,2
COMBINE_TIMEPOINTS = 1,2

With these settings:

  • Separate phenotype-methylation objects are created for timepoints 1 and 2.
  • A combined longitudinal object is created for timepoints 1 and 2.
  • GLM uses the first timepoint by default.
  • LME uses the combined longitudinal object.

The GLM timepoint can be selected explicitly:

GLM_TIMEPOINT = 1

Model variables

PHENOTYPES_GLM = TreatmentGroup
PHENOTYPES_LME = TreatmentGroup

COVARIATES = Sex,Age
FACTOR_VARS = Sex,Timepoint,TreatmentGroup,Participant
SCALE_VARS = Age
  • PHENOTYPES_GLM lists the variables tested by GLM.
  • PHENOTYPES_LME lists the variables tested by LME.
  • COVARIATES lists adjustment variables.
  • FACTOR_VARS lists categorical variables.
  • SCALE_VARS lists numeric model variables to centre and divide by their standard deviations.
  • Use NULL when no variables should be scaled.

GLM settings

GLM_LIBS = glm2
INTERACTION_GLM = NULL

To fit a GLM interaction:

INTERACTION_GLM = Timepoint

Longitudinal lmerTest/lme4 settings

LME_LIBS = lme4,lmerTest
PERSON_VAR = Participant
INTERACTION_LME = Timepoint

LME_CORRELATION_STRUCTURE = none
LME_CORRELATION_VAR = NULL

With PHENOTYPES_LME = TreatmentGroup, the model contains:

TreatmentGroup × Timepoint + covariates + (1 | Participant)

The participant-level term (1 | Participant) is a random intercept.

Omnibus interaction test

To calculate one joint p-value for the complete phenotype-by-timepoint interaction:

LME_OMNIBUS_TEST = TRUE
LME_OMNIBUS_DDF = Satterthwaite

For a three-level phenotype, the omnibus test evaluates the interaction coefficients together and provides one p-value per CpG.

Available denominator degrees-of-freedom methods are:

LME_OMNIBUS_DDF = Satterthwaite

or:

LME_OMNIBUS_DDF = Kenward-Roger

Kenward–Roger requires the pbkrtest package:

install.packages("pbkrtest")

Omnibus testing is currently available for the lmerTest/lme4 engine, not the nlme engine.

Parallel processing and resumable summaries

N_CORES = 32
RESUME_FROM_SUMMARY = TRUE
CHUNK_SIZE = 10000
  • N_CORES is the maximum number of worker processes.
  • dnaEPICO may use fewer workers when the available CPUs, memory, or workload do not support the requested number.
  • RESUME_FROM_SUMMARY = TRUE reuses a compatible completed phenotype summary.
  • An interrupted phenotype without a complete summary restarts from its first CpG.
  • Full lmer, nlme, and glm2 model objects are not retained for every CpG.

Step 5: Run the pipeline

Run commands from the directory containing the Makefile.

Check configured outputs

make status MODEL=model1

Preprocessing, quality control, SVA, and phenotype preparation

make f3 MODEL=model1

This route generates:

  • RGSet and other preprocessing objects
  • beta, M-value, and copy-number matrices
  • quality-control plots and metrics
  • sex-prediction and mismatch information
  • probe-filtering results
  • cell-composition estimates
  • surrogate variable or PC columns
  • separate and combined phenotype-methylation objects
  • processing logs
  • an updated web report

The SVA step updates the single model-specific phenoLC.csv with the calculated PC columns.

Include the GLM analysis

make f4 MODEL=model1

This runs the first three stages, fits the CpG-wise GLM analysis, and refreshes the report.

The principal workbook is:

data/model1/methylationGLM/annotatedGLM.xlsx

The workbook contains the sheets:

  1. annotatedGLM
  2. metadata
  3. dictionary

Its result columns follow this general structure:

IlmnID
<phenotype coefficient p-values>
Name
chr
pos
UCSC_RefGene_Group
UCSC_RefGene_Name
Relation_to_Island
<phenotype>_Model.Message

The exact coefficient columns depend on whether the phenotype is numeric, binary, multilevel, or used in an interaction.

Include the longitudinal LME analysis

make f3lme MODEL=model1

This runs preprocessing and phenotype preparation, fits the longitudinal mixed-effects analysis, and refreshes the report.

The principal workbook is:

data/model1/methylationLME/annotatedLME.xlsx

The workbook contains:

  1. annotatedLME
  2. metadata
  3. dictionary

When omnibus testing is enabled, the result columns follow this general structure:

IlmnID
<interaction coefficient p-values>
<interaction omnibus p-value>
Name
chr
pos
UCSC_RefGene_Group
UCSC_RefGene_Name
Relation_to_Island
<remaining omnibus statistics>
<phenotype>_Model.Message

The omnibus statistics include:

  • F-statistic
  • numerator degrees of freedom
  • denominator degrees of freedom
  • raw p-value
  • adjusted p-value
  • denominator degrees-of-freedom method

Run the complete workflow

make all MODEL=model1

This runs preprocessing, phenotype preparation, GLM, LME, and report generation.

Run several configured models

If the Makefile contains:

MODELS = modelA modelB modelC

the available multi-model commands include:

make f3_models
make f4_models
make f3lme_models
make models

Each model must have an appropriate PHENO_FILE entry in the per-model configuration block.

Step 6: Inspect the outputs

The main generated directories are:

my_project/
├── data/model1/
├── rData/model1/
├── results/model1/
├── figures/model1/
├── logs/model1/
└── reports/model1/

Model workbooks

data/model1/methylationGLM/annotatedGLM.xlsx
data/model1/methylationLME/annotatedLME.xlsx

Compact phenotype summaries

rData/model1/methylationGLM/models/<Phenotype>SummaryGLM.rds
rData/model1/methylationLME/models/<Phenotype>SummaryLME.rds

These files contain compact result summaries used for reproducible output generation and phenotype-level restart. They do not contain one complete native model object per CpG.

Text summaries and significant-CpG files

results/model1/summary/
results/model1/cpgs/

Interactive report

reports/model1/docs/index.html

The report presents the GLM and LME results through paged viewers without embedding the entire result table in the HTML page.

The supporting TSV files are stored under:

reports/model1/assets/results/glm_results/
reports/model1/assets/results/lme_results/

Keep the complete reports/model1 directory together. Copying only index.html will separate the page from the result files that it needs.

Step 7: Monitor or clean a run

Check which expected outputs are present:

make status MODEL=model1

Remove generated outputs for one model:

make clean MODEL=model1

Remove outputs for every model listed in MODELS:

make clean_models

The clean commands retain the shared input directory containing the original IDAT and phenotype files.

Cluster/HPC option

dnaEPICO can be submitted as a PBS job. The following is an example and must be adapted to the local HPC environment:

#!/bin/bash -l
#PBS -N dnaEPICO_model1
#PBS -q batch_name
#PBS -l select=1:ncpus=32:mem=512gb
#PBS -l walltime=48:00:00
#PBS -j oe

module purge
module load GCC/13.2.0
module load R/4.4.1

cd "$PBS_O_WORKDIR"

Rscript -e "library(dnaEPICO); print(packageVersion('dnaEPICO'))"
quarto --version

make all MODEL=model1 N_CORES=32

The resource values are examples. Required memory and runtime depend on:

  • Number of samples
  • Number of CpGs
  • Number and type of covariates
  • GLM, lmerTest/lme4, or nlme engine
  • Use of interactions
  • Use of omnibus tests
  • Number of requested workers

N_CORES is treated as a maximum. dnaEPICO records the selected worker count and memory-related limits in the model log.

Submit the job with:

qsub pipeline.pbs

Monitor it with:

qstat -u "$USER"

Clone this wiki locally