Skip to main content

migec

PyPI CI docs python C++ license

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 — checkout, refine and assemble — with cell barcodes, whitelists, dual-end and positional (10x) layouts, cell calling, QC figures, and suggest/subsample/plot. Index hopping, .mig bucket output and 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. Collapsing them into a consensus removes essentially all sequencing error — which is what makes rare-variant detection and error-free repertoire profiling possible. The difficulty is entirely in the details:

  • Barcodes acquire errors too. Distinguishing an error-child barcode from a genuine collision needs the birthday bound, the base qualities, and the fact that a polymerase error in an early PCR cycle carries high quality in every read that inherits it. Treating that as a sequencing error is the dominant residual mistake in UMI counting.
  • A molecule seen three times is still information. Cutting at a coverage threshold throws away real sequence. migec keeps low-coverage molecules that have no plausible parent and reports the uncertainty rather than deleting the data.
  • Consensus cannot fix an error made before amplification. An RT or first-cycle PCR error is in every read. Any quality above that floor is a fiction, so migec measures the floor from the data and refuses to claim more.

Pipeline

the migec pipeline

Three stages, plus three tools that read no reads: suggest says where the barcode is, subsample cuts a fixture that is still a library, and plot draws every QC panel with gnuplot from the TSVs the stages already wrote. Regenerate the figure with dot -Tsvg assets/pipeline.dot -o assets/pipeline.svg.

Output is ordinary FASTQ. 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 (for tools that keep them), so this works and was measured (docs/downstream.rst):

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

minibwa spells the comment flag -y on map and -C on the legacy mem subcommand — the wrong one exits with an error rather than dropping the tags quietly. Whether to align before grouping (position + UMI, the fgbio/UMIErrorCorrect order) or after collapsing is a real choice, and docs/downstream.rst works through what each one buys.

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.

Install

pip install migec

Wheels for CPython 3.10–3.13 on Linux x86-64 and macOS arm64. From source: bash setup.sh.

Usage

There is exactly one thing migec has to be told: where the barcode is. Most libraries put it 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
migec checkout R1.fq.gz R2.fq.gz --bc-pattern '^XXXXXXXXXXXXXXXXNNNNNNNNNN' -o out/   # the same

N is a UMI base, X a cell-barcode base, and slices are half-open and 0-based like Python's, so 0:8 is eight bases and the next slice may start at 8. A leading ^ — and every slice list, since a position is only a position if it is measured from somewhere — anchors the barcode at the first base, which is what --max-offset 0 used to have to say by hand. Getting that wrong is not a tuning mistake: a layout with no constant sequence gives a free scan no evidence to choose an offset with, and migec refuses rather than picking one.

Or name the chemistry (migec sheet --presets prints all of them, and where each layout is written down):

preset layout
umi ^NNNNNNNN generic inline UMI
migec cagtggtatcaacgcagagtNNNNtNNNNtNNNN MIGEC 5'-RACE RepSeq
primerid NNNNNNNNNcagtttaacttttgggccatcca HIV-1 Primer ID, as used by MAGERI
duplex ^NNNNNNNNNNNN..... on both mates duplex sequencing
10x ^XXXXXXXXXXXXXXXXNNNNNNNNNNNN 10x Chromium 3' v3
10x-v2 ^XXXXXXXXXXXXXXXXNNNNNNNNNN 10x Chromium 3' v2 and 5'
tso500 ^NNNNN..... on R1 Illumina TSO500 ctDNA — read the warning in docs/layouts.rst
smarter-umi ^NNNNNNNNNN... SMARTer template-switching RNA-seq
migec checkout R1.fq.gz R2.fq.gz --preset 10x-v2 -o out/

fgbio, Picard, samtools and TSO500 write the same thing as a read structure, taken verbatim:

migec checkout R1.fq.gz R2.fq.gz --read-structure 5M5S+T -o out/    # TSO500: `5M5S+T +T`

Many samples in one file

