🧬 DenseCall2

An open platform for nanopore signal classification, basecalling with native modification detection, and methylation-aware genome language modelling.

✨ Highlights

🦠 Signal classification (label reads by species / barcode / taxon) now ships as a separate package: denseclass. It was previously densecall classify; see Signal Classification for the migration note.

📑 Table of Contents

• Installation
• Part 1 — Basecalling
• Part 2 — Training
• Part 3 — Genome Language Modelling
• Part 4 — RNA m6A Benchmarking
• Part 5 — DNA CpG Methylation Benchmarking
• Signal Classification
• Repository Layout
• Supported Subcommands
• Acknowledgements
• Citation

🛠️ Installation

⚠️ Platform Note: DenseCall2 has been tested on Linux Ubuntu. macOS is not supported due to CUDA/NVIDIA toolchain dependencies.

Quick Install

Install from PyPI:

conda create -n densecall python=3.10
conda activate densecall
pip install densecall

Or install from a local checkout during development:

pip install -e .

🛠 CUDA Toolkit

DenseCall2 may need to compile CUDA code at install time and requires the NVIDIA CUDA Toolkit (which provides the nvcc compiler and runtime libraries). Install it via conda (requires conda-forge + nvidia channels):

conda create -n densecall python=3.10
conda activate densecall
conda install -c nvidia cuda-toolkit=12.8
pip install densecall

💡 Make sure the toolkit version matches the CUDA version your PyTorch was built against (e.g. torch==2.9.1 uses CUDA 12.8) and that nvcc is on $PATH before running the install.

⚡ flash-attn (Optional)

The Conformer encoder uses flash-attention for accelerated inference. It is not declared as a dependency because it must be compiled against your exact CUDA / PyTorch build:

pip install flash-attn --no-build-isolation

💡 Fallback: If flash-attn is not installed, the code automatically falls back to PyTorch's built-in SDPA (memory-efficient / eager attention). The pipeline still runs without the flash-attention speedup.

Dependencies: See densecall/requirements.txt for the full pinned set. CTC decoding uses the external fast_ctc_decode package — no Cython compilation needed.

🧬 Part 1 — Basecalling

Basecalling converts raw nanopore signal (POD5) into nucleotide sequences using a hybrid architecture that combines CTC with an attention-based autoregressive decoder.

Decoding Pathways

📦 Model Directories

A model directory contains a config.toml and one or more weights_*.tar checkpoints. The basecaller CLI auto-selects the highest-numbered checkpoint.

⌨️ Command-line Reference

📌 --alignment-threads is only used when the output is a BAM file (i.e., when --reference is provided).

🧪 RNA004 Basecalling with m6A Detection (CTC Pathway)

densecall basecaller rna004 reads.pod5 \
  --recursive --batchsize 32 --rna --mod --max-reads 20000 \
  --chunksize 12000 --overlap 600 > calls.fq

Without --mod, the attention (ASR) decoder is used:

densecall basecaller rna004 reads.pod5 \
  --recursive --batchsize 32 --rna --max-reads 20000 \
  --chunksize 12000 --overlap 0 > calls_asr.fq

🧬 DNA R9/R10 Basecalling

# DNA R9.4.1
densecall basecaller ../models/dna_r9.4.1_hac_CG@v1 flowcell_dir/ \
  --recursive --batchsize 32 --max-reads 20000 \
  --chunksize 12000 --overlap 600 > calls_r9.fq

# DNA R10.4.1
densecall basecaller ../models/dna_r10.4.1_e8_400bps_hac_CG@v1 flowcell_dir/ \
  --recursive --batchsize 32 --max-reads 20000 \
  --chunksize 12000 --overlap 600 > calls_r10.fq

☸️ DNA CpG Methylation Basecalling (Aligned to Reference)

densecall basecaller ../models/dna_r10.4.1_e8_400bps_hac_CG@v1 sample.pod5 \
  --read-ids sample_readids.txt --reference hg38.fa --mod \
  --alignment-threads 12 --batchsize 256 --chunksize 5000 > sample.bam

🧭 RNA004 Basecalling with Reference-Guided Alignment

densecall basecaller rna004 reads.pod5 --recursive \
  --chunksize 5000 --overlap 200 --batchsize 256 --mod --rna \
  --reference hg38.fa --mm2-preset splice --alignment-threads 12 > calls.bam

📏 Basecalling Accuracy Evaluation

densecall accuracy calls.fq reference.fa -d outdir/

🏋️ Part 2 — Training

Training a basecalling / modification model follows a three-stage pipeline:

sigmap → convert → train

Pipeline Overview

1. 🔗 sigmap — Emit Signal Features from BAM + POD5

--bam-and-pod5 takes pairs of <name.sorted.bam> <pod5_dir>. --modified-bases-models takes a remora-style .pt model for modification labels; its motif must be declared with --motif <motif> <offset> --base <code> --alphabet <alphabet>. Use --reverse for direct-RNA (reverse-complement mapped), and --mod_th to set the modification probability threshold.

densecall sigmap \
  --bam-and-pod5 sample.sorted.bam pod5_dir \
  --save_name sample.hdf5 \
  --max-reads 200000 \
  --levels rna004_9mer_levels_v1.txt \
  --modified-bases-models mod_model.pt \
  --motif DRACH 2 --base Y --alphabet AYCGT --reverse --mod_th 0.9

2. 📦 convert — Chunk the HDF5 into Training Data

densecall convert sample.hdf5 m6A/ --chunksize 12000

3. 🏗️ train — Train the Model

train writes checkpoints and a config.toml into <training_directory> (the directory passed to basecaller later). --directory points at the combined chunk directory containing the training/validation chunks. --compile uses torch.compile; --new starts a fresh training run instead of resuming.

densecall train model_0.9 \
  --directory combined/ \
  --batch 32 --epochs 30 --grad-accum-split 2 \
  --no-quantile-grad-clip --lr 0.002 \
  --alphabet NAYCGT --chunks 1000000 \
  --config models/configs/hybrid_transformer.toml \
  --new --compile

The trained model_0.9/ directory (containing config.toml + weights_*.tar) can then be used directly with densecall basecaller model_0.9 reads.pod5 ....

🧠 Part 3 — Genome Language Modelling

Two language-model families are shipped as subcommand namespaces:

🧬 DNA Language Model (densecall lm)

Methylation-aware DNA language modelling on NA12878. The architecture is a conv + transformer backbone whose hyperparameters live in densecall/lm/config.toml.

Pipeline: pretrain → posttrain → predict

1. Pretrain — Masked-LM Training Over the Genome

--vocab standard (ACGTN) trains the plain model; --vocab z trains the methylation-aware variant on a genome whose modified bases are written as the Z token (ATCGNZ); --vocab y does the same for m6A with the Y token (ATCGNY). The masking pool and the FASTA fetch whitelist are both derived from the selected vocabulary, so Y survives instead of being folded into N. Up-sampling of mod-containing windows is controlled by mod_upsampling in the config (falls back to the legacy z_upsampling). Splits are --train-chroms / --val-chroms (default: train chr1–chr21+X, val chr22).

# Standard vocabulary (ACGTN)
densecall lm pretrain --fasta ./data/human/genome.fasta --vocab standard \
  --output-dir ./pretrain_standard --config ./config.toml

# Methylation-aware vocabulary (ATCGNZ)
densecall lm pretrain --fasta ./data/human/genome_methylated.fasta --vocab z \
  --output-dir ./pretrain_z --config ./config.toml

# m6A-aware vocabulary (ATCGNY)
densecall lm pretrain --fasta ./data/human/genome_y.fasta --vocab y \
  --output-dir ./pretrain_y --config ./config.toml

2. Posttrain — DNA-to-BigWig Fine-Tune

--genome is the FASTA, --bigwig/--bigwig-name are functional tracks (RNA-seq, methylation, accessibility), --bed is a split (train/val/test) BED, and --pretrained-backbone points at the backbone_pretrained.pt from step 1. --freeze-backbone / --freeze-layers N freeze weights to fit only the head (or the top N transformer layers).

densecall lm posttrain --genome ./data/human/genome_methylated.fasta \
  --bigwig ./NA12878/ENCFF808QGQ_RNA-seq.bigWig --bigwig-name "RNA-SEQ" \
  --bed ./data/human/splits.bed --pretrained-backbone ./pretrain_z/backbone_pretrained.pt \
  --vocab z --output ./posttrain_z --config ./config.toml

3. Predict — Generate Predicted BigWig / BED Track

densecall lm predict --genome ./data/human/genome_methylated.fasta \
  --checkpoint ./posttrain_z/checkpoint_latest.pt --out-bed ./lmz.bigWig \
  --bigwig ./NA12878/ENCFF808QGQ_RNA-seq.bigWig --bigwig-name NA12878 \
  --fig-name track.png --config ./config.toml

The train subcommand runs the same posttrain model in --mode train (no --checkpoint needed) when you want to train from a bare posttraining entry rather than loading a backbone:

densecall lm train --genome ./data/human/genome.fasta \
  --bigwig track.bigWig --bigwig-name RNA-SEQ --bed ./data/human/splits.bed \
  --vocab standard --output ./train_out --config ./config.toml

🎯 RNA Splice-Site Prediction (densecall rnalm)

Investigates whether an m6A "Y token" improves splice-site prediction over a plain ACGT baseline. A U-Net + transformer model (UNetSplicePredictor) scores splice site and per-tissue usage in a center-cropped window.

Pipeline: data_maker → data_maker_m6a → train → analysis

💡 Pretrained-DNA arms. --dna_vocab lets the RNA model consume the same inputs as a densecall lm checkpoint instead of the legacy 5/6-token encoding. dna_vocab=none keeps legacy ids unchanged; standard/y remap them through a vectorised LUT to the 11-/12-token DNA vocabulary (densecall/rnalm/dna_vocab.py) — the h5 files are never rewritten. Pair it with --vocab_size 11/12, and warm-start from a pretrained backbone with --load_pretrained (backbone learning rate is multiplied by --backbone_lr_mult, default 0.1).

1. Data Generation

data_maker / data_maker_m6a build HDF5 windows of fixed context (2*context_len + center_len, default 10000 bp) with an all-token integer encoding (0=N, 1=A, 2=C, 3=G, 4=T, 5=Y).

data_maker_m6a injects a Y token at filtered m6A sites (cov>5, ratio>10%, strand-matched, from a modkit pileup BED).

densecall rnalm data_maker     --mode train -c dataset/splice_data/configs_15tissue.yaml
densecall rnalm data_maker     --mode test  -c dataset/splice_data/configs_15tissue.yaml
densecall rnalm data_maker_m6a --mode train -c dataset/splice_data/configs_15tissue_m6a.yaml
densecall rnalm data_maker_m6a --mode test  -c dataset/splice_data/configs_15tissue_m6a.yaml

2. Train

Train baseline (--vocab_size 5, no-Y h5s) and y_aware (--vocab_size 6, m6A h5s) with identical flags so analysis can compare them fairly.

⚠️ --seq_len (default 8960) must be divisible by 2 ** num_downsamples (8960 % 128 == 0 with the default num_downsamples=7).

# Baseline (vocab 5)
densecall rnalm train \
  --train_h5 dataset/splice_data/gtex_500_15tis/dataset_train.h5 \
  --test_h5  dataset/splice_data/gtex_500_15tis/dataset_test.h5 \
  --vocab_size 5 --task both --seq_len 8960 --epochs 10 --batch_size 32 \
  --lr 1e-4 --output_dir ./checkpoints/baseline --amp --num_workers 4

# Y-aware (vocab 6)
densecall rnalm train \
  --train_h5 dataset/splice_data/gtex_500_15tis_m6a/dataset_train.h5 \
  --test_h5  dataset/splice_data/gtex_500_15tis_m6a/dataset_test.h5 \
  --vocab_size 6 --task both --seq_len 8960 --epochs 10 --batch_size 32 \
  --lr 1e-4 --output_dir ./checkpoints/y_aware --amp --num_workers 4

# 768-dim arms warm-started from the DNA LM pretrain checkpoints
# (legacy h5 ids remapped to the DNA vocabulary; h5 untouched)
densecall rnalm train \
  --train_h5 dataset/splice_data/gtex_500_15tis/dataset_train.h5 \
  --test_h5  dataset/splice_data/gtex_500_15tis/dataset_test.h5 \
  --vocab_size 11 --dna_vocab standard \
  --load_pretrained /path/to/pretrain_standard/pretrain_best.pt \
  --backbone_lr_mult 0.1 --task both --seq_len 8960 --epochs 30 --batch_size 32 \
  --lr 1e-4 --output_dir ./checkpoints/baseline_pt768 --amp --num_workers 4

3. Analysis

Compare the two checkpoints on their respective test h5s. Writes 7 PDFs + 4 CSVs to --output_dir. By default all metrics are restricted to Y-bearing windows near the splice site; --all_windows removes that restriction. --bootstrap_b (default 30) controls the cluster-bootstrap CI iterations; --no_bootstrap skips them. vocab_size and dna_vocab are read back from each checkpoint's saved args, so mixed 5/6- and 11/12-token arms compare directly. When the bootstrap runs, metrics_summary.csv gains ci_low, ci_high and p_delta_gt0 columns for every site metric — including the F1 rows (f1_acceptor, f1_donor, f1_macro_site).

densecall rnalm analysis \
  --baseline_ckpt checkpoints/baseline/best_model.pt \
  --yaware_ckpt   checkpoints/y_aware/best_model.pt \
  --test_h5_baseline dataset/splice_data/gtex_500_15tis/dataset_test.h5 \
  --test_h5_yaware   dataset/splice_data/gtex_500_15tis_m6a/dataset_test.h5 \
  --device cuda --batch_size 64 --config dataset/splice_data/configs_15tissue.yaml \
  --output_dir ./visualization/15tissue

💡 Quick smoke tests: Add --max_samples 2000 to train, or --max_windows 3000 to analysis. To sanity-check m6A injection, open the m6A h5 and confirm the integer encoding reaches token 5 (assert (X == 5).any()).

🔬 Part 4 — RNA m6A Benchmarking

Standalone evaluation scripts for RNA m6A (N6-methyladenosine) detection, reproducing the analysis of Zou et al. 2025 (Briefings in Bioinformatics).

⚠️ These are plain Python scripts — run them directly (python m6a_benchmark_complete.py ...), they do not import the densecall package.

They take WT vs IVT (in-vitro-transcribed control) modBAMs produced by densecall basecaller --mod (or any modBAM caller), plus a GLORI ground-truth site table, and output the paper's Fig 2–6 data and comparison figures.

m6a_benchmark_complete.py — Per-Caller Full Evaluation

Evaluates one basecaller against the GLORI ground truth and writes a full set of CSVs + figures (threshold sweep, WT-predicted vs GT modification-ratio density, IVT-vs-WT per-site scatter, per-read/per-site distributions, motif-specific performance, per-site detail, all exonic DRACH WT sites).

python m6a_benchmark_complete.py \
    --wt-bam densecall_wt.bam \
    --ivt-bam densecall_ivt.bam \
    --gt-excel 41587_2022_1487_MOESM3_ESM.xlsx \
    --gt-exonic-bed glori_exonic.bed \
    --fasta hg38.fa \
    --name DenseCall2 --fixed-read-thr 0.6 \
    --output-dir ./densecall2_results --calibrate

📌 Requires modkit on $PATH (used to extract full from the modBAMs). --calibrate calibrates the WT vs IVT probability scales before thresholding.

plot_m6a_comparison.py — Side-by-Side Comparison

Plots two callers (e.g. Dorado vs DenseCall2) on their respective evaluation output directories, with a consistent per-model color scheme matching the article layout:

python plot_m6a_comparison.py \
    --dir1 ./dorado_results --name1 Dorado \
    --dir2 ./densecall2_results --name2 DenseCall2 \
    --rec-thr1 0.6 --rec-thr2 0.6 \
    --output-dir ./comparison_figures

🧫 Part 5 — DNA CpG Methylation Benchmarking

Evaluation of DNA CpG (5mC) modification basecalling against per-CpG ground truth derived from whole-genome bisulfite (WGBS) data. A --mod BAM is validated with modkit validate --bam-and-bed against bisulfite-derived CpG truth BEDs.

Ground-Truth CpG BEDs

Ground-truth CpG BEDs are derived from bisulfite data (CpG.gz.bismark.zero.cov.gz), keeping only autosomal CpG sites with sequencing coverage ≥ 10, then tri-labeled by the bisulfite methylation fraction:

modkit validate --bam-and-bed reads the fifth BED column as the per-CpG label (0 = negative, 1 = positive) and compares it to the basecaller's predicted calls, reporting sensitivity / specificity / precision / F1 etc. for each run.

1. 📥 Produce a --mod BAM

densecall basecaller model_dir sample.pod5 \
  --mod --batchsize 128 --chunksize 5000 --reference reference.fa \
  --recursive --alignment-threads 12 > test.sam
samtools sort --write-index test.sam -o test_sort.bam -@ 12

2. ✅ Validate Against the Ground-Truth Tiers

for tier in 10_90_10 10_80_20; do
  modkit validate --bam-and-bed test_sort.bam \
    ../data/sequences/cpg_${tier}.bed -t 12 \
    -o validate_${tier}.txt
done

3. 📊 Per-Site Bisulfite Correlation

For a direct comparison of the predicted methylation fraction against bisulfite at every covered CpG:

modkit pileup test_sort.bam predict_persite.bed \
  --ref reference.fa --preset traditional --region chr22 -t 12 --no-filtering
bedtools intersect -a predict_persite.bed \
  -b CpG.gz.bismark.zero.cov.gz -wa -wb > merge.bed

🦠 Signal Classification

Signal classification — labelling each raw read by species, barcode or taxon for microbiome abundance profiling — ships as its own package since v0.0.2.9.24. It used to be densecall classify, wrapping densecall/microbiome/.

pip install denseclass        # or: pip install -e ../denseclass
denseclass --help

It is trainable on real reads or on signal simulated from reference genomes, and the trained model/ directory is a normal model directory that the basecalling train/predict tooling can also use.

Migration is a command rename — every flag and every subcommand keeps its signature:

See the denseclass README for the full examples.

📂 Repository Layout

densecall/
├── cli/                  # basecaller, accuracy, convert, sigmap, train, lm, rnalm
├── asr/                  # RNA hybrid CTC + attention-decoder model (densecall.asr)
├── conformer/            # Conformer encoder + mod-call pipeline (densecall.conformer)
├── lm/                   # methylation-aware DNA language model (densecall.lm subcommands)
├── rnalm/                # RNA splice-site prediction + m6A Y token (densecall.rnalm subcommands)
├── scripts/              # standalone m6A benchmark tools (run directly, see Part 4)
├── models/configs/       # conformer.toml, hybrid_transformer.toml
└── *.py                  # core modules (io, util, mod_util, reader, ...)

📦 Model weights are not bundled in this repository. Download the trained weights for rna004, dna_r9.4.1_hac_CG@v1 and dna_r10.4.1_e8_400bps_hac_CG@v1 from the release assets (or the storage link provided), and point densecall basecaller <model_dir> reads.pod5 ... at the directory containing config.toml + weights_*.tar.

🔧 Supported Subcommands

Signal classification moved out to the separate denseclass package — see Signal Classification.

💙 Acknowledgements

• The RNA language model dataset generation code (data_maker / data_maker_m6a) and the training datasets are based on the curated data package from ShenLab-Genomics/biombenchmark.
• Parts of the basecalling training and decoding code are derived from nanoporetech/bonito.
• The U-Net architecture used in the DNA and RNA language models is inspired by NTv3 in instadeepai/nucleotide-transformer.

📖 Citation

Please cite the following publication if you use DenseCall2 in your work:

@article{densecall2026,
  title   = {DenseCall2: An open platform for nanopore signal classification, basecalling with native modification detection, and methylation-aware genome language modelling},
  author  = {Linlian and colleagues},
  year    = {2026}
}