Fit CpG-wise GLMs for one-timepoint methylation analyses

Source: R/methylationGLM.R

methylationGLM() prepares phenotype-plus-methylation data, fits one Gaussian GLM per CpG and phenotype, summarizes and annotates the results, and creates optional coefficient tables and diagnostic plots. It writes outputs only when saveOutputs = TRUE. Numeric CpG columns are passed to glm2::glm2() without a separate methylation-domain filter. Native model messages, warnings, and errors are recorded in one phenotype-specific Model.Message field. Annotated outputs contain CpGs with at least one returned coefficient or omnibus p-value; aggregate availability and condition counts are recorded in workbook metadata.

methylationGLM(
  inputPheno = "rData/preprocessingPheno/mergeData/phenoBT1.RData",
  outputLogs = "logs",
  outputRData = "rData/methylationGLM/models",
  outputPlots = "figures/methylationGLM",
  phenotypes = c("DASS_Depression", "DASS_Anxiety", "DASS_Stress", "PCL5_TotalScore",
    "MHCSF_TotalScore", "BRS_TotalScore"),
  covariates = paste0("Sex,Age,Ethnicity,TraumaDefinition,Leukocytes,",
    "Epithelial.cells"),
  factorVars = "Sex,Ethnicity,TraumaDefinition",
  scaleVars = NULL,
  cpgPrefix = "cg",
  cpgLimit = NA,
  methylationScale = "beta",
  nCores = 32,
  plotWidth = 2000,
  plotHeight = 1000,
  plotDPI = 150,
  interactionTerm = NULL,
  omnibusTest = FALSE,
  vennDPhenotypes = NULL,
  vennDLabels = NULL,
  vennDOmnibusPhenotypes = NULL,
  vennDOmnibusLabels = NULL,
  libPath = NULL,
  glmLibs = "glm2",
  prsMap = NULL,
  summaryPval = NA,
  summaryResidualSD = TRUE,
  saveSignificantCpGs = FALSE,
  significantCpGDir = "preliminaryResults/cpgs/methylationGLM",
  significantCpGPval = 0.05,
  saveTxtSummaries = TRUE,
  chunkSize = NULL,
  summaryTxtDir = "preliminaryResults/summary/methylationGLM",
  fdrThreshold = 0.05,
  padjmethod = "fdr",
  annotationPackage = "IlluminaHumanMethylationEPICv2anno.20a1.hg38",
  annotationCols = c("Name", "chr", "pos", "UCSC_RefGene_Group", "UCSC_RefGene_Name",
    "Relation_to_Island", "GencodeV41_Group"),
  gencodeHub = FALSE,
  annotatedGLMOut = "data/methylationGLM",
  reportAssetsDir = NULL,
  display = FALSE,
  verbose = FALSE,
  logs = FALSE,
  saveOutputs = FALSE,
  resumeFromSummary = TRUE
)

Arguments

inputPheno

Character. Path to the merged phenotype-plus-methylation .RData or .rds object created by preprocessingPheno(). The default points to the timepoint-1 object produced by the package workflow.

outputLogs

Character. Directory used for optional log files.

outputRData

Character. Directory used for compact, resumable phenotype summaries.

outputPlots

Character. Directory used for optional TIFF plots.

phenotypes

Character vector or comma-separated phenotype variables to model.

covariates

Character. Comma-separated covariate variables included in each GLM.

factorVars

Character. Comma-separated variables to convert to factors before modeling.

scaleVars

Character vector, comma-separated variable names, or NULL. Numeric fixed-effect variables to centre and divide by their sample standard deviations before model fitting.

cpgPrefix

Character. Prefix used to identify methylation columns in the merged phenotype-plus-methylation input object. The default is 'cg'.

cpgLimit

Integer or NA. Maximum number of CpGs to analyse. Use NA to keep all CpGs matching cpgPrefix.

methylationScale

Character. Methylation metric represented by the CpG columns. One of 'Beta', 'M', or 'CN', in any combination of upper- and lower-case letters. The default is 'beta'.

nCores

Integer. Maximum number of worker processes to use while fitting models. Automatic fitting remains serial below the glm2 crossover and caps workers by the CpG workload, available CPUs, and detected memory.

plotWidth

Integer. TIFF width in pixels when plots are written to disk.

plotHeight

Integer. TIFF height in pixels when plots are written to disk.

plotDPI

Integer. TIFF resolution in DPI when plots are written to disk.

interactionTerm

Character or NULL. Optional interaction term. When supplied and present in the input data, the phenotype is modeled together with its interaction against this variable.

omnibusTest

Logical. If TRUE, use car::linearHypothesis() to test the complete phenotype-by-interaction term, or the phenotype main effect when interactionTerm = NULL, once per CpG. One-degree-of-freedom terms are tested and therefore reproduce the corresponding coefficient p-value.

vennDPhenotypes

Character vector, comma-separated phenotype names, or NULL. Selected phenotypes are expanded to all coefficient p-value columns for model-level Venn diagrams and workbook tables.

vennDLabels

Character vector, comma-separated display labels, or NULL. Labels follow the resolved coefficient order and preserve case.

vennDOmnibusPhenotypes

Character vector, comma-separated phenotype names, or NULL. These use only omnibus p-value columns and require omnibusTest = TRUE.

vennDOmnibusLabels

Character vector, comma-separated display labels, or NULL, supplied in the resolved omnibus phenotype order.

libPath

Character vector or NULL. Optional library paths forwarded to worker processes. By default, the current .libPaths() are used.

glmLibs

Character. Comma-separated package names to validate on worker processes. The default is 'glm2'.

