Spaces:
Running
Running
| """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 | |
| # -------------------------------------------------------------------------- | |
| 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") | |