diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 4de7c09c..e4cfa5c0 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -84,4 +84,47 @@ jobs: files: ./coverage.xml fail_ci_if_error: false + min-deps: + # Installs every runtime dependency pinned to its declared lower bound + # (e.g. typer==0.12.4) and runs the CLI smoke tests. Guards against the + # floor being too low. + runs-on: ubuntu-latest + + steps: + - name: Checkout repository + uses: actions/checkout@v7 + + - name: Set up Python + uses: actions/setup-python@v7 + with: + python-version: '3.11' + cache: pip + + - name: Build minimum-dependency constraints file + shell: python + run: | + import re, tomllib + from pathlib import Path + + pyproject = tomllib.loads(Path('pyproject.toml').read_text()) + deps = pyproject['project']['dependencies'] + # Each entry looks like "typer>=0.12.4" (some carry a trailing comment). + # Pin to the declared lower bound: "typer>=0.12.4" -> "typer==0.12.4". + pins = [] + for dep in deps: + m = re.match(r"^\s*([A-Za-z0-9_.-]+)\s*>=\s*([0-9][0-9A-Za-z.\-]*)", dep) + if not m: + raise SystemExit(f"could not parse lower bound: {dep!r}") + pins.append(f"{m.group(1)}=={m.group(2)}") + Path('min-constraints.txt').write_text("\n".join(pins) + "\n") + print("\n".join(pins)) + + - name: Install minimum runtime dependencies + run: | + python -m pip install --upgrade pip + python -m pip install -c min-constraints.txt -e ".[dev]" + + - name: Run CLI smoke tests against minimum dependencies + run: pytest tests/test_cli_smoke.py -v + diff --git a/.gitignore b/.gitignore index 3dbd2725..7a27a4dd 100644 --- a/.gitignore +++ b/.gitignore @@ -46,4 +46,5 @@ graphify-out coverage.xml TODO.md production/ -monitoring/ \ No newline at end of file +monitoring/ +docs/RELEASE_SUMMARY_RESPRO.md \ No newline at end of file diff --git a/README.md b/README.md index 3fc083ae..b97e9693 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,15 @@ conda create -n respro conda activate respro conda install bioconda::respro ``` + +Install via Docker (BioContainers): + +```bash +docker pull quay.io/biocontainers/respro +``` + Install via pip: + ```bash git clone https://github.com/the-foxlab/ResistanceProfiler pip install -e ".[dev]" diff --git a/docs/docs/cli-reference.md b/docs/docs/cli-reference.md index c91fff15..9c0ee3e3 100644 --- a/docs/docs/cli-reference.md +++ b/docs/docs/cli-reference.md @@ -56,6 +56,21 @@ respro init \ If your dataset only contains atomic mutation rules, omit `--formula-rules`. +Optionally ship a per-database example consensus FASTA that users can profile +with a single command (see [Profile FASTA input](#profile-fasta-input)) and that +the webapp exposes as an "Example" button. The example must be a single-record +FASTA: + +```bash +respro init \ + --name "Docs Demo" \ + --genbank some_reference.gb \ + --rules rules.tsv \ + --formula-rules combinatorial_rules.tsv \ + --example example_consensus.fasta \ + --output myrespro.db +``` + For metadata and interpretation algorithm options, see [Database Preparation](database-preparation.md). ## Extend or validate rules in an existing project @@ -79,6 +94,16 @@ respro add \ --formula-rules combinatorial_rules.tsv ``` +Use `--example example_consensus.fasta` to store or overwrite a per-database +example FASTA, or `--no-example` to clear a previously stored example: + +```bash +respro add \ + --project myrespro.db \ + --rules rules.tsv \ + --example example_consensus.fasta +``` + ## Profile FASTA input ```bash @@ -91,6 +116,19 @@ respro fasta \ --export pdf ``` +To profile the example consensus FASTA stored in the project database (set via +`respro init --example` or `respro add --example`), use `--example` instead of +`--fasta`. The two options are mutually exclusive: + +```bash +respro fasta \ + --project myrespro.db \ + --example \ + --output my_output +``` + +If no example is stored, the command fails with a message. + ## Profile VCF input ```bash diff --git a/docs/docs/database-preparation.md b/docs/docs/database-preparation.md index da55e7a0..85732b77 100644 --- a/docs/docs/database-preparation.md +++ b/docs/docs/database-preparation.md @@ -33,6 +33,7 @@ Notes: - `--genbank` can be repeated for multiple files. - `--no-additional-info` skips network lookups for extra metadata. - `--metadata` accepts a JSON file with curated project metadata. See the section below for the supported keys and value rules. +- `--example` optionally stores a single-record consensus FASTA shipped with the database. Users can then profile it via `respro fasta --example` or via the webapp "Example" button. Use `respro add --example` to overwrite and `respro add --no-example` to clear it. - After initialization, later profiling runs use this database as the internal coordinate and rule source. After this command succeeds, the file `myrespro.db` should exist. diff --git a/docs/docs/how-it-works.md b/docs/docs/how-it-works.md index bc658efe..7d5a8395 100644 --- a/docs/docs/how-it-works.md +++ b/docs/docs/how-it-works.md @@ -54,6 +54,19 @@ For VCF input, variants are defined in user-provided reference coordinates. ResP - For reverse-strand features, alleles are reverse-complemented to the internal forward strand. - Anchor-changed indels (where the VCF anchor base differs from the internal reference) are automatically split into a substitution plus a canonical indel before annotation. +#### Spliced genes and unspliced queries + +When a project feature is a spliced CDS (more than one segment in `feature_segment`) and the user supplies an unspliced whole-genome query, minimap2 reports the inter-exon intron as a single large `I` (FASTA) or `D` (raw mappy) operation in the CIGAR. Without special handling this intron would be misinterpreted as a giant coding insertion (e.g. a multi-kilobase frameshift) and would crash per-exon identity. + +ResPro classifies such intron operations: + +- The exon-junction CDS offsets are derived from `feature.segments` in genomic 5'→3' order (matching the normalized CIGAR's walking order for both strands). +- A CIGAR `I` op is classified as an intron when **both** hold: (1) its CDS position coincides with a known junction offset within `alignment.intron_junction_tolerance` (default 5 nt; configurable in `defaults.toml`), and (2) its length is strictly greater than that same tolerance. +- Classified intron `I` ops are removed from the CIGAR stored on the `FeatureMatch` (producing an exon-only CIGAR) and recorded as `IntronInterval`s (CDS junction position, query span, length). +- Identity and CDS coverage are recomputed over exons only, so a perfect exons match reports ~99% identity rather than the genomic span including the intron. +- In FASTA mode the intron query span is stripped from the region before codon walking, so no intron insertion variant is emitted. Real coding insertions (any `I` op not at a junction within tolerance, or of length ≤ tolerance) are still emitted. +- In VCF mode the query-to-CDS coordinate map skips intron query positions and offsets exon-2 (and later) CDS positions past the intron span, so variants inside the intron are excluded and exon-2 variants remap to the correct CDS offset. + ### Amino-acid consequence interpretation Per-feature nucleotide changes are translated into amino-acid consequences. Supported classes include: diff --git a/docs/docs/install.md b/docs/docs/install.md index 6382892f..ac5cd506 100644 --- a/docs/docs/install.md +++ b/docs/docs/install.md @@ -31,6 +31,14 @@ respro --version Bioconda is the recommended install path if you already use conda/mamba — it handles the `mappy` native dependency automatically without a C compiler. +### via BioContainers (Docker) + +```bash +docker pull quay.io/biocontainers/respro +``` + +The BioContainers image provides a containerized CLI environment. Mount your data directories and run `respro` commands inside the container. + ## Web app ### Docker (recommended) diff --git a/docs/docs/webapp.md b/docs/docs/webapp.md index 45a001cd..1b8f9cfd 100644 --- a/docs/docs/webapp.md +++ b/docs/docs/webapp.md @@ -73,7 +73,7 @@ ships pre-built `.db` files does **not** need access to any of them. | Phase | Host | Purpose | |---|---|---| -| Maintained-DB bootstrap & updates | `raw.githubusercontent.com` | Fetch `manifest.json`, `rules.tsv`, `metadata.json`, `formula-rules.tsv` from `the-foxlab/respro-databases` | +| Maintained-DB bootstrap & updates | `raw.githubusercontent.com` | Fetch `manifest.json`, `rules.tsv`, `metadata.json`, `formula-rules.tsv`, and (when declared) `example.fasta` from `the-foxlab/respro-databases` | | Maintained-DB bootstrap | `eutils.ncbi.nlm.nih.gov` | Fetch GenBank reference records (`efetch.fcgi?db=nuccore`) referenced by a database's rules | | `respro init` enrichment | `pubchem.ncbi.nlm.nih.gov` | Resolve drug names to CIDs, descriptions, titles, structure images (PUG REST) | | `respro init` enrichment | `eutils.ncbi.nlm.nih.gov` | PubMed article summaries (`esummary.fcgi`), PMC ID conversion | @@ -216,6 +216,15 @@ In batch VCF mode, an optional BAM can be attached to each sample for per-sample | `RESPRO_WEB_MAINTAINED_BOOTSTRAP` | `false` | When `true` (or `1`/`yes`/`on`), missing maintained databases are downloaded into `data/project_databases/` at startup, and a weekly background thread checks for updates. | | `RESPRO_WEB_MAINTAINED_DB_UPDATE_INTERVAL_SECONDS` | `604800` (7 days) | Interval between maintained-database update checks. Set to `0` to disable the weekly thread (a one-time check still runs at startup when bootstrap is enabled). Invalid values fall back to the default. | +A maintained database may declare an optional `example_fasta_path` field in its +`manifest.json` entry. When present, the referenced single-record consensus +FASTA is downloaded as `example.fasta` alongside the rules and stored in the +project database. The webapp then shows an **Example** button next to the +"One sample" / "Multiple samples" choice in the Analyze tab for that database; +clicking it switches to FASTA mode and profiles the stored example in one step, +without requiring any user-supplied input file. The button is only shown for +databases that actually ship an example. + ### Minimal `.env` example For local development you can usually leave everything at defaults. For a public deployment, start from: diff --git a/pyproject.toml b/pyproject.toml index 7cac5dfb..b51570eb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,7 @@ classifiers = [ ] dependencies = [ 'click>=8.1', - 'typer>=0.12', + 'typer>=0.12.4', 'rich>=13.0', 'numpy>=1.24', 'biopython>=1.81', diff --git a/respro/cli/fasta.py b/respro/cli/fasta.py index 5c0c7526..97516e4b 100644 --- a/respro/cli/fasta.py +++ b/respro/cli/fasta.py @@ -5,6 +5,8 @@ from __future__ import annotations import logging +import os +import tempfile from pathlib import Path from typing import Annotated @@ -25,6 +27,7 @@ from respro.core.annotation import annotate_variants from respro.core.fasta_to_vcf import fasta_to_vcf from respro.core.query import resolve_fasta_query +from respro.db.rules_queries import get_project_example_fasta from respro.db.schema import open_project_db from respro.utils.cli_errors import cli_error, render_click_exception from respro.utils.logging import err_console @@ -36,9 +39,16 @@ def _profile_fasta_command( typer.Option('--project', '-p', exists=True, help='Project database.') ], consensus_fasta: Annotated[ - Path, + Path | None, typer.Option('--fasta', '-f', exists=True, help='Input consensus FASTA sequence.') - ], + ] = None, + use_example: Annotated[ + bool, + typer.Option( + '--example', + help='Profile the example consensus FASTA stored in the project database.', + ) + ] = False, sample: Annotated[ str, typer.Option('--sample', '-s', help='Sample name for the report.') @@ -87,6 +97,12 @@ def _profile_fasta_command( project_conn = None results_conn = None + if use_example and consensus_fasta is not None: + cli_error('Provide either --fasta or --example, not both.') + if not use_example and consensus_fasta is None: + cli_error('Provide --fasta or --example to specify the input consensus sequence.') + + example_temp_path: Path | None = None try: export_formats = _parse_export_formats(export) @@ -95,6 +111,25 @@ def _profile_fasta_command( if project_row is None: cli_error('No project found in the database') + # When --example is used, materialise the stored example FASTA text into a temp file so the + # existing read/align pipeline is reused unchanged. + if use_example: + example_text = get_project_example_fasta(project_conn) + if example_text is None: + cli_error( + f'No example FASTA is stored in project database {project!s}. ' + 'Add one with `respro init --example `.' + ) + assert example_text is not None + temp_fd, temp_name = tempfile.mkstemp(prefix='respro_example_', suffix='.fasta') + example_temp_path = Path(temp_name) + os.close(temp_fd) + example_temp_path.write_text(example_text) + consensus_fasta = example_temp_path + + # cli_error raises typer.Exit, but mypy cannot infer that; narrow explicitly. + assert consensus_fasta is not None + results_conn = _init_results_db_connection( str(results_db) if results_db else None, project_conn, logger, ) @@ -163,6 +198,8 @@ def _profile_fasta_command( project_conn.close() if results_conn is not None: results_conn.close() + if example_temp_path is not None: + example_temp_path.unlink(missing_ok=True) def register(app: typer.Typer) -> None: diff --git a/respro/cli/init.py b/respro/cli/init.py index effa7a7e..52995e9b 100644 --- a/respro/cli/init.py +++ b/respro/cli/init.py @@ -25,6 +25,7 @@ from respro.db.project_metadata import load_metadata_json, store_project_metadata from respro.db.schema import PROJECT_SCHEMA_VERSION, create_schema, open_project_db from respro.io.genbank import ParsedGenBankReference, parse_genbank_sources +from respro.io.reference import read_fasta from respro.utils.cli_errors import cli_error from respro.utils.files import require_file, resolve_output_file from respro.utils.logging import err_console @@ -32,6 +33,31 @@ logger = logging.getLogger(__name__) +def _read_example_fasta_text(example_fasta: Path) -> str: + """ + Read and validate a per-database example consensus FASTA file. + + The example must contain exactly one record. Returns the original file text so the stored + blob round-trips to the same file the curator provided. + + :param example_fasta: path to the example FASTA file + :return: the raw file text + :raises ValueError: if the file is empty, unparseable, or contains more than one record + """ + require_file(example_fasta, 'Example FASTA') + text = example_fasta.read_text() + try: + records = read_fasta(example_fasta) + except Exception as exc: # Biopython raises various errors for malformed FASTA + raise ValueError(f'Example FASTA could not be parsed: {exc}') from exc + if len(records) != 1: + raise ValueError( + f'Example FASTA must contain a single record; got {len(records)}. ' + 'The example is profiled as a single consensus sequence.' + ) + return text + + # ────────────────────────────────────────────────────────────────────── # Orchestration # ────────────────────────────────────────────────────────────────────── @@ -46,6 +72,7 @@ def init_project( metadata_json: Path | None = None, overwrite: bool = False, additional_info: bool = True, + example_fasta: Path | None = None, ) -> Path: """ Create and populate a new project database. @@ -58,6 +85,7 @@ def init_project( :param overwrite: if True, delete an existing database at db_path before creating a fresh one :param additional_info: if True (default), query PubChem for drug metadata and resolve publication DOIs/titles via NCBI and CrossRef; failures are non-fatal + :param example_fasta: optional single-record consensus FASTA stored with the database as an example :return: path to the created database """ if not genbank_paths: @@ -68,6 +96,7 @@ def init_project( require_file(rules_tsv, 'Rules TSV') if formula_rules_tsv is not None: require_file(formula_rules_tsv, 'Formula rules TSV') + example_text = _read_example_fasta_text(example_fasta) if example_fasta else '' metadata_payload, algorithms = load_metadata_json(metadata_json) if metadata_json else ({}, []) genbank_records = parse_genbank_sources(genbank_paths) @@ -90,6 +119,11 @@ def init_project( additional_info=additional_info, ) store_project_metadata(conn, project_id, metadata_payload) + if example_text: + conn.execute( + 'UPDATE project SET example_fasta = ? WHERE id = ?', + (example_text, project_id), + ) if algorithms: algorithms = _sanitize_effect_as_resistant_algorithms(conn, project_id, algorithms) if algorithms: @@ -121,6 +155,8 @@ def add_to_project( genbank_paths: list[Path] | None = None, additional_info: bool = True, validate_only: bool = False, + example_fasta: Path | None = None, + clear_example: bool = False, ) -> Path: """ Add curated rules and optional GenBank annotations to an existing project. @@ -130,12 +166,17 @@ def add_to_project( :param genbank_paths: optional GenBank files with additional references/features :param additional_info: if True, query PubChem for new drugs and resolve publication metadata :param validate_only: if True, run full rules validation/import path and roll back all DB changes + :param example_fasta: optional single-record consensus FASTA replacing the stored example + :param clear_example: if True, clear any stored example FASTA :return: path to the updated database """ require_file(db_path, 'Project database') require_file(rules_tsv, 'Rules TSV') if formula_rules_tsv is not None: require_file(formula_rules_tsv, 'Formula rules TSV') + if example_fasta is not None and clear_example: + raise ValueError('Provide either --example or --no-example, not both.') + example_text = _read_example_fasta_text(example_fasta) if example_fasta else '' records: list[ParsedGenBankReference] = [] for genbank_path in genbank_paths or []: @@ -180,6 +221,13 @@ def add_to_project( apply_ic50_threshold_classification(conn, project_id, ic50_config) if additional_info: _get_drugs_from_pubchem(conn, project_id) + if example_text: + conn.execute( + 'UPDATE project SET example_fasta = ? WHERE id = ?', + (example_text, project_id), + ) + elif clear_example: + conn.execute('UPDATE project SET example_fasta = ? WHERE id = ?', ('', project_id)) conn.execute( "UPDATE project SET updated_at = datetime('now') WHERE id = ?", (project_id,), @@ -306,6 +354,12 @@ def _init_command( help='Query PubChem for drug metadata and resolve publications via NCBI/CrossRef.', ) ] = True, + example: Annotated[ + Path | None, typer.Option( + '--example', '-ex', exists=True, + help='Optional single-record consensus FASTA stored with the database as an example.', + ) + ] = None, ) -> None: """ Initialise a project database from one or more GenBank reference records and resistance rules provided in TSV. @@ -327,6 +381,7 @@ def _init_command( metadata_json=metadata, overwrite=overwrite, additional_info=additional_info, + example_fasta=example, ) except (FileExistsError, FileNotFoundError, ValueError) as exc: cli_error(str(exc)) @@ -359,6 +414,15 @@ def _init_add_command( validate: Annotated[ bool, typer.Option('--validate', '-v', help='Validate rules and exit without writing DB changes.') ] = False, + example: Annotated[ + Path | None, typer.Option( + '--example', '-ex', exists=True, + help='Replace the stored example consensus FASTA with this single-record file.', + ) + ] = None, + clear_example: Annotated[ + bool, typer.Option('--no-example', help='Clear any stored example consensus FASTA.') + ] = False, ) -> None: """ Add curated rules and optional GenBank annotations to an existing project database. @@ -373,6 +437,8 @@ def _init_add_command( formula_rules_tsv=formula_rules, additional_info=additional_info, validate_only=validate, + example_fasta=example, + clear_example=clear_example, ) except (FileExistsError, FileNotFoundError, ValueError) as exc: cli_error(str(exc)) diff --git a/respro/cli/maintained_db.py b/respro/cli/maintained_db.py index 37d698c9..69312670 100644 --- a/respro/cli/maintained_db.py +++ b/respro/cli/maintained_db.py @@ -167,6 +167,7 @@ def _download_command( rules_tsv=files['rules'], # type: ignore[arg-type] formula_rules_tsv=files['formula_rules'], # type: ignore[arg-type] metadata_json=files['metadata'], # type: ignore[arg-type] + example_fasta=files['example'], # type: ignore[arg-type] overwrite=True, additional_info=additional_info, ) diff --git a/respro/config/cli_settings.py b/respro/config/cli_settings.py index 3cb03a53..0a078ae3 100644 --- a/respro/config/cli_settings.py +++ b/respro/config/cli_settings.py @@ -86,6 +86,7 @@ class CliAlignmentConfig: gap_extension_penalty_1: int gap_open_penalty_2: int gap_extension_penalty_2: int + intron_junction_tolerance: int @dataclass(frozen=True) @@ -163,6 +164,7 @@ def _bins(section: dict) -> CliAfBinsConfig: gap_extension_penalty_1=int(alignment['gap_extension_penalty_1']), gap_open_penalty_2=int(alignment['gap_open_penalty_2']), gap_extension_penalty_2=int(alignment['gap_extension_penalty_2']), + intron_junction_tolerance=int(alignment['intron_junction_tolerance']), ), ) diff --git a/respro/config/defaults.toml b/respro/config/defaults.toml index 4e90bd37..e1163856 100644 --- a/respro/config/defaults.toml +++ b/respro/config/defaults.toml @@ -66,3 +66,10 @@ mismatch_penalty = 4 gap_extension_penalty_1 = 2 gap_open_penalty_2 = 24 gap_extension_penalty_2 = 1 +# Maximum CDS-position distance (in nt) between a CIGAR `I` op and a known +# exon-junction offset for the `I` op to be classified as an intron rather +# than a coding insertion. Only applies to features with >1 segment (spliced +# CDS). mappy may place the junction a few nt off the stored segment boundary, +# so a small tolerance accommodates soft clipping near exon edges without +# classifying real in-exon insertions as introns. +intron_junction_tolerance = 5 diff --git a/respro/core/alignment.py b/respro/core/alignment.py index 4450febc..980eeecc 100644 --- a/respro/core/alignment.py +++ b/respro/core/alignment.py @@ -16,12 +16,22 @@ from respro.config.cli_settings import CLI_CONFIG from respro.db.features import load_feature_segments_by_feature_id -from respro.db.models import FeatureMatch, FeatureRecord +from respro.db.models import FeatureMatch, FeatureRecord, FeatureSegment, IntronInterval logger = logging.getLogger(__name__) _RE_CIGAR = re.compile(r'(\d+)([MIDNSHP=X])') +__all__ = [ + 'IntronInterval', + 'classify_introns', + 'exon_junction_cds_offsets', + 'match_query_to_features', + 'load_features', + 'parse_cigar', + 'cigar_to_coordinate_map', +] + # ────────────────────────────────────────────────────────────────────── # Public API # ────────────────────────────────────────────────────────────────────── @@ -186,13 +196,25 @@ def _match_with_mappy( continue h = primary[0] - identity = h.mlen / h.blen if h.blen else 0.0 - cds_coverage = (h.q_en - h.q_st) / len(cds) - query_coverage = (h.r_en - h.r_st) / len(query_upper) - strand = '+' if h.strand == 1 else '-' cigar = _normalize_mappy_cigar(h.cigar_str, strand) + # Classify intron I-ops by exon-junction position for spliced features. + # Produces an exon-only CIGAR (intron I ops removed) plus intron + # intervals carried on the FeatureMatch. query_start in each interval is + # relative to the alignment's coding-orientation query region start + # (i.e. relative to the first base the CIGAR consumes in the query); + # downstream strand-aware consumers convert to forward-strand coords. + junctions = exon_junction_cds_offsets(feature) + intron_tolerance = cfg.intron_junction_tolerance + cigar, intron_intervals = classify_introns(cigar, junctions, intron_tolerance) + intron_lengths = [iv.length for iv in intron_intervals] + + identity, cds_coverage = _recompute_exon_metrics( + h, cds, feature, intron_lengths, cigar, + ) + query_coverage = (h.r_en - h.r_st) / len(query_upper) + matches.append(FeatureMatch( feature=feature, identity=identity, @@ -203,11 +225,53 @@ def _match_with_mappy( strand=strand, cigar=cigar, cds_start=h.q_st, + intron_intervals=tuple(intron_intervals), )) return matches +def _recompute_exon_metrics( + h: mappy.Alignment, + cds: str, + feature: FeatureRecord, + intron_lengths: list[int], + exon_only_cigar: str, +) -> tuple[float, float]: + """ + Recompute identity and CDS coverage over exons only for spliced features. + + When introns were classified, the intron ``D`` ops (deletions in the CDS + query = intron spans in the genome reference) would normally be counted in + ``h.blen`` and crash identity. However mappy reports ``h.blen`` as the + query (CDS) block length, which already excludes ``D`` ops, so the intron + does not inflate it. Identity is therefore computed as ``h.mlen`` over the + aligned exon span (sum of ``M`` and ``D`` ops in the exon-only CIGAR), + which equals the coding bases assessed. CDS coverage is the aligned exon + length over the sum of exon lengths (the true coding span), clamped to + [0.0, 1.0]. Using the exon-only CIGAR avoids mappy soft-clipping inflating + coverage above 1.0. + + :param h: mappy primary alignment hit + :param cds: feature CDS sequence (spliced, coding orientation) + :param feature: feature record (for exon lengths) + :param intron_lengths: lengths of classified intron I ops + :param exon_only_cigar: exon-only CIGAR (intron I ops removed) + :return: ``(identity, cds_coverage)`` both in [0.0, 1.0] + """ + aligned_cds = sum(int(n) for n, op in re.findall(r'(\d+)([MD])', exon_only_cigar)) + if intron_lengths: + # Exon-only identity: matches over the assessed coding span. mlen + # already excludes intron D ops; aligned_cds is the exon M+D span. + identity = h.mlen / aligned_cds if aligned_cds > 0 else 0.0 + else: + identity = h.mlen / h.blen if h.blen else 0.0 + + exon_span = sum(seg.end - seg.start for seg in feature._coding_segments) or len(cds) + cds_coverage = aligned_cds / exon_span if exon_span else 0.0 + return identity, min(cds_coverage, 1.0) + + def _swap_cigar_indels(cigar: str) -> str: """ Swap I and D operations in a CIGAR string. @@ -270,6 +334,125 @@ def parse_cigar(cigar: str) -> list[tuple[int, str]]: return [(int(m.group(1)), m.group(2)) for m in _RE_CIGAR.finditer(cigar)] +def exon_junction_cds_offsets(feature: FeatureRecord) -> list[int]: + """ + Return the 0-based CDS offsets at which each exon junction occurs. + + A junction offset is the cumulative coding length up to the end of an exon; + it is the CDS position at which the next exon begins. For a single-segment + feature the list is empty. + + Junctions are derived from :attr:`FeatureRecord.segments` in **genomic + 5'->3' order** (segment_index order), not from ``_coding_segments``. This + matches the normalized CIGAR's walking order: for '+' strand the CIGAR + walks genomic order (== coding order), and for '-' strand + :func:`_normalize_mappy_cigar` reverses the CIGAR so it also walks genomic + 5'->3' order. Using ``_coding_segments`` (which reverses for '-' strand) + would mismatch the CIGAR and misclassify introns. + + :param feature: feature record with optional ``segments`` + :return: ordered list of exon-junction CDS offsets (one per internal boundary) + """ + segments = feature.segments if feature.segments else ( + FeatureSegment(segment_index=0, start=feature.start, end=feature.end), + ) + offsets: list[int] = [] + cumulative = 0 + for segment in segments[:-1]: + cumulative += segment.end - segment.start + offsets.append(cumulative) + return offsets + + +def classify_introns( + cigar: str, + junction_offsets: list[int], + tolerance: int, +) -> tuple[str, list[IntronInterval]]: + """ + Classify CIGAR ``I`` ops as introns by exon-junction position. + + Walks the normalized CIGAR (CDS=reference, genome=query) tracking the CDS + position consumed by ``M`` and ``D`` ops. An ``I`` op (insertion in the + query/genome) is classified as an intron iff **both** hold: + + 1. Its CDS position coincides with a known exon-junction offset within + ``tolerance`` nt (the primary classifier — no standalone length heuristic). + 2. Its length is strictly greater than ``tolerance``. + + The length guard is derived from ``tolerance`` rather than carried as a + separate knob: an insertion of length ≤ tolerance can never be + misclassified as an intron, regardless of where it lands relative to the + junction, while real introns (typically hundreds/thousands of nt, always + > tolerance) always satisfy the guard. Such intron ``I`` ops are removed + from the returned exon-only CIGAR and recorded as :class:`IntronInterval` + entries; real coding insertions (any ``I`` op not near a junction, or of + length ≤ tolerance) are kept. + + Adjacent same-operation runs in the exon-only CIGAR are merged so downstream + consumers receive a clean CIGAR with no spurious length-1 adjacency. + + :param cigar: normalized CIGAR string (CDS=reference, genome=query) + :param junction_offsets: exon-junction CDS offsets (from + :func:`exon_junction_cds_offsets`); empty for single-exon features + :param tolerance: maximum CDS-position distance for an ``I`` op to be + classified as an intron, and the strict lower bound on intron length + (from ``alignment.intron_junction_tolerance``) + :return: ``(exon_only_cigar, intron_intervals)`` where intron_intervals is + ordered by CDS junction position + """ + if not junction_offsets: + return cigar, [] + + junctions = sorted(junction_offsets) + introns: list[IntronInterval] = [] + exon_ops: list[tuple[int, str]] = [] + + cds_pos = 0 + query_pos = 0 # tracked for IntronInterval.query_start (relative to cigar start) + for length, op in parse_cigar(cigar): + if op == 'I': + near_junction = any( + abs(cds_pos - j) <= tolerance for j in junctions + ) + if near_junction and length > tolerance: + introns.append(IntronInterval( + cds_junction_pos=cds_pos, + query_start=query_pos, + length=length, + )) + else: + exon_ops.append((length, op)) + query_pos += length + continue + + if op == 'D': + cds_pos += length + elif op == 'M': + cds_pos += length + query_pos += length + + exon_ops.append((length, op)) + + merged = _merge_adjacent_ops(exon_ops) + exon_only_cigar = ''.join(f'{length}{op}' for length, op in merged) + return exon_only_cigar, introns + + +def _merge_adjacent_ops(ops: list[tuple[int, str]]) -> list[tuple[int, str]]: + """Merge adjacent CIGAR operations of the same type into single runs.""" + if not ops: + return [] + merged: list[tuple[int, str]] = [ops[0]] + for length, op in ops[1:]: + last_len, last_op = merged[-1] + if op == last_op: + merged[-1] = (last_len + length, op) + else: + merged.append((length, op)) + return merged + + def cigar_to_coordinate_map(cigar: str, query_start: int) -> dict[int, int | None]: """ Convert a CIGAR string to a CDS-position → query-position map. diff --git a/respro/core/fasta_to_vcf.py b/respro/core/fasta_to_vcf.py index 36c4fc3f..2e13b3ff 100644 --- a/respro/core/fasta_to_vcf.py +++ b/respro/core/fasta_to_vcf.py @@ -23,7 +23,7 @@ from Bio.Seq import Seq from respro.core.vcf_coverage import _merge_codon_gaps -from respro.db.models import CoverageGap, FeatureMatch, FeatureRecord, VariantCall +from respro.db.models import CoverageGap, FeatureMatch, FeatureRecord, IntronInterval, VariantCall logger = logging.getLogger(__name__) @@ -72,6 +72,13 @@ def _profile_feature_to_variants( if match.strand == '-': region = str(Seq(region).reverse_complement()) + # Spliced features aligned against an unspliced query carry intron intervals + # (in coding-orientation, relative to the region start). The exon-only CIGAR + # already excludes the intron I ops, so the region must also drop the intron + # spans — otherwise the CIGAR's M ops would walk exon1 straight into the + # intron bases, misaligning exon 2 and emitting spurious variants. + region = _strip_introns_from_region(region, match.intron_intervals) + aligned_cds_len = sum(int(n) for n, op in re.findall(r'(\d+)([MD])', match.cigar)) covered_cds_start = match.cds_start covered_cds_end = match.cds_start + aligned_cds_len @@ -91,6 +98,34 @@ def _profile_feature_to_variants( ) +def _strip_introns_from_region( + region: str, + intron_intervals: tuple[IntronInterval, ...], +) -> str: + """ + Remove intron spans from the coding-orientation query region. + + ``IntronInterval.query_start`` is relative to the start of the + coding-orientation region (the first base the CIGAR consumes), and + ``length`` is the intron length in that orientation. This drops those + spans so the exon-only CIGAR's ``M`` ops walk exon 1 directly into exon 2 + without crossing the intron. + + :param region: coding-orientation query slice (already RC'd for '-' strand) + :param intron_intervals: intron intervals from the FeatureMatch + :return: region with intron spans removed + """ + if not intron_intervals: + return region + # Sort by start and drop each [start, start+length) span. Build a keep-mask + # to avoid index shifts from repeated slicing. + drop = [False] * len(region) + for iv in intron_intervals: + for i in range(iv.query_start, min(iv.query_start + iv.length, len(region))): + drop[i] = True + return ''.join(base for base, d in zip(region, drop) if not d) + + def _variants_from_alignment( aligned_ref: str, aligned_query: str, diff --git a/respro/core/vcf_remap.py b/respro/core/vcf_remap.py index 153016db..9f491238 100644 --- a/respro/core/vcf_remap.py +++ b/respro/core/vcf_remap.py @@ -10,7 +10,7 @@ from respro.core.alignment import cigar_to_coordinate_map from respro.core.query import QueryRecord, pick_best_reference_id, select_matches_for_reference -from respro.db.models import FeatureMatch, FeatureRecord, VariantCall +from respro.db.models import FeatureMatch, FeatureRecord, IntronInterval, VariantCall logger = logging.getLogger(__name__) @@ -46,7 +46,7 @@ def remap_variants( for match in matches: q2c = _build_query_to_cds_map( match.cigar, match.query_start, match.query_end, - match.strand, query_len, match.cds_start, + match.strand, query_len, match.cds_start, match.intron_intervals, ) c2q = {cds_pos: query_pos for query_pos, cds_pos in q2c.items()} match_maps.append((match, q2c, c2q)) @@ -410,6 +410,7 @@ def _build_query_to_cds_map( strand: str, query_len: int, cds_start: int = 0, + intron_intervals: tuple[IntronInterval, ...] = (), ) -> dict[int, int]: """ Invert a CIGAR-based coordinate map to query-position → CDS-position. @@ -417,21 +418,52 @@ def _build_query_to_cds_map( For '-' strand matches the CIGAR was built against the reverse-complement query, so positions are first converted back to forward-strand coordinates. - :param cigar: CIGAR string from alignment + For spliced features aligned against an unspliced query, ``cigar`` is the + exon-only CIGAR (intron ``I`` ops removed) and ``intron_intervals`` carries + the intron query spans. Because the exon-only CIGAR consumes query + positions sequentially with no intron gap, exon-2 (and later) CDS positions + would otherwise map into the intron query span. Each mapped query position + is therefore shifted past the cumulative length of all introns whose + coding-orientation start precedes it, and intron query positions themselves + are excluded from the inverted map. + + :param cigar: exon-only CIGAR string from alignment :param query_start: 0-based forward-strand start (from FeatureMatch) :param query_end: 0-based forward-strand end (from FeatureMatch) :param strand: alignment strand ('+' or '-') :param query_len: total query sequence length :param cds_start: 0-based CDS offset where this alignment starts + :param intron_intervals: intron intervals (query_start relative to the + coding-orientation region start) carried on the FeatureMatch :return: mapping {forward_query_pos: cds_pos} """ - cds_to_query: dict[int, int | None] + # Intron spans in the orientation cigar_to_coordinate_map produces. + # iv.query_start is relative to the CIGAR's coding-orientation region start. + # For '-' strand the CIGAR is in RC-orientation, so the shift must happen + # in RC space BEFORE converting to forward-strand coordinates — otherwise + # the RC-orientation intron spans would be matched against forward-strand + # query positions, silently dropping exon-1 variants and mis-remapping + # exon-2 variants. if strand == '+': + orient_start = query_start + intron_spans = [ + (orient_start + iv.query_start, orient_start + iv.query_start + iv.length) + for iv in intron_intervals + ] cds_to_query = cigar_to_coordinate_map(cigar, query_start) + if intron_spans: + cds_to_query = _shift_past_introns(cds_to_query, intron_spans) else: - # Recover RC-space start from the stored forward-strand end rc_start = query_len - query_end + orient_start = rc_start + intron_spans = [ + (orient_start + iv.query_start, orient_start + iv.query_start + iv.length) + for iv in intron_intervals + ] cds_to_query_rc = cigar_to_coordinate_map(cigar, rc_start) + if intron_spans: + cds_to_query_rc = _shift_past_introns(cds_to_query_rc, intron_spans) + # Convert RC-orientation positions to forward-strand AFTER the shift. cds_to_query = {} for cds_pos, rc_pos in cds_to_query_rc.items(): if rc_pos is not None: @@ -455,3 +487,45 @@ def _build_query_to_cds_map( return query_to_cds +def _shift_past_introns( + cds_to_query: dict[int, int | None], + intron_spans: list[tuple[int, int]], +) -> dict[int, int | None]: + """ + Shift each mapped query position past the cumulative length of preceding introns. + + The exon-only CIGAR maps CDS positions to query offsets that ignore intron + query spans. For a CDS position whose raw mapped offset falls at or after + an intron's coding-orientation start, add the cumulative intron length so + the position lands in the correct exon in the full (unspliced) query. + + Mapped positions that fall *inside* an intron span (should not happen for + true exon positions, but defensively) are set to None. + + Operates in a single orientation space (the CIGAR's coding orientation: + forward for '+' strand, RC for '-' strand). The caller converts to + forward-strand coordinates afterward for '-' strand. + + :param cds_to_query: CDS-position → raw orientation-absolute query position + :param intron_spans: intron spans in the same orientation as the raw positions + :return: CDS-position → corrected query position (same orientation as input) + """ + if not intron_spans: + return cds_to_query + spans = sorted(intron_spans) + corrected: dict[int, int | None] = {} + for cds_pos, qpos in cds_to_query.items(): + if qpos is None: + corrected[cds_pos] = None + continue + # Cumulative intron length for introns whose start is at or before qpos. + shift = sum(end - start for start, end in spans if start <= qpos) + new_pos = qpos + shift + # Exclude positions that land inside an intron span after shifting. + if any(start <= new_pos < end for start, end in spans): + corrected[cds_pos] = None + continue + corrected[cds_pos] = new_pos + return corrected + + diff --git a/respro/db/_rules_features.py b/respro/db/_rules_features.py index 43b2031e..f5a31c12 100644 --- a/respro/db/_rules_features.py +++ b/respro/db/_rules_features.py @@ -187,28 +187,33 @@ def _validate_reference_amino_acids( if 0 <= pos_0based and end_pos <= len(aa_seq): actual = aa_seq[pos_0based:end_pos].upper() if actual != ref_block: - mismatch_keys.add((feature_name, position_raw, reference_identifier, ref_aa)) - mismatch_details.append( + key = (feature_name, position_raw, reference_identifier, ref_aa) + if key not in mismatch_keys: + mismatch_details.append( + f' {reference_identifier} feature {feature_name!r} pos {pos} ({coord_base}-based): ' + f'rule says {ref_aa!r}, feature sequence has {actual!r} — rule will be skipped' + ) + mismatch_keys.add(key) + else: + out_of_range_key = (feature_name, position_raw, reference_identifier, ref_aa) + if out_of_range_key not in mismatch_keys: + out_of_range.append( f' {reference_identifier} feature {feature_name!r} pos {pos} ({coord_base}-based): ' - f'rule says {ref_aa!r}, feature sequence has {actual!r} — rule will be skipped' + f'out of range (aa_sequence length = {len(aa_seq)}) — rule will be skipped' ) - else: - out_of_range.append( - f' {reference_identifier} feature {feature_name!r} pos {pos} ({coord_base}-based): ' - f'out of range (aa_sequence length = {len(aa_seq)}) — rule will be skipped' - ) + mismatch_keys.add(out_of_range_key) if out_of_range: logger.warning( - '%d rule(s) reference positions beyond the end of the annotated protein ' - 'and will be skipped:\n%s', + '%d position(s) with rules reference coordinates beyond the end of the ' + 'annotated protein and will be skipped:\n%s', len(out_of_range), '\n'.join(out_of_range), ) if mismatch_details: logger.warning( - '%d rule(s) have reference AA mismatches and will be skipped:\n%s', + '%d position(s) with rules have reference AA mismatches and will be skipped:\n%s', len(mismatch_details), '\n'.join(mismatch_details), ) diff --git a/respro/db/cache.py b/respro/db/cache.py index 95456677..e1305931 100644 --- a/respro/db/cache.py +++ b/respro/db/cache.py @@ -5,11 +5,12 @@ from __future__ import annotations import hashlib +import json import logging import sqlite3 from respro.db.features import load_feature_segments_by_feature_id -from respro.db.models import FeatureMatch, FeatureRecord +from respro.db.models import FeatureMatch, FeatureRecord, IntronInterval logger = logging.getLogger(__name__) @@ -24,6 +25,42 @@ def sequence_checksum(sequence: str) -> str: return hashlib.sha256(sequence.upper().encode()).hexdigest() +def _serialize_intron_intervals(intervals: tuple[IntronInterval, ...]) -> str: + """Serialize intron intervals to a compact JSON array string.""" + return json.dumps( + [ + { + 'cds_junction_pos': iv.cds_junction_pos, + 'query_start': iv.query_start, + 'length': iv.length, + } + for iv in intervals + ] + ) + + +def _deserialize_intron_intervals(blob: str | None) -> tuple[IntronInterval, ...]: + """Deserialize intron intervals from a JSON array string. + + Tolerates ``None`` or malformed blobs (legacy rows) by returning an empty + tuple so that unspliced-feature mappings load without intron spans. + """ + if not blob: + return tuple() + try: + records = json.loads(blob) + except (json.JSONDecodeError, TypeError): + return tuple() + return tuple( + IntronInterval( + cds_junction_pos=int(rec['cds_junction_pos']), + query_start=int(rec['query_start']), + length=int(rec['length']), + ) + for rec in records + ) + + def load_cached_mappings( conn: sqlite3.Connection, checksum: str, @@ -45,6 +82,7 @@ def load_cached_mappings( rows = conn.execute( 'SELECT qgm.feature_id, qgm.identity, qgm.cds_coverage, qgm.query_coverage, ' 'qgm.cds_start, qgm.query_start, qgm.query_end, qgm.strand, qgm.cigar, ' + 'qgm.intron_intervals, ' 'g.reference_id, g.name, g.protein, g.start, g.end, g.strand AS feature_strand, ' 'g.codon_start, g.nt_sequence, g.aa_sequence, ' 'r.accession AS reference_accession ' @@ -84,6 +122,7 @@ def load_cached_mappings( query_end=row['query_end'], strand=row['strand'], cigar=row['cigar'], + intron_intervals=_deserialize_intron_intervals(row['intron_intervals']), )) logger.info('Loaded %d cached mapping(s) for checksum %s…', len(matches), checksum[:12]) @@ -119,8 +158,8 @@ def store_mappings( conn.execute( 'INSERT OR REPLACE INTO query_feature_mapping ' '(query_ref_id, feature_id, identity, cds_coverage, query_coverage, ' - 'cds_start, query_start, query_end, strand, cigar) ' - 'VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)', + 'cds_start, query_start, query_end, strand, cigar, intron_intervals) ' + 'VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?)', ( qref_id, match.feature.id, @@ -132,6 +171,7 @@ def store_mappings( match.query_end, match.strand, match.cigar, + _serialize_intron_intervals(match.intron_intervals), ), ) diff --git a/respro/db/models.py b/respro/db/models.py index ed6e03a4..9f88422c 100644 --- a/respro/db/models.py +++ b/respro/db/models.py @@ -147,6 +147,23 @@ class FeatureSegment: end: int +@dataclass(frozen=True) +class IntronInterval: + """One intron detected between two CDS exons in a spliced-feature alignment. + + Carried on :class:`FeatureMatch` so downstream consumers (FASTA variant + emission, VCF remap, cache, visualization) can skip the intron without + re-deriving its position from the CIGAR. + """ + + cds_junction_pos: int # 0-based CDS offset where the intron begins + query_start: int # 0-based query offset relative to the alignment's + # coding-orientation region start (RC-relative for + # '-' strand); downstream consumers add the region's + # absolute start to obtain strand-specific coordinates + length: int # intron length in query nt + + @dataclass class ResistanceRule: """A single resistance rule loaded from the database.""" @@ -370,6 +387,7 @@ class FeatureMatch: strand: str cigar: str cds_start: int = 0 # first aligned CDS position (0-based); used to reconstruct gapped strings + intron_intervals: tuple[IntronInterval, ...] = field(default_factory=tuple) @dataclass(frozen=True) diff --git a/respro/db/rules_queries.py b/respro/db/rules_queries.py index f381e292..b3c91c35 100644 --- a/respro/db/rules_queries.py +++ b/respro/db/rules_queries.py @@ -671,10 +671,25 @@ def get_project_summary_for_display(conn: sqlite3.Connection) -> dict: 'SELECT id, name, uuid, created_at, schema_version, ' 'metadata_maintainers, metadata_contact, metadata_publication_pmid, ' 'metadata_publication_doi, metadata_website, metadata_description, ' - 'metadata_maintainer_update, metadata_license, metadata_tsv_checksum ' + 'metadata_maintainer_update, metadata_license, metadata_tsv_checksum, ' + 'example_fasta ' 'FROM project ORDER BY id LIMIT 1' ).fetchone() if row is None: raise ValueError('Project DB contains no project metadata.') return dict(row) + + +def get_project_example_fasta(conn: sqlite3.Connection) -> str | None: + """ + Return the stored example consensus FASTA text, or ``None`` when none is stored. + + :param conn: open project DB connection + :return: FASTA text (single record) or ``None`` if the ``example_fasta`` column is empty + """ + row = conn.execute('SELECT example_fasta FROM project ORDER BY id LIMIT 1').fetchone() + if row is None: + raise ValueError('Project DB contains no project metadata.') + value = row['example_fasta'] + return value if value else None diff --git a/respro/db/schema.py b/respro/db/schema.py index 1edabca5..9e2211ce 100644 --- a/respro/db/schema.py +++ b/respro/db/schema.py @@ -9,7 +9,7 @@ from respro.utils.files import require_file -PROJECT_SCHEMA_VERSION = 1 +PROJECT_SCHEMA_VERSION = 2 RESULTS_SCHEMA_VERSION = 1 PROJECT_SCHEMA_SQL = """\ @@ -29,6 +29,7 @@ metadata_maintainer_update TEXT NOT NULL DEFAULT '', metadata_license TEXT NOT NULL DEFAULT '', metadata_tsv_checksum TEXT NOT NULL DEFAULT '', + example_fasta TEXT NOT NULL DEFAULT '', -- optional per-database example consensus FASTA (single record) schema_version INTEGER NOT NULL DEFAULT 1 ); @@ -193,6 +194,7 @@ query_end INTEGER NOT NULL, strand TEXT NOT NULL DEFAULT '+', cigar TEXT NOT NULL, + intron_intervals TEXT NOT NULL DEFAULT '[]', UNIQUE(query_ref_id, feature_id) ); @@ -404,6 +406,7 @@ 'metadata_maintainer_update': "TEXT NOT NULL DEFAULT ''", 'metadata_license': "TEXT NOT NULL DEFAULT ''", 'metadata_tsv_checksum': "TEXT NOT NULL DEFAULT ''", + 'example_fasta': "TEXT NOT NULL DEFAULT ''", }, 'reference': { 'accession': "TEXT DEFAULT ''", @@ -462,6 +465,7 @@ 'query_feature_mapping': { 'query_coverage': 'REAL NOT NULL DEFAULT 0', 'cds_start': 'INTEGER NOT NULL DEFAULT 0', + 'intron_intervals': "TEXT NOT NULL DEFAULT '[]'", }, } diff --git a/respro/io/maintained_db.py b/respro/io/maintained_db.py index dc5f3b44..4bec2f87 100644 --- a/respro/io/maintained_db.py +++ b/respro/io/maintained_db.py @@ -122,6 +122,14 @@ def download_database_files(db_name: str, dest_dir: Path) -> dict[str, object]: f'{db_name}/formula-rules.tsv', ) + example_path: Path | None = None + if 'example.fasta' in file_map: + example_path = _download_file( + file_map['example.fasta'], + dest_dir / 'example.fasta', + f'{db_name}/example.fasta', + ) + accessions = _parse_reference_identifiers(rules_path) genbank_paths = _fetch_genbank_records(accessions, dest_dir) @@ -129,6 +137,7 @@ def download_database_files(db_name: str, dest_dir: Path) -> dict[str, object]: 'rules': rules_path, 'metadata': metadata_path, 'formula_rules': formula_rules_path, + 'example': example_path, 'genbank': genbank_paths, } @@ -161,6 +170,15 @@ def list_output_files(db_name: str) -> list[dict]: } ) + example_fasta_path = entry.get('example_fasta_path', '') + if example_fasta_path: + output_files.append( + { + 'name': 'example.fasta', + 'download_url': _resolve_manifest_path_to_download_url(example_fasta_path), + } + ) + return output_files @@ -213,6 +231,10 @@ def _validate_manifest_entry(entry: object, idx: int) -> None: if not isinstance(formula_rules_path, str): raise RuntimeError(f'Malformed manifest: databases[{idx}].formula_rules_path must be a string') + example_fasta_path = entry.get('example_fasta_path', '') + if not isinstance(example_fasta_path, str): + raise RuntimeError(f'Malformed manifest: databases[{idx}].example_fasta_path must be a string') + metadata = entry.get('metadata') if not isinstance(metadata, dict): raise RuntimeError(f'Malformed manifest: databases[{idx}].metadata must be an object') diff --git a/respro/report/alignment_visualization.py b/respro/report/alignment_visualization.py index 57074c68..8e2276d4 100644 --- a/respro/report/alignment_visualization.py +++ b/respro/report/alignment_visualization.py @@ -7,7 +7,7 @@ from markupsafe import Markup, escape from respro.core.annotation import reverse_complement -from respro.core.fasta_to_vcf import _gapped_strings_from_cigar +from respro.core.fasta_to_vcf import _gapped_strings_from_cigar, _strip_introns_from_region from respro.db.models import AnnotatedVariant, FeatureMatch @@ -89,6 +89,13 @@ def build_feature_alignments( if match.strand == '-': region = reverse_complement(region) + # Spliced features aligned against an unspliced query carry intron + # intervals (coding-orientation, relative to the region start). The + # exon-only CIGAR excludes intron I ops, so the region must drop the + # intron spans too — otherwise the CIGAR's M ops walk exon 1 into the + # intron bases and the snippet renders a giant intron gap row. + region = _strip_introns_from_region(region, match.intron_intervals) + aligned_ref_coding, aligned_query_coding = _gapped_strings_from_cigar( feature.nt_sequence.upper(), region, match.cigar, match.cds_start, ) diff --git a/respro/report/plots.py b/respro/report/plots.py index 0d2703ce..32bc734a 100644 --- a/respro/report/plots.py +++ b/respro/report/plots.py @@ -719,12 +719,14 @@ def _draw_genome_overview( label_x, label_y, feature.display_name, - ha='center', - va='bottom', + ha='left', + va='baseline', fontsize=7, fontweight='bold', color=FEATURE_HIGHLIGHTED_EDGE, zorder=200 + y, + rotation_mode='anchor', + rotation=25, ) ax.set_title('Genome overview', fontsize=8, loc='left', fontweight='bold') diff --git a/tests/test_cli_smoke.py b/tests/test_cli_smoke.py new file mode 100644 index 00000000..7cf262cc --- /dev/null +++ b/tests/test_cli_smoke.py @@ -0,0 +1,27 @@ +""" +Smoke tests for the top-level ``respro`` CLI. +""" + +from __future__ import annotations + +from typer.testing import CliRunner + +from respro import __version__ +from respro.cli.main import app + + +class TestCliSmoke: + """Minimal end-to-end checks that the CLI group builds and runs.""" + + def test_version_flag_prints_version_and_exits_zero(self): + """``respro --version`` should print the package version and exit 0.""" + result = CliRunner().invoke(app, ['--version']) + assert result.exit_code == 0, result.output + assert __version__ in result.output + + def test_help_flag_exits_zero_and_lists_commands(self): + """``respro --help`` should exit 0 and mention the application name.""" + result = CliRunner().invoke(app, ['--help']) + assert result.exit_code == 0, result.output + # The app help string contains "ResistanceProfiler". + assert 'respro' in result.output.lower() or 'resistance' in result.output.lower() diff --git a/tests/test_config.py b/tests/test_config.py new file mode 100644 index 00000000..3e0c31e5 --- /dev/null +++ b/tests/test_config.py @@ -0,0 +1,346 @@ +""" +Tests for configuration loading and validation. + +Covers: respro/config/cli_settings.py +- _load_cli_config() +- All dataclass configuration types +- defaults.toml structure validation +""" + +from __future__ import annotations + +import tomllib +from importlib.resources import files +from urllib.parse import urlparse + +from respro.config.cli_settings import ( + CliConfig, + _load_cli_config, +) + + +class TestLoadCliConfig: + """Tests for _load_cli_config().""" + + def test_loads_without_error(self): + """Should load configuration without errors.""" + config = _load_cli_config() + assert isinstance(config, CliConfig) + + def test_returns_cli_config_instance(self): + """Should return CliConfig dataclass.""" + config = _load_cli_config() + assert hasattr(config, 'timeouts') + assert hasattr(config, 'urls') + assert hasattr(config, 'parsing') + assert hasattr(config, 'matching') + assert hasattr(config, 'similarity') + assert hasattr(config, 'af_bins') + assert hasattr(config, 'af_bins_fasta') + assert hasattr(config, 'alignment') + + def test_all_sections_populated(self): + """Should populate all configuration sections.""" + config = _load_cli_config() + assert config.timeouts is not None + assert config.urls is not None + assert config.parsing is not None + assert config.matching is not None + assert config.similarity is not None + assert config.af_bins is not None + assert config.af_bins_fasta is not None + assert config.alignment is not None + + +class TestCliTimeoutConfig: + """Tests for CliTimeoutConfig.""" + + def test_pubchem_timeout_positive(self): + """Should have positive pubchem timeout.""" + config = _load_cli_config() + assert config.timeouts.pubchem > 0 + assert isinstance(config.timeouts.pubchem, int) + + def test_pubmed_timeout_positive(self): + """Should have positive pubmed timeout.""" + config = _load_cli_config() + assert config.timeouts.pubmed > 0 + assert isinstance(config.timeouts.pubmed, int) + + def test_crossref_timeout_positive(self): + """Should have positive crossref timeout.""" + config = _load_cli_config() + assert config.timeouts.crossref > 0 + assert isinstance(config.timeouts.crossref, int) + + +class TestCliUrlConfig: + """Tests for CliUrlConfig.""" + + def test_pubchem_urls_have_placeholders(self): + """Should have {cid} placeholder in pubchem URLs.""" + config = _load_cli_config() + assert '{cid}' in config.urls.pubchem_compound_page + assert '{cid}' in config.urls.pubchem_structure_png + assert '{name}' in config.urls.pubchem_cid_lookup + + def test_ncbi_urls_have_placeholders(self): + """Should have {pmid}/{identifier} placeholders in NCBI URLs.""" + config = _load_cli_config() + assert '{pmid}' in config.urls.ncbi_pubmed_esummary + assert '{identifier}' in config.urls.ncbi_pmc_idconv + assert '{accession}' in config.urls.ncbi_nuccore_efetch + + def test_crossref_url_has_placeholder(self): + """Should have {doi} placeholder in CrossRef URL.""" + config = _load_cli_config() + assert '{doi}' in config.urls.crossref_works + + def test_protein_url_has_placeholder(self): + """Should have {protein_id} placeholder in protein URL.""" + config = _load_cli_config() + assert '{protein_id}' in config.urls.ncbi_protein_page + + def test_github_raw_url_present(self): + """Should have GitHub raw URL.""" + config = _load_cli_config() + parsed = urlparse(config.urls.github_respro_db_raw) + host = parsed.hostname + assert host is not None + assert host == 'githubusercontent.com' or host.endswith('.githubusercontent.com') + + def test_all_urls_are_strings(self): + """Should all be string type.""" + config = _load_cli_config() + for field_name in dir(config.urls): + if not field_name.startswith('_'): + value = getattr(config.urls, field_name) + assert isinstance(value, str) + + +class TestCliParsingConfig: + """Tests for CliParsingConfig.""" + + def test_doi_prefixes_is_tuple(self): + """Should have DOI prefixes as tuple.""" + config = _load_cli_config() + assert isinstance(config.parsing.doi_prefixes, tuple) + + def test_doi_prefixes_not_empty(self): + """Should have at least one DOI prefix.""" + config = _load_cli_config() + assert len(config.parsing.doi_prefixes) > 0 + + def test_doi_prefixes_are_strings(self): + """Should all be strings.""" + config = _load_cli_config() + for prefix in config.parsing.doi_prefixes: + assert isinstance(prefix, str) + + +class TestCliMatchingConfig: + """Tests for CliMatchingConfig.""" + + def test_af_threshold_is_float(self): + """Should have float AF threshold.""" + config = _load_cli_config() + assert isinstance(config.matching.combination_member_af_threshold, float) + + def test_af_threshold_in_valid_range(self): + """Should be between 0 and 1.""" + config = _load_cli_config() + assert 0 < config.matching.combination_member_af_threshold < 1 + + +class TestCliSimilarityConfig: + """Tests for CliSimilarityConfig.""" + + def test_high_threshold_positive(self): + """Should have positive high threshold.""" + config = _load_cli_config() + assert config.similarity.high > 0 + assert isinstance(config.similarity.high, int) + + def test_moderate_threshold_non_negative(self): + """Should have non-negative moderate threshold.""" + config = _load_cli_config() + assert config.similarity.moderate >= 0 + assert isinstance(config.similarity.moderate, int) + + def test_high_greater_than_moderate(self): + """High threshold should be greater than moderate.""" + config = _load_cli_config() + assert config.similarity.high > config.similarity.moderate + + +class TestCliAfBinsConfig: + """Tests for CliAfBinsConfig.""" + + def test_high_bin_valid_range(self): + """High bin should be in valid range.""" + config = _load_cli_config() + high = config.af_bins.high + assert 0 <= high[0] <= high[1] <= 1 + + def test_intermediate_bin_valid_range(self): + """Intermediate bin should be in valid range.""" + config = _load_cli_config() + intermediate = config.af_bins.intermediate + assert 0 <= intermediate[0] <= intermediate[1] <= 1 + + def test_low_bin_valid_range(self): + """Low bin should be in valid range.""" + config = _load_cli_config() + low = config.af_bins.low + assert 0 <= low[0] <= low[1] <= 1 + + def test_bins_dont_overlap(self): + """Bins should not overlap.""" + config = _load_cli_config() + assert config.af_bins.low[1] < config.af_bins.intermediate[0] + assert config.af_bins.intermediate[1] < config.af_bins.high[0] + + def test_as_dict_method(self): + """Should return dict from as_dict().""" + config = _load_cli_config() + result = config.af_bins.as_dict() + assert isinstance(result, dict) + assert 'high' in result + assert 'intermediate' in result + assert 'low' in result + + def test_fasta_bins_different_from_vcf(self): + """FASTA bins should differ from VCF bins.""" + config = _load_cli_config() + assert config.af_bins.intermediate != config.af_bins_fasta.intermediate + + +class TestCliAlignmentConfig: + """Tests for CliAlignmentConfig.""" + + def test_preset_is_string(self): + """Should have string preset.""" + config = _load_cli_config() + assert isinstance(config.alignment.preset, str) + assert len(config.alignment.preset) > 0 + + def test_k_positive(self): + """Should have positive k (k-mer size).""" + config = _load_cli_config() + assert config.alignment.k > 0 + assert isinstance(config.alignment.k, int) + + def test_w_positive(self): + """Should have positive w (minimizer window).""" + config = _load_cli_config() + assert config.alignment.w > 0 + assert isinstance(config.alignment.w, int) + + def test_best_n_positive(self): + """Should have positive best_n.""" + config = _load_cli_config() + assert config.alignment.best_n > 0 + assert isinstance(config.alignment.best_n, int) + + def test_gap_open_penalty_positive(self): + """Should have positive gap open penalty.""" + config = _load_cli_config() + assert config.alignment.gap_open_penalty > 0 + assert isinstance(config.alignment.gap_open_penalty, int) + + def test_match_score_positive(self): + """Should have positive match score.""" + config = _load_cli_config() + assert config.alignment.match_score > 0 + assert isinstance(config.alignment.match_score, int) + + def test_mismatch_penalty_positive(self): + """Should have positive mismatch penalty.""" + config = _load_cli_config() + assert config.alignment.mismatch_penalty > 0 + assert isinstance(config.alignment.mismatch_penalty, int) + + def test_gap_extension_penalty_1_positive(self): + """Should have positive gap extension penalty 1.""" + config = _load_cli_config() + assert config.alignment.gap_extension_penalty_1 > 0 + assert isinstance(config.alignment.gap_extension_penalty_1, int) + + def test_gap_open_penalty_2_positive(self): + """Should have positive gap open penalty 2.""" + config = _load_cli_config() + assert config.alignment.gap_open_penalty_2 > 0 + assert isinstance(config.alignment.gap_open_penalty_2, int) + + def test_gap_extension_penalty_2_positive(self): + """Should have positive gap extension penalty 2.""" + config = _load_cli_config() + assert config.alignment.gap_extension_penalty_2 > 0 + assert isinstance(config.alignment.gap_extension_penalty_2, int) + + def test_intron_junction_tolerance_non_negative(self): + """Should have non-negative intron junction tolerance.""" + config = _load_cli_config() + assert config.alignment.intron_junction_tolerance >= 0 + assert isinstance(config.alignment.intron_junction_tolerance, int) + + +class TestDefaultsTomlStructure: + """Tests for defaults.toml file structure.""" + + def test_file_exists_and_valid_toml(self): + """Should exist and be valid TOML.""" + defaults_path = files('respro.config').joinpath('defaults.toml') + content = defaults_path.read_text(encoding='utf-8') + # Should not raise tomllib.TOMLDecodeError + payload = tomllib.loads(content) + assert payload is not None + + def test_all_required_sections_present(self): + """Should have all required sections.""" + defaults_path = files('respro.config').joinpath('defaults.toml') + payload = tomllib.loads(defaults_path.read_text(encoding='utf-8')) + required_sections = [ + 'timeouts', + 'urls', + 'parsing', + 'matching', + 'similarity', + 'af_bins', + 'af_bins_fasta', + 'alignment', + ] + for section in required_sections: + assert section in payload, f"Missing section: {section}" + + def test_no_duplicate_keys(self): + """Should not have duplicate keys in sections.""" + defaults_path = files('respro.config').joinpath('defaults.toml') + payload = tomllib.loads(defaults_path.read_text(encoding='utf-8')) + # tomllib would fail on duplicate keys, so if we got here, it's valid + assert len(payload) > 0 + + def test_timeouts_section_structure(self): + """Should have correct timeouts structure.""" + defaults_path = files('respro.config').joinpath('defaults.toml') + payload = tomllib.loads(defaults_path.read_text(encoding='utf-8')) + timeouts = payload['timeouts'] + assert 'pubchem' in timeouts + assert 'pubmed' in timeouts + assert 'crossref' in timeouts + + def test_urls_section_structure(self): + """Should have correct urls structure.""" + defaults_path = files('respro.config').joinpath('defaults.toml') + payload = tomllib.loads(defaults_path.read_text(encoding='utf-8')) + urls = payload['urls'] + assert 'pubchem_compound_page' in urls + assert 'ncbi_pubmed_esummary' in urls + assert 'crossref_works' in urls + + def test_alignment_section_has_intron_junction_tolerance(self): + """Should have intron_junction_tolerance in alignment section.""" + defaults_path = files('respro.config').joinpath('defaults.toml') + payload = tomllib.loads(defaults_path.read_text(encoding='utf-8')) + alignment = payload['alignment'] + assert 'intron_junction_tolerance' in alignment diff --git a/tests/test_example_fasta.py b/tests/test_example_fasta.py new file mode 100644 index 00000000..f45be91e --- /dev/null +++ b/tests/test_example_fasta.py @@ -0,0 +1,275 @@ +""" +Tests for the per-database example FASTA feature (Feature: db-example-fasta). + +Covers the schema storage layer: the ``project.example_fasta`` column, its auto-migration +on open, the ``get_project_example_fasta`` accessor, and the ``get_project_summary_for_display`` +SELECT surfacing the column. Also covers ``init_project``/``add_to_project`` example wiring. +""" + +from __future__ import annotations + +import textwrap +from pathlib import Path + +import pytest +from conftest import TINY_REF_SEQ, write_genbank + +from respro.cli.init import add_to_project, init_project +from respro.db.rules_queries import get_project_example_fasta, get_project_summary_for_display +from respro.db.schema import PROJECT_SCHEMA_SQL, create_schema, open_project_db + +_EXAMPLE_FASTA = '>example_sample\nATGAAAGCTTTTGGCCCC\n' + + +class TestExampleFastaSchema: + def test_new_db_has_example_fasta_column(self, tmp_path: Path) -> None: + db_path = tmp_path / 'fresh.db' + conn = create_schema(db_path) + columns = { + row['name'] + for row in conn.execute('PRAGMA table_info(project)').fetchall() + } + conn.close() + + assert 'example_fasta' in columns + + def test_example_fasta_defaults_empty_for_new_db(self, tmp_path: Path) -> None: + db_path = tmp_path / 'fresh.db' + conn = create_schema(db_path) + conn.execute("INSERT INTO project (name, schema_version, uuid) VALUES ('P', 1, 'u')") + value = conn.execute('SELECT example_fasta FROM project').fetchone()[0] + conn.close() + + assert value == '' + + def test_open_project_db_auto_adds_example_fasta_column(self, tmp_path: Path) -> None: + # Build a complete DB, then strip the example_fasta column to simulate a legacy DB. + db_path = tmp_path / 'legacy.db' + conn = create_schema(db_path) + conn.execute("INSERT INTO project (name, schema_version, uuid) VALUES ('Legacy', 1, 'u')") + conn.execute( + 'CREATE TABLE project_legacy AS ' + 'SELECT id, name, uuid, created_at, updated_at, schema_version, ' + 'metadata_maintainers, metadata_contact, metadata_publication_pmid, ' + 'metadata_publication_doi, metadata_website, metadata_description, ' + 'metadata_maintainer_update, metadata_license, metadata_tsv_checksum ' + 'FROM project' + ) + conn.execute('DROP TABLE project') + conn.execute('ALTER TABLE project_legacy RENAME TO project') + conn.commit() + conn.close() + + migrated_conn = open_project_db(db_path) + columns = { + row['name'] + for row in migrated_conn.execute('PRAGMA table_info(project)').fetchall() + } + value = migrated_conn.execute('SELECT example_fasta FROM project').fetchone()[0] + migrated_conn.close() + + assert 'example_fasta' in columns + assert value == '' + + def test_get_project_example_fasta_returns_none_when_empty(self, tmp_path: Path) -> None: + db_path = tmp_path / 'fresh.db' + conn = create_schema(db_path) + conn.execute("INSERT INTO project (name, schema_version, uuid) VALUES ('P', 1, 'u')") + + assert get_project_example_fasta(conn) is None + conn.close() + + def test_get_project_example_fasta_returns_stored_text(self, tmp_path: Path) -> None: + db_path = tmp_path / 'fresh.db' + conn = create_schema(db_path) + conn.execute("INSERT INTO project (name, schema_version, uuid) VALUES ('P', 1, 'u')") + conn.execute('UPDATE project SET example_fasta = ? WHERE id = 1', (_EXAMPLE_FASTA,)) + + assert get_project_example_fasta(conn) == _EXAMPLE_FASTA + conn.close() + + def test_get_project_example_fasta_raises_when_no_project_row(self, tmp_path: Path) -> None: + db_path = tmp_path / 'fresh.db' + conn = create_schema(db_path) + + with pytest.raises(ValueError): + get_project_example_fasta(conn) + conn.close() + + def test_summary_for_display_includes_example_fasta(self, tmp_path: Path) -> None: + db_path = tmp_path / 'fresh.db' + conn = create_schema(db_path) + conn.execute("INSERT INTO project (name, schema_version, uuid) VALUES ('P', 1, 'u')") + conn.execute('UPDATE project SET example_fasta = ? WHERE id = 1', (_EXAMPLE_FASTA,)) + + summary = get_project_summary_for_display(conn) + conn.close() + + assert summary['example_fasta'] == _EXAMPLE_FASTA + + def test_example_fasta_column_in_schema_sql(self) -> None: + assert 'example_fasta' in PROJECT_SCHEMA_SQL + + +# ─── init_project / add_to_project ─────────────────────────────────── + +_VALID_EXAMPLE = '>example_sample\nATGAAAGCTTTTGGCCCC\n' +_MULTI_RECORD_EXAMPLE = '>a\nATGAAA\n>b\nGGGCCC\n' +_INVALID_FASTA = 'not a fasta at all\nno header line\n' + + +def _tiny_genbank(tmp_path: Path) -> Path: + gb = tmp_path / 'tiny.gb' + write_genbank(gb, [ + { + 'id': 'tiny_ref', + 'accession': 'tiny_ref', + 'sequence': TINY_REF_SEQ, + 'features': [{'feature': 'gag', 'protein': 'Gag', 'start': 1, 'end': 87, 'strand': '+'}], + } + ]) + return gb + + +def _rules_tsv(tmp_path: Path) -> Path: + tsv = tmp_path / 'rules.tsv' + tsv.write_text(textwrap.dedent("""\ + feature\treference_identifier\tposition\treference\tmutation\tantiviral\tphenotype + gag\ttiny_ref\t2\tK\tE\tDrugA\tresistant + """)) + return tsv + + +class TestInitProjectExample: + def test_init_stores_example_fasta_text(self, tmp_path: Path) -> None: + example = tmp_path / 'example.fasta' + example.write_text(_VALID_EXAMPLE) + db = tmp_path / 'proj.db' + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + example_fasta=example, + additional_info=False, + ) + + conn = open_project_db(db) + assert get_project_example_fasta(conn) == _VALID_EXAMPLE + conn.close() + + def test_init_without_example_leaves_empty(self, tmp_path: Path) -> None: + db = tmp_path / 'proj.db' + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + additional_info=False, + ) + + conn = open_project_db(db) + assert get_project_example_fasta(conn) is None + conn.close() + + def test_init_rejects_multi_record_example(self, tmp_path: Path) -> None: + example = tmp_path / 'example.fasta' + example.write_text(_MULTI_RECORD_EXAMPLE) + db = tmp_path / 'proj.db' + with pytest.raises(ValueError, match='single record'): + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + example_fasta=example, + additional_info=False, + ) + assert not db.exists() + + def test_init_rejects_invalid_fasta(self, tmp_path: Path) -> None: + example = tmp_path / 'example.fasta' + example.write_text(_INVALID_FASTA) + db = tmp_path / 'proj.db' + with pytest.raises(ValueError, match='(?i)example fasta'): + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + example_fasta=example, + additional_info=False, + ) + assert not db.exists() + + +class TestAddToProjectExample: + def test_add_overwrites_existing_example(self, tmp_path: Path) -> None: + example = tmp_path / 'example.fasta' + example.write_text(_VALID_EXAMPLE) + db = tmp_path / 'proj.db' + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + example_fasta=example, + additional_info=False, + ) + + new_example = tmp_path / 'example2.fasta' + new_example.write_text('>other\nGGGCCC\n') + add_to_project( + db_path=db, + rules_tsv=_rules_tsv(tmp_path), + example_fasta=new_example, + additional_info=False, + ) + + conn = open_project_db(db) + assert get_project_example_fasta(conn) == '>other\nGGGCCC\n' + conn.close() + + def test_add_clears_example_when_none(self, tmp_path: Path) -> None: + example = tmp_path / 'example.fasta' + example.write_text(_VALID_EXAMPLE) + db = tmp_path / 'proj.db' + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + example_fasta=example, + additional_info=False, + ) + + add_to_project( + db_path=db, + rules_tsv=_rules_tsv(tmp_path), + example_fasta=None, + clear_example=True, + additional_info=False, + ) + + conn = open_project_db(db) + assert get_project_example_fasta(conn) is None + conn.close() + + def test_add_rejects_multi_record_example(self, tmp_path: Path) -> None: + db = tmp_path / 'proj.db' + init_project( + db_path=db, + name='test', + genbank_paths=[_tiny_genbank(tmp_path)], + rules_tsv=_rules_tsv(tmp_path), + additional_info=False, + ) + example = tmp_path / 'example.fasta' + example.write_text(_MULTI_RECORD_EXAMPLE) + with pytest.raises(ValueError, match='single record'): + add_to_project( + db_path=db, + rules_tsv=_rules_tsv(tmp_path), + example_fasta=example, + additional_info=False, + ) diff --git a/tests/test_features.py b/tests/test_features.py new file mode 100644 index 00000000..9efb8236 --- /dev/null +++ b/tests/test_features.py @@ -0,0 +1,215 @@ +""" +Tests for feature and reference loading. + +Covers: respro/db/features.py +- load_features_for_reference() +- load_feature_segments_by_feature_id() +- _is_ncbi_protein_accession() +- _resolve_ncbi_protein_url() +""" + +from __future__ import annotations + +import sqlite3 +from urllib.parse import urlparse + +import pytest + +from respro.db.features import ( + _is_ncbi_protein_accession, + _resolve_ncbi_protein_url, + load_feature_segments_by_feature_id, + load_features_for_reference, +) +from respro.db.models import FeatureRecord, FeatureSegment + + +@pytest.fixture +def in_memory_db(): + """Create in-memory database with minimal schema.""" + conn = sqlite3.connect(':memory:') + conn.row_factory = sqlite3.Row + # Create only the tables we need for these tests + conn.execute(''' + CREATE TABLE reference ( + id INTEGER PRIMARY KEY AUTOINCREMENT, + project_id INTEGER NOT NULL, + name TEXT NOT NULL, + accession TEXT, + organism TEXT, + taxonomy TEXT, + length INTEGER NOT NULL + ) + ''') + conn.execute(''' + CREATE TABLE feature ( + id INTEGER PRIMARY KEY AUTOINCREMENT, + reference_id INTEGER NOT NULL, + name TEXT NOT NULL, + protein TEXT, + start INTEGER NOT NULL, + end INTEGER NOT NULL, + strand TEXT NOT NULL, + codon_start INTEGER, + nt_sequence TEXT, + aa_sequence TEXT, + feature_type TEXT DEFAULT 'CDS', + parent_feature_name TEXT, + FOREIGN KEY (reference_id) REFERENCES reference(id) + ) + ''') + conn.execute(''' + CREATE TABLE feature_segment ( + feature_id INTEGER NOT NULL, + segment_index INTEGER NOT NULL, + start INTEGER NOT NULL, + end INTEGER NOT NULL, + PRIMARY KEY (feature_id, segment_index), + FOREIGN KEY (feature_id) REFERENCES feature(id) + ) + ''') + conn.commit() + return conn + + +class TestLoadFeatureSegmentsByFeatureId: + """Tests for load_feature_segments_by_feature_id().""" + + def test_returns_empty_dict_for_empty_list(self, in_memory_db): + """Should return empty dict for empty feature ID list.""" + result = load_feature_segments_by_feature_id(in_memory_db, []) + assert result == {} + + def test_returns_empty_dict_for_nonexistent_ids(self, in_memory_db): + """Should return empty dict for nonexistent IDs.""" + result = load_feature_segments_by_feature_id(in_memory_db, [999]) + assert result == {} + + def test_loads_segments_for_existing_features(self, in_memory_db): + """Should load segments for existing features.""" + # Insert test data + in_memory_db.execute( + 'INSERT INTO reference (id, project_id, name, length) VALUES (1, 1, "test", 1000)' + ) + in_memory_db.execute( + 'INSERT INTO feature (id, reference_id, name, start, end, strand) VALUES (1, 1, "gene1", 100, 500, "+")' + ) + in_memory_db.execute( + 'INSERT INTO feature_segment (feature_id, segment_index, start, end) VALUES (1, 0, 100, 300)' + ) + in_memory_db.execute( + 'INSERT INTO feature_segment (feature_id, segment_index, start, end) VALUES (1, 1, 400, 500)' + ) + in_memory_db.commit() + + result = load_feature_segments_by_feature_id(in_memory_db, [1]) + assert 1 in result + assert len(result[1]) == 2 + assert isinstance(result[1][0], FeatureSegment) + assert result[1][0].segment_index == 0 + assert result[1][1].segment_index == 1 + + +class TestLoadFeaturesForReference: + """Tests for load_features_for_reference().""" + + def test_returns_empty_list_for_no_features(self, in_memory_db): + """Should return empty list when no features exist.""" + in_memory_db.execute( + 'INSERT INTO reference (id, project_id, name, length) VALUES (1, 1, "test", 1000)' + ) + in_memory_db.commit() + result = load_features_for_reference(in_memory_db, 1) + assert result == [] + + def test_loads_features_with_segments(self, in_memory_db): + """Should load features with their segments.""" + # Insert test data + in_memory_db.execute( + 'INSERT INTO reference (id, project_id, name, length) VALUES (1, 1, "test", 1000)' + ) + in_memory_db.execute( + 'INSERT INTO feature (id, reference_id, name, start, end, strand, nt_sequence, aa_sequence, feature_type) ' + 'VALUES (1, 1, "gene1", 100, 500, "+", "ATGC", "MK", "CDS")' + ) + in_memory_db.execute( + 'INSERT INTO feature_segment (feature_id, segment_index, start, end) VALUES (1, 0, 100, 300)' + ) + in_memory_db.commit() + + result = load_features_for_reference(in_memory_db, 1) + assert len(result) == 1 + assert isinstance(result[0], FeatureRecord) + assert result[0].name == 'gene1' + assert len(result[0].segments) == 1 + + +class TestIsNcbiProteinAccession: + """Tests for _is_ncbi_protein_accession().""" + + def test_matches_xxx_format(self): + """Should match XXX12345.1 format.""" + assert _is_ncbi_protein_accession('AAA12345.1') is True + assert _is_ncbi_protein_accession('NP_123456.2') is True + + def test_matches_xx_underscore_format(self): + """Should match XX_123456.1 format.""" + assert _is_ncbi_protein_accession('YP_009137097.1') is True + assert _is_ncbi_protein_accession('NP_123456.2') is True + + def test_matches_xxxx_format(self): + """Should match XXXX12345678.1 format.""" + assert _is_ncbi_protein_accession('KAFS00000001.1') is True + + def test_rejects_no_version(self): + """Should reject accessions without version suffix.""" + assert _is_ncbi_protein_accession('AAA12345') is False + + def test_rejects_invalid_format(self): + """Should reject invalid formats.""" + assert _is_ncbi_protein_accession('invalid') is False + assert _is_ncbi_protein_accession('12345') is False + + def test_case_insensitive(self): + """Should be case-insensitive.""" + assert _is_ncbi_protein_accession('aaa12345.1') is True + assert _is_ncbi_protein_accession('yp_009137097.1') is True + + +class TestResolveNcbiProteinUrl: + """Tests for _resolve_ncbi_protein_url().""" + + def test_returns_url_for_valid_accession(self): + """Should return URL for valid accession.""" + cache: dict[str, str] = {} + url = _resolve_ncbi_protein_url('AAA12345.1', cache) + parsed = urlparse(url) + assert parsed.hostname == 'www.ncbi.nlm.nih.gov' + assert 'AAA12345.1' in url + + def test_returns_empty_for_invalid_accession(self): + """Should return empty string for invalid accession.""" + cache: dict[str, str] = {} + url = _resolve_ncbi_protein_url('invalid', cache) + assert url == '' + + def test_returns_empty_for_empty_input(self): + """Should return empty string for empty input.""" + cache: dict[str, str] = {} + url = _resolve_ncbi_protein_url('', cache) + assert url == '' + + def test_uses_cache(self): + """Should use cached URL.""" + cache: dict[str, str] = {} + url1 = _resolve_ncbi_protein_url('AAA12345.1', cache) + assert 'AAA12345.1' in cache + url2 = _resolve_ncbi_protein_url('AAA12345.1', cache) + assert url1 == url2 + assert len(cache) == 1 # Only one entry + + def test_strips_whitespace(self): + """Should strip whitespace from protein_id.""" + cache: dict[str, str] = {} + url = _resolve_ncbi_protein_url(' AAA12345.1 ', cache) + assert 'AAA12345.1' in url diff --git a/tests/test_formula_parsing.py b/tests/test_formula_parsing.py new file mode 100644 index 00000000..a0ea3fef --- /dev/null +++ b/tests/test_formula_parsing.py @@ -0,0 +1,416 @@ +""" +Tests for boolean formula expression parsing. + +Covers: respro/db/_rules_formula.py +- _tokenize_formula_expression() +- _parse_formula_expression() +- _canonicalize_formula_ast() +- _formula_ast_to_string() +- _formula_sort_key() +""" + +from __future__ import annotations + +import pytest + +from respro.db._rules_formula import ( + _canonicalize_formula_ast, + _formula_ast_to_string, + _formula_sort_key, + _parse_formula_expression, + _tokenize_formula_expression, +) + + +class TestTokenizeFormulaExpression: + """Tests for _tokenize_formula_expression().""" + + def test_tokenizes_simple_atom(self): + """Should tokenize a single atomic rule ID.""" + tokens = _tokenize_formula_expression('mut_1') + assert tokens == ['mut_1'] + + def test_tokenizes_multiple_atoms(self): + """Should tokenize multiple atomic rule IDs.""" + tokens = _tokenize_formula_expression('mut_1 mut_2 mut_3') + assert tokens == ['mut_1', 'mut_2', 'mut_3'] + + def test_tokenizes_operators_uppercase(self): + """Should convert operators to uppercase.""" + tokens = _tokenize_formula_expression('mut1 and mut2 or mut3') + assert tokens == ['mut1', 'AND', 'mut2', 'OR', 'mut3'] + + def test_tokenizes_parentheses(self): + """Should tokenize parentheses.""" + tokens = _tokenize_formula_expression('(mut1 OR mut2) AND mut3') + assert tokens == ['(', 'mut1', 'OR', 'mut2', ')', 'AND', 'mut3'] + + def test_tokenizes_not_operator(self): + """Should tokenize NOT operator.""" + tokens = _tokenize_formula_expression('NOT mut1') + assert tokens == ['NOT', 'mut1'] + + def test_tokenizes_xor_operator(self): + """Should tokenize XOR operator.""" + tokens = _tokenize_formula_expression('mut1 XOR mut2') + assert tokens == ['mut1', 'XOR', 'mut2'] + + def test_handles_complex_expression(self): + """Should tokenize complex expressions.""" + tokens = _tokenize_formula_expression('(A AND B) OR (NOT C)') + assert tokens == ['(', 'A', 'AND', 'B', ')', 'OR', '(', 'NOT', 'C', ')'] + + def test_rejects_empty_expression(self): + """Should raise error for empty expression.""" + with pytest.raises(ValueError, match='empty expression'): + _tokenize_formula_expression('') + + def test_rejects_whitespace_only(self): + """Should raise error for whitespace-only expression.""" + with pytest.raises(ValueError, match='empty expression'): + _tokenize_formula_expression(' ') + + def test_rejects_unsupported_characters(self): + """Should raise error for unsupported characters.""" + with pytest.raises(ValueError, match='unsupported characters'): + _tokenize_formula_expression('mut1 + mut2') + + def test_rejects_special_characters(self): + """Should raise error for special characters.""" + with pytest.raises(ValueError, match='unsupported characters'): + _tokenize_formula_expression('mut1 & mut2') + + def test_allows_underscores_in_ids(self): + """Should allow underscores in rule IDs.""" + tokens = _tokenize_formula_expression('rule_123 AND rule_456') + assert tokens == ['rule_123', 'AND', 'rule_456'] + + def test_allows_dashes_in_ids(self): + """Should allow dashes in rule IDs.""" + tokens = _tokenize_formula_expression('rule-1 AND rule-2') + assert tokens == ['rule-1', 'AND', 'rule-2'] + + def test_allows_colons_in_ids(self): + """Should allow colons in rule IDs.""" + tokens = _tokenize_formula_expression('gene:A AND gene:B') + assert tokens == ['gene:A', 'AND', 'gene:B'] + + def test_allows_dots_in_ids(self): + """Should allow dots in rule IDs.""" + tokens = _tokenize_formula_expression('v1.0 AND v2.0') + assert tokens == ['v1.0', 'AND', 'v2.0'] + + def test_allows_numeric_ids(self): + """Should allow purely numeric rule IDs.""" + tokens = _tokenize_formula_expression('123 AND 456') + assert tokens == ['123', 'AND', '456'] + + def test_case_insensitive_operators(self): + """Should handle operators in any case.""" + tokens = _tokenize_formula_expression('mut1 And mut2 Or mut3') + assert tokens == ['mut1', 'AND', 'mut2', 'OR', 'mut3'] + + +class TestParseFormulaExpression: + """Tests for _parse_formula_expression().""" + + def test_parses_single_atom(self): + """Should parse a single atomic rule ID.""" + canonical_str, referenced_ids = _parse_formula_expression('mut_1') + assert canonical_str == 'mut_1' + assert referenced_ids == ['mut_1'] + + def test_parses_and_expression(self): + """Should parse AND expression.""" + canonical_str, referenced_ids = _parse_formula_expression('A AND B') + assert canonical_str == '(A AND B)' + assert referenced_ids == ['A', 'B'] + + def test_parses_or_expression(self): + """Should parse OR expression.""" + canonical_str, referenced_ids = _parse_formula_expression('A OR B') + assert canonical_str == '(A OR B)' + assert referenced_ids == ['A', 'B'] + + def test_parses_xor_expression(self): + """Should parse XOR expression.""" + canonical_str, referenced_ids = _parse_formula_expression('A XOR B') + assert canonical_str == '(A XOR B)' + assert referenced_ids == ['A', 'B'] + + def test_parses_not_expression(self): + """Should parse NOT expression.""" + canonical_str, referenced_ids = _parse_formula_expression('NOT A') + assert canonical_str == '(NOT A)' + assert referenced_ids == ['A'] + + def test_operator_precedence_not_over_and(self): + """NOT should have higher precedence than AND.""" + canonical_str, referenced_ids = _parse_formula_expression('NOT A AND B') + # Structure: AND with (NOT A) and B - canonical form sorts operands + assert 'NOT A' in canonical_str + assert 'B' in canonical_str + assert 'AND' in canonical_str + assert referenced_ids == ['A', 'B'] + + def test_operator_precedence_and_over_or(self): + """AND should have higher precedence than OR.""" + canonical_str, referenced_ids = _parse_formula_expression('A AND B OR C') + # Structure: OR of (A AND B) and C + assert 'A AND B' in canonical_str + assert 'OR' in canonical_str + assert referenced_ids == ['A', 'B', 'C'] + + def test_operator_precedence_xor_between_and_or(self): + """XOR should have higher precedence than OR, lower than AND.""" + canonical_str, referenced_ids = _parse_formula_expression('A AND B XOR C OR D') + assert 'A AND B' in canonical_str + assert 'XOR' in canonical_str + assert 'OR' in canonical_str + assert referenced_ids == ['A', 'B', 'C', 'D'] + + def test_left_associative_and(self): + """AND should be left-associative.""" + canonical_str, referenced_ids = _parse_formula_expression('A AND B AND C') + assert canonical_str == '(A AND B AND C)' + assert referenced_ids == ['A', 'B', 'C'] + + def test_left_associative_or(self): + """OR should be left-associative.""" + canonical_str, referenced_ids = _parse_formula_expression('A OR B OR C') + assert canonical_str == '(A OR B OR C)' + assert referenced_ids == ['A', 'B', 'C'] + + def test_left_associative_xor(self): + """XOR should be left-associative.""" + canonical_str, referenced_ids = _parse_formula_expression('A XOR B XOR C') + assert canonical_str == '(A XOR B XOR C)' + assert referenced_ids == ['A', 'B', 'C'] + + def test_right_associative_not(self): + """NOT should be right-associative.""" + canonical_str, referenced_ids = _parse_formula_expression('NOT NOT A') + assert canonical_str == '(NOT (NOT A))' + assert referenced_ids == ['A'] + + def test_parentheses_override_precedence(self): + """Parentheses should override default precedence.""" + canonical_str, referenced_ids = _parse_formula_expression('A AND (B OR C)') + assert canonical_str == '(A AND (B OR C))' + assert referenced_ids == ['A', 'B', 'C'] + + def test_nested_parentheses(self): + """Should handle nested parentheses.""" + canonical_str, referenced_ids = _parse_formula_expression('(A OR (B AND C)) XOR D') + assert 'B AND C' in canonical_str + assert 'OR' in canonical_str + assert 'XOR' in canonical_str + assert referenced_ids == ['A', 'B', 'C', 'D'] + + def test_complex_expression(self): + """Should parse complex nested expressions.""" + expr = '(A AND B) OR (NOT C AND D)' + canonical_str, referenced_ids = _parse_formula_expression(expr) + assert 'A AND B' in canonical_str + assert 'NOT C' in canonical_str + assert 'D' in canonical_str + assert 'OR' in canonical_str + assert referenced_ids == ['A', 'B', 'C', 'D'] + + def test_not_in_parentheses(self): + """Should handle NOT inside parentheses.""" + canonical_str, referenced_ids = _parse_formula_expression('(NOT A) AND B') + assert 'NOT A' in canonical_str + assert 'B' in canonical_str + assert 'AND' in canonical_str + assert referenced_ids == ['A', 'B'] + + def test_multiple_not_operators(self): + """Should handle multiple NOT operators.""" + canonical_str, referenced_ids = _parse_formula_expression('NOT A AND NOT B') + assert 'NOT A' in canonical_str + assert 'NOT B' in canonical_str + assert 'AND' in canonical_str + assert referenced_ids == ['A', 'B'] + + def test_rejects_unbalanced_opening_paren(self): + """Should raise error for unbalanced opening parenthesis.""" + with pytest.raises(ValueError, match='unexpected end of expression'): + _parse_formula_expression('(A AND B') + + def test_rejects_unbalanced_closing_paren(self): + """Should raise error for unbalanced closing parenthesis.""" + with pytest.raises(ValueError, match=r"unexpected token '\)'"): + _parse_formula_expression('A AND B)') + + def test_rejects_operator_after_paren(self): + """Should raise error for operator after opening paren.""" + with pytest.raises(ValueError, match="unexpected token 'AND'"): + _parse_formula_expression('(AND A)') + + def test_rejects_double_operator(self): + """Should raise error for consecutive operators.""" + with pytest.raises(ValueError): + _parse_formula_expression('A AND OR B') + + def test_tracks_duplicate_rule_ids(self): + """Should track duplicate rule IDs in referenced_ids.""" + canonical_str, referenced_ids = _parse_formula_expression('A AND A') + assert referenced_ids == ['A', 'A'] + + def test_canonical_form_flattens_and(self): + """Should flatten nested AND operators in canonical form.""" + canonical_str, _ = _parse_formula_expression('A AND (B AND C)') + assert canonical_str == '(A AND B AND C)' + + def test_canonical_form_sorts_operands(self): + """Should sort operands for deterministic canonical form.""" + canonical_str1, _ = _parse_formula_expression('B AND A') + canonical_str2, _ = _parse_formula_expression('A AND B') + assert canonical_str1 == canonical_str2 + + def test_canonical_form_distinguishes_not(self): + """Canonical form should distinguish NOT from positive atoms.""" + canonical_str1, _ = _parse_formula_expression('A AND B') + canonical_str2, _ = _parse_formula_expression('NOT A AND B') + assert canonical_str1 != canonical_str2 + + +class TestCanonicalizeFormulaAst: + """Tests for _canonicalize_formula_ast().""" + + def test_atom_unchanged(self): + """Should return atom nodes unchanged.""" + atom = ('ATOM', 'mut_1') + result = _canonicalize_formula_ast(atom) + assert result == atom + + def test_not_wrapped(self): + """Should wrap NOT nodes but canonicalize child.""" + not_node = ('NOT', ('ATOM', 'A')) + result = _canonicalize_formula_ast(not_node) + assert result == ('NOT', ('ATOM', 'A')) + + def test_flattens_nested_and(self): + """Should flatten nested AND nodes.""" + nested = ('AND', [('AND', [('ATOM', 'A'), ('ATOM', 'B')]), ('ATOM', 'C')]) + result = _canonicalize_formula_ast(nested) + assert result == ('AND', [('ATOM', 'A'), ('ATOM', 'B'), ('ATOM', 'C')]) + + def test_flattens_nested_or(self): + """Should flatten nested OR nodes.""" + nested = ('OR', [('OR', [('ATOM', 'A'), ('ATOM', 'B')]), ('ATOM', 'C')]) + result = _canonicalize_formula_ast(nested) + assert result == ('OR', [('ATOM', 'A'), ('ATOM', 'B'), ('ATOM', 'C')]) + + def test_sorts_and_operands(self): + """Should sort AND operands for deterministic ordering.""" + unsorted = ('AND', [('ATOM', 'B'), ('ATOM', 'A')]) + result = _canonicalize_formula_ast(unsorted) + assert result == ('AND', [('ATOM', 'A'), ('ATOM', 'B')]) + + def test_sorts_or_operands(self): + """Should sort OR operands for deterministic ordering.""" + unsorted = ('OR', [('ATOM', 'C'), ('ATOM', 'A'), ('ATOM', 'B')]) + result = _canonicalize_formula_ast(unsorted) + assert result == ('OR', [('ATOM', 'A'), ('ATOM', 'B'), ('ATOM', 'C')]) + + def test_sorts_xor_operands(self): + """Should sort XOR operands for deterministic ordering.""" + unsorted = ('XOR', [('ATOM', 'Z'), ('ATOM', 'A')]) + result = _canonicalize_formula_ast(unsorted) + assert result == ('XOR', [('ATOM', 'A'), ('ATOM', 'Z')]) + + def test_keeps_not_after_positive(self): + """Should sort NOT nodes after positive atoms.""" + mixed = ('AND', [('NOT', ('ATOM', 'B')), ('ATOM', 'A')]) + result = _canonicalize_formula_ast(mixed) + assert result == ('AND', [('ATOM', 'A'), ('NOT', ('ATOM', 'B'))]) + + def test_handles_complex_nested_structure(self): + """Should canonicalize complex nested structures.""" + ast = ( + 'OR', + [ + ('AND', [('NOT', ('ATOM', 'C')), ('ATOM', 'A')]), + ('ATOM', 'B'), + ], + ) + result = _canonicalize_formula_ast(ast) + assert result[0] == 'OR' + assert len(result[1]) == 2 + + +class TestFormulaAstToString: + """Tests for _formula_ast_to_string().""" + + def test_atom_to_string(self): + """Should convert atom to string.""" + result = _formula_ast_to_string(('ATOM', 'mut_1')) + assert result == 'mut_1' + + def test_not_to_string(self): + """Should convert NOT to string.""" + result = _formula_ast_to_string(('NOT', ('ATOM', 'A'))) + assert result == '(NOT A)' + + def test_and_to_string(self): + """Should convert AND to string.""" + result = _formula_ast_to_string(('AND', [('ATOM', 'A'), ('ATOM', 'B')])) + assert result == '(A AND B)' + + def test_or_to_string(self): + """Should convert OR to string.""" + result = _formula_ast_to_string(('OR', [('ATOM', 'A'), ('ATOM', 'B')])) + assert result == '(A OR B)' + + def test_xor_to_string(self): + """Should convert XOR to string.""" + result = _formula_ast_to_string(('XOR', [('ATOM', 'A'), ('ATOM', 'B')])) + assert result == '(A XOR B)' + + def test_nested_expression_to_string(self): + """Should convert nested expressions to string.""" + ast = ('OR', [('AND', [('ATOM', 'A'), ('ATOM', 'B')]), ('ATOM', 'C')]) + result = _formula_ast_to_string(ast) + assert result == '((A AND B) OR C)' + + def test_triple_and_to_string(self): + """Should convert triple AND to string.""" + ast = ('AND', [('ATOM', 'A'), ('ATOM', 'B'), ('ATOM', 'C')]) + result = _formula_ast_to_string(ast) + assert result == '(A AND B AND C)' + + def test_not_of_and_to_string(self): + """Should convert NOT of AND to string.""" + ast = ('NOT', ('AND', [('ATOM', 'A'), ('ATOM', 'B')])) + result = _formula_ast_to_string(ast) + assert result == '(NOT (A AND B))' + + +class TestFormulaSortKey: + """Tests for _formula_sort_key().""" + + def test_atom_comes_before_not(self): + """Atoms should sort before NOT nodes.""" + atom_key = _formula_sort_key(('ATOM', 'A')) + not_key = _formula_sort_key(('NOT', ('ATOM', 'B'))) + assert atom_key < not_key + + def test_compound_same_rank_as_atom(self): + """AND/OR/XOR should have same rank (1).""" + and_key = _formula_sort_key(('AND', [('ATOM', 'A')])) + or_key = _formula_sort_key(('OR', [('ATOM', 'A')])) + xor_key = _formula_sort_key(('XOR', [('ATOM', 'A')])) + assert and_key[0] == or_key[0] == xor_key[0] == 1 + + def test_not_has_highest_rank(self): + """NOT should have highest rank (2).""" + not_key = _formula_sort_key(('NOT', ('ATOM', 'A'))) + assert not_key[0] == 2 + + def test_sorts_by_string_representation(self): + """Should use string representation as secondary sort key.""" + a_key = _formula_sort_key(('ATOM', 'A')) + b_key = _formula_sort_key(('ATOM', 'B')) + assert a_key < b_key diff --git a/tests/test_init_project.py b/tests/test_init_project.py index 362ae34d..c699e7a4 100644 --- a/tests/test_init_project.py +++ b/tests/test_init_project.py @@ -73,6 +73,7 @@ def _rule( reference: str, mutation: str = 'E', reference_identifier: str = '', + antiviral: str = 'DrugX', ) -> dict[str, str]: """Build a minimal TSV row dict for a single resistance rule.""" row = { @@ -80,7 +81,7 @@ def _rule( 'position': str(position), 'reference': reference, 'mutation': mutation, - 'antiviral': 'DrugX', + 'antiviral': antiviral, } if reference_identifier: row['reference_identifier'] = reference_identifier @@ -253,6 +254,40 @@ def test_collects_all_mismatch_keys(self, caplog) -> None: assert len(mismatch_keys) == 3 assert any('3' in r.message and 'mismatch' in r.message.lower() for r in caplog.records) + def test_emits_each_mismatch_only_once_per_position(self, caplog) -> None: + import logging + + # Same feature/position/ref_aa repeated across multiple antivirals — + # the mismatch should be logged exactly once and the key returned once. + rows = [ + _rule('gX', position=1, reference='Z', antiviral='DrugA'), + _rule('gX', position=1, reference='Z', antiviral='DrugB'), + _rule('gX', position=1, reference='Z', antiviral='DrugC'), + ] + with caplog.at_level(logging.WARNING, logger='respro.db.rules_import'): + mismatch_keys = _validate_reference_amino_acids(rows, _features('gX', _AA_SEQ), coord_base=1) + assert mismatch_keys == {('gX', '1', '', 'Z')} + mismatch_messages = [ + r.message for r in caplog.records if 'mismatch' in r.message.lower() + ] + assert len(mismatch_messages) == 1 + assert mismatch_messages[0].count("pos 1 (1-based)") == 1 + + def test_emits_out_of_range_only_once_per_position(self, caplog) -> None: + import logging + + rows = [ + _rule('gX', position=10, reference='M', antiviral='DrugA'), + _rule('gX', position=10, reference='M', antiviral='DrugB'), + ] + with caplog.at_level(logging.WARNING, logger='respro.db.rules_import'): + _validate_reference_amino_acids(rows, _features('gX', _AA_SEQ), coord_base=1) + out_of_range_messages = [ + r.message for r in caplog.records if 'out of range' in r.message + ] + assert len(out_of_range_messages) == 1 + assert out_of_range_messages[0].count("pos 10 (1-based)") == 1 + def test_skips_rows_without_ref_aa(self) -> None: # Row has no 'reference'/'ref_aa' key — must not raise rows = [ diff --git a/tests/test_intron_classification.py b/tests/test_intron_classification.py new file mode 100644 index 00000000..2855cb32 --- /dev/null +++ b/tests/test_intron_classification.py @@ -0,0 +1,741 @@ +""" +Tests for intron-aware CIGAR classification in spliced-feature alignment. + +When a user supplies an unspliced (whole-genome) query, minimap2 reports the +intron between two CDS exons as a single large ``I`` op in the normalized CIGAR +(CDS=reference, genome=query). These tests pin the behaviour that classifies +such ``I`` ops as introns by exon-junction CDS position (not by length) so that +downstream consumers never see the intron as a coding insertion. +""" + +from __future__ import annotations + +import pytest + +from respro.core.alignment import ( + classify_introns, + exon_junction_cds_offsets, + match_query_to_features, +) +from respro.core.fasta_to_vcf import fasta_to_vcf +from respro.core.vcf_remap import _build_query_to_cds_map, remap_variants +from respro.db.models import FeatureRecord, FeatureSegment, IntronInterval, VariantCall + +# ────────────────────────────────────────────────────────────────────── +# exon_junction_cds_offsets +# ────────────────────────────────────────────────────────────────────── + +class TestExonJunctionCdsOffsets: + def test_single_segment_feature_has_no_junctions(self) -> None: + """A feature with one segment (or no segments) has no exon junctions.""" + from respro.db.models import FeatureRecord + + feature = FeatureRecord( + id=1, reference_id=1, name='g', protein='', + start=0, end=30, strand='+', codon_start=0, + ) + assert exon_junction_cds_offsets(feature) == [] + + def test_two_segment_forward_strand(self) -> None: + """Forward strand: junction offset = length of first coding segment.""" + from respro.db.models import FeatureRecord, FeatureSegment + + feature = FeatureRecord( + id=1, reference_id=1, name='g', protein='', + start=0, end=18, strand='+', codon_start=0, + segments=( + FeatureSegment(segment_index=0, start=0, end=6), + FeatureSegment(segment_index=1, start=12, end=18), + ), + ) + # First exon is 6 nt (0..6), so the junction is at CDS offset 6. + assert exon_junction_cds_offsets(feature) == [6] + + def test_two_segment_reverse_strand(self) -> None: + """Reverse strand: junction offsets use genomic (segment_index) order, + NOT _coding_segments reversed order, because the normalized CIGAR for + '-' strand is reversed and walks genomic 5'->3' order. With unequal + segment lengths the two orderings yield different junction offsets, so + this test discriminates the chosen convention.""" + from respro.db.models import FeatureRecord, FeatureSegment + + feature = FeatureRecord( + id=1, reference_id=1, name='g', protein='', + start=0, end=26, strand='-', codon_start=0, + segments=( + FeatureSegment(segment_index=0, start=0, end=6), # 6 nt (genomic 5') + FeatureSegment(segment_index=1, start=18, end=26), # 8 nt (genomic 3') + ), + ) + # Genomic order (segment_index): first exon 6 nt → junction at CDS offset 6. + # Coding order (_coding_segments reversed): first exon 8 nt → junction 8. + # The normalized '-' CIGAR walks genomic order, so junction must be 6. + assert exon_junction_cds_offsets(feature) == [6] + + def test_three_segment_forward_strand(self) -> None: + """Three exons of lengths 4, 5, 6 → junctions at 4 and 9.""" + from respro.db.models import FeatureRecord, FeatureSegment + + feature = FeatureRecord( + id=1, reference_id=1, name='g', protein='', + start=0, end=40, strand='+', codon_start=0, + segments=( + FeatureSegment(segment_index=0, start=0, end=4), + FeatureSegment(segment_index=1, start=10, end=15), + FeatureSegment(segment_index=2, start=30, end=36), + ), + ) + assert exon_junction_cds_offsets(feature) == [4, 9] + + +# ────────────────────────────────────────────────────────────────────── +# classify_introns +# ────────────────────────────────────────────────────────────────────── + +class TestClassifyIntrons: + def test_no_junctions_returns_cigar_unchanged(self) -> None: + """A single-exon feature: no I op is ever classified as an intron.""" + cigar, introns = classify_introns( + cigar='10M3I5M', junction_offsets=[], tolerance=5, + ) + assert cigar == '10M3I5M' + assert introns == [] + + def test_intron_at_junction_is_removed(self) -> None: + """A 3902-nt I op exactly at the exon junction is classified as intron.""" + # 6M (exon 1) + 3902I (intron) + 7M (exon 2), junction at CDS offset 6. + cigar, introns = classify_introns( + cigar='6M3902I7M', junction_offsets=[6], tolerance=5, + ) + # Adjacent M runs merge after intron removal → clean 13M. + assert cigar == '13M' + assert len(introns) == 1 + assert introns[0].cds_junction_pos == 6 + assert introns[0].length == 3902 + + def test_intron_within_tolerance_is_removed(self) -> None: + """An I op whose CDS position is within ±tolerance of the junction is + classified as an intron (mappy may place the junction a few nt off).""" + # Junction at 100; the I op sits at CDS position 102 (2 nt inside tolerance 5). + # 102M consumed before I then 3902I then 50M. + cigar, introns = classify_introns( + cigar='102M3902I50M', junction_offsets=[100], tolerance=5, + ) + assert cigar == '152M' + assert len(introns) == 1 + + def test_coding_insertion_inside_exon_is_kept(self) -> None: + """A small I op NOT near any junction is a real coding insertion.""" + # Junction at 100; 3I at CDS position 5 (deep in exon 1, |5-100|=95 > 5). + # 6I at CDS position 108 (|108-100|=8 > 5) → also kept (beyond tolerance). + cigar, introns = classify_introns( + cigar='5M3I103M6I50M', junction_offsets=[100], tolerance=5, + ) + assert '3I' in cigar + assert '6I' in cigar + assert introns == [] + + def test_coding_insertion_far_from_junction_is_kept(self) -> None: + """A 3-nt insertion deep inside exon 1 (far from junction) is kept.""" + # Junction at 100; 3I at CDS position 5. + cigar, introns = classify_introns( + cigar='5M3I95M2000I50M', junction_offsets=[100], tolerance=5, + ) + assert '3I' in cigar + assert '2000I' not in cigar + assert len(introns) == 1 + assert introns[0].length == 2000 + + def test_multiple_introns_multiple_junctions(self) -> None: + """Three exons → two junctions; two intron I ops both removed.""" + # Exons 40, 50, 60 nt. Junctions at 40 and 90. + # 40M 100I 50M 200I 60M + cigar, introns = classify_introns( + cigar='40M100I50M200I60M', junction_offsets=[40, 90], tolerance=5, + ) + assert cigar == '150M' + assert len(introns) == 2 + assert introns[0].cds_junction_pos == 40 + assert introns[0].length == 100 + assert introns[1].cds_junction_pos == 90 + assert introns[1].length == 200 + + def test_i_op_beyond_tolerance_is_kept(self) -> None: + """An I op at a position more than `tolerance` from any junction is kept.""" + # Junction at 6; I op at CDS position 20 (tolerance 5 → 20-6=14 > 5). + cigar, introns = classify_introns( + cigar='20M50I10M', junction_offsets=[6], tolerance=5, + ) + assert cigar == '20M50I10M' + assert introns == [] + + def test_deletion_ops_advance_cds_position(self) -> None: + """D ops consume CDS positions; an intron I op after a D must still be + detected at the right junction offset.""" + # 30M 1D 10M (consumes 41 CDS nt: 30+1+10) then 3902I then 50M, junction at 41. + cigar, introns = classify_introns( + cigar='30M1D10M3902I50M', junction_offsets=[41], tolerance=5, + ) + assert '3902I' not in cigar + assert len(introns) == 1 + assert introns[0].cds_junction_pos == 41 + + def test_small_insertion_within_tolerance_is_kept(self) -> None: + """A real coding insertion whose CDS position is within ±tolerance of a + junction must NOT be classified as an intron if its length is ≤ tolerance. + This is the edge-case guard: a 3-nt insertion near an exon boundary is a + real variant, not an intron.""" + # Junction at 62; 3I at CDS position 60 (|60-62|=2 ≤ tolerance 5), and + # length 3 ≤ tolerance 5 → kept as a real insertion. + cigar, introns = classify_introns( + cigar='60M3I2M200I50M', junction_offsets=[62], tolerance=5, + ) + assert '3I' in cigar, '3-nt insertion near junction must be kept (length ≤ tolerance)' + # The 200I is at CDS pos 62 (60M + 2M = 62; the 3I does not advance cds_pos). + # Its query_start is 65 (60M + 3I consumes 3 query nt + 2M = 65). + assert introns == [IntronInterval(cds_junction_pos=62, query_start=65, length=200)] + + def test_large_insertion_within_tolerance_is_classified_as_intron(self) -> None: + """An I op within tolerance AND longer than tolerance is an intron.""" + # Junction at 62; 200I at CDS position 60 (within tolerance 5), length 200 > 5. + cigar, introns = classify_introns( + cigar='60M200I2M50M', junction_offsets=[62], tolerance=5, + ) + assert '200I' not in cigar + assert len(introns) == 1 + assert introns[0].length == 200 + + def test_insertion_exactly_at_tolerance_boundary_is_kept(self) -> None: + """An I op of length exactly equal to tolerance is kept (guard is strict >). + Only length strictly greater than tolerance is classified as an intron.""" + cigar, introns = classify_introns( + cigar='62M5I50M', junction_offsets=[62], tolerance=5, + ) + assert '5I' in cigar + assert introns == [] + + def test_insertion_one_above_tolerance_boundary_is_intron(self) -> None: + """An I op of length tolerance+1 at the junction is classified as an intron.""" + cigar, introns = classify_introns( + cigar='62M6I50M', junction_offsets=[62], tolerance=5, + ) + assert '6I' not in cigar + assert len(introns) == 1 + assert introns[0].length == 6 + + +# ────────────────────────────────────────────────────────────────────── +# match_query_to_features — end-to-end intron classification +# ────────────────────────────────────────────────────────────────────── + +def _spliced_feature() -> FeatureRecord: + """A two-exon forward-strand CDS: exon1 0..30, intron 30..60, exon2 60..90. + + The stored nt_sequence is the spliced CDS (exon1+exon2 = 60 nt). + """ + exon1 = 'ATGAAAGCTTTTGGCCCCAAATTTGGGCCC' # 30 nt + exon2 = 'AAGCTTTTGGCCCCAAATTTGGGCCCAAT'[:30] # 30 nt (pad if needed) + # Ensure exactly 30 nt each. + exon1 = (exon1 + 'A' * 30)[:30] + exon2 = (exon2 + 'A' * 30)[:30] + return FeatureRecord( + id=1, reference_id=1, name='spl', protein='Spl', + start=0, end=90, strand='+', codon_start=0, + nt_sequence=exon1 + exon2, + segments=( + FeatureSegment(segment_index=0, start=0, end=30), + FeatureSegment(segment_index=1, start=60, end=90), + ), + ) + + +class TestMatchQueryToFeaturesSpliced: + def test_unspliced_query_classifies_intron_and_reports_exon_identity(self) -> None: + """An unspliced query (exon1+intron+exon2) yields an exon-only CIGAR, + one intron interval, identity≈1.0 (not crashed by the intron), and + cds_coverage=1.0 over the exon span.""" + feature = _spliced_feature() + intron = 'N' * 30 # 30-nt intron between the two exons + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + + # The intron I op must be removed from the exon-only CIGAR. + assert 'I' not in m.cigar or '30I' not in m.cigar + # One intron interval recorded at the exon junction (CDS offset 30). + assert len(m.intron_intervals) == 1 + assert m.intron_intervals[0].cds_junction_pos == 30 + assert m.intron_intervals[0].length == 30 + # Identity reflects exons only (perfect match → ~1.0), not crashed by intron. + assert m.identity > 0.95 + # CDS coverage is 1.0 over the exon span (60 nt). + assert m.cds_coverage == 1.0 + + def test_spliced_query_no_intron_op_classified(self) -> None: + """A spliced query (exons only, no intron) leaves the CIGAR with no + intron intervals — the happy path is unchanged.""" + feature = _spliced_feature() + query = feature.nt_sequence # already spliced + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + assert m.intron_intervals == () + assert m.identity > 0.95 + + def test_reverse_strand_spliced_feature_classifies_intron(self) -> None: + """A '-' strand two-exon CDS aligned against an unspliced query must + classify the intron at the correct exon junction. This is the HCMV + UL89 scenario: the normalized CIGAR is reversed for '-' strand, so + junction offsets must be in genomic (segment_index) order, not coding + order, to match the CIGAR's walking order.""" + from Bio.Seq import Seq + + exon1 = 'ATGAAAGCTTTTGGCCCCAAATTTGGGCCC' # 30 nt, genomic 5' segment + exon2 = 'AAGCTTTTGGCCCCAAATTTGGGCCCAAT'[:30] + exon2 = (exon2 + 'A' * 30)[:30] + # '-' strand: the spliced CDS is revcomp(exon1+exon2) read 5'->3'. + spliced_cds = str(Seq(exon1 + exon2).reverse_complement()) + feature = FeatureRecord( + id=2, reference_id=1, name='spl_minus', protein='SplM', + start=0, end=90, strand='-', codon_start=0, + nt_sequence=spliced_cds, + segments=( + FeatureSegment(segment_index=0, start=0, end=30), + FeatureSegment(segment_index=1, start=60, end=90), + ), + ) + # Unspliced query: exon1 + intron + exon2 (genomic 5'->3'). + intron = 'N' * 30 + query = exon1 + intron + exon2 + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + # The intron must be classified (one interval), not left as a giant I op. + assert len(m.intron_intervals) == 1, ( + f'expected 1 intron interval for - strand spliced feature, got {m.intron_intervals}' + ) + assert m.intron_intervals[0].length == 30 + # Exon-only CIGAR must not contain the 30I intron op. + assert '30I' not in m.cigar + # Identity reflects exons only. + assert m.identity > 0.95 + assert m.cds_coverage == 1.0 + + +# ────────────────────────────────────────────────────────────────────── +# fasta_to_vcf — intron must not become a giant insertion variant +# ────────────────────────────────────────────────────────────────────── + +class TestFastaToVcfIntron: + def test_unspliced_query_emits_no_intron_insertion_variant(self) -> None: + """The prove-it test for the real bug: an unspliced query (exon1 + + intron + exon2) must NOT emit a giant ~intron-length frameshift + insertion variant. Only genuine exon SNPs (if any) are emitted.""" + feature = _spliced_feature() + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + variants, gaps = fasta_to_vcf(query, matches) + + # No variant may have an ALT longer than a few bases (no 30-nt intron insertion). + intron_insertions = [v for v in variants if len(v.alt) - len(v.ref) >= 30] + assert intron_insertions == [], ( + f'Expected no intron-length insertion variant, got {intron_insertions}' + ) + # A perfect unspliced query (exons match exactly) yields zero variants. + assert variants == [] + + def test_unspliced_query_reports_exon_only_identity_and_coverage(self) -> None: + """For a perfect unspliced query, identity≈1.0 and cds_coverage=1.0 + (computed over exons, not the genomic span including the intron).""" + feature = _spliced_feature() + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + assert m.identity == pytest.approx(1.0, abs=1e-6) + assert m.cds_coverage == pytest.approx(1.0, abs=1e-6) + assert m.intron_intervals, 'spliced alignment should carry intron intervals' + + def test_spliced_query_happy_path_unchanged(self) -> None: + """A spliced query (exons only, no intron) follows the unchanged happy + path: no intron I op is present, so nothing is classified as an intron + and identity/coverage reflect a perfect match.""" + feature = _spliced_feature() + query = feature.nt_sequence # spliced CDS, no intron + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + assert m.intron_intervals == () + assert m.identity == pytest.approx(1.0, abs=1e-6) + variants, _ = fasta_to_vcf(query, matches) + assert variants == [] + + def test_unspliced_query_with_exon_snp_emits_only_the_snp(self) -> None: + """A SNP in exon 2 of an unspliced query is emitted; the intron is not.""" + feature = _spliced_feature() + intron = 'N' * 30 + # Introduce one SNP at exon 2 position 45 (CDS offset 45). + exon2_snp = ( + feature.nt_sequence[:30] + + intron + + feature.nt_sequence[30:45] + + ('T' if feature.nt_sequence[45] != 'T' else 'A') + + feature.nt_sequence[46:] + ) + matches = match_query_to_features(exon2_snp, [feature]) + assert len(matches) == 1 + variants, _ = fasta_to_vcf(exon2_snp, matches) + + # Exactly one variant (the exon-2 SNP), no intron insertion. + intron_insertions = [v for v in variants if len(v.alt) - len(v.ref) >= 30] + assert intron_insertions == [] + assert len(variants) == 1 + + def test_real_coding_insertion_inside_exon_is_still_emitted(self) -> None: + """A genuine 3-nt coding insertion inside exon 1 (not at the junction) + is still emitted as a variant — intron classification must not suppress + real in-exon insertions. Tested on a spliced query (no intron present) + so the insertion is the only indel; the intron-classification path is + a no-op here (no junction I op), proving it never suppresses real I ops.""" + feature = _spliced_feature() + # Insert 3 nt 'AAA' at CDS offset 10 (inside exon 1, far from junction 30). + query = feature.nt_sequence[:10] + 'AAA' + feature.nt_sequence[10:] + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + # No intron classified (the 3I is at offset 10, not near junction 30). + assert matches[0].intron_intervals == () + variants, _ = fasta_to_vcf(query, matches) + + # The 3-nt insertion should appear as an insertion variant (alt longer than ref). + insertions = [v for v in variants if len(v.alt) > len(v.ref)] + assert any(len(v.alt) - len(v.ref) == 3 for v in insertions), ( + f'Expected a 3-nt coding insertion variant, got {variants}' + ) + + +# ────────────────────────────────────────────────────────────────────── +# VCF remap — intron-aware query-to-CDS coordinate map +# ────────────────────────────────────────────────────────────────────── + +class TestVcfRemapIntron: + def test_exon2_variant_in_unspliced_query_remaps_to_second_segment(self) -> None: + """A VCF variant in exon 2 of an unspliced query (exon1+intron+exon2) + must remap to the correct second-segment genomic position, not be + shifted by the intron length or land in the intron query span.""" + feature = _spliced_feature() # exon1 0..30, exon2 60..90 (genomic) + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + assert len(m.intron_intervals) == 1 + + # A variant at the first base of exon 2 in the query. + # Query layout: exon1 [0..30) + intron [30..60) + exon2 [60..90). + # First exon-2 base is at query position 60. + exon2_first_base = query[60] + var = VariantCall( + chrom='c', pos=60, ref=exon2_first_base, alt='T', allele_freq=1.0, depth=50, + ) + remapped, warnings = remap_variants([var], [m], query) + + assert not warnings + assert len(remapped) == 1 + # Exon 2 starts at genomic position 60; first exon-2 base → genomic 60. + assert remapped[0].pos == 60 + + def test_variant_inside_intron_query_span_is_excluded(self) -> None: + """A VCF variant whose query position falls inside the intron span is + not remapped (it is non-coding).""" + feature = _spliced_feature() + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + + matches = match_query_to_features(query, [feature]) + m = matches[0] + # Variant at query position 45 (middle of the intron span 30..60). + var = VariantCall(chrom='c', pos=45, ref='N', alt='A', allele_freq=1.0, depth=50) + remapped, warnings = remap_variants([var], [m], query) + + assert remapped == [] + + def test_build_query_to_cds_map_skips_intron_query_span(self) -> None: + """The query-to-CDS map must skip intron query positions and offset + exon-2 CDS positions past the intron span.""" + feature = _spliced_feature() + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + + matches = match_query_to_features(query, [feature]) + m = matches[0] + + q2c = _build_query_to_cds_map( + m.cigar, m.query_start, m.query_end, m.strand, len(query), m.cds_start, + m.intron_intervals, + ) + # Exon 1: query 0..30 → CDS 0..30. + assert q2c[0] == 0 + assert q2c[29] == 29 + # Intron query positions 30..59 must NOT be in the map. + for qpos in range(30, 60): + assert qpos not in q2c, f'intron query pos {qpos} should not map to a CDS pos' + # Exon 2: query 60 → CDS 30 (first base of exon 2 in the spliced CDS). + assert q2c[60] == 30 + assert q2c[89] == 59 + + def test_reverse_strand_exon_variants_remap_correctly(self) -> None: + """'-' strand spliced feature against an unspliced query: exon-1 and + exon-2 VCF variants must remap to the correct genomic segments, and + intron-span variants must be excluded. This is the HCMV UL89 scenario + (the motivating bug).""" + from Bio.Seq import Seq + + exon1 = 'ATGAAAGCTTTTGGCCCCAAATTTGGGCCC' # 30 nt, genomic seg0 [0..30) + exon2 = 'AAGCTTTTGGCCCCAAATTTGGGCCCAAT'[:30] + exon2 = (exon2 + 'A' * 30)[:30] + spliced_cds = str(Seq(exon1 + exon2).reverse_complement()) + feature = FeatureRecord( + id=3, reference_id=1, name='splm', protein='SplM', + start=0, end=90, strand='-', codon_start=0, + nt_sequence=spliced_cds, + segments=( + FeatureSegment(segment_index=0, start=0, end=30), + FeatureSegment(segment_index=1, start=60, end=90), + ), + ) + intron = 'N' * 30 + query = exon1 + intron + exon2 # unspliced, genomic 5'->3' + + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + assert len(m.intron_intervals) == 1, ( + f'expected 1 intron interval, got {m.intron_intervals}' + ) + + q2c = _build_query_to_cds_map( + m.cigar, m.query_start, m.query_end, m.strand, len(query), m.cds_start, + m.intron_intervals, + ) + + # Intron query positions [30..60) must NOT be in the map. + for qpos in range(30, 60): + assert qpos not in q2c, f'intron query pos {qpos} should not map to a CDS pos' + + # For '-' strand the spliced CDS is revcomp(exon1+exon2); coding 5'->3' + # order is seg1 (genomic exon2) then seg0 (genomic exon1). The normalized + # CIGAR walks genomic 5'->3' (seg0 then seg1), so: + # query [0..30) = genomic exon1 (seg0) → CDS offsets [30..60) + # query [60..90) = genomic exon2 (seg1) → CDS offsets [0..30) + exon1_cds = {q2c[q] for q in range(0, 30) if q in q2c} + exon2_cds = {q2c[q] for q in range(60, 90) if q in q2c} + assert exon1_cds, 'exon-1 query positions must map to CDS positions' + assert exon2_cds, 'exon-2 query positions must map to CDS positions' + assert all(30 <= c < 60 for c in exon1_cds), ( + f'exon-1 (genomic seg0) CDS offsets must be in [30,60), got {exon1_cds}' + ) + assert all(0 <= c < 30 for c in exon2_cds), ( + f'exon-2 (genomic seg1) CDS offsets must be in [0,30), got {exon2_cds}' + ) + + # End-to-end remap: an exon-2 variant (query pos 60, first base of seg1) + # maps to CDS offset 29 (last coding base of seg1) → genomic 60, inside + # segment 1 [60..90). A regression that mis-remaps to seg0 must fail. + var_exon2 = VariantCall( + chrom='c', pos=60, ref=query[60], alt='T', allele_freq=1.0, depth=50, + ) + remapped, warnings = remap_variants([var_exon2], [m], query) + assert not warnings, f'unexpected warnings: {warnings}' + assert len(remapped) == 1, f'expected 1 remapped variant, got {remapped}' + genomic = remapped[0].pos + assert genomic == 60, ( + f'exon-2 variant should remap to genomic 60 (seg1), got {genomic}' + ) + + # An intron-span variant is excluded. + var_intron = VariantCall( + chrom='c', pos=45, ref='N', alt='A', allele_freq=1.0, depth=50, + ) + remapped_intron, _ = remap_variants([var_intron], [m], query) + assert remapped_intron == [] + + def test_adjacent_m_runs_merged_in_exon_only_cigar(self) -> None: + """After removing an intron I op, adjacent M runs are merged.""" + cigar, introns = classify_introns( + cigar='6M3902I7M', junction_offsets=[6], tolerance=5, + ) + # 6M + 7M should merge to 13M for a clean exon-only CIGAR. + assert cigar == '13M' + + +# ────────────────────────────────────────────────────────────────────── +# Cache round-trip for intron_intervals +# ────────────────────────────────────────────────────────────────────── + +class TestCacheIntronRoundTrip: + """Storing and reloading a FeatureMatch with intron_intervals must preserve + the intron intervals and the exon-only CIGAR exactly.""" + + def test_store_load_preserves_intron_intervals(self, tmp_path) -> None: # type: ignore[no-untyped-def] + from respro.db.cache import load_cached_mappings, sequence_checksum, store_mappings + from respro.db.schema import create_schema, open_project_db + + db_path = tmp_path / 'test_project.db' + conn = create_schema(db_path) + conn.execute( + 'INSERT INTO project (name, schema_version) VALUES (?, ?)', ('Test', 1), + ) + conn.execute( + 'INSERT INTO reference (project_id, name, length) VALUES (?, ?, ?)', + (1, 'ref1', 90), + ) + feature = _spliced_feature() + conn.execute( + 'INSERT INTO feature (reference_id, name, protein, start, end, strand, ' + 'codon_start, nt_sequence, aa_sequence) VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?)', + ( + feature.reference_id, feature.name, feature.protein, + feature.start, feature.end, feature.strand, + feature.codon_start, feature.nt_sequence, feature.aa_sequence, + ), + ) + for seg in feature.segments: + conn.execute( + 'INSERT INTO feature_segment (feature_id, segment_index, start, end) ' + 'VALUES (?, ?, ?, ?)', + (feature.id, seg.segment_index, seg.start, seg.end), + ) + conn.commit() + + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + assert matches[0].intron_intervals, 'spliced alignment should carry intron intervals' + + chk = sequence_checksum(query) + store_mappings(conn, 'unspliced_ref', query, chk, matches) + conn.close() + + # Re-open via open_project_db to exercise the optional-column migration. + conn2 = open_project_db(db_path) + loaded = load_cached_mappings(conn2, chk) + conn2.close() + + assert loaded is not None + assert len(loaded) == 1 + m = loaded[0] + assert m.cigar == matches[0].cigar + assert m.intron_intervals == matches[0].intron_intervals + assert len(m.intron_intervals) == 1 + iv = m.intron_intervals[0] + assert iv.cds_junction_pos == 30 + assert iv.query_start == 30 + assert iv.length == 30 + + def test_legacy_rows_without_intron_column_backfill_as_empty(self, tmp_path) -> None: # type: ignore[no-untyped-def] + """A row inserted before the intron_intervals column existed must load + with an empty intron_intervals tuple after migration.""" + from respro.db.cache import load_cached_mappings, sequence_checksum, store_mappings + from respro.db.schema import create_schema, open_project_db + + db_path = tmp_path / 'test_project.db' + conn = create_schema(db_path) + conn.execute( + 'INSERT INTO project (name, schema_version) VALUES (?, ?)', ('Test', 1), + ) + conn.execute( + 'INSERT INTO reference (project_id, name, length) VALUES (?, ?, ?)', + (1, 'ref1', 30), + ) + conn.execute( + 'INSERT INTO feature (reference_id, name, protein, start, end, strand, ' + 'nt_sequence, aa_sequence) VALUES (?, ?, ?, ?, ?, ?, ?, ?)', + (1, 'gag', 'Gag', 0, 30, '+', 'ATGAAAGCTTTTGGCCCCAAATTTGGGCCC', 'MKAFGPKFGP'), + ) + conn.commit() + + feature = FeatureRecord( + id=1, reference_id=1, name='gag', protein='Gag', + start=0, end=30, strand='+', codon_start=0, + nt_sequence='ATGAAAGCTTTTGGCCCCAAATTTGGGCCC', aa_sequence='MKAFGPKFGP', + segments=tuple(), + ) + from respro.db.models import FeatureMatch + legacy_match = FeatureMatch( + feature=feature, identity=1.0, cds_coverage=1.0, query_coverage=1.0, + cds_start=0, query_start=0, query_end=30, strand='+', cigar='30M', + intron_intervals=tuple(), + ) + query = 'ATGAAAGCTTTTGGCCCCAAATTTGGGCCC' + chk = sequence_checksum(query) + store_mappings(conn, 'legacy_ref', query, chk, [legacy_match]) + conn.close() + + conn2 = open_project_db(db_path) + loaded = load_cached_mappings(conn2, chk) + conn2.close() + + assert loaded is not None + assert len(loaded) == 1 + assert loaded[0].intron_intervals == tuple() + + +# ────────────────────────────────────────────────────────────────────── +# Alignment visualization for spliced features (verification) +# ────────────────────────────────────────────────────────────────────── + +class TestAlignmentVisualizationSpliced: + """The per-codon alignment snippet must render exon1→exon2 with no intron + gap row, because the exon-only CIGAR + intron-stripped region leave only + coding bases.""" + + def test_snippet_has_no_intron_gap_and_exon2_aligns(self) -> None: + from respro.report.alignment_visualization import build_feature_alignments + + feature = _spliced_feature() + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + m = matches[0] + assert m.intron_intervals, 'spliced alignment should carry intron intervals' + + alns = build_feature_alignments(query, [m]) + assert 'spl' in alns + aln = alns['spl'] + # No gap row spanning ~30 intron bases: the aligned strings must be + # exactly the spliced CDS length (60) with no long '-' run. + assert len(aln.aligned_ref) == len(aln.aligned_query) + # The longest run of '-' in aligned_ref must be far shorter than the + # 30-nt intron (at most a few nt for any unaligned CDS tail). + import re as _re + longest_gap = max((len(run) for run in _re.findall(r'-+', aln.aligned_ref)), default=0) + assert longest_gap < 10, f'intron gap leaked into snippet: {longest_gap}-nt gap' + # Exon 2 query bases must align to exon 2 CDS bases (first exon-2 base + # is the last char of the spliced CDS region after stripping). + # The aligned query must contain exon2 sequence contiguously. + exon2 = feature.nt_sequence[30:] + assert exon2 in aln.aligned_query.replace('-', ''), 'exon 2 not aligned contiguously' + + def test_identity_is_exon_only_not_genomic(self) -> None: + """match.identity carried into the report must reflect exon-only identity + (~1.0 for a perfect exons match), not the genomic 33.96%.""" + feature = _spliced_feature() + intron = 'N' * 30 + query = feature.nt_sequence[:30] + intron + feature.nt_sequence[30:] + matches = match_query_to_features(query, [feature]) + assert len(matches) == 1 + assert matches[0].identity > 0.99 diff --git a/tests/test_maintained_bootstrap.py b/tests/test_maintained_bootstrap.py index 49f0d059..8daa5337 100644 --- a/tests/test_maintained_bootstrap.py +++ b/tests/test_maintained_bootstrap.py @@ -46,6 +46,7 @@ def _make_files(tmp_path: Path, *, tsv_checksum: str = _MANAGER_TSV_CHECKSUM) -> 'rules': rules, 'metadata': metadata, 'formula_rules': None, + 'example': None, 'genbank': [gb], } @@ -80,7 +81,7 @@ def _install_init_checksum_writer(remote_checksum: str) -> None: """ def fake_init(*, db_path, name, genbank_paths, rules_tsv, formula_rules_tsv, metadata_json, - overwrite=False, additional_info=True): + overwrite=False, additional_info=True, example_fasta=None): # Build a minimal valid DB then stamp the remote checksum. from respro.db.schema import create_schema @@ -262,7 +263,7 @@ def test_downloads_missing_database_and_initializes_atomically(self, tmp_path: P init_calls: list[Path] = [] def init_side_effect(*, db_path, name, genbank_paths, rules_tsv, formula_rules_tsv, - metadata_json, overwrite=False, additional_info=True): + metadata_json, overwrite=False, additional_info=True, example_fasta=None): from respro.db.schema import create_schema conn = create_schema(db_path) diff --git a/tests/test_maintained_db.py b/tests/test_maintained_db.py index 7bbc248a..ef4edb53 100644 --- a/tests/test_maintained_db.py +++ b/tests/test_maintained_db.py @@ -84,6 +84,16 @@ def _manifest_with_formula_rules_path(formula_rules_path: str) -> dict: manifest['databases'][0]['formula_rules_path'] = formula_rules_path return manifest + +def _manifest_with_example_fasta_path(example_fasta_path: str) -> dict: + """Return a manifest copy with a custom example_fasta_path for the hsv entry.""" + manifest = json.loads(json.dumps(_MANIFEST)) + manifest['databases'][0]['example_fasta_path'] = example_fasta_path + return manifest + + +_EXAMPLE_FASTA_CONTENT = b'>example_sample\nATGAAAGCTTTTGGCCCC\n' + _RULES_TSV_CONTENT = ( 'feature\tposition\tref_aa\tmut_aa\tphenotype\treference_identifier\n' 'UL23\t168\tA\tT\tresistant\tX04770\n' @@ -189,6 +199,54 @@ def test_raises_on_http_error(self) -> None: list_output_files('hsv_daehne_jaki') +# ── list_output_files: example_fasta_path ───────────────────────────────────── + +class TestListOutputFilesExampleFasta: + def test_includes_example_fasta_when_path_present(self) -> None: + manifest = _manifest_with_example_fasta_path( + 'databases/hsv_daehne_jaki/output/example.fasta' + ) + with patch('urllib.request.urlopen', return_value=_json_mock(manifest)): + result = list_output_files('hsv_daehne_jaki') + names = {entry['name'] for entry in result} + assert 'example.fasta' in names + + def test_omits_example_fasta_when_path_absent(self) -> None: + with patch('urllib.request.urlopen', return_value=_json_mock(_MANIFEST)): + result = list_output_files('hsv_daehne_jaki') + names = {entry['name'] for entry in result} + assert 'example.fasta' not in names + + def test_omits_example_fasta_when_path_empty(self) -> None: + manifest = _manifest_with_example_fasta_path('') + with patch('urllib.request.urlopen', return_value=_json_mock(manifest)): + result = list_output_files('hsv_daehne_jaki') + names = {entry['name'] for entry in result} + assert 'example.fasta' not in names + + +# ── _validate_manifest_entry: example_fasta_path ────────────────────────────── + +class TestValidateManifestEntryExampleFasta: + def test_accepts_absent_example_fasta_path(self) -> None: + # The default _MANIFEST has no example_fasta_path; listing must not raise. + with patch('urllib.request.urlopen', return_value=_json_mock(_MANIFEST)): + result = list_maintained_databases() + assert 'hsv_daehne_jaki' in result + + def test_accepts_empty_example_fasta_path(self) -> None: + manifest = _manifest_with_example_fasta_path('') + with patch('urllib.request.urlopen', return_value=_json_mock(manifest)): + result = list_maintained_databases() + assert 'hsv_daehne_jaki' in result + + def test_rejects_non_string_example_fasta_path(self) -> None: + manifest = _manifest_with_example_fasta_path(123) # type: ignore[arg-type] + with patch('urllib.request.urlopen', return_value=_json_mock(manifest)): + with pytest.raises(RuntimeError, match='example_fasta_path must be a string'): + list_maintained_databases() + + # ── _parse_reference_identifiers ───────────────────────────────────────────── class TestParseReferenceIdentifiers: @@ -404,6 +462,45 @@ def se(*args, **kwargs): # noqa: ANN001, ANN202 assert result['formula_rules'] is None + def test_downloads_example_fasta_when_path_present(self, tmp_path: Path) -> None: + rules_bytes = _RULES_TSV_CONTENT.encode() + meta_bytes = json.dumps(_METADATA).encode() + manifest = _manifest_with_example_fasta_path( + 'databases/hsv_daehne_jaki/output/example.fasta' + ) + + responses = [ + _json_mock(manifest), # manifest fetch + _bytes_mock(rules_bytes), # rules.tsv + _bytes_mock(meta_bytes), # metadata.json + _bytes_mock(b'feature\tformula\n'), # formula-rules.tsv + _bytes_mock(_EXAMPLE_FASTA_CONTENT), # example.fasta + _bytes_mock(_GENBANK_CONTENT), # X04770.gb + _bytes_mock(_GENBANK_CONTENT), # X04771.gb + ] + call_count = [0] + + def se(*args, **kwargs): # noqa: ANN001, ANN202 + idx = call_count[0] + call_count[0] += 1 + return responses[min(idx, len(responses) - 1)] + + with patch('urllib.request.urlopen', side_effect=se): + result = download_database_files('hsv_daehne_jaki', tmp_path) + + assert isinstance(result['example'], Path) + assert result['example'].read_bytes() == _EXAMPLE_FASTA_CONTENT + + def test_example_is_none_when_path_absent(self, tmp_path: Path) -> None: + rules_bytes = _RULES_TSV_CONTENT.encode() + meta_bytes = json.dumps(_METADATA).encode() + + se = self._make_download_side_effect(rules_bytes, meta_bytes) + with patch('urllib.request.urlopen', side_effect=se): + result = download_database_files('hsv_daehne_jaki', tmp_path) + + assert result['example'] is None + # ── CLI: databases --list ──────────────────────────────────────────────── @@ -466,6 +563,7 @@ def test_download_calls_init_project(self, tmp_path: Path) -> None: 'rules': tmp_path / 'rules.tsv', 'metadata': tmp_path / 'metadata.json', 'formula_rules': None, + 'example': None, 'genbank': [fake_gb], } @@ -493,6 +591,7 @@ def test_download_fails_when_no_genbank(self, tmp_path: Path) -> None: 'rules': tmp_path / 'rules.tsv', 'metadata': tmp_path / 'metadata.json', 'formula_rules': None, + 'example': None, 'genbank': [], } with patch('respro.cli.maintained_db.download_database_files', return_value=fake_files): @@ -524,6 +623,7 @@ def test_download_outputs_success_path(self, tmp_path: Path) -> None: 'rules': tmp_path / 'rules.tsv', 'metadata': tmp_path / 'metadata.json', 'formula_rules': None, + 'example': None, 'genbank': [fake_gb], } with ( @@ -545,6 +645,7 @@ def test_download_uses_name_db_default_when_output_omitted(self, tmp_path: Path) 'rules': tmp_path / 'rules.tsv', 'metadata': tmp_path / 'metadata.json', 'formula_rules': None, + 'example': None, 'genbank': [fake_gb], } @@ -568,6 +669,7 @@ def test_download_treats_output_without_suffix_as_directory(self, tmp_path: Path 'rules': tmp_path / 'rules.tsv', 'metadata': tmp_path / 'metadata.json', 'formula_rules': None, + 'example': None, 'genbank': [fake_gb], } @@ -592,6 +694,7 @@ def test_download_disables_additional_info_with_flag(self, tmp_path: Path) -> No 'rules': tmp_path / 'rules.tsv', 'metadata': tmp_path / 'metadata.json', 'formula_rules': None, + 'example': None, 'genbank': [fake_gb], } diff --git a/tests/test_profile_fasta.py b/tests/test_profile_fasta.py index 4e4aa220..ee08eb32 100644 --- a/tests/test_profile_fasta.py +++ b/tests/test_profile_fasta.py @@ -1451,3 +1451,102 @@ def test_divergent_sequence_real_codon_deletion_detected(self) -> None: f'3-nt deletion, got {len(deletions)}: ' f'{[(a.codon_pos, a.ref_aa, a.alt_aa) for a in deletions]}' ) + + +class TestFastaExampleFlag: + """End-to-end CLI tests for ``respro fasta --example`` (Feature: db-example-fasta).""" + + @pytest.fixture() + def fasta_db_with_example(self, fasta_db: Path) -> Path: + """Return a fasta_db that has a single-record example consensus FASTA stored.""" + mutant = 'ATG' + 'GAA' + TINY_REF_SEQ[6:] # K→E at codon 1 → resistance hit + example_text = f'>stored_example\n{mutant}\n' + conn = open_project_db(fasta_db) + conn.execute('UPDATE project SET example_fasta = ? WHERE id = 1', (example_text,)) + conn.commit() + conn.close() + return fasta_db + + def test_example_profiles_stored_fasta( + self, fasta_db_with_example: Path, tmp_path: Path, + ) -> None: + """``--example`` profiles the stored example and produces the same hit as ``--fasta``.""" + output_dir = tmp_path / 'out' + result = CliRunner().invoke(app, [ + 'fasta', + '--project', str(fasta_db_with_example), + '--example', + '--output', str(output_dir), + ]) + + assert result.exit_code == 0, result.output + assert '1 total database hits' in result.output or '1 hit' in result.output + + def test_example_matches_explicit_fasta_result( + self, fasta_db_with_example: Path, tmp_path: Path, + ) -> None: + """``--example`` and ``--fasta`` with the same sequence yield identical hit counts.""" + mutant = 'ATG' + 'GAA' + TINY_REF_SEQ[6:] + fasta_path = tmp_path / 'explicit.fasta' + fasta_path.write_text(f'>explicit\n{mutant}\n') + + out_example = tmp_path / 'out_example' + out_explicit = tmp_path / 'out_explicit' + runner = CliRunner() + res_example = runner.invoke(app, [ + 'fasta', '--project', str(fasta_db_with_example), + '--example', '--output', str(out_example), + ]) + res_explicit = runner.invoke(app, [ + 'fasta', '--project', str(fasta_db_with_example), + '--fasta', str(fasta_path), '--output', str(out_explicit), + ]) + + assert res_example.exit_code == 0, res_example.output + assert res_explicit.exit_code == 0, res_explicit.output + assert '1 total database hits' in res_example.output + assert '1 total database hits' in res_explicit.output + + def test_example_errors_when_db_has_no_example(self, fasta_db: Path, tmp_path: Path) -> None: + """``--example`` on a DB without a stored example fails with a clear message.""" + output_dir = tmp_path / 'out' + result = CliRunner().invoke(app, [ + 'fasta', + '--project', str(fasta_db), + '--example', + '--output', str(output_dir), + ]) + + assert result.exit_code != 0 + clean = _strip_ansi(result.output) + assert 'example' in clean.lower() + + def test_example_and_fasta_are_mutually_exclusive( + self, fasta_db_with_example: Path, tmp_path: Path, + ) -> None: + """Providing both ``--example`` and ``--fasta`` is an error.""" + fasta_path = tmp_path / 'extra.fasta' + fasta_path.write_text(f'>x\n{TINY_REF_SEQ}\n') + output_dir = tmp_path / 'out' + result = CliRunner().invoke(app, [ + 'fasta', + '--project', str(fasta_db_with_example), + '--fasta', str(fasta_path), + '--example', + '--output', str(output_dir), + ]) + + assert result.exit_code != 0 + + def test_neither_example_nor_fasta_errors( + self, fasta_db_with_example: Path, tmp_path: Path, + ) -> None: + """Providing neither ``--example`` nor ``--fasta`` is an error.""" + output_dir = tmp_path / 'out' + result = CliRunner().invoke(app, [ + 'fasta', + '--project', str(fasta_db_with_example), + '--output', str(output_dir), + ]) + + assert result.exit_code != 0 diff --git a/tests/test_report_queries.py b/tests/test_report_queries.py new file mode 100644 index 00000000..89b8330e --- /dev/null +++ b/tests/test_report_queries.py @@ -0,0 +1,177 @@ +""" +Tests for report query helpers. + +Covers: respro/db/report_queries.py +- load_numeric_metric_thresholds() +- load_drug_class_map() +- load_drug_alias_map() +- _parse_numeric_value() +""" + +from __future__ import annotations + +import json +import sqlite3 + +import pytest + +from respro.db.report_queries import ( + load_drug_alias_map, + load_drug_class_map, + load_numeric_metric_thresholds, +) + + +@pytest.fixture +def memory_db(): + """Create in-memory database.""" + conn = sqlite3.connect(':memory:') + conn.row_factory = sqlite3.Row + return conn + + +class TestLoadNumericMetricThresholds: + """Tests for load_numeric_metric_thresholds().""" + + def test_returns_empty_dict_for_none_connection(self): + """Should return empty dict for None connection.""" + result = load_numeric_metric_thresholds(None) + assert result == {} + + def test_returns_dict_with_none_when_no_data(self, memory_db): + """Should return dict with None values when no data exists.""" + self._create_rule_tables(memory_db) + result = load_numeric_metric_thresholds(memory_db) + assert result == {'ic50': None, 'fold_ic50': None, 'score': None} + + def test_returns_none_for_insufficient_values(self, memory_db): + """Should return None when fewer than 2 values.""" + self._create_rule_tables(memory_db) + memory_db.execute('INSERT INTO resistance_rule (score) VALUES ("5.0")') + memory_db.commit() + result = load_numeric_metric_thresholds(memory_db) + assert result.get('score') is None + + def test_computes_mean_and_std(self, memory_db): + """Should compute mean and std for numeric fields.""" + self._create_rule_tables(memory_db) + memory_db.execute('INSERT INTO resistance_rule (score) VALUES ("10.0")') + memory_db.execute('INSERT INTO resistance_rule (score) VALUES ("20.0")') + memory_db.commit() + result = load_numeric_metric_thresholds(memory_db) + assert 'score' in result + assert result['score'] is not None + mean, std = result['score'] + assert mean == 15.0 + + def test_handles_null_values(self, memory_db): + """Should handle NULL values gracefully.""" + self._create_rule_tables(memory_db) + memory_db.execute('INSERT INTO resistance_rule (score) VALUES (NULL)') + memory_db.execute('INSERT INTO resistance_rule (score) VALUES ("10.0")') + memory_db.commit() + result = load_numeric_metric_thresholds(memory_db) + assert result.get('score') is None + + def test_parses_numeric_values_with_prefixes(self, memory_db): + """Should parse values with >, <, ~ prefixes.""" + self._create_rule_tables(memory_db) + memory_db.execute('INSERT INTO resistance_rule (score) VALUES (">10.0")') + memory_db.execute('INSERT INTO resistance_rule (score) VALUES ("~20.0")') + memory_db.commit() + result = load_numeric_metric_thresholds(memory_db) + assert result.get('score') is not None + + @staticmethod + def _create_rule_tables(conn): + """Create resistance_rule and resistance_formula_rule tables.""" + conn.execute(''' + CREATE TABLE resistance_rule ( + id INTEGER PRIMARY KEY, feature_id INTEGER, drug_id INTEGER, + external_id TEXT, position INTEGER, reference TEXT, mutation TEXT, + phenotype TEXT, clinical_phenotype TEXT, ic50 TEXT, fold_ic50 TEXT, + score TEXT, publication TEXT, source TEXT, comment TEXT + ) + ''') + conn.execute(''' + CREATE TABLE resistance_formula_rule ( + id INTEGER PRIMARY KEY, drug_id INTEGER, formula_id TEXT, + normalized_expression TEXT, phenotype TEXT, clinical_phenotype TEXT, + ic50 TEXT, fold_ic50 TEXT, score TEXT, source TEXT, comment TEXT + ) + ''') + conn.commit() + + +class TestLoadDrugClassMap: + """Tests for load_drug_class_map().""" + + def test_returns_empty_dict_for_none_connection(self): + """Should return empty dict for None connection.""" + result = load_drug_class_map(None) + assert result == {} + + def test_returns_empty_dict_when_no_algorithm(self, memory_db): + """Should return empty dict when no drug_groups algorithm.""" + memory_db.execute('CREATE TABLE interpretation_algorithm (id INTEGER PRIMARY KEY, algorithm_name TEXT, config_json TEXT)') + memory_db.commit() + result = load_drug_class_map(memory_db) + assert result == {} + + def test_builds_drug_to_class_mapping(self, memory_db): + """Should build drug to class mapping.""" + memory_db.execute('CREATE TABLE interpretation_algorithm (id INTEGER PRIMARY KEY, algorithm_name TEXT, config_json TEXT)') + config = { + 'groups': { + 'NRTI': ['tenofovir', 'emtricitabine'], + 'NNRTI': ['efavirenz', 'rilpivirine'], + } + } + memory_db.execute( + 'INSERT INTO interpretation_algorithm (algorithm_name, config_json) VALUES (?, ?)', + ('drug_groups', json.dumps(config)) + ) + memory_db.commit() + result = load_drug_class_map(memory_db) + assert result['tenofovir'] == 'NRTI' + assert result['emtricitabine'] == 'NRTI' + assert result['efavirenz'] == 'NNRTI' + + def test_lowercases_drug_names(self, memory_db): + """Should lowercase drug names.""" + memory_db.execute('CREATE TABLE interpretation_algorithm (id INTEGER PRIMARY KEY, algorithm_name TEXT, config_json TEXT)') + config = {'groups': {'Class': ['DrugA', 'DrugB']}} + memory_db.execute( + 'INSERT INTO interpretation_algorithm (algorithm_name, config_json) VALUES (?, ?)', + ('drug_groups', json.dumps(config)) + ) + memory_db.commit() + result = load_drug_class_map(memory_db) + assert 'druga' in result + assert 'drugb' in result + + +class TestLoadDrugAliasMap: + """Tests for load_drug_alias_map().""" + + def test_returns_empty_dict_for_none_connection(self): + """Should return empty dict for None connection.""" + result = load_drug_alias_map(None) + assert result == {} + + def test_returns_empty_dict_when_no_drugs(self, memory_db): + """Should return empty dict when no drugs in database.""" + memory_db.execute('CREATE TABLE drug (id INTEGER PRIMARY KEY, name TEXT, alias TEXT)') + memory_db.commit() + result = load_drug_alias_map(memory_db) + assert result == {} + + def test_builds_alias_mapping(self, memory_db): + """Should build drug alias mapping.""" + memory_db.execute('CREATE TABLE drug (id INTEGER PRIMARY KEY, name TEXT, alias TEXT)') + memory_db.execute('INSERT INTO drug (name, alias) VALUES ("Tenofovir", "TDF")') + memory_db.execute('INSERT INTO drug (name, alias) VALUES ("Emtricitabine", "FTC")') + memory_db.commit() + result = load_drug_alias_map(memory_db) + assert 'tenofovir' in result + assert result['tenofovir'] == 'TDF' diff --git a/tests/test_rules_alleles.py b/tests/test_rules_alleles.py new file mode 100644 index 00000000..4019abe4 --- /dev/null +++ b/tests/test_rules_alleles.py @@ -0,0 +1,306 @@ +""" +Tests for mutation and allele normalization in rule imports. + +Covers: respro/db/_rules_alleles.py +- _is_noop_mutation() +- _is_supported_mutation_token() +- _normalize_rule_alleles_for_storage() +- Anchorless deletion detection +- FSX/frameshift handling +""" + +from __future__ import annotations + +from respro.db._rules_alleles import ( + _RE_ANCHORLESS_DEL, + _RE_REWRITE_TOKEN, + _is_noop_mutation, + _is_supported_mutation_token, + _normalize_rule_alleles_for_storage, +) + + +class TestIsNoopMutation: + """Tests for _is_noop_mutation().""" + + def test_same_aa_returns_true(self): + """Should return True when reference and mutation are the same.""" + assert _is_noop_mutation('A', 'A') is True + assert _is_noop_mutation('G', 'G') is True + + def test_case_insensitive(self): + """Should be case-insensitive.""" + assert _is_noop_mutation('A', 'a') is True + assert _is_noop_mutation('g', 'G') is True + assert _is_noop_mutation('Leu', 'LEU') is True + + def test_different_aa_returns_false(self): + """Should return False when reference and mutation differ.""" + assert _is_noop_mutation('A', 'G') is False + assert _is_noop_mutation('K', 'R') is False + + def test_empty_strings(self): + """Should handle empty strings.""" + assert _is_noop_mutation('', '') is True + assert _is_noop_mutation('A', '') is False + + +class TestIsSupportedMutationToken: + """Tests for _is_supported_mutation_token().""" + + def test_standard_aa_mutations_supported(self): + """Should support standard amino acid mutations.""" + assert _is_supported_mutation_token('A') is True + assert _is_supported_mutation_token('K') is True + assert _is_supported_mutation_token('L') is True + assert _is_supported_mutation_token('ACD') is True + + def test_stop_codon_supported(self): + """Should support stop codon symbol.""" + assert _is_supported_mutation_token('*') is True + + def test_fsx_mutations_supported(self): + """Should support FSX (frameshift) mutations.""" + assert _is_supported_mutation_token('FSX') is True + assert _is_supported_mutation_token('KFSX') is True + assert _is_supported_mutation_token('RFSX') is True + + def test_ins_any_supported(self): + """Should support INS_any insertion.""" + assert _is_supported_mutation_token('INS_any') is True + + def test_any_not_supported(self): + """Should NOT support ANY token.""" + assert _is_supported_mutation_token('ANY') is False + + def test_invalid_aa_not_supported(self): + """Should reject invalid amino acid letters.""" + assert _is_supported_mutation_token('Z') is False + assert _is_supported_mutation_token('X') is False + assert _is_supported_mutation_token('BJZ') is False + + def test_numbers_not_supported(self): + """Should reject numeric tokens.""" + assert _is_supported_mutation_token('123') is False + + def test_mixed_case_handled(self): + """Should handle mixed case input (converted to uppercase internally).""" + assert _is_supported_mutation_token('abc') is False # B is not valid AA + assert _is_supported_mutation_token('AKC') is True # Valid AAs + assert _is_supported_mutation_token('KFSX') is True + + +class TestAnchorlessDeletionRegex: + """Tests for _RE_ANCHORLESS_DEL pattern.""" + + def test_matches_single_aa_del(self): + """Should match single AA deletion.""" + assert _RE_ANCHORLESS_DEL.match('Q35del') is not None + assert _RE_ANCHORLESS_DEL.match('A123del') is not None + + def test_matches_multiple_aa_del(self): + """Should match multiple AA deletion.""" + assert _RE_ANCHORLESS_DEL.match('DD676del') is not None + assert _RE_ANCHORLESS_DEL.match('KK34del') is not None + + def test_case_insensitive(self): + """Should match case-insensitively.""" + assert _RE_ANCHORLESS_DEL.match('q35DEL') is not None + assert _RE_ANCHORLESS_DEL.match('Q35DEL') is not None + assert _RE_ANCHORLESS_DEL.match('q35del') is not None + + def test_rejects_no_del_suffix(self): + """Should reject tokens without 'del' suffix.""" + assert _RE_ANCHORLESS_DEL.match('Q35') is None + assert _RE_ANCHORLESS_DEL.match('Q35deletion') is None + + def test_rejects_no_position(self): + """Should reject tokens without position.""" + assert _RE_ANCHORLESS_DEL.match('Qdel') is None + assert _RE_ANCHORLESS_DEL.match('del') is None + + def test_rejects_empty(self): + """Should reject empty string.""" + assert _RE_ANCHORLESS_DEL.match('') is None + + +class TestRewriteTokenRegex: + """Tests for _RE_REWRITE_TOKEN pattern.""" + + def test_matches_standard_mutation(self): + """Should match standard mutation format.""" + match = _RE_REWRITE_TOKEN.match('A123G') + assert match is not None + assert match.groups() == ('A', '123', 'G') + + def test_matches_stop_codon(self): + """Should match stop codon mutation.""" + match = _RE_REWRITE_TOKEN.match('K100*') + assert match is not None + assert match.groups() == ('K', '100', '*') + + def test_matches_from_stop(self): + """Should match mutation from stop codon.""" + match = _RE_REWRITE_TOKEN.match('*50K') + assert match is not None + assert match.groups() == ('*', '50', 'K') + + def test_matches_multi_aa(self): + """Should match multi-AA mutations.""" + match = _RE_REWRITE_TOKEN.match('DD100KK') + assert match is not None + assert match.groups() == ('DD', '100', 'KK') + + def test_rejects_no_position(self): + """Should reject mutations without position.""" + assert _RE_REWRITE_TOKEN.match('AG') is None + + def test_rejects_invalid_format(self): + """Should reject invalid formats.""" + assert _RE_REWRITE_TOKEN.match('123AG') is None + assert _RE_REWRITE_TOKEN.match('A123') is None + + +class TestNormalizeRuleAllelesForStorage: + """Tests for _normalize_rule_alleles_for_storage().""" + + def test_returns_tuple_for_valid_direct_aa(self): + """Should return normalized tuple for direct AA token.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='A', + mutation_raw='G', + position_0based=10, + context='test row', + errors=errors, + ) + assert result is not None + pos, ref, mut = result + assert pos == 10 + assert ref == 'A' + assert mut == 'G' + assert errors == [] + + def test_normalizes_stop_codon(self): + """Should normalize stop codon to *.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='K', + mutation_raw='STOP', + position_0based=50, + context='test', + errors=errors, + ) + assert result is not None + _, _, mut = result + assert mut == '*' + + def test_normalizes_frameshift_to_fsx(self): + """Should normalize frameshift to fsX format (case preserved from reference).""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='R', + mutation_raw='FS', + position_0based=100, + context='test', + errors=errors, + ) + assert result is not None + _, ref, mut = result + assert mut == 'RfsX' # Reference case preserved, 'fsX' suffix + + def test_handles_fsx_suffix(self): + """Should handle existing FSX suffix.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='K', + mutation_raw='KFSX', + position_0based=25, + context='test', + errors=errors, + ) + assert result is not None + _, _, mut = result + assert mut == 'KfsX' # Case preserved from input + + def test_returns_none_for_empty_reference(self): + """Should return None for empty reference.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='', + mutation_raw='G', + position_0based=10, + context='test', + errors=errors, + ) + assert result is None + + def test_returns_none_for_empty_mutation(self): + """Should return None for empty mutation.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='A', + mutation_raw='', + position_0based=10, + context='test', + errors=errors, + ) + assert result is None + + def test_appends_error_for_unrecognized_mutation(self): + """Should append error for unrecognized mutation format.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='A', + mutation_raw='invalid_mutation_123', + position_0based=10, + context='row 5', + errors=errors, + ) + assert result is None + assert len(errors) == 1 + assert 'row 5: unrecognised mutation' in errors[0] + + def test_strips_whitespace(self): + """Should strip whitespace from inputs.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa=' A ', + mutation_raw=' G ', + position_0based=10, + context='test', + errors=errors, + ) + assert result is not None + _, ref, mut = result + assert ref == 'A' + assert mut == 'G' + + def test_case_insensitive_reference(self): + """Should uppercase reference.""" + errors: list[str] = [] + result = _normalize_rule_alleles_for_storage( + reference_aa='a', + mutation_raw='g', + position_0based=10, + context='test', + errors=errors, + ) + assert result is not None + _, ref, mut = result + assert ref == 'A' + assert mut == 'G' + + def test_preserves_position(self): + """Should preserve 0-based position.""" + errors: list[str] = [] + for pos in [0, 10, 100, 1000]: + result = _normalize_rule_alleles_for_storage( + reference_aa='A', + mutation_raw='G', + position_0based=pos, + context='test', + errors=errors, + ) + assert result is not None + assert result[0] == pos diff --git a/tests/test_rules_normalize.py b/tests/test_rules_normalize.py new file mode 100644 index 00000000..babd6dd2 --- /dev/null +++ b/tests/test_rules_normalize.py @@ -0,0 +1,465 @@ +""" +Tests for TSV row normalization helpers in rule imports. + +Covers: respro/db/_rules_normalize.py +- _get_value() +- _parse_ic50_value() +- _parse_single_ic50() +- _normalize_ic50_from_row() +- _normalize_fold_ic50_from_row() +- _normalize_score_from_row() +- _normalize_phenotype_token() +- _append_contradictory_comment() +""" + +from __future__ import annotations + +from respro.db._rules_normalize import ( + _CONTRADICTORY_COMMENT, + _append_contradictory_comment, + _get_value, + _normalize_fold_ic50_from_row, + _normalize_ic50_from_row, + _normalize_phenotype_token, + _normalize_score_from_row, + _parse_ic50_value, + _parse_single_ic50, +) + + +class TestGetValue: + """Tests for _get_value() helper.""" + + def test_returns_first_non_empty_value(self): + """Should return the first non-empty value from the provided keys.""" + row = {'a': '', 'b': 'value_b', 'c': 'value_c'} + result = _get_value(row, 'a', 'b', 'c') + assert result == 'value_b' + + def test_returns_stripped_value(self): + """Should strip whitespace from the returned value.""" + row = {'key': ' value with spaces '} + result = _get_value(row, 'key') + assert result == 'value with spaces' + + def test_returns_empty_string_when_all_keys_missing(self): + """Should return empty string when no keys have values.""" + row = {'a': '', 'b': '', 'c': ''} + result = _get_value(row, 'a', 'b', 'c') + assert result == '' + + def test_returns_empty_string_when_row_empty(self): + """Should return empty string for empty row.""" + row: dict[str, str] = {} + result = _get_value(row, 'a', 'b', 'c') + assert result == '' + + def test_handles_none_values(self): + """Should skip None values and continue to next key.""" + row = {'a': None, 'b': 'found', 'c': 'ignored'} + result = _get_value(row, 'a', 'b', 'c') + assert result == 'found' + + def test_single_key_exists(self): + """Should work with single key that exists.""" + row = {'key': 'value'} + result = _get_value(row, 'key') + assert result == 'value' + + def test_single_key_missing(self): + """Should return empty string for single missing key.""" + row = {'other': 'value'} + result = _get_value(row, 'key') + assert result == '' + + +class TestParseIc50Value: + """Tests for _parse_ic50_value() numeric parsing.""" + + def test_parses_integer(self): + """Should parse integer values.""" + assert _parse_ic50_value('10') == 10.0 + assert _parse_ic50_value('0') == 0.0 + + def test_parses_float(self): + """Should parse decimal values.""" + assert _parse_ic50_value('10.5') == 10.5 + assert _parse_ic50_value('0.001') == 0.001 + + def test_parses_scientific_notation(self): + """Should parse scientific notation.""" + assert _parse_ic50_value('1.5e3') == 1500.0 + assert _parse_ic50_value('2.0E-3') == 0.002 + assert _parse_ic50_value('1e10') == 1e10 + + def test_parses_negative_values(self): + """Should parse negative values (though biologically unusual).""" + assert _parse_ic50_value('-5.0') == -5.0 + assert _parse_ic50_value('-1e2') == -100.0 + + def test_parses_with_whitespace(self): + """Should strip whitespace before parsing.""" + assert _parse_ic50_value(' 10.5 ') == 10.5 + assert _parse_ic50_value('\t5.0\n') == 5.0 + + def test_returns_none_for_empty(self): + """Should return None for empty string.""" + assert _parse_ic50_value('') is None + + def test_returns_none_for_whitespace_only(self): + """Should return None for whitespace-only string.""" + assert _parse_ic50_value(' ') is None + + def test_returns_none_for_none_literal(self): + """Should return None for 'none' string.""" + assert _parse_ic50_value('none') is None + assert _parse_ic50_value('None') is None + assert _parse_ic50_value('NONE') is None + + def test_returns_none_for_invalid_string(self): + """Should return None for non-numeric strings.""" + assert _parse_ic50_value('resistant') is None + assert _parse_ic50_value('N/A') is None + assert _parse_ic50_value('unknown') is None + + def test_extracts_number_from_text(self): + """Should extract numeric portion from mixed text.""" + assert _parse_ic50_value('>10.5') == 10.5 + assert _parse_ic50_value('~5.0') == 5.0 + assert _parse_ic50_value('approximately 10') == 10.0 + + def test_parses_leading_decimal_point(self): + """Should parse numbers starting with decimal point.""" + assert _parse_ic50_value('.5') == 0.5 + assert _parse_ic50_value('.123') == 0.123 + + +class TestParseSingleIc50: + """Tests for _parse_single_ic50() with error handling.""" + + def test_returns_canonical_integer(self): + """Should return canonical format for integers.""" + errors: list[str] = [] + result = _parse_single_ic50('10', errors=errors, context='test') + assert result == '10' + assert errors == [] + + def test_returns_canonical_float(self): + """Should return canonical format for floats (removes trailing zeros).""" + errors: list[str] = [] + result = _parse_single_ic50('10.500', errors=errors, context='test') + assert result == '10.5' + assert errors == [] + + def test_returns_canonical_integer_format(self): + """Should return integer format without decimal point for whole numbers.""" + errors: list[str] = [] + result = _parse_single_ic50('20.0', errors=errors, context='test') + assert result == '20' # :g format removes trailing .0 + assert errors == [] + + def test_returns_empty_for_none(self): + """Should return empty string for 'none' value.""" + errors: list[str] = [] + result = _parse_single_ic50('none', errors=errors, context='test') + assert result == '' + assert errors == [] + + def test_returns_empty_for_empty(self): + """Should return empty string for empty input.""" + errors: list[str] = [] + result = _parse_single_ic50('', errors=errors, context='test') + assert result == '' + assert errors == [] + + def test_appends_error_for_invalid_value(self): + """Should append error message for invalid values.""" + errors: list[str] = [] + result = _parse_single_ic50('invalid', errors=errors, context='test context') + assert result == '' + assert len(errors) == 1 + assert "test context: invalid ic50 value 'invalid'" in errors[0] + + def test_strips_whitespace(self): + """Should strip whitespace before parsing.""" + errors: list[str] = [] + result = _parse_single_ic50(' 10.5 ', errors=errors, context='test') + assert result == '10.5' + assert errors == [] + + +class TestNormalizeIc50FromRow: + """Tests for _normalize_ic50_from_row() column resolution.""" + + def test_reads_ic50_column(self): + """Should read from 'ic50' column.""" + row = {'ic50': '10.5', 'ic_50': '20.0'} + errors: list[str] = [] + result = _normalize_ic50_from_row(row, errors=errors, context='test') + assert result == '10.5' + + def test_reads_ic_50_column_as_fallback(self): + """Should read from 'ic_50' column when 'ic50' missing.""" + row = {'ic_50': '20.0'} + errors: list[str] = [] + result = _normalize_ic50_from_row(row, errors=errors, context='test') + assert result == '20' # :g format removes trailing .0 + + def test_prefers_ic50_over_ic_50(self): + """Should prefer 'ic50' over 'ic_50' when both present.""" + row = {'ic50': '10.5', 'ic_50': '20.0'} + errors: list[str] = [] + result = _normalize_ic50_from_row(row, errors=errors, context='test') + assert result == '10.5' + + def test_returns_empty_for_missing_columns(self): + """Should return empty string when columns missing.""" + row = {'other': 'value'} + errors: list[str] = [] + result = _normalize_ic50_from_row(row, errors=errors, context='test') + assert result == '' + assert errors == [] + + def test_appends_error_for_invalid_value(self): + """Should append error for invalid IC50 value.""" + row = {'ic50': 'not_a_number'} + errors: list[str] = [] + result = _normalize_ic50_from_row(row, errors=errors, context='row 5') + assert result == '' + assert len(errors) == 1 + assert 'row 5: invalid ic50 value' in errors[0] + + +class TestNormalizeFoldIc50FromRow: + """Tests for _normalize_fold_ic50_from_row() column resolution.""" + + def test_reads_fold_ic50_column(self): + """Should read from 'fold_ic50' column.""" + row = {'fold_ic50': '5.5', 'fold_ic_50': '10.0'} + errors: list[str] = [] + result = _normalize_fold_ic50_from_row(row, errors=errors, context='test') + assert result == '5.5' + + def test_reads_fold_ic_50_column_as_fallback(self): + """Should read from 'fold_ic_50' column when 'fold_ic50' missing.""" + row = {'fold_ic_50': '10.0'} + errors: list[str] = [] + result = _normalize_fold_ic50_from_row(row, errors=errors, context='test') + assert result == '10' # :g format removes trailing .0 + + def test_prefers_fold_ic50_over_fold_ic_50(self): + """Should prefer 'fold_ic50' over 'fold_ic_50'.""" + row = {'fold_ic50': '5.5', 'fold_ic_50': '10.0'} + errors: list[str] = [] + result = _normalize_fold_ic50_from_row(row, errors=errors, context='test') + assert result == '5.5' + + def test_returns_empty_for_missing_columns(self): + """Should return empty string when columns missing.""" + row = {'other': 'value'} + errors: list[str] = [] + result = _normalize_fold_ic50_from_row(row, errors=errors, context='test') + assert result == '' + + +class TestNormalizeScoreFromRow: + """Tests for _normalize_score_from_row().""" + + def test_parses_integer_score(self): + """Should parse integer scores.""" + row = {'score': '5'} + errors: list[str] = [] + result = _normalize_score_from_row(row, errors=errors, context='test') + assert result == '5' + + def test_parses_float_score(self): + """Should parse float scores.""" + row = {'score': '5.75'} + errors: list[str] = [] + result = _normalize_score_from_row(row, errors=errors, context='test') + assert result == '5.75' + + def test_returns_empty_for_none(self): + """Should return empty string for 'none' value.""" + row = {'score': 'none'} + errors: list[str] = [] + result = _normalize_score_from_row(row, errors=errors, context='test') + assert result == '' + + def test_returns_empty_for_missing_column(self): + """Should return empty string when column missing.""" + row = {'other': 'value'} + errors: list[str] = [] + result = _normalize_score_from_row(row, errors=errors, context='test') + assert result == '' + + def test_appends_error_for_invalid_score(self): + """Should append error for invalid score value.""" + row = {'score': 'not_a_number'} + errors: list[str] = [] + result = _normalize_score_from_row(row, errors=errors, context='row 10') + assert result == '' + assert len(errors) == 1 + assert 'row 10: invalid score value' in errors[0] + + def test_uses_score_column_only(self): + """Should only read from 'score' column.""" + row = {'score': '5.0', 'other_score': '10.0'} + errors: list[str] = [] + result = _normalize_score_from_row(row, errors=errors, context='test') + assert result == '5' # :g format removes trailing .0 + + +class TestNormalizePhenotypeToken: + """Tests for _normalize_phenotype_token().""" + + def test_maps_resistant_variants(self): + """Should map resistant variants to 'resistant'.""" + assert _normalize_phenotype_token('resistant') == 'resistant' + assert _normalize_phenotype_token('resistance') == 'resistant' + assert _normalize_phenotype_token('res') == 'resistant' + + def test_maps_sensitive_variants(self): + """Should map sensitive variants to 'sensitive'.""" + assert _normalize_phenotype_token('sensitive') == 'sensitive' + assert _normalize_phenotype_token('susceptible') == 'sensitive' + assert _normalize_phenotype_token('sens') == 'sensitive' + + def test_maps_intermediate_variants(self): + """Should map intermediate variants.""" + assert _normalize_phenotype_token('intermediate') == 'intermediate' + assert _normalize_phenotype_token('interm') == 'intermediate' + assert _normalize_phenotype_token('i') == 'intermediate' + + def test_maps_contradictory(self): + """Should map contradictory variants.""" + assert _normalize_phenotype_token('contradictory') == 'contradictory' + assert _normalize_phenotype_token('contra') == 'contradictory' + assert _normalize_phenotype_token('conflict') == 'contradictory' + assert _normalize_phenotype_token('conflicting') == 'contradictory' + + def test_maps_unknown_to_unknown(self): + """Should map unknown/empty to 'unknown'.""" + assert _normalize_phenotype_token('unknown') == 'unknown' + assert _normalize_phenotype_token('none') == 'unknown' + assert _normalize_phenotype_token('') == 'unknown' + assert _normalize_phenotype_token(' ') == 'unknown' + + def test_case_insensitive(self): + """Should be case-insensitive.""" + assert _normalize_phenotype_token('RESISTANT') == 'resistant' + assert _normalize_phenotype_token('Resistant') == 'resistant' + assert _normalize_phenotype_token('ReSiStAnT') == 'resistant' + + def test_strips_whitespace(self): + """Should strip whitespace.""" + assert _normalize_phenotype_token(' resistant ') == 'resistant' + assert _normalize_phenotype_token('\tsensitive\n') == 'sensitive' + + def test_returns_none_for_unmapped(self): + """Should return None for unmapped phenotypes.""" + assert _normalize_phenotype_token('invalid') is None + assert _normalize_phenotype_token('mixed') is None + + +class TestAppendContradictoryComment: + """Tests for _append_contradictory_comment().""" + + def test_returns_original_when_not_contradictory(self): + """Should return original comment when phenotype is not contradictory.""" + result = _append_contradictory_comment( + 'Some comment', + phenotype='resistant', + clinical_phenotype='resistant', + ) + assert result == 'Some comment' + + def test_returns_standard_comment_when_empty(self): + """Should return standard comment when input is empty.""" + result = _append_contradictory_comment( + '', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == _CONTRADICTORY_COMMENT + + def test_appends_when_comment_present(self): + """Should append standard comment to existing comment.""" + result = _append_contradictory_comment( + 'Multiple studies show different results', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == 'Multiple studies show different results. Publications have contradictory phenotype associations.' + + def test_adds_period_before_appending(self): + """Should add period if comment doesn't end with punctuation.""" + result = _append_contradictory_comment( + 'Results vary', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == 'Results vary. Publications have contradictory phenotype associations.' + + def test_preserves_existing_punctuation(self): + """Should preserve existing punctuation.""" + result = _append_contradictory_comment( + 'Results vary.', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == 'Results vary. Publications have contradictory phenotype associations.' + + def test_handles_exclamation_mark(self): + """Should handle exclamation mark as existing punctuation.""" + result = _append_contradictory_comment( + 'Results vary!', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == 'Results vary! Publications have contradictory phenotype associations.' + + def test_handles_question_mark(self): + """Should handle question mark as existing punctuation.""" + result = _append_contradictory_comment( + 'Results vary?', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == 'Results vary? Publications have contradictory phenotype associations.' + + def test_skips_if_already_contains_standard_comment(self): + """Should not duplicate if comment already contains standard text.""" + result = _append_contradictory_comment( + _CONTRADICTORY_COMMENT, + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == _CONTRADICTORY_COMMENT + + def test_case_insensitive_check_for_existing_comment(self): + """Should check for existing comment case-insensitively.""" + result = _append_contradictory_comment( + 'publications have contradictory phenotype associations.', + phenotype='contradictory', + clinical_phenotype='contradictory', + ) + assert result == 'publications have contradictory phenotype associations.' + + def test_works_with_clinical_phenotype_only(self): + """Should work when only clinical_phenotype is contradictory.""" + result = _append_contradictory_comment( + 'Comment', + phenotype='resistant', + clinical_phenotype='contradictory', + ) + assert 'Publications have contradictory' in result + + def test_works_with_phenotype_only(self): + """Should work when only phenotype is contradictory.""" + result = _append_contradictory_comment( + 'Comment', + phenotype='contradictory', + clinical_phenotype='resistant', + ) + assert 'Publications have contradictory' in result diff --git a/tests/test_utils_cli.py b/tests/test_utils_cli.py new file mode 100644 index 00000000..4334f67c --- /dev/null +++ b/tests/test_utils_cli.py @@ -0,0 +1,150 @@ +""" +Tests for CLI utilities - error handling and logging. + +Covers: respro/utils/cli_errors.py, respro/utils/logging.py +- cli_error() +- render_click_exception() +- setup_logging() +- _LogHighlighter +""" + +from __future__ import annotations + +import logging + +import click +import pytest +import typer + +from respro.utils.cli_errors import cli_error, render_click_exception +from respro.utils.logging import _LogHighlighter, setup_logging + + +class TestCliError: + """Tests for cli_error().""" + + def test_raises_typer_exit(self): + """Should raise typer.Exit.""" + with pytest.raises(typer.Exit) as exc_info: + cli_error('test error') + assert exc_info.value.exit_code == 1 + + def test_default_exit_code_is_1(self): + """Should exit with code 1 by default.""" + with pytest.raises(typer.Exit) as exc_info: + cli_error('error') + assert exc_info.value.exit_code == 1 + + def test_custom_exit_code(self): + """Should support custom exit codes.""" + with pytest.raises(typer.Exit) as exc_info: + cli_error('error', exit_code=2) + assert exc_info.value.exit_code == 2 + + def test_message_includes_error_prefix(self, capsys): + """Should prepend 'Error: ' to message.""" + try: + cli_error('something went wrong') + except typer.Exit: + pass + captured = capsys.readouterr() + assert 'Error: something went wrong' in captured.err + + +class TestRenderClickException: + """Tests for render_click_exception().""" + + def test_raises_typer_exit(self): + """Should raise typer.Exit.""" + exc = click.ClickException('test') + with pytest.raises(typer.Exit): + render_click_exception(exc) + + def test_preserves_exit_code(self): + """Should preserve exception's exit code.""" + exc = click.ClickException('test') + exc.exit_code = 3 + with pytest.raises(typer.Exit) as exc_info: + render_click_exception(exc) + assert exc_info.value.exit_code == 3 + + def test_works_with_usage_error(self, capsys): + """Should work with UsageError.""" + exc = click.UsageError('usage problem') + try: + render_click_exception(exc) + except typer.Exit: + pass + captured = capsys.readouterr() + assert 'usage problem' in captured.err.lower() + + +class TestSetupLogging: + """Tests for setup_logging().""" + + def test_returns_logger(self): + """Should return Logger instance.""" + logger = setup_logging() + assert isinstance(logger, logging.Logger) + assert logger.name == 'respro' + + def test_default_level_is_warning(self): + """Should default to WARNING level.""" + logger = setup_logging() + assert logger.level == logging.WARNING + + def test_verbosity_0_is_warning(self): + """Should set WARNING for verbosity 0.""" + logger = setup_logging(verbosity=0) + assert logger.level == logging.WARNING + + def test_verbosity_1_is_info(self): + """Should set INFO for verbosity 1.""" + logger = setup_logging(verbosity=1) + assert logger.level == logging.INFO + + def test_verbosity_2_is_debug(self): + """Should set DEBUG for verbosity 2+.""" + logger = setup_logging(verbosity=2) + assert logger.level == logging.DEBUG + logger2 = setup_logging(verbosity=10) + assert logger2.level == logging.DEBUG + + def test_handler_attached(self): + """Should attach a handler.""" + logger = setup_logging() + assert len(logger.handlers) > 0 + + def test_rich_handler_used(self): + """Should use RichHandler.""" + from rich.logging import RichHandler + logger = setup_logging() + assert any(isinstance(h, RichHandler) for h in logger.handlers) + + def test_idempotent(self): + """Should be idempotent (not add multiple handlers).""" + logger1 = setup_logging() + handler_count = len(logger1.handlers) + logger2 = setup_logging() + assert len(logger2.handlers) == handler_count + + +class TestLogHighlighter: + """Tests for _LogHighlighter.""" + + def test_instantiates(self): + """Should instantiate without error.""" + highlighter = _LogHighlighter() + assert highlighter is not None + + def test_has_regex_patterns(self): + """Should have regex patterns defined.""" + highlighter = _LogHighlighter() + assert len(highlighter.highlights) > 0 + + def test_patterns_match_strings(self): + """Should have patterns that match expected strings.""" + highlighter = _LogHighlighter() + # Just verify the patterns compile and exist + assert any('quoted' in str(p) or 'str_val' in str(p) for p in highlighter.highlights) + assert any('module' in str(p) or 'respro' in str(p) for p in highlighter.highlights) diff --git a/tests/test_web_api.py b/tests/test_web_api.py index 0c684a4f..d69e5ddb 100644 --- a/tests/test_web_api.py +++ b/tests/test_web_api.py @@ -623,6 +623,31 @@ def test_databases_endpoint_includes_project_metadata( assert database['metadata']['contact'] == 'team@example.org' assert database['metadata']['license'] == 'MIT' + def test_databases_endpoint_has_example_false_by_default(self, client: TestClient) -> None: + response = client.get('/api/databases') + assert response.status_code == 200 + database = response.json()['data']['items'][0] + assert database['has_example'] is False + + def test_databases_endpoint_has_example_true_when_stored( + self, + client: TestClient, + startup_config: StartupConfig, + ) -> None: + project_db = sorted(startup_config.project_databases_dir.glob('*.db'))[0] + conn = sqlite3.connect(project_db) + conn.execute( + "UPDATE project SET example_fasta = ? WHERE id = 1", + ('>example\nATGAAAGCT\n',), + ) + conn.commit() + conn.close() + + response = client.get('/api/databases') + assert response.status_code == 200 + database = response.json()['data']['items'][0] + assert database['has_example'] is True + def test_mutations_endpoint_alias(self, client: TestClient) -> None: project_db = sorted(client.app.state.startup_config.project_databases_dir.glob('*.db'))[0] conn = sqlite3.connect(project_db) @@ -713,6 +738,56 @@ def test_profile_fasta( report_payload = _download_artifact_json(client, json_id) assert report_payload['run']['vcf_path'] == 'original-upload.fasta' + def test_profile_fasta_use_example_profiles_stored_example( + self, + client: TestClient, + startup_config: StartupConfig, + ) -> None: + """use_example=true enqueues a fasta --example job without an uploaded FASTA.""" + project_db = sorted(startup_config.project_databases_dir.glob('*.db'))[0] + conn = sqlite3.connect(project_db) + # Store the tiny reference itself as the example so it profiles cleanly (0 hits). + conn.execute( + "UPDATE project SET example_fasta = ? WHERE id = 1", + (f'>stored_example\n{TINY_REF_SEQ}\n',), + ) + conn.commit() + conn.close() + + submit = client.post( + '/api/profile/fasta', + json={ + 'use_example': True, + 'sample': 'web-example', + }, + ) + assert submit.status_code == 200, submit.text + job_id = submit.json()['job_id'] + assert job_id + + payload = _poll_job(client, job_id) + assert payload['status'] == 'succeeded', payload + result = payload['result'] + assert result['database_id'] == project_db.name + assert result['sample_name'] == 'web-example' + + def test_profile_fasta_use_example_rejects_missing_example( + self, + client: TestClient, + startup_config: StartupConfig, + ) -> None: + """use_example=true on a DB without an example fails the job with a clear error.""" + submit = client.post( + '/api/profile/fasta', + json={'use_example': True, 'sample': 'web-example'}, + ) + assert submit.status_code == 200 + job_id = submit.json()['job_id'] + + payload = _poll_job(client, job_id) + assert payload['status'] == 'failed' + assert 'example' in (payload.get('error') or '').lower() + def test_profile_fasta_path_outside_uploads_rejected( self, client: TestClient, diff --git a/web/backend/jobs.py b/web/backend/jobs.py index 26b445aa..415e50ce 100644 --- a/web/backend/jobs.py +++ b/web/backend/jobs.py @@ -25,16 +25,17 @@ def run_profile_fasta( *, project_db: str, output_dir: str, - fasta_path: str, + fasta_path: str | None, sample: str, threads: int, input_display_name: str | None = None, artifact_base_name: str | None = None, + use_example: bool = False, ) -> dict: """RQ job wrapper for FASTA profiling.""" output_html_path = _build_web_output_html_path( output_dir=Path(output_dir), - input_name=artifact_base_name or input_display_name or Path(fasta_path).name, + input_name=artifact_base_name or input_display_name or (Path(fasta_path).name if fasta_path else 'example'), ) return _run_job_with_logging( mode='fasta', @@ -43,10 +44,11 @@ def run_profile_fasta( job_func=lambda: _run_profile_fasta_subprocess( project_db=Path(project_db), output_html_path=output_html_path, - fasta_path=Path(fasta_path), + fasta_path=Path(fasta_path) if fasta_path else None, sample=sample, threads=threads, input_display_name=input_display_name, + use_example=use_example, ), ) @@ -116,10 +118,11 @@ def _run_profile_fasta_subprocess( *, project_db: Path, output_html_path: Path, - fasta_path: Path, + fasta_path: Path | None, sample: str, threads: int, input_display_name: str | None, + use_example: bool = False, ) -> dict: """Execute FASTA profiling through the respro CLI and return the web API payload.""" command = [ @@ -127,8 +130,12 @@ def _run_profile_fasta_subprocess( 'fasta', '--project', str(project_db), - '--fasta', - str(fasta_path), + ] + if use_example: + command.append('--example') + else: + command.extend(['--fasta', str(fasta_path)]) + command.extend([ '--sample', sample, '--output', @@ -140,7 +147,7 @@ def _run_profile_fasta_subprocess( 'json', '--export', 'pdf', - ] + ]) if input_display_name: command.extend(['--input-display-name', input_display_name]) @@ -154,7 +161,7 @@ def _run_profile_fasta_subprocess( 'sample_name': run_payload.get('sample_name', sample), 'created_at': run_payload.get('created_at', ''), 'reference_name': run_payload.get('reference_name', ''), - 'query_name': _read_first_fasta_name(fasta_path), + 'query_name': _read_first_fasta_name(fasta_path) if fasta_path else 'example', 'report_html_path': str(artifacts['html']), 'report_json_path': str(artifacts['json']), 'report_pdf_path': str(artifacts['pdf']), diff --git a/web/backend/models.py b/web/backend/models.py index 71b0df6c..569294e2 100644 --- a/web/backend/models.py +++ b/web/backend/models.py @@ -15,7 +15,8 @@ class ProfileFastaPayload(BaseModel): ownership-checks them server-side. """ - fasta_id: str = Field(max_length=WEB_BACKEND_CONFIG.defaults.opaque_id_max_length) + fasta_id: str | None = Field(default=None, max_length=WEB_BACKEND_CONFIG.defaults.opaque_id_max_length) + use_example: bool = False input_display_name: str | None = Field(default=None, max_length=WEB_BACKEND_CONFIG.defaults.display_name_max_length) database_id: str | None = Field(default=None, max_length=WEB_BACKEND_CONFIG.defaults.opaque_id_max_length) sample: str | None = Field(default=None, max_length=WEB_BACKEND_CONFIG.defaults.sample_name_max_length) diff --git a/web/backend/routes/profile.py b/web/backend/routes/profile.py index 5631561a..a302a79b 100644 --- a/web/backend/routes/profile.py +++ b/web/backend/routes/profile.py @@ -55,12 +55,22 @@ def profile_fasta_route( queue: Queue = Depends(get_queue), ) -> JobSubmitResponse: session = get_session(request) - fasta_path = _resolve_upload_path( - payload.fasta_id, - session.session_hash, - config.allowed_roots, - label='FASTA', - ) + if payload.use_example: + fasta_path = None + fasta_id = None + else: + if not payload.fasta_id: + raise HTTPException( + status_code=422, + detail='fasta_id is required when use_example is false.', + ) + fasta_path = _resolve_upload_path( + payload.fasta_id, + session.session_hash, + config.allowed_roots, + label='FASTA', + ) + fasta_id = payload.fasta_id consume_sample_quota( request, sample_count=1, @@ -73,13 +83,18 @@ def profile_fasta_route( run_profile_fasta, project_db=str(project_db), output_dir=str(config.results_dir), - fasta_path=str(fasta_path), + fasta_path=str(fasta_path) if fasta_path else None, + use_example=payload.use_example, sample=payload.sample or defaults.profile_sample_name, threads=payload.threads if payload.threads is not None else defaults.profile_threads, input_display_name=payload.input_display_name, **enqueue_options, ) - record_job(session_hash=session.session_hash, upload_ids=[payload.fasta_id], job_id=job.id) + record_job( + session_hash=session.session_hash, + upload_ids=[fasta_id] if fasta_id else [], + job_id=job.id, + ) logger.info('Queue job enqueued: job_id=%s mode=fasta database_id=%s', job.id, project_db.name) return JobSubmitResponse(job_id=job.id) diff --git a/web/backend/services/browse.py b/web/backend/services/browse.py index 9478a5d9..d6732f2a 100644 --- a/web/backend/services/browse.py +++ b/web/backend/services/browse.py @@ -67,6 +67,7 @@ def list_databases(project_databases_dir: Path) -> dict: 'schema_version': project_row['schema_version'], 'supported_organisms': organisms, 'mutation_count': mutation_count, + 'has_example': bool(project_row.get('example_fasta')), 'metadata': metadata, 'algorithms': algorithms, } diff --git a/web/backend/services/maintained_bootstrap.py b/web/backend/services/maintained_bootstrap.py index b1b9784e..d44bc702 100644 --- a/web/backend/services/maintained_bootstrap.py +++ b/web/backend/services/maintained_bootstrap.py @@ -130,6 +130,7 @@ def _rebuild_database_files(name: str, dest_db_path: Path) -> None: rules_tsv=files['rules'], # type: ignore[arg-type] formula_rules_tsv=files['formula_rules'], # type: ignore[arg-type] metadata_json=files['metadata'], # type: ignore[arg-type] + example_fasta=files['example'], # type: ignore[arg-type] overwrite=False, additional_info=True, ) diff --git a/web/frontend/src/assets/example.svg b/web/frontend/src/assets/example.svg new file mode 100644 index 00000000..802f400c --- /dev/null +++ b/web/frontend/src/assets/example.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/web/frontend/src/components/DashboardView.jsx b/web/frontend/src/components/DashboardView.jsx index cb5c8fcb..1566eab0 100644 --- a/web/frontend/src/components/DashboardView.jsx +++ b/web/frontend/src/components/DashboardView.jsx @@ -75,6 +75,7 @@ export function DashboardView({ isCancelingJob, cancelActiveJob, runSelectedProfile, + runExampleProfile, buildReportUrl, buildArtifactUrl, uploadFastaFile, @@ -207,6 +208,7 @@ export function DashboardView({ isCancelingJob={isCancelingJob} cancelActiveJob={cancelActiveJob} runSelectedProfile={runSelectedProfile} + runExampleProfile={runExampleProfile} statusError={statusError} selectedProfileReportPath={selectedProfileReportPath} setSelectedProfileReportPath={setSelectedProfileReportPath} diff --git a/web/frontend/src/components/tabs/AnalyzeTab.jsx b/web/frontend/src/components/tabs/AnalyzeTab.jsx index f5cac841..4a9d9098 100644 --- a/web/frontend/src/components/tabs/AnalyzeTab.jsx +++ b/web/frontend/src/components/tabs/AnalyzeTab.jsx @@ -2,6 +2,7 @@ import { useEffect, useRef, useState } from 'react'; import analyzeIconSrc from '../../assets/icon-analyze.svg'; import batchIconSrc from '../../assets/batch.svg'; +import exampleIconSrc from '../../assets/example.svg'; import infoIconSrc from '../../assets/info.svg'; import { Spinner } from '../Spinner'; @@ -34,6 +35,7 @@ export function AnalyzeTab({ isCancelingJob, cancelActiveJob, runSelectedProfile, + runExampleProfile, statusError, selectedProfileReportPath, setSelectedProfileReportPath, @@ -134,6 +136,22 @@ export function AnalyzeTab({

+ {selectedDatabase?.has_example ? ( + + ) : null} {ANALYZE_SUBMODES.map((subMode) => (