Skip to content

Build a Custom MDB Atlas from bedMethyl

Jump to bottom
Fu-Yilei edited this page Aug 27, 2026 · 1 revision

Build a Custom SniffCell MDB Atlas from bedMethyl Files

This tutorial builds a cell-type methylation atlas from multiple ONT or PacBio bedMethyl files using MDB, then uses the atlas with sniffcell find. MDB stores 5mC and 5hmC as separate assay views, so SniffCell can discover ctDMRs that are hidden when both modifications are collapsed into total modified C.

Requirements

  • All samples use the same reference assembly and chromosome naming.
  • One bedMethyl file per purified/reference sample.
  • Unique sample IDs and a sample-to-cell-type table.
  • For donor-aware filtering, donor and cell-type labels for every sample.
  • mdb, sniffcell, modkit, and tabix available in the environment.

Install the current MDB source:

git clone https://github.com/Fu-Yilei/mdb.git
cd mdb
pip install .

Verify both CLIs:

mdb --help
sniffcell find --help

1. Generate modification-resolved bedMethyl files

If bedMethyl files already exist, skip this step. For ONT BAMs containing 5mC and 5hmC calls, keep the modifications separate:

modkit pileup sample.bam sample.combined.bed.gz \
  --reference GRCh38_no_alt.fa \
  --modified-bases 5mC 5hmC \
  --cpg \
  --combine-strands \
  --bgzf \
  --bgzf-threads 4 \
  --threads 16

tabix -p bed sample.combined.bed.gz

Do not use --combine-mods for a dual 5mC/5hmC atlas. In modkit bedMethyl, column 4 should identify m rows as 5mC and h rows as 5hmC. MDB detects these rows and writes separate assay views.

2. Build the reference CpG index once

Use the exact FASTA used for pileup:

mdb index \
  --reference GRCh38_no_alt.fa \
  --output GRCh38.cpg_index.npz

Add --sex if chrX and chrY should be included. Every sample bundle merged into a cohort must use the same index.

3. Convert every bedMethyl into an MDB sample bundle

mdb create \
  --platform ont \
  --npz GRCh38.cpg_index.npz \
  --bed sample.combined.bed.gz \
  --output sample.smdb \
  --sample-id SAMPLE_ID \
  --min_coverage 5

Use --platform pacbio for PacBio bedMethyl. Repeat this command for every purified/reference sample. For an HPC array, cap concurrent tasks explicitly; for example:

#SBATCH --array=1-100%30

Make a one-path-per-line bundle manifest:

/path/neuron_donor1.smdb
/path/oligo_donor1.smdb
/path/neuron_donor2.smdb
/path/oligo_donor2.smdb

4. Merge sample bundles while preserving modifications

mdb merge \
  --inputs sample_bundles.txt \
  --output reference_samples.mmdb \
  --cohort-backend zarr \
  --workers 4

Do not add -m/--modifiedc when building a dual atlas. That option creates a summed modifiedC = 5mC + 5hmC view and is intended only for a total-modified-C control.

5. Aggregate CpGs into SniffCell/Loyfer blocks

sniffcell find uses the Loyfer block coordinate system:

mdb group \
  --input reference_samples.mmdb \
  --output reference_samples.loyfer.mmdb \
  --grouping loyfer \
  --min-observed-cpgs 2 \
  --min-group-cpgs 2 \
  --threads 8

The grouped database contains groups.npz, which defines the genomic rows used by sniffcell find.

6A. Use the new samples as a standalone atlas

Create celltypes.json with exact MDB sample IDs:

{
  "brain_reference": {
    "Neuron": ["neuron_donor1", "neuron_donor2", "neuron_donor3", "neuron_donor4"],
    "Oligodendrocyte": ["oligo_donor1", "oligo_donor2", "oligo_donor3", "oligo_donor4"]
  }
}

The grouped reference_samples.loyfer.mmdb can now be supplied directly to sniffcell find. Set:

atlas_path=reference_samples.loyfer.mmdb

6B. Optionally extend the legacy Loyfer atlas

To retain an existing SniffCell/Loyfer NPY atlas as unresolved modifiedC while adding modification-resolved samples:

mdb atlas \
  --legacy-npy atlas/all_celltypes_blocks.npy \
  --legacy-index atlas/all_celltypes_blocks.index.gz \
  --legacy-samples atlas/all_celltypes.txt \
  --legacy-assay modifiedC \
  --cohorts reference_samples.loyfer.mmdb \
  --output combined_loyfer_ont.mmdb \
  --cohort-backend zarr

The output keeps separate views:

  • modifiedC__combined__combined: legacy bisulfite signal;
  • 5mC__combined__combined: modification-resolved 5mC;
  • 5hmC__combined__combined: modification-resolved 5hmC.

mdb atlas requires each grouped custom row to match the legacy Loyfer row system exactly. Alignment totals and input provenance are written to atlas_manifest.json. For the following commands, set:

atlas_path=combined_loyfer_ont.mmdb

7. Call single-channel ctDMRs

sniffcell find \
  --mdb "$atlas_path" \
  --assay 5hmC \
  --celltypes_file celltypes.json \
  --celltypes_keys brain_reference \
  --output brain.5hmC.ctdmr.tsv \
  --diff_threshold 0.40 \
  --min_rows 1 \
  --min_cpgs 3 \
  --max_gap_bp 3000

Use --assay 5mC or --assay modifiedC for the other views.

8. Call 5mC and 5hmC together without collapsing them

--assay dual runs separate 5mC and 5hmC calls and concatenates the results with a modification column:

sniffcell find \
  --mdb "$atlas_path" \
  --assay dual \
  --celltypes_file celltypes.json \
  --celltypes_keys brain_reference \
  --output brain.dual.ctdmr.tsv \
  --diff_threshold 0.40 \
  --min_rows 1 \
  --min_cpgs 3 \
  --max_gap_bp 3000

The two signals are never summed or averaged. The downstream locus identity is (chr, start, end, modification), not coordinates alone.

9. Recommended paired-donor filtering for brain 5hmC

Brain 5hmC can vary substantially among individuals and disease states. For paired neuron/oligodendrocyte reference samples, create metadata.tsv:

sample_id	donor	cell_type
neuron_donor1	donor1	Neuron
oligo_donor1	donor1	Oligodendrocyte
neuron_donor2	donor2	Neuron
oligo_donor2	donor2	Oligodendrocyte
neuron_donor3	donor3	Neuron
oligo_donor3	donor3	Oligodendrocyte
neuron_donor4	donor4	Neuron
oligo_donor4	donor4	Oligodendrocyte

Then add:

  --paired-stability-filter \
  --paired-metadata metadata.tsv \
  --paired-min-donors 4 \
  --paired-support-effect 0.30 \
  --paired-min-support 3 \
  --paired-min-effect 0.15

With --diff_threshold 0.40, this requires:

  • absolute median paired neuron-minus-oligodendrocyte effect at least 0.40;
  • the same effect direction in every complete donor pair;
  • at least three donors with absolute effect at least 0.30;
  • no donor with absolute effect below 0.15.

These are optional dual/paired-mode safeguards, not global SniffCell defaults.

10. Validate the atlas before using it

Check the following:

  1. atlas_manifest.json lists the intended assay views and sample counts.
  2. Every sample ID in celltypes.json matches columns.tsv exactly.
  3. 5mC and 5hmC appear as separate views.
  4. No donor or cell type is missing from the paired metadata.
  5. ctDMR tables contain the modification and atlas_source columns.
  6. Inspect effect distributions and donor-level agreement before using new 5hmC ctDMRs for deconvolution.

The legacy --npy/--index/--meta SniffCell interface remains available for reproducibility, but MDB is recommended for new multi-modification atlases.

Clone this wiki locally