SEnsible Step-wise Analysis of DNA MEthylation

latest cli ↗ Cite ↗

Infinium analysis made simpler — the same betas as the R package, ~20× faster.

GitHub· R / Bioconductor· YAME· InfiniumAnnotation· cinderplot· EPIC · EPICv2 · HM450 · MSA

Install & fetch — then fully offline

# build from source (conda install is up top)
git clone --recurse-submodules \
  https://github.com/zwdzwd/sesame-cli
cd sesame-cli && make

# one-time: annotation into the shared store (yame downloads)
yame fetch InfiniumAnnotation/EPICv2
yame fetch genomes/hg38

…builds a local, digest-verified store ($YAME_DATA_HOME, default ~/.local/share/yame) — offline from here on:

EPICv2/  # probe annotation
├ EPICv2.ordering.tsv.gz
├ EPICv2.hg38.mask.cm
├ EPICv2.hg38.coord.tsv.gz
└ SHA256SUMS
genome/hg38/  # genome tracks
├ seqinfo.tsv.gz  gaps.tsv.gz
├ cytoband.tsv.gz  genes.bed.gz(+.tbi)
└ SHA256SUMS

Preprocess — IDAT → betas

QCDPB (openSesame default) over a cohort in one process → YAME .cg + qc.tsv.

# a folder of IDATs, searched recursively
sesame preprocess --out out/ in/

# or explicit prefixes (just betas)
sesame preprocess --output beta --out out/ \
  idat_prefix_1 idat_prefix_2

…writes one .cg per output over the cohort (each with a .cg.idx of sample names) + qc.tsv:

out/
├ beta.cg (+.idx)          # betas (fmt4)
├ intensity.cg (+.idx)     # M/U (fmt3)
├ pval.cg (+.idx)          # detection p
└ qc.tsv                   # QC metrics

qc.tsv holds 65 per-sample metrics. frac_dt — the fraction of probes above detection — is the most critical for assessing data quality (also per type, frac_dt_cg/_ch/_rs).

# head qc.tsv  (key columns)
sample    frac_dt  frac_dt_cg  frac_dt_rs
sampleA   0.978    0.979       1.00
sampleB   0.994    0.995       1.00

More knobs: --prep QCEPB swaps in linear dye bias (E = dyeBiasL); --collapse averages EPICv2/MSA replicate probes to their cg-prefix (betasCollapseToPfx); --detection pneg emits the negative-control-ECDF detection p; and GCT (bisulfite conversion) joins qc.tsv on every platform.

Bimodal distribution of GM12878 beta values — GM12878 betas — the classic bimodal shape: unmethylated (~0) and methylated (~1).

Read a .cg back — attach-probe

Betas and signals are stored positionally — one value per probe, in a fixed order, with no Probe IDs inside. The IDs live once in the platform's ordering table, which is what keeps the files small and byte-for-byte alignable: stack N samples into a matrix by column with no key-matching. attach-probe joins the two back into a labelled TSV.

$ sesame attach-probe out/beta.cg
Probe_ID          206909630040_R03C01
cg00000029_TC21   0.9124
cg00000769_BC11   0.0331
# the ordering is inferred from the row count -- each
# platform's is a distinct length. Name it if you'd rather:
$ sesame attach-probe --platform EPICv2 out/beta.cg
# --all for every sample column (matrix)

Peek at a raw IDAT — idat-dump

The bead-level record as Illumina wrote it, before any of it becomes a probe.

$ sesame idat-dump --tsv sample_Grn.idat.gz
address    mean  sd   nbeads
10600313   4821  118  12
10600322   287   34   9

Array QC plots — preprocess → cinderplot

The two diagnostics from sesame's supplemental QC, with no R in the loop: preprocess writes the signals, attach-probe labels them, tabl reshapes them, and cinderplot draws them.

