Metadata-Version: 2.4
Name: densecall
Version: 0.0.2.9.9
Home-page: https://github.com/n-damo/densecall
Author: linlian
Author-email: 21620151153308@stu.xmu.edu.cn
License: MIT
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Operating System :: POSIX :: Linux
Classifier: License :: OSI Approved :: MIT License
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
Requires-Dist: ont-bonito==1.0.0
Requires-Dist: torch==2.9.1
Requires-Dist: sageattention==1.0.6
Requires-Dist: torchaudio==2.9.1
Requires-Dist: torchvision==0.24.1
Requires-Dist: ont-koi==0.6.5
Requires-Dist: Cython==3.2.9
Requires-Dist: setuptools==77.0.3
Requires-Dist: numpy==1.23.5
Requires-Dist: pyranges==0.1.4
Requires-Dist: pybigwig==0.3.25
Requires-Dist: pyfaidx==0.9.0.4
Requires-Dist: torchmetrics==1.9.0
Requires-Dist: mappy
Requires-Dist: fitter==1.7.1
Requires-Dist: gtfparse==2.0.1
Requires-Dist: pyslow5==1.3.1
Requires-Dist: blue-crab==0.4.0
Requires-Dist: toml==0.10.2
Requires-Dist: tqdm<5,>4
Requires-Dist: scipy==1.14.1
Requires-Dist: pysam==0.21.0
Requires-Dist: parasail==1.3.4
Requires-Dist: pandas==1.5.3
Requires-Dist: requests~=2.32.2
Requires-Dist: ont-fast5-api==4.1.3
Requires-Dist: fast-ctc-decode==0.3.5
Requires-Dist: python-dateutil==2.8.2
Requires-Dist: pytorch-ranger==0.1.1
Requires-Dist: dtaidistance==2.3.11
Requires-Dist: torch_pruning==1.3.7
Requires-Dist: ruptures==1.1.9
Requires-Dist: tslearn==0.6.3
Requires-Dist: patchworklib==0.6.3
Requires-Dist: plotnine==0.12.4
Requires-Dist: matplotlib==3.10.8
Requires-Dist: seaborn==0.13.2
Requires-Dist: biopython==1.86
Requires-Dist: igv==0.9.6
Requires-Dist: deap==1.4.1
Requires-Dist: tabulate==0.9.0
Requires-Dist: onnx-simplifier==0.4.36
Requires-Dist: onnx2torch==1.5.14
Requires-Dist: onnxconverter-common==1.14.0
Requires-Dist: onnxmltools==1.12.0
Requires-Dist: onnxruntime-gpu==1.22.0
Requires-Dist: onnxscript==0.5.1
Requires-Dist: memory-efficient-attention-pytorch==0.1.6
Requires-Dist: edlib==1.3.9
Requires-Dist: ont-remora==3.3.0
Requires-Dist: rotary-embedding-torch==0.8.3
Requires-Dist: debugpy==1.8.8
Requires-Dist: statsmodels
Requires-Dist: transformers
Requires-Dist: uncalled4
Requires-Dist: lhotse
Requires-Dist: nltk==3.9.1
Requires-Dist: muon_optimizer==0.1.0
Requires-Dist: accelerate
Requires-Dist: flashlight==0.1.1
Requires-Dist: flashlight-text==0.0.7
Requires-Dist: pyabf
Requires-Dist: hyperopt
Requires-Dist: tensorboard
Dynamic: author
Dynamic: author-email
Dynamic: classifier
Dynamic: description
Dynamic: description-content-type
Dynamic: home-page
Dynamic: license
Dynamic: requires-dist
Dynamic: requires-python

<div align="center">

# 🧬 DenseCall2

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

<table align="center" style="border: none;">
<tr style="border: none;">
<td style="border: none; padding: 0 8px;">
  <a href="https://pypi.org/project/densecall">
    <img src="https://img.shields.io/pypi/v/densecall" alt="PyPI version">
  </a>
</td>
<td style="border: none; padding: 0 8px;">
  <a href="https://pypi.org/project/densecall">
    <img src="https://img.shields.io/pypi/pyversions/densecall" alt="Python version">
  </a>
</td>
<td style="border: none; padding: 0 8px;">
  <img src="https://img.shields.io/badge/License-MIT-blue.svg" alt="License: MIT">
</td>
<td style="border: none; padding: 0 8px;">
  <img src="https://img.shields.io/badge/platform-Linux%20x86_64-lightgrey.svg" alt="Platform: Linux x86_64">
</td>
</tr>
</table>

</div>

---

## ✨ Highlights

| Feature | Description |
|:--|:--|
| 🧬 **Basecalling + Modification Detection** | Hybrid CTC + attention (ASR) architecture. CTC pathway emits per-read modification probabilities as BAM MM/ML tags; attention decoder delivers higher canonical accuracy with key-value caching for long reads. |
| 🦠 **Signal Classification** | Labels each raw read by species / barcode / taxon for microbiome abundance profiling. Trainable on real reads or signal simulated from reference genomes. |
| 🧠 **Genome Language Modelling** | Methylation-aware DNA language model (`densecall lm`) and RNA splice-site predictor with optional m6A "Y token" (`densecall rnalm`). |
| 🔬 **Modification Benchmarking** | Standalone evaluation scripts for RNA m6A (against GLORI ground truth) and DNA CpG (against whole-genome bisulfite data). |

---

## 📑 Table of Contents

- [Installation](#-installation)
- [Part 1 — Basecalling](#-part-1--basecalling)
- [Part 2 — Training](#️-part-2--training)
- [Part 3 — Signal Classification](#-part-3--signal-classification)
- [Part 4 — Genome Language Modelling](#-part-4--genome-language-modelling)
- [Part 5 — RNA m6A Benchmarking](#-part-5--rna-m6a-benchmarking)
- [Part 6 — DNA CpG Methylation Benchmarking](#-part-6--dna-cpg-methylation-benchmarking)
- [Repository Layout](#-repository-layout)
- [Supported Subcommands](#-supported-subcommands)
- [Acknowledgements](#-acknowledgements)
- [Citation](#-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:

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

Or install from a local checkout during development:

```bash
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:

```bash
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

| Pathway | Flag | Use Case |
|:--|:--|:--|
| **CTC** | `--mod` | Frame-level sequence alignment; emits per-read modification probabilities as BAM MM/ML tags. Best for simultaneous basecalling + modification detection. |
| **Attention (ASR)** | *(default)* | Higher canonical accuracy with key-value caching for low-latency long reads. |

### 📦 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.

| Model | Architecture | Alphabet | Chunksize | Overlap |
|:--|:--|:--|--:|--:|
| `models/rna004/` | RNA004 hybrid (CTC + ASR) | `NAYCGT` | 12000 | 0 |
| `models/dna_r9.4.1_hac_CG@v1/` | DNA R9.4.1 Conformer | `NACZGT` (Z = methylated C) | 12000 | 600 |
| `models/dna_r10.4.1_e8_400bps_hac_CG@v1/` | DNA R10.4.1 Conformer | `NACZGT` (Z = methylated C) | 12000 | 600 |

### ⌨️ 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)

```bash
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:

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

#### 🧬 DNA R9/R10 Basecalling

```bash
# 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)

```bash
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

```bash
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

```bash
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

| Stage | Command | Description |
|:--|:--|:--|
| **1. sigmap** | `densecall sigmap` | Maps POD5 reads to a reference using an existing basecall + minimap2 alignment. Writes per-read signal features (with modification states from a remora-style model) to an HDF5 file. |
| **2. convert** | `densecall convert` | Splits the HDF5 into fixed-size chunks (with a held-out validation set) in a training directory. |
| **3. train** | `densecall train` | Trains the architecture declared in a `models/configs/*.toml` on the chunked data. Uses the `muon` optimizer by default. |

### 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.

```bash
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

```bash
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.

```bash
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

| Mode | Description |
|:--|:--|
| **Real-data** | `gendata` extracts signal chunks from reads labeled by species, then `train` (no `--simu`) fits the CNN directly. |
| **Simulation** | `train --simu <refs.fa> --kmer-model <levels>` synthesizes chunks on the fly from a FASTA of reference genomes plus a k-mer level table. Presets: `dna-r9-min`, `dna-r10-min`, `rna-r9-min`, `rna004-min`, ... |

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

```bash
# 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)

```bash
# 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:

```bash
# 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:

| Namespace | Description |
|:--|:--|
| **`densecall lm`** | Methylation-aware DNA language model (`pretrain`, `posttrain`, `train`, `predict`). |
| **`densecall rnalm`** | RNA splice-site prediction with optional m6A **"Y token"** (`data_maker`, `data_maker_m6a`, `train`, `analysis`). |

### 🧬 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).

```bash
# 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).

```bash
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

```bash
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:

```bash
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.

| Vocabulary | Token Set | `vocab_size` |
|:--|:--|--:|
| Baseline | ACGTN | 5 |
| m6A-aware | ACGTNY | 6 |

**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).

| Mode | Chromosomes |
|:--|:--|
| **test** | chr1, 3, 5, 7, 9 |
| **train** | All others |

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

```bash
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`).

```bash
# 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.

```bash
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).

```bash
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:

```bash
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:

| Tier | Methylated | Unmethylated | Ambiguous |
|:--|:--|:--|:--|
| `*_10_90_10.bed` | ≥ 90% | < 10% | 10–90% (discarded) |
| `*_10_80_20.bed` | ≥ 80% | < 10% | 10–80% (discarded) |

`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

```bash
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

```bash
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:

```bash
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

| Subcommand | Purpose |
|:--|:--|
| `basecaller` | Basecall POD5/FAST5 to FASTQ/BAM/SAM (CTC or attention decoder, `--mod` for modification calls) |
| `accuracy` | Alignment-based basecalling accuracy evaluation |
| `convert` | Chunk HDF5 signal datasets into training data (`sigmap` output) |
| `sigmap` | Map reads to reference and emit per-read signal features with modification states |
| `train` | Train basecalling / modification models (`sigmap` → `convert` → `train` pipeline) |
| `classify` | Signal classification (`label` / `split` / `gendata` / `train` / `predict` / `precision` / `evaluation` / `comparison`) |
| `lm` | Methylation-aware DNA language model (`pretrain` / `posttrain` / `train` / `predict`) |
| `rnalm` | RNA splice-site prediction with optional m6A Y token (`data_maker` / `data_maker_m6a` / `train` / `analysis`) |

---

## 💙 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](https://github.com/ShenLab-Genomics/biombenchmark).
- Parts of the basecalling training and decoding code are derived from [nanoporetech/bonito](https://github.com/nanoporetech/bonito).
- The U-Net architecture used in the DNA and RNA language models is inspired by NTv3 in [instadeepai/nucleotide-transformer](https://github.com/instadeepai/nucleotide-transformer).

---

## 📖 Citation

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

```bibtex
@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}
}
```