Then it is a barcode table — MIGEC's, read verbatim. Uppercase is matched exactly (IUPAC degeneracy allowed), lowercase is the fuzzy adapter region, and UMI runs need not be contiguous:

S1	aaACTcagtggtatcaacgcagagtNNNNtNNNNtNNNN
S2	aaAGAcagtggtatcaacgcagagtNNNNtNNNNtNNNN

Column 3 is MIGEC's slave pattern — a second pattern on the other mate whose captured positions extend the UMI, which is how a 24 nt dual-end barcode is declared:

S1	NNNNNNNNNNNNtgact	agtcaNNNNNNNNNNNN
migec suggest reads.fq.gz                            # where is the barcode? read it off the data
migec sheet barcodes.txt                             # what will each row extract?
migec checkout reads.fq.gz -b barcodes.txt -o out/
migec checkout R1.fq.gz R2.fq.gz -b barcodes.txt -o out/ -t 8
migec refine out/S1.fq.gz -o ref/                    # correct barcode errors
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
migec plot out/                                      # QC figures from the tables just written

Every stage takes -t/--threads (one per core by default) and --limit-read N / --limit-umi N, which stop the intake after N reads or N distinct barcodes. Limits are for getting an answer out of a 400 GB run in a minute. Never a sample: the first N reads of a FASTQ are one corner of one flowcell, so a limited run reports that it was limited and nothing measured under one transfers to the library. subsample is the sampler.

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

2.2 s (903,599 reads/s) = 1.5 s matching on 8 threads + 0.7 s UMI statistics, serial
peak RSS 131.0 MB of which UMI counters 21.2 MB

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

Paired input searches both mates for the tag and swaps the pair so R1 always carries it — an amplicon library sequenced in both orientations otherwise loses half of each MIG at consensus, and nothing upstream reports it.

Reads come out trimmed of adapter, sample tag and UMI, with the barcode carried in SAM-style tags that survive bwa mem -C and minimap2 -y into the BAM:

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

alongside checkout.summary.tsv, checkout.coverage.tsv (the MIG size histogram) and checkout.umi_composition.tsv (per-position base usage, entropy, information content).

Note: umi_tools spells a cell barcode C, which is cytosine here; pasting one is refused with the translation rather than compiled into a pattern that matches nothing. And on a barcode-only read — 10x R1 is 26 nt of barcode and nothing else — refine and assemble take R2.

On sc5p_v2_hs_PBMC_1k VDJ-T: 100% of 3,155,166 reads assigned, 221,024 barcodes at 14.28 reads each, 813 cells called.

It corrects the barcodes, with the evidence that survives at one read

migec refine out/S1.fq.gz -o ref/

A barcode one substitution from another is either an error child of it or an independent molecule. The count ratio separates them on a deep amplicon and is worth nothing at 1–3 reads per UMI, so refine also uses the barcode's own base quality at the position that differs — checkout already writes it to QX — and payload agreement, since an error child is a read of the parent's molecule. Agreement is worth log(1/clonality), and the clonality is measured rather than assumed:

barcodes    23,910 distinct
  merged    3,855 (16.1%) into a parent, 3,889 reads moved
molecules   20,055 after correction          <- 20,000 were simulated

barcode error   2.87e-03 per base            <- 3.0e-03 injected
clonality       0.0100 of random barcode pairs carry the same payload anyway
                -- payload agreement is worth about 100x odds towards the same molecule here

Note: At ~1 read per UMI 80% of barcode errors cannot be fixed by anyone — the parent barcode was never sequenced. migec corrects a tenth of the rest and destroys no real molecule at any depth measured, which is the side to err on: a wrong merge deletes a molecule and nothing downstream can tell, while a missed correction only inflates a count. Every corrected read keeps what it was in an OX:Z: tag. See docs/refine.rst.

It collapses each molecule, and caps what it claims

