🧬 DenseCall2

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

✨ Highlights

📑 Table of Contents

• Installation
• Part 1 — Basecalling
• Part 2 — Training
• Part 3 — Signal Classification
• Part 4 — Genome Language Modelling
• Part 5 — RNA m6A Benchmarking
• Part 6 — DNA CpG Methylation Benchmarking
• 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 .

⚡ 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 — Signal Classification

Signal classification labels each raw nanopore read with a species / barcode / taxon, typically for microbiome abundance profiling. It is driven by the densecall classify subcommand namespace, which wraps the densecall/microbiome/ subpackage.

Available Subcommands

densecall classify label        Build label map / inverse label map from a FASTA
densecall classify split        Split read→label pairs into train / test sets
densecall classify gendata      Convert POD5/FAST5 + labels into train.hdf5 / valid.hdf5
densecall classify train        Train the classification CNN (real or simulated data)
densecall classify predict      Classify reads with a trained model
densecall classify precision    One-class precision / recall / F1 / ROC evaluation
densecall classify evaluation   Full multiclass panel (confusion matrix, PR/ROC, ...)
densecall classify comparison   Side-by-side comparison of two models

Training Modes

The CNN architecture is set by a small TOML (e.g. dna.toml with strides=[2,2,2,2,2] for DNA, rna.toml with strides=[5,2,2,2,2] for RNA) passed via -config. A label_map.json/inverse_label_map.json pair must be produced first with classify label.

🧪 Example: Barcode / Species Classification on Real Reads

# Build label map
densecall classify label species.fa
densecall classify split readid_species.txt --train_file train_labels.txt --test_file test_labels.txt

# Generate training data
densecall classify gendata train.pod5 \
  --label_path train_labels.txt --save_path ./ --seqlen 6000 --do-trim

# Train the model
densecall classify train model_real \
  -config dna.toml --directory ./ --label-map label_map.json \
  --chunksize 6000 -f --new -y

# Predict
densecall classify predict model_real test.pod5 inverse_label_map.json \
  --chunksize 6000 > result_real.txt

# Evaluate
densecall classify evaluation result_real.txt test_labels.txt \
  -m label_map.json -o ./real

🧫 Example: Simulation-Based Training (No Real Per-Species Data)

# Build label map
densecall classify label species.fa

# Train with simulated data
densecall classify train model_simu \
  -config dna.toml --directory ./ --label-map label_map.json \
  --chunksize 25000 -f --new --simu species.fa --preset dna-r10-min -y

# Predict
densecall classify predict model_simu reads_dir/ inverse_label_map.json \
  --chunksize 25000 > result_simu.txt

# Evaluate precision
densecall classify precision result_simu.txt test_labels.txt -m label_map.json

🧬 RNA Direct-Read Classification

RNA reads are reverse-complement mapped and sequenced after the poly(A) tail, so gendata/predict take --rna --tail, and RNA models use the rna.toml config:

# Generate training data
densecall classify gendata train.pod5 \
  --label_path train_labels.txt --save_path ./ --rna --seqlen 50000 --tail

# Train
densecall classify train model \
  -config rna.toml --directory ./ --label-map label_map.json \
  -f --new -y --chunksize 50000 --rna

# Predict
densecall classify predict model test.pod5 inverse_label_map.json \
  --rna --tail --chunksize 50000 > result_real.txt

The trained model/ directory is a normal densecall model directory and can also be used with the basecalling train/predict tooling where applicable.

🧠 Part 4 — 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). 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

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

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

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.

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 5 — 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 6 — 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

📂 Repository Layout

densecall/
├── cli/                  # basecaller, accuracy, convert, sigmap, train, classify, lm, rnalm
├── asr/                  # RNA hybrid CTC + attention-decoder model (densecall.asr)
├── conformer/            # Conformer encoder + mod-call pipeline (densecall.conformer)
├── microbiome/           # signal classification (densecall.classify subcommands)
├── 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 5)
├── 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

💙 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}
}