Layouts: where the barcode is ============================= There is exactly one thing migec has to be told, and everything after it -- correction, consensus, the quality cap -- is the same three commands whatever the answer. Four ways to say it, in the order you should reach for them. If the barcode has already been extracted -- a capture, exome or ctDNA kit puts it in the index read, so it reaches you in the ``RX`` tag of a BAM rather than inside R1 -- then there is no layout to declare and ``checkout`` is not the stage you want. See :doc:`Bring your own UMI `. 1. A position ------------- Most libraries put the barcode at a fixed offset in one read. That is the primary mode, and it needs neither a sample sheet nor an anchor: .. code-block:: bash 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 slice Two spellings, one meaning: .. list-table:: :header-rows: 1 :widths: 30 30 40 * - pattern - slice - what it is * - ``^NNNNNNNN`` - ``0:8`` - an 8 nt UMI at the first base * - ``^NNNN.NNNNN`` - ``0:4,5:10`` - a 9 nt UMI split by one skipped base * - ``^XXXXXXXXXXXXXXXXNNNNNNNNNN`` - ``cell:0:16,16:26`` - a 16 nt cell barcode then a 10 nt UMI (10x) * - ``^NNNNXNNN`` - ``0:4,cell:4:5,5:8`` - a UMI interrupted by one cell-barcode base In a pattern, ``N`` is a UMI base, ``X`` a cell-barcode base, ``.`` a base that is skipped -- neither scored nor captured -- and anything else (``ACGT`` or any IUPAC symbol) is constant sequence that gets scored. Slices are **half-open and 0-based**, like Python's: ``0:8`` is eight bases and the next slice may start at 8. Each is a UMI slice unless it is prefixed ``cell:``. Slices must be in increasing order and must not overlap, because one base belongs to one barcode. Never: ``umi_tools`` spells a cell-barcode position ``C``, and ``C`` is **cytosine** here. Pasting one of its patterns would otherwise compile into a demand for a run of literal cytosines that matches nothing, so it is refused by name with the translation to ``X`` rather than accepted. Anchoring ~~~~~~~~~ A leading ``^`` says the barcode starts at the first base. Every slice list says the same thing by construction, since a position is only a position if it is measured from somewhere. Both set ``--max-offset 0``, and a layout with nothing to score is anchored automatically even without the caret -- which is why the flag no longer appears in any of the examples. Never: this is not a convenience. Placement is a hypothesis test, and a pattern with no constant sequence supplies no evidence for it. Asked to scan freely, ``compile()`` refuses rather than picking an offset; asked to scan a 5 nt dual-end handle, it refuses too, because ``TGACT`` occurs by chance about every kilobase and the bar for a free scan is ``log2(offsets/alpha)`` bits, which five bases cannot pay. Anchored, there is only one place to be and the bar does not apply. See :doc:`checkout` for the arithmetic. 2. A named preset ----------------- .. code-block:: bash migec sheet --presets # all of them, and their sources migec checkout R1.fq.gz R2.fq.gz --preset 10x-v2 -o out/ A preset places the barcode and stops there. What the *experiment* implies -- how many reads a consensus needs, which pre-amplification floor applies, whether the reads under one barcode are even co-terminal -- is the other axis, and ``migec sheet --assay`` prints it as a paste-ready recipe. Eight profiles: ``airr``, ``amplicon`` (a targeted PCR panel, *not* an alias of ``airr``), ``exome``, ``ctdna``, ``mrd``, ``rnaseq``, ``10x-gex``, ``10x-vdj``. See :doc:`assays`. .. list-table:: :header-rows: 1 :widths: 14 34 52 * - preset - layout - what it is, and where the layout is written down * - ``umi`` - ``^NNNNNNNN`` - generic inline UMI. Change the run length, or write the slice. * - ``migec`` - ``cagtggtatcaacgcagagtNNNNtNNNNtNNNN`` - MIGEC 5'-RACE RepSeq: the SMART adapter then a 12 nt UMI split by two spacers. ``misc/barcodes.txt`` of MIGEC 1.2.9, tag ``v1-final``. Prefix a sample tag per row to demultiplex. * - ``primerid`` - ``NNNNNNNNNcagtttaacttttgggccatcca`` - HIV-1 Primer ID amplicon as used by MAGERI. Recovered by ``migec suggest`` from ``SRR1763769``; the primer places it, so the scan stays free. * - ``duplex`` - ``^NNNNNNNNNNNN.....`` on both mates - duplex sequencing: a 12 nt UMI and a 5 nt spacer per mate, 24 nt together. * - ``10x`` - ``^XXXXXXXXXXXXXXXXNNNNNNNNNNNN`` - 10x Chromium 3' v3/v3.1: 16 nt cell barcode, 12 nt UMI, on R1. * - ``10x-v2`` - ``^XXXXXXXXXXXXXXXXNNNNNNNNNN`` - 10x Chromium 3' v2 and 5' v1/v2: 16 nt cell barcode, 10 nt UMI. * - ``tso500`` - ``^NNNNN.....`` on R1 only - Illumina TSO500 ctDNA. The fgbio read structure is ``5M5S+T +T`` -- R2 is all template. See the warning below. * - ``smarter-umi`` - ``^NNNNNNNNNNGGG`` - SMARTer template-switching RNA-seq: a 10 nt inline UMI, then the ``GGG`` the template switch leaves behind. ``ncgr/UMI-analysis``, whose quality filter reads offset 0 length 10. Measured on 10x's own ``sc5p_v2_hs_PBMC_1k`` VDJ-T run with ``--preset 10x-v2``: **100% of 3,155,166 reads assigned**, 221,024 barcodes at 14.28 reads each, an effective UMI length of 9.97 of 10, and 813 cells called by OrdMag from the 305,702 molecules on R2. ``SOURCES.md`` carries the fetch command and the rest of the run. Warning: the ``duplex`` preset extracts the tags and emits **single-strand** consensuses. Pairing the two strands of a molecule into a duplex consensus is not implemented, so no duplex error rate should be quoted from this output. Never: **a 5 nt UMI does not identify a molecule, and TSO500's does not pretend to.** 4^5 is 1,024 barcodes against the tens of thousands of fragments a ctDNA panel region carries, so the space is saturated by construction and the birthday bound says most barcodes are shared. TSO500's own pipeline resolves that by grouping on the UMI **and the mapping position** (``fgbio GroupReadsByUmi``, which runs after alignment); migec groups on the barcode, before any alignment exists, so it cannot. It will report the space as saturated, set ``err_unreliable``, and warn -- and on this chemistry that warning is the correct answer, not a threshold to raise. Use migec here to extract and tag; do the grouping position-aware, downstream. 3. A read structure ------------------- fgbio, Picard, samtools and the TSO500 pipelines describe a layout as a *read structure*, and migec takes them verbatim: ``M`` a molecular barcode, ``B`` a sample/cell barcode, ``S`` a skip, ``T`` template. .. code-block:: bash migec checkout R1.fq.gz R2.fq.gz --read-structure 5M5S+T -o out/ # TSO500: `5M5S+T +T` migec checkout R1.fq.gz R2.fq.gz --read-structure 12M5S+T --read-structure2 12M5S+T -o out/ .. list-table:: :header-rows: 1 :widths: 20 44 36 * - structure - pattern - platform * - ``5M5S+T`` - ``NNNNN.....`` - TSO500 * - ``16B10M+T`` - ``XXXXXXXXXXXXXXXXNNNNNNNNNN`` - 10x 5' * - ``8M+T`` - ``NNNNNNNN`` - a plain inline UMI The pattern stops at the first template segment, because everything after it is payload and migec trims to exactly that point. ``+`` means "the rest of the read" and is valid only on the last segment; ``+M`` is refused, since an unbounded barcode has no length for the collision arithmetic to use. A read structure is positional by definition, so it carries its own anchor. Translating from zUMIs ~~~~~~~~~~~~~~~~~~~~~~ zUMIs -- and so NASC-seq2, which drives it -- writes the layout as a ``base_definition``: .. code-block:: yaml file1: base_definition: - UMI(12-19) - cDNA(23-200) find_pattern: ATTGCGCAATG;2 file2: base_definition: - cDNA(1-200) - BC(201-220) Never: **zUMIs ranges are 1-based and inclusive; migec slices are 0-based and half-open.** They are not the same numbers and they are not off by a constant either -- ``UMI(12-19)`` is eight bases at offsets 11 through 18, which is ``11:19`` here. Subtract one from the start and leave the stop alone: .. list-table:: :header-rows: 1 :widths: 30 30 40 * - zUMIs - migec - note * - ``UMI(12-19)`` - ``11:19`` - 8 nt; start-1, stop unchanged * - ``BC(1-6,20-26)`` - ``cell:0:6,cell:19:26`` - two ranges, one barcode -- migec concatenates them the same way * - ``cDNA(23-200)`` - (nothing) - the payload; migec trims to the end of the pattern and keeps the rest ``find_pattern: ATTGCGCAATG;2`` is a constant anchor with two mismatches allowed. Write it into the pattern in lowercase -- ``attgcgcaatg`` is scored at half weight, which is the soft-match region -- and drop the slice form, since a pattern that carries its own anchor does not need to be anchored at a position. 4. A barcode table ------------------ For many samples in one file. This is MIGEC's own ``barcodes.txt``, read verbatim, which is why the published tables run unchanged: .. code-block:: text S1 aaACTcagtggtatcaacgcagagtNNNNtNNNNtNNNN S2 aaAGAcagtggtatcaacgcagagtNNNNtNNNNtNNNN Uppercase is matched exactly, lowercase is the fuzzy adapter (scored at half weight), and the ``t``\ s between UMI runs are pattern bases, not barcode -- the UMI is 12 nt, so the barcode space is ``4^12`` and not ``4^14``. Column 3 is the *slave* pattern, on the other mate, whose captured positions **extend** the UMI: .. code-block:: text S1 NNNNNNNNNNNNtgact agtcaNNNNNNNNNNNN 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. Rows may share a sample id -- that is how a sample sequenced with more than one tag is declared, and one output file is written per *sample id*, never per row. .. code-block:: bash migec sheet barcodes.txt # what will each row extract, before anything runs If you do not know the layout ----------------------------- Do not guess: .. code-block:: bash migec suggest reads.fq.gz It segments the per-cycle base composition into UMI, constant and payload runs and prints a paste-ready pattern. It recovered the 9 nt + ``CAGTTTAACTTTTGGGCCAT`` layout of ``SRR1763769`` unaided. Note: it stops at the last *constant* run -- composition alone cannot tell a UMI from diverse payload, only the anchor can, and it says so when that is all it found. See :doc:`suggest`.