Ultra-fast analysis of sparse DNA methylomes via MRMP (Most Recurrent Methylation Pattern) encoding — cell-type annotation, deconvolution, and genome-wide CpG upscaling.
GitHub· models· For coding agents· hg38 · mm10 · WGBS · single-cell · spatial
Every model is one self-contained bundle carrying its own MRMP feature definition (and
labels), so a query .cg runs directly — no separate annotation
files.
# conda — installs the `methscope` binary
conda install -c zhou-lab -c conda-forge methscope
# optional, linux-64 only: adds a `methscope-cuda` binary with the CUDA
# backend for `upscale-train`. Everything else — upscale, classify,
# deconv — is pure C and needs no GPU, so most users can skip this.
conda install -c zhou-lab -c conda-forge methscope-cuda
# or build from source (bundles YAME)
git clone --recurse-submodules https://github.com/zhou-lab/methscope-cli
cd methscope-cli && make # add CUDA=1 for the training backend
# install yame to fetch prebuilt model and test data
conda install -c zhou-lab -c conda-forge yameData and models are fetched with yame fetch; each workflow below
opens by fetching just the file it needs into the working directory
(yame fetch -c, digest-checked). It never prompts, so it is safe
inside a container build.
Atomic downloads. Downloads land on a .part
sibling and are renamed only once the bytes verify, so an interrupted
yame fetch can never leave a short file that later reads as a valid
model. Fetching now lives in yame, so keep it on
PATH alongside the conda env that provided
libxgboost.
An MRMP is one recurrent methylation pattern across a reference panel; the mask that
labels every genomic CpG with its pattern is what every model is built on.
mrmp-build constructs it natively from a discretized reference
.cg — no YAME + awk + sort pipeline, ~576 MB peak, ~50 s
for a whole-genome reference. The fetchable example below is chr20 only, so it finishes in
about a second.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
# a 40-cell-type chr20 reference (39.5 MB)
yame fetch -c hg38/data/human_hg38_40_celltypes_chr20.cg
methscope mrmp-build human_hg38_40_celltypes_chr20.cg chr20_40celltypes.mrmproot 40 classes | 491 patterns | 7,218 CpGs
split threshold: 20,000 CpGs -> 1 group(s)
! unsplittable (40 classes): GSM5652176_Adipocytes-Z000000T7 ... (34 more)
1 node(s), 29,904 bytes -> chr20_40celltypes.mrmp40 samples is the ceiling: a pattern packs as a base-3 uint64 and 3^40 < 2^64 < 3^41. The whole-genome models use 35 Zhou pseudo-bulks.
What it is: dimensions, the binstring parameters, the PNA share.
methscope inspect chr20_40celltypes.mrmpMRMP chr20_40celltypes.mrmp
format MRMPIDX1 v1, one set
name root
reference human_hg38_40_celltypes_chr20.cg
classes 40
patterns 491
covered 7,218 of 773,477 CpGs (0.93%)
PNA 766,259 CpGs (99.07%) match no pattern
block 29,904 bytes
membership 18,266 bytes, RLE + BGZF (169x vs dense)
resolution --mincov 1 --beta-threshold 0.500
checksum 8cf2026fa6de8d46 (over the per-CpG key stream)--patterns adds the ranked head; --top sets
how many (default 20).
methscope inspect chr20_40celltypes.mrmp --patterns --top 5#pattern label count (top 5 of 491)
1111111101111111111111111111111111111111 P1 698
1111101111111111111111111111111111111111 P2 527
1111111111011111111111111111111111111111 P3 478
1111111110111111111111111111111111111111 P4 429
1111111111111101111111111111111111111111 P5 337For a multi-set routing tree, --tree shows each MRMP
node and the classes it decides. It works on the tree bundle written by
mrmp-build.
methscope inspect chr20_40celltypes.mrmp --treeTREE chr20_40celltypes.mrmp
1 node, 40 classes, 491 patterns over 7,218 CpGs
root 40 classes 491 patterns 7,218 CpGs
|-- GSM5652176_Adipocytes-Z000000T7
|-- GSM5652181_Saphenous-Vein-Endothel-Z000000RM
...The runtime .cm mask is inspection / interop only, not a pipeline
step: upscale-train materializes it into the bundle from the
.mrmp itself.
methscope mrmp-export --top 1000 chr20_40celltypes.mrmp chr20_40celltypes.cmRanking is deterministic (count desc, key asc), so a build is byte-reproducible from the
reference alone. The artifact is the build pipeline's currency:
upscale-featurize, upscale-set-units, and
upscale-train all read it, so the msur's group map and the mask a
model ships cannot drift apart.
hg38 · 35 samples · 29,401,795 CpGs · 2,345,190 patterns · 155 MB header 128 B reference, sample names, binstring params, checksum patterns 16 B × 2,345,190 base-3 key + count, ranked 11111111111111111111111111111111111 10,736,116 P1 11111111111011111111111111111111111 2,859,193 P2 00000000000000000000000000000000000 2,214,235 P3 membership 4 B × 29,401,795 rank per CpG CpG 0 → 0 P1 CpG 7 → 4294967295 PNA
Classify a methylome against a bundled classifier. Each step fetches the file it needs with
yame fetch -c (digest-checked), so the block runs as-is on a clean
machine.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
# 4 Loyfer cells vs the shipped 33-type classifier
yame fetch -c hg38/data/human_hg38_celltypes.cg
yame fetch -c hg38/models/hg38_celltype.clfx
methscope classify hg38_celltype.clfx human_hg38_celltypes.cgcell prediction_label confidence certainty
Oligodendrocyte Oligodendrocyte 0.999738 0.999186
Pancreas-Beta Pancreatic.Islet.Cell 0.999662 0.998899
Blood-NK NK 0.978498 0.894815
Blood-Monocytes Monocyte 0.999126 0.993053Loyfer truth → Zhou2025 label: 4/4 concordant. The same call runs a sex / age / any bundled model. Linear models store per-feature training means, so missing/sparse betas are imputed at classify time — accurate even at <1% genome coverage.
yame fetch -c hg38/models/hg38_sex.clfx
methscope inspect hg38_sex.clfxcontainer MSBNDL1 (MethScope BuNDLe v1) - 3 sections, 30,931 bytes
# section offset size content
--- --------------- ----------- ----------- --------------------------------
0 mrmp 0 30,653 YAME .cm mask, 3 states
1 kind 30,761 8 "logistic" framework mark
2 model 30,769 154 methscope-linear text spec
[2] model logistic linear classifier - run via `classify`
method logistic
labels Female, Male
bias -7.50339
scale 1
features 2
Xa_lo weight=-13.9125 mean=0.24253
Xa_hi weight=13.6301 mean=0.752643hg38_celltype.clfx · 24.5 MB · 62 human cell types mrmp fmt2 MRMP definition — the feature set the model was fit on kind framework mark: xgboost | logistic | threshold model the booster / linear spec, labels embedded labels 62 cell types (linear models also store per-feature means)
.ubjx): one self-describing file
with the inner model, its MRMP feature definition, and the labels, so a query
.cg runs directly. Three frameworks (xgboost / logistic /
threshold) are supported and the kind section records which, so
classify needs no flag to know how to run it.The 40-cell-type chr20 reference doubles as training data: label each cell by lineage, straight from the sample names.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
yame fetch -c hg38/data/human_hg38_40_celltypes_chr20.cg
methscope mrmp-build --flat human_hg38_40_celltypes_chr20.cg chr20.mrmp
methscope mrmp-export --top 1000 chr20.mrmp chr20.cm
cut -f1 human_hg38_40_celltypes_chr20.cg.idx | sed 's/^GSM[0-9]*_//; s/-Z[0-9A-Z]*$//' | sed -E 's/.*Endothel.*|.*Endothelium.*/Endothelial/; s/.*Epithelial.*|.*Keratinocytes.*|.*Podocytes.*/Epithelial/; s/^Blood-.*|.*Macrophages.*/Immune/; s/.*Smooth-Muscle.*|Skeletal-Muscle|.*Cardiomyocyte.*/Muscle/; s/^Neuron$|.*Oligodendrocytes.*/Neural/; s/.*Fibroblasts.*|Adipocytes|Lung-Pleura/Mesenchymal/; s/^(Liver-Hepatocytes|Pancreas-Acinar|Pancreas-Beta)$/Secretory/' > labels.txtFeaturize once, then train from the artifact. The labels ride inside the
.msfm, so a later train cannot be handed a mismatched label
file, and the chain travels with it — classify-train --data
therefore takes no positional MRMP.
methscope classify-featurize -l labels.txt -o chr20.msfm \
human_hg38_40_celltypes_chr20.cg chr20.mrmp
methscope classify-train --data chr20.msfm -o chr20_lineage.clfx
methscope classify chr20_lineage.clfx human_hg38_40_celltypes_chr20.cgShape only — this particular example does not train. The default
xgboost framework needs --data; the flat
single-booster form these commands used to take was removed. What is still missing here is
the labels: a node trains over the classes its MRMP declares, which for
chr20.mrmp are the 40 reference samples, so a 7-lineage
labels.txt leaves it with no training rows for a
node's classes. Give a label set that matches the MRMP's own classes, or build the
MRMP from a store already labelled the way you intend to train.
--framework logistic | threshold are unfitted and still take the
<query.cg> <ref.cm> positional form.
The same wrapping by hand: unpack a bundle into its parts, then re-wrap
(unbundle writes chr20_lineage.ubj +
chr20_lineage.mrmp; -k re-stamps the kind).
Any .clfx works — the shipped models round-trip too.
methscope unbundle chr20_lineage.clfx
methscope bundle -m chr20_lineage.mrmp -k tree -o rewrapped.clfx chr20_lineage.ubjSolve a mixture against a signature panel by non-negative least squares. The solver rebuilds its pattern set per query from the CpGs that query measured, so it is robust to extreme sparsity and works at cfDNA-scale coverage.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
# nine known mixtures vs the shipped 62-type reference (1.5 GB)
yame fetch -c hg38/data/human_hg38_immune_mixture.cg
yame fetch -c hg38/models/hg38_62celltypes.msdref
methscope deconv hg38_62celltypes.msdref human_hg38_immune_mixture.cgcell class fraction
mac70_mono30_2pow22 Macrophage 0.696336
mac70_mono30_2pow22 Mono 0.303664
cd4_70_cd8_30_2pow22 Tmem_CD4 0.498298
cd4_70_cd8_30_2pow22 Tmem_CD8 0.268361
…Only the classes at or above --min-frac are listed, largest
first. --report gives the same thing to read rather than to
parse, and --wide the full 62-column matrix for a join.
methscope deconv --report hg38_62celltypes.msdref human_hg38_immune_mixture.cgmac70_mono30_2pow22: Macrophage 69.6%; Mono 30.4%
cd4_70_cd8_30_2pow22: Tmem_CD4 49.8%; Tmem_CD8 26.8%; Tnaive_CD4 20.7%; Tnaive_CD8 2.6%
tnaive_cd4_70_cd8_30_2pow22: Tnaive_CD4 69.9%; Tnaive_CD8 30.1%Every record has a known truth, so this checks an answer rather than merely producing one. Each name gives both the composition and the sparsity: three 70/30 pairs (macrophage/monocyte, CD4/CD8 T cell, naive CD4/naive CD8) at 222, 218 and 216 binarized CpGs. The pooled T-cell mixtures split across memory and naive sub-types, so sum those to score them. The 216 rung sits at the depth where a single draw becomes seed-dependent, which is why error grows down the ladder. The 62-type reference is built on the Zhou single-cell atlas plus Loyfer hepatocytes and BLUEPRINT granulocytes.
Per-cell panel, so each cell must deconvolve to itself.
yame fetch -c hg38/data/human_hg38_celltypes.cg
yame fetch -c hg38/models/hg38_celltype.clfx
methscope deconv-build-ref -o self.msdref human_hg38_celltypes.cg
methscope deconv self.msdref human_hg38_celltypes.cgcell class fraction
Oligodendrocyte Oligodendrocyte 1.000000
Pancreas-Beta Pancreas-Beta 1.000000
Blood-NK Blood-NK 1.000000
Blood-Monocytes Blood-Monocytes 1.000000Impute dense CpG methylation from a sparse methylome with a bundled decoder — pure-C
inference, no CUDA or BLAS at run time. The default output keeps the predicted
fraction per CpG (format 4); --binary thresholds it at 0.5
into 0/1 calls (format 6), and --probs writes the same values as
a TSV.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
# what goes in: a sparse methylome
yame fetch -c hg38/data/human_hg38_test.cg
yame summary human_hg38_test.cgQFile Query MFile Mask N_univ N_query N_mask N_overlap Log2OddsRatio Beta Depth
human_hg38_test.cg 1 NA global 29000 23857 29000 23857 NA 0.823 NA# impute the whole genome from it — fractions, the default
yame fetch -c hg38/models/hg38_10k1.updecx
methscope upscale -o human_hg38_test_reconstructed.cg hg38_10k1.updecx human_hg38_test.cg
yame summary human_hg38_test_reconstructed.cgQFile Query MFile Mask N_univ N_query N_mask N_overlap Log2OddsRatio Beta Depth
human_hg38_test_reconstructed.cg 1 NA global 29401795 10000 0 10000 NA 0.617 NA# or 0/1 calls, for anything that wants a set rather than a level
methscope upscale --binary -o human_hg38_test_calls.cg hg38_10k1.updecx human_hg38_test.cg
yame summary human_hg38_test_calls.cgQFile Query MFile Mask N_univ N_query N_mask N_overlap Log2OddsRatio Beta Depth
human_hg38_test_calls.cg 1 NA global 10000 6239 10000 6239 NA 0.624 NAThis decoder writes one 10,000-CpG block. The two rows describe the same prediction from different sides: the fraction file spans the genome row space and reports 10,000 predicted CpGs at mean beta 0.617, while the binarized one is a set — 6,239 of the block's 10,000 CpGs called methylated, so its beta 0.624 is that fraction rather than a mean level. Keep the default unless a downstream tool demands format 6: thresholding is lossy, and it is what the accuracy check below needs.
Score against the truth, which is a set of 0/1 calls — so compare the
--binary output, not the fractions. Against the fraction file every
row differs and this prints 0/9794 (0.0%), which looks like a
broken model rather than a mismatched format.
yame fetch -c hg38/data/human_hg38_test.truth.cg
paste <(yame rowsub -I 1_10000 human_hg38_test_calls.cg | yame unpack -a - | cut -f4) <(yame rowsub -I 1_10000 human_hg38_test.truth.cg | yame unpack -a - | cut -f4) | awk '$2!=2{t++;if($1==$2)c++} END{printf "%d/%d (%.1f%%)\n",c,t,100*c/t}'9211/9794 (94.0%)See the whole block first. The view is coordinate-aware — hprint
resolves the hg38 coordinate track from the store (fetched below), window-averaging 10,000 CpGs
into 60 columns; -g shows deciles.
cat human_hg38_test.truth.cg human_hg38_test.cg human_hg38_test_reconstructed.cg > three.cg
yame index -s <(printf 'truth\ninput\nreconstructed\n') three.cg
yame fetch hg38/cpg_nocontig.cr
yame hprint -g -r chr1:921649-1151482 -w 60 three.cg chr1:921649-1151482 (10000 CpGs, win=167)
|.........|.........|.........|.........|.........|.........
921649 958945 995856 1033422 1067326 1112421
truth 104523059802945587792005953275337999952335758999995168975689
input ...0.....5.0...........09...............9..0....9...9....9..
reconstructed 1066340498029456878930069532854479999525367589999951689757890 unmethylated through 9 methylated, . no coverage, coloured blue to red. Each column averages 167 CpGs. The input row is mostly gaps, and where it does show a digit that digit is one or two observed CpGs standing in for the whole window — so it reads lower or higher than the truth because it is a far thinner sample of the same window, not because it disagrees with it. Truth and reconstruction, both dense, track each other across the block.
Zoom to single CpGs over the same block the score used — the same
rowsub -I as above. -c turns colour off,
which cut needs since it counts bytes, not columns.
cat human_hg38_test.truth.cg human_hg38_test.cg human_hg38_test_reconstructed.cg | yame rowsub -I 1_10000 - | yame hprint -c - | cut -c1621-1680truth 111111010111111111111110111100000010000000000000000000000000
input 222222222222222221222222222222222222222222220222222222222222
recon 11111101011111111111111011110000001000000000000000000000000060 CpGs from the middle of the block. Two were observed, and the reconstruction matches the truth at all 60 — including the isolated 0 at position 7 and the lone 1 sitting in the unmethylated run.
Training the decoder is a separate, three-step pipeline, and every step reads the same
.mrmp.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
# the whole pipeline from scratch, on the 40-cell-type chr20 reference
# 1. define the features
yame fetch -c hg38/data/human_hg38_40_celltypes_chr20.cg
methscope mrmp-build --flat --force human_hg38_40_celltypes_chr20.cg chr20.mrmproot 40 classes | 491 patterns | 7,218 CpGs
1 node(s), 29,904 bytes -> chr20.mrmp# 2. featurize the truth atlas -> sampling/truth msur
methscope upscale-featurize --reps 20 --sample 8000 human_hg38_40_celltypes_chr20.cg chr20.mrmp chr20.msur[methscope] upscale-featurize: inflating 40 truth cells (1 worker, about 0.2 GB resident)
[methscope] upscale-featurize: wrote chr20.msur and chr20.msur.tsv# 3. group CpGs into processing units
methscope upscale-set-units --unit-cpgs 4096 human_hg38_40_celltypes_chr20.cg chr20.msuiunits=109 pure=13 mixed=67 PNA=29 oversized=7
unit CpGs min=1380 mean=7096.1 max=261816inspect reads every build artifact, detected from its magic.
methscope inspect chr20.msuiformat MSUIDX1 v1
cpgs 773477 (657409 real + 116068 PNA)
pattern_length 40
units 109 (29 PNA)methscope inspect chr20.msurformat MSURAW2 v2
cells 40
replicates 20
embedded_truth yes (trainable)hg38 · 207 cells × 100 reps · 29,401,795 CpGs · 15.0 GB offset section size 0 header 72 B 72 groups u16 × 29.4M 58.8 MB CpG → 1..1000, 0 = PNA 58,803,662 truth u16 × 207 × 29.4M 12.2 GB beta×65534, 65535 = missing (always embedded) 12,231,146,792 records 124,000 B × 20,700 2.6 GB 1000 beta + 1000 count + 29,000 CpG ids
.cg is never copied.hg38 · 16,384 CpGs/unit · 700 units · 2,345,190 memberships · 174 MB header 26,841,955 real CpGs + 2,559,840 PNA; checksum ties it to .mrmp units 24 B × 700 {offset, first_member, n_member, n_cpg, flags} u0 10,736,116 CpGs pure (P1, all-1) u1 2,859,193 CpGs pure (P2) u2 2,214,235 CpGs pure (P3, all-0) cpg 4 B × 29,401,795 output slot → genomic CpG index membership 24 B × 2,345,190 {key, offset, count, unit}
hg38 · 1000 MRMP in · 700 units · 29,401,795 CpGs out · 2.8 GB header patterns, input_dim, n_units, activation, section offsets mean/scale f32 × input_dim train-cell standardization units 32 B × 700 {out offset, param offset, bytes, n_cpg, rank} cpg 4 B × 29,401,795 output slot → genomic CpG index params f32 per unit: z = act(A·x + a); logit = E·z + b pure rank 16 · mixed/PNA rank 32
.cm mask.Training runs on CPU by default and threads over units, which are independent — that is what makes the run resumable. The 109 chr20 units take ~17 s on 8 threads.
methscope upscale-train -i chr20.msur --units chr20.msui --mrmp chr20.mrmp -o chr20.updecx --work-dir work/ --threads 8P=491 features=beta+count input=982 trunk=none pure=16 mixed=factor/32
backend=cpu
units 109/109 (resumed=0) last_val_mae=0.193608
wrote bare UPDEC2 work//model.updec2 (95924368 bytes), resumed=0The validation MAE moves a little between runs — 0.171 and 0.194 on two passes over this same input — so treat it as a scale, not a checksum. The unit count, the pattern count and the model size are stable.
make CUDA=1 CUDA_HOME=/path/to/cuda CUDA_ARCH=sm_80 builds the GPU
backend, which is picked automatically when a device answers; checkpoints are interchangeable,
so you can start on CPU and finish on a GPU node. The chr20 demo is small; a whole-genome
decoder is memory- and GPU-heavy, so submit it to a GPU node:
sbatch -q gpu --gres=gpu:1 --mem=64G -c 8 -t 8:00:00 --wrap 'methscope upscale-train -i wg.msur --units wg.msui --mrmp wg.mrmp -o hg38_wg.updecx --work-dir work/'--split rows: <cell_index>TAB<train|val|test>,
one per source cell; without it, cells are cut 70/15/15 by --seed.
--device cpu forces the portable backend even on a GPU node.