From ecb6aeb257e136133693e1999104ea2c9e173d80 Mon Sep 17 00:00:00 2001 From: Ben Woodcroft Date: Wed, 28 Oct 2026 12:30:42 +1000 Subject: [PATCH 1/4] Archive OTU table v6: record frameshift repairs as per-read variants Archive OTU tables now store each read after frameshift repair, and a new read_unaligned_sequence_variants field records the repairs as VCF-style [POS, REF, ALT] variants (the stored read being the reference), so the read as sequenced can be recreated. With --context-window, shortening the read to the region around the hit is not recorded. This replaces reads_with_repaired_deletions; version 5 tables are rejected. Variants are computed once in the DIAMOND prefilter and travel with each window sequence, replacing the scattered repaired-deletion bookkeeping in pipe. This also fixes reverse reads of a pair being given their mate's sequence in the paired DIAMOND-taxonomy path. renew keeps the variants, and summarise --resolve-ambiguous-windows uses them to find repair-inserted Ns. Extended TSV OTU tables gain the column, written as POS:REF>ALT per read, or '.' for an unmodified read. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019URUjTJM4Se2CmgGzPmPs5 --- ARCHITECTURE.md | 12 ++- CHANGELOG.md | 4 +- docs/Glossary.md | 18 +++- singlem/archive_otu_table.py | 119 ++++++++++++-------- singlem/diamond_spkg_searcher.py | 58 ++++++---- singlem/frameshift_repair.py | 168 ++++++++++++++++++++++++++--- singlem/genome_fasta_mux.py | 4 +- singlem/otu_table.py | 13 ++- singlem/pipe.py | 148 +++++++++---------------- singlem/pipe_sequence_extractor.py | 7 ++ singlem/renew.py | 14 ++- singlem/sequence_classes.py | 11 +- singlem/summariser.py | 77 +++++++------ test/test_appraise.py | 12 +-- test/test_archive_otu_table.py | 41 +++++-- test/test_diamond_spkg_searcher.py | 34 ++++-- test/test_frameshift_repair.py | 102 ++++++++++++++++-- test/test_pipe.py | 137 +++++++++++++++++++---- test/test_summariser.py | 49 ++++++--- 19 files changed, 724 insertions(+), 304 deletions(-) diff --git a/ARCHITECTURE.md b/ARCHITECTURE.md index 509c9de2..f8ae400b 100644 --- a/ARCHITECTURE.md +++ b/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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 3ea7b23b..6590815e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,8 +2,8 @@ * `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`. diff --git a/docs/Glossary.md b/docs/Glossary.md index 5e7762bd..774ed88f 100644 --- a/docs/Glossary.md +++ b/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. diff --git a/singlem/archive_otu_table.py b/singlem/archive_otu_table.py index 10f2688b..c918b5fa 100644 --- a/singlem/archive_otu_table.py +++ b/singlem/archive_otu_table.py @@ -5,26 +5,32 @@ On-disk encodings ----------------- Versions 1-4 store `otus` as a list of rows, each row a list of values in -`fields` order. Version 5 stores the same rows column-wise, hoists fields whose +`fields` order. Version 6 stores the same rows column-wise, hoists fields whose value is identical in every row into `constant_fields`, and dereplicates read sequences into a shared `reads` list that the rows refer to. Those three changes only remove repetition; nothing is added or lost. -Version 5 also adds one field, `reads_with_repaired_deletions`, which records -which of an OTU's reads had an ambiguous base inserted into their window by -frameshift repair. `pipe` needs that to tell such a base apart from one already -present in a raw read, and until version 5 it only existed in memory, so the -resolution of ambiguous windows had to happen inside the same `pipe` run. It is -None whenever no read of the OTU has one, which is every OTU of a run without -frameshift repair, so constant hoisting collapses the whole column to a single -entry there. - -The gain is real but modest: on test/data/SRR8653040.json (95 OTUs, 1197 reads) -gzip -9 goes from 36.6 KB to 35.6 KB and zstd -19 from 32.1 KB to 30.4 KB. Most -of what remains is beyond the reach of any layout - the read sequences -themselves, and the per-package sha256s in the header, which are incompressible -and alone account for 13% of that gzipped file. Storing fewer or shorter reads is -the only thing that would shrink an archive substantially. +Version 6 also adds one field, `read_unaligned_sequence_variants`, recording the +modifications made to each read by frameshift repair. The unaligned sequence +stored is the read after repair, and for each read the field holds either None +(not modified) or the repairs as VCF-style variants that turn the stored +sequence back into the read as sequenced (see +frameshift_repair.repair_frameshifts_with_variants). `pipe` also needs this to +tell an ambiguous base that repair inserted apart from one already present in a +raw read, when resolving such bases is deferred to `summarise`. The field is +None for an OTU none of whose reads were repaired, which is every OTU of a run +without frameshift repair, so constant hoisting collapses the whole column to a +single entry there. + +Version 5 was the same layout with a different record of frameshift repair, and +was only ever written by development versions of SingleM. It is not read. + +The gain from the layout is real but modest: on test/data/SRR8653040.json (95 +OTUs, 1197 reads) gzip -9 goes from 36.6 KB to 35.6 KB and zstd -19 from 32.1 KB +to 30.4 KB. Most of what remains is beyond the reach of any layout - the read +sequences themselves, and the per-package sha256s in the header, which are +incompressible and alone account for 13% of that gzipped file. Storing fewer or +shorter reads is the only thing that would shrink an archive substantially. In memory the table is the same either way: `self.data` is a list of rows in `self.fields` order, so nothing downstream of `read()` needs to know which @@ -37,7 +43,7 @@ class ArchiveOtuTable: - version = 5 + version = 6 FIELDS_VERSION1 = str.split( 'gene sample sequence num_hits coverage taxonomy read_names nucleotides_aligned taxonomy_by_known?' @@ -51,18 +57,21 @@ class ArchiveOtuTable: FIELDS_VERSION4 = str.split( 'gene sample sequence num_hits coverage taxonomy read_names nucleotides_aligned taxonomy_by_known? read_unaligned_sequences equal_best_hit_taxonomies taxonomy_assignment_method' ) - FIELDS_VERSION5 = FIELDS_VERSION4 + str.split('reads_with_repaired_deletions') - FIELDS_OF_EACH_VERSION = [ - FIELDS_VERSION1, - FIELDS_VERSION2, - FIELDS_VERSION3, - FIELDS_VERSION4, - FIELDS_VERSION5, - ] - FIELDS = FIELDS_OF_EACH_VERSION[version - 1] + FIELDS_VERSION6 = FIELDS_VERSION4 + str.split('read_unaligned_sequence_variants') + # Version 5 is absent, being unsupported (see above). + FIELDS_OF_EACH_VERSION = { + 1: FIELDS_VERSION1, + 2: FIELDS_VERSION2, + 3: FIELDS_VERSION3, + 4: FIELDS_VERSION4, + 6: FIELDS_VERSION6, + } + FIELDS = FIELDS_OF_EACH_VERSION[version] # The first version written column-wise rather than row-wise. - FIRST_COLUMNAR_VERSION = 5 + FIRST_COLUMNAR_VERSION = 6 + + UNSUPPORTED_VERSION = 5 # The one field whose values are dereplicated into the shared 'reads' list, # being both by far the largest (~85% of an uncompressed archive) and the @@ -80,11 +89,7 @@ class ArchiveOtuTable: # the common case into a run of identical values, which does compress away. NEXT_NEW_READ = 0 - # Fields that exist for SingleM's own bookkeeping, and so are left out of the - # extended (TSV) OTU table. reads_with_repaired_deletions indexes into an - # OTU's read names, which is of no use to someone reading a TSV, and adding a - # column there would change that format to no-one's benefit. - FIELDS_NOT_IN_EXTENDED_OTU_TABLE = frozenset(['reads_with_repaired_deletions']) + UNALIGNED_SEQUENCE_VARIANTS_FIELD_NAME = 'read_unaligned_sequence_variants' READ_NAME_FIELD_INDEX = 6 SAMPLE_ID_FIELD_INDEX = FIELDS_VERSION4.index('sample') @@ -95,7 +100,7 @@ class ArchiveOtuTable: TAXONOMY_FIELD_INDEX = FIELDS_VERSION4.index('taxonomy') NUCLEOTIDES_ALIGNED_FIELD_INDEX = FIELDS_VERSION4.index('nucleotides_aligned') TAXONOMY_BY_KNOWN_FIELD_INDEX = FIELDS_VERSION4.index('taxonomy_by_known?') - READS_WITH_REPAIRED_DELETIONS_FIELD_INDEX = FIELDS_VERSION5.index('reads_with_repaired_deletions') + UNALIGNED_SEQUENCE_VARIANTS_FIELD_INDEX = FIELDS_VERSION6.index(UNALIGNED_SEQUENCE_VARIANTS_FIELD_NAME) def __init__(self, singlem_packages=None): self.singlem_packages = singlem_packages @@ -161,7 +166,7 @@ def _is_constant_column(values): return all(type(value) is type(first) and value == first for value in values) def _columnar_body(self): - '''The version 5 encoding of self.data: columns, minus the fields that + '''The version 6 encoding of self.data: columns, minus the fields that are constant across every row, with read sequences dereplicated.''' rows = self._rows_padded_to_fields() columns = {field: [row[i] for row in rows] for (i, field) in enumerate(self.fields)} @@ -250,7 +255,12 @@ def _rows_from_columnar_body(j): def read(input_io, min_version=None): otus = ArchiveOtuTable() j = json.load(input_io) - if j['version'] not in range(1, ArchiveOtuTable.version + 1): + if j['version'] == ArchiveOtuTable.UNSUPPORTED_VERSION: + raise Exception( + "Archive OTU tables of version {} were only written by development " + "versions of SingleM, and are not supported. Regenerate it with the " + "current version of SingleM.".format(j['version'])) + if j['version'] not in ArchiveOtuTable.FIELDS_OF_EACH_VERSION: raise Exception("Wrong OTU table version detected") otus.version = j['version'] if min_version is not None and otus.version < min_version: @@ -261,7 +271,7 @@ def read(input_io, min_version=None): otus.singlem_package_sha256s = j['singlem_package_sha256s'] otus.fields = j['fields'] - if otus.fields != ArchiveOtuTable.FIELDS_OF_EACH_VERSION[j['version'] - 1]: + if otus.fields != ArchiveOtuTable.FIELDS_OF_EACH_VERSION[j['version']]: raise Exception("Unexpected archive OTU table format detected") if otus.version >= ArchiveOtuTable.FIRST_COLUMNAR_VERSION: @@ -270,6 +280,14 @@ def read(input_io, min_version=None): otus.data = j['otus'] return otus + @staticmethod + def version_with_fields(fields): + '''The version whose field list is fields, or None if there is none.''' + for (version, version_fields) in ArchiveOtuTable.FIELDS_OF_EACH_VERSION.items(): + if version_fields == fields: + return version + return None + def sort(self): self.data.sort() @@ -308,15 +326,28 @@ def nucleotides_aligned(self): def taxonomy_by_known(self): return self.data[ArchiveOtuTable.TAXONOMY_BY_KNOWN_FIELD_INDEX] - def reads_with_repaired_deletions(self): - '''Indices into read_names() of the reads whose window carries an - ambiguous base that frameshift repair inserted, as opposed to one that - was already in the raw read. None when there are none, and when the OTU - came from an archive written before version 5, which did not record - it.''' - if len(self.data) <= ArchiveOtuTable.READS_WITH_REPAIRED_DELETIONS_FIELD_INDEX: + def read_unaligned_sequence_variants(self): + '''For each read of read_unaligned_sequences(), the modifications that + frameshift repair made to it, as a list of VCF-style [POS, REF, ALT] + variants that turn the stored sequence back into the read as sequenced, + or None if it was not modified (see + frameshift_repair.repair_frameshifts_with_variants). None rather than a + list when no read was modified, and when the OTU came from an archive + written before version 6, which did not record it.''' + if len(self.data) <= ArchiveOtuTable.UNALIGNED_SEQUENCE_VARIANTS_FIELD_INDEX: return None - return self.data[ArchiveOtuTable.READS_WITH_REPAIRED_DELETIONS_FIELD_INDEX] + return self.data[ArchiveOtuTable.UNALIGNED_SEQUENCE_VARIANTS_FIELD_INDEX] + + def read_original_unaligned_sequences(self): + '''read_unaligned_sequences(), but as the reads were before any + modification by frameshift repair.''' + from .frameshift_repair import revert_variants + sequences = self.read_unaligned_sequences() + variants = self.read_unaligned_sequence_variants() + if variants is None: + return sequences + return [sequence if read_variants is None else revert_variants(sequence, read_variants) + for (sequence, read_variants) in zip(sequences, variants)] def coverage(self): return self.data[ArchiveOtuTable.COVERAGE_FIELD_INDEX] diff --git a/singlem/diamond_spkg_searcher.py b/singlem/diamond_spkg_searcher.py index 1ac0a0ef..93397a1c 100755 --- a/singlem/diamond_spkg_searcher.py +++ b/singlem/diamond_spkg_searcher.py @@ -9,7 +9,8 @@ from subprocess import Popen, PIPE from .utils import FastaNameToSampleName -from .frameshift_repair import walk_btop, repair_frameshifts as repair_frameshifts_in_sequence +from .frameshift_repair import walk_btop, repair_frameshifts as repair_frameshifts_in_sequence, \ + repair_frameshifts_with_variants class DiamondSpkgSearcher: def __init__(self, num_threads, working_directory): @@ -137,7 +138,7 @@ def _prefilter(self, diamond_database, read_files, is_reverse_reads, performance query_sequence_lengths = {} num_frameshifts_repaired = 0 num_reads_with_frameshifts = 0 - repaired_deletion_qseqids = set() + unaligned_sequence_variants = {} do_logging = not sys.stderr.isatty() or logging.getLogger().level != logging.INFO if not do_logging: # Start animation thread @@ -229,6 +230,15 @@ def _prefilter(self, diamond_database, read_files, is_reverse_reads, performance # to the prefilter database to the fasta file. truncated_qseq = full_qseq[int(qstart_int)-1:int(qend_int)] + # The read as it is to be recorded in the archive OTU + # table. With a context window that is only the region + # around the hit, and cutting it down to that is not + # recorded as a modification of the read. + if context_window is not None: + recorded_qseq = truncated_qseq + else: + recorded_qseq = full_qseq + # Restore the reading frame of the region that is about to # be translated and aligned to the HMM. Only frameshifts # inside the region matter, and they are re-based onto it. @@ -244,24 +254,32 @@ def _prefilter(self, diamond_database, read_files, is_reverse_reads, performance truncated_qseq, frameshifts_here) num_frameshifts_repaired += len(frameshifts_here) num_reads_with_frameshifts += 1 - # Only deletions insert an ambiguous base; track - # those hits so resolve_ambiguous_windows() can - # later tell a repair-inserted 'N' apart from an - # 'N' that was already in the raw read. - if any(kind == 'deletion' for (_, kind) in frameshifts_here): - repaired_deletion_qseqids.add(unique_qseqid) + # The same repairs, but of the read as recorded, + # described so that the archive OTU table can + # hold the repaired read and still allow the + # original to be recreated. This also lets + # resolve_ambiguous_windows() tell an 'N' that + # repair inserted apart from one that was already + # in the raw read. + recorded_offset = region_start if context_window is not None else 0 + (_, variants) = repair_frameshifts_with_variants( + recorded_qseq, + [(position - recorded_offset, kind) + for (position, kind) in frameshifts + if region_start <= position < region_end]) + if len(variants) > 0: + unaligned_sequence_variants[unique_qseqid] = variants fasta_file.write(f'>{unique_qseqid}\n{truncated_qseq}\n') # Write the full read sequence to the fasta file as well # so it can be included in the archive OTU table later # on. Write without the unique_qseqid modification, to save space. + # This is the read before repair, since a read with + # hits to several markers can be repaired differently + # for each; pipe applies each hit's repairs to it. if qseqid not in seen_full_qseqs: - if context_window is not None: - qseq_full_seq = full_qseq[int(qstart_int)-1:int(qend_int)] - else: - qseq_full_seq = full_qseq - full_qseq_f.write(f'>{qseqid}\n{qseq_full_seq}\n') + full_qseq_f.write(f'>{qseqid}\n{recorded_qseq}\n') seen_full_qseqs.add(qseqid) # Stop animation and clear the line @@ -288,23 +306,23 @@ def _prefilter(self, diamond_database, read_files, is_reverse_reads, performance diamond_results.append(DiamondSearchResult( fasta_path, full_qseq_fasta_path, best_hits, query_sequence_lengths, - repaired_deletion_qseqids)) + unaligned_sequence_variants)) logging.info(f"Found {len(best_hits)} hits for {os.path.basename(file)}") return diamond_results class DiamondSearchResult: def __init__(self, query_sequence_file, full_query_sequences_file, best_hits, query_sequence_lengths, - repaired_deletion_qseqids=None): + unaligned_sequence_variants=None): self.query_sequences_file = query_sequence_file self.full_query_sequences_file = full_query_sequences_file self.best_hits = best_hits self.query_sequence_lengths = query_sequence_lengths - # unique_qseqid values (see _prefilter) of hits where frameshift repair - # inserted an ambiguous base for a deletion, i.e. whose 'N' is safe to - # resolve in resolve_ambiguous_windows(), as opposed to an 'N' that was - # already present in the raw read. - self.repaired_deletion_qseqids = repaired_deletion_qseqids if repaired_deletion_qseqids is not None else set() + # unique_qseqid (see _prefilter) of each hit whose frameshifts were + # repaired, to the repairs as VCF-style variants (see + # frameshift_repair.repair_frameshifts_with_variants) relative to the + # read as written to full_query_sequences_file. + self.unaligned_sequence_variants = unaligned_sequence_variants if unaligned_sequence_variants is not None else {} def sample_name(self): return FastaNameToSampleName().fasta_to_name(self.query_sequences_file) diff --git a/singlem/frameshift_repair.py b/singlem/frameshift_repair.py index 010ea712..cbefd94a 100644 --- a/singlem/frameshift_repair.py +++ b/singlem/frameshift_repair.py @@ -131,8 +131,7 @@ def repair_frameshifts(sequence, frameshifts, ambiguous_char=AMBIGUOUS_CHAR): frame. Deletions have ambiguous_char inserted (the deleted base is unknown); - insertions have the extra base removed. Edits are applied from the 3' end - backwards so that earlier positions stay valid. + insertions have the extra base removed. Parameters ---------- @@ -145,17 +144,155 @@ def repair_frameshifts(sequence, frameshifts, ambiguous_char=AMBIGUOUS_CHAR): ------- the repaired sequence. ''' - for (position, kind) in sorted(frameshifts, key=lambda f: -f[0]): - if position < 0 or position > len(sequence): + return repair_frameshifts_with_variants(sequence, frameshifts, ambiguous_char)[0] + + +def repair_frameshifts_with_variants(sequence, frameshifts, ambiguous_char=AMBIGUOUS_CHAR): + '''As repair_frameshifts(), but also describe the repairs as VCF-style + variants that turn the repaired sequence back into the original one. + + The repaired sequence plays the part of the VCF reference, and the original + the alternate allele, because it is the repaired sequence that is kept (e.g. + in the archive OTU table) and the original that needs recreating from it. + Each variant is [POS, REF, ALT]: POS is 1-based in the repaired sequence, + REF is the repaired sequence from POS onwards, and ALT what the original has + in its place. As in VCF, each carries the base before the edit, or the base + after it when the edit is at the very start. A deleted base repaired as 'N' + after a 'C' at POS 7 is [7, 'CN', 'C']; an extra 'T' removed after that 'C' + is [7, 'C', 'CT']. Edits close enough for those bases to overlap are merged + into one variant, so variants never overlap and are in order of POS. + + Returns + ------- + (repaired sequence, list of variants), the list being empty when no + frameshift could be applied. + ''' + # Edits as (position in the repaired sequence, bases there in the repaired + # sequence, bases there in the original), without the anchoring base. + edits = [] + repaired_parts = [] + repaired_length = 0 + consumed = 0 + for (position, kind) in sorted(frameshifts, key=lambda f: f[0]): + if position < consumed or position > len(sequence) or \ + (kind != 'deletion' and position == len(sequence)): logging.warning( "Ignoring frameshift at out-of-range position %i of a %i bp sequence", position, len(sequence)) continue + repaired_parts.append(sequence[consumed:position]) + repaired_length += position - consumed if kind == 'deletion': - sequence = sequence[:position] + ambiguous_char + sequence[position:] + edits.append((repaired_length, ambiguous_char, '')) + repaired_parts.append(ambiguous_char) + repaired_length += 1 + consumed = position + else: + edits.append((repaired_length, '', sequence[position])) + consumed = position + 1 + repaired_parts.append(sequence[consumed:]) + repaired = ''.join(repaired_parts) + + # Give each edit its anchoring base, merging those whose spans then overlap. + clusters = [] + for edit in edits: + (position, ref, _) = edit + if position > 0: + (start, end) = (position - 1, position + len(ref)) else: - sequence = sequence[:position] + sequence[position + 1:] - return sequence + (start, end) = (0, min(len(repaired), len(ref) + 1)) + if len(clusters) > 0 and start < clusters[-1][1]: + clusters[-1][1] = max(clusters[-1][1], end) + clusters[-1][2].append(edit) + else: + clusters.append([start, end, [edit]]) + + variants = [] + for (start, end, cluster_edits) in clusters: + alt_parts = [] + cursor = start + for (position, ref, alt) in cluster_edits: + alt_parts.append(repaired[cursor:position]) + alt_parts.append(alt) + cursor = position + len(ref) + alt_parts.append(repaired[cursor:end]) + variants.append([start + 1, repaired[start:end], ''.join(alt_parts)]) + return (repaired, variants) + + +def apply_variants(original, variants): + '''Recreate the repaired sequence from the original one, given the variants + repair_frameshifts_with_variants() returned for it.''' + parts = [] + cursor = 0 + shift = 0 + for (position, ref, alt) in variants: + start = position - 1 + shift + if start < cursor or original[start:start + len(alt)] != alt: + raise Exception( + "Variant {} does not match the sequence it is applied to".format( + [position, ref, alt])) + parts.append(original[cursor:start]) + parts.append(ref) + cursor = start + len(alt) + shift += len(alt) - len(ref) + parts.append(original[cursor:]) + return ''.join(parts) + + +def revert_variants(repaired, variants): + '''Recreate the original sequence from the repaired one, given the variants + repair_frameshifts_with_variants() returned for it. This is how a reader of + an archive OTU table recovers a read as it was before its frameshifts were + repaired.''' + parts = [] + cursor = 0 + for (position, ref, alt) in variants: + start = position - 1 + if start < cursor or repaired[start:start + len(ref)] != ref: + raise Exception( + "Variant {} does not match the sequence it is applied to".format( + [position, ref, alt])) + parts.append(repaired[cursor:start]) + parts.append(alt) + cursor = start + len(ref) + parts.append(repaired[cursor:]) + return ''.join(parts) + + +def original_length(repaired_length, variants): + '''Length of the original sequence, given that of the repaired one.''' + if variants is None: + return repaired_length + return repaired_length + sum(len(alt) - len(ref) for (_, ref, alt) in variants) + + +def variants_insert_ambiguous_base(variants, ambiguous_char=AMBIGUOUS_CHAR): + '''True if the variants record a repaired deletion, i.e. the repaired + sequence has an ambiguous_char that the original did not, as opposed to one + that was already in the raw read.''' + if variants is None: + return False + return any(ref.count(ambiguous_char) > alt.count(ambiguous_char) + for (_, ref, alt) in variants) + + +NO_VARIANTS_TEXT = '.' + + +def format_variants(otu_variants, num_reads): + '''Text form of an OTU's read_unaligned_sequence_variants, for a TSV cell. + + Reads are separated by spaces, as other per-read fields are, and each read + is NO_VARIANTS_TEXT if unmodified, or else its variants as POS:REF>ALT + separated by commas, e.g. "121:AN>A,150:T>TT". otu_variants may be None, + meaning that none of the num_reads reads were modified.''' + if otu_variants is None: + otu_variants = [None] * num_reads + return ' '.join( + NO_VARIANTS_TEXT if variants is None + else ','.join('%i:%s>%s' % (pos, ref, alt) for (pos, ref, alt) in variants) + for variants in otu_variants) def shift_position(position, frameshifts): @@ -335,7 +472,7 @@ def resolve_windows(ambiguous_windows, abundances, def resolve_ambiguous_windows(window_sequences, max_divergence=DEFAULT_MAX_FRAMESHIFT_REPAIR_DIVERGENCE, ambiguous_char=AMBIGUOUS_CHAR, - repaired_names=None): + only_repair_inserted=False): '''Fill in ambiguous bases in window sequences from other windows in the same set. @@ -357,12 +494,12 @@ def resolve_ambiguous_windows(window_sequences, modified in place (their aligned_sequence attribute is updated). max_divergence: int maximum number of mismatches at unambiguous positions. - repaired_names: set of str, or None - if given, only windows whose .name is in this set are eligible to be - resolved, since an ambiguous_char elsewhere is not from frameshift - repair (e.g. an 'N' already present in the raw read) and inventing a - base for it would misrepresent the read. If None, every ambiguous - window is eligible, as when no such provenance is available. + only_repair_inserted: bool + if True, only windows whose unaligned_sequence_variants record a + repaired deletion are eligible to be resolved, since an ambiguous_char + elsewhere is not from frameshift repair (e.g. an 'N' already present in + the raw read) and inventing a base for it would misrepresent the read. + If False, every ambiguous window is eligible. Returns ------- @@ -370,7 +507,8 @@ def resolve_ambiguous_windows(window_sequences, ''' ambiguous = [s for s in window_sequences if ambiguous_char in s.aligned_sequence - and (repaired_names is None or s.name in repaired_names)] + and (not only_repair_inserted + or variants_insert_ambiguous_base(s.unaligned_sequence_variants))] if len(ambiguous) == 0: return 0 diff --git a/singlem/genome_fasta_mux.py b/singlem/genome_fasta_mux.py index 90712c8d..9afb032d 100644 --- a/singlem/genome_fasta_mux.py +++ b/singlem/genome_fasta_mux.py @@ -24,7 +24,7 @@ def demux_otu_table(self, otu_table_object): archive = otu_table_object.archive(None) new_otus = [] - if ArchiveOtuTable.version != 5: + if ArchiveOtuTable.version != 6: raise ValueError("Unsupported OTU table version provided to genome demux - likely a programming bug: %d" % ArchiveOtuTable.version) for otu in archive: @@ -56,7 +56,7 @@ def demux_otu_table(self, otu_table_object): otu.equal_best_hit_taxonomies(), otu.taxonomy_assignment_method(), # Genome input does not go through the DIAMOND prefilter, so - # no window of it has a frameshift-repaired base. + # no sequence of it is frameshift-repaired. None ]) for genome_name, otus in genome_to_otus.items(): diff --git a/singlem/otu_table.py b/singlem/otu_table.py index da9dca5c..609fac2b 100644 --- a/singlem/otu_table.py +++ b/singlem/otu_table.py @@ -2,6 +2,7 @@ from .archive_otu_table import ArchiveOtuTable from .otu_table_entry import OtuTableEntry +from .frameshift_repair import format_variants class OtuTable: @@ -146,9 +147,15 @@ def write_otus_to(otu_table_entries, output_io, fields_to_print=None): if fields_to_print is None: fields_to_print = OtuTable.DEFAULT_OUTPUT_FIELDS output_io.write("\t".join(fields_to_print)+"\n") for o in otu_table_entries: - field_indices_to_print = [o.fields.index(f) for f in fields_to_print] - output_io.write("\t".join([OtuTable._to_printable(cell) for cell in [ - o.data[i] for i in field_indices_to_print]])+"\n") + output_io.write("\t".join([ + OtuTable._field_to_printable(o, f) for f in fields_to_print])+"\n") + + @staticmethod + def _field_to_printable(otu, field): + cell = otu.data[otu.fields.index(field)] + if field == ArchiveOtuTable.UNALIGNED_SEQUENCE_VARIANTS_FIELD_NAME: + return format_variants(cell, len(otu.data[otu.fields.index('read_names')])) + return OtuTable._to_printable(cell) @staticmethod def _to_printable(e): diff --git a/singlem/pipe.py b/singlem/pipe.py index 09419cbe..b9deb8cc 100644 --- a/singlem/pipe.py +++ b/singlem/pipe.py @@ -27,7 +27,7 @@ from .taxonomy import * from .otu_table_collection import StreamingOtuTableCollection, OtuTableCollection from .genome_fasta_mux import GenomeFastaMux -from .frameshift_repair import resolve_ambiguous_windows, DEFAULT_MAX_FRAMESHIFT_REPAIR_DIVERGENCE +from .frameshift_repair import resolve_ambiguous_windows, apply_variants, DEFAULT_MAX_FRAMESHIFT_REPAIR_DIVERGENCE from graftm.sequence_extractor import SequenceExtractor from graftm.greengenes_taxonomy import GreenGenesTaxonomy @@ -424,11 +424,6 @@ def return_cleanly(): if not diamond_prefilter: diamond_package_assignment = False - # Populated below when the DIAMOND prefilter runs with frameshift - # repair on; used later to restrict resolve_ambiguous_windows() to - # windows whose 'N' actually came from a repaired deletion. - repaired_deletion_qseqids = set() - #### Extract diamond_prefilter_performance_parameters from metapackage (v5 metapackages only) if diamond_prefilter: # Set the min ORF length in DIAMOND, as this saves CPU time and @@ -532,24 +527,10 @@ def finish_sra_extraction_processes(sra_extraction_processes, sra_extraction_com if any([len(r.best_hits)>0 for r in diamond_forward_search_results]): found_a_hit = True forward_read_files = list([r.query_sequences_file for r in diamond_forward_search_results]) - repaired_deletion_qseqids = set() - for r in diamond_forward_search_results: - repaired_deletion_qseqids.update(r.repaired_deletion_qseqids) if analysing_pairs: reverse_read_files = list([r.query_sequences_file for r in diamond_reverse_search_results]) if any([len(r.best_hits)>0 for r in diamond_reverse_search_results]): found_a_hit = True - for r in diamond_reverse_search_results: - repaired_deletion_qseqids.update(r.repaired_deletion_qseqids) - if input_sra_files: - # unknown_sequences are renamed by KingfisherSra to strip the - # SRA .0/.1/.2 direction suffix and the ••gene suffix (see - # KingfisherSra._split_extracted_reads); normalise these - # qseqids the same way so they still match by .name below. - sra_regex = KingfisherSra()._split_regex() - repaired_deletion_qseqids = { - (m.group(1) if (m := sra_regex.match(qseqid)) else qseqid) - for qseqid in repaired_deletion_qseqids} logging.info("Finished DIAMOND prefilter phase") if not found_a_hit: logging.info("No reads identified in any samples, stopping") @@ -651,7 +632,7 @@ def finish_sra_extraction_processes(sra_extraction_processes, sra_extraction_com num_resolved += resolve_ambiguous_windows( window_sequences, max_divergence=max_frameshift_repair_divergence, - repaired_names=repaired_deletion_qseqids) + only_repair_inserted=True) logging.info( "Resolved the ambiguous base(s) of {} window sequence(s)".format( num_resolved)) @@ -676,7 +657,6 @@ def finish_sra_extraction_processes(sra_extraction_processes, sra_extraction_com assignment_singlem_db=assignment_singlem_db, diamond_forward_qseqs = diamond_forward_qseqs if diamond_prefilter else None, diamond_reverse_qseqs = diamond_reverse_qseqs if (diamond_prefilter and analysing_pairs) else None, - repaired_deletion_read_names = repaired_deletion_qseqids, ) if genome_fasta_files: @@ -701,8 +681,6 @@ def assign_taxonomy_and_process(self, **kwargs): assignment_singlem_db = kwargs.pop('assignment_singlem_db') diamond_forward_qseqs = kwargs.pop('diamond_forward_qseqs') diamond_reverse_qseqs = kwargs.pop('diamond_reverse_qseqs') - # renew has no frameshift repair information, so it does not pass this. - repaired_deletion_read_names = kwargs.pop('repaired_deletion_read_names', set()) if len(kwargs) > 0: raise Exception("Unexpected arguments detected: %s" % kwargs) @@ -741,7 +719,6 @@ def assign_taxonomy_and_process(self, **kwargs): known_sequence_tax if known_sequence_taxonomy else None, diamond_forward_qseqs, diamond_reverse_qseqs, - repaired_deletion_read_names, # outputs otu_table_object, package_to_taxonomy_bihash, @@ -787,7 +764,6 @@ def _process_taxonomically_assigned_reads( known_sequence_tax, diamond_forward_qseqs, diamond_reverse_qseqs, - repaired_deletion_read_names, # outputs otu_table_object, package_to_taxonomy_bihash): @@ -809,12 +785,12 @@ def add_info(infos, otu_table_object, known_tax): # before printing otu table names = [name.split('••')[0] for name in info.names] - (names, unaligned_sequences, aligned_lengths, repaired_deletions) = \ + (names, unaligned_sequences, aligned_lengths, unaligned_sequence_variants) = \ self._read_columns_sorted_by_name( names, info.unaligned_sequences, info.aligned_lengths, - info.repaired_deletions) + info.unaligned_sequence_variants) to_print = [ singlem_package.graftm_package_basename(), @@ -829,10 +805,11 @@ def add_info(infos, otu_table_object, known_tax): unaligned_sequences, info.equal_best_taxonomies, info.taxonomy_assignment_method, - # None rather than an empty list when no read has one, so the - # column collapses to a single constant in the archive (see - # ArchiveOtuTable). - repaired_deletions if len(repaired_deletions) > 0 else None] + # None rather than a list of Nones when no read was + # repaired, so the column collapses to a single constant in + # the archive (see ArchiveOtuTable). + unaligned_sequence_variants + if any(v is not None for v in unaligned_sequence_variants) else None] otu_table_object.data.append(to_print) def extract_placement_parser( @@ -871,9 +848,7 @@ def process_readset( known_sequence_taxonomy, None, None, - None, - None, - repaired_deletion_read_names) + None) add_info(known_infos, otu_table_object, True) if not analysing_pairs and len(readset.unknown_sequences) == 0: @@ -884,9 +859,35 @@ def process_readset( return [] else: # if any sequences were aligned (not just already known) if diamond_forward_qseqs: + # Up to now each read's unaligned_sequence has been just the + # part of it that hit in the prefilter. The archive OTU table + # records the whole read (or with a context window, the region + # around the hit), with its frameshifts repaired. This is the + # last point at which it is known which reads are forward and + # which reverse. + # + # The names here are e.g. + # 'HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482••2524614704', + # whereas the full sequences are named without the ••marker. forward_full_qseqs = SeqReader().read_nucleotide_sequences(diamond_forward_qseqs[sample_name]) - if diamond_reverse_qseqs: - reverse_full_qseqs = SeqReader().read_nucleotide_sequences(diamond_reverse_qseqs[sample_name]) + if analysing_pairs: + reverse_full_qseqs = SeqReader().read_nucleotide_sequences(diamond_reverse_qseqs[sample_name]) + sides = [(readset[0], forward_full_qseqs), (readset[1], reverse_full_qseqs)] + else: + sides = [(readset, forward_full_qseqs)] + for (side_readset, full_qseqs) in sides: + for s in side_readset.unknown_sequences: + original = full_qseqs[s.name.split('••')[0]] + if s.unaligned_sequence_variants is None: + s.unaligned_sequence = original + else: + s.unaligned_sequence = apply_variants(original, s.unaligned_sequence_variants) + # Coverage is calculated from the length of the read as + # sequenced, not as repaired. With a context window the + # recorded read is only part of it, and the length + # already set is that of the whole read. + if not self._context_window: + s.full_nucleotide_sequence_length = len(original) if analysing_pairs: aligned_seqs = list(itertools.chain( @@ -906,12 +907,6 @@ def process_readset( equal_best_hit_hash = assignment_result.get_equal_best_hits(singlem_package, sample_name) equal_best_taxonomies = {} if analysing_pairs: - if diamond_forward_qseqs: - # Include all reads, not just those with taxonomy hits, - # since aligned_seqs may contain reads found by HMMER - # but not by DIAMOND. Reverse first, then forward overwrites. - read_name_to_fullseq = dict(reverse_full_qseqs) - read_name_to_fullseq.update(forward_full_qseqs) for (name, best_hits) in best_hit_hash[1].items(): taxonomies[name] = best_hits for (name, best_hits) in best_hit_hash[0].items(): @@ -923,8 +918,6 @@ def process_readset( # Overwrite reverse hit with the forward hit equal_best_taxonomies[name] = equal_best_hits else: - if diamond_forward_qseqs: - read_name_to_fullseq = forward_full_qseqs for (name, best_hits) in best_hit_hash.items(): taxonomies[name] = best_hits for (name, equal_best_hits) in equal_best_hit_hash.items(): @@ -979,31 +972,6 @@ def process_readset( taxonomies = known_sequence_tax else: taxonomies = {} - # We need to go through each - # unalignedalignednucleotidesequence and put in the actual - # unaligned sequence, because this is the last time we have - # access to knowing what is the forward and what is the - # reverse read sequence. - read_name_to_fullseq = None - # In [8]: readset[0].unknown_sequences[0].name - # Out[8]: 'HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482••2524614704' - # - # In [9]: forward_full_qseqs - # Out[9]: {'HWI-ST1243:156:D1K8 - if analysing_pairs: - for s in readset[0].unknown_sequences: - s.unaligned_sequence = forward_full_qseqs[s.name.split('••')[0]] - if not self._context_window: # Do not reset length if context windowing is used, because the full_qseq is truncated - s.full_nucleotide_sequence_length = None - for s in readset[1].unknown_sequences: - s.unaligned_sequence = reverse_full_qseqs[s.name.split('••')[0]] - if not self._context_window: - s.full_nucleotide_sequence_length = None - else: - for s in readset.unknown_sequences: - s.unaligned_sequence = forward_full_qseqs[s.name.split('••')[0]] - if not self._context_window: - s.full_nucleotide_sequence_length = None new_infos = list(self._seqs_to_counts_and_taxonomy( aligned_seqs, singlem_assignment_method, known_sequence_tax if known_sequence_taxonomy else {}, @@ -1013,8 +981,6 @@ def process_readset( SMAFA_NAIVE_THEN_DIAMOND_ASSIGNMENT_METHOD) else None, placement_parser if singlem_assignment_method == PPLACER_ASSIGNMENT_METHOD else None, assignment_methods, - read_name_to_fullseq if diamond_forward_qseqs else None, - repaired_deletion_read_names, )) if output_jplace: @@ -1078,9 +1044,7 @@ def _seqs_to_counts_and_taxonomy(self, sequences, per_read_taxonomies, per_read_equal_best_taxonomies, placement_parser, - taxonomy_assignment_methods, - read_name_to_fullseq, - repaired_deletion_read_names): + taxonomy_assignment_methods): '''Given an array of UnalignedAlignedNucleotideSequence objects, and taxonomic assignment-related results, yield over 'Info' objects that contain e.g. the counts of the aggregated sequences and corresponding median @@ -1099,10 +1063,6 @@ def _seqs_to_counts_and_taxonomy(self, sequences, Used only if assignment_method is PPLACER_ASSIGNMENT_METHOD. taxonomy_assignment_methods: SingleAnswerAssignmentMethodStore or MultiAnswerAssignmentMethodStore What method(s) were used to assign taxonomy to individual OTUs. - repaired_deletion_read_names: set of str - names of the reads whose window has an ambiguous base that frameshift - repair inserted, recorded per OTU so that the base can still be - resolved after this run finishes (see Summariser). ''' class CollectedInfo: def __init__(self): @@ -1115,7 +1075,7 @@ def __init__(self): self.aligned_lengths = [] self.orf_names = [] self.known_sequence_taxonomies = [] - self.repaired_deletions = [] + self.unaligned_sequence_variants = [] seq_to_collected_info = {} for s in sequences: @@ -1156,18 +1116,14 @@ def __init__(self): if per_read_equal_best_taxonomies is not None and equal_best_tax is not None: collected_info.equal_best_taxonomies.append(equal_best_tax) collected_info.names.append(s.name) - if read_name_to_fullseq: - collected_info.unaligned_sequences.append( - read_name_to_fullseq[s.name.split('••')[0]]) - else: - collected_info.unaligned_sequences.append(s.unaligned_sequence) + collected_info.unaligned_sequences.append(s.unaligned_sequence) collected_info.coverage += s.coverage_increment() collected_info.aligned_lengths.append(s.aligned_length) collected_info.orf_names.append(s.orf_name) - collected_info.repaired_deletions.append(s.name in repaired_deletion_read_names) + collected_info.unaligned_sequence_variants.append(s.unaligned_sequence_variants) class Info: - def __init__(self, seq, count, taxonomy, equal_best_taxonomies, names, unaligned_sequences, coverage, aligned_lengths, taxonomy_assignment_method, repaired_deletions): + def __init__(self, seq, count, taxonomy, equal_best_taxonomies, names, unaligned_sequences, coverage, aligned_lengths, taxonomy_assignment_method, unaligned_sequence_variants): self.seq = seq self.count = count self.taxonomy = taxonomy @@ -1177,7 +1133,7 @@ def __init__(self, seq, count, taxonomy, equal_best_taxonomies, names, unaligned self.coverage = coverage self.aligned_lengths = aligned_lengths self.taxonomy_assignment_method = taxonomy_assignment_method - self.repaired_deletions = repaired_deletions + self.unaligned_sequence_variants = unaligned_sequence_variants for seq, collected_info in seq_to_collected_info.items(): # All seqs in OTU have the same assignment method, except in rare @@ -1227,7 +1183,7 @@ def __init__(self, seq, count, taxonomy, equal_best_taxonomies, names, unaligned collected_info.coverage, collected_info.aligned_lengths, otu_taxonomy_assignment_method, - collected_info.repaired_deletions) + collected_info.unaligned_sequence_variants) def _median_taxonomy(self, taxonomies): levels_to_counts = [] @@ -1906,7 +1862,7 @@ def _diamond_assign_taxonomy_paired_output_directory( "read1" if is_forward else "read2") @staticmethod - def _read_columns_sorted_by_name(names, unaligned_sequences, aligned_lengths, repaired_deletions): + def _read_columns_sorted_by_name(names, unaligned_sequences, aligned_lengths, unaligned_sequence_variants): '''Put an OTU's per-read data in read name order. These lists are parallel - index i describes the same read in each - and @@ -1917,18 +1873,12 @@ def _read_columns_sorted_by_name(names, unaligned_sequences, aligned_lengths, re towards nucleotides_aligned without appearing in the window itself, so two reads can share a window and still differ there. - Returns the four lists in name order, except that the last becomes the - indices of the reads whose window has an ambiguous base inserted by - frameshift repair, which is how the OTU table records it. + Returns the four lists in name order. ''' reads = sorted( - zip(names, unaligned_sequences, aligned_lengths, repaired_deletions), + zip(names, unaligned_sequences, aligned_lengths, unaligned_sequence_variants), key=lambda read: read[0]) - return ( - [read[0] for read in reads], - [read[1] for read in reads], - [read[2] for read in reads], - [i for (i, read) in enumerate(reads) if read[3]]) + return tuple([read[i] for read in reads] for i in range(4)) def _remove_single_sequence_duplicates(self, readset): '''In extremely rare circumstances, a single read can have >1 OTU diff --git a/singlem/pipe_sequence_extractor.py b/singlem/pipe_sequence_extractor.py index fe5043d6..ca69f49f 100644 --- a/singlem/pipe_sequence_extractor.py +++ b/singlem/pipe_sequence_extractor.py @@ -346,6 +346,13 @@ def _extract_reads_by_diamond_for_package_and_sample(prefilter_result, spkg, logging.debug("[PID: {}] Found {} window sequences for spkg {}".format(pid, len(window_seqs),spkg.base_directory())) + + # Attached to the window sequences here, while their names still match the + # prefilter's, so that the record of each read's frameshift repairs travels + # with it through any renaming (e.g. of SRA reads) to the archive OTU table. + for window_seq in window_seqs: + window_seq.unaligned_sequence_variants = \ + prefilter_result.unaligned_sequence_variants.get(window_seq.name) return ExtractedReadSet( sample_name, spkg, sequences, [], window_seqs diff --git a/singlem/renew.py b/singlem/renew.py index a91fa968..d71c0651 100644 --- a/singlem/renew.py +++ b/singlem/renew.py @@ -9,6 +9,7 @@ from .archive_otu_table import ArchiveOtuTable from .pipe_sequence_extractor import ExtractedReads, ExtractedReadSet from .sequence_classes import Sequence, UnalignedAlignedNucleotideSequence +from .frameshift_repair import original_length from .otu_table_collection import StreamingOtuTableCollection class Renew: @@ -139,11 +140,20 @@ def process_otu_batch(state): nucleotides_aligned = otu.data[nucleotides_aligned_field] if len(read_names) != len(nucleotides_aligned) or len(seqs) != len(nucleotides_aligned): raise Exception("Unexpected format of otu found: {}".format(otu)) - for (name, seq, num_aligned) in zip(read_names, seqs, nucleotides_aligned): + # Carried through so that the output archive still records the + # frameshift repairs made by the pipe run that the input came from. + all_variants = otu.read_unaligned_sequence_variants() + if all_variants is None: + all_variants = [None] * len(read_names) + for (name, seq, num_aligned, variants) in zip(read_names, seqs, nucleotides_aligned, all_variants): state.current_sequences.append(Sequence(name, seq)) state.current_unaligned_aligned_nuc_seqs.append( UnalignedAlignedNucleotideSequence( - name, None, otu.sequence, seq, num_aligned)) + name, None, otu.sequence, seq, num_aligned, + # As in pipe, coverage is from the length of the read as + # sequenced rather than as repaired. + full_nucleotide_sequence_length=None if variants is None else original_length(len(seq), variants), + unaligned_sequence_variants=variants)) if ignore_missing_singlem_packages: logging.info("Ignored {} OTUs that were not found in the specified singlem packages".format(num_ignored_otus)) diff --git a/singlem/sequence_classes.py b/singlem/sequence_classes.py index 91dcc623..657bdc81 100644 --- a/singlem/sequence_classes.py +++ b/singlem/sequence_classes.py @@ -48,7 +48,8 @@ class UnalignedAlignedNucleotideSequence: ''' - def __init__(self, name, orf_name, aligned_sequence, unaligned_sequence, aligned_length, full_nucleotide_sequence_length=None): + def __init__(self, name, orf_name, aligned_sequence, unaligned_sequence, aligned_length, full_nucleotide_sequence_length=None, + unaligned_sequence_variants=None): ''' Parameters --------- @@ -63,6 +64,13 @@ def __init__(self, name, orf_name, aligned_sequence, unaligned_sequence, aligned aligned_length: the number of nucleotides used in the alignment, including columns that were removed as not aligned + unaligned_sequence_variants: list of [POS, REF, ALT], or None + the frameshift repairs made to the read, as VCF-style variants that + turn the repaired read back into the original one (see + frameshift_repair.repair_frameshifts_with_variants). They are + relative to the whole read as recorded in the archive OTU table, + which is not necessarily unaligned_sequence until pipe puts that + read there. None if the read was not repaired. ''' self.name = name self.orf_name = orf_name @@ -70,6 +78,7 @@ def __init__(self, name, orf_name, aligned_sequence, unaligned_sequence, aligned self.unaligned_sequence = unaligned_sequence self.aligned_length = aligned_length self.full_nucleotide_sequence_length = full_nucleotide_sequence_length + self.unaligned_sequence_variants = unaligned_sequence_variants def coverage_increment(self): '''Given the alignment came from a read of length diff --git a/singlem/summariser.py b/singlem/summariser.py index 3d1e5883..75cb78cf 100644 --- a/singlem/summariser.py +++ b/singlem/summariser.py @@ -14,9 +14,10 @@ from .taxonomy import QUERY_BASED_ASSIGNMENT_METHOD, DIAMOND_ASSIGNMENT_METHOD, NO_ASSIGNMENT_METHOD from .condense import CondensedCommunityProfile, CondensedCommunityProfileCamiIiiGtdbWriter from .frameshift_repair import ( - resolve_windows, AMBIGUOUS_CHAR, DEFAULT_MAX_FRAMESHIFT_REPAIR_DIVERGENCE) + resolve_windows, variants_insert_ambiguous_base, AMBIGUOUS_CHAR, + DEFAULT_MAX_FRAMESHIFT_REPAIR_DIVERGENCE) -REPAIRED_DELETIONS_FIELD = 'reads_with_repaired_deletions' +VARIANTS_FIELD = ArchiveOtuTable.UNALIGNED_SEQUENCE_VARIANTS_FIELD_NAME class Summariser: @staticmethod @@ -166,9 +167,7 @@ def write_otu_table(**kwargs): if output_extras: OtuTable.write_otus_to( table_collection, output_table_io, - fields_to_print=[ - field for field in table_collection.example_field_names() - if field not in ArchiveOtuTable.FIELDS_NOT_IN_EXTENDED_OTU_TABLE]) + fields_to_print=table_collection.example_field_names()) else: OtuTable.write_otus_to(table_collection, output_table_io) @@ -354,8 +353,7 @@ def _resolve_ambiguous_windows_in_dataframe(df, max_divergence): genes = df.get_column('gene').to_list() sequences = df.get_column('sequence').to_list() num_hits = df.get_column('num_hits').to_list() - read_names = df.get_column('read_names').to_list() - repaired = df.get_column(REPAIRED_DELETIONS_FIELD).to_list() + variants = df.get_column(VARIANTS_FIELD).to_list() # Abundances are pooled across samples within a gene, as in pipe, so that # a low coverage sample can still draw on a donor. @@ -368,9 +366,12 @@ def _resolve_ambiguous_windows_in_dataframe(df, max_divergence): resolvable = [] num_partial = 0 for (i, sequence) in enumerate(sequences): - if AMBIGUOUS_CHAR not in sequence or repaired[i] is None or len(repaired[i]) == 0: + if AMBIGUOUS_CHAR not in sequence or variants[i] is None: continue - if len(repaired[i]) != len(read_names[i]): + repaired = [variants_insert_ambiguous_base(v) for v in variants[i]] + if not any(repaired): + continue + if not all(repaired): # Only some of this OTU's reads have a repair-inserted base, so # resolving it would mean splitting the OTU in two, and the # per-read coverage needed to divide up its coverage exactly is @@ -394,14 +395,9 @@ def _resolve_ambiguous_windows_in_dataframe(df, max_divergence): resolved = resolutions[sequences[i]] if resolved is not None: sequences[i] = resolved - # No ambiguous base is left to resolve, so the OTU no longer - # needs marking as having one. - repaired[i] = None num_resolved += 1 - df = df.with_columns( - pl.Series('sequence', sequences), - pl.Series(REPAIRED_DELETIONS_FIELD, repaired, dtype=pl.Object)) + df = df.with_columns(pl.Series('sequence', sequences)) return (df, num_resolved, num_partial) @staticmethod @@ -435,7 +431,7 @@ def archive_schema(fields): "nucleotides_aligned", "read_unaligned_sequences", "equal_best_hit_taxonomies", - "reads_with_repaired_deletions", + VARIANTS_FIELD, }: schema[field] = pl.Object elif field == "num_hits": @@ -462,6 +458,15 @@ def rows_of_width(otus, width): return [row if len(row) == width else list(row) + [None] * (width - len(row)) for row in otus] + def archive_dataframe(rows, fields): + '''Build a column at a time, since polars infers a type for the + values of an Object column when building row-wise, and fails on the + mixed integers and strings of read_unaligned_sequence_variants.''' + schema = archive_schema(fields) + return pl.DataFrame([ + pl.Series(field, [row[i] for row in rows], dtype=schema[field]) + for (i, field) in enumerate(fields)]) + def read_archive_table(df, f, prev_ar): logging.debug("Reading archive table {} into RAM ..".format(a)) ar = ArchiveOtuTable.read(f) @@ -470,9 +475,7 @@ def read_archive_table(df, f, prev_ar): # alignment_hmm_sha256s = ar.alignment_hmm_sha256s # singlem_package_sha256s = ar.singlem_package_sha256s fields = ar.fields - df = pl.DataFrame( - ar.data, schema=archive_schema(fields), orient="row" - ) + df = archive_dataframe(ar.data, fields) else: fields = wider_fields(prev_ar.fields, ar.fields) if fields is None: @@ -487,15 +490,12 @@ def read_archive_table(df, f, prev_ar): for field in fields[len(df.columns):]: df = df.with_columns( pl.Series(field, [None] * len(df), dtype=schema[field])) - df2 = pl.DataFrame( - rows_of_width(ar.data, len(fields)), schema=archive_schema(fields), - orient="row" - ) + df2 = archive_dataframe(rows_of_width(ar.data, len(fields)), fields) df = pl.concat([df, df2], how="vertical") # The output is written at whichever version has these fields, so that # combining an old archive with a new one gives a new one. ar.fields = fields - ar.version = ArchiveOtuTable.FIELDS_OF_EACH_VERSION.index(fields) + 1 + ar.version = ArchiveOtuTable.version_with_fields(fields) return df, ar for a in archive_otu_tables: @@ -536,11 +536,11 @@ def read_archive_table(df, f, prev_ar): # taxonomy - which is the taxonomy the whole dataset would have given it, # since pipe resolves windows before assigning taxonomy. if resolve_ambiguous_windows: - if REPAIRED_DELETIONS_FIELD not in ar.fields: + if VARIANTS_FIELD not in ar.fields: raise Exception( "Cannot resolve ambiguous windows in an archive OTU table of version " - "{}, which does not record which reads have a frameshift-repaired " - "base. Regenerate it with a current version of pipe.".format(ar.version)) + "{}, which does not record which reads were frameshift-repaired. " + "Regenerate it with a current version of pipe.".format(ar.version)) logging.info("Resolving ambiguous bases in frameshift-repaired windows ..") (df, num_resolved, num_partial) = Summariser._resolve_ambiguous_windows_in_dataframe( df, max_frameshift_repair_divergence) @@ -593,19 +593,18 @@ def read_archive_table(df, f, prev_ar): tax_assignment_method, ] - if REPAIRED_DELETIONS_FIELD in ar.fields: - # The indices point into each OTU's own read names, and those are - # concatenated above in group order, so each is shifted by the - # number of reads the earlier OTUs of the group contributed. - merged_repaired = [] - offset = 0 - for (names, indices) in zip( + if VARIANTS_FIELD in ar.fields: + # Concatenated in the same order as the read names above. An OTU + # with no repaired read holds None rather than a list, so it + # contributes one None per read. + merged_variants = [] + for (names, otu_variants) in zip( grouped.get_column("read_names").to_list(), - grouped.get_column(REPAIRED_DELETIONS_FIELD).to_list()): - if indices is not None: - merged_repaired.extend(offset + i for i in indices) - offset += len(names) - collapsed_row.append(merged_repaired if len(merged_repaired) > 0 else None) + grouped.get_column(VARIANTS_FIELD).to_list()): + merged_variants.extend( + [None] * len(names) if otu_variants is None else otu_variants) + collapsed_row.append( + merged_variants if any(v is not None for v in merged_variants) else None) collapsed_rows.append(collapsed_row) diff --git a/test/test_appraise.py b/test/test_appraise.py index 95f85186..21951ace 100644 --- a/test/test_appraise.py +++ b/test/test_appraise.py @@ -55,11 +55,11 @@ def _shas_of(self, alignment_shas_json, package_shas_json): def assertArchiveEqual(self, expected_version, expected_shas, expected_otu_jsons, observed_json): '''Assert an archive OTU table holds the given OTUs. The rows are compared - as parsed by ArchiveOtuTable rather than as text, since from version 5 on + as parsed by ArchiveOtuTable rather than as text, since from version 6 on the same rows have more than one valid on-disk spelling.''' observed = ArchiveOtuTable.read(StringIO(observed_json)) self.assertEqual(expected_version, observed.version) - self.assertEqual(ArchiveOtuTable.FIELDS_OF_EACH_VERSION[expected_version - 1], observed.fields) + self.assertEqual(ArchiveOtuTable.FIELDS_OF_EACH_VERSION[expected_version], observed.fields) (alignment_shas, package_shas) = expected_shas self.assertEqual(alignment_shas, observed.alignment_hmm_sha256s) self.assertEqual(package_shas, observed.singlem_package_sha256s) @@ -710,8 +710,8 @@ def test_print_appraisal_archive_input_output(self): ) expected_shas = self._shas_of(alignment_shas, package_shas) - self.assertArchiveEqual(5, expected_shas, [binned_otu], found_otu_table_io.getvalue()) - self.assertArchiveEqual(5, expected_shas, [unbinned_otu], not_found_otu_table_io.getvalue()) + self.assertArchiveEqual(6, expected_shas, [binned_otu], found_otu_table_io.getvalue()) + self.assertArchiveEqual(6, expected_shas, [unbinned_otu], not_found_otu_table_io.getvalue()) def test_print_appraisal_archive_input_output_cli(self): with in_tempdir(): @@ -741,8 +741,8 @@ def test_print_appraisal_archive_input_output_cli(self): observed_not_found = "".join(f.readlines()) expected_shas = self._shas_of(alignment_shas, package_shas) - self.assertArchiveEqual(5, expected_shas, [binned_otu], observed_found) - self.assertArchiveEqual(5, expected_shas, [unbinned_otu], observed_not_found) + self.assertArchiveEqual(6, expected_shas, [binned_otu], observed_found) + self.assertArchiveEqual(6, expected_shas, [unbinned_otu], observed_not_found) def test_contamination(self): metagenome_otu_table = [ diff --git a/test/test_archive_otu_table.py b/test/test_archive_otu_table.py index 93b795c2..766284d7 100644 --- a/test/test_archive_otu_table.py +++ b/test/test_archive_otu_table.py @@ -34,14 +34,14 @@ def otu(gene='4.11.22seqs', sample='sample1', sequence='ACGT', num_hits=1, coverage=1.0, taxonomy='Root', read_names=None, nucleotides_aligned=None, taxonomy_by_known=False, read_unaligned_sequences=None, equal_best_hit_taxonomies=None, - taxonomy_assignment_method='diamond', reads_with_repaired_deletions=None): + taxonomy_assignment_method='diamond', read_unaligned_sequence_variants=None): return [gene, sample, sequence, num_hits, coverage, taxonomy, read_names if read_names is not None else ['read1'], nucleotides_aligned if nucleotides_aligned is not None else [60], taxonomy_by_known, read_unaligned_sequences if read_unaligned_sequences is not None else ['AAAACCCC'], equal_best_hit_taxonomies, taxonomy_assignment_method, - reads_with_repaired_deletions] + read_unaligned_sequence_variants] class Tests(unittest.TestCase): @@ -69,7 +69,7 @@ def assertRoundTrips(self, data): def test_writes_the_current_version(self): j = json.loads(self.write([otu()])) - self.assertEqual(5, j['version']) + self.assertEqual(6, j['version']) self.assertEqual(ArchiveOtuTable.FIELDS, j['fields']) def test_round_trip_single_otu(self): @@ -175,13 +175,34 @@ def test_reads_version_4(self): self.assertEqual(95, len(observed.data)) self.assertEqual(ArchiveOtuTable.FIELDS_VERSION4, observed.fields) - def test_version_4_and_version_5_hold_the_same_otus(self): + def test_version_4_and_version_6_hold_the_same_otus(self): with open(os.path.join(path_to_data, 'SRR8653040.json')) as f: version_4 = ArchiveOtuTable.read(f) - version_5 = ArchiveOtuTable.read(StringIO(self.write(version_4.data))) - # Version 5 adds reads_with_repaired_deletions, which a version 4 archive - # does not record, so its OTUs gain a None for it. - self.assertEqual([row + [None] for row in version_4.data], version_5.data) + version_6 = ArchiveOtuTable.read(StringIO(self.write(version_4.data))) + # Version 6 adds read_unaligned_sequence_variants, which a version 4 + # archive does not record, so its OTUs gain a None for it. + self.assertEqual([row + [None] for row in version_4.data], version_6.data) + + def test_version_5_rejected(self): + j = json.loads(self.write([otu()])) + j['version'] = 5 + with self.assertRaisesRegex(Exception, 'version 5'): + ArchiveOtuTable.read(StringIO(json.dumps(j))) + + def test_round_trip_variants(self): + variants = [None, [[3, 'GN', 'G'], [8, 'C', 'CT']]] + observed = self.round_trip([otu( + read_names=['read1', 'read2'], nucleotides_aligned=[60, 60], + read_unaligned_sequences=['AAAACCCC', 'ACGNTTTCAA'], + read_unaligned_sequence_variants=variants)]) + entry = list(observed)[0] + self.assertEqual(variants, entry.read_unaligned_sequence_variants()) + self.assertEqual(['AAAACCCC', 'ACGTTTCTAA'], entry.read_original_unaligned_sequences()) + + def test_original_unaligned_sequences_without_variants(self): + entry = list(self.round_trip([otu()]))[0] + self.assertEqual(None, entry.read_unaligned_sequence_variants()) + self.assertEqual(['AAAACCCC'], entry.read_original_unaligned_sequences()) def test_version_4_tables_are_still_written_row_wise(self): table = ArchiveOtuTable() @@ -199,8 +220,8 @@ def test_version_4_tables_are_still_written_row_wise(self): def test_min_version(self): with self.assertRaises(InsufficientArchiveOtuTableVersionException): - ArchiveOtuTable.read(StringIO(self.write([otu()])), min_version=6) - ArchiveOtuTable.read(StringIO(self.write([otu()])), min_version=5) + ArchiveOtuTable.read(StringIO(self.write([otu()])), min_version=7) + ArchiveOtuTable.read(StringIO(self.write([otu()])), min_version=6) def test_unknown_version_rejected(self): j = json.loads(self.write([otu()])) diff --git a/test/test_diamond_spkg_searcher.py b/test/test_diamond_spkg_searcher.py index b4b73bf1..eaffbce1 100644 --- a/test/test_diamond_spkg_searcher.py +++ b/test/test_diamond_spkg_searcher.py @@ -51,7 +51,7 @@ def test_malformed_diamond_output_suggests_validating_input_reads(self): repair_frameshifts=False, ) - def _prefilter_on_one_hit(self, tempdir, line, repair_frameshifts): + def _prefilter_on_one_hit(self, tempdir, line, repair_frameshifts, context_window=None): searcher = DiamondSpkgSearcher(num_threads=1, working_directory=tempdir) with patch( 'singlem.diamond_spkg_searcher.Popen', @@ -64,7 +64,7 @@ def _prefilter_on_one_hit(self, tempdir, line, repair_frameshifts): performance_parameters='', sample_names=['sample.fna'], min_orf_length=0, - context_window=None, + context_window=context_window, repair_frameshifts=repair_frameshifts, ) @@ -83,15 +83,19 @@ def test_frameshift_repair_inserts_ambiguous_base_for_a_deletion(self): results = self._prefilter_on_one_hit(tempdir, line, True) with open(results[0].query_sequences_file) as f: output = f.read() + with open(results[0].full_query_sequences_file) as f: + full = f.read() # Only the aligned region (nt 1-26) is written, with an 'N' inserted at # the frameshift, so it is one base longer than the region. repaired = output.split('\n')[1] self.assertEqual('ATGAAANCCCGGGTTTAAACCCGGGTT', repaired) self.assertEqual(27, len(repaired)) self.assertEqual(0, len(repaired) % 3) - # The hit's deletion should be recorded so resolve_ambiguous_windows() - # can later tell this 'N' apart from one in the raw read. - self.assertEqual({'read1••gene'}, results[0].repaired_deletion_qseqids) + # The repair is recorded against the whole read, which is what goes + # into the archive OTU table, and the read is written there unrepaired. + self.assertEqual( + {'read1••gene': [[6, 'AN', 'A']]}, results[0].unaligned_sequence_variants) + self.assertEqual('>read1\n{}\n'.format(sequence), full) def test_frameshift_repair_removes_extra_base_for_an_insertion(self): sequence = 'ATGAAACCCGGGTTTAAACCCGGGTTTAAA' @@ -107,9 +111,9 @@ def test_frameshift_repair_removes_extra_base_for_an_insertion(self): self.assertEqual('ATGAAACCGGGTTTAAACCCGGGTTTA', repaired) self.assertEqual(27, len(repaired)) self.assertEqual(0, len(repaired) % 3) - # An insertion does not leave an ambiguous base, so nothing needs - # tracking for resolve_ambiguous_windows(). - self.assertEqual(set(), results[0].repaired_deletion_qseqids) + # The removed base is recorded so that the read can be recreated. + self.assertEqual( + {'read1••gene': [[6, 'A', 'AC']]}, results[0].unaligned_sequence_variants) def test_frameshift_repair_inserts_ambiguous_base_for_a_deletion_on_reverse_strand(self): # Same deletion as test_frameshift_repair_inserts_ambiguous_base_for_a_deletion, @@ -125,6 +129,20 @@ def test_frameshift_repair_inserts_ambiguous_base_for_a_deletion_on_reverse_stra self.assertEqual(27, len(repaired)) self.assertEqual(0, len(repaired) % 3) + def test_frameshift_repair_with_context_window_is_relative_to_the_window(self): + # With a context window only the region around the hit (here nt 4-30) is + # recorded, and cutting the read down to it is not itself a variant. + sequence = 'GGGATGAAACCCGGGTTTAAACCCGGGTTTAAA' + line = 'read1\t{}\tgene~ref1\t7\t29\t2/-6\n'.format(sequence) + with tempfile.TemporaryDirectory() as tempdir: + results = self._prefilter_on_one_hit(tempdir, line, True, context_window=3) + with open(results[0].full_query_sequences_file) as f: + full = f.read() + name = 'read1:4-32,33••gene' + self.assertEqual('>read1:4-32,33\n{}\n'.format(sequence[3:32]), full) + # The read is missing a base after nt 12, i.e. nt 9 of the region. + self.assertEqual({name: [[9, 'CN', 'C']]}, results[0].unaligned_sequence_variants) + def test_no_frameshift_repair_when_not_requested(self): # Without --repair-frameshifts, DIAMOND is not asked for BTOP and the # sequence is written through unchanged. diff --git a/test/test_frameshift_repair.py b/test/test_frameshift_repair.py index f0d56c07..0333dcf9 100644 --- a/test/test_frameshift_repair.py +++ b/test/test_frameshift_repair.py @@ -28,8 +28,9 @@ sys.path = [os.path.join(os.path.dirname(os.path.realpath(__file__)), '..')] + sys.path from singlem.frameshift_repair import ( - parse_btop, walk_btop, repair_frameshifts, shift_position, - resolve_ambiguous_windows) + parse_btop, walk_btop, repair_frameshifts, repair_frameshifts_with_variants, + apply_variants, revert_variants, original_length, variants_insert_ambiguous_base, + format_variants, shift_position, resolve_ambiguous_windows) class ParseBtopTests(unittest.TestCase): @@ -138,6 +139,87 @@ def test_restores_reading_frame(self): self.assertEqual(0, len(repaired) % 3) +class RepairFrameshiftsWithVariantsTests(unittest.TestCase): + def assertRecreates(self, sequence, frameshifts): + (repaired, variants) = repair_frameshifts_with_variants(sequence, frameshifts) + self.assertEqual(repair_frameshifts(sequence, frameshifts), repaired) + self.assertEqual(sequence, revert_variants(repaired, variants)) + self.assertEqual(repaired, apply_variants(sequence, variants)) + self.assertEqual(len(sequence), original_length(len(repaired), variants)) + return (repaired, variants) + + def test_no_frameshifts(self): + self.assertEqual(('ACGT', []), self.assertRecreates('ACGT', [])) + + def test_deletion(self): + # The base before the edit anchors it, as in VCF. + self.assertEqual( + ('ACNGTACGT', [[2, 'CN', 'C']]), + self.assertRecreates('ACGTACGT', [(2, 'deletion')])) + + def test_insertion(self): + self.assertEqual( + ('ACTACGT', [[2, 'C', 'CG']]), + self.assertRecreates('ACGTACGT', [(2, 'insertion')])) + + def test_edit_at_start_anchored_by_following_base(self): + self.assertEqual( + ('NACGT', [[1, 'NA', 'A']]), + self.assertRecreates('ACGT', [(0, 'deletion')])) + self.assertEqual( + ('CGT', [[1, 'C', 'AC']]), + self.assertRecreates('ACGT', [(0, 'insertion')])) + + def test_deletion_at_end(self): + self.assertEqual( + ('ACGTN', [[4, 'TN', 'T']]), + self.assertRecreates('ACGT', [(4, 'deletion')])) + + def test_several_frameshifts(self): + self.assertEqual( + ('ACNGTCGT', [[2, 'CN', 'C'], [5, 'T', 'TA']]), + self.assertRecreates('ACGTACGT', [(2, 'deletion'), (4, 'insertion')])) + + def test_adjacent_frameshifts_merged(self): + # Anchored separately, these would both start at the C, which VCF does + # not allow. + (repaired, variants) = self.assertRecreates( + 'ACGTACGT', [(2, 'insertion'), (3, 'deletion')]) + self.assertEqual('ACNTACGT', repaired) + self.assertEqual([[2, 'CN', 'CG']], variants) + + def test_random_frameshifts_recreate_the_original(self): + import random + rng = random.Random(42) + for _ in range(2000): + sequence = ''.join(rng.choice('ACGTN') for _ in range(rng.randint(1, 12))) + frameshifts = list({ + rng.randint(0, len(sequence)): rng.choice(['deletion', 'insertion']) + for _ in range(rng.randint(0, 4))}.items()) + (_, variants) = self.assertRecreates(sequence, frameshifts) + for (first, second) in zip(variants, variants[1:]): + self.assertLessEqual(first[0] + len(first[1]), second[0]) + + def test_mismatched_variant_rejected(self): + with self.assertRaises(Exception): + revert_variants('ACGT', [[2, 'GN', 'G']]) + with self.assertRaises(Exception): + apply_variants('ACGT', [[2, 'GN', 'G']]) + + def test_variants_insert_ambiguous_base(self): + self.assertTrue(variants_insert_ambiguous_base([[2, 'CN', 'C']])) + self.assertFalse(variants_insert_ambiguous_base([[2, 'C', 'CG']])) + self.assertFalse(variants_insert_ambiguous_base(None)) + # An 'N' that was in the read all along is not one repair inserted. + self.assertFalse(variants_insert_ambiguous_base([[2, 'N', 'NA']])) + + def test_format_variants(self): + self.assertEqual('. .', format_variants(None, 2)) + self.assertEqual( + '121:AN>A,150:T>TT .', + format_variants([[[121, 'AN', 'A'], [150, 'T', 'TT']], None], 2)) + + class ShiftPositionTests(unittest.TestCase): def test_no_frameshifts(self): self.assertEqual(10, shift_position(10, [])) @@ -156,9 +238,10 @@ class FakeWindowSequence: '''Stands in for UnalignedAlignedNucleotideSequence, which needs more setup than these tests require.''' - def __init__(self, aligned_sequence, name=None): + def __init__(self, aligned_sequence, name=None, unaligned_sequence_variants=None): self.aligned_sequence = aligned_sequence self.name = name + self.unaligned_sequence_variants = unaligned_sequence_variants class ResolveAmbiguousWindowsTests(unittest.TestCase): @@ -173,23 +256,22 @@ def test_resolved_from_identical_window(self): self.assertEqual(1, resolve_ambiguous_windows(sequences)) self.assertEqual('ACGTACGT', ambiguous.aligned_sequence) - def test_repaired_names_restricts_resolution(self): + def test_only_repair_inserted_restricts_resolution(self): # An 'N' that came from the raw read rather than frameshift repair must # not be filled in, even if a perfect donor is available. raw_n = FakeWindowSequence('ACGNACGT', name='raw_read') - repaired = FakeWindowSequence('ACGNACGT', name='repaired_read') + repaired = FakeWindowSequence( + 'ACGNACGT', name='repaired_read', unaligned_sequence_variants=[[10, 'GN', 'G']]) sequences = [raw_n, repaired, FakeWindowSequence('ACGTACGT', name='donor')] self.assertEqual( - 1, resolve_ambiguous_windows(sequences, repaired_names={'repaired_read'})) + 1, resolve_ambiguous_windows(sequences, only_repair_inserted=True)) self.assertEqual('ACGNACGT', raw_n.aligned_sequence) self.assertEqual('ACGTACGT', repaired.aligned_sequence) - def test_repaired_names_none_resolves_everything(self): - # Without provenance information, every ambiguous window is eligible, - # matching the previous behaviour. + def test_not_only_repair_inserted_resolves_everything(self): ambiguous = FakeWindowSequence('ACGNACGT', name='unknown') sequences = [ambiguous, FakeWindowSequence('ACGTACGT', name='donor')] - self.assertEqual(1, resolve_ambiguous_windows(sequences, repaired_names=None)) + self.assertEqual(1, resolve_ambiguous_windows(sequences, only_repair_inserted=False)) self.assertEqual('ACGTACGT', ambiguous.aligned_sequence) def test_most_abundant_donor_wins(self): diff --git a/test/test_pipe.py b/test/test_pipe.py index d9ff270b..ce4ad5e6 100644 --- a/test/test_pipe.py +++ b/test/test_pipe.py @@ -21,6 +21,7 @@ # If not, see . #======================================================================= +import csv import unittest import os.path import tempfile @@ -45,7 +46,7 @@ class Tests(unittest.TestCase): headers = str.split('gene sample sequence num_hits coverage taxonomy') - headers_with_extras = headers + str.split('read_names nucleotides_aligned taxonomy_by_known? read_unaligned_sequences equal_best_hit_taxonomies taxonomy_assignment_method') + headers_with_extras = headers + str.split('read_names nucleotides_aligned taxonomy_by_known? read_unaligned_sequences equal_best_hit_taxonomies taxonomy_assignment_method read_unaligned_sequence_variants') maxDiff = None two_packages = '%s %s' % ( os.path.join(path_to_data, '4.11.22seqs.gpkg.spkg'), @@ -756,22 +757,22 @@ def test_read_columns_sorted_by_name_keeps_reads_together(self): # sorting. The aligned lengths here differ between reads of the one OTU, # which happens when an insert column inside the window's span is a gap in # some reads and not others. - (names, unaligned_sequences, aligned_lengths, repaired_deletions) = \ + (names, unaligned_sequences, aligned_lengths, variants) = \ SearchPipe._read_columns_sorted_by_name( ['read_c', 'read_a', 'read_b'], ['CCC', 'AAA', 'BBB'], [63, 60, 57], - [False, True, False]) + [None, [[1, 'AN', 'A']], None]) self.assertEqual(['read_a', 'read_b', 'read_c'], names) self.assertEqual(['AAA', 'BBB', 'CCC'], unaligned_sequences) self.assertEqual([60, 57, 63], aligned_lengths) - # read_a is the one with a frameshift-repaired base, and sorts first. - self.assertEqual([0], repaired_deletions) + # read_a is the one that was frameshift-repaired, and sorts first. + self.assertEqual([[[1, 'AN', 'A']], None, None], variants) def test_read_columns_sorted_by_name_keeps_duplicate_names_in_order(self): (names, unaligned_sequences, aligned_lengths, _) = \ SearchPipe._read_columns_sorted_by_name( - ['read_a', 'read_a'], ['first', 'second'], [60, 57], [False, False]) + ['read_a', 'read_a'], ['first', 'second'], [60, 57], [None, None]) self.assertEqual(['read_a', 'read_a'], names) self.assertEqual(['first', 'second'], unaligned_sequences) self.assertEqual([60, 57], aligned_lengths) @@ -791,11 +792,11 @@ def test_archive_otu_groopm_compatibility(self): path_to_script, n.name, os.path.join(path_to_data,'4.11.22seqs.gpkg.spkg')) # Read the JSON directly rather than through ArchiveOtuTable, since - # that is what an external consumer does. From version 5 the OTUs are + # that is what an external consumer does. From version 6 the OTUs are # stored column-wise, with the fields that are the same in every OTU # hoisted into constant_fields. j = json.loads(extern.run(cmd)) - self.assertEqual(5, j['version']) + self.assertEqual(6, j['version']) def column(field): if field in j['constant_fields']: return [j['constant_fields'][field]] * j['num_otus'] @@ -807,6 +808,102 @@ def column(field): for name in read_names] ) + def _frameshifted_reads(self): + '''Reads with a base deleted from, or inserted into, the window region + of one that hits 4.11.22seqs, which DIAMOND reports as frameshifts.''' + with open(os.path.join(path_to_data, 'nanopore', 'fastas', 'one_hit.fasta')) as f: + read = f.read().split('\n')[1] + return { + 'deletion': read[:120] + read[121:], + 'insertion': read[:80] + 'T' + read[80:], + } + + def _pipe_frameshifted_reads(self, tempdir, reads, extra_args=''): + fasta = os.path.join(tempdir, 'reads.fa') + with open(fasta, 'w') as f: + for (name, seq) in reads.items(): + f.write('>{}\n{}\n'.format(name, seq)) + archive = os.path.join(tempdir, 'archive.json') + extern.run('{} pipe --sequences {} --archive-otu-table {} --singlem-packages {} ' + '--assignment-method diamond {}'.format( + path_to_script, fasta, archive, + os.path.join(path_to_data, '4.11.22seqs.gpkg.spkg'), extra_args)) + with open(archive) as f: + return (archive, ArchiveOtuTable.read(f)) + + def _read_data(self, archive): + return { + name: (stored, variants, original) + for otu in archive + for (name, stored, variants, original) in zip( + otu.read_names(), otu.read_unaligned_sequences(), + otu.read_unaligned_sequence_variants() or [None] * len(otu.read_names()), + otu.read_original_unaligned_sequences())} + + def test_archive_records_frameshift_repairs_as_variants(self): + reads = self._frameshifted_reads() + with tempfile.TemporaryDirectory() as tempdir: + (_, archive) = self._pipe_frameshifted_reads(tempdir, reads) + observed = self._read_data(archive) + self.assertEqual(sorted(reads), sorted(observed)) + # The stored read is the repaired one, and the variants turn it back into + # the read as sequenced. + (stored, variants, original) = observed['deletion'] + self.assertEqual([[121, 'AN', 'A']], variants) + self.assertEqual(len(reads['deletion']) + 1, len(stored)) + self.assertEqual(reads['deletion'], original) + (stored, variants, original) = observed['insertion'] + self.assertEqual([[79, 'T', 'TT']], variants) + self.assertEqual(len(reads['insertion']) - 1, len(stored)) + self.assertEqual(reads['insertion'], original) + + def test_otu_table_extras_include_frameshift_variants(self): + reads = self._frameshifted_reads() + with tempfile.TemporaryDirectory() as tempdir: + otu_table = os.path.join(tempdir, 'otu_table.tsv') + self._pipe_frameshifted_reads( + tempdir, reads, '--otu-table {} --output-extras'.format(otu_table)) + with open(otu_table) as f: + rows = list(csv.DictReader(f, delimiter='\t')) + observed = {} + for row in rows: + observed.update(zip( + row['read_names'].split(' '), + row['read_unaligned_sequence_variants'].split(' '))) + self.assertEqual({'deletion': '121:AN>A', 'insertion': '79:T>TT'}, observed) + + def test_archive_frameshift_variants_with_context_window(self): + # Only the region around the hit is stored, and that it is shorter than + # the read is not recorded as a variant, so the variants recreate just + # that region of the read. + reads = self._frameshifted_reads() + with tempfile.TemporaryDirectory() as tempdir: + (_, archive) = self._pipe_frameshifted_reads( + tempdir, reads, '--context-window 20') + observed = {name.split(':')[0]: data for (name, data) in self._read_data(archive).items()} + self.assertEqual(sorted(reads), sorted(observed)) + for (name, (stored, variants, original)) in observed.items(): + self.assertEqual(1, len(variants)) + self.assertLess(len(original), len(reads[name])) + self.assertIn(original, reads[name]) + self.assertNotEqual(stored, original) + + def test_renew_keeps_frameshift_variants(self): + reads = self._frameshifted_reads() + with tempfile.TemporaryDirectory() as tempdir: + (archive_path, archive) = self._pipe_frameshifted_reads(tempdir, reads) + renewed_path = os.path.join(tempdir, 'renewed.json') + extern.run('{} renew --input-archive-otu-table {} --singlem-packages {} ' + '--assignment-method diamond --archive-otu-table {}'.format( + path_to_script, archive_path, + os.path.join(path_to_data, '4.11.22seqs.gpkg.spkg'), renewed_path)) + with open(renewed_path) as f: + renewed = ArchiveOtuTable.read(f) + self.assertEqual(self._read_data(archive), self._read_data(renewed)) + self.assertEqual( + sorted(round(otu.coverage, 6) for otu in archive), + sorted(round(otu.coverage, 6) for otu in renewed)) + def test_protein_package_non60_length(self): expected = [ "\t".join(self.headers), @@ -899,7 +996,7 @@ def test_full_seqs_reported_in_extras(self): ''' expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 1 2.11 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 60 False CTAAATAAATAAATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA [\'2524614704\'] diamond', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 1 2.11 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 60 False CTAAATAAATAAATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA [\'2524614704\'] diamond .', ''] with tempfile.NamedTemporaryFile(mode='w',suffix='.fa') as n: n.write(inseqs) @@ -917,7 +1014,7 @@ def test_paired_reads_one_read_each(self): # Reads should be merged expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 4.88 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT [\'2524614704\'] [\'2524614704\'] diamond', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 4.88 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT [\'2524614704\'] [\'2524614704\'] diamond . .', ''] inseqs = '''>HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 1:N:0:TAAGGCGACTAAGCCT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA @@ -949,7 +1046,7 @@ def test_paired_reads_two_reads_each(self): # Reads should be merged expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 3 7.32 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 seq3 60 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA [\'2524614704\'] [\'2524614704\'] [\'2524614704\'] diamond', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 3 7.32 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 seq3 60 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA [\'2524614704\'] [\'2524614704\'] [\'2524614704\'] diamond . . .', ''] inseqs = '''>HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 1:N:0:TAAGGCGACTAAGCCT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA @@ -986,7 +1083,7 @@ def test_paired_reads_one_read_each_diamond(self): # Reads should be merged expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 4.88 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT [\'2524614704\'] [\'2524614704\'] diamond', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 4.88 Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT [\'2524614704\'] [\'2524614704\'] diamond . .', ''] inseqs = '''>HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 1:N:0:TAAGGCGACTAAGCCT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA @@ -1019,7 +1116,7 @@ def test_paired_reads_one_read_each_diamond_example(self): # Reads should be merged expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 4.88 2524614704 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT None diamond_example', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 4.88 2524614704 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAAT None diamond_example . .', ''] inseqs = '''>HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 1:N:0:TAAGGCGACTAAGCCT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA @@ -1087,7 +1184,7 @@ def test_two_orfs_in_same_read(self): # more frequently) expected = [ "\t".join(self.headers_with_extras), - 'S1.12.ribosomal_protein_S12_S23 CGTGGTGTCTGCACCCGGGTGTACACCACCACCNCGAAGAAGCCGAACTCGGCGTTGCGC 1 1.64 Root; d__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Propionibacteriales; f__Nocardioidaceae; g__Aeromicrobium A00178:38:H5NYYDSXX:2:1552:32524:1517 60 False CGGGATGTAGGCAGTGACCTCCACGCCTGAGGAGAGCCGGACGCGTGCGACCTTGCGCAACGCCGAGTTCGGCTTCTTCGGGTGGTGGTGTACACCCGGGTGCAGACACCACGGCGCTGGGGCGAACCCTTGAGCGCAGGGGTGTTGGTCT [\'GCA_001426755.1\', \'GCA_900167475.1\'] diamond', + 'S1.12.ribosomal_protein_S12_S23 CGTGGTGTCTGCACCCGGGTGTACACCACCACCNCGAAGAAGCCGAACTCGGCGTTGCGC 1 1.64 Root; d__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Propionibacteriales; f__Nocardioidaceae; g__Aeromicrobium A00178:38:H5NYYDSXX:2:1552:32524:1517 60 False CGGGATGTAGGCAGTGACCTCCACGCCTGAGGAGAGCCGGACGCGTGCGACCTTGCGCAACGCCGAGTTCGGCTTCTTCGNGGTGGTGGTGTACACCCGGGTGCAGACACCACGGCGCTGGGGCGAACCCTTGAGCGCAGGGGTGTTGGTCT [\'GCA_001426755.1\', \'GCA_900167475.1\'] diamond 80:GN>G', ''] inseqs = '''>A00178:38:H5NYYDSXX:2:1552:32524:1517 1:N:0:CAACGGA+ATCCGTT CGGGATGTAGGCAGTGACCTCCACGCCTGAGGAGAGCCGGACGCGTGCGACCTTGCGCAACGCCGAGTTCGGCTTCTTCGGGTGGTGGTGTACACCCGGGTGCAGACACCACGGCGCTGGGGCGAACCCTTGAGCGCAGGGGTGTTGGTCT @@ -1346,8 +1443,8 @@ def test_annoy_only_assignment_paired(self): @unittest.skipIf(not TEST_ANNOY, "annoy not installed") def test_annoy_only_assignment_output_extras_single(self): - expected = ['gene sample sequence num_hits coverage taxonomy read_names nucleotides_aligned taxonomy_by_known? read_unaligned_sequences equal_best_hit_taxonomies taxonomy_assignment_method', - '4.11.22seqs 4.11.22seqs.gpkg.spkg_inseqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 1 2.44 Root; part_of_sdb HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA Root; part_of_sdb; Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 Root; part_of_sdb; novel_domain singlem_query_based'] + expected = ['gene sample sequence num_hits coverage taxonomy read_names nucleotides_aligned taxonomy_by_known? read_unaligned_sequences equal_best_hit_taxonomies taxonomy_assignment_method read_unaligned_sequence_variants', + '4.11.22seqs 4.11.22seqs.gpkg.spkg_inseqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 1 2.44 Root; part_of_sdb HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 60 False ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA Root; part_of_sdb; Root; d__Bacteria; p__Firmicutes; c__Clostridia; o__Clostridiales; f__Lachnospiraceae; g__[Lachnospiraceae_bacterium_NK4A179]; s__Lachnospiraceae_bacterium_NK4A179 Root; part_of_sdb; novel_domain singlem_query_based .'] cmd = '{} pipe --sequences {} --otu-table /dev/stdout --singlem-packages {} --assignment-singlem-db {} --assignment-method annoy --output-extras'.format( path_to_script, os.path.join(path_to_data, '4.11.22seqs.gpkg.spkg_inseqs.fna'), @@ -1367,7 +1464,7 @@ def test_annoy_only_assignment_output_archive_single(self): os.path.join(path_to_data, '4.11.22seqs.gpkg.spkg_inseqs_and_inseqs2.manually_different_species.otu_table.csv.sdb'), ) # Compared as OTUs rather than as JSON, since the same OTUs have more than - # one valid on-disk spelling from archive version 5 on. + # one valid on-disk spelling from archive version 6 on. observed = ArchiveOtuTable.read(StringIO(extern.run(cmd))) self.assertEqual(ArchiveOtuTable.version, observed.version) self.assertEqual(["4b0bf5b3d7fd2ca16e54eed59d3a07eab388f70f7078ac096bf415f1c04731d9"], @@ -1532,7 +1629,7 @@ def test_read_name_parsing(self): all_ascii_chars = ''.join(chr(i) for i in list(range(33, 127)) + list(range(161, 255))) expected = [ "\t".join(self.headers_with_extras), - f'S1.12.ribosomal_protein_S12_S23 CGTGGTGTCTGCACCCGGGTGTACACCACCACCNCGAAGAAGCCGAACTCGGCGTTGCGC 1 1.64 Root; d__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Propionibacteriales; f__Nocardioidaceae; g__Aeromicrobium {all_ascii_chars} 60 False CGGGATGTAGGCAGTGACCTCCACGCCTGAGGAGAGCCGGACGCGTGCGACCTTGCGCAACGCCGAGTTCGGCTTCTTCGGGTGGTGGTGTACACCCGGGTGCAGACACCACGGCGCTGGGGCGAACCCTTGAGCGCAGGGGTGTTGGTCT [\'GCA_001426755.1\', \'GCA_900167475.1\'] diamond', + f'S1.12.ribosomal_protein_S12_S23 CGTGGTGTCTGCACCCGGGTGTACACCACCACCNCGAAGAAGCCGAACTCGGCGTTGCGC 1 1.64 Root; d__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Propionibacteriales; f__Nocardioidaceae; g__Aeromicrobium {all_ascii_chars} 60 False CGGGATGTAGGCAGTGACCTCCACGCCTGAGGAGAGCCGGACGCGTGCGACCTTGCGCAACGCCGAGTTCGGCTTCTTCGNGGTGGTGGTGTACACCCGGGTGCAGACACCACGGCGCTGGGGCGAACCCTTGAGCGCAGGGGTGTTGGTCT [\'GCA_001426755.1\', \'GCA_900167475.1\'] diamond 80:GN>G', ''] inseqs = inseqs = ( f'>{all_ascii_chars}\n' @@ -1555,7 +1652,7 @@ def test_no_assign_taxonomy_keeps_full_qseq(self): truncated, and the wrong coverage value to be reported.''' expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 2.54 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 True ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA None no_assign_taxonomy', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 2.54 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 seq2 60 60 True ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA TTCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA None no_assign_taxonomy . .', ''] inseqs = '''>HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 1:N:0:TAAGGCGACTAAGCCT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA @@ -1588,7 +1685,7 @@ def test_context_window_no_assign_taxonomy(self): truncated, and the wrong coverage value to be reported.''' expected = [ "\t".join(self.headers_with_extras), - '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 2.54 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482:1-99,253 seq2:2-104,304 60 60 True ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGA TCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAATAAAA None no_assign_taxonomy', + '4.11.22seqs TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA 2 2.54 HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482:1-99,253 seq2:2-104,304 60 60 True ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGA TCAGCTGCACGACGTACCATAGTGTTTTTGTATACTTTATACTCAACACCAGCTTCACGTAATTGTGAACGTAAGTCAGTAACTTCAGCTACTGTTAATAAAA None no_assign_taxonomy . .', ''] inseqs = '''>HWI-ST1243:156:D1K83ACXX:7:1106:18671:79482 1:N:0:TAAGGCGACTAAGCCT ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA diff --git a/test/test_summariser.py b/test/test_summariser.py index 9dca2ab9..9fd430bb 100644 --- a/test/test_summariser.py +++ b/test/test_summariser.py @@ -361,16 +361,32 @@ def _chunk_archive(self, path, otus, version=None): table.singlem_package_sha256s = ['spkg_sha'] if version is not None: table.version = version - table.fields = ArchiveOtuTable.FIELDS_OF_EACH_VERSION[version - 1] + table.fields = ArchiveOtuTable.FIELDS_OF_EACH_VERSION[version] otus = [otu[:len(table.fields)] for otu in otus] table.data = otus with open(path, 'w') as f: table.write_to(f) + # How frameshift repair records a deletion it inserted an 'N' for. + REPAIRED_DELETION = [[2, 'AN', 'A']] + def _otu(self, sequence, read_names, repaired=None): + '''repaired gives the indices of the reads that had a deletion repaired.''' + if repaired is None: + variants = None + else: + variants = [self.REPAIRED_DELETION if i in repaired else None + for i in range(len(read_names))] return ['gene1', 'chunked_sample', sequence, len(read_names), float(len(read_names)), 'Root; d__Bacteria', read_names, [60] * len(read_names), False, - ['AAAA'] * len(read_names), None, 'no_assign_taxonomy', repaired] + ['AANAA'] * len(read_names), None, 'no_assign_taxonomy', variants] + + def _repaired_reads(self, otu): + '''Names of the reads of otu recorded as having a deletion repaired.''' + variants = otu.read_unaligned_sequence_variants() + if variants is None: + return [] + return sorted(name for (name, v) in zip(otu.read_names(), variants) if v is not None) def _collapse(self, otus_per_chunk, extra_args=''): with in_tempdir(): @@ -399,8 +415,9 @@ def test_resolve_ambiguous_windows_across_chunks(self): self.assertEqual(self.DONOR_WINDOW, otu.sequence) self.assertEqual(4, otu.count) self.assertEqual(['r1', 'r2', 'r3', 'r4'], sorted(otu.read_names())) - # Nothing ambiguous is left, so nothing is left marked as repaired. - self.assertEqual(None, otu.reads_with_repaired_deletions()) + # The window is resolved, but the read itself is still as repaired, so + # its record of the repair stays. + self.assertEqual(['r1'], self._repaired_reads(otu)) def test_chunks_combine_unresolved_without_the_flag(self): combined = self._collapse( @@ -429,7 +446,7 @@ def test_ambiguous_window_with_no_donor_is_left_alone(self): otu = list(combined)[0] self.assertEqual(self.AMBIGUOUS_WINDOW, otu.sequence) # Still marked, so a later combination that does have a donor can resolve it. - self.assertEqual([0], otu.reads_with_repaired_deletions()) + self.assertEqual(['r1'], self._repaired_reads(otu)) def test_partially_repaired_otu_is_left_alone(self): # Two reads share the window but only one got its ambiguous base from a @@ -442,19 +459,23 @@ def test_partially_repaired_otu_is_left_alone(self): sorted([self.AMBIGUOUS_WINDOW, self.DONOR_WINDOW]), sorted(otu.sequence for otu in combined)) - def test_repaired_deletion_indices_are_reindexed_when_otus_merge(self): - # Two chunks each hold part of the same unresolvable OTU; the indices are - # into each chunk's own read names, so they shift when the reads join up. + def test_variants_stay_with_their_reads_when_otus_merge(self): + # Two chunks each hold part of the same unresolvable OTU, one of them + # with no repaired reads at all, which is recorded as a single None. combined = self._collapse( - [[self._otu(self.AMBIGUOUS_WINDOW, ['r1', 'r2'], repaired=[0, 1])], - [self._otu(self.AMBIGUOUS_WINDOW, ['r3'], repaired=[0])]]) + [[self._otu(self.AMBIGUOUS_WINDOW, ['r1', 'r2'], repaired=[1])], + [self._otu(self.AMBIGUOUS_WINDOW, ['r3'])]]) self.assertEqual(1, len(combined.data)) otu = list(combined)[0] self.assertEqual(3, len(otu.read_names())) - self.assertEqual([0, 1, 2], otu.reads_with_repaired_deletions()) + self.assertEqual(3, len(otu.read_unaligned_sequence_variants())) + self.assertEqual(['r2'], self._repaired_reads(otu)) + self.assertEqual( + {'r1': 'AANAA', 'r2': 'AAAA', 'r3': 'AANAA'}, + dict(zip(otu.read_names(), otu.read_original_unaligned_sequences()))) def test_collapse_archives_of_different_versions(self): - # An archive made before reads_with_repaired_deletions existed can still + # An archive made before read_unaligned_sequence_variants existed can still # be combined with one made after, since the fields only ever grow. with in_tempdir(): self._chunk_archive('old.json', [self._otu('AAAA', ['r1'])], version=4) @@ -470,8 +491,8 @@ def test_collapse_archives_of_different_versions(self): self.assertEqual(ArchiveOtuTable.version, combined.version) by_sequence = {otu.sequence: otu for otu in combined} self.assertEqual(['AAAA', 'CCCC'], sorted(by_sequence)) - self.assertEqual(None, by_sequence['AAAA'].reads_with_repaired_deletions()) - self.assertEqual([0], by_sequence['CCCC'].reads_with_repaired_deletions()) + self.assertEqual(None, by_sequence['AAAA'].read_unaligned_sequence_variants()) + self.assertEqual(['r2'], self._repaired_reads(by_sequence['CCCC'])) def test_collapse_archives_of_the_same_old_version(self): # Combining two old archives still gives that old version out. From 23bc193e7855ec9bf94aed5bd022824b1bd6ccab Mon Sep 17 00:00:00 2001 From: Ben Woodcroft Date: Wed, 28 Oct 2026 12:51:04 +1000 Subject: [PATCH 2/4] Address PR review: chunked weebill, resolved window taxonomy, weebill-only samples * pipe: reject --read-chunk-number when weebill would run, since weebill would profile the whole of the reads rather than the chunk. * summarise --resolve-ambiguous-windows: a window resolved to a donor takes the donor's taxonomy, equal best hits and assignment method, as the whole dataset would have given it. * condense: profile named weebill samples that have no marker OTUs from their weebill coverage, rather than dropping them. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019URUjTJM4Se2CmgGzPmPs5 --- CHANGELOG.md | 4 +++- singlem/condense.py | 29 +++++++++++++++++++++++++---- singlem/main.py | 2 +- singlem/pipe.py | 6 ++++++ singlem/summariser.py | 29 ++++++++++++++++++++++++----- test/test_condense.py | 40 ++++++++++++++++++++++++++++++++++++++++ test/test_pipe.py | 20 ++++++++++++++++++++ test/test_summariser.py | 24 ++++++++++++++++++++++++ 8 files changed, 143 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6590815e..a4f5a89b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,7 +6,9 @@ * 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. diff --git a/singlem/condense.py b/singlem/condense.py index 5d2c8947..cd5e25de 100644 --- a/singlem/condense.py +++ b/singlem/condense.py @@ -151,16 +151,37 @@ def _condense_to_otu_table(self, metapackage, **kwargs): if target_domains[domain] in [1, 2]: raise Exception("Number of markers for all domains must either be >= 3 or equal to 0. Only {} markers for domain '{}' found".format(target_domains[domain], domain)) - for sample, sample_otus in input_otu_table.each_sample_otus(generate_archive_otu_table=True): - - logging.debug("Processing sample {} ..".format(sample)) + def condense_sample(sample, sample_otus): apply_diamond_expectation_maximisation = True weebill_hits = self._weebill_hits_for_sample(sample, weebill_sample_to_hits) if weebill_profile is not None else None - yield self._condense_a_sample(sample, sample_otus, markers, target_domains, trim_percent, min_taxon_coverage, + return self._condense_a_sample(sample, sample_otus, markers, target_domains, trim_percent, min_taxon_coverage, True, apply_diamond_expectation_maximisation, metapackage, output_after_em_otu_table, viral_mode, weebill_hits, alpha, joint, joint_l1_penalty, joint_absence_weight, joint_min_markers, joint_adaptive_weebill_weight, joint_pin_weebill_species, joint_novel_budget) + condensed_samples = set() + for sample, sample_otus in input_otu_table.each_sample_otus(generate_archive_otu_table=True): + logging.debug("Processing sample {} ..".format(sample)) + condensed_samples.add(sample) + yield condense_sample(sample, sample_otus) + + # A sample in which weebill detected species but no marker window was + # extracted, as can happen at low depth, is absent from the OTU table, but + # its weebill species still belong in its profile. + for sample in self._weebill_only_samples(weebill_sample_to_hits, condensed_samples): + logging.info("Sample {} has no marker OTUs, so is profiled from its weebill coverage alone".format(sample)) + empty_otus = ArchiveOtuTable() + empty_otus.fields = ArchiveOtuTable.FIELDS + yield condense_sample(sample, empty_otus) + + @staticmethod + def _weebill_only_samples(weebill_sample_to_hits, condensed_samples): + '''Named weebill samples with hits that are not among condensed_samples. A + profile without a sample column (the None key) applies to every SingleM + sample, so it names no sample of its own.''' + return sorted(sample for (sample, hits) in weebill_sample_to_hits.items() + if sample is not None and len(hits) > 0 and sample not in condensed_samples) + def _domain_coverage_estimates(self, sample_otus, target_domains, trim_percent): '''Per domain, the total genome-equivalent coverage of everything in that domain, estimated the way the standard condense estimates any taxon's: the trimmed mean diff --git a/singlem/main.py b/singlem/main.py index 2ea398f4..7dc128a3 100755 --- a/singlem/main.py +++ b/singlem/main.py @@ -143,7 +143,7 @@ def add_less_common_pipe_arguments(argument_group, extra_args=False): argument_group.add_argument('--read-chunk-number', type=int, metavar='chunk_number', - help='Process only this specific chunk number (1-based index). Requires --read-chunk-size.') + help='Process only this specific chunk number (1-based index). Requires --read-chunk-size. Incompatible with running weebill, so requires --no-weebill when the metapackage bundles a weebill database and a taxonomic profile or weebill sketch is requested.') argument_group.add_argument('--output-jplace', metavar='filename', help='Output a jplace format file for each singlem package to a file starting with this string, each with one entry per OTU. Requires \'%s\' as the --assignment_method [default: unused]' % pipe.PPLACER_ASSIGNMENT_METHOD) argument_group.add_argument('--singlem-packages', nargs='+', help='SingleM packages to use [default: use the set from the default metapackage]') argument_group.add_argument('--assignment-singlem-db', '--assignment_singlem_db', help='Use this SingleM DB when assigning taxonomy [default: not set, use the default]') diff --git a/singlem/pipe.py b/singlem/pipe.py index b9deb8cc..f1fb9f76 100644 --- a/singlem/pipe.py +++ b/singlem/pipe.py @@ -98,6 +98,12 @@ def run(self, **kwargs): # --input-weebill-sketch' produce one without needing the raw reads again. sketching_weebill = output_weebill_sketch is not None and weebill_available and not running_weebill if (running_weebill or sketching_weebill) and weebill_forward_reads and not weebill_unsupported_input: + # weebill would profile the whole of the reads rather than just the + # chunk the markers were searched in. + if kwargs.get('read_chunk_number') is not None: + raise Exception("weebill cannot currently be run on a single chunk of reads, so " + "--read-chunk-number is incompatible with a weebill profile or sketch. Pass --no-weebill " + "to profile from the markers alone.") # Check weebill is installed before the marker search rather than after # it, which is an hour or more of work to throw away. if shutil.which('weebill') is None: diff --git a/singlem/summariser.py b/singlem/summariser.py index 75cb78cf..28167370 100644 --- a/singlem/summariser.py +++ b/singlem/summariser.py @@ -18,6 +18,8 @@ DEFAULT_MAX_FRAMESHIFT_REPAIR_DIVERGENCE) VARIANTS_FIELD = ArchiveOtuTable.UNALIGNED_SEQUENCE_VARIANTS_FIELD_NAME +# The fields describing an OTU's taxonomic assignment, which depend on its window. +ASSIGNMENT_FIELDS = ['taxonomy', 'equal_best_hit_taxonomies', 'taxonomy_assignment_method'] class Summariser: @staticmethod @@ -358,9 +360,14 @@ def _resolve_ambiguous_windows_in_dataframe(df, max_divergence): # Abundances are pooled across samples within a gene, as in pipe, so that # a low coverage sample can still draw on a donor. abundances_by_gene = {} - for (gene, sequence, hits) in zip(genes, sequences, num_hits): + # The row of each donor window with the most hits, whose taxonomic + # assignment a window resolving to it takes, as the grouping below does. + donor_rows = {} + for (i, (gene, sequence, hits)) in enumerate(zip(genes, sequences, num_hits)): if AMBIGUOUS_CHAR not in sequence: abundances_by_gene.setdefault(gene, Counter())[sequence] += hits + if (gene, sequence) not in donor_rows or hits > num_hits[donor_rows[(gene, sequence)]]: + donor_rows[(gene, sequence)] = i ambiguous_by_gene = {} resolvable = [] @@ -384,6 +391,12 @@ def _resolve_ambiguous_windows_in_dataframe(df, max_divergence): resolvable.append(i) ambiguous_by_gene.setdefault(genes[i], set()).add(sequence) + # The ambiguous window was assigned a taxonomy of its own, but the whole + # dataset would have resolved it before assigning taxonomy, and so given + # it the donor's assignment. + assignment_fields = [field for field in ASSIGNMENT_FIELDS if field in df.columns] + assignments = {field: df.get_column(field).to_list() for field in assignment_fields} + num_resolved = 0 for (gene, windows) in ambiguous_by_gene.items(): resolutions = resolve_windows( @@ -395,9 +408,15 @@ def _resolve_ambiguous_windows_in_dataframe(df, max_divergence): resolved = resolutions[sequences[i]] if resolved is not None: sequences[i] = resolved + donor = donor_rows[(gene, resolved)] + for values in assignments.values(): + values[i] = values[donor] num_resolved += 1 - df = df.with_columns(pl.Series('sequence', sequences)) + df = df.with_columns( + [pl.Series('sequence', sequences)] + + [pl.Series(field, values, dtype=df.schema[field]) + for (field, values) in assignments.items()]) return (df, num_resolved, num_partial) @staticmethod @@ -532,9 +551,9 @@ def read_archive_table(df, f, prev_ar): raise Exception("Multiple taxonomy_by_known found: {}".format(', '.join(map(str, taxonomy_by_known)))) # Done before the OTUs are grouped below, so that a window which resolves - # to one that is already present is grouped with it, taking that OTU's - # taxonomy - which is the taxonomy the whole dataset would have given it, - # since pipe resolves windows before assigning taxonomy. + # to one that is already present is grouped with it. It also takes that + # OTU's taxonomic assignment, which is what the whole dataset would have + # given it, since pipe resolves windows before assigning taxonomy. if resolve_ambiguous_windows: if VARIANTS_FIELD not in ar.fields: raise Exception( diff --git a/test/test_condense.py b/test/test_condense.py index 470dd52d..1927f385 100644 --- a/test/test_condense.py +++ b/test/test_condense.py @@ -365,6 +365,46 @@ def test_weebill_hits_for_sample_no_sample_column_used_for_every_sample(self): self.assertEqual(hits, condenser._weebill_hits_for_sample('sample1', weebill_sample_to_hits)) self.assertEqual(hits, condenser._weebill_hits_for_sample('sample2', weebill_sample_to_hits)) + def test_weebill_only_samples(self): + hits = {_canonical_species_key(self.S_ECOLI): WeebillHit(self.S_ECOLI, 9.0)} + self.assertEqual(['sample2', 'sample3'], Condenser._weebill_only_samples( + {'sample3': hits, 'sample1': hits, 'sample2': hits, 'empty': {}}, {'sample1'})) + # Hits without a sample name are not a sample of their own. + self.assertEqual([], Condenser._weebill_only_samples({None: hits}, set())) + + def test_sample_detected_only_by_weebill_is_profiled(self): + # sample1 has weebill hits but no marker OTUs, so it is not in the OTU + # table, but its profile should still have its weebill species. + import tempfile + from singlem.metapackage import Metapackage + from singlem.otu_table_collection import StreamingOtuTableCollection + from singlem.condense import DEFAULT_JOINT_CONDENSE_ARGUMENTS + species = 'Root; d__Archaea; p__X; c__X; o__X; f__X; g__X; s__X' + with tempfile.TemporaryDirectory() as td: + extern.run("{} metapackage --singlem-packages {} --no-nucleotide-sdb --no-taxon-genome-lengths " + "--metapackage {}".format( + path_to_script, + os.path.join(path_to_data, '..', '4.11.22seqs.v3_archaea_targetted.gpkg.spkg'), + os.path.join(td, 'a.smpkg'))) + metapackage = Metapackage.acquire(os.path.join(td, 'a.smpkg')) + for joint in (False, True): + otus = ArchiveOtuTable(metapackage.singlem_packages) + otus.data = [['4.11.22seqs', 'other', 'A' * 60, 1, 2.0, species, ['r1'], [60], False, + ['A' * 60], [species], QUERY_BASED_ASSIGNMENT_METHOD, None]] + collection = StreamingOtuTableCollection() + collection.add_archive_otu_table_object(otus) + profile = os.path.join(td, 'profile.tsv') + Condenser().condense( + input_streaming_otu_table=collection, output_otu_table=profile, krona=None, + metapackage=metapackage, + weebill_profile=os.path.join(path_to_data, 'small_weebill_profile.tsv'), + **(DEFAULT_JOINT_CONDENSE_ARGUMENTS if joint else {})) + with open(profile) as f: + rows = [line.rstrip('\n').split('\t') for line in f][1:] + self.assertIn( + ['sample1', '9.0', 'Root; ' + '; '.join(self.S_ECOLI.split(';'))], rows, + 'joint' if joint else 'standard') + def test_weebill_profile_is_unknown_corrected(self): # The column name is the only record of whether weebill was run with -u, and # condense keys alpha off it: True_cov is already on SingleM's coverage diff --git a/test/test_pipe.py b/test/test_pipe.py index ce4ad5e6..65229244 100644 --- a/test/test_pipe.py +++ b/test/test_pipe.py @@ -1356,6 +1356,26 @@ def test_read_chunk_size_forward_gzip(self): list([line.split("\t") for line in expected]), extern.run(cmd).replace('chunked_reads','')) + def test_read_chunk_number_incompatible_with_weebill(self): + # weebill would see the whole of the reads, not just the chunk. + with tempfile.TemporaryDirectory() as td: + metapackage = os.path.join(td, 'a.smpkg') + extern.run("{} metapackage --singlem-packages {} --no-nucleotide-sdb --no-taxon-genome-lengths " + "--weebill-db {} --weebill-c 100 --metapackage {}".format( + path_to_script, + os.path.join(path_to_data, '4.11.22seqs.v3_archaea_targetted.gpkg.spkg'), + os.path.join(path_to_data, 'dummy.syl2db'), + metapackage)) + cmd = "{} pipe --sequences {} --otu-table /dev/null --metapackage {} --read-chunk-size 2 " \ + "--read-chunk-number 1 --output-weebill-sketch {}".format( + path_to_script, + os.path.join(path_to_data, '1_pipe', 'small.fa'), + metapackage, + os.path.join(td, 'sketch')) + with self.assertRaises(extern.ExternCalledProcessError) as context: + extern.run(cmd) + self.assertIn('--read-chunk-number is incompatible with a weebill', str(context.exception)) + def test_read_chunk_size_paired(self): target_seq_aln = 'TTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTA' target_seq = 'ATTAACAGTAGCTGAAGTTACTGACTTACGTTCACAATTACGTGAAGCTGGTGTTGAGTATAAAGTATACAAAAACACTATGGTACGTCGTGCAGCTGAA' diff --git a/test/test_summariser.py b/test/test_summariser.py index 9fd430bb..4680b407 100644 --- a/test/test_summariser.py +++ b/test/test_summariser.py @@ -419,6 +419,30 @@ def test_resolve_ambiguous_windows_across_chunks(self): # its record of the repair stays. self.assertEqual(['r1'], self._repaired_reads(otu)) + def test_resolved_window_takes_the_donors_assignment(self): + # The ambiguous window has more reads than its donor, so without care the + # combined OTU would keep the taxonomy assigned to the ambiguous window, + # rather than the one the whole dataset would have given its reads. + def assigned(otu, taxonomy, equal_best_hit_taxonomies): + otu[5] = taxonomy + otu[10] = equal_best_hit_taxonomies + otu[11] = 'diamond' + return otu + genus = 'Root; d__Bacteria; g__A' + species = 'Root; d__Bacteria; g__A; s__A b' + combined = self._collapse( + [[assigned(self._otu(self.AMBIGUOUS_WINDOW, ['r1', 'r2', 'r3', 'r4', 'r5'], + repaired=range(5)), genus, [genus])], + [assigned(self._otu(self.DONOR_WINDOW, ['r6']), species, [species])]], + extra_args='--resolve-ambiguous-windows') + self.assertEqual(1, len(combined.data)) + otu = list(combined)[0] + self.assertEqual(self.DONOR_WINDOW, otu.sequence) + self.assertEqual(6, otu.count) + self.assertEqual(species, otu.taxonomy) + self.assertEqual([species, species], otu.equal_best_hit_taxonomies()) + self.assertEqual('diamond', otu.taxonomy_assignment_method()) + def test_chunks_combine_unresolved_without_the_flag(self): combined = self._collapse( [[self._otu(self.AMBIGUOUS_WINDOW, ['r1'], repaired=[0])], From bf04812962b9d039455879ceeade02596fb1126c Mon Sep 17 00:00:00 2001 From: Ben Woodcroft Date: Wed, 28 Oct 2026 12:59:37 +1000 Subject: [PATCH 3/4] Fix polars 2.0 incompatibilities polars 2.0 removed the 'columns' argument of DataFrame.pivot (now 'on') and GroupBy.count (now GroupBy.len), breaking summarise's species-by-site and taxonomic-level-coverage outputs, and supplement. Both replacements are available in polars 1.x too. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019URUjTJM4Se2CmgGzPmPs5 --- singlem/condense.py | 2 +- singlem/summariser.py | 4 ++-- singlem/supplement.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/singlem/condense.py b/singlem/condense.py index cd5e25de..0494818c 100644 --- a/singlem/condense.py +++ b/singlem/condense.py @@ -1381,7 +1381,7 @@ def taxonomic_level_coverage_table(self, assume_8_levels=False): ((pl.col('coverage') / pl.col('coverage').sum()).alias('relative_abundance') * 100).round(2), ) - if assume_8_levels or len(result.select(pl.col('level')).group_by('level').count()) in [7, 8]: + if assume_8_levels or len(result.select(pl.col('level')).group_by('level').len()) in [7, 8]: # If there's 7 or 8 (including 0) levels, then assume that this is a regular taxonomy going on. levels = ['root', 'domain', 'phylum', 'class', 'order', 'family', 'genus', 'species'] level_id_to_level_name = {i: levels[i] for i in range(len(levels))} diff --git a/singlem/summariser.py b/singlem/summariser.py index 28167370..2b401b91 100644 --- a/singlem/summariser.py +++ b/singlem/summariser.py @@ -717,7 +717,7 @@ def write_species_by_site_table(**kwargs): if target_level is not None: profiles = pl.concat(all_profiles[level_index]) num_samples = len(profiles) - profiles.pivot(index=['taxonomy'], columns=['sample'], values=['relative_abundance']).fill_null(0).write_csv(output_species_by_site_relative_abundance_table, separator='\t') + profiles.pivot(index=['taxonomy'], on=['sample'], values=['relative_abundance']).fill_null(0).write_csv(output_species_by_site_relative_abundance_table, separator='\t') logging.info("Wrote site by species table for {} samples".format(num_samples)) else: # Write out a CSV for each level, with the name being the prefix, and the level being the suffix @@ -725,7 +725,7 @@ def write_species_by_site_table(**kwargs): num_samples = len(profiles) profiles = pl.concat(profiles) output_filename = output_species_by_site_relative_abundance_table+'-'+levels[level]+'.tsv' - profiles.pivot(index=['taxonomy'], columns=['sample'], values=['relative_abundance']).fill_null(0).write_csv(output_filename, separator='\t') + profiles.pivot(index=['taxonomy'], on=['sample'], values=['relative_abundance']).fill_null(0).write_csv(output_filename, separator='\t') logging.info("Wrote site by species table {} for {} sample(s) at level {}".format( output_filename, num_samples, levels[level])) diff --git a/singlem/supplement.py b/singlem/supplement.py index a63a23e7..20099690 100755 --- a/singlem/supplement.py +++ b/singlem/supplement.py @@ -392,7 +392,7 @@ def run_hmmsearch_on_one_genome(lock, data, matched_transcripts_fna, output_matc df2 = df.sort('bitscore', descending=True).group_by('transcript').first() # Remove genes where there is >1 hit # If any HMM is associated with multiple genes, remove it as being potentially dubious - ok_genes = list(df2.group_by('hmm').count().filter(pl.col('count') == 1).select('hmm'))[0] + ok_genes = list(df2.group_by('hmm').len().filter(pl.col('len') == 1).select('hmm'))[0] df3 = df2.filter(pl.col('hmm').is_in(ok_genes)) logging.debug("After filtering, {} hits remain, from {} total hits".format(len(df3), len(df2))) From e841c415ac0dfd5ba89f14ac7828fb4793d8fe93 Mon Sep 17 00:00:00 2001 From: Ben Woodcroft Date: Wed, 28 Oct 2026 15:43:28 +1000 Subject: [PATCH 4/4] Test that pipe profiles weebill species when no marker is hit pipe already continues to condense with an empty OTU table when no read hits a marker, and condense now profiles weebill-only samples, so a run whose reads hit no marker still gets its weebill species. Test this end to end with a stand-in weebill. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_019URUjTJM4Se2CmgGzPmPs5 --- test/test_pipe.py | 52 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/test/test_pipe.py b/test/test_pipe.py index 65229244..6e3dace2 100644 --- a/test/test_pipe.py +++ b/test/test_pipe.py @@ -1356,6 +1356,58 @@ def test_read_chunk_size_forward_gzip(self): list([line.split("\t") for line in expected]), extern.run(cmd).replace('chunked_reads','')) + def test_weebill_profile_without_marker_hits(self): + # The reads hit no marker, but weebill detects a species in them, which + # should still be profiled rather than pipe stopping for lack of hits. + with tempfile.TemporaryDirectory() as td: + metapackage = os.path.join(td, 'a.smpkg') + shutil.copytree(os.path.join(path_to_data, 'archaeal_small.v4.smpkg'), metapackage) + contents_path = os.path.join(metapackage, 'CONTENTS.json') + with open(contents_path) as f: + contents = json.load(f) + contents.update({ + 'singlem_metapackage_version': 7, + 'taxonomy_database_name': 'test', 'taxonomy_database_version': None, + 'diamond_prefilter_performance_parameters': None, + 'diamond_taxonomy_assignment_performance_parameters': None, + 'makeidx_sensitivity_params': None, 'avg_num_genes_per_species': None, + 'weebill_dbs': [{'db': 'fake.syl2db', 'c': 100}]}) + with open(contents_path, 'w') as f: + json.dump(contents, f) + open(os.path.join(metapackage, 'fake.syl2db'), 'w').close() + + # Stands in for weebill, reporting one genome the metapackage knows. + bin_dir = os.path.join(td, 'bin') + os.mkdir(bin_dir) + fake_weebill = os.path.join(bin_dir, 'weebill') + with open(fake_weebill, 'w') as f: + f.write('#!/bin/bash\n' + 'while [ $# -gt 0 ]; do\n' + ' if [ "$1" = -o ]; then out=$2; fi\n' + ' if [ "$1" = -r ]; then reads=$2; fi\n' + ' shift\n' + 'done\n' + 'printf "Sample_file\\tGenome_file\\tTrue_cov\\n%s\\tGCA_011042755.1_genomic.fna.gz\\t3.5\\n" ' + '"$reads" > "$out"\n') + os.chmod(fake_weebill, 0o755) + + reads = os.path.join(td, 'reads.fa') + with open(reads, 'w') as f: + f.write('>r1\n{}\n>r2\n{}\n'.format('ACGT' * 50, 'TTGCA' * 40)) + otu_table = os.path.join(td, 'otu_table.tsv') + profile = os.path.join(td, 'profile.tsv') + extern.run('PATH={}:$PATH {} pipe --sequences {} --metapackage {} --otu-table {} -p {}'.format( + bin_dir, path_to_script, reads, metapackage, otu_table, profile)) + with open(otu_table) as f: + self.assertEqual(['gene\tsample\tsequence\tnum_hits\tcoverage\ttaxonomy\n'], f.readlines()) + with open(profile) as f: + rows = [line.rstrip('\n').split('\t') for line in f] + self.assertEqual(['sample', 'coverage', 'taxonomy'], rows[0]) + self.assertIn( + ['reads', '3.5', 'Root; d__Archaea; p__Thermoproteota; c__Korarchaeia; o__Korarchaeales; ' + 'f__Korarchaeaceae; g__DRBY01; s__DRBY01 sp011042755'], + rows[1:]) + def test_read_chunk_number_incompatible_with_weebill(self): # weebill would see the whole of the reads, not just the chunk. with tempfile.TemporaryDirectory() as td: