mopdepth

A fast BAM/CRAM depth calculation tool written in Crystal, inspired by mosdepth.

This is an experiment to see if well-known tools can be ported to Crystal using “vibe coding”.

Features

• Fast depth calculation for BAM/CRAM files
• Multiple processing modes (fast mode, fragment mode, CIGAR-based)
• Per-base and region-based depth analysis
• BED file support for custom regions
• Plain, gzip, and BGZF BED input
• Explicit FASTA references for CRAM input
• D4 per-base output via the pure-Crystal d4.cr library
• Window-based analysis
• Comprehensive filtering options (MAPQ, fragment length, flags)

Installation

Prerequisites

• Crystal
• hts-lib (for BAM/CRAM support)
• Make

Build from source

git clone https://github.com/kojix2/mopdepth
cd mopdepth
make

The default target always builds an optimized release binary at bin/mopdepth. Use make debug only when a non-optimized development binary is needed.

D4 support is included in the regular build through bio-cr/d4.cr. It requires neither Rust nor the native d4binding library. Run make test-d4 to execute the D4 write/read-back integration tests.

Library API

Applications can use mopdepth without invoking its CLI:

require "mopdepth"

options = Mopdepth::Options.new(
  mode: Mopdepth::Mode::Standard,
  minimum_mapping_quality: 20
)

Mopdepth::Reader.open("sample.bam", threads: 4) do |reader|
  reader.each_coverage(options: options) do |coverage|
    puts "#{coverage.region.chromosome}: #{coverage.values.size} bases"
  end
end

Mopdepth::Reader owns the indexed BAM/CRAM handle and exposes coverage, coverage_into, and each_coverage. Coordinates are always zero-based and half-open. coverage_into reuses a caller-owned Array(Int32); returned arrays contain one value per base and never expose the internal difference-array sentinel. A reader is serial: use one reader per concurrent worker.

Long operations support cooperative cancellation and progress callbacks:

control = Mopdepth::Control.new { |progress| puts progress }
# Another fiber or thread may call control.cancel.
coverage = reader.coverage("chr1", options: options, control: control)

Usage

./mopdepth [options] <prefix> <BAM-or-CRAM>

Basic example

./mopdepth output sample.bam

Options

• -t, --threads THREADS: BAM decompression threads
• -c, --chrom CHROM: Restrict to chromosome
• -b, --by BY: BED file or numeric window size
• -f, --fasta FASTA: FASTA reference for CRAM input (defaults to REF_PATH)
• -n, --no-per-base: Skip per-base output
• --d4: Write per-base depth as D4
• -Q, --mapq MAPQ: MAPQ threshold
• -l, --min-frag-len MIN: Minimum fragment length
• -u, --max-frag-len MAX: Maximum fragment length
• -x, --fast-mode: Fast mode (read start/end positions only)
• -a, --fragment-mode: Count full fragment (proper pairs only)
• -m, --use-median: Use median for region stats instead of mean
• -q, --quantize QUANTIZE: Write quantized output (for example, 0:1:4:)
• -T, --thresholds THRESHOLDS: Comma-separated thresholds for region coverage
• -F, --flag FLAG: Exclude reads with FLAG bits set
• -i, --include-flag FLAG: Include only reads with FLAG bits set
• -R, --read-groups GROUPS: Comma-separated read group IDs
• -M, --mos: Use mosdepth-compatible filenames (mosdepth.); default is mopdepth.
• -v, --version: Show version
• -h, --help: Show help message

Processing modes

• Default mode: CIGAR-based depth calculation (most accurate)
• Fast mode (-x): Uses read start/end positions (faster but less accurate)
• Fragment mode (-a): Counts full fragments for paired-end reads

Note: Fast mode and fragment mode cannot be used together.

Window sizes must be positive. BED intervals are validated against the alignment header: coordinates must be 0-based half-open, non-empty, non-negative, and contained in the named reference. References absent from the alignment header are reported and ignored. Overlapping BED intervals remain independent observations.

With -n -b regions.bed, chromosomes absent from the BED are skipped entirely when neither -q nor -T is requested. They are therefore absent from the chromosome summaries and global distributions, and do not contribute to total. For chromosomes present in the BED, total and the global distribution still cover each entire processed chromosome (or the selected --chrom range), while total_region and the region distribution cover the BED intervals. This skip is a mosdepth-compatible performance optimization.

--chrom also accepts a 1-based inclusive range such as chr1:100-200. Unlike mosdepth 0.3.x, mopdepth applies that range to calculation and output. In fragment mode it still examines the complete chromosome so fragments spanning the selected range are not missed.

Output files

• Summary: <prefix>.(mopdepth|mosdepth).summary.txt
• Per-base: <prefix>.per-base.bed.gz (unless -n)
• D4 per-base: <prefix>.per-base.d4 (--d4, replacing per-base BED.gz)
• Global dist: <prefix>.(mopdepth|mosdepth).global.dist.txt
• Regions: <prefix>.regions.bed.gz (when --by)
• Region dist: <prefix>.(mopdepth|mosdepth).region.dist.txt (when --by)
• Quantized: <prefix>.quantized.bed.gz (when --quantize)
• Thresholds: <prefix>.thresholds.bed.gz (when --thresholds and --by)

By default, files are named with the mopdepth.* label. Use -M/--mos to switch to mosdepth.*.

Histogram-based outputs have bounded depth bins, matching mosdepth's performance conventions: -m reports 65,535 for medians above 65,535; numeric-window region.dist places rounded means or medians above 511 in bin 511; and global or BED-region distributions place per-base depths above 400,000 in bin 399,990. These bounds do not clip per-base depth, threshold counts, quantized output, or mean and maximum depth in the summary.

mosdepth compatibility

mopdepth follows mosdepth's coverage, filtering, NoData, threshold, quantization, and window-distribution conventions where those conventions are well-defined. Compressed output and CSI files are compared by decoded content and query results, not by their binary bytes.

Two known mosdepth edge-case errors are intentionally not reproduced: contained read pairs use the true intersection when removing mate overlap, and quantization examines the final real base of every reference. Region means are accumulated as integer depth sums before formatting to avoid floating-point accumulation drift.

D4-enabled builds support mosdepth-style --d4 per-base output and create an embedded secondary frame index after successfully closing the file. -M/--mos changes compatible output names; it does not enable bug-for-bug compatibility.

Summary file format

The summary file contains the following columns:

• chrom: Chromosome name
• length: Chromosome length
• bases: Total depth (sum of all depths)
• mean: Mean depth
• min: Minimum depth
• max: Maximum depth

Examples

Basic depth calculation

./mopdepth output sample.bam

With BED regions

./mopdepth -b regions.bed output sample.bam

Window-based analysis (1kb windows)

./mopdepth -b 1000 output sample.bam

Fast mode with MAPQ filtering

./mopdepth -x -Q 20 output sample.bam

Fragment mode for paired-end data

./mopdepth -a -l 100 -u 1000 output sample.bam

License

MIT License