checkout – find the barcode and cut it out#
migec checkout finds the barcode pattern in each read, extracts the sample tag and the UMI,
trims the synthetic sequence away, and puts the barcode in the read header.
migec checkout reads.fq.gz --barcodes barcodes.txt --out out/
migec checkout R1.fq.gz R2.fq.gz --barcodes barcodes.txt --out out/ -t 8
Barcode patterns#
The pattern grammar is MIGEC’s, so published barcode tables work unchanged:
S1 aaACTcagtggtatcaacgcagagtNNNNtNNNNtNNNN
S2 aaAGAcagtggtatcaacgcagagtNNNNtNNNNtNNNN
symbol |
meaning |
|---|---|
|
a scored position, matched exactly |
IUPAC ( |
a scored position with a degenerate set: |
lowercase |
a scored position at half weight — the adapter, where a mismatch is expected |
|
a UMI position: captured, never scored |
|
a wildcard: neither scored nor captured |
N always means UMI here, never IUPAC “any base” — use . for an uncaptured wildcard.
UMI runs need not be contiguous. NNNNtNNNNtNNNN captures twelve bases into one UMI, with the
lowercase t spacers scored at half weight — which is how the published MIGEC tables are
written.
migec sheet barcodes.txt prints what each row will extract without running anything.
How a match is accepted#
Not by counting mismatches. For each scored position, with IUPAC set \(S\) of size \(m\), observed base \(b\) and error probability \(e\):
which is the log-likelihood ratio of “the tag is here” against “this is random sequence”. \(w\) is 1.0 for uppercase and 0.5 for lowercase. At \(m = 1\) a match is worth +2.00 bits, a mismatch −9.55 bits at Q30 and −0.60 bits at Q2.
So a mismatch on a bad base is nearly free and a mismatch on a good base is fatal. That is what
MIGEC’s “good mismatch / bad mismatch” counting was reaching for, done continuously — and without
its two defects, both of which are reproduced as tests here: v1 read the quality string from the
start of the read rather than from the match offset, and a dangling else meant low-quality
mismatches were never counted at all.
The acceptance threshold defaults to a Bonferroni bound over the offsets actually scanned, \(\log_2(n_\text{offsets} \cdot n_\text{patterns} / \alpha)\).
Note
Reads are not i.i.d. uniform ACGT — shared primers and composition bias violate the null badly. Treat the default threshold as a starting point and calibrate it against shuffled decoy patterns on your own data.
Ambiguous is not the same as unmatched#
A read whose best sample tag does not beat the runner-up by the placement margin (min_margin,
5 bits, an internal parameter rather than a flag) is reported as
ambiguous rather than assigned to one of them. That is a different diagnosis from unmatched:
ambiguous means the barcodes in the sheet are too close together, unmatched means the pattern is
wrong or absent. One counter cannot say both, so there are two.
Paired input and strand normalisation#
Give a second mate and checkout looks for the tag in R1 first, then — only if R1 came up empty, so the cost falls on reads that would otherwise be discarded — in R2. When it turns up there the pair is swapped, so the output R1 always carries the tag. For single-end input the same fallback reverse-complements the read.
This is not a convenience. Amplicon libraries are sequenced in both orientations, and a MIG holding
both orientations of one molecule loses half its reads at consensus while nothing upstream reports
it. The count of reads that were flipped is in the report and in normalised.
The mate is passed through whole, including when a dual-end (slave) pattern matched inside it:
trimming it would need its own payload_begin carried alongside, and by then the mate’s barcode
bases are already in the UMI. Both mates carry the RX/QX/BC tags, because a
downstream tool that sees only one of them cannot group the pair.
Trimming#
--trim pattern (the default) drops everything up to and including the matched pattern — the
adapter, the sample tag and the UMI. What is left is exactly the biological payload. This is what
you want before alignment: the tag is synthetic and will soft-clip at best, mismap at worst.
--trim none keeps the read whole and puts the UMI in the header only.
Barcode transfer to the header#
The barcode travels on in SAM-style tags:
@r0 RX:Z:GCTAAAGACAAT QX:Z:IIIIIIIIIIII BC:Z:S1
TACATAACATACACGTCAGCACGAAACTTGTTGGCCCAGTGTGAATCGCTT
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
RX/QX are the UMI and its qualities (read by fgbio, Picard and umi_tools); BC is the
sample barcode. Tags are separated by TAB, with one space after the read name.
Warning
The TABs are not cosmetic. bwa mem -C and minimap2 -y copy the FASTQ comment verbatim
into the SAM record, so it has to be SAM-conformant there or the resulting BAM is malformed.
Warning
dnaio — used by arda’s rnaseq module — drops the comment entirely. Anything a
downstream Python tool must see has to be in the read name.
Quality filtering of the UMI#
--min-umi-quality defaults to 0, i.e. nothing is dropped. MIGEC used 15 and MAGERI 20, both
as hard drops. The default here is deliberately different: a low-quality UMI base is a reason to be
less certain about which molecule a read belongs to, not a reason to throw the read away. The
correction step can usually recover it, and a molecule seen three times is information.
Output#
file |
content |
|---|---|
|
trimmed reads with barcodes in the header (single-end) |
|
the same, for paired input |
|
reads matching no pattern ( |
|
one row per sample: yields, UMI statistics, correction |
|
reads and distinct UMIs per power-of-two MIG size |
|
per-position base usage, entropy, information, collision |
|
reads per observed i7 x i5 combination, and which were ordered |
|
reads per lane and tile: the run’s yield map |
|
everything above, machine-readable |
See UMI statistics for what the last two contain and how to read them, and
Speed and memory for --threads and what the run costs in time and memory.
Writing the partition instead: --mig#
--mig writes <sample>.<bbb>.mig buckets in place of the per-sample FASTQ. The reads come
out range-partitioned on the barcode — the same partition, on the same key, that
assemble – one consensus per molecule otherwise builds for itself in its first pass — so assemble reads them
directly and skips that pass:
migec checkout reads.fq.gz -b barcodes.txt -o out --mig
migec assemble out/S1.000.mig -o asm # one bucket names the whole partition
Measured on the 500,000-read benchmark corpus, four samples, four threads: 1.16 s → 0.98 s end
to end for the identical 124,878 molecules. checkout pays 0.06 s of it and assemble saves
0.25 s, which is its whole partition pass.
FASTQ stays the default, and should. A .mig file is a migec intermediate that nothing else
reads: every aligner, every pipeline in Downstream: what consumes the consensus, and every example here speak FASTQ.
Reach for --mig when the reads are going straight to assemble and nowhere else.
Note: one bucket file names the whole partition — assemble collects that sample’s siblings
beside it. A directory holding more than one sample’s buckets is refused by name rather than
assembled together, because a UMI repeats across samples by design and merging two samples’ reads
would create molecules that never existed.
Note: the open-file budget is for the run rather than for each sample, so the bucket count falls as the sheet grows — 256 buckets for one sample, 64 each for four, two each for a 96-plex. A sample of a 96-plex sheet also holds a 96th of the reads, so this is proportionate rather than a compromise.
Note: the .mig header carries no quality calibration. It is fitted from the whole run and the
buckets are opened while the run is still going; checkout.json carries the fit, and a table
written before it was measured would be worse than an absent one.
Note: paired input puts both mates in one bucket file rather than in two files, since the
record holds them together. assemble consenses mate 1, which is what it does from FASTQ as
well — there the second mate is simply a file you did not hand it.
Note: unmatched reads stay FASTQ (--write-unmatched). They carry no barcode, so there is no
bucket to put them in.
Index hopping, from the header#
On a patterned flowcell a free index primer can prime a neighbouring cluster, so a molecule from one sample is read carrying its own i7 and another sample’s i5. The read then lands in that other sample and looks exactly like one of its reads: nothing in the sequence, the barcode or the count says otherwise, which is why a per-sample yield can never see it.
The instrument already wrote the evidence. checkout reads the index pair out of the read
header — the last field of 1:N:0:ATCACG+CGTGAT — for every read, matched or not, and
writes the contingency table:
i7 i5 reads share_of_i7 share_of_i5 declared
ATCACG CGTGAT 9803 0.980 0.980 1
CGATGT ACATCG 9797 0.980 0.980 1
ATCACG ACATCG 200 0.020 0.020 0
CGATGT CGTGAT 200 0.020 0.020 0
A combination counts as declared when it holds at least 5% of the reads of its own i7 and of its own i5. The sample sheet migec is given carries the in-line barcode, not the index pair, so the declared set is inferred — and the gap is wide enough to infer it from: hopping runs at 0.1–2% while a declared combination is the bulk of its own index. The raw counts are in the table so the inference can be disagreed with.
Note: it matters most where it is smallest. At 0.1% hopping, a 1% variant in a deeply sequenced sample contaminates its neighbour at 1e-5 — which is exactly the level a rare-variant caller is asked to believe.
Never: a single-indexed run is not estimable, and that is not the same as zero. With one index
there are no combinations, so nothing can be off-diagonal; estimable is false and the rate is
not reported as a finding.
Where on the flowcell#
The same header carries the lane and the tile — instrument:run:flowcell:lane:tile:x:y — and
checkout.tiles.tsv is one row per (lane, tile) with its share of its lane. A tile is a physical
patch of the flowcell, so a bubble, a dead tile, an edge effect and an underloaded lane all show
here and in no other table migec writes: a run with one starved tile and a healthy run have the same
read count, the same barcode statistics and the same molecule count.
Note: the deeper spatial question is whether two reads of one molecule are the same cluster read
twice. Optical duplicates on an unpatterned flowcell, and ExAmp pad-hopping duplicates on a
patterned one, are not independent observations, and the consensus posterior adds one
log-likelihood per read as though they were. That scan needs the reads of a molecule together and
the pixel coordinates with them, which would cost 12 bytes a read in the one stage whose memory
bound is the claim — so it lives in scripts/diagnose.py, at Picard’s own distances (100 px
unpatterned, 2500 px patterned).
Never: an SRA-normalised header has no coordinates at all (@SRR1763769.1 1/2), and a run whose
headers did not survive reports no map rather than a map of one tile that reads as a
single-tile flowcell.
What the reported Phred is actually worth#
The pattern’s own constant bases are known sequence — the adapter and the sample tag — so a
disagreement there is an instrument error and nothing else. That makes them the only free
calibration standard in the read, and checkout counts match/mismatch at every unambiguous
scored position, indexed by the reported Phred.
reported Phred is worth 1.04x its nominal error, measured on 46,289,536 constant pattern bases
the fit's intercept is 3.9e-03 per base -- the SYNTHESISED anchor's own defect rate,
not a sequencing floor
The table is fitted as ê(q) = ε_qi + a · 10^(−q/10), weighted by how many bases carried each
Q — which matters, because RTA3 emits about four distinct Q values and an unweighted fit would let
one seen a hundred times outvote one seen a billion.
Warning
The intercept is not a sequencing floor. The standard being measured against is a
synthesised oligo, and oligo synthesis carries roughly one defect per 200–500 bases. On
SRR1763769 the intercept comes out at 3.9·10⁻³, spread evenly over all 23 anchor positions
with none polymorphic — against an independently measured 0.55% rate of one-base-short
barcodes from failed couplings in the same oligo. Same order, same cause.
So it is reported as a diagnostic of the primer and deliberately left out of the calibrated
error. error(q) applies the slope only; folding the intercept in would add 4·10⁻³ to every
base likelihood in the pipeline on the strength of the primer’s quality.
A pattern position whose mismatch rate sits far above the median is dropped before fitting: it is
polymorphic, or the pattern is wrong about it, and either way it is not measuring the instrument.
checkout.pattern_positions.tsv shows every position and whether it was used;
checkout.quality_calibration.tsv is the per-Q table.
Dual-end barcodes#
Column 3 of a MIGEC barcode table is the slave barcode: a second pattern on the other mate whose captured positions extend the UMI rather than starting a new one. MAGERI’s design, quoted from its Methods, is twelve bases at each end of the molecule:
S1 NNNNNNNNNNNNTGACT AGTCANNNNNNNNNNNN
migec checkout R1.fq.gz R2.fq.gz -b barcodes.txt -o out/
Never: Both halves must match or the read is unmatched. Accepting the master alone would emit 12 nt UMIs beside 24 nt ones, and every collision estimate downstream would then be computed over two barcode spaces at once.
--max-offset: where the pattern may start#
The default is automatic and this flag should not be passed. A leading ^, a slice list, a
read structure and a pattern with nothing to score all anchor at the first base; anything with an
adapter to place it gets a free scan. -1 forces the free scan and 0 forces the anchor.
A dual-end design needs the anchor, and the reason is not convenience:
A five-base handle like TGACT is worth 10 bits. The acceptance bar is a Bonferroni bound over
the offsets scanned — log2(offsets × patterns / α) — which over a 77 nt read is 12.6 bits. So a
free scan correctly refuses it: TGACT occurs by chance about every kilobase and cannot be
placed. Anchored at offset 0 the bar is 6.6 bits, the handle clears it, and the placement is
determined by the chemistry rather than by the sequence.
Note
The bound is charged for the offsets actually scanned. Billing an anchored scan for the sixty offsets it never tries refuses every read of a design that is perfectly well determined — which is what it did until this was fixed.