prsMap

Character or NULL. Optional phenotype-to-PRS mapping in the form 'Phenotype1:PRS_1,Phenotype2:PRS_2'.

summaryPval

Numeric or NA. Optional p-value threshold applied to the returned CpG summary tables. Use NA to keep all summary rows.

summaryResidualSD

Logical. If TRUE, append residual standard deviations to the CpG summary tables and residual diagnostic plots.

saveSignificantCpGs

Logical. If TRUE, collect significant CpG coefficient tables in the returned object and optionally write them to disk when saveOutputs = TRUE.

significantCpGDir

Character. Directory used for optional significant CpG coefficient tables.

significantCpGPval

Numeric. P-value threshold used to collect or write significant CpG coefficient tables. The threshold is applied to omnibus p-values when omnibusTest = TRUE, and to target coefficient p-values otherwise.

saveTxtSummaries

Logical. If TRUE and saveOutputs = TRUE, write tab-delimited summary tables to summaryTxtDir.

chunkSize

Integer or NULL. Number of CpGs processed per summary extraction chunk. NULL chooses a value automatically.

summaryTxtDir

Character. Directory used for optional tab-delimited GLM summary tables.

fdrThreshold

Numeric. False-discovery-rate threshold used to highlight CpGs in the residual-significance diagnostic plots.

padjmethod

Character. P-value adjustment method passed to stats::p.adjust(). The default is 'fdr'.

annotationPackage

Character. Annotation package or object name passed to minfi::getAnnotation(), for example 'IlluminaHumanMethylationEPICv2anno.20a1.hg38'.

annotationCols

Character vector or comma-separated annotation columns to append to the combined GLM summary table. Available columns depend on the selected annotation package.

gencodeHub

Logical. If TRUE, append release-aware GENCODE gene-body and nearest-TSS annotations obtained through AnnotationHub. The selected array annotation must use GRCh38 coordinates.

annotatedGLMOut

Character. Directory used for the optional annotated GLM summary XLSX workbook.

reportAssetsDir

Character or NULL. Report results directory used for the compressed TSV table and its compact metadata sidecars. NULL writes only the model outputs and annotated workbook.

display

Logical. If TRUE, draw exploratory and diagnostic plots on the active graphics device.

verbose

Logical. If TRUE, emit progress messages with message(). The default is FALSE.

logs

Logical. If TRUE, write the same progress messages to file.path(outputLogs, 'log_methylationGLM.txt').

saveOutputs

Logical. If TRUE, write compact phenotype summaries, text summaries, significant-CpG tables, annotated results, and TIFF plots. The default is FALSE.

resumeFromSummary

Logical. If TRUE and saveOutputs = TRUE, reuse a complete phenotype summary when its input file and model configuration match the current analysis. If processing stops before a phenotype summary is complete, that phenotype is fitted again from its first CpG.

Value

A list with class 'dnaEPICO_methylationGLM'.

preparedData

Object returned by prepareMethylationGLMData() containing the merged phenotype-plus-methylation analysis table and modeling metadata.

distributionPlots

Object returned by plotMethylationGLMDistributions() describing any exploratory plots that were generated or written.

designPlots

Missingness and numeric-correlation plots for the variables used in the model.

modelFits

Object returned by fitMethylationGLMModels() containing compact per-phenotype coefficient, omnibus, and condition results.

modelSummaries

Object returned by summarizeMethylationGLMModels() containing the combined CpG summary tables used for reporting and annotation.

significantCpGs

Object returned by collectSignificantCpGsMethylationGLM() containing optional phenotype-specific significant-CpG tables.

diagnosticPlots

Object returned by plotMethylationGLMDiagnostics() describing the diagnostic plot objects and any written TIFF files.

annotation

Object returned by annotateMethylationGLMSummaries() containing the annotated combined summary table.

manhattanPlots

Versioned circular and rectangular Manhattan plots for every raw p-value column in the annotated results.

vennDPlots

Requested coefficient and omnibus model-level Venn plots, worksheet tables, and label mappings.

savedFiles

Object returned by writeMethylationGLMOutputs() when saveOutputs = TRUE, otherwise NULL.

runSettings

High-level run metadata including the generic analysis label, methylation scale, display label, selected merged-object prefix, and internal response-column name.

See dnaEPICO_methylationGLM for a class-level overview.

Examples

if (requireNamespace(
    "IlluminaHumanMethylation450kanno.ilmn12.hg19",
    quietly = TRUE
)) {
    tmp <- tempdir()
    toy_path <- file.path(tmp, "phenoBT1.RData")
    phenoBT1 <- data.frame(
        Sample_Name = c("S1", "S2", "S3", "S4"),
        status = factor(c("Case", "Case", "Control", "Control")),
        sex = factor(c("F", "M", "F", "M")),
        cg00000029 = c(0.20, 0.25, 0.22, 0.27),
        cg00000108 = c(0.60, 0.55, 0.52, 0.58),
        check.names = FALSE
    )
    save(phenoBT1, file = toy_path)

    result <- methylationGLM(
        inputPheno = toy_path,
        phenotypes = "status",
        covariates = "sex",
        factorVars = "status,sex",
        cpgLimit = 2,
        nCores = 1,
        summaryPval = 1,
        annotationPackage = "IlluminaHumanMethylation450kanno.ilmn12.hg19",
        annotationCols = "Name,chr,pos",
        display = FALSE,
        verbose = FALSE,
        logs = FALSE,
        saveOutputs = FALSE
    )

    class(result)
}
#> [1] "dnaEPICO_methylationGLM"