Dye bias Q-Q plot, Infinium-I red versus green signal — **Dye bias Q-Q** — green quantiles below the diagonal: red reads hotter across the range. Re-run after `--prep CD` and the curve should land on the line. Ports `sesameQC_plotRedGrnQQ()`.

Total intensity versus beta, as a binned density — **Intensity vs beta** — two lobes at beta ~0.05 and ~0.9, funnelling to 0.5 as signal dies into background below ~2¹⁰. That funnel is where a detection cutoff belongs. Ports `sesameQC_plotIntensVsBetas()`.

Differential methylation — dml

Per-probe OLS with t/F tests + BH, matching R's lm to ~1e-9.

sesame dml --betas out/beta.cg \
  --index EPICv2.ordering.tsv.gz \
  --meta samples.tsv \
  --formula '~ group + age' > dml.tsv

De-identify — deidentify

Strip the genetic fingerprint: zero (or reversibly scramble) the SNP (rs) probe intensities in an IDAT — every other bead is byte-identical, so betas are unchanged except the SNPs.

# zero SNP beads (a prefix does both channels)
sesame deidentify sample_prefix

# or reversibly scramble, then restore
sesame deidentify --randomize --seed 42 s.idat
sesame deidentify -r --seed 42 s_noid.idat

Copy number — cnv

OLS vs a normal panel → log2 → genome bins → CBS (DNAcopy port). Two output files.

sesame preprocess --prep "" --raw-signal \
  --output total_intensity --out t/ tumor
sesame cnv --platform EPICv2 \
  t/total_intensity.cg segments.tsv bins.tsv
# --probes: per-probe; --exclude chrY

cnv">

# plot the two files with cinderplot
cinderplot 'out/bins.tsv
  + aes(chrom=chrom, x=start, xend=end, y=log2ratio, colour=log2ratio)
  + geom_point(size=.5)
  + geom_segment(data="out/segments.tsv", y=seg.mean, color="black")
  + scale_x_genome("hg38/seqinfo.tsv.gz")
  + scale_colour_gradient2(low="#762a83", high="#b2182b")
  + ideogram("hg38/cytoband.tsv.gz")'

K562 (EPICv2/hg38): chr9p deletion, chr13/22 amps.

Genotype SNPs → VCF — vcf

SNP (rs) probes → VCF (formatVCF); GT + fractions exact vs R.

sesame vcf tumor --platform EPICv2 \
  --snp EPICv2.hg38.snp.tsv.gz > geno.vcf
# --variants: only rs + channel-switch

Region & gene views — region → cinderplot

A locus's betas as long-form TSV → cinderplot over the gene models in zhou-lab/genomes.

DNMT3A methylation on a probe-anchored heatmap — **--gene DNMT3A** — gene models + probe-anchored heatmap.

HOXA cluster methylation on a probe-anchored heatmap — **chr:beg-end** — the HOXA cluster.

Speed — vs the R package

One EPICv2 sample, IDAT → betas, median of 3. R is timed twice: the openSesame() call alone, and the whole invocation — startup and library(sesame) included — which is what a pipeline pays, since it starts a fresh R per sample.

Runtime comparison: sesame C versus R openSesame — **2.5 s** vs **30.7 s** in-session (**12×**) and **46.9 s** per invocation (**19×**). One machine, one array — regenerate rather than quoting these.

Betas utilities — mliftover · impute

Post-process a beta.cg → beta.cg, positional to the platform ordering.

# lift betas across platforms (mLiftOver)
sesame mliftover --to EPIC --platform EPICv2 \
  beta.cg beta.EPIC.cg

# fill NA betas (mean / genomic neighbours)
sesame impute --method mean beta.cg out.cg
sesame impute --method neighbors \
  --platform EPICv2 beta.cg out.cg

mliftover prefix-joins to the target ordering (EPICv2↔EPIC↔HM450); impute fills NAs by the probe's cross-sample mean, or the nearest same-strand genomic neighbours (imputeBetas*).

