Bring your own UMI#

migec checkout is one way of producing migec’s input. It is not the only one, and for a capture assay it is usually the wrong one: on an exome, a hybrid-capture panel, a ctDNA or an MRD library the UMI is in the index read, so it never appears inside R1 or R2 and there is no pattern for checkout to find. What the sequencing facility hands over is a BAM from fgbio, Picard or the kit vendor’s pipeline, with the UMI already in the RX tag.

That file is a migec input. refine, assemble and subsample group on RX, and where that tag came from is not their business.

The contract#

A record is migec input when it carries these tags. Only RX is required.

tag

what it is

if it is missing

RX

the UMI

refine and assemble refuse the file by name. Nothing to group on

QX

the UMI’s own base qualities, one Phred per base

the correction posterior falls back to the library’s global error rate, exactly as it does for a v1 .mig file. Correction still runs; it just has one fewer piece of evidence

CB, CY

cell barcode and its qualities

the library is treated as one cell, which is what a bulk library is

BC

sample id

the sample is named sample, or whatever --sample says

In a FASTQ those tags live in the comment — everything after the first space or tab in the header — separated by tabs, which is what makes the line a valid SAM record once an aligner appends it (File formats). In a BAM they are ordinary tags. The two are the same input.

From a BAM#

migec refine   umi_tagged.bam -o ref/
migec assemble ref/S1.fq.gz   -o asm/

BAM, SAM and CRAM are recognised from the file itself, not from the name — there is no flag. Behind the scenes migec converts once with samtools into a temporary FASTQ inside the output directory, which is deleted when the stage returns:

samtools collate -u -O in.bam | \
  samtools fastq -n -T RX,QX,CB,CY,BC -1 R1.fq -2 R2.fq -0 R0.fq -s S.fq -

Two things about that command are load-bearing.

Warning

Collate is not optional on an aligned file. samtools fastq -1/-2 pairs by adjacency, so on a coordinate-sorted BAM it writes one molecule’s mate 1 beside another molecule’s mate 2. assemble matches mates strictly by position, so the wrong pair would be consensed and nothing downstream could tell. An unaligned BAM cannot be coordinate-sorted, so migec skips collate there and the record order is preserved; @SQ in the header is what decides, and a header that cannot be read counts as aligned.

Note

The temporary FASTQ is roughly 4x the BAM on disk while the stage runs, and it is written next to your output, not in /tmp — the same siting as the range-partition buckets, for the same reason.

Starting from an index read#

This is the common case and it needs no BAM at all to begin with: the facility gives you R1, R2 and an I1/UMI FASTQ, and the UMI is in the third file. One samtools command puts it where migec reads it, and samtools is already required for the BAM path:

samtools import -1 R1.fq.gz -2 R2.fq.gz --i1 UMI.fq.gz \
    --barcode-tag RX --quality-tag QX -o tagged.bam
migec refine   tagged.bam  -o ref/
migec assemble ref/S1.fq.gz -o asm/

--barcode-tag RX is the whole trick: samtools import writes an index read into BC/QT by default, which is the sample barcode, and the two are not the same thing.

Measured, on a 2,000-molecule simulated library where the UMI was split out into its own FASTQ: this route and checkout → refine on the un-split reads agree exactly — 9,297 reads, 2,193 distinct barcodes, 193 merged, 2,000 molecules from 2,000 injected, on both.

Getting the UMI into RX some other way#

If it is already in a BAM from the vendor’s pipeline, it is already RX and there is nothing to do. Otherwise these write it; all are the vendor’s own tool and migec reimplements none of them.

# a BAM plus the UMI read
fgbio AnnotateBamWithUmis -i in.bam -f UMI.fq.gz -o tagged.bam

# from FASTQ, with a read structure
fgbio FastqToBam -i R1.fq.gz R2.fq.gz UMI.fq.gz --read-structures +T +T +M -o tagged.bam

# Picard
picard FastqToSam F1=R1.fq.gz F2=R2.fq.gz UMI_FASTQ=UMI.fq.gz ... O=tagged.bam

Note: the samtools import recipe above was run here; these three are quoted from their own documentation. fgbio needs a JDK 17+, and on an older JDK it fails with an UnsupportedClassVersionError from htsjdk rather than a version message.

Warning

umi_tools extract puts the UMI in the read name, not in a tag, and migec does not parse read names — the name is only ever copied through. This is deliberate: a name is free-form and a tag is not, so guessing at a name is how one pipeline’s _ separator becomes another’s silent truncation. samtools import -U converts a name-borne UMI into RX, but only for the Illumina/CASAVA convention — on a read_ACGTACGT name it silently writes no tag, which is why migec refuses a file with no RX by name rather than reporting zero molecules.

From a FASTQ somebody else tagged#

Nothing special is needed. If the comment carries RX:Z: after a tab, it works:

@A00123:45:HXXX:1:1101:1000:1000 RX:Z:ACGTACGT<TAB>QX:Z:IIIIIIII
TATCAGAGTAGTGGTATTTCACAGGCGGCCAGCAGGG
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

samtools fastq -T RX,QX writes exactly this, so a round trip through a BAM and back is a no-op as far as migec is concerned.

What this does not change#

Grouping is still on sample + cell + UMI and never on the alignment position. That is the difference between migec and a map-first tool such as fgbio’s GroupReadsByUmi or umi_tools group, and it is measured rather than asserted — see Grouping accuracy: Calib, UMI-tools, fgbio. Reading their file format is not the same as adopting their model: a tool that transports the barcode composes with migec, and a tool that deduplicates on it replaces a stage of migec (Downstream: what consumes the consensus).