migec assemble out/S1.fq.gz -o cons/
migec assemble out/S1.fq.gz -o cons/ --contig     # random-primed reads that tile the molecule

A molecule is sample + cell barcode + UMI, never the UMI alone — the same UMI in two cells is two molecules, and that is the design rather than a defect. Reads are range partitioned on the packed key into .mig buckets and one bucket is sorted at a time, so nothing scales with the library: 531,365 reads/s, and 121 MB at 16 buckets against 203 MB at one.

The per-column posterior is LL[j][b] = Σ_i (r==b ? log(1−e) : log(e/3)), and then the number that matters:

Q(j) = −10 log10( p_cons(j) + p_floor )

The floor is added, not compared, and it is named rather than guessed — an RT or first-cycle-PCR error is in every read of the molecule and no consensus removes it:

migec assemble ... --rt-error rt        # 1e-4, caps at Q40. Anything with an RT step (default)
migec assemble ... --rt-error medium    # 1e-5, caps at Q50. No RT, an ordinary polymerase
migec assemble ... --rt-error high      # 1e-6, caps at Q60. No RT, a proofreading polymerase
migec assemble ... --rt-error 7.37e-5   # or the rate itself, e.g. TruSight Oncology 500 v2

This is the one-molecule floor and every record here is one molecule. 10x state it exactly: "The estimated error rate for the V(D)J RT reaction is 1e-4 per base. Therefore, assembled bases that are covered by a single UMI are assigned Q40, and bases covered by at least two UMIs are assigned Q60." The Q60 branch needs two molecules to agree — an RT error is common-mode within a molecule and independent between them — and combining molecules is arda's job. 1e-4 is also what X2 measured here independently (1.54e-4 on SRR1763769), and the polymerase classes come from McInerney et al. 2014 (Taq 4.3e-5, Pfu 2.8e-6, Phusion 2.6e-6, Pwo 2.4e-6 per bp per duplication) — with the first cycle worth ~5x an ordinary one (Shagin et al. 2017).

Which is exactly what comes out: one read carries the payload's own error and the curve flattens at the floor rather than at the instrument — drawn as a box per depth bin, further down.

Very deep barcodes are capped at 10,000 reads into the consensus — past that the column posterior has long since saturated while the group still costs time and memory, which is 10x's rule and their reasoning. The cap applies to the reads that are consensed, never to the reads that are counted: cD stays the true depth of the molecule.

--contig is for random-primed libraries, where reads sharing a barcode tile the molecule instead of starting at the same base. They are placed against each other by seed matching, cut into overlap components, and one consensus is emitted per component — a component is never extended across a gap, because 27.3% of 10x groups hold more than one and a single consensus over those asserts sequence no read covers. Assembling a cell's full receptor and calling doublets is arda's job, not this one.

Note: Contig assembly needs a barcode that is not saturated: two fragments of two different molecules sharing a barcode have no sequence in common, which is exactly what two fragments of one look like. assemble runs the same birthday arithmetic on the barcodes it saw and reports how many molecules a group holds on average — above 1, the warning says so. See docs/assemble.rst.

When the deliverable is a count, not a sequence

migec assemble ref/S1.fq.gz -o cons/ --fast

Counting mode: the group's most frequent exact sequence, with every base carrying the best quality any read of that sequence reported. No column model, so no per-base error correction and no sub-clustering — and the RT floor still caps what it claims. Use it for expression and clonotype abundance, where a molecule count is the answer and the sequence only has to be right often enough to assign. Measured against the full path on 8-read molecules at 5e-3 per base: the column posterior removes essentially every sequencing error, and the majority string keeps whatever it carried. --fast is refused with --contig, whose tiling reads share no exact sequence to take a majority over.

It draws its own QC

migec plot out/                       # every panel whose table is in out/
migec plot cons/ -o figs/ --format pdf

