DMSA models probes as the
multivariate response, so every engine wants the same two objects: a
samples x probes matrix and a phenotype
data.frame with rows in the same order. No mainstream
preprocessing pipeline ends there. minfi ends in a
GenomicRatioSet (probes in rows), sesame in a beta matrix
(probes in rows), meffil and ewastools in matrices (probes in rows),
ChAMP in myNorm plus a pd sheet, RnBeads in an
RnBSet, methylprep in a CSV export. The gap between those
shapes and the analysis is where silent mistakes live: a transposed
matrix, M-values treated as betas, a sample sheet in a different order
than the arrays.
dmsa_import() is the one door in. It takes whichever end
object a pipeline produced and returns an object
dmsa_frame() accepts directly:
(For a curated SummarizedExperiment where you want to
pick the assay yourself - a recountmethylation download,
say - the lower-level dmsa_se() does the split with no
further opinions.)
library(dmsa)
set.seed(1)
## the shape most pipelines produce: probes in rows, samples in columns
B <- matrix(runif(60), nrow = 6,
dimnames = list(sprintf("cg%08d", 1:6), paste0("S", 1:10)))
ph <- data.frame(Sample_Name = paste0("S", 1:10), age = rnorm(10))
imp <- dmsa_import(B, pheno = ph)
imp
#> dmsa_import: 10 samples x 6 probes
#> source: matrix
#> values: betas
#> alignment: matched on pheno column 'Sample_Name'
#> alpha map: 0 of 6 probes are Alpha-mapped
#> note: orientation: probe ids in rownames
#> next: dmsa_frame(<this object>, outcome = ..., ...)Three things happened, none silently. The probe axis was identified
from the dimension names and the matrix transposed; the value scale was
detected (betas here - M-values are converted with a note, because
dmsa_frame() converts betas to M-values itself and arriving
twice-converted would corrupt the frame); and the sheet was matched to
the sample names by identifier. Each of those is either unambiguous or
an error that says what to pass:
orientation = when neither axis has probe-looking names,
values = to override the scale, align_by = to
name the identifier column (minfi Basename paths and
Sentrix_ID/Sentrix_Position pairs are tried
automatically; align_by = ".position" declares same-order
alignment explicitly). A raw container - an RGChannelSet, a
list of sesame SigDFs - is refused with the preprocessing
recipe rather than half-handled.
None of these packages are dmsa dependencies; the chunks are not run here.
## minfi - any preprocessed object works via the SummarizedExperiment route
rg <- minfi::read.metharray.exp(targets = targets)
gr <- minfi::ratioConvert(minfi::mapToGenome(minfi::preprocessNoob(rg)))
imp <- dmsa_import(gr) # colData is the sheet
imp <- dmsa_import(minfi::preprocessRaw(rg)) # Meth/Unmeth handled too
## sesame - openSesame() returns the beta matrix directly
imp <- dmsa_import(sesame::openSesame(idat_dir), pheno = sheet)
## ewastools
mt <- ewastools::read_idats(sheet$Basename)
imp <- dmsa_import(ewastools::dont_normalize(mt), pheno = sheet)
## meffil
qc <- meffil::meffil.qc(samplesheet, verbose = FALSE)
norm <- meffil::meffil.normalize.quantiles(qc, number.pcs = 5)
imp <- dmsa_import(meffil::meffil.normalize.samples(norm), pheno = samplesheet)
## ChAMP
load <- ChAMP::champ.load(directory, arraytype = "EPIC")
imp <- dmsa_import(ChAMP::champ.norm(beta = load$beta), pheno = load$pd)
## RnBeads - extract, then import
imp <- dmsa_import(RnBeads::meth(rnb.set, row.names = TRUE),
pheno = RnBeads::pheno(rnb.set))
## methylprep - export betas to CSV, then import the file
imp <- dmsa_import("beta_values.csv", pheno = "sample_sheet.csv")If preprocessing choice changed DMSA’s conclusions, the adapter would be a convenience and a hazard. So we put one plate of real arrays through five pipelines and asked the same two pre-registered questions of each result.
Data. 94 EPIC v1 arrays (one plate; 48 couples’ families, 12 chips) from Project Alpha’s first batch - and, importantly, a different physical scan of these participants than the arrays behind the package’s bundled maps: the two sets share zero barcodes. Person- and sample-level covariates (sex, age, epithelial and fibroblast fractions) were joined from the canonical study table by participant id; array-level structure (chip) entered as a random intercept over this plate’s own 12 chips; the canonical arrays’ control-probe surrogate variables were omitted, since they describe other physical arrays. Across pipelines, everything was therefore identical except the preprocessing.
Questions. (A) The pre-registered positive control:
the immune, inflammation & HLA system against the immune-cell
fraction IC (with the registered deviation - the epithelial fraction
leaves the covariate set, because IC = 1 - Epi - Fib makes the full set
rank-deficient). Cell composition is a signal known to dominate blood
methylation; a method that misses it is broken. (B) A theory pairing
expected to stay quiet at this sample size: the HPA axis against a
symptom inventory (BSI). All runs:
system > module > gene > probe, B = 1999
permutations, seed 1.
Result. The control fired at the permutation floor in all five pipelines - p = 0.0005 is 1/(B+1): not one of 1,999 permutations reached the observed statistic - including on ewastools’ completely unnormalized betas:
| pipeline | preprocessing | usable probes | system p | aligned direction |
|---|---|---|---|---|
| minfi noob | noob background/dye | 34 | 0.0005 | + |
| minfi raw | none (Meth/Unmeth ratio) | 35 | 0.0005 | + |
| meffil | functional normalization | 34 | 0.0005 | + |
| sesame | openSesame default | 26 | 0.0005 | + |
| ewastools | none | 34 | 0.0005 | + |
The agreement runs the whole cascade down. All 8 immune modules were selected in every pipeline. At gene level, 18-20 genes were significant per pipeline and 17 in all five (Jaccard 0.81): the HLA class I/II battery (HLA-DMB, -DOB, -DPA1, -DPB1, -DQA2, -DRA, -E, -F, HLA-F-AS1) with ICAM1, IL10, IL12RB1, IL12RB2, IL18, NLRP3, TLR4 and TNF; gene-level aligned directions were identical across all five pipelines for 92% of the shared genes. And the null behaved: HPA x BSI came back non-significant in every pipeline (p = 0.41-0.81), with the cascade never descending - zero modules and zero genes selected, in any pipeline.
The one visible pipeline difference is coverage, not conclusions: sesame’s default masking left 26 usable probes where the others kept 34-35 - a deliberate conservatism, reported by the frame, changing nothing downstream.
Why this is expected, not lucky. DMSA’s families are declared in the map rather than harvested from whatever the analyst loaded, the probe directions are fixed by the direction catalogue before any data arrive, and the permutation null is computed under the same transform as the observed statistic. Preprocessing can move individual betas and drop probes; what it does not get to do is redefine the hypothesis, the alignment, or the reference distribution. The pipeline choice moves probe counts; the answer stays the answer.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [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: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] dmsa_0.99.20 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.39 R6_2.6.1 fastmap_1.2.0 xfun_0.62
#> [5] maketools_1.3.2 cachem_1.1.0 knitr_1.52 htmltools_0.5.9
#> [9] buildtools_1.0.0 lifecycle_1.0.5 cli_3.6.6 sass_0.4.10
#> [13] jquerylib_0.1.4 compiler_4.6.1 sys_3.4.3 tools_4.6.1
#> [17] evaluate_1.0.5 bslib_0.12.0 yaml_2.3.12 otel_0.2.0
#> [21] jsonlite_2.0.0 rlang_1.3.0