mdb

mdb builds, merges, queries, and analyzes CpG-by-sample methylation matrices from ONT and PacBio BED inputs. It also provides cohort PCA/UMAP, summary statistics, interactive methylation visualization, ONT ASM PCA, and strand-bias hotspot detection. Cohorts can also be reduced to predefined or data-driven CpG groups for downstream analysis.

• PyPI package: methdb
• CLI command: mdb

Install

pip install methdb

To install from a local checkout instead:

pip install .

Verify:

mdb --help
mdb --version

Core Concepts

• Sample bundle (.smdb): one sample, multiple track views (assay/haplotype/strand).
• Cohort store (.mmdb): merged sample bundles for population-scale queries.
• Grouped cohort store (.mmdb): the same cohort format, with rows reduced to predefined or de novo CpG groups and an additional group-coordinate index.
• Backends:
  • zarr (default, compressed, block-aligned merge writes)
  • npy (optional compatibility backend)

Quick Start

1) Build CpG index

mdb index -r GRCh38_no_alt.fa -o GRCh38.cpg_index.npz

Include chrX/chrY:

mdb index -r GRCh38_no_alt.fa -o GRCh38.cpg_index.npz --sex

2) Create sample bundle

ONT (modkit output file, directory, or prefix):

mdb create \
  -p ont \
  -n GRCh38.cpg_index.npz \
  -b /path/to/ont_input \
  -o sample_ont.smdb \
  -c 5 \
  --sample-id SAMPLE_ONT

Accepted ONT layouts are:

• a BED/BED.GZ file, treated as the combined track;
• a directory containing any of combined.bed(.gz), hp1.bed(.gz), or hp2.bed(.gz); or
• a prefix resolving to <prefix>.combined.bed.gz, <prefix>.hap1.bed.gz, and/or <prefix>.hap2.bed.gz.

When a modkit BED contains both + and - records, mdb create automatically builds separate strand tracks.

PacBio (prefix or a specific BED file):

mdb create \
  -p pacbio \
  -n GRCh38.cpg_index.npz \
  -b /path/to/pacbio_prefix \
  -o sample_pb.smdb \
  -c 5 \
  --sample-id SAMPLE_PB

A PacBio prefix resolves to <prefix>.combined.bed(.gz), <prefix>.hap1.bed(.gz), and <prefix>.hap2.bed(.gz). Pointing directly to a BED/BED.GZ file creates a combined track. Prefix or direct-file input is recommended. Do not rely on PacBio directory discovery in the current implementation; pass the prefix or a specific BED path instead.

3) Merge sample bundles into a cohort

Default backend (zarr):

mdb merge \
  -i sample_ont.smdb sample_pb.smdb \
  -o cohort.mmdb \
  --workers 2 \
  --block-size 64 \
  --zarr-row-chunk 65536 \
  --zarr-codec zstd \
  --zarr-clevel 5 \
  --zarr-shuffle bitshuffle \
  --zarr-codec-threads 4

NPY backend (explicit):

mdb merge \
  -i sample_ont.smdb sample_pb.smdb \
  -o cohort_npy.mmdb \
  --cohort-backend npy \
  --workers 2 \
  --block-size 64

Build modifiedC view (5mC + 5hmC where available):

mdb merge -i sample_ont.smdb sample_pb.smdb -o cohort_modifiedc.mmdb -m

4) Append new samples to existing cohort

mdb append \
  -c cohort.mmdb \
  -i new_sample1.smdb new_sample2.smdb

5) Query values

Point query:

mdb query \
  -i cohort.mmdb \
  --sample-id SAMPLE_PB \
  --assay 5mC \
  --haplotype combined \
  --strand combined \
  --locus chr1:10469

Range query:

mdb query \
  -i cohort.mmdb \
  --sample-id SAMPLE_PB \
  --assay 5mC \
  --haplotype combined \
  --strand combined \
  --region chr1:10469-12000

6) Build a CpG-grouped cohort database