Both are local fills — a neighbouring CpG or the column mean — and they only reach probes the array already carries. For imputation as an inference problem, use methscope: methscope upscale predicts genome-wide CpG methylation from a sparse methylome with a trained model, which is both more accurate and not limited to the manifest.

A platform is its ordering table — --platform takes one by path as readily as by name.

The ordering table — IDs live here, once

Probe_ID           M         U         col
cg00000029_TC21    NA        41791408  2
cg00000769_BC11    26756338  59799309  G
ctl_1060944_NEG    NA        10609447  2

Four columns and nothing else — not the Illumina/sesame manifest with tabs instead of commas. M is the methylated bead address (Illumina AddressB_ID, NA for Infinium-II), U the unmethylated one (AddressA_ID), and col the design: G / R Infinium-I in green or red, 2 Infinium-II. Sequences and coordinates live in separate files a custom array can do without. A legacy fifth mask column is still read; published orderings keep the mask in a companion .cm.

Tab- or comma-separated, plain or gzipped — the separator is sniffed from the header, so a table exported straight out of a CSV manifest needs no conversion. A header line is required.

Row order is load-bearing and permanent. A .cg stores one value per row with no Probe_IDs inside, so row i is probe i forever — sorting or inserting a row re-labels every .cg written against the old file. Freeze it when you publish, and version it if it changes.

Your own array — hand it the ordering

--platform takes a name or a path, so a custom or pre-release array needs no special flag. Every companion asset has its own override, inferred only when the platform is a name you published.

# named platform: ordering, mask, coords all inferred
sesame preprocess --platform EPICv2 --out out/ idats/

# custom array: the ordering IS the platform
sesame preprocess --platform MyArray.ordering.csv \
  --mask MyArray.mask.cm --prep QCDPB --out out/ idats/

# without --mask, prep C or "" work; Q/P/B need one
sesame preprocess --platform MyArray.ordering.csv \
  --prep C --out out/ idats/

Overrides: --mask --coords --snp --normals. From an Illumina manifest the ordering is a projection: IlmnID→Probe_ID, AddressB_ID→M, AddressA_ID→U, design + colour channel→col. Full spec in the README.

Making the mask — a .cm bitset

A .cm is a bitset with one bit per ordering row, in the same order — 1 means drop the probe. Write the column, pack it, name the track:

# one 0/1 per ordering row, in the ordering's order
awk -F'\t' 'NR>1{print (bad) ? 1 : 0}' \
  MyArray.ordering.csv > m.txt

yame pack -f b m.txt MyArray.mask.cm   # b = bitset
echo M_custom > names.txt
yame index -s names.txt MyArray.mask.cm  # name it

Given directly, --mask wants a single-track .cm — a published platform's mask holds many (EPICv2 has 25) and sesame picks the recommended subset by platform name, so handing one over without a name would be ambiguous. It says so rather than guessing.

Citation — the SeSAMe papers

sesame-cli is a C reimplementation of the SeSAMe method family. Please cite the work behind the parts you use:

Infinium annotation & masking — Zhou, Laird & Shen. Nucleic Acids Research 2017. 10.1093/nar/gkw967
SeSAMe / pOOBAH — Zhou, Triche, Laird & Shen. Nucleic Acids Research 2018. 10.1093/nar/gky691
mLiftOver — Chen & Zhou. Bioinformatics 2024. 10.1093/bioinformatics/btae423
KnowYourCG (KYCG) — Goldberg et al. Science Advances 2025. 10.1126/sciadv.adw3027
MSA array — Goldberg et al. Cell Genomics 2025. 10.1016/j.xgen.2025.100929
ELBAR / low-input — Lee et al. Nucleic Acids Research 2024. 10.1093/nar/gkae127
Cross-species analysis — Ding, Kaur, Horvath & Zhou. Briefings in Bioinformatics 2023. 10.1093/bib/bbac617
MM285 mouse array — Zhou et al. Cell Genomics 2022. 10.1016/j.xgen.2022.100144