migec

UMI barcode extraction, correction and consensus assembly for barcoded sequencing data.

A complete C++20 rewrite of MIGEC (Shugay et al., Nature Methods 2014) and MAGERI (Shugay et al., PLoS Computational Biology 2017).

Version 2 is under construction. All three stages work today, with cell barcodes, whitelists, dual-end and positional (10x) layouts, mate merging, cell calling, index-hopping estimation, QC figures, and suggest/subsample/plot. The published benchmark comparisons are what remain; see ROADMAP.md. The Groovy MIGEC 1.2.9 is archived on branch legacy-v1 and at tag v1-final — Java users want the jars on the 1.2.9 release.

Why

Tag each molecule with a random barcode before amplification and every read carrying that barcode descends from one original molecule, so collapsing them into a consensus removes essentially all sequencing error. The difficulty is entirely in the details: barcodes acquire errors of their own and an error child has to be told from a genuine collision; a molecule seen three times is still information rather than something to threshold away; and no consensus can repair an error made before the first amplification cycle, because it is in every read. migec measures that floor from the data and refuses to claim a quality above it.

Install

uv tool install migec      # or: uv pip install migec, or: pip install migec
migec info                 # prints the three version strings; they must agree

Wheels for CPython 3.10–3.13 on Linux x86-64 and macOS arm64, so nothing compiles. On a cluster whose system Python is older than 3.10, bring your own:

uv venv --python 3.12 ~/envs/migec && source ~/envs/migec/bin/activate && uv pip install migec

From source, for development: git clone https://github.com/antigenomics/migec && cd migec && bash setup.sh. See installation.

Where the barcode is

This is the one thing migec has to be told. Most libraries put the barcode at a fixed offset in one read, so that is the primary way to say it — a position, no sheet, no anchor:

migec checkout reads.fq.gz --bc-pattern '^NNNNNNNN'  -o out/    # 8 nt UMI at the read start
migec checkout reads.fq.gz --bc-pattern '0:8'        -o out/    # the same, as a half-open slice
migec checkout reads.fq.gz --bc-pattern '0:4,5:10'   -o out/    # 9 nt UMI split by one spacer base
migec checkout R1.fq.gz R2.fq.gz --bc-pattern 'cell:0:16,16:26' -o out/     # 10x

N is a UMI base, X a cell-barcode base, and slices are half-open and 0-based like Python's. A leading ^, and every slice list, anchors the barcode at the first base, so --max-offset never has to be passed. Or name the chemistry — migec sheet --presets prints all of them with the source each layout is written down in:

migec checkout R1.fq.gz R2.fq.gz --preset 10x-v2 -o out/          # a named chemistry
migec checkout R1.fq.gz R2.fq.gz --read-structure 5M5S+T -o out/  # fgbio, Picard, samtools, TSO500
migec checkout reads.fq.gz -b barcodes.txt -o out/                # many samples: MIGEC's own table
migec suggest  reads.fq.gz                                        # if you do not know: read it off the data

A preset says where the barcode is. It does not say what a consensus is worth, and that matters more: the same 12 nt UMI serves a repertoire census and an MRD assay, and the right settings are opposite. migec sheet --assay ctdna prints the second half — the --min-reads and the pre-amplification floor the experiment implies, for eight profiles from airr to mrd (assays, layouts).

The pipeline

migec checkout  reads.fq.gz -b barcodes.txt -o out/   # find the barcode, cut it out, demultiplex
migec refine    out/S1.fq.gz -o ref/                  # correct the errors IN the barcode
migec assemble  ref/S1.fq.gz -o cons/                 # one consensus per molecule
migec subsample out/S1.fq.gz -o small.fq.gz --keep 1  # a fixture that is still a library

Those four and suggest are the pipeline; plot, sheet and info read no reads at all (commands).

Every stage takes -t/--threads (one per core by default) and --limit-read N / --limit-umi N, which stop the intake early — for getting an answer out of a 400 GB run in a minute, never as a sample. -t changes the wall clock and nothing else: the output is byte-identical at any thread count, which is what makes a retry on different cores comparable. Every run says what it did, what it cost, and what the barcode was worth:

reads       2,000,000
  assigned  2,000,000 (100.0%)
  unmatched 0 (0.0%)
  ambiguous 0 (0.0%)

1.6 s (1,243,801 reads/s) = 1.5 s matching on 8 threads + 0.1 s UMI statistics
peak RSS 136.0 MB of which UMI counters 11.5 MB

sample             reads        UMIs  reads/UMI  UMI len  eff len
S1               500,000     125,000       4.00       12    12.00

What comes out

Ordinary FASTQ, trimmed of adapter, sample tag and UMI. One record is one molecule, and its identity is carried twice — in the read name (<sample>.<cell>.<umi>, for tools that drop FASTQ comments) and in tab-separated SAM tags that survive into a BAM:

@r0 RX:Z:GCTAAAGACAAT	QX:Z:IIIIIIIIIIII	BC:Z:S1
TACATAACATACACGTCAGCACGAAACTTGTTGGCCCAGTGTGAATCGCTT

migec plot cons/          # twenty QC panels with gnuplot, straight off those TSVs

How deeply was each molecule sequenced, and is the barcode any good? Two panels answer the questions you ask first, both drawn straight off checkout's own tables on a real 10x library (sc5p_v2_hs_PBMC_1k VDJ-T, 3.16 M read pairs, 311,421 molecules):

Left, MIG size in doubling bins — bars are molecules, the line is the reads inside them. Most molecules are seen once or twice while most reads pile onto a few very deep ones, which is the distribution that decides whether a consensus is worth assembling at all.

Right, the barcode's own PWM, every position of cell barcode and UMI against the 1/4 a free synthesiser mix would give. The two halves read differently on purpose: the 16 cell positions (marked on the axis) swing between 0.198 and 0.305 because they are drawn from 10x's whitelist, while the 10 UMI positions sit tight around 1/4. A UMI position that wandered like a cell position would be a synthesis defect, and effective_length in checkout.summary.tsv is what it costs you — 25.85 usable bases of 26 here.

Each of these was run against real output (downstream):

minimap2 -ax sr -y ref.fa cons/S1.consensus.fq.gz | samtools sort -o S1.bam   # RX, CB, MI in the BAM
minibwa map -y -t8 ref.fa cons/S1.consensus.fq.gz | samtools sort -o S1.bam   # `-y`, not bwa's `-C`
bwa mem -C     ref.fa cons/S1.consensus.fq.gz     | samtools sort -o S1.bam
arda amplicon --r1 cons/S1.consensus.fq.gz -p S1      # AIRR sequence_id IS the molecule id
salmon quant -i tx.idx -l A -r cons/S1.consensus.fq.gz -o quant/   # NumReads are molecule counts

Never run alevin, bustools or STARsolo on a consensus FASTQ. They read the barcode out of a raw barcode read and deduplicate themselves; migec already did, and that read no longer exists.

What it is worth: rare variants, measured

Commercial cfDNA reference material with certified allele frequencies, including a 0%-certified arm that is a true negative by construction. Three replicates per arm, one panel, one aligner, matched molecule-support thresholds, substitutions only — migec emits no indels by design and 56% of UMIErrorCorrect's calls are deletions (post-processing, assets/ctdna_callers.tsv).

migec + LoFreq matches the best sensitivity at every arm and reports 3.8× fewer false positives on the true negative; migec + Mutect2 reports the fewest of anything measured and pays for it at 0.25%. Against MAGERI, the other descendant of MIGEC 1, on a simulated shallow library where both tools emit the same consensuses at the same accuracy: MAGERI reports 142 variants of which 137 are at positions nothing was injected at, migec + LoFreq reports 5 and is right about all five (validation).

The no consensus row is the same reads, trimming, barcode correction, aligner and caller — only a record is a read rather than a molecule. Collapsing cuts false positives 2.8×, makes the measured frequency right (1.02× of certified against 1.27×), and detects more from 38× less depth: at 0.125% the consensus finds the hotspot in 1 of 3 replicates and a 197,772× read pileup finds it in none. A read count is not a molecule count.

Two flags decide more than the choice of caller, and both are measured rather than argued:

gatk Mutect2 --max-reads-per-alignment-start 0 ...   # or it sees 1.5% of your molecules
migec assemble ... --min-reads 3                     # 0% arm: 10.0 -> 2.0 calls per sample

What makes it different

• Barcode correction uses the evidence that survives at one read per UMI — the barcode's own base quality and payload agreement, not only the count ratio, which reports zero there (refine).
• Emitted quality is capped at the measured RT/first-cycle floor, Q40 by default and fitted from the data with --pre-amp-error auto, never taken from the instrument's (quality floor).
• Every model-derived number has something model-free beside it: collisions, the barcode error rate, the split threshold (nulls, barcode space).
• Nothing scales with the library — a range partition into buckets and a sorted counter array, 22 B per distinct UMI against a hash map's 48 (performance).
• Twenty QC panels, each a gnuplot script over a TSV a stage already wrote, so a figure can never disagree with the report (plots).
• For rare variants the molecule count decides, not the caller, and it is fixed before any software runs (variants, detection).
• BAM, SAM and CRAM are inputs, so a capture, exome or ctDNA kit that puts the UMI in the index read never needs checkout at all (bring your own UMI).

Documentation

https://antigenomics.github.io/migec/

Citing

Until the v2 paper exists, cite the original methods:

• Shugay M et al. Towards error-free profiling of immune repertoires. Nat Methods 11:653–655 (2014). doi:10.1038/nmeth.2960
• Shugay M et al. MAGERI: Computational pipeline for molecular-barcoded targeted resequencing. PLoS Comput Biol 13(5):e1005480 (2017). doi:10.1371/journal.pcbi.1005480

License

GPL-3.0-or-later. The archived v1 code on legacy-v1 remains under its original MiLaboratory non-commercial license.