Roadmap

August 14, 2026 · View on GitHub

checkout, refine and assemble all work. Implemented: the .mig intermediate format (reader, writer, range partitioning, CRC-checked blocks, provenance and quality-calibration in the header), FASTQ IO for plain and gzipped input with strict validation, barcode packing/unpacking and the IUPAC/Phred primitives, the pybind11 module, the read simulator with ground truth, CI (C++ on ubuntu+macos, Python 3.10/3.12 matrix, pinned ruff), sphinx docs with a zero-warning gate, and the PyPI publish workflow.

migec checkout works: degenerate barcode patterns in the MIGEC dialect read verbatim from published barcode tables, quality-aware log-likelihood acceptance, sample assignment with an ambiguity verdict distinct from "unmatched", UMI extraction from non-contiguous runs, trimming of adapter/tag/UMI, barcode transfer into SAM-style FASTQ headers, the power-of-two coverage histogram, per-position base composition with Shannon entropy and information content, collision entropy and effective barcode length, data-driven UMI error-rate estimation, and count correction with a sequencing/polymerase/independent-molecule mixture plus the collision-corrected molecule count. scripts/spikein_ratio.py computes the published spike-in validation metric, and scripts/compare_calib.py scores UMI grouping against Calib by adjusted Rand index with splitting and merging reported separately.

Nothing in the pipeline scales with the library any more. checkout's counters and refine's barcode table both range-partition themselves past a byte budget, and correction follows both into the partition in two passes with the key rotated.

Throughput and footprint: every stage threads and every stage is byte-identical at any -t -- checkout 1.55 M reads/s end to end (1.70 M matching), refine 1.55 M, assemble 2.47 M, all at 16 threads; ~22 bytes per distinct UMI against a hash map's ~48. Verified under the thread sanitizer as well as by comparison at 1..16 threads. The UMI counters range-partition themselves past a 1 GB budget, and correction follows them into the partition in two passes with the key rotated, so nothing in checkout scales with the library any more. checkout --mig writes the reads into that same partition, which is the pass assemble otherwise builds for itself: 1.16 s -> 0.98 s end to end on 500 k reads for identical molecules.

QC: twenty gnuplot panels over the tables the stages write, including the four figures a user already knows how to read -- Cell Ranger's barcode rank plot on unique UMIs, the MIG size spectrum with molecules and reads, the rank/Zipf curve, and consensus quality as a box over an exact (depth, quality) grid.

Milestones are ordered by risk, not by pipeline order: the consensus quality model is the scientific claim and is validated before any throughput work.

What is left, in the order it should be done

Everything below is either open or half-open. Nothing else on this page is.

#	item	milestone	why it is next
1	The RT-vs-first-cycle-PCR split inside the floor	M3	the sequencing/pre-amplification split now has a standard — the deep-MIG consensus residual, which is what `--pre-amp-error auto` fits (`docs/quality_floor.rst`). What is still unsplit is what MADE the floor, and that needs two chemistries to compare rather than a better estimator: the same template through an RT protocol and through a DNA one. Data, not code
2	The rest of the published comparisons	M5	MIGEC v1, MAGERI, UMIErrorCorrect, UMI-tools, fgbio and Cell Ranger are done and scored (`docs/postprocessing.rst`, `docs/validation.rst`, `docs/grouping.rst`, `docs/single_cell.rst`). Open: UMI-VarCal, which cannot run on the certified arms (single-end). This is what the version number is waiting on, not the code
3	The last two rows of the ctDNA table	M5	the ground truth is found and scored: `PRJNA788522` / `PRJNA507366`, certified VAF, real 12 nt inline UMI. Done: LoFreq and Mutect2 on the identical consensus BAM, a MAPQ 20 floor, adapter trimming, and UMIErrorCorrect end to end (`assets/ctdna_callers.tsv`). Running: the no-consensus baseline -- same reads, same aligner, same callers, reads instead of molecules -- and UMI-VarCal
4	Bit-parallel matcher	M2	last, deliberately: the scan is O(offsets x pattern) and is not the bottleneck. It goes in when a benchmark says so

Nothing is blocked on data any more. Britanova et al ageing (bulk TCR, shallow -- the real 1-3 reads/UMI dataset) has been run: one HiSeq lane, ten donors, 149,588,907 read pairs at a 16 nt UMI and 2.44-2.57 reads per UMI, in assets/shallow_repertoire.tsv and docs/validation.rst. The ctDNA ground truth was found rather than built, by screening SRA read structure instead of trusting the published claim that none exists (scripts/sra_fetch.py probe).

M0 — skeleton, format, simulator

Archive the Groovy implementation on legacy-v1 / v1-final, start master fresh
.mig format frozen in docs/formats.rst, with a round-trip and truncation test
FASTQ reader/writer; a record straddling a buffer refill is covered by a test
Read simulator with truth files, and tests for the simulator itself
CMake + scikit-build-core + CI + docs shell
Push to origin

X — experiments that must precede the milestones they inform

X1 — read-start dispersion within (CB,UMI) on one 10x run. Done, 2026-08-13, on pbmc_1k_v3 at MALAT1/ACTB/B2M — 8,326 multi-read groups. The co-terminal assumption is false: 7.8% of groups overall, and 0.3% of those with ≥6 reads, so what co-terminality there is is a small-group coincidence rather than a property of the chemistry. 92% of groups are wider than one 91 nt read. But 72.7% still form a single overlap component (rising to 81.3% at ≥6 reads), so partitioning by overlap first reduces three quarters of the work to the ungapped problem MIGEC already solves — and the remaining 27.3% is why the partition is mandatory: a single consensus over them asserts sequence across a gap no read covers. Also: only 1.5% of groups hold more than one read at this depth, so for shallow 3' GEX the UMI buys counting, not error correction. Written up in docs/fragmented.rst, script in scripts/read_start_dispersion.py.
X2 — emitted-quality calibration on a clonal control, stratified by MIG size. Done, 2026-08-13, on SRR1763769 (2.12 M reads, HIV-1 Primer ID). The floor is of order 1e-4, not 1e-6: 1.54e-4 [1.36e-4, 1.74e-4] at MIGs of ≥80 reads, so no emitted quality above Q40 is supportable by default, and a blanket 1e-6 is excluded for an RT protocol by two orders of magnitude. Matches the ~1 in 10,000 the source paper reports (doi:10.1128/JVI.00522-15). Note: The curve is still declining at 80 reads, so this is an upper bound: the 9 nt Primer ID puts the library at 49.6% occupancy and checkout flags it saturated, so collided MIGs contribute mismatches counted here as error. Not the p_floor + a/c least-squares fit the plan specified — that model is wrong for a majority vote and returned a negative probability on simulated data. docs/quality_floor.rst, scripts/quality_floor.py.
X3 — three permutation nulls on one deep real dataset. Done, 2026-08-13, on SRR1763769 — 125,369 distinct 9 nt barcodes at 47.8% occupancy. (1) Position independence holds to ~1% of the collision rate. Tested as a distribution — JSD against the product measure q(u) = Π_j p_j(u_j), floored by a same-size draw from q — not by one functional. Invisible on all distinct barcodes (saturation), clear at reads ≥ 2 (z up to 33), and entirely nearest-neighbour: every adjacent pair positive, every distant pair zero. Cause measured: 0.55% of reads carry a barcode one base short, a coupling that did not fire, which frameshifts everything after it. Π_j m_j stays; the 1.86x collision excess is the read threshold. Note: The first version reported 1.04x and all of it was artefact — N as a fifth base, and the plug-in Σ p̂² whose bias grows with k. (2) 97% of distance-1 pairs are chance (844,243 observed, 817,358 under a column shuffle), and a size-preserving count shuffle over the fixed graph finds ~18,000 genuine error children, plateauing from a count ratio of 5 upward. The permutation background puts the barcode error at 1.4e-3, 0.70x of the Phred + polymerase prediction, against checkout's analytic 8.0e-4 at 0.39x — M3 takes the permuted background. (3) The split threshold is 8.68, not 2.00: a curveball randomisation preserving both margins of the reads x positions minor-allele matrix puts the 1% false-positive point 19x above the nominal p < 0.01, because a low-quality read carries minor bases at many positions at once and mimics a linked subclone. Note: bootstrap 95% CI [8.42, 9.14]** over 82,800 randomisations; a tenth as many gave 9.61 and 11.66, so the interval is the number. (4) Shallow libraries (1-3 reads/UMI) are a separate regime and three of these results do not transfer to them — written up rather than quoted. docs/nulls.rst, scripts/permutation_nulls.py.

M1 — `assemble`, the consensus and quality model

M2 — `checkout`

M3 — `refine` (the stage works; cell calling and the FDR threshold are open)

Measured before building (scripts/correction_accuracy.py, 2026-08-13): the existing count-ratio correction works from 3.1 reads/UMI upward (recall ≥0.8, precision ≥0.9) and collapses below it — 0.02 recall / 0.25 precision at 1.1 reads/UMI. ε follows: 0.98x of injected at 3.1 reads/UMI, 0.20x at 1.1. Molecules are never lost (≥0.99 kept at every depth), so the failure is missed correction, not destroyed data. That is the whole shallow regime, so M3's error model has to use the evidence that survives at one read:

reads/UMI	recall	precision	ε / true
1.11	0.022	0.254	0.20
1.51	0.235	0.567	0.58
2.32	0.587	0.846	0.86
3.12	0.800	0.936	0.98
7.12	0.959	0.977	0.96

Barcode base quality as evidence (BarcodeEvidence::position_error). QX was carried through checkout and unread. A sequencing miscall in the barcode has low Phred at the base it changed; an early-PCR child has a high one in every read. Works at one read.
Payload agreement as evidence (BarcodeEvidence::payload). A barcode error child is a read of the parent's molecule, so its payload matches. Worth log(1/clonality), and the clonality is measured from the data by sampling random barcode pairs — decisive in a diverse repertoire, worth nothing in a clonal library, and the number says which this is. It also lifts the count gates, which is what makes a singleton-vs-singleton merge possible at all, and it refuses merges the count ratio would have made on a disagreeing payload.
The error likelihood is a rate, not a conditional. The zero-truncated Poisson divided out (1 − e^−λ) — precisely the term saying whether a child should exist — so for a singleton child it tended to 1 for every λ and the error rate stopped mattering at exactly the coverage where nothing else was available. Untruncated, both sides are expected counts of neighbouring barcodes and the comparison is like for like.

Measured after (scripts/correction_accuracy.py), against the achievable ceiling rather than against all children — a child whose parent barcode was never sequenced has nothing to merge into:

reads/UMI reachable recall of those precision molecules kept
1.11 0.204 0.108 0.818 1.000
2.32 0.904 0.816 0.830 0.987
3.12 0.975 0.914 0.926 0.991
7.12 1.000 0.979 0.997 0.999
13.30 1.000 0.983 0.999 1.000

Note: At ~1 read/UMI 80% of barcode errors are unfixable in principle — the parent was never sequenced — and of the rest migec fixes 11% while destroying no real molecule. Precision is the side to err on: a wrong merge deletes a molecule and nothing downstream can tell, a missed correction only inflates the count.
migec refine works: barcode table, correction, read rewrite with OX:Z: preserving the original, <sample>.barcodes.tsv, coverage histogram after correction. Recovered 20,055 molecules from 20,000 simulated with ε at 0.96x of injected. Holds the table, never the reads; three streaming passes
Quality calibration measured against the pattern's own constant bases, fitted as ê(q) = ε_qi + a·10^(−q/10) weighted by bases per Q. On SRR1763769 the slope is 1.04 over 46.3 M bases, so the reported Phred's scaling is right. Never: The intercept (3.9e-3) is not a sequencing floor -- the standard is a synthesised oligo, and synthesis runs ~1 defect per 200-500 bases. It is spread evenly over all 23 anchor positions, none polymorphic, and matches the independently measured 0.55% one-base-short rate. Reported as a diagnostic of the primer and left out of error()
Sequencing vs quality-independent separation of the template (2026-08-14). The pattern bases can only calibrate the primer, so the standard is the template itself: at 20+ reads a consensus has suppressed sequencing error to nothing, and what it still gets wrong against the library's modal sequence is what was in the molecule before amplification. --pre-amp-error auto (M1) is that measurement. Still open, and it is data rather than an estimator: splitting the floor into RT and first-cycle PCR needs the same template through both chemistries
Correction posterior: birthday prior with Rényi-2 collision entropy, phred, and a polymerase mixture component for early-cycle PCR children. The distance-1 background comes from X3's column shuffle, not from C(n,2)·P_coll·shell
MIG-size threshold at a target FDR; keep-orphan retention. The residual is measured, not derived: a surviving barcode that still looks like a child of a surviving neighbour, by count or by its reads agreeing on the molecule. Note: Count alone reports zero residual at 1-3 reads per UMI, which is where it is worst. On a 1.23 reads/barcode library: 5.25% of 1-read molecules, threshold ≥2; on 4.62 reads/barcode: 0%, threshold 1. Never: Reported, never applied — every molecule stays in the output
Bucketed correction (2026-08-14). refine's table range-partitions itself past a byte budget and correction follows it, in two passes with the key rotated -- a plain partition splits a barcode from its neighbour for the top b/2 positions and would bound the memory while silently not correcting them. Two things it needed that the counter version did not. The table carries the evidence: a BarcodeEvidence indexed against entries() cannot survive a partition, and dropping it leaves the bucketed run on the count ratio alone, which reports nothing at 1-3 reads/UMI. And the two passes scan rather than merge, with one global apply afterwards: a barcode can have a plausible parent on each side of the boundary, and merging inside a pass takes the first rather than the best -- 2 barcodes in 6,591 landed elsewhere than the resident run put them. With both, a partitioned run and a resident one agree on every scalar and on every output file byte for byte (tests/synthetic/test_refine_bucketed.py). It also folded the evidence pass into the table pass, so refine now streams the reads twice rather than three times
Cell calling (OrdMag + knee) and the QC tables. Molecules per cell, never reads; the knee reported next to the call and a warning when they disagree by more than 3x. On 500 real cells over 20,000 ambient barcodes it calls exactly the 500. <sample>.cells.tsv
QC tables: <sample>.rank.tsv (molecule rank + CDF, log-spaced), <sample>.bins.tsv (per MIG size: barcodes, reads, fraction merged as error, payload entropy), drawn by notebooks/refine_diagnostics.py
<sample>.cell_rank.tsv (2.1.0) — Cell Ranger's barcode rank plot, cells sorted by distinct UMIs with the call on the curve. Never reads: one over-amplified molecule would put an empty droplet high up, which is the artefact the plot exists to show
<sample>.sizes.tsv (2.1.0) — the MIG size spectrum at exact sizes, one row per distinct depth. Power-of-two bins turn the rank/Zipf curve into four steps
Gate: estimated ε within 20% of injected at 1–3 reads/UMI, not only at 7; ≥95% of no-parent 3–5-read MIGs retained (already ≥99% at every depth measured)

reads/UMI	reachable	recall of those	precision	molecules kept
1.11	0.204	0.108	0.818	1.000
2.32	0.904	0.816	0.830	0.987
3.12	0.975	0.914	0.926	0.991
7.12	1.000	0.979	0.997	0.999
13.30	1.000	0.983	0.999	1.000

M4 — end to end

suggest, subsample, marimo notebooks, full docs
Gate met, measured (docs/downstream.rst): minimap2 -ax sr -y and bwa mem -C carry RX/CB/MI into a valid sorted BAM on 600/600 records; arda amplicon reads the consensus directly and its AIRR sequence_id is the molecule id; salmon and kallisto quantify it plainly. STAR unverified — the brew arm64 build reads 0 reads from any FASTQ
plot — twenty QC panels drawn with gnuplot from the tables the stages already write, and assets/ holds the pipeline figure (graphviz) and the example panels the README shows, regenerated by scripts/example_figures.py
The four familiar figures (2.1.0) — barcode rank on Cell Ranger's axes, the MIG size spectrum on log1p with molecules and reads, the rank/Zipf curve, unique UMIs per sample barcode. Consensus quality is a box over an exact (depth, quality) grid, never a thinned scatter: quality is discrete and capped, so a cloud draws a flat line whatever the bin holds and the thinning removes the only thing it could have added
Publication defaults (2.1.0) — transparent background, one ink colour that reads on light and dark, key inside the plot box rather than in a gutter that widens every figure
Overrepresented k-mers in suggest — exact counts in a flat $4^{8}$ array, measured against the reads' own composition, stitched back into the sequence they came from. Run on a stage's OUTPUT it answers "did the trim remove the primer", which nothing else here could
checkout.trimming.tsv — the payload length distribution after trimming. A pattern matched one base off still matches; this is where that shows
Never: sort is not a command. Partitioning happens inside assemble and exposing it would document a third format with no independent meaning (project/design-io-interop.md).

M5 — benchmarks and release

Deliberately not doing

Alignment and variant calling — MAGERI's job; the pipeline hands off to arda/minimap2/bwa-meme.
Indels — Illumina rates ~1e-6/base and no dataset to verify against. Substitutions only.
Duplex consensus (DCS) — v2.0 extracts duplex tags and emits single-strand consensuses. No error-suppression claim in this repo is based on duplex data until DCS exists.
Full-length receptor assembly, doublet calling and contaminating-chain filtering — arda's job. --contig assembles one molecule's fragments into a contig and stops there.
EmptyDrops-style cell rescue — Cell Ranger's job. Ours is OrdMag plus a knee, and the benchmark gate is written against that rather than against a Jaccard we cannot reach.
An external merge sort — range partitioning plus an in-RAM sort per bucket covers it, and nbuckets == 1 is the in-memory case, so there is one code path rather than two.