Examples#

Three commands, whatever the platform. Only checkout needs to know where the barcode is.

migec checkout  reads.fq.gz -b barcodes.txt -o co/   # find and cut out the barcode
migec refine    co/S1.fq.gz -o ref/                  # fix errors IN the barcode
migec assemble  ref/S1.fq.gz -o asm/                 # collapse each molecule

Declaring the layout#

Four equivalent ways to say where the barcode is, in the order to reach for them:

form	looks like	when
a position	`^NNNNNNNN` or `0:8`	the primary mode – the chemistry fixes the barcode at an offset
`--preset`	`10x-v2`, `tso500`, `duplex`, …	a chemistry with a name; `migec sheet --presets` lists them
`--read-structure`	`5M5S+T`	fgbio, Picard, samtools and TSO500 all speak this
barcode table	`S1<TAB>aaACTcagtgg...NNNNtNNNNtNNNN`	many samples in one file; MIGEC’s own format, read verbatim

In a pattern, N is a UMI base, X a cell-barcode base, uppercase is matched exactly, lowercase is the fuzzy adapter region, and . is skipped. Slices are half-open and 0-based like Python’s, each a UMI slice unless prefixed cell:. A leading ^, a slice list and a read structure all anchor the barcode at the first base, so --max-offset never has to be passed.

Layouts: where the barcode is has all of it in one place. Do not guess the layout – suggest – read the layout off the data reads it off the data, and Downstream: what consumes the consensus is what happens after.

By platform#

platform	how to declare it	notes
Bulk amplicon (MIGEC, RepSeq)	`aaACTcagtggtatcaacgcagagtNNNNtNNNNtNNNN`	the primer anchors it, so a free scan places it
HIV-1 Primer ID	`NNNNNNNNNcagtttaacttttgggccatcca`	recovered from the data by `migec suggest`
10x droplet	`--preset 10x-v2`, or `--bc-pattern 'cell:0:16,16:26'`	positional; `refine`/`assemble` take R2
TSO500	`--preset tso500`, or `--read-structure 5M5S+T`	`5M5S+T +T`: a 5 nt UMI on R1 only. Warning: 1,024 barcodes does not identify a molecule – TSO500 groups position-aware, downstream
UMI RNA-seq (SMARTer)	`--preset smarter-umi`	10 nt inline UMI, then the `GGG` the template switch leaves
Dual-end (MAGERI)	column 3 of the sheet, `NNNNNNNNNNNNtgact` / `agtcaNNNNNNNNNNNN`	both halves must match
Shallow bulk (1-3 reads/UMI)	any of the above	nothing changes; what the numbers can mean does
Capture, exome, ctDNA, MRD	nothing – the UMI is already in `RX`	the kit puts it in the index read, so there is no layout and no `checkout`. See Bring your own UMI

Starting from a BAM#

A hybrid-capture kit reads the UMI on the index, so what arrives is an fgbio, Picard or vendor BAM with the UMI in RX — never a FASTQ with a barcode inside the read. Two commands, no layout, no checkout:

migec refine   tagged.bam  -o ref/   # correct errors in the UMI
migec assemble ref/S1.fq.gz -o asm/  # one consensus per molecule

BAM, SAM and CRAM are recognised from the file, not the name; samtools does the conversion and the temporary FASTQ is deleted when the stage returns. To see the round trip on data you already have:

migec checkout reads.fq.gz --bc-pattern '^NNNNNNNNNNNN' -o co/
samtools import -T '*' -s co/S1.fq.gz -o S1.bam     # the tags become real BAM tags
migec refine   S1.bam -o from_bam/                  # identical to refining co/S1.fq.gz

Runnable notebooks#

Six marimo notebooks, each a plain Python file with its own PEP 723 dependency header, so uv builds the environment and nothing has to be installed first. Two of them need the network: platforms.py downloads its fixtures from isalgo/umi_data on first run, and ctdna_variants.py fetches its runs straight from SRA with scripts/sra_fetch.py, because anything with a public accession is regenerated rather than mirrored. The other four simulate a library whose true molecule and clone counts are known, so every number they print can be checked rather than admired.

uv run marimo edit notebooks/platforms.py

notebook	what it answers
`notebooks/platforms.py`	every layout above, with a real end-to-end run of two of them
`notebooks/barcode_space.py`	is my barcode long enough? collisions, occupancy, the error budget
`notebooks/refine_diagnostics.py`	the coverage curve, the barcode-rank plot, and where the errors are
`notebooks/exome_capture.py`	duplicates or real molecules? why coordinate deduplication undercounts a capture panel
`notebooks/airr_repertoire.py`	how much of a repertoire is PCR? clonotype counts from reads against molecules
`notebooks/ctdna_variants.py`	how many molecules a variant caller actually gets, on cell-free DNA reference material at known allele frequency (Variant calling: which caller, and what it can possibly see)

Looking at the run#

migec plot co/ -o figs/          # every panel whose table is in co/
migec plot ref/ -o figs/
migec plot asm/ -o figs/ --format pdf

Four of the twenty panels are figures you have already read somewhere else, and they are the ones to look at first:

panel	the question	the failure it shows
`cell_rank`	is my cell calling right?	Cell Ranger’s barcode rank plot, on unique UMIs. No knee means no cells.
`mig_size_spectrum`	is the library over-sequenced?	most molecules shallow while most reads sit in the deep ones
`mig_size_zipf`	is amplification even?	a bent rank curve where Zipf would be straight
`umi_error_rate`	is the barcode error rate believable?	the two estimators part company, which means the barcode space filled
`consensus_quality`	what quality am I allowed to claim?	the boxes flatten at the RT floor, not at the instrument

Every panel is a gnuplot script over a TSV the stage already wrote, so a figure can be redrawn without the FASTQ. Without gnuplot installed the .gp scripts are still written. See plot.

Then what#

The consensus is ordinary FASTQ, so an aligner or a quantifier takes it directly. All four were run against real assemble output; Downstream: what consumes the consensus has the table and the record counts.

minimap2 -ax sr -y  ref.fa asm/S1.consensus.fq.gz | samtools sort -o S1.bam
minibwa map -y -t8  ref.fa asm/S1.consensus.fq.gz | samtools sort -o S1.bam
bwa mem -C          ref.fa asm/S1.consensus.fq.gz | samtools sort -o S1.bam
salmon quant -i tx.idx -l A -r asm/S1.consensus.fq.gz -o quant/   # NumReads = molecules

Never: not alevin, bustools or STARsolo. They deduplicate from a raw barcode read that no longer exists, so running them on a consensus collapses the library twice. Downstream: what consumes the consensus also covers when to align before collapsing instead.

Pipelines#

Two, for two different situations (Pipelines: Nextflow and SLURM has both in full):

Nextflow – integrations/nextflow/ runs the three stages and then continues: --mode ctdna aligns and calls variants, --mode airr calls clonotypes with arda, --mode consensus stops at the consensus. It drops into nf-core/airrflow or any pipeline that hands you FASTQ pairs.
SLURM – integrations/slurm/ is two sbatch templates and a sample sheet, for a cohort where the deliverable is the consensus. Both run as ordinary bash without SLURM, which is how to check a layout before queueing anything.

All three stages are byte-identical at any thread count, so a retry with more cores cannot change a result – which is what makes an escalating retry safe in either.

nextflow run integrations/nextflow --mode ctdna --input 'd/*_R{1,2}.fq.gz' --fasta ref.fa
sbatch --array=1-11 integrations/slurm/migec_array.sbatch samples.tsv