Ultra-fast analysis of sparse DNA methylomes via MRMP (Most Recurrent Methylation Pattern) encoding — cell-type annotation, deconvolution, and genome-wide CpG upscaling.
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
# or build from source (bundles YAME)
git clone --recurse-submodules https://github.com/zhou-lab/methscope-cli
cd methscope-cli && make
# 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.mrmpsamples=40 CpGs=773477 distinct patterns=116450 (+PNA) candidates=116450
PNA CpGs=116068 (15.01%) checksum=9e7f6d292a2ba94640 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.mrmpformat MRMPIDX1 v1
reference human_hg38_40_celltypes_chr20.cg
samples 40
cpgs 773477
distinct_candidates 116450
selectable_patterns 116450
pna_cpgs 116068 15.0060%
binstring mincov=1 beta=0.500 max_ambig_frac=1.000 min_major_fold=10.000
content_checksum 9e7f6d292a2ba946--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 116450)
1111111111111111111111111111111111111111 P1 261816
0000000000000000000000000000000000000000 P2 58021
1111101111111111111111111111111111111111 P3 22839
1111111101111111111111111111111111111111 P4 11536
1111111111111011111111111111111111111111 P5 7027The 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 Zhou2025 62-type atlas
yame fetch -c hg38/data/human_hg38_celltypes.cg
yame fetch -c hg38/models/hg38_celltype.ubjx
methscope classify human_hg38_celltypes.cg hg38_celltype.ubjxcell prediction_label confidence
Oligodendrocyte ODC 0.915134 # oligodendrocyte
Pancreas-Beta Beta 0.930708 # islet beta
Blood-NK NK CD16 0.812749 # natural killer
Blood-Monocytes Mono 0.835970 # monocyteLoyfer 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.ubjx
methscope inspect hg38_sex.ubjxcontainer MSBNDL1 (MethScope BuNDLe v1) - 3 sections, 30,931 bytes
# section offset size content
0 mrmp 0 30,653 fmt2 MRMP definition, 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.ubjx · 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 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.txtFeatures + labels → one runnable bundle, in a single step (the positional takes the
runtime .cm, not the .mrmp artifact).
methscope classify-train -l labels.txt -o chr20_lineage.ubjx human_hg38_40_celltypes_chr20.cg chr20.cmtrained 7-class xgboost model on 40 cells x 1000 feature(s), 6 rounds
-> chr20_lineage.ubjx (booster+MRMP bundle)--framework xgboost | logistic | threshold.
methscope classify human_hg38_40_celltypes_chr20.cg chr20_lineage.ubjxcell prediction_label confidence
GSM5652176_Adipocytes-Z000000T7 Mesenchymal 0.229483
GSM5652181_Saphenous-Vein-Endothel-Z000000RM Endothelial 0.377817
GSM5652187_Kidney-Tubular-Endothel-Z000000PX Endothelial 0.36265640/40 here, but this is the training set — a sanity check, not a validation.
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).
methscope unbundle chr20_lineage.ubjx
methscope bundle -m chr20_lineage.mrmp -k xgboost -o rewrapped.ubjx chr20_lineage.ubjSolve a mixture against a signature panel by non-negative least squares. Uses all MRMP patterns and is robust to extreme sparsity, so it works on cfDNA-scale coverage.
# a scratch dir for the outputs
mkdir -p ~/tmp/methscope && cd ~/tmp/methscope
# a whole-body mixture vs the shipped 65-type reference
yame fetch -c hg38/data/human_hg38_immune_mixture.cg
yame fetch -c hg38/models/hg38_65celltypes.refx
methscope deconv human_hg38_immune_mixture.cg hg38_65celltypes.refxcell … Macrophage Mono NK_CD56 …
1 … 0.7486 0.2370 0.0000 …Truth 70/30; the other 63 types come back ~0. The 65-type reference spans 58 Zhou single-cell + 7 Loyfer organ/blood types for whole-body / cfDNA deconvolution.
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.ubjx
methscope deconv-build-ref -o self.refx human_hg38_celltypes.cg hg38_celltype.ubjx
methscope deconv human_hg38_celltypes.cg self.refx# Oligodendrocyte Oligodendrocyte=100%
# Pancreas-Beta Pancreas-Beta=100% ... (each 100% itself)--matrix emits the raw record x pattern beta TSV instead of a
.refx.
methscope deconv-build-ref --matrix query.cg 2025_Zhou_MajorPseudo.mrmpImpute dense CpG methylation from a sparse methylome with a bundled decoder — pure-C
inference, no CUDA or BLAS at run time. Output is a .cg of binary
0/1 calls; --probs emits per-CpG probabilities instead.
# 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
yame fetch -c hg38/models/hg38_10k1.updecx
methscope upscale -o human_hg38_test_reconstructed.cg hg38_10k1.updecx human_hg38_test.cg
# what comes out
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 10000 6239 10000 6239 NA 0.624 NAThis decoder writes one 10,000-CpG block, so N_univ is 10,000. Inside that block the input carried 13 covered CpGs; the truth has 9,794 at beta 0.593, and the reconstruction fills all 10,000 at beta 0.624 — it recovers the region's methylation level, not just its average.
Score the calls against the truth (the .cg holds binary 0/1 calls,
so this is exact-match accuracy; --probs emits probabilities
instead).
yame fetch -c hg38/data/human_hg38_test.truth.cg
paste <(yame rowsub -I 1_10000 human_hg38_test_reconstructed.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 human_hg38_40_celltypes_chr20.cg chr20.mrmpsamples=40 CpGs=773477 distinct patterns=116450 (+PNA) candidates=116450# 2. featurize the truth atlas -> sampling/truth msur
methscope upscale-featurize --reps 20 --sample 8000 human_hg38_40_celltypes_chr20.cg chr20.mrmp chr20.msurwrote chr20.msur and chr20.msur.tsv # 40 cells x 20 reps, 91 MB# 3. group CpGs into processing units
methscope upscale-set-units --unit-cpgs 4096 chr20.mrmp 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 8backend=cpu; cells=40 reps=20 patterns=1000 input=2000 units=109 threads=8
units 109/109 (resumed=0) last_val_mae=0.156606
wrote bare UPDEC2 work/model.updec2 (109288672 bytes)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.