"""Reading what the user submits: sequences and an optional count table. Everything here is pure Python and independent of the model, so it is tested directly (test_inputs.py). """ from __future__ import annotations import csv import io import math import re import zlib from dataclasses import dataclass, field from typing import Dict, List, Optional, Sequence, Tuple from taxonova.data.fasta import parse_sequence_id, parse_taxonomy LEVEL_ORDER = ["domain", "phylum", "class", "order", "family", "genus", "species"] #: Most sequences one submission may hold. Measured on 2026-10-30 with two #: threads, 500 V4 sequences took 5.5 s end to end (screen plus model), so #: 20,000 take roughly 4 minutes; long reads take longer. Larger sets belong #: on the command line, which has no limit. MAX_SEQUENCES = 20_000 #: Largest upload accepted, per file, as sent (compressed or not). MAX_UPLOAD_BYTES = 50 * 1024 * 1024 #: Largest text accepted after gzip decompression. MAX_TEXT_BYTES = 120 * 1024 * 1024 #: Most samples in a count table. MAX_SAMPLES = 5_000 _SIZE_RE = re.compile(r"(?:^|;)size=(\d+(?:\.\d+)?)(?:;|$)", re.IGNORECASE) #: Characters a sequence may hold: IUPAC nucleotide codes in either case, #: plus what taxonova's FASTA reader removes (alignment gaps '-' and '.', #: whitespace). The package reads any other character as N; the web app #: refuses it instead, because such input (protein, a pasted table, a #: binary file) is not nucleotide sequence and would otherwise come back #: as a page of Unclassified rows with no reason given. _NOT_NUCLEOTIDE = re.compile(r"[^ACGTURYSWKMBDHVNacgturyswkmbdhvn.\-\s]") #: Records named in the alphabet error. MAX_BAD_RECORDS_SHOWN = 5 class InputError(ValueError): """A problem with the user's input, worded for the user.""" Triple = Tuple[str, str, dict] # -------------------------------------------------------------------------- # Bytes -> text # -------------------------------------------------------------------------- def decode_upload(data: bytes, what: str = "file") -> str: """Decode an upload to text, inflating gzip if needed.""" if len(data) > MAX_UPLOAD_BYTES: raise InputError(f"The {what} is larger than {MAX_UPLOAD_BYTES // (1024 * 1024)} MB.") if data[:2] == b"\x1f\x8b": data = _gunzip(data, what) elif data[:3] == b"BZh" or data[:4] == b"PK\x03\x04": raise InputError(f"The {what} is compressed in a format other than gzip. Use .gz or plain text.") try: return data.decode("utf-8-sig") except UnicodeDecodeError: try: return data.decode("latin-1") except UnicodeDecodeError: # pragma: no cover - latin-1 decodes anything raise InputError(f"The {what} is not a text file.") def _gunzip(data: bytes, what: str) -> bytes: """Inflate every gzip member (bgzip and concatenated .gz files have several).""" limit = MAX_TEXT_BYTES + 1 parts: List[bytes] = [] total = 0 rest = data while rest: if rest[:2] != b"\x1f\x8b": if rest.strip(b"\x00"): raise InputError(f"The {what} has data after the end of the gzip stream.") break d = zlib.decompressobj(16 + zlib.MAX_WBITS) try: out = d.decompress(rest, limit - total) except zlib.error: raise InputError(f"The {what} looks like gzip but could not be decompressed.") total += len(out) parts.append(out) if total > MAX_TEXT_BYTES or d.unconsumed_tail: raise InputError( f"The {what} is larger than {MAX_TEXT_BYTES // (1024 * 1024)} MB once decompressed.") if not d.eof: raise InputError(f"The {what} is a truncated gzip file.") rest = d.unused_data return b"".join(parts) # -------------------------------------------------------------------------- # Sequences # -------------------------------------------------------------------------- def read_fasta_string(text: str) -> List[Triple]: """Parse FASTA text into (seq_id, seq, taxonomy) triples. Headers are parsed with taxonova's own functions, so ids are the first word of the header, as on the command line. """ sequences = [] seq_id = header = None seq_lines: List[str] = [] for line in text.splitlines(): line = line.strip() if not line: continue if line.startswith(">"): if seq_id is not None: sequences.append((seq_id, "".join(seq_lines), parse_taxonomy(header, expect_taxonomy=False))) header = line seq_id = parse_sequence_id(line) seq_lines = [] elif line.startswith(";"): continue # old-style FASTA comment line, skipped as taxonova's reader does elif seq_id is not None: seq_lines.append(line) else: raise InputError("The FASTA text has sequence lines before the first '>' header.") if seq_id is not None: sequences.append((seq_id, "".join(seq_lines), parse_taxonomy(header, expect_taxonomy=False))) return sequences def prepare_plain_sequences(text: str) -> List[Triple]: """One sequence per line, no headers: name them seq1, seq2, ...""" unknown = {lvl: "unknown" for lvl in LEVEL_ORDER} return [(f"seq{i}", line.strip(), dict(unknown)) for i, line in enumerate(text.splitlines(), start=1) if line.strip()] def parse_sequences(text: str) -> List[Triple]: """FASTA or one-sequence-per-line text -> triples, with clear errors.""" stripped = text.lstrip() if not stripped: return [] if stripped.startswith("@"): raise InputError( "This looks like FASTQ. Taxonova classifies sequences in FASTA format, " "for example the ASV or OTU sequences from DADA2, QIIME 2 or VSEARCH.") seqs = read_fasta_string(text) if ">" in stripped[:1] or "\n>" in text else prepare_plain_sequences(text) empty = [s for s, q, _ in seqs if not q] if empty: raise InputError(f"{len(empty)} sequence(s) have a header but no sequence, " f"for example '{empty[0][:60]}'.") check_alphabet(seqs) if len(seqs) > MAX_SEQUENCES: raise InputError( f"The input holds {len(seqs):,} sequences; the web app accepts up to " f"{MAX_SEQUENCES:,} per run. Split the file, or use the command-line tool " f"(pip install taxonova), which has no limit.") return seqs def check_alphabet(sequences: Sequence[Triple]) -> None: """Refuse the whole input if any record holds a character that is not a nucleotide code. A file that is mostly valid is refused too, rather than classified with the bad records left out or read as N: a stray character usually means the wrong file, a protein FASTA or a damaged export, and a silent partial run would give abundances over a set the user did not submit. The error names up to ``MAX_BAD_RECORDS_SHOWN`` records and the first offending character, so the file can be fixed. """ bad = [] first = None for sid, seq, _ in sequences: m = _NOT_NUCLEOTIDE.search(seq) if m: bad.append(sid) if first is None: first = (sid, m.group(0), m.start() + 1) if not bad: return shown = ", ".join(f"'{s[:40]}'" for s in bad[:MAX_BAD_RECORDS_SHOWN]) more = f" and {len(bad) - MAX_BAD_RECORDS_SHOWN:,} more" if len(bad) > MAX_BAD_RECORDS_SHOWN else "" sid, ch, pos = first shown_ch = ch if ch.isprintable() else ch.encode("unicode_escape").decode("ascii") raise InputError( f"{len(bad):,} sequence(s) contain characters that are not nucleotide codes " f"(IUPAC A C G T U R Y S W K M B D H V N; gaps '-' and '.' are allowed): {shown}{more}. " f"The first is '{sid[:40]}': '{shown_ch}' at position {pos:,}. Nothing was classified. " "Taxonova classifies nucleotide sequences (DNA or RNA) only; correct or remove these " "records and submit again.") def rna_to_dna(sequences: List[Triple]) -> Tuple[List[Triple], int]: """Replace U with T. Taxonova's models were trained on DNA, so RNA input otherwise gets different, weaker calls. Returns the triples and how many sequences were changed.""" out, changed = [], 0 for sid, seq, tax in sequences: if "U" in seq or "u" in seq: seq = seq.replace("U", "T").replace("u", "t") changed += 1 out.append((sid, seq, tax)) return out, changed def size_annotation(seq_id: str) -> Optional[float]: """The N of a ';size=N' header annotation, if present.""" m = _SIZE_RE.search(seq_id) return float(m.group(1)) if m else None def strip_size(seq_id: str) -> str: """'ASV1;size=120;' -> 'ASV1'.""" return _SIZE_RE.sub(";", seq_id).strip(";") def normalise_sequence(seq: str) -> str: return re.sub(r"[\s.\-]", "", seq).upper().replace("U", "T") # -------------------------------------------------------------------------- # Count tables # -------------------------------------------------------------------------- @dataclass class CountTable: """Counts per feature and sample, after matching to the submitted sequences.""" samples: List[str] #: Input position -> counts per sample (only matched features). counts: Dict[int, List[float]] orientation: str # "features-as-rows" | "features-as-columns" | "size-annotation" keyed_by: str # "id" | "sequence" warnings: List[str] = field(default_factory=list) def _sniff_delimiter(lines: Sequence[str], filename: str) -> str: if re.sub(r"\.gz$", "", filename.lower()).endswith(".csv"): return "," sample = "\n".join(lines[:5]) return "\t" if sample.count("\t") >= sample.count(",") else "," def _as_number(cell: str) -> Optional[float]: """A count, 0.0 for empty or NA, None if not a finite number.""" cell = cell.strip() if cell == "" or cell.upper() in ("NA", "NAN"): return 0.0 try: v = float(cell) except ValueError: return None return v if math.isfinite(v) else None def _is_non_finite(cell: str) -> bool: try: return not math.isfinite(float(cell.strip())) except ValueError: return False def read_table(text: str, filename: str = "") -> Tuple[str, List[str], List[str], List[List[str]]]: """Text -> (corner, header cells, row keys, row cells). Comments dropped. Accepts BIOM-converted TSV ('# Constructed from biom file' then '#OTU ID\\t...'), QIIME 2 exports, DADA2 seqtab written by R (empty corner cell) and plain TSV/CSV. """ lines = [ln for ln in text.splitlines() if ln.strip()] lines = [ln for ln in lines if not (ln.startswith("#") and not re.match(r"#\s*(OTU|ASV|Feature|Sample)", ln, re.I))] if len(lines) < 2: raise InputError("The count table needs a header row and at least one data row.") delim = _sniff_delimiter(lines, filename) rows = list(csv.reader(io.StringIO("\n".join(lines)), delimiter=delim)) header = [c.strip() for c in rows[0]] corner = header[0].lstrip("#").strip() names = header[1:] keys, cells = [], [] for r in rows[1:]: if not r or not any(c.strip() for c in r): continue keys.append(r[0].strip()) cells.append(r[1:]) # R's write.table without col.names=NA leaves the header one short. if cells and len(names) == len(cells[0]) - 1: names = [header[0]] + names corner = "" return corner, names, keys, cells def _numeric_matrix(names: List[str], cells: List[List[str]], warnings: List[str]): """Drop wholly non-numeric columns (e.g. a taxonomy column); error on mixed ones.""" ncol = len(names) keep = [] for j in range(ncol): vals = [_as_number(r[j]) if j < len(r) else 0.0 for r in cells] inf = [i for i, r in enumerate(cells) if j < len(r) and vals[i] is None and _is_non_finite(r[j])] if inf: raise InputError( f"Column '{names[j]}' has a value that is not a count on data row {inf[0] + 1}: " f"'{cells[inf[0]][j][:40]}'.") if all(v is None for v in vals): warnings.append(f"Column '{names[j]}' is not numeric and was ignored.") continue bad = [i for i, v in enumerate(vals) if v is None] if bad: raise InputError( f"Column '{names[j]}' has a value that is not a count on data row {bad[0] + 1}: " f"'{cells[bad[0]][j][:40]}'.") if any(v < 0 for v in vals): raise InputError(f"Column '{names[j]}' has negative counts.") keep.append((names[j], vals)) if not keep: raise InputError("The count table has no numeric columns.") return [k for k, _ in keep], [list(col) for _, col in keep] def match_count_table(text: str, sequences: Sequence[Triple], filename: str = "") -> CountTable: """Parse a count table and match it to the submitted sequences. Features may be rows or columns, and may be keyed by sequence id (with or without ';size=') or by the sequence itself (DADA2 seqtab). The orientation and key type with the most matches win. """ warnings: List[str] = [] corner, names, keys, cells = read_table(text, filename) by_id: Dict[str, int] = {} by_seq: Dict[str, int] = {} dup_ids = 0 for i, (sid, seq, _) in enumerate(sequences): for k in {sid, strip_size(sid)}: if k in by_id and by_id[k] != i: dup_ids += 1 by_id.setdefault(k, i) by_seq.setdefault(normalise_sequence(seq), i) def score(labels): ids = sum(1 for x in labels if x in by_id) seqs = sum(1 for x in labels if normalise_sequence(x) in by_seq) if labels else 0 return max(ids, seqs), ("id" if ids >= seqs else "sequence") row_score, row_key = score(keys) col_score, col_key = score(names) if row_score == 0 and col_score == 0: ex_tab = ", ".join(f"'{k[:30]}'" for k in keys[:3]) ex_fa = ", ".join(f"'{s[0][:30]}'" for s in sequences[:3]) raise InputError( "None of the feature IDs in the count table match the submitted sequences. " f"The table starts with {ex_tab}; the FASTA starts with {ex_fa}. Use the same " "IDs (or the sequences themselves) in both files.") if row_score >= col_score: orientation, keyed_by = "features-as-rows", row_key samples, columns = _numeric_matrix(names, cells, warnings) feature_labels = keys per_feature = [[col[i] for col in columns] for i in range(len(keys))] else: orientation, keyed_by = "features-as-columns", col_key feature_labels = names sample_names = keys # Transpose: rows are samples here. mat = [] for r_i, r in enumerate(cells): vals = [] for j in range(len(names)): v = _as_number(r[j]) if j < len(r) else 0.0 if v is None: raise InputError( f"Sample '{sample_names[r_i]}' has a value that is not a count for " f"'{names[j][:30]}': '{r[j][:40]}'.") if v < 0: raise InputError(f"Sample '{sample_names[r_i]}' has negative counts.") vals.append(v) mat.append(vals) samples = sample_names per_feature = [[mat[s][f] for s in range(len(samples))] for f in range(len(names))] if len(samples) > MAX_SAMPLES: raise InputError(f"The count table has {len(samples):,} samples; the limit is {MAX_SAMPLES:,}.") if len(set(samples)) != len(samples): raise InputError("Sample names in the count table are not unique.") lookup = by_id if keyed_by == "id" else by_seq counts: Dict[int, List[float]] = {} unmatched_reads = 0.0 unmatched = 0 seen_keys: set = set() duplicates: List[str] = [] for label, vals in zip(feature_labels, per_feature): key = label if keyed_by == "id" else normalise_sequence(label) if key in seen_keys: duplicates.append(label) seen_keys.add(key) pos = lookup.get(key) if pos is None: unmatched += 1 unmatched_reads += sum(vals) continue if pos in counts: counts[pos] = [a + b for a, b in zip(counts[pos], vals)] else: counts[pos] = list(vals) total = sum(sum(v) for v in counts.values()) + unmatched_reads if unmatched: pct = 100 * unmatched_reads / total if total else 0 warnings.append( f"{unmatched:,} feature(s) in the count table are not among the submitted " f"sequences and were left out ({pct:.1f}% of counts).") missing = len(sequences) - len(counts) if missing: warnings.append(f"{missing:,} submitted sequence(s) have no counts in the table " f"and are shown with zero abundance.") if duplicates: warnings.append( f"{len(duplicates):,} feature(s) occur more than once in the count table, for " f"example '{duplicates[0][:40]}'; their counts were added together.") if dup_ids and keyed_by == "id": warnings.append("Some sequence IDs occur more than once in the FASTA; their counts " "were assigned to the first occurrence.") return CountTable(samples=list(samples), counts=counts, orientation=orientation, keyed_by=keyed_by, warnings=warnings) def counts_from_size_annotations(sequences: Sequence[Triple]) -> Optional[CountTable]: """A one-sample table from ';size=N' headers, when every header has one.""" sizes = [size_annotation(s[0]) for s in sequences] if not sequences or any(v is None for v in sizes): return None return CountTable(samples=["size"], counts={i: [v] for i, v in enumerate(sizes)}, orientation="size-annotation", keyed_by="id")