After mdb merge, mdb group writes a separate database containing group-level mean methylation values while retaining every assay, haplotype, strand, and sample column. It never modifies the source .mmdb. Only true multi-CpG groups are materialized by default (--min-group-cpgs 2); singleton CpGs remain available exclusively through their original source-cohort views. The grouped output is marked with "representation": "cpg_groups" in manifest.json.

Use the exact Loyfer atlas-block index consumed by sniffcell find:

mdb group -i cohort.mmdb -o cohort.loyfer.mmdb --grouping loyfer

This loads the packaged, versioned loyfer_grch38_v1.npz index containing 7,990,609 blocks from the exact atlas index used by sniffcell find. Loyfer intervals are G-anchored, so the stored metadata fixes a one-base left shift when mapping them onto the C-anchored .mmdb CpG index. This is the block representation from the Loyfer methylation atlas.

Use the DECODE/Nanopolish definition, which joins adjacent CpGs separated by at most 10 bp:

mdb group -i cohort.mmdb -o cohort.decode.mmdb --grouping decode

This loads the packaged, versioned decode_grch38_10bp_v1.npz index containing 22,179,326 autosomal CpG units. Its parameters are fixed: adjacent CpGs at most 10 bp apart belong to the same unit. This applies the published DECODE/Nanopolish CpG-unit rule to the complete GRCh38 autosomal CpG index; it does not attempt to recreate the paper's cohort-specific Nanopolish and quality-control exclusions. There is no runtime DECODE distance option.

The packaged definitions are immutable inputs to named modes. Set MDB_GROUP_INDEX_DIR or use --predefined-index-dir only to relocate the same versioned files. The file method, version, assembly, parameters, source, and SHA-256 are written into grouping_summary.json, and the mapped group index is copied into the output database.

Or learn groups from the cohort. Adjacent CpGs are joined when their cross-sample methylation profiles meet the correlation threshold and genomic distance limit:

mdb group \
  -i cohort.mmdb \
  -o cohort.denovo.mmdb \
  --grouping denovo \
  --denovo-min-correlation 0.8 \
  --denovo-max-gap-bp 200 \
  --denovo-min-shared-samples 3 \
  --threads 4

The de novo method is adjacent-correlation segmentation. It uses the requested source view (5mC__combined__combined by default), computes Pearson correlation for every neighboring CpG pair using samples observed at both positions, and links a pair when correlation is at least 0.8, distance is at most 200 bp, and at least three samples are shared. Connected runs of passing adjacent links are the groups. A missing/constant/under-supported pair creates a boundary. The learned boundaries are then reused to aggregate every view in the cohort. This is deterministic for a given input and parameter set, but it is a greedy local segmentation rather than a statistical change-point model. --threads scans chromosomes concurrently during de novo discovery and aggregates independent assay/haplotype/strand views concurrently when writing the grouped store.

The ordinary index.npz indexes reduced group rows at their GRCh38 start coordinates, so existing cohort readers continue to work. Each grouped store also includes groups.npz and grouping_summary.json. The additional index maps every group row to its GRCh38 start/end, source CpG row interval, CpG count, and predefined-reference interval. Add --write-group-table only when the redundant, human-readable groups.tsv.gz is useful; omitting it avoids serializing millions of text rows during routine builds.

The summary records the method, parameters, citation/provenance, aggregation rule, achieved reduction, omitted singleton count, and unchanged source cohort. Group values are arithmetic means of observed CpG beta values and require at least --min-observed-cpgs observations per sample.

Grouped outputs retain the normal cohort-store kind and backend, so they can be passed directly to mdb pca, mdb stats, mdb viz, and other cohort readers. Point and range queries use each group's GRCh38 start in index.npz; use groups.npz for the complete group interval and source-CpG mapping, and query the unchanged source cohort when individual-CpG resolution is required.

As a reference benchmark, grouping a 29,152,456-CpG by 294-sample Zarr cohort containing six 5mC/5hmC views produced the following results. These runs also wrote groups.tsv.gz; elapsed time depends on storage and CPU performance.

For the de novo run, group discovery took 108 seconds; the remaining time was spent aggregating and compressing the six output views.

7) Run PCA on cohort view

mdb pca \
  -i cohort.mmdb \
  -o cohort_pca \
  --assay 5mC \
  --haplotype combined \
  --strand combined \
  --n_pcs 10 \
  --frac_cpgs 0.1 \
  --plot_style studio \
  --plot_style_variants \
  --outlier_detect \
  --outlier_alpha 0.999 \
  --outlier_n_pcs 10

When --outlier_detect is enabled, pca also writes:

• outlier_report.tsv (all samples with outlier columns),
• outliers_only.tsv (flagged samples),
• pca_with_outliers_marked.html,
• pca_no_outliers.html,
• pca_pairplot_no_outliers.html,
• pca_pairplot_no_outliers.png.

The no-outlier plots require at least one inlier, and the no-outlier pairplot requires at least two. Plotly PNG exports are written only when a compatible image engine is available.

When --plot_style_variants is enabled, extra style comparison HTML files are written (for example pca_studio.html, pca_sunrise.html, pca_paper.html depending on selected primary style).

Standard PCA outputs include:

• embedding.tsv,
• params.json,
• pca_umap.log,
• pca.html,
• pca_pairplot.html,
• pca_pairplot.png.

Add -m metadata.tsv to make aligned metadata columns available for plot coloring. With metadata coloring enabled, --pairplot-hue selects the initial color field for both pca.html and the pairplot; all usable metadata columns remain available through the interactive dropdown.

For example, a grouped cohort is analyzed exactly like the original cohort:

mdb pca \
  -i cohort.denovo.mmdb \
  -o cohort_denovo_pca \
  -m metadata.tsv \
  --assay 5mC \
  --haplotype combined \
  --strand combined \
  --frac_cpgs 0.001 \
  --n_pcs 10 \
  --pairplot-mode metadata \
  --pairplot-hue tissue_name

Add --umap to calculate a UMAP embedding and write umap.html.

Restrict PCA to CpGs overlapping BED intervals:

mdb pca \
  -i cohort.mmdb \
  -o cohort_pca_regions \
  --cpg-bed regions.bed

Alternatively, average all CpGs within each BED interval and run PCA on the resulting region-by-sample matrix:

mdb pca \
  -i cohort.mmdb \
  -o cohort_pca_region_means \
  --cpg-bed regions.bed \
  --cpg-bed-agg

The aggregation mode also writes region_avg.tsv. Both BED modes require a current cohort store (.mmdb) because legacy flat NPY folders do not contain genomic positions.

8) Summarize cohort stats

mdb stats \
  -i cohort.mmdb \
  -o cohort_stats \
  -m metadata.tsv \
  --assay all \
  --haplotype all \
  --strand all \
  --plot_style studio

stats writes:

• sample_stats.tsv: one row per sample per cohort track with observed CpG counts and fractions.
• track_stats.tsv: per-track summary across assay / haplotype / strand subclasses.
• metadata_group_stats.tsv: per-track grouped summaries for plotted metadata columns.
• cpg_count_scatter.html: interactive sample-rank scatter with PCA-style color dropdown.
• frac_cpg_scatter.html: corresponding sample-rank scatter for the observed CpG fraction.
• cpg_count_by_track.html: interactive box/point plot of observed CpGs stratified by metadata across tracks.
• frac_cpg_by_track.html: corresponding box/point plot for the observed CpG fraction.
• params.json and stats.log: run parameters and detailed logging.

metadata_group_stats.tsv is written only when the aligned metadata provides usable grouping columns. PNG versions of the plots are written when Plotly PNG export is available.

9) Build binned methylation profile HTML

mdb viz \
  -i /path/to/cohort.mmdb \
  -o cohort_viz \
  -m metadata.tsv \
  --assay 5mC,5hmC \
  --haplotype combined \
  --strand combined \
  --bin-length 100000

viz also accepts a directory of .smdb bundles, bundle globs, or a .txt manifest of bundle paths. It writes:

• methylation_viz.html: interactive heatmap HTML for switching tracks, chromosomes, metadata-group rows, and selected sample rows.
• binned_profiles.npz: compressed sample-by-bin matrices for each selected track.
• sample_metadata_aligned.tsv: aligned metadata used by the HTML.
• tissue_name_group_profiles.tsv.gz: precomputed tissue-level mean profiles when tissue_name metadata is available.
• viz_manifest.json: run parameters, resolved inputs, selected tracks, and output paths.
• viz.log: detailed run log.

10) Plot a locus-focused methylation profile from a cohort

mdb plot \
  -i cohort.mmdb \
  -o locus_plot \
  -r chr15:80150000-80200000 \
  -m metadata.tsv \
  --sample-id sampleA,sampleB \
  --assay 5mC,5hmC \
  --combine-tracks sum \
  --window-size 20 \
  --color-by tissue_name

plot is separate from viz: it focuses on one genomic interval from a cohort store and writes an interactive line-profile HTML with controls for:

• selected track view(s) within the requested assay / haplotype / strand filters
• sample-line mode or grouped-mean mode using the chosen metadata color field
• browser-side sliding-window smoothing
• metadata-driven recoloring when aligned metadata is provided

It writes:

• methylation_plot.html: interactive locus-focused profile viewer with track, color, grouping, and window controls.
• region_profiles.npz: compressed CpG-level matrices for the plotted region and selected tracks.
• smoothed_profiles.tsv.gz: smoothed profile coordinates and values using the initial command-line window size.
• sample_metadata_aligned.tsv: aligned metadata used by the HTML controls.
• plot_manifest.json: run parameters and output paths.
• plot.log: detailed run log.

11) Run PCA on ONT ASM segments

asmpca expects modkit DMR segment BEDs containing fields such as name, effect_size, and cohen_h:

mdb asmpca \
  -i sample1.segments.bed.gz sample2.segments.bed.gz \
  -o asm_pca \
  -m metadata.tsv \
  --feature-mode dmr_location \
  --min-region-samples 2 \
  --n-pcs 10

Inputs may be BED paths, directories, globs, or a text manifest of paths. By default, cohort DMR regions are the merged union of input rows whose name is different. Supply --dmr-regions regions.bed to use an external region set.

The default dmr_location mode uses binary DMR presence as the PCA features. To project a segment statistic onto the regions instead:

mdb asmpca \
  -i asm_segments/*.bed.gz \
  -o asm_metric_pca \
  --feature-mode segment_metric \
  --metric effect_size

Core outputs are embedding.tsv, dmr_regions_used.bed, params.json, pca.html, pca_pairplot.html, pca_pairplot.png, and pca_umap.log.

12) Detect strand-biased methylation hotspots

strand requires a cohort store containing matching plus and minus views for the requested assay and haplotype:

mdb strand \
  -i cohort.mmdb \
  -o strand_hotspots \
  -m metadata.tsv \
  --group-by tissue_broad \
  --assay 5hmC \
  --haplotype combined \
  --min-paired-frac 0.8 \
  --min-mean-total 0.005 \
  --cluster-gap-bp 1000 \
  --top-n-hotspots 500 \
  --workers 4

The command scores strand imbalance genome-wide, clusters adjacent same-direction CpGs, and optionally repeats hotspot calling within metadata groups. Core outputs include:

• per_sample_metrics.tsv.gz,
• hotspots_global.tsv and hotspots_global.bed,
• hotspot_sample_profiles.tsv.gz when hotspot profiles are available,
• per-group hotspot TSV/BED and summary files when --group-by is used,
• interactive strand-bias and hotspot HTML reports,
• params.json and strand.log.

Important Notes

• create --reader currently defaults to scan and the active create path uses scan-based reading.
• create --workers is currently accepted by the CLI but is not used by the active scan-based create path.
• merge and append require sample bundles created by current mdb create (manifest-based .smdb layout).
• pca now supports both current cohort stores (.mmdb) and legacy flat merged .npy folders.
• PCA color categories are sorted before plotting so colors stay stable across full vs no-outlier comparisons.

License

MIT (LICENSE).