Archive OTU table v6: record frameshift repairs as per-read variants

ARCHITECTURE.md

@@ -49,8 +49,11 @@ Ordered roughly by where a reader should start.
   of lost windows on Nanopore reads. Resolution is the one step whose answer depends on the rest
   of the dataset, so a chunked run (`--read-chunk-size`) skips it and `summarise
   --resolve-ambiguous-windows` does it over the combined chunks instead; both go through the
-  shared `resolve_windows()`, and the archive's `reads_with_repaired_deletions` field is what
-  carries the provenance that far.
+  shared `resolve_windows()`. The repairs to each read are recorded as VCF-style variants
+  (`repair_frameshifts_with_variants()`), computed in the prefilter, attached to each window
+  sequence as `unaligned_sequence_variants`, and written to the archive's
+  `read_unaligned_sequence_variants` field; that is also what carries the provenance of a
+  repair-inserted `N` as far as `summarise`.
 - **`prefilter_pad.py`** — `PrefilterPadder` plus the pure helpers `window_alignment_positions()`
   and `pad_aligned_sequence()`. Emits a prefilter FASTA where every on-target sequence is
   padded to a fixed length (30aa + 20aa window + 30aa, `X`-padded) with the window in a
@@ -62,9 +65,10 @@ Ordered roughly by where a reader should start.
   The OTU data model. `OtuTableEntry` is one (marker, sample, window sequence, count,
   coverage, taxonomy) row; `OtuTable` is the plain TSV form; `ArchiveOtuTable` is the richer
   JSON form (retains per-read data, enabling `renew`); collection classes provide streaming
-  iteration over many tables (`StreamingOtuTableCollection`). Archive version 5 writes the
+  iteration over many tables (`StreamingOtuTableCollection`). Archive version 6 writes the
   OTUs column-wise, hoists fields that are constant across every OTU into `constant_fields`,
-  and dereplicates read sequences into a shared `reads` list; versions 1-4 are still read.
+  and dereplicates read sequences into a shared `reads` list; versions 1-4 are still read,
+  but not version 5, which only development versions wrote.
   In memory the table is a list of rows either way, so only `ArchiveOtuTable` itself cares.
 - **`condense.py`** — `Condenser` turns per-marker OTUs into a `CondensedCommunityProfile`
   (a coverage tree of `WordNode`s), applying trimmed means and genus/species

CHANGELOG.md

@@ -2,11 +2,13 @@
 
 * `pipe`/`condense`: Add joint SingleM + weebill taxonomic profiling (`--joint`), which uses NNLS deconvolution against a weebill coverage profile instead of the default EM-based condense, with sylph-only species injection, species-level pinning (`--joint-pin-weebill-species`) and a novel-lineage budget (`--joint-novel-budget`) to avoid fabricating novel lineages in communities weebill already fully explains. `pipe`/`renew` run weebill automatically whenever the metapackage bundles a weebill database; disable with `--no-weebill`. Add `--output-weebill-sketch`/`renew --input-weebill-sketch` to save/reuse a weebill read sketch without keeping the raw reads. Metapackage format bumped to v7 to bundle a two-stage weebill database (`metapackage --weebill-db`/`--weebill-c`).
 * `pipe`: Remove `--weebill-injection`; the joint SingleM + weebill deconvolution is now the only way weebill is integrated. `pipe` now raises an error if a taxonomic profile is requested and the metapackage bundles no weebill database, unless `--no-weebill` is given.
-* Archive OTU tables are now written at version 5, which stores the OTUs column-wise, stores fields that are constant across every OTU once rather than per OTU, and dereplicates read sequences into a shared list. Versions 1-4 are still read, and `ArchiveOtuTable` presents the same list of rows whichever version it read, but code that parses the JSON itself rather than going through `ArchiveOtuTable` needs updating - see the Archive OTU table entry in the glossary for the layout.
-* Archive OTU tables version 5 also add a `reads_with_repaired_deletions` field, recording which of an OTU's reads have an ambiguous base that frameshift repair inserted.
+* Archive OTU tables are now written at version 6, which stores the OTUs column-wise, stores fields that are constant across every OTU once rather than per OTU, and dereplicates read sequences into a shared list. Versions 1-4 are still read, and `ArchiveOtuTable` presents the same list of rows whichever version it read, but code that parses the JSON itself rather than going through `ArchiveOtuTable` needs updating - see the Archive OTU table entry in the glossary for the layout.
+* Archive OTU tables version 6 also add a `read_unaligned_sequence_variants` field. `read_unaligned_sequences` now holds each read as modified by frameshift repair, and this field records the modifications as VCF-style `[POS, REF, ALT]` variants, so that the read as sequenced can be recreated. Extended (`--output-extras`) OTU tables gain the same column, written as `POS:REF>ALT`. With `--context-window`, the shortening of the read to the region around the hit is not recorded as a variant. Version 5 archive OTU tables, written only by development versions, are not read.
 * `pipe`: When processing one chunk of a dataset with `--read-chunk-size`/`--read-chunk-number`, the ambiguous bases left by frameshift repair are no longer resolved within the chunk. Which window a read ended up with otherwise depended on what else was in that chunk, so the chunks did not combine into the result the whole dataset would have given.
 * `pipe`: Fix the `nucleotides_aligned` of an archive OTU table being ordered independently of `read_names` and `read_unaligned_sequences`, which were sorted by read name while it was not. The two orders only differ when the reads of one OTU do not all have the same aligned length, which takes an insert column inside the window's span that is a gap in some of them, so this was rarely visible - but `renew` reads the three as parallel.
-* `summarise`: Add `--resolve-ambiguous-windows` (and `--max-frameshift-repair-divergence`), which resolves those bases while archive OTU tables are combined, where every chunk is present. Use it when combining the chunks of a dataset processed with `pipe --read-chunk-size`.
+* `summarise`: Add `--resolve-ambiguous-windows` (and `--max-frameshift-repair-divergence`), which resolves those bases while archive OTU tables are combined, where every chunk is present. Use it when combining the chunks of a dataset processed with `pipe --read-chunk-size`. A window resolved to one already present takes that window's taxonomic assignment, as it would have had the whole dataset been processed at once.
+* `pipe`: `--read-chunk-number` is incompatible with running weebill, which would profile the whole of the reads rather than the chunk. Pass `--no-weebill`.
+* `condense`/`pipe`: A sample in which weebill detects species but which has no marker OTUs is now profiled from its weebill coverage, rather than left out of the taxonomic profile.
 * `summarise`: `--collapse-to-sample-name` and `--collapse-paired-with-unpaired-archive-otu-table` now combine archive OTU tables of different versions, rather than rejecting them. The output takes the newest version's fields, and OTUs from an older archive have no value for the fields that version lacked. Previously an archive made before this release could not be combined with one made after it.
 * `pipe`: Frameshift repair, which uses the frameshifts DIAMOND reports during the prefilter to restore the reading frame of each read before it is aligned to the HMM, so reads carrying single base indels still yield a window, is now on by default. Substantially improves window recovery on Nanopore data, where indels rather than substitutions are the dominant error. Where a repaired deletion leaves a base of unknown identity, it is taken from the most abundant window within `--max-frameshift-repair-divergence` mismatches, and only for windows where repair actually inserted that base, not a pre-existing ambiguous base in the raw read. Disable with `--no-repair-frameshifts`. 
 * `condense`: Fix a bug where, in a multi-sample run, a weebill profile with hits for exactly one named sample had those hits reused for every other sample instead of being skipped for them.

docs/Glossary.md

@@ -38,24 +38,32 @@ The extended OTU table form generated with the `--output-extras` option to the [
   1. read_names - the names of the reads which encode the OTU sequence
   2. nucleotides_aligned - the number of nucleotides which aligned to the window (usually 60, but can be more or less if there are gaps or inserts)
   3. taxonomy_by_known? - whether the taxonomy of the OTU was determined by known genomes (TRUE) or by the reads themselves (FALSE). Currently this is a disused column and is always marked FALSE.
-  4. read_unaligned_sequences - the raw sequences of the reads which encode the OTU sequence
+  4. read_unaligned_sequences - the sequences of the reads which encode the OTU sequence, after any frameshift repair
   5. equal_best_hit_taxonomies - the taxonomies of the best hits to the OTU sequence, if there are multiple equally good hits. This is a JSON array of strings.
+  6. taxonomy_assignment_method - the method used to assign the taxonomy of the OTU
+  7. read_unaligned_sequence_variants - the modifications made to each read by frameshift repair (see archive OTU table, below), with one space-separated entry per read in the order of read_names. An entry is `.` for a read that was not modified, and otherwise the read's variants written as `POS:REF>ALT` and separated by commas, e.g. `121:AN>A,150:T>TT`.
 
 
 ## **Archive OTU table**
 Similar to an extended form OTU table, but in JSON form for machine readability and with formatting version recorded. The [renew](/tools/renew) subcommand which re-analyses a dataset requires this format of OTU table rather than the default tab-separated OTU table format. The canonical file extension for SingleM packages is `.json`.
 
-From version 5, archive OTU tables are laid out to reduce the repetition that dominates the file, so the same OTUs take less space on disk:
+From version 6, archive OTU tables are laid out to reduce the repetition that dominates the file, so the same OTUs take less space on disk:
 
 * `otus` maps each field name to the list of that field's values across the OTUs, rather than being a list of one row per OTU.
 * Fields whose value is identical in every OTU are stored once in `constant_fields` and left out of `otus`. `num_otus` gives the number of OTUs, since the columns alone do not when every field is constant.
 * `reads` holds the unaligned read sequences, dereplicated: a read hitting two markers is stored once rather than once per marker. The `read_unaligned_sequences` field then holds, for each read of an OTU, either `0` for "the next read of `reads` not yet used" or `n` for a repeat of `reads[n - 1]`.
 
-Version 5 also records frameshift repairs. A missing read base is stored as `N` in the OTU `sequence`; `reads_with_repaired_deletions` stores the zero-based indices of the affected reads in `read_names`. This distinguishes repair-inserted `N`s from ambiguous bases already present in the reads. The value is `null` when no read is affected. Repairs that remove an extra base need no additional metadata.
+Version 6 also records the modifications that frameshift repair makes to each read. The `read_unaligned_sequences` are the reads *after* repair: where a base is missing from a read, an `N` is inserted, and where a read has an extra base, it is removed. The `read_unaligned_sequence_variants` field records, for each read of an OTU, how to turn that stored sequence back into the read as sequenced. It is `null` for a read that was not modified, and is `null` for the whole OTU when none of its reads were. Otherwise it is a list of VCF-style variants `[POS, REF, ALT]`, in which the stored sequence plays the part of the VCF reference:
 
-This field is internal and is omitted from extended TSV OTU tables. [summarise](/tools/summarise) uses it with `--resolve-ambiguous-windows` to resolve deferred `N`s after combining chunks produced by [pipe](/tools/pipe) `--read-chunk-size`.
+* `POS` is the 1-based position in the stored sequence at which the variant starts.
+* `REF` is the stored sequence from `POS` onwards, and `ALT` what the read as sequenced has in its place.
+* As in VCF, each variant includes the base before the edit (or the base after it, when the edit is at the very start of the read). So `[7, "CN", "C"]` means an `N` was inserted after the `C` at position 7, and `[7, "C", "CT"]` means a `T` after that `C` was removed. Edits close enough for these bases to overlap are combined into one variant, so variants never overlap and are listed in order of `POS`.
 
-Versions 1 to 4, which store `otus` as a list of rows with every field present in each, are still read. Code that reads archive OTU tables through `singlem.archive_otu_table.ArchiveOtuTable` sees a list of rows regardless of the version on disk; code that parses the JSON itself needs to handle the layout of the version it is given.
+Replacing each `REF` with its `ALT` recreates the read as sequenced. When reading archive OTU tables through `singlem.archive_otu_table.ArchiveOtuTable`, `read_original_unaligned_sequences()` does this. With `pipe --context-window`, only the region around each hit is stored rather than the whole read, and this shortening is not recorded as a variant.
+
+The variants also show which reads had an `N` inserted by repair (a variant whose `REF` has more `N`s than its `ALT`), as opposed to ambiguous bases already present in the read. [summarise](/tools/summarise) uses that with `--resolve-ambiguous-windows` to resolve deferred `N`s in OTU sequences after combining chunks produced by [pipe](/tools/pipe) `--read-chunk-size`.
+
+Version 5 archive OTU tables were only ever written by development versions of SingleM, and are not supported. Versions 1 to 4, which store `otus` as a list of rows with every field present in each, are still read. Code that reads archive OTU tables through `singlem.archive_otu_table.ArchiveOtuTable` sees a list of rows regardless of the version on disk; code that parses the JSON itself needs to handle the layout of the version it is given.
 
 ## **Coverage, unfilled coverage and filled coverage**
 First, to define *coverage*. Coverage of a single species is defined in the standard way as in e.g. [CoverM](https://github.com/wwood/CoverM). It is the average number of reads covering each position in that species' genome. SingleM (and Lyrebird) do not directly calculate this via mapping, but instead assume that the coverage of a species' single copy marker genes is representative of the coverage of the whole genome. In particular, it is estimated from the reads that encode sequences recruited to the highly conserved marker gene windows.