Twenty panels. Every one is a gnuplot script over a TSV a stage already wrote, so a figure can be redrawn from the table next to it long after the FASTQ is gone, and a figure can never disagree with the number in the report. gnuplot is not a Python dependency: without it the .gp scripts are still written. Every SVG is transparent and mid-grey, so one file serves a light README, a dark README and a printed page — and the legend sits inside the plot box, not in a gutter that makes every figure wider than its data. See docs/plots.rst.

Four of them are the ones you already know how to read.

The barcode rank plot, on Cell Ranger's axes, because that is the figure every user of a droplet protocol has seen. Barcodes sorted by how many distinct UMIs they carry — never by reads, since one over-amplified molecule would otherwise put an empty droplet high on the curve, which is the artefact the plot exists to show. The call is drawn on the curve.

barcode rank plot

The MIG size spectrum, molecules and the reads they account for, on log(1 + size). Both series, on their own axes, because they peak in different places the moment a library is over-sequenced: most molecules are shallow, most reads are in the deep ones. A figure with only the first says the library is fine; a figure with only the second says it is saturated. log1p, so a molecule seen once still has a place on the axis.

MIG size spectrum

The rank/Zipf curve — molecule size against rank, log-log. A straight line is Zipf, and amplification bias bends it. This is why refine writes the size spectrum at exact sizes rather than in power-of-two bins: four bins make four steps, and a straight line cannot be told from a bent one.

molecule size against rank

Consensus quality against depth, as a box, never a thinned scatter. Emitted quality is discrete and capped at the floor, so at any real depth every molecule sits on one or two integers: a cloud of dots draws that as a flat line whether the bin holds ten molecules or ten million, and thinning it throws away the tails that were the only thing the cloud could have shown. These are exact order statistics over every molecule, read off a (depth, quality) count grid.

consensus quality against depth

checkout.summary.tsv carries unique UMIs and reads per sample barcode, drawn the same way by the sample_umis panel — the multiplexed analogue of the same question.

It finds what the trim left behind

migec suggest profiles any FASTQ, so point it at the output of a stage rather than its input:

migec suggest out/S1.fq.gz -o qc/       # did the trim actually remove the primer?

An 8-mer occurs by chance about every 65 kb, so synthetic sequence that survived shows up as a run of k-mers each shifted one base from the last, hundreds of times more often than the reads' own base composition predicts — and the run is stitched back into the sequence it came from:

kmer           count   obs/exp    reads  mean pos
GGGCCATC      20,023     701.2  100.0%      22.1
TGGGCCAT      20,018     656.4  100.0%      21.0
TTGGGCCA      20,018     656.4  100.0%      20.0

overlapping into: CAGTTTAACTTTTGGGCCATCCA

Overrepresentation is against the reads' own base composition, never a flat 1/4: a 70% AT library makes every AT-rich k-mer look enriched against uniform, and the table would then be a description of the GC content rather than a finding.

It tells you where the barcode is

migec suggest reads the layout off the reads rather than off the protocol. A UMI cycle is one the synthesiser mixed — all four bases near 1/4, ~2 bits. A constant cycle is one base near 100%. Everything else is payload.

 cycle      A      C      G      T  1/4 dev     Q  layout
     0  0.271  0.205  0.257  0.267    0.045    33  N  UMI
     ...
     9  0.020  0.971  0.004  0.006    0.721    37  |  constant

segments:
    0-8   umi         9 nt  (mean 1/4 deviation 0.038)
    9-31  constant   23 nt  (mean 1/4 deviation 0.718)  CAGTTTAACTTTTGGGCCATCCA

pattern  NNNNNNNNNcagtttaacttttgggccatcca

That is a real HIV Primer ID library (SRR1763769) with nothing supplied but the FASTQ. The pattern pastes straight into a barcode table, and checking it out assigns 95.0% of reads.

It tells you whether the barcode was big enough

A 12 nt UMI is 4¹² = 16,777,216 sequences — if the synthesiser delivered exactly 25% of each base. It never does, so the usable space is the collision (Rényi-2) entropy, 1 / Π_j Σ_a p_j(a)², never Shannon: H₂ ≤ H₁, so Shannon overstates the space and understates collisions, which is the direction that silently merges molecules.

From there the birthday problem, in the form that survives a full space:

occupied = S·(1 − e^−λ)     molecules = S·λ     P(k>1 | k≥1) = (1 − e^−λ − λe^−λ)/(1 − e^−λ)
sample              space  occupancy  MIGs >1 mol   molecules   err pred   err est
CTRL              250,902      49.9%        30.6%     173,482    2.0e-03   2.7e-04

warning: CTRL: 31% of MIGs hold more than one molecule (50% of a 250,902 barcode space is
  occupied). Their consensus is a mixture of templates, not a molecule
warning: CTRL: the barcode error estimate (2.7e-04) is not reliable here -- 50% of each barcode's
  1-substitution neighbourhood is itself occupied ...

scripts/collision_check.py checks that prediction against something model-free — two molecules sharing a barcode with different sequences are visible in the reads — and finds 1.86× more collisions than predicted. That is not the position-independence assumption in Π_j m_j, which a permutation puts at 1.01× — it is the read threshold, since a collided barcode carries two molecules' reads and is over-represented among the MIGs big enough to show a split.

The barcode error rate is estimated from the distance-1 excess and reported next to what the reported Phred and the polymerase predict. Note: The estimator has a working range: it recovers 0.92× of an injected rate at 0.3% occupancy and 0.23× at 50%, always collapsing downward, so it is flagged unreliable past 5% neighbourhood occupancy rather than quietly believed.

Full derivations in docs/barcode_space.rst; notebooks/barcode_space.py draws it.

Every derivation has a permutation that checks it

scripts/permutation_nulls.py measures three quantities the pipeline otherwise derives, assuming nothing. On the same HIV library — 125,369 distinct 9 nt barcodes, 47.8% occupancy:

derived from a model measured by permutation
positions are independent 1.01× excess, purely nearest-neighbour holds, to ~1%
distance-1 pairs are error children 97% are chance; ~18,000 are real permute the background
split a MIG at nominal p < 0.01 1% false positives at -log10 p = 8.68 19× over-call

The independence null is a distribution — the product measure q(u) = Π_j p_j(u_j) — so it is tested with Jensen-Shannon divergence against a same-size draw from q, not with one functional of each. The residual dependence is entirely between adjacent positions, and it has a cause: 0.55% of reads carry a barcode one base short, a coupling step that did not fire, which frameshifts everything after it.

The last one is the one that mattered. Reads are not exchangeable — a low-quality read carries a minor base at many positions at once and looks exactly like a linked subclone — so the null has to preserve both margins of the reads × positions matrix, per-position error count and per-read error load. The nominal threshold calls 30.62% of MIGs as two molecules; the permutation calls 1.60%. The threshold is a Monte Carlo estimate and its error is quoted: 8.68, bootstrap 95% CI [8.42, 9.14] over 82,800 randomisations — a tail quantile from a tenth as many gave 9.61 and 11.66, so the interval is the number, not the point. See docs/nulls.rst.

Most UMIs have 1–3 reads, and that is the normal case

Bulk repertoire profiling and shallow 3' single-cell both put the MIG size histogram's mass at 1–3 reads. migec runs there and says what it can support rather than quoting a number calibrated on a deep library: the split threshold is inert (a pair of columns can carry at most log10 C(n, n/2), so it needs ~30 reads), the count-ratio error-child null has no dynamic range, and nothing is thresholded away--min-reads defaults to 1, because a molecule seen once is still a molecule and the answer to a barcode error is to correct it, not to delete it.

    MIG size      groups    share
           1      31,888    79.4%
         2-3       8,176    20.4%
         4-7         112     0.3%

warning: 79.4% of molecules were seen once. A consensus over one read is that read --
  the UMI is buying counting here, not error correction

It is also the memory-hostile case, because distinct barcodes are what everything scales with, so it is what the benchmarks use: 190,595 reads/s at 1.02 reads/UMI, 259 B resident per distinct barcode, still bounded by the bucket rather than the library.

Speed and memory are reported, not assumed

--threads defaults to one per core on all three stages, and the output is byte-identical whatever it is set to — checkout matches in chunks and writes them back in input order, refine scans the barcode neighbourhood in parallel and applies the merges it finds serially, and assemble gives each worker its own bucket and concatenates them in bucket order. -t changes the wall clock and nothing else, which is asserted per stage in C++, at the CLI, and under the thread sanitizer.

stage 1 thread 16 threads bound by
checkout 202,717 1,056,472 reads
refine 605,611 1,012,368 distinct barcodes
assemble 549,745 2,051,937 reads, then the largest bucket

reads/s. Two of those used to be 222,017 and 202,977, and the first fix was not a thread: zlib at its default level 6 was 83% of refine's wall clock, compressing an intermediate the next stage decompresses immediately. Level 1 costs 21% more bytes and gave 3x before a single thread was added. Measure the stage before parallelising it.

The same lesson twice in assemble: after the consensus was threaded, the partition was 2.07 s of a 2.69 s run against a 0.23 s gzip -dc floor for the same file — so it was not the inflate. It runs on the workers now, by ownership rather than locking (worker w owns every bucket with bucket % threads == w), and half the win was in the reader: the chunk is assigned into rather than cleared, because clear() destroys four std::string per record and the reader ends up in malloc instead of inflate.

4 M reads, -t 16 before after
wall clock 2.70 s 1.95 s
reads/s 1,481,946 2,051,937
peak RSS 1,479 MB 789 MB

The memory fell with it, because the estimate deciding how finely to cut the input said a gzipped FASTQ goes resident at 8x its on-disk size and it is really 19x — a resident record is two heap std::string with their allocator headers, not the 180 bytes of payload. Guessing low picks too few buckets, and pass 2 holds sixteen of them at once.

threads reads/s matching reads/s peak RSS
1 202,717 217,506 60 MB
2 362,486 412,855 78 MB
4 603,773 758,801 97 MB
8 910,500 1,335,791 136 MB
16 1,056,472 1,684,654 217 MB

checkout thread scaling

2 M single-end 129 nt reads, four barcode patterns, 4 reads per molecule, M-series laptop. python scripts/benchmark_threads.py --reads 2000000 -o assets/ writes that table, and the figure is drawn from it — the two cannot drift apart. Two things had to be true for the matching to scale. zlib compresses random DNA at 7 MB/s at its default level 6, so compression runs on the workers (concatenated gzip members are a valid gzip stream) at level 1 — 137 MB/s for 13% more bytes. And the log-likelihood score tabulates into 1.2 kB, because the log2 in the inner loop was 90% of runtime.

Both columns are reported, because they scale with different things and only one of them threads: matching scales with reads, the UMI statistics with distinct UMIs at ~1.5–2 µs each, on one thread. That serial tail is why the end-to-end column flattens at 16 while matching is still climbing.

The UMI counters are a sorted (key, count) array rather than a hash map: ~22 bytes per distinct UMI against ~48, which at the 4·10⁸ distinct UMIs of an ordinary NovaSeq run is 8.8 GB against 19 GB. That still does not fit a laptop — the counters are not yet partitioned, so checkout warns when they pass 1 GB rather than letting you find out from the OOM killer. See docs/performance.rst.

Grouping accuracy is measured against Calib

scripts/compare_calib.py scores both tools' read partitions against a known truth with the adjusted Rand index, and reports splitting and merging separately — splitting inflates the molecule count and is recoverable, merging mixes molecules and destroys real variants.

UMI UMI error ARI reads split reads merged
12 nt 0 1.0000 0.0000 0.0000
12 nt 5·10⁻³ 0.9348 0.5165 0.0004
6 nt 0 0.8877 0.0000 0.3982

Calib clusters on barcode and sequence; migec today groups on the barcode alone, and the gap is exactly the collision rate — which eff len predicts before any clustering runs. A clean 12 nt barcode needs nothing cleverer; a 6 nt one cannot be rescued by any amount of barcode cleverness, only by sequence, which is what assemble adds.

eff len is the number to look at

A 12 nt UMI is not worth 12 nt unless its bases are uniform. effective_length is $-\sum_j \log_4 \sum_a p_j(a)^2$ — what the barcode is actually worth. A 12 nt UMI with eight fixed positions has an effective length of 4, a usable space of 256, and will collide constantly.

The distinction matters more than it looks: a sequence logo draws Shannon entropy, but the probability two molecules collide is the Rényi-2 (collision) entropy. Since H₂ ≤ H₁, using Shannon overstates the usable space and understates collisions — the direction that silently merges distinct molecules. Both are reported; only the collision form feeds any decision.

Pipelines

integrations/nextflow/ is an nf-core-style local module set — four processes, a subworkflow, meta.yml, nextflow.config, environment.yml — that drops into nf-core/airrflow or anything else that hands you FASTQ pairs. SLURM is the pipeline's business, not the module's: it declares label and task.cpus and nothing more.

All three stages thread, and each is byte-identical at any -t, so a retry with different cores cannot change a result — which is what makes an escalating errorStrategy 'retry' safe here. Per-sample keys in meta beat the params.migec_* defaults, so one run can mix chemistries. See docs/nextflow.rst; nextflow is not installed on the machine these docs were measured on, so the modules are reviewed against the nf-core spec rather than verified by a run.

Documentation

https://antigenomics.github.io/migec/ — see docs/formats.rst for the on-disk format, and ROADMAP.md for what is implemented.

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.

Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Source Distribution

migec-2.1.0.tar.gz (616.4 kB view details)

Uploaded Source

Built Distributions

If you're not sure about the file name format, learn more about wheel file names.

migec-2.1.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (561.1 kB view details)

Uploaded CPython 3.13manylinux: glibc 2.17+ x86-64

migec-2.1.0-cp313-cp313-macosx_11_0_arm64.whl (378.7 kB view details)

Uploaded CPython 3.13macOS 11.0+ ARM64

migec-2.1.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (561.2 kB view details)

Uploaded CPython 3.12manylinux: glibc 2.17+ x86-64

migec-2.1.0-cp312-cp312-macosx_11_0_arm64.whl (378.2 kB view details)

Uploaded CPython 3.12macOS 11.0+ ARM64

migec-2.1.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (560.1 kB view details)

Uploaded CPython 3.11manylinux: glibc 2.17+ x86-64

migec-2.1.0-cp311-cp311-macosx_11_0_arm64.whl (376.3 kB view details)

Uploaded CPython 3.11macOS 11.0+ ARM64

migec-2.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (559.6 kB view details)

Uploaded CPython 3.10manylinux: glibc 2.17+ x86-64

migec-2.1.0-cp310-cp310-macosx_11_0_arm64.whl (374.6 kB view details)

Uploaded CPython 3.10macOS 11.0+ ARM64

File details

Details for the file migec-2.1.0.tar.gz.

File metadata

  • Download URL: migec-2.1.0.tar.gz
  • Upload date:
  • Size: 616.4 kB
  • Tags: Source
  • Uploaded using Trusted Publishing? Yes
  • Uploaded via: twine/7.0.0 CPython/3.13.14

File hashes

Hashes for migec-2.1.0.tar.gz
Algorithm Hash digest
SHA256 f4b3fc57c7363648ff95b1d6d8a6dc4484006a347491e1ed36d8b67e6aaa28b7
MD5 3a351982f58daec115f8eb63fee6f51a
BLAKE2b-256 c914c8a03efc75eec1ac0c0a33f15926a36544b52da31c247eaeafbc557945cb

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0.tar.gz:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm Hash digest
SHA256 aaef5a2832f2948a0030bf022d6cf5515d93396e7fe61c5b7441225fa5d3ae33
MD5 1ab109517678c6e44ad5f35d3fd0cee7
BLAKE2b-256 9ac117269d8b012473ea0a095d8350b9a2a23973e46445bea33f96fa834e0231

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp313-cp313-macosx_11_0_arm64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp313-cp313-macosx_11_0_arm64.whl
Algorithm Hash digest
SHA256 50dabeb27891d7c0b380dc0c1a1ce9c6cfa92013961be60cfbe1df05ca07261b
MD5 4d9e0ef57442b35527cb3668943feed4
BLAKE2b-256 c692ebaeae3b313c4710381defe6d5e0c4a6c760cc3d310709db06923e6be88f

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp313-cp313-macosx_11_0_arm64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm Hash digest
SHA256 78e2ebd03b0b58cce4fca5f4c735c288c3f17ea3989bbd5f04f86e5a2f71ff8a
MD5 f01217cc2cd967de5675af2e45a321c4
BLAKE2b-256 3f5aee47894b19e7efcb68296da5eb287a8bdf49023e1dada67445482cce3e87

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp312-cp312-macosx_11_0_arm64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp312-cp312-macosx_11_0_arm64.whl
Algorithm Hash digest
SHA256 ae94ab14c80f26573e27b8c03330a2acdf7f3e643dfa3bdb01c89dfee9b38cc0
MD5 de725170122da847e8444785dd931260
BLAKE2b-256 5e2d034214c3e510aad485107bdfbad8b69425d75f40c0538653ab04b7473419

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp312-cp312-macosx_11_0_arm64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm Hash digest
SHA256 1df009ca8a9895868f6624faddfcd2b491351175a3d66be6e62798c9ccf65529
MD5 7be1def3fc455949fb2bce723cb60b2d
BLAKE2b-256 fc09bd99d06802828a1768197cac13ffaffe18a7227ce869dc8f847633f4cd3b

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp311-cp311-macosx_11_0_arm64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp311-cp311-macosx_11_0_arm64.whl
Algorithm Hash digest
SHA256 0c5831a4ead9b6be85cbaf22e2885f418fece3551b2f8f1b9b23a30b0a2d31d5
MD5 304ac01a677821d2e06018a4ca329472
BLAKE2b-256 e645115f4b4740aa823d53206c2a69aaa47f509451d9c933281b0c716195b819

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp311-cp311-macosx_11_0_arm64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm Hash digest
SHA256 306a8ee82469a941be0ef828bbf92151d76f409c71961ce011973dea9e7316d2
MD5 b3f8eb66d7bc62bcf711fc2817c0fb6e
BLAKE2b-256 e7d522761b41f1328a0eceac59cfe7e27565cd4db8815869cd38c11a1af46972

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

File details

Details for the file migec-2.1.0-cp310-cp310-macosx_11_0_arm64.whl.

File metadata

File hashes

Hashes for migec-2.1.0-cp310-cp310-macosx_11_0_arm64.whl
Algorithm Hash digest
SHA256 2cfc32151c8a8ff42ac88672ce643510359e030873344006803863f54d3fba74
MD5 2d1d2f0cd24c2cec445213731511b2d7
BLAKE2b-256 08bd93656cc4f785d644b33a1d16c17b9ed648272021fb780665a3bc7d0edc66

See more details on using hashes here.

Provenance

The following attestation bundles were made for migec-2.1.0-cp310-cp310-macosx_11_0_arm64.whl:

Publisher: publish.yml on antigenomics/migec

Attestations: Values shown here reflect the state when the release was signed and may no longer be current.

Supported by

AWS Cloud computing and Security Sponsor Datadog Monitoring Depot Continuous Integration Fastly CDN Google Download Analytics Sentry Error logging StatusPage Status page