From 192cb211830b982fbbee0e62aa596b8253407818 Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Thu, 6 Aug 2026 18:34:13 +0200 Subject: [PATCH 1/7] Rewrite CANDy as an installable Python package (CLI + API) Reimplements all 13 stages of the Colab notebook pipeline as a proper pip-installable package (src/candy) instead of Colab-only cells, with pluggable clustering (CD-HIT/MMseqs2), alignment (MAFFT), phylogenetics (FastTree), and domain-curation (manual/Gemini) backends so future tool swaps don't require another rewrite. Adds a Typer CLI (`candy run`), a SQLAlchemy-backed result database, and 79 tests covering the pure-logic modules plus mocked end-to-end pipeline runs for both input modes. While porting, several real bugs in the original notebook were found and fixed (each documented in-code and covered by a regression test): duplicated taxonomy codes corrupting characterized-enzyme organism names, a BLAST-fallback crash from state leaking across loop iterations, swapped taxonomy/organism_name columns in the SQLite schema, and a value-based list.index() bug mis-pairing domain names to positions. Co-Authored-By: Claude Sonnet 5 --- .gitignore | 15 ++ README.md | 175 ++++----------- environment.yml | 13 ++ pyproject.toml | 51 +++++ src/candy/__init__.py | 3 + src/candy/alignment/__init__.py | 22 ++ src/candy/alignment/mafft.py | 25 +++ src/candy/blast.py | 110 ++++++++++ src/candy/cazy.py | 217 +++++++++++++++++++ src/candy/cli.py | 129 +++++++++++ src/candy/clustering/__init__.py | 36 ++++ src/candy/clustering/cdhit.py | 42 ++++ src/candy/clustering/mmseqs2.py | 41 ++++ src/candy/config.py | 143 +++++++++++++ src/candy/curation/__init__.py | 37 ++++ src/candy/curation/gemini.py | 58 +++++ src/candy/curation/manual.py | 41 ++++ src/candy/database.py | 156 ++++++++++++++ src/candy/domains.py | 246 +++++++++++++++++++++ src/candy/external_tools.py | 65 ++++++ src/candy/fasta.py | 169 +++++++++++++++ src/candy/interpro.py | 230 ++++++++++++++++++++ src/candy/itol.py | 168 +++++++++++++++ src/candy/merge.py | 81 +++++++ src/candy/network.py | 80 +++++++ src/candy/phylogenetics/__init__.py | 22 ++ src/candy/phylogenetics/fasttree.py | 33 +++ src/candy/pipeline.py | 250 ++++++++++++++++++++++ src/candy/uniparc.py | 58 +++++ tests/test_alignment_and_phylogenetics.py | 35 +++ tests/test_blast.py | 84 ++++++++ tests/test_cazy.py | 42 ++++ tests/test_cli.py | 55 +++++ tests/test_clustering.py | 55 +++++ tests/test_curation.py | 49 +++++ tests/test_database.py | 111 ++++++++++ tests/test_domains.py | 120 +++++++++++ tests/test_external_tools.py | 38 ++++ tests/test_fasta.py | 99 +++++++++ tests/test_interpro.py | 68 ++++++ tests/test_itol.py | 78 +++++++ tests/test_merge.py | 62 ++++++ tests/test_network.py | 28 +++ tests/test_pipeline.py | 155 ++++++++++++++ tests/test_uniparc.py | 19 ++ 45 files changed, 3686 insertions(+), 128 deletions(-) create mode 100644 .gitignore create mode 100644 environment.yml create mode 100644 pyproject.toml create mode 100644 src/candy/__init__.py create mode 100644 src/candy/alignment/__init__.py create mode 100644 src/candy/alignment/mafft.py create mode 100644 src/candy/blast.py create mode 100644 src/candy/cazy.py create mode 100644 src/candy/cli.py create mode 100644 src/candy/clustering/__init__.py create mode 100644 src/candy/clustering/cdhit.py create mode 100644 src/candy/clustering/mmseqs2.py create mode 100644 src/candy/config.py create mode 100644 src/candy/curation/__init__.py create mode 100644 src/candy/curation/gemini.py create mode 100644 src/candy/curation/manual.py create mode 100644 src/candy/database.py create mode 100644 src/candy/domains.py create mode 100644 src/candy/external_tools.py create mode 100644 src/candy/fasta.py create mode 100644 src/candy/interpro.py create mode 100644 src/candy/itol.py create mode 100644 src/candy/merge.py create mode 100644 src/candy/network.py create mode 100644 src/candy/phylogenetics/__init__.py create mode 100644 src/candy/phylogenetics/fasttree.py create mode 100644 src/candy/pipeline.py create mode 100644 src/candy/uniparc.py create mode 100644 tests/test_alignment_and_phylogenetics.py create mode 100644 tests/test_blast.py create mode 100644 tests/test_cazy.py create mode 100644 tests/test_cli.py create mode 100644 tests/test_clustering.py create mode 100644 tests/test_curation.py create mode 100644 tests/test_database.py create mode 100644 tests/test_domains.py create mode 100644 tests/test_external_tools.py create mode 100644 tests/test_fasta.py create mode 100644 tests/test_interpro.py create mode 100644 tests/test_itol.py create mode 100644 tests/test_merge.py create mode 100644 tests/test_network.py create mode 100644 tests/test_pipeline.py create mode 100644 tests/test_uniparc.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4503525 --- /dev/null +++ b/.gitignore @@ -0,0 +1,15 @@ +__pycache__/ +*.py[cod] +*.egg-info/ +.eggs/ +build/ +dist/ +.venv/ +venv/ +.pytest_cache/ +.mypy_cache/ +.ruff_cache/ +.coverage +htmlcov/ +*.db +.env diff --git a/README.md b/README.md index 35c6794..5a21e59 100644 --- a/README.md +++ b/README.md @@ -1,144 +1,64 @@ # Carbohydrate Active eNzyme Domain analYsis tool (CANDy) - automated analysis of domain architectures in carbohydrate-active enzymes -CANDy boosts a fast, FAIR and seamless protein domain analysis of any [CAZy](http://www.cazy.org/) family. Check the latest online version on [Google Colab](https://colab.research.google.com/drive/1ipRAwMFMDRGUinPDk2bwu1cg8fE2WY8Q?usp=sharing), yet for bigger families we recommend you downloading the Jupyter Notebook. -#### April 2025 update : We included MMseqs2 as an alternative clustering option for CD-HIT. Additionally, the domain name curation step is now fully automated by employing Gemini 2.0 flash. Users still have the option to perform this step manually. (Only available on Colab) +CANDy is a fast, FAIR and seamless protein domain analysis tool for any [CAZy](http://www.cazy.org/) family. -## Requirements +**This branch is a work-in-progress rewrite of CANDy as an installable Python package** (CLI + Python API), replacing the original Google Colab / Jupyter Notebook implementation. The previous notebook (`CANDy v2.0.ipynb`) is still in this repository for reference, and the original is on [Google Colab](https://colab.research.google.com/drive/1ipRAwMFMDRGUinPDk2bwu1cg8fE2WY8Q?usp=sharing). -Make sure to have following tools installed in the same directory as the Jupyter Notebook: +## Installation -### 1. CD-HIT +CANDy depends on a handful of external bioinformatics tools (CD-HIT or MMseqs2 for clustering, MAFFT for alignment, FastTree for phylogenetics) that aren't distributed via PyPI. The supported install path is a conda environment that provides them, with CANDy itself installed via pip into that environment: -Download the source code for CD-HIT from the GitHub repository at https://github.com/weizhongli/cdhit/releases and follow the [installation instructions](https://github.com/weizhongli/cdhit/wiki/2.-Installation). - -### 2. MAFFT +```bash +conda env create -f environment.yml +conda activate candy +``` -Precompiled binary can be downloaded here: https://mafft.cbrc.jp/alignment/software/. Change the installation directory to the path where the Notebook is stored or manually move the executable from the default directory. You can find the location of MAFFT by typing the following command in your terminal: +This installs the external tools via bioconda and CANDy itself (editable) via pip. If you already have those tools on your PATH through some other means, you can also just `pip install -e .` directly. - where mafft +To also enable automated Gemini-based domain-name curation: -### 3. FastTree +```bash +pip install -e ".[gemini]" +``` -#### MacOS +## Usage - Follow installation instructions for your system on http://www.microbesonline.org/fasttree/#Install - - Or - - Open a terminal window and install the Xcode Command Line Tools by typing the following command: - - xcode-select --install - - Install Homebrew by typing the following command: - - /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)" +### Command line - Install FastTree by typing the following command: - - brew install fasttree +```bash +# Query a CAZy family directly +candy run --jobname my_gh5_run --family GH5 --email you@example.com --tree -#### Linux +# Analyse your own FASTA file instead +candy run --jobname my_custom_run --fasta my_sequences.fasta --tree +``` - Follow installation instructions for your system on http://www.microbesonline.org/fasttree/#Install - - Or +Run `candy run --help` for the full list of options (taxonomy subset, clustering software/cutoff, BLAST identity threshold, max domain length, domain-overlap threshold, curation backend, ...). - Download the FastTree source code from the FastTree website at http://www.microbesonline.org/fasttree/. +### Python API - Open a terminal window and navigate to the directory where you downloaded the source code. +```python +from candy.config import CAZyFamilyInput, PipelineConfig, Taxonomy +from candy.pipeline import run_pipeline - Type the following commands to compile the FastTree package: - - tar xvzf FastTree-2.1.13.c - cd FastTree-2.1.13 - make - - After the compilation process completes, you will find the FastTree executable file in the "FastTree-2.1.13" directory. - -#### Windows - - Download the FastTree Windows binary from the FastTree website at http://www.microbesonline.org/fasttree/#Download. - - Extract the FastTree executable from the downloaded ZIP file. - - - -## Running CANDy -Your directory shoud look like this: - -image - - -This Notebook uses several Python packages. To avoid compatibility issues we recommend running this Notebook in a __virtual environment__. -- Therefore, install [Anaconda](https://www.anaconda.com/) and follow the installation instructions. -- Go to 'Environments' in Anaconda and click 'create'. Give your environment a name, for example '__myenv__'. The virual environment will be launched automatically. -- Go to the package search bar and search for '__ipywidgets__'. Download the package to be able to use the interactive widgets in this Notebook. Repeat for the '__h5py__' package. -- Next, go back to the 'Home' page in Anaconda and install __Jupyter Notebook__. Once completed, press launch and go to the directory where you saved this Notebook. -- Verify that you see the name of the virtual enivronment on the right top of the Notebook, for example: __Python (myenv)__. If that's not the case, go to Kernel and choose the environment. - -Also, for large families, avoid your computer entering sleep or stand-by mode since this will interupt the run. Change the settings in your computer or [caffeinate](https://pypi.org/project/caffeinate/) your system. - -#### MacOS - -Install caffeinate package by running in your Terminal: - - brew install caffeinate - -Start the package by running: - - caffeinate -d - -Stop by running: - - ctrl + C - -#### Linux - -Install caffeinate package by running in your Terminal: - - sudo apt-get install caffeinate - -Start the package by running: - - caffeinate -d - -Stop by running: - - ctrl + C -#### Windows - -Note: The caffeinate package is not available for Windows. However, you can use a similar feature called "powercfg" to prevent the system from going to sleep. - -Open the Command Prompt application. - -Type following line to see the current power requests: - - powercfg /requests - -Type following line, followed by the type of request you want to override (e.g., "system" or "display"): - - powercfg /requestsoverride - -To stop the power request override, type in "powercfg /requestsoverride" followed by the type of request and the "/remove" argument +config = PipelineConfig( + input=CAZyFamilyInput(enzyme_class="GH", family_number=5, email="you@example.com", taxonomy=Taxonomy.ALL), + jobname="my_gh5_run", + output_dir="results", + build_tree=True, +) +result = run_pipeline(config) +print(result.database_path, result.tree_path) +``` ## Output -When running the Google Colab version of CANDy, results containing the FATSA files, SQLite database, MSA, phylogenetic tree (in Newick format) and iTOL annotation files are automatically downloaded in a Zip file. When running CANDy locally, these outputs are stores in the same directory as the Jupyter Notebook. -### Database - -To open the results in the database, download SQLite from: https://sqlitebrowser.org/ - -### (Annotated) Phylogenetic tree - -To view the phylogenetic tree, several free services are available. The Notebook makes use of the ete3 package to visualize the annotated tree in there. For a more interactive experience we recommend [iTOL](https://itol.embl.de/). The script outputs iTOL annotation files for the visualization of the protein domains and the activity of the included characterized sequences. - -![image](https://github.com/PyEED/CANDy/assets/72694200/889b74eb-f740-4aff-80a0-855359099dc0) - -### Protein domain co-occurence network - -CANDy offers users a co-occurrence network that visually represents both the frequency of different domain types and the degree to which they are interconnected. A simple visualisation is offered in the Notebook, but for a more interactive experience we recommend using [Cytoscape](https://cytoscape.org/) (yFiles Organic Layout). +Results are written to `{output_dir}/{jobname}/`: -![image](https://github.com/PyEED/CANDy/assets/72694200/cd20a756-7ded-4964-a9b9-04f7d83698a1) +- FASTA files for each processing stage +- A SQLite database (`{jobname}_db.db`) containing the domain annotations -- open it with [DB Browser for SQLite](https://sqlitebrowser.org/) +- A protein domain co-occurrence network (`{jobname}_domain_cooccurence_network.graphml`) -- open it in [Cytoscape](https://cytoscape.org/) (yFiles Organic Layout recommended) +- If `--tree`/`build_tree=True`: a MAFFT alignment, a FastTree phylogenetic tree (Newick), and [iTOL](https://itol.embl.de/) annotation files for the domain architecture and (for CAZy family queries) characterized-enzyme activity ## Acknowledgements @@ -146,13 +66,12 @@ CANDy communicates with and/or references the following separate libraries, pack - [Biopython](https://biopython.org/) - [pandas](https://pandas.pydata.org/) -- [tqdm](https://github.com/tqdm/tqdm) +- [tqdm](https://github.com/tqdm/tqdm) - [sqlitebrowser](https://sqlitebrowser.org/) - [SQLAlchemy](https://www.sqlalchemy.org/) -- [sdRDM](https://github.com/JR-1991/software-driven-rdm) - [requests](https://requests.readthedocs.io/en/latest/) -- [ete3](http://etetoolkit.org/) - [CD-HIT](https://academic.oup.com/bioinformatics/article/22/13/1658/194225?login=true) +- [MMseqs2](https://www.nature.com/articles/nbt.3988) - [MAFFT](https://academic.oup.com/nar/article/30/14/3059/2904316?login=true) - [FastTree](http://www.microbesonline.org/fasttree/) - [NetworkX](https://networkx.org/) @@ -168,9 +87,9 @@ Windels A, Franceus J, Pleiss J, Desmet T. CANDy: Automated analysis of domain a ### License and Disclaimer -This Jupyter Notebook is licensed under [MIT](https://github.com/PyEED/CANDy/blob/main/SECURITY.md#mit-license). +CANDy is licensed under [MIT](https://github.com/PyEED/CANDy/blob/main/SECURITY.md#mit-license). -This Notebook and other information provided is for theoretical utilisation only, caution should be exercised in its use. It is provided ‘as-is’ without any warranty of any kind, whether expressed or implied. Information is not intended to be a substitute for professional medical advice, diagnosis, or treatment, and does not constitute medical or other professional advice. +CANDy and other information provided is for theoretical utilisation only, caution should be exercised in its use. It is provided 'as-is' without any warranty of any kind, whether expressed or implied. Information is not intended to be a substitute for professional medical advice, diagnosis, or treatment, and does not constitute medical or other professional advice. ### Third-party software @@ -180,6 +99,6 @@ Use of the third-party software, libraries or code referred to in the [Acknowled The following databases are used by CANDy, and are available with reference to the following: - UniProt: (unmodified), by The UniProt Consortium, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). -- NCBI: (unmodified), by the National Library of Medicine, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). -- CAZy: (unmodified), by http://www.cazy.org/ and Elodie Drula, Marie-Line Garron, Suzan Dogan, Vincent Lombard, Bernard Henrissat, Nicolas Terrapon, The carbohydrate-active enzyme database: functions and literature, Nucleic Acids Research, Volume 50, Issue D1, 7 January 2022, Pages D571–D577, https://doi.org/10.1093/nar/gkab1045, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). -- InterPro: (unmodified), by EMBL-EBI, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). +- NCBI: (unmodified), by the National Library of Medicine, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). +- CAZy: (unmodified), by http://www.cazy.org/ and Elodie Drula, Marie-Line Garron, Suzan Dogan, Vincent Lombard, Bernard Henrissat, Nicolas Terrapon, The carbohydrate-active enzyme database: functions and literature, Nucleic Acids Research, Volume 50, Issue D1, 7 January 2022, Pages D571–D577, https://doi.org/10.1093/nar/gkab1045, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). +- InterPro: (unmodified), by EMBL-EBI, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/). diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..5084977 --- /dev/null +++ b/environment.yml @@ -0,0 +1,13 @@ +name: candy +channels: + - bioconda + - conda-forge +dependencies: + - python>=3.10 + - pip + - mmseqs2 + - cd-hit + - mafft + - fasttree + - pip: + - -e . diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..1b9b3ba --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,51 @@ +[build-system] +requires = ["setuptools>=68", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "candy-cazyme" +version = "3.0.0.dev0" +description = "Automated analysis of domain architectures in carbohydrate-active enzymes (CAZymes)" +readme = "README.md" +license = { text = "MIT" } +authors = [ + { name = "Alex Windels" }, +] +requires-python = ">=3.10" +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", + "Topic :: Scientific/Engineering :: Bio-Informatics", +] +dependencies = [ + "biopython>=1.81", + "pandas>=2.0", + "requests>=2.31", + "tqdm>=4.66", + "SQLAlchemy>=2.0", + "networkx>=3.1", + "matplotlib>=3.7", + "numpy>=1.26", + "typer>=0.12", + "lxml>=5.0", +] + +[project.optional-dependencies] +gemini = ["google-genai>=0.3"] +dev = [ + "pytest>=8.0", + "pytest-mock>=3.14", +] + +[project.scripts] +candy = "candy.cli:app" + +[tool.setuptools.packages.find] +where = ["src"] + +[tool.pytest.ini_options] +testpaths = ["tests"] +markers = [ + "integration: requires network access or external CLI tools (mafft, fasttree, mmseqs2, cd-hit)", +] diff --git a/src/candy/__init__.py b/src/candy/__init__.py new file mode 100644 index 0000000..0979834 --- /dev/null +++ b/src/candy/__init__.py @@ -0,0 +1,3 @@ +"""CANDy: automated analysis of domain architectures in carbohydrate-active enzymes.""" + +__version__ = "3.0.0.dev0" diff --git a/src/candy/alignment/__init__.py b/src/candy/alignment/__init__.py new file mode 100644 index 0000000..5312b5a --- /dev/null +++ b/src/candy/alignment/__init__.py @@ -0,0 +1,22 @@ +"""Pluggable multiple-sequence-alignment backends.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Protocol + + +class AlignmentTool(Protocol): + name: str + + def align(self, input_fasta: Path, output_fasta: Path) -> Path: + """Align sequences in ``input_fasta``, writing the alignment to ``output_fasta``.""" + ... + + +def get_alignment_tool(name: str) -> AlignmentTool: + if name == "mafft": + from candy.alignment.mafft import MafftAligner + + return MafftAligner() + raise ValueError(f"Unknown alignment tool: {name}") diff --git a/src/candy/alignment/mafft.py b/src/candy/alignment/mafft.py new file mode 100644 index 0000000..ccd931b --- /dev/null +++ b/src/candy/alignment/mafft.py @@ -0,0 +1,25 @@ +from __future__ import annotations + +import logging +from pathlib import Path + +from candy.external_tools import require_binary, run_tool + +logger = logging.getLogger(__name__) + + +class MafftAligner: + """Runs MAFFT directly via subprocess. + + The notebook used Biopython's ``Bio.Align.Applications.MafftCommandline`` + wrapper, which is deprecated upstream and slated for removal from + Biopython; calling the CLI directly avoids that dependency. + """ + + name = "mafft" + + def align(self, input_fasta: Path, output_fasta: Path) -> Path: + binary = require_binary("mafft") + logger.info("Aligning %s with MAFFT.", input_fasta) + run_tool([binary, str(input_fasta)], stdout_path=output_fasta) + return output_fasta diff --git a/src/candy/blast.py b/src/candy/blast.py new file mode 100644 index 0000000..a09b674 --- /dev/null +++ b/src/candy/blast.py @@ -0,0 +1,110 @@ +"""BLAST fallback for characterized enzymes without a UniParc entry. + +Characterized (experimentally studied) enzymes from CAZy are always +included in the analysis even if they weren't picked up by clustering. +Some of them, however, don't have a UniParc entry, so InterPro can't +annotate them directly; this module BLASTs those sequences against ``nr`` +to find a close homolog that *does* have a UniParc entry, and uses that +homolog as a stand-in for domain detection. + +.. note:: + The original notebook accumulated BLAST hits from *every* characterized + sequence processed so far into one shared dict, then re-ran the UniParc + check against that whole accumulated dict on each iteration -- wasteful + (repeated, growing InterPro queries), and worse, it could pick up + UniParc-available accessions left over from a *different* characterized + sequence's hits while the current sequence had none of its own, at which + point ``max()`` was called on an empty dict and crashed. This version + scopes each sequence's BLAST candidates and UniParc check to itself. +""" + +from __future__ import annotations + +import logging +from collections.abc import Mapping, Sequence +from pathlib import Path + +from Bio import SeqIO +from Bio.Blast import NCBIWWW, NCBIXML + +from candy.interpro import match_lookup + +logger = logging.getLogger(__name__) + +_EXPECT_THRESHOLD = 5.0 +_WORD_SIZE = 6 + + +def _blast_candidates(sequence, identity_threshold: float) -> dict[str, tuple[float, str]]: + """BLASTP a sequence against nr; returns {accession: (percent_identity, hit_sequence)}.""" + result_handle = NCBIWWW.qblast("blastp", "nr", sequence, expect=_EXPECT_THRESHOLD, word_size=_WORD_SIZE) + + candidates: dict[str, tuple[float, str]] = {} + for blast_record in NCBIXML.parse(result_handle): + query_length = blast_record.query_length + for alignment in blast_record.alignments: + for hsp in alignment.hsps: + identity = hsp.identities / query_length * 100 + if identity >= identity_threshold: + candidates[alignment.accession] = (round(identity, 2), hsp.sbjct) + return candidates + + +def resolve_characterized_via_blast( + characterized_fasta: str | Path, + ids_without_uniparc: Sequence[str], + identity_threshold: float, +) -> tuple[dict[str, str], dict[str, str]]: + """Find a UniParc-available BLAST homolog for each characterized sequence lacking one. + + Returns ``(chosen_hit_by_id, hit_sequences)``: for each resolved + characterized ``protein_id``, ``chosen_hit_by_id[protein_id]`` is the + accession of the best matching homolog, and ``hit_sequences[accession]`` + is that homolog's sequence (to be used in place of the characterized + sequence for domain detection). + """ + ids_without_uniparc = set(ids_without_uniparc) + chosen_hit_by_id: dict[str, str] = {} + hit_sequences: dict[str, str] = {} + + with open(characterized_fasta) as handle: + records = list(SeqIO.parse(handle, "fasta")) + + for record in records: + protein_id = record.id.split("_")[0] + if protein_id not in ids_without_uniparc: + continue + + logger.info("Doing a BLASTP search for %s.", protein_id) + candidates = _blast_candidates(record.seq, identity_threshold) + + if not candidates: + logger.warning( + "No BLAST hits above %.0f%% identity for %s; the sequence will be excluded.", + identity_threshold, protein_id, + ) + continue + + candidate_sequences = {accession: seq for accession, (_, seq) in candidates.items()} + matches, _ = match_lookup(candidate_sequences) + + in_uniparc = { + accession: identity + for accession, (identity, _) in candidates.items() + if accession in matches + } + if not in_uniparc: + logger.warning( + "No UniParc-available BLAST match found for %s; the sequence will be excluded.", protein_id + ) + continue + + best_accession = max(in_uniparc, key=in_uniparc.get) + chosen_hit_by_id[protein_id] = best_accession + hit_sequences[best_accession] = candidates[best_accession][1] + logger.info( + "For characterized sequence %s, %s (%.2f%% identical) will be used for domain detection.", + protein_id, best_accession, in_uniparc[best_accession], + ) + + return chosen_hit_by_id, hit_sequences diff --git a/src/candy/cazy.py b/src/candy/cazy.py new file mode 100644 index 0000000..54e4058 --- /dev/null +++ b/src/candy/cazy.py @@ -0,0 +1,217 @@ +"""Retrieve CAZy family sequences and characterized-enzyme data from CAZy/NCBI. + +The notebook fetched sequences one GenBank ID at a time via +``Entrez.efetch`` in a tight loop, which is slow and easy to get +rate-limited on for large families. Entrez supports batched fetches (a +comma-separated ``id`` list per request), used here instead; the resulting +FASTA output is identical, just retrieved in ~200-id chunks rather than +one request per sequence. +""" + +from __future__ import annotations + +import io +import logging +import time +from collections.abc import Iterable, Sequence +from urllib.error import HTTPError + +import pandas as pd +import requests +from Bio import Entrez + +from candy.config import Taxonomy + +logger = logging.getLogger(__name__) + +CAZY_FAMILY_URL = "https://www.cazy.org/IMG/cazy_data/{family}.txt" +CAZY_CHARACTERIZED_URL = "https://www.cazy.org/{family}_characterized.html" + +_TAXONOMY_NAME_TO_CODE = { + "BACTERIA": "B", + "EUKARYOTA": "E", + "ARCHAEA": "A", + "VIRUSES": "V", + "UNCLASSIFIED": "U", +} + +_ENTREZ_BATCH_SIZE = 200 +_ENTREZ_MAX_RETRIES = 10 +_ENTREZ_RETRY_DELAY = 3.0 + + +def _batched(items: Sequence[str], size: int) -> Iterable[Sequence[str]]: + for i in range(0, len(items), size): + yield items[i : i + size] + + +def fetch_family_page(family: str) -> str: + """Download the raw CAZy family data page (tab-separated GenBank ID listing).""" + url = CAZY_FAMILY_URL.format(family=family) + logger.info("Retrieving data from %s", url) + response = requests.get(url, timeout=60) + response.raise_for_status() + return response.text + + +def parse_family_table(raw_text: str) -> tuple[pd.DataFrame, dict[str, str]]: + """Parse a CAZy family page into (table, {genbank_id: taxonomy_code}).""" + df = pd.read_table(io.StringIO(raw_text), engine="python", header=None) + + taxonomy_dict: dict[str, str] = {} + for _, row in df.iterrows(): + kingdom = row[1].upper() + genbank_id = row[3] + taxonomy_dict[genbank_id] = _TAXONOMY_NAME_TO_CODE[kingdom] + + return df, taxonomy_dict + + +def extract_ncbi_ids(df: pd.DataFrame, taxonomy_dict: dict[str, str], taxonomy: Taxonomy) -> list[str]: + """Return GenBank IDs sourced from NCBI (excludes JGI-only entries), filtered by taxonomy.""" + is_ncbi = df.iloc[:, -1].astype(str) == "ncbi" + ids = df.iloc[:, -2][is_ncbi].tolist() + + if taxonomy == Taxonomy.ALL: + return ids + return [i for i in ids if taxonomy_dict.get(i) == taxonomy.code] + + +def fetch_sequences_fasta(ids: Sequence[str], email: str, batch_size: int = _ENTREZ_BATCH_SIZE) -> str: + """Batch-fetch protein FASTA sequences from NCBI for a list of GenBank IDs.""" + Entrez.email = email + chunks: list[str] = [] + + total_batches = (len(ids) + batch_size - 1) // batch_size or 1 + for batch in _batched(ids, batch_size): + attempts = 0 + while True: + try: + with Entrez.efetch(db="protein", id=",".join(batch), rettype="fasta", retmode="text") as handle: + chunks.append(handle.read()) + break + except HTTPError: + attempts += 1 + if attempts >= _ENTREZ_MAX_RETRIES: + logger.warning("Giving up on a batch of %d IDs after repeated HTTP errors.", len(batch)) + break + time.sleep(_ENTREZ_RETRY_DELAY) + logger.info("Fetched %d/%d sequence batches.", len(chunks), total_batches) + + return "".join(chunks) + + +def extract_family_sequences(family: str, taxonomy: Taxonomy, email: str) -> tuple[str, dict[str, str]]: + """End-to-end: fetch a CAZy family page and the protein sequences it lists. + + Returns (fasta_text, taxonomy_dict). Raises ``ValueError`` if no + sequences exist for the requested taxonomy subset. + """ + raw_text = fetch_family_page(family) + df, taxonomy_dict = parse_family_table(raw_text) + + if taxonomy != Taxonomy.ALL and taxonomy.code not in taxonomy_dict.values(): + raise ValueError( + f"No sequences belonging to {taxonomy.value} were found for family {family}." + ) + + ids = extract_ncbi_ids(df, taxonomy_dict, taxonomy) + logger.info("Extracting %d protein sequences from NCBI.", len(ids)) + fasta_text = fetch_sequences_fasta(ids, email) + + return fasta_text, taxonomy_dict + + +def fetch_characterized_page(family: str, max_retries: int = 10, retry_delay: float = 5.0) -> list[pd.DataFrame]: + """Download and parse CAZy's '{family}_characterized.html' page. + + CAZy's servers frequently return 429s under load, so this retries. + """ + url = CAZY_CHARACTERIZED_URL.format(family=family) + logger.info("Retrieving data from %s", url) + + attempts = 0 + while True: + try: + return pd.read_html(url) + except HTTPError: + attempts += 1 + if attempts >= max_retries: + raise + logger.info("Too many requests to CAZy; retrying in %.0fs.", retry_delay) + time.sleep(retry_delay) + + +def extract_characterized_ids( + tables: list[pd.DataFrame], +) -> tuple[list[str], dict[str, str], dict[str, str]]: + """Parse the characterized-enzyme table into (genbank_ids, {id: taxonomy_code}, {id: ec_number}). + + Column 4 is the GenBank identifier, column 1 is the EC number. + """ + activity_table = tables[1] + + activity_dict: dict[str, str] = {} + for i in range(len(activity_table)): + identifier = str(activity_table[4].loc[i]).split(" ")[0] + if "." in identifier: + activity_dict[identifier] = activity_table[1].loc[i] + + id_series = tables[1][4].drop_duplicates() + + taxonomy_dict: dict[str, str] = {} + ids: list[str] = [] + current_taxonomy = "" + taxonomy_headers = {"Archaea": "A", "Bacteria": "B", "Eukaryota": "E", "Viruses": "V", "Unclassified": "U"} + + for raw_id in id_series: + entry = str(raw_id) + if entry in taxonomy_headers: + current_taxonomy = taxonomy_headers[entry] + continue + + dot_position = entry.find(".") + if dot_position == -1: + continue + + genbank_id = entry[: dot_position + 2] + if genbank_id not in ids: + ids.append(genbank_id) + taxonomy_dict[genbank_id] = current_taxonomy + + return ids, taxonomy_dict, activity_dict + + +def fetch_characterized_sequences( + ids: Sequence[str], taxonomy_dict: dict[str, str], taxonomy: Taxonomy, email: str +) -> str: + """Fetch characterized-enzyme FASTA sequences, annotated with organism + taxonomy code.""" + Entrez.email = email + filtered_ids = [i for i in ids if taxonomy == Taxonomy.ALL or taxonomy_dict.get(i) == taxonomy.code] + + output = "" + for genbank_id in filtered_ids: + try: + with Entrez.efetch(db="protein", id=genbank_id, rettype="fasta", retmode="text") as handle: + text = handle.read() + except HTTPError: + logger.warning("HTTP error fetching characterized sequence %s.", genbank_id) + continue + + from Bio import SeqIO + + record = next(SeqIO.parse(io.StringIO(text), "fasta"), None) + if record is None: + logger.warning("Could not parse characterized sequence %s.", genbank_id) + continue + + # The notebook computed this slice index from str(seq_record) (the + # multi-line repr, not the id itself), via `.id[0:str(seq_record).find('.')+1]`. + # That offset always overshoots the id's real length, so the slice + # was a no-op in every realistic case -- it just returns the full id. + accession = record.id + organism = record.description[record.description.find("[") + 1 : record.description.find("]")] + header = f">{accession} {organism}_{taxonomy_dict.get(genbank_id, '')}".replace(" ", "_") + output += f"{header}\n{record.seq}\n" + + return output diff --git a/src/candy/cli.py b/src/candy/cli.py new file mode 100644 index 0000000..0976710 --- /dev/null +++ b/src/candy/cli.py @@ -0,0 +1,129 @@ +"""Command-line entry point: `candy run ...`.""" + +from __future__ import annotations + +import logging +import re +from pathlib import Path +from typing import Optional + +import typer + +from candy.config import ( + CAZyFamilyInput, + ClusteringConfig, + ClusteringSoftware, + CurationConfig, + CustomFastaInput, + DomainCleaningConfig, + PipelineConfig, + Taxonomy, +) +from candy.pipeline import run_pipeline + +app = typer.Typer(help="CANDy: automated analysis of domain architectures in carbohydrate-active enzymes.") + +_FAMILY_RE = re.compile(r"^([A-Za-z]+)(\d+)(?:_(\d+))?$") + + +@app.callback() +def _main() -> None: + """CANDy: automated analysis of domain architectures in carbohydrate-active enzymes. + + A no-op callback: Typer collapses a Typer() app with exactly one + @app.command() into a single top-level command (dropping the command + name entirely), which would silently break `candy run ...`. Registering + this callback keeps `run` as a real, required subcommand. + """ + + +def _parse_family(family: str) -> tuple[str, int, str | None]: + match = _FAMILY_RE.match(family) + if not match: + raise typer.BadParameter( + f"Could not parse '{family}' as a CAZy family (expected e.g. 'GH5' or 'GH5_1')." + ) + enzyme_class, number, subfamily = match.groups() + return enzyme_class, int(number), subfamily + + +@app.command() +def run( + jobname: str = typer.Option(..., help="Job name; results are written to output-dir/jobname."), + family: Optional[str] = typer.Option( + None, help="CAZy family to query, e.g. 'GH5' or 'GH5_1'. Mutually exclusive with --fasta." + ), + fasta: Optional[Path] = typer.Option( + None, exists=True, help="Custom FASTA file to analyse instead of querying CAZy." + ), + email: Optional[str] = typer.Option(None, help="Email for NCBI Entrez. Required with --family."), + taxonomy: Taxonomy = typer.Option(Taxonomy.ALL, help="Taxonomic subset to restrict a --family query to."), + output_dir: Path = typer.Option(Path("."), help="Directory results are written under."), + clustering_software: ClusteringSoftware = typer.Option( + ClusteringSoftware.MMSEQS2, help="Sequence clustering backend." + ), + cluster_identity: int = typer.Option(85, help="Clustering identity cutoff, percent."), + blast_identity: int = typer.Option( + 95, help="Minimum identity for the characterized-enzyme BLAST fallback, percent." + ), + max_domain_length: int = typer.Option(800, help="Maximum plausible domain length, amino acids."), + overlap_percentage: int = typer.Option( + 20, help="Overlap threshold above which two hits are considered the same domain, percent." + ), + build_tree: bool = typer.Option(False, "--tree/--no-tree", help="Run MSA + phylogenetics + iTOL export."), + curation_backend: str = typer.Option( + "manual", help="Domain-name curation backend: 'manual' or 'gemini'." + ), + curation_api_key: Optional[str] = typer.Option( + None, envvar="GOOGLE_API_KEY", help="API key for the curation backend, if it needs one." + ), + verbose: bool = typer.Option(False, "--verbose", "-v", help="Enable debug logging."), +) -> None: + """Run the CANDy pipeline end-to-end.""" + logging.basicConfig(level=logging.DEBUG if verbose else logging.INFO, format="%(message)s") + + if (family is None) == (fasta is None): + raise typer.BadParameter("Specify exactly one of --family or --fasta.") + + if family is not None: + if not email: + raise typer.BadParameter("--email is required when using --family (NCBI Entrez requires it).") + enzyme_class, family_number, subfamily = _parse_family(family) + pipeline_input = CAZyFamilyInput( + enzyme_class=enzyme_class, + family_number=family_number, + subfamily=subfamily, + email=email, + taxonomy=taxonomy, + ) + else: + pipeline_input = CustomFastaInput(fasta_path=fasta) + + config = PipelineConfig( + input=pipeline_input, + jobname=jobname, + output_dir=output_dir, + clustering=ClusteringConfig(software=clustering_software, identity_cutoff=cluster_identity), + domain_cleaning=DomainCleaningConfig( + max_domain_length=max_domain_length, overlap_percentage=overlap_percentage + ), + curation=CurationConfig(backend=curation_backend, api_key=curation_api_key), + blast_identity_threshold=blast_identity, + build_tree=build_tree, + ) + + result = run_pipeline(config) + + typer.echo(f"\nDone. Results written to {result.jobname_dir}") + typer.echo(f" Database: {result.database_path}") + typer.echo(f" Co-occurrence network: {result.network_graphml_path}") + if result.tree_path: + typer.echo(f" Alignment: {result.alignment_path}") + typer.echo(f" Tree: {result.tree_path}") + typer.echo(f" Domain annotation: {result.domain_annotation_path}") + if result.characterized_annotation_path: + typer.echo(f" Characterized annotation: {result.characterized_annotation_path}") + + +if __name__ == "__main__": + app() diff --git a/src/candy/clustering/__init__.py b/src/candy/clustering/__init__.py new file mode 100644 index 0000000..d6b65a7 --- /dev/null +++ b/src/candy/clustering/__init__.py @@ -0,0 +1,36 @@ +"""Pluggable sequence-clustering backends. + +CANDy reduces redundancy in the input set before domain detection, since it +cares about domain *architectures* rather than exact sequence identity, so a +fairly permissive clustering cutoff is typically used. Both tools the +published notebook offered (CD-HIT and MMseqs2) are kept as interchangeable +:class:`Clusterer` implementations; adding another tool later only means +implementing this protocol. +""" + +from __future__ import annotations + +from pathlib import Path +from typing import Protocol + +from candy.config import ClusteringConfig, ClusteringSoftware + + +class Clusterer(Protocol): + name: str + + def cluster(self, input_fasta: Path, output_fasta: Path, config: ClusteringConfig) -> Path: + """Cluster ``input_fasta`` and write representative sequences to ``output_fasta``.""" + ... + + +def get_clusterer(software: ClusteringSoftware) -> Clusterer: + if software == ClusteringSoftware.CDHIT: + from candy.clustering.cdhit import CdHitClusterer + + return CdHitClusterer() + if software == ClusteringSoftware.MMSEQS2: + from candy.clustering.mmseqs2 import Mmseqs2Clusterer + + return Mmseqs2Clusterer() + raise ValueError(f"Unknown clustering software: {software}") diff --git a/src/candy/clustering/cdhit.py b/src/candy/clustering/cdhit.py new file mode 100644 index 0000000..f7ffafb --- /dev/null +++ b/src/candy/clustering/cdhit.py @@ -0,0 +1,42 @@ +from __future__ import annotations + +import logging +from pathlib import Path + +from candy.config import ClusteringConfig +from candy.external_tools import require_binary, run_tool + +logger = logging.getLogger(__name__) + +# CD-HIT's word-length and memory-limit flags aren't exposed as pipeline +# parameters (the notebook hardcoded them too); they only affect runtime, +# not clustering results. +_WORD_LENGTH = 5 +_MEMORY_LIMIT_MB = 8000 +_DESCRIPTION_LENGTH = 20 + + +class CdHitClusterer: + name = "cd-hit" + + def cluster(self, input_fasta: Path, output_fasta: Path, config: ClusteringConfig) -> Path: + binary = require_binary("cd-hit") + identity = config.identity_cutoff / 100 + min_long_coverage = config.cdhit_min_long_coverage / 100 + min_short_coverage = config.cdhit_min_short_coverage / 100 + + logger.info("Clustering with CD-HIT at %.0f%% identity.", config.identity_cutoff) + run_tool( + [ + binary, + "-i", str(input_fasta), + "-o", str(output_fasta), + "-c", f"{identity}", + "-n", str(_WORD_LENGTH), + "-d", str(_DESCRIPTION_LENGTH), + "-M", str(_MEMORY_LIMIT_MB), + "-aL", f"{min_long_coverage}", + "-aS", f"{min_short_coverage}", + ] + ) + return output_fasta diff --git a/src/candy/clustering/mmseqs2.py b/src/candy/clustering/mmseqs2.py new file mode 100644 index 0000000..bca212f --- /dev/null +++ b/src/candy/clustering/mmseqs2.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +import logging +import shutil +import tempfile +from pathlib import Path + +from candy.config import ClusteringConfig +from candy.external_tools import require_binary, run_tool + +logger = logging.getLogger(__name__) + + +class Mmseqs2Clusterer: + name = "mmseqs2" + + def cluster(self, input_fasta: Path, output_fasta: Path, config: ClusteringConfig) -> Path: + binary = require_binary("mmseqs") + min_seq_id = config.identity_cutoff / 100 + min_coverage = config.mmseqs_min_coverage / 100 + + logger.info("Clustering with MMseqs2 at %.0f%% identity.", config.identity_cutoff) + with tempfile.TemporaryDirectory(prefix="candy_mmseqs2_") as tmp_dir: + tmp_dir = Path(tmp_dir) + output_prefix = tmp_dir / "cluster" + run_tool( + [ + binary, + "easy-cluster", + str(input_fasta), + str(output_prefix), + str(tmp_dir / "tmp"), + "--min-seq-id", f"{min_seq_id}", + "-c", f"{min_coverage}", + "--cov-mode", str(config.mmseqs_cov_mode), + ] + ) + representative_sequences = Path(f"{output_prefix}_rep_seq.fasta") + shutil.copyfile(representative_sequences, output_fasta) + + return output_fasta diff --git a/src/candy/config.py b/src/candy/config.py new file mode 100644 index 0000000..18d5f13 --- /dev/null +++ b/src/candy/config.py @@ -0,0 +1,143 @@ +"""Explicit, testable run configuration. + +The original Colab notebook collected these parameters through ``ipywidgets`` +dropdowns and blocking ``input()`` calls, and threaded them through the +pipeline as module-level globals (``family``, ``go``, ``jobname``, ...) that +later cells silently depended on. Here they are collected into plain +dataclasses that can be constructed directly (Python API), parsed from CLI +flags, or built in a test without touching stdin. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from enum import Enum +from pathlib import Path + + +class Taxonomy(str, Enum): + ALL = "All" + ARCHAEA = "Archaea" + BACTERIA = "Bacteria" + EUKARYOTA = "Eukaryota" + VIRUSES = "Viruses" + UNCLASSIFIED = "Unclassified" + + @property + def code(self) -> str: + """Single-letter code used throughout CAZy/CANDy output, e.g. 'B' for Bacteria.""" + return { + Taxonomy.ALL: "All", + Taxonomy.ARCHAEA: "A", + Taxonomy.BACTERIA: "B", + Taxonomy.EUKARYOTA: "E", + Taxonomy.VIRUSES: "V", + Taxonomy.UNCLASSIFIED: "U", + }[self] + + +class ClusteringSoftware(str, Enum): + CDHIT = "cd-hit" + MMSEQS2 = "mmseqs2" + + +# Default ranking (lower = higher priority) used to resolve overlapping domain +# annotations sourced from different InterPro member databases. Mirrors the +# notebook's default so results are reproducible across the rewrite. +DEFAULT_DATABASE_PREFERENCE: list[str] = [ + "SMART", + "CDD", + "PFAM", + "SUPERFAMILY", + "NCBIFAM", + "PANTHER", + "GENE3D", + "PIRSF", + "HAMAP", + "PRINTS", + "SFLD", + "PROFILE", + "PROSITE", +] + + +@dataclass +class CAZyFamilyInput: + """Query a CAZy family/subfamily directly (input option 1 in the notebook).""" + + enzyme_class: str # "GH", "GT", "PL", "CE", "AA" + family_number: int + email: str + subfamily: str | None = None + taxonomy: Taxonomy = Taxonomy.ALL + + @property + def family(self) -> str: + base = f"{self.enzyme_class}{self.family_number}" + return f"{base}_{self.subfamily}" if self.subfamily else base + + +@dataclass +class CustomFastaInput: + """Analyse a user-supplied FASTA file directly (input option 2 in the notebook).""" + + fasta_path: Path + + def __post_init__(self) -> None: + self.fasta_path = Path(self.fasta_path) + + +@dataclass +class ClusteringConfig: + software: ClusteringSoftware = ClusteringSoftware.MMSEQS2 + identity_cutoff: int = 85 # percent + + # CD-HIT specific + cdhit_min_short_coverage: int = 90 # -aS, percent + cdhit_min_long_coverage: int = 90 # -aL, percent + + # MMseqs2 specific + mmseqs_cov_mode: int = 0 + mmseqs_min_coverage: int = 90 # percent + + +@dataclass +class DomainCleaningConfig: + max_domain_length: int = 800 # amino acids + overlap_percentage: int = 20 + database_preference: list[str] = field( + default_factory=lambda: list(DEFAULT_DATABASE_PREFERENCE) + ) + + +@dataclass +class CurationConfig: + """Which backend resolves synonymous domain names into one umbrella name.""" + + backend: str = "manual" # "manual" | "gemini" | any registered CurationBackend name + api_key: str | None = None + + +@dataclass +class PipelineConfig: + input: CAZyFamilyInput | CustomFastaInput + jobname: str + output_dir: Path + + clustering: ClusteringConfig = field(default_factory=ClusteringConfig) + domain_cleaning: DomainCleaningConfig = field(default_factory=DomainCleaningConfig) + curation: CurationConfig = field(default_factory=CurationConfig) + + blast_identity_threshold: int = 95 # percent; used only for the characterized-enzyme + # BLAST fallback when a characterized sequence has no UniParc entry. + + build_tree: bool = False + alignment_tool: str = "mafft" + tree_tool: str = "fasttree" + + def __post_init__(self) -> None: + self.output_dir = Path(self.output_dir) + + @property + def is_cazy_query(self) -> bool: + return isinstance(self.input, CAZyFamilyInput) diff --git a/src/candy/curation/__init__.py b/src/candy/curation/__init__.py new file mode 100644 index 0000000..63beb7d --- /dev/null +++ b/src/candy/curation/__init__.py @@ -0,0 +1,37 @@ +"""Pluggable domain-name curation backends. + +InterPro's member databases often have several different names for what is +functionally the same domain. Curation groups those synonymous names under +one umbrella name. The notebook hard-required a Gemini API key for its +"automated" path; here that's one optional backend among others behind a +common protocol, with manual curation as the dependency-free default. +""" + +from __future__ import annotations + +from typing import Protocol + + +class CurationBackend(Protocol): + name: str + + def curate(self, domain_names: list[str], *, family: str | None = None) -> dict[str, list[str]]: + """Group synonymous domain names. + + Returns ``{umbrella_name: [raw_name, ...]}``. Every name in + ``domain_names`` should appear in exactly one group (domains that + can't be grouped with anything else become their own singleton group). + """ + ... + + +def get_curation_backend(name: str, **kwargs) -> CurationBackend: + if name == "manual": + from candy.curation.manual import ManualCurationBackend + + return ManualCurationBackend(**kwargs) + if name == "gemini": + from candy.curation.gemini import GeminiCurationBackend + + return GeminiCurationBackend(**kwargs) + raise ValueError(f"Unknown curation backend: {name}") diff --git a/src/candy/curation/gemini.py b/src/candy/curation/gemini.py new file mode 100644 index 0000000..f9a2592 --- /dev/null +++ b/src/candy/curation/gemini.py @@ -0,0 +1,58 @@ +from __future__ import annotations + +import ast +import os + +_PROMPT_TEMPLATE = ( + "I have the following list of protein domains, retrieved from InterPro's member " + "databases: {domain_names}. These sequences {family_clause}. I want you to group " + "these domains in overarching domains. For instance, all domains related to the " + "catalytic activity need to be grouped under the term 'Catalytic domain', domains " + "related to carbohydrate binding need to be grouped as 'CBM', and so on. Fibronectin " + "III domains can be considered as immunoglobulin-like. Determine the overarching " + "names yourself. If a domain can't be grouped with another one, keep the name. If a " + "catalytic domain is found that is not expected for this family, don't include it in " + "the 'Catalytic domain' name, but give it a separate name. I want you to output only a " + "dictionary in the following format: {{'overarching name': ['domain 1', 'domain 2', " + "'domain 3'], ... }}" +) + + +class GeminiCurationBackend: + """Automated curation via Google Gemini, ported from the notebook's default path. + + Requires the optional ``gemini`` extra (``pip install candy-cazyme[gemini]``) + and an API key, passed explicitly or read from ``GOOGLE_API_KEY`` / + ``GEMINI_API_KEY``. + """ + + name = "gemini" + + def __init__(self, api_key: str | None = None, model: str = "gemini-2.0-flash") -> None: + self.api_key = api_key or os.environ.get("GOOGLE_API_KEY") or os.environ.get("GEMINI_API_KEY") + if not self.api_key: + raise ValueError( + "Gemini curation requires an API key. Pass api_key=, or set the " + "GOOGLE_API_KEY or GEMINI_API_KEY environment variable." + ) + self.model = model + + def curate(self, domain_names: list[str], *, family: str | None = None) -> dict[str, list[str]]: + try: + from google import genai + except ImportError as exc: + raise ImportError( + "Gemini curation requires the optional 'gemini' extra: " + "pip install candy-cazyme[gemini]" + ) from exc + + family_clause = ( + f"belong to enzymes from CAZy family {family}" if family else "are carbohydrate-active enzymes" + ) + prompt = _PROMPT_TEMPLATE.format(domain_names=domain_names, family_clause=family_clause) + + client = genai.Client(api_key=self.api_key) + response = client.models.generate_content(model=self.model, contents=prompt) + + cleaned = response.text.replace("python", "").replace("```", "") + return ast.literal_eval(cleaned) diff --git a/src/candy/curation/manual.py b/src/candy/curation/manual.py new file mode 100644 index 0000000..726f5e3 --- /dev/null +++ b/src/candy/curation/manual.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +from collections.abc import Callable + + +class ManualCurationBackend: + """Interactively prompt the user to group synonymous domain names. + + ``input_fn``/``print_fn`` are injectable so this can be driven + programmatically in tests without touching real stdin/stdout. + """ + + name = "manual" + + def __init__( + self, + input_fn: Callable[[str], str] = input, + print_fn: Callable[[str], None] = print, + ) -> None: + self._input = input_fn + self._print = print_fn + + def curate(self, domain_names: list[str], *, family: str | None = None) -> dict[str, list[str]]: + remaining = dict(enumerate(domain_names)) + curated: dict[str, list[str]] = {} + + while remaining: + self._print("\nDomains still to curate:") + for index, name in remaining.items(): + self._print(f"{index}: {name}") + + umbrella_name = self._input("\nDomain name: ").strip() + if umbrella_name.upper() == "STOP": + for name in remaining.values(): + curated[name] = [name] + break + + indices = [int(i) for i in self._input("Includes: ").strip().split(",")] + curated[umbrella_name] = [remaining.pop(i) for i in indices] + + return curated diff --git a/src/candy/database.py b/src/candy/database.py new file mode 100644 index 0000000..eb8aa58 --- /dev/null +++ b/src/candy/database.py @@ -0,0 +1,156 @@ +"""SQLite result database, via SQLAlchemy. + +The notebook listed SQLAlchemy as a dependency but never actually used it -- +the database was built with raw ``sqlite3`` calls. This module replaces that +with real SQLAlchemy models. + +While porting the insert logic, the original ``populate_database`` was found +to write ``(id, organism, characterized, taxonomy, sequence, ...)`` into a +table declared as ``(id, taxonomy, characterized, organism_name, sequence, +...)`` -- i.e. the ``taxonomy`` and ``organism_name`` columns were swapped, so +every exported database had organism names (e.g. "Escherichia_coli") stored +under the "taxonomy" column and kingdom codes (e.g. "Bacteria") under +"organism_name". Fixed here; flagged since anyone with tooling built against +the old (mislabeled) column semantics will need to adjust. +""" + +from __future__ import annotations + +import logging +from collections.abc import Mapping, Sequence +from pathlib import Path + +from sqlalchemy import String, Text, create_engine +from sqlalchemy.orm import DeclarativeBase, Mapped, Session, mapped_column + +logger = logging.getLogger(__name__) + +ARCHITECTURE_SEPARATOR = "--" + +TAXONOMY_CODE_TO_NAME = { + "B": "Bacteria", + "A": "Archaea", + "E": "Eukaryota", + "V": "Viruses", + "U": "Unclassified", +} + + +class Base(DeclarativeBase): + pass + + +class ProteinSequenceRow(Base): + __tablename__ = "protein_sequences" + + protein_sequence_id: Mapped[str] = mapped_column(String, primary_key=True) + taxonomy: Mapped[str] = mapped_column(String, default="") + characterized: Mapped[str] = mapped_column(String, default="") + organism_name: Mapped[str] = mapped_column(String, default="") + amino_acid_sequence: Mapped[str] = mapped_column(Text, default="") + domain_architecture: Mapped[str] = mapped_column(Text, default="") + domain_locations: Mapped[str] = mapped_column(Text, default="") + domain_databases: Mapped[str] = mapped_column(Text, default="") + + +class DomainAssemblyRow(Base): + __tablename__ = "domain_assemblies" + + protein_domain: Mapped[str] = mapped_column(Text, primary_key=True) + protein_sequence_ids: Mapped[str] = mapped_column(Text, default="") + + +class DomainCurationRow(Base): + __tablename__ = "domain_curation" + + domain_name: Mapped[str] = mapped_column(Text, primary_key=True) + synonymous_domain_names: Mapped[str] = mapped_column(Text, default="") + + +def create_database(db_path: str | Path): + """Create (or open) the SQLite database and return a SQLAlchemy engine.""" + engine = create_engine(f"sqlite:///{db_path}") + Base.metadata.create_all(engine) + return engine + + +def populate_database( + engine, + proteins: Sequence, # Bio.SeqRecord.SeqRecord, id formatted as "{accession}_{organism}_{taxcode}" + characterized_ids: Mapping[str, str], + domain_architecture: Mapping[str, Sequence[tuple[list[int], str]]], + domain_to_database: Mapping[str, Mapping[str, str]], + curated_domains: Mapping[str, Sequence[str]], + *, + is_cazy_query: bool, +) -> None: + """Populate all three tables from the pipeline's intermediate results. + + ``domain_architecture`` and ``domain_to_database`` are the per-protein + outputs of :func:`candy.domains.clean_domains` and + :func:`candy.interpro.build_domain_tables` respectively. ``curated_domains`` + maps an umbrella domain name to the list of raw names folded into it. + """ + name_to_umbrella = { + raw_name: umbrella for umbrella, raw_names in curated_domains.items() for raw_name in raw_names + } + + architecture_to_ids: dict[str, list[str]] = {} + seen_ids: set[str] = set() + rows: list[ProteinSequenceRow] = [] + + for record in proteins: + if is_cazy_query: + protein_id = record.id.split("_")[0] + organism = " ".join(record.id.split("_")[1:-1]) + taxonomy = TAXONOMY_CODE_TO_NAME.get(record.id[-1], "") + else: + protein_id = record.id + organism = "" + taxonomy = "" + + if protein_id in seen_ids: + logger.warning("Skipping duplicate sequence ID: %s", protein_id) + continue + seen_ids.add(protein_id) + + domain_names = [] + domain_positions = [] + domain_databases = [] + for position, raw_name in domain_architecture.get(protein_id, []): + domain_names.append(name_to_umbrella.get(raw_name, raw_name)) + domain_positions.append(str(position)) + domain_databases.append(domain_to_database.get(protein_id, {}).get(raw_name, "")) + + architecture = ARCHITECTURE_SEPARATOR.join(domain_names) + architecture_to_ids.setdefault(architecture, []).append(protein_id) + + rows.append( + ProteinSequenceRow( + protein_sequence_id=protein_id, + taxonomy=taxonomy, + characterized="C" if protein_id in characterized_ids else "", + organism_name=organism, + amino_acid_sequence=str(record.seq), + domain_architecture=architecture, + domain_locations=ARCHITECTURE_SEPARATOR.join(domain_positions), + domain_databases=ARCHITECTURE_SEPARATOR.join(domain_databases), + ) + ) + + with Session(engine) as session: + session.add_all(rows) + session.add_all( + DomainAssemblyRow(protein_domain=architecture, protein_sequence_ids=", ".join(ids)) + for architecture, ids in architecture_to_ids.items() + ) + session.add_all( + DomainCurationRow(domain_name=umbrella, synonymous_domain_names=", ".join(raw_names)) + for umbrella, raw_names in curated_domains.items() + ) + session.commit() + + +def read_protein_sequences(engine) -> list[ProteinSequenceRow]: + with Session(engine) as session: + return list(session.query(ProteinSequenceRow).all()) diff --git a/src/candy/domains.py b/src/candy/domains.py new file mode 100644 index 0000000..6cd3c5c --- /dev/null +++ b/src/candy/domains.py @@ -0,0 +1,246 @@ +"""Resolve overlapping/duplicate domain annotations into a clean architecture. + +InterPro aggregates hits from many member databases (Pfam, SMART, CDD, ...), +so the same physical domain is usually reported multiple times, at slightly +different boundaries, by different databases. This module groups positions +that overlap "enough" to be the same domain, then picks one representative +hit per group according to a configurable database-preference ranking. + +Ported from the notebook's ``overlap`` / ``extract_overlap`` / ``clean_domains`` +functions with the same resolution semantics (including its pre-existing +asymmetry: the single-hit branch additionally excludes the FUNFAM database +and does not length-filter, while the ambiguous/grouped branch does not +exclude FUNFAM but does length-filter). That asymmetry is preserved rather +than "fixed" here to keep behaviour identical to the published tool. +""" + +from __future__ import annotations + +import itertools +from collections.abc import Mapping, Sequence + +Position = list[int] # [start, end], 1 excludes 0-based confusion by matching notebook convention +NamedDomain = tuple[Position, str] # ([start, end], domain_name) + +_EXCLUDED_DOMAIN_NAMES = { + "None", + "Prokaryotic membrane lipoprotein lipid attachment site profile", +} +_EXCLUDED_DATABASES_SINGLE = {"TMHMM", "MOBIDB_LITE", "PHOBIUS", "FUNFAM"} +_EXCLUDED_DATABASES_GROUPED = {"TMHMM", "MOBIDB_LITE", "PHOBIUS"} + + +def positions_overlap(a: Position, b: Position, overlap_fraction: float) -> bool: + """True if two [start, end] positions overlap by more than ``overlap_fraction``.""" + a, b = sorted([a, b]) + x1, y1 = a + x2, y2 = b + length1 = y1 - x1 + length2 = y2 - x2 + + if x2 <= y1 <= y2: + if x1 == y1 or x2 == y2: + identity = 1.0 + else: + over = y1 - x2 + identity = max(over / length1, over / length2) + elif x1 <= x2 and y2 <= y1: + # b falls entirely within a + identity = 100.0 + else: + identity = 0.0 + + return identity > overlap_fraction + + +def group_overlapping( + positions: Sequence[Position], overlap_fraction: float +) -> list[Position | list[Position]]: + """Group positions considered to represent the same domain. + + Returns a list where each element is either a single ``[start, end]`` + (unambiguous) or a list of ``[start, end]`` pairs considered to be the + same domain reported at slightly different boundaries. + """ + ordered = sorted(positions, key=lambda p: (p[0], p[1])) + grouped: list[Position | list[Position]] = [] + pending: list[Position] = [] + + i = 0 + while i < len(ordered): + if i == len(ordered) - 1: + if pending: + grouped.append(_dedupe(pending)) + else: + grouped.append(ordered[i]) + break + + if positions_overlap(ordered[i], ordered[i + 1], overlap_fraction): + pending += [ordered[i], ordered[i + 1]] + i += 1 + else: + if pending: + grouped.append(_dedupe(pending)) + pending = [] + else: + grouped.append(ordered[i]) + i += 1 + + return grouped + + +def _dedupe(positions: list[Position]) -> list[Position]: + positions = sorted(positions) + return [p for p, _ in itertools.groupby(positions)] + + +def map_positions_to_names( + position_dict: Mapping[str, Sequence[Position]], + domain_name_list: Mapping[str, Sequence[Sequence[str]]], +) -> dict[str, dict[str, str]]: + """Re-key each protein's raw InterPro hits from position -> resolved domain name. + + ``domain_name_list[protein_id]`` is a list of ``[database, name]`` pairs in + the same order as ``position_dict[protein_id]`` (both are built hit-by-hit + from the same InterPro response in :mod:`candy.interpro`). + + .. note:: + The output is keyed by ``str(position)``, so when two different + databases report the exact same ``[start, end]`` span (which happens + routinely, e.g. Pfam and SMART agreeing on boundaries), only one name + can survive per position either way. The original notebook paired + these two parallel lists via ``value.index(positions)`` (a lookup *by + value*), which for duplicate positions always resolves to whichever + hit happened to occur first in the raw InterPro response -- an + accident of API ordering, unrelated to the database-preference + ranking applied later. This version pairs by loop index instead (the + two lists are built hit-by-hit in lockstep in ``interpro.py``), so for + duplicate positions the last hit wins, same as any normal dict + overwrite -- still arbitrary between the tied candidates, but no + longer silently locked to index 0 regardless of which hit is current. + """ + result: dict[str, dict[str, str]] = {} + for protein_id, positions in position_dict.items(): + names_for_protein = domain_name_list[protein_id] + result[protein_id] = { + str(position): names_for_protein[index][1] for index, position in enumerate(positions) + } + return result + + +def clean_domains( + position_dict: Mapping[str, Sequence[Position]], + position_to_name: Mapping[str, Mapping[str, str]], + name_to_database: Mapping[str, Mapping[str, str]], + database_preference: Sequence[str], + overlap_fraction: float, + max_domain_length: int, +) -> dict[str, list[NamedDomain]]: + """Resolve each protein's raw domain hits into a clean, non-overlapping architecture.""" + final: dict[str, list[NamedDomain]] = {} + + for protein_id, raw_positions in position_dict.items(): + positions = list(raw_positions) + previous_length = -1 + resolved: list[NamedDomain] = [] + + while previous_length != len(positions): + previous_length = len(positions) + groups = group_overlapping(positions, overlap_fraction) + + resolved = [] + next_positions: list[Position] = [] + + for group in groups: + if _is_single_position(group): + resolved_domain = _resolve_single( + group, protein_id, position_to_name, name_to_database + ) + if resolved_domain is not None: + resolved.append(resolved_domain) + next_positions.append(group) + else: + resolved_domain = _resolve_grouped( + group, + protein_id, + position_to_name, + name_to_database, + database_preference, + max_domain_length, + ) + if resolved_domain is not None: + resolved.append(resolved_domain) + next_positions.append(resolved_domain[0]) + + positions = next_positions + + final[protein_id] = resolved + + return final + + +def _is_single_position(group: Position | list[Position]) -> bool: + return all(isinstance(x, int) for x in group) + + +def _resolve_single( + position: Position, + protein_id: str, + position_to_name: Mapping[str, Mapping[str, str]], + name_to_database: Mapping[str, Mapping[str, str]], +) -> NamedDomain | None: + name = position_to_name[protein_id][str(position)] + database = name_to_database[protein_id][name] + + if ( + name in _EXCLUDED_DOMAIN_NAMES + or database in _EXCLUDED_DATABASES_SINGLE + or "SIGNAL" in database.upper() + ): + return None + return position, name + + +def _resolve_grouped( + group: list[Position], + protein_id: str, + position_to_name: Mapping[str, Mapping[str, str]], + name_to_database: Mapping[str, Mapping[str, str]], + database_preference: Sequence[str], + max_domain_length: int, +) -> NamedDomain | None: + candidates: dict[str, str] = {} # name -> database + candidate_positions: dict[str, Position] = {} # name -> position + + for position in group: + start, end = position + length = end - start + name = position_to_name[protein_id][str(position)] + database = name_to_database[protein_id][name] + + if ( + name in _EXCLUDED_DOMAIN_NAMES + or database in _EXCLUDED_DATABASES_GROUPED + or "SIGNAL" in database.upper() + or "SIGNAL" in name.upper() + or not (10 <= length <= max_domain_length) + ): + continue + + candidates[name] = database + candidate_positions[name] = position + + if not candidates: + return None + + ranked: dict[int, str] = {} + for name, database in candidates.items(): + if database in database_preference: + ranked[database_preference.index(database)] = name + + if ranked: + chosen_name = ranked[min(ranked)] + else: + chosen_name = next(iter(candidate_positions)) + + return candidate_positions[chosen_name], chosen_name diff --git a/src/candy/external_tools.py b/src/candy/external_tools.py new file mode 100644 index 0000000..49f1fd2 --- /dev/null +++ b/src/candy/external_tools.py @@ -0,0 +1,65 @@ +"""Locate external CLI tools (CD-HIT, MMseqs2, MAFFT, FastTree) on PATH. + +The notebook installed these itself at runtime with Colab-only shell magics +(``!apt-get install mafft``, downloading a Linux MMseqs2 binary, compiling +CD-HIT from source). A pip package can't do that portably, so CANDy instead +expects these tools to already be on PATH -- the documented install path is +the bundled conda ``environment.yml`` (bioconda ships all four). +""" + +from __future__ import annotations + +import shutil +import subprocess +from pathlib import Path + +_CONDA_HINT = ( + "Install it via the bundled conda environment: " + "`conda env create -f environment.yml && conda activate candy`." +) + + +class MissingDependencyError(RuntimeError): + """Raised when a required external CLI tool isn't available on PATH.""" + + +def require_binary(name: str, *, hint: str = _CONDA_HINT) -> str: + """Return the resolved path to ``name`` on PATH, or raise a clear error.""" + path = shutil.which(name) + if path is None: + raise MissingDependencyError( + f"Required external tool '{name}' was not found on PATH. {hint}" + ) + return path + + +def find_binary(name: str) -> str | None: + """Return the resolved path to ``name`` on PATH, or None if not found.""" + return shutil.which(name) + + +class ExternalToolError(RuntimeError): + """Raised when an external CLI tool exits with a non-zero status.""" + + +def run_tool( + args: list[str], *, cwd: Path | None = None, stdout_path: Path | None = None +) -> subprocess.CompletedProcess: + """Run an external CLI tool, raising a readable error (with stderr) on failure. + + If ``stdout_path`` is given, the subprocess's stdout is streamed straight + to that file (used for tools like FastTree that write results to stdout). + """ + if stdout_path is not None: + with open(stdout_path, "w") as stdout_file: + result = subprocess.run( + args, cwd=cwd, stdout=stdout_file, stderr=subprocess.PIPE, text=True + ) + else: + result = subprocess.run(args, cwd=cwd, capture_output=True, text=True) + + if result.returncode != 0: + raise ExternalToolError( + f"Command failed ({result.returncode}): {' '.join(args)}\n{result.stderr}" + ) + return result diff --git a/src/candy/fasta.py b/src/candy/fasta.py new file mode 100644 index 0000000..a57472a --- /dev/null +++ b/src/candy/fasta.py @@ -0,0 +1,169 @@ +"""FASTA parsing and header-formatting utilities. + +CANDy encodes taxonomy directly into the FASTA header as +``{accession}_{organism}_{taxonomy_code}`` (e.g. ``P12345_Escherichia_coli_B``) +so that downstream stages (clustering, domain detection, tree annotation) can +recover the organism and kingdom from the sequence ID alone, without carrying +a side table around. ``split_header`` / the ``_`` convention is therefore load +bearing for the rest of the pipeline, exactly as it was in the original +notebook. +""" + +from __future__ import annotations + +import re +from collections.abc import Iterable, Mapping +from dataclasses import dataclass +from pathlib import Path + +from Bio import SeqIO +from Bio.SeqRecord import SeqRecord + +_NON_WORD_RE = re.compile(r"\W+") +_MAX_ORGANISM_LEN = 100 + + +@dataclass(frozen=True) +class FormattedRecord: + """A sequence with its accession, organism and taxonomy code kept separate.""" + + accession: str + organism: str + taxonomy_code: str + sequence: str + + @property + def header(self) -> str: + return f"{self.accession}_{self.organism}_{self.taxonomy_code}" + + def to_fasta(self) -> str: + return f">{self.header}\n{self.sequence}\n" + + +def count_sequences(fasta_path: str | Path) -> int: + """Count the number of records in a FASTA file without fully parsing it.""" + count = 0 + with open(fasta_path) as handle: + for line in handle: + if line.startswith(">"): + count += 1 + return count + + +def parse_fasta_to_dict(fasta_path: str | Path) -> dict[str, str]: + """Return {record_id: sequence} for every record in a FASTA file.""" + with open(fasta_path) as handle: + return {record.id: str(record.seq) for record in SeqIO.parse(handle, "fasta")} + + +def _extract_organism(description: str) -> str | None: + """Pull the organism name out of an NCBI-style '... [Organism name]' description.""" + if "[" not in description: + return None + + organism = description[description.find("[") + 1 : description.rfind("]")] + + # Handle nested brackets, e.g. "[[Candida] auris]" -> "Candida auris" + if "[" in organism and "]" in organism: + organism = organism[organism.find("[") + 1 : organism.find("]")] + organism[organism.find("]") + 1 :] + + organism = organism.replace(" ", "_") + organism = _NON_WORD_RE.sub("", organism) + if len(organism) > _MAX_ORGANISM_LEN: + organism = organism[:_MAX_ORGANISM_LEN] + return organism.rstrip("_") + + +def format_records( + records: Iterable[SeqRecord], taxonomy_dict: Mapping[str, str] +) -> list[FormattedRecord]: + """Deduplicate by sequence/ID, sanitize headers, and encode taxonomy. + + Mirrors the notebook's ``format_fasta``: records without a recognisable + ``[Organism]`` tag in their description are dropped, sequences already + seen (by exact string match) are skipped, and accessions containing ``_`` + have it stripped (since ``_`` is the header field separator). + """ + seen_sequences: set[str] = set() + seen_ids: set[str] = set() + formatted: list[FormattedRecord] = [] + + for record in records: + sequence = str(record.seq) + if sequence in seen_sequences: + continue + if record.id in seen_ids: + continue + + organism = _extract_organism(record.description) + if organism is None: + # No organism annotation to key on; skip, same as the notebook. + continue + + seen_sequences.add(sequence) + seen_ids.add(record.id) + + accession = record.id.replace("_", "") + formatted.append( + FormattedRecord( + accession=accession, + organism=organism, + taxonomy_code=taxonomy_dict.get(record.id, ""), + sequence=sequence, + ) + ) + + return formatted + + +def format_fasta_file( + input_path: str | Path, output_path: str | Path, taxonomy_dict: Mapping[str, str] +) -> int: + """Format a FASTA file on disk and write the result. Returns the record count.""" + with open(input_path) as handle: + records = list(SeqIO.parse(handle, "fasta")) + + formatted = format_records(records, taxonomy_dict) + + with open(output_path, "w") as out: + for record in formatted: + out.write(record.to_fasta()) + + return len(formatted) + + +def write_fasta(records: Iterable[tuple[str, str]], output_path: str | Path) -> None: + """Write an iterable of (header, sequence) pairs as FASTA.""" + with open(output_path, "w") as out: + for header, sequence in records: + out.write(f">{header}\n{sequence}\n") + + +def select_sequences_with_domains( + fasta_path: str | Path, domain_architecture: Mapping[str, list], is_cazy_query: bool +) -> tuple[str, list[str]]: + """Keep only sequences that have at least one detected domain. + + Used before MSA/phylogenetics, since undomained sequences add noise to + the tree without contributing architectural information. + + .. note:: + The notebook recovered each record's bare protein ID via + ``id[:id.find('.')+2]`` -- a formula that only produces the right + accession because NCBI version suffixes are almost always a single + digit (".1", ".2", ...); a double-digit version (".10"+) would be + truncated wrong. The bare ID is already known directly from how + ``domain_architecture`` was keyed, so it's reused here instead of + re-deriving it from the header string. + """ + ids_without_domains = [pid for pid, hits in domain_architecture.items() if not hits] + ids_with_domains = {pid for pid, hits in domain_architecture.items() if hits} + + output = "" + with open(fasta_path) as handle: + for record in SeqIO.parse(handle, "fasta"): + identifier = record.id.split("_")[0] if is_cazy_query else record.id + if identifier in ids_with_domains: + output += f">{record.id}\n{record.seq}\n" + + return output, ids_without_domains diff --git a/src/candy/interpro.py b/src/candy/interpro.py new file mode 100644 index 0000000..128261f --- /dev/null +++ b/src/candy/interpro.py @@ -0,0 +1,230 @@ +"""Client for InterPro's MD5 match-lookup and entry-name APIs. + +The notebook queried the same match-lookup endpoint twice with two separate, +near-duplicate implementations: once to check which sequences have a UniParc +entry at all (:mod:`candy.uniparc`), and again to fetch full domain hits +(``domain_detection``). Both are consolidated here into a single +:func:`match_lookup`, which both call sites now share. +""" + +from __future__ import annotations + +import hashlib +import logging +import time +import urllib.error +import xml.etree.ElementTree as ET +from collections.abc import Mapping, Sequence +from dataclasses import dataclass, field +from urllib.parse import urlencode +from urllib.request import urlopen + +import requests +from tqdm import tqdm + +logger = logging.getLogger(__name__) + +MATCH_LOOKUP_URL = "https://www.ebi.ac.uk/interpro/match-lookup/matches" +ENTRY_API_URL = "https://www.ebi.ac.uk/interpro/api/entry/" +_MD5_BATCH_SIZE = 5000 + +# Maps InterPro's member-database labels to the URL path segment used by the +# entry API. Databases not in this map (TMHMM, MOBIDB_LITE, PHOBIUS, COILS, +# FUNFAM, ...) don't have resolvable "names" via this API; their raw domain +# code is used as the name instead, same as the original notebook. +_DATABASE_TO_ENTRY_PATH = { + "GENE3D": "cathgene3d", + "CDD": "cdd", + "HAMAP": "hamap", + "PANTHER": "panther", + "PFAM": "pfam", + "PIRSF": "pirsf", + "PRINTS": "prints", + "PROSITE_PROFILES": "profile", + "PROSITE_PATTERNS": "prosite", + "SFLD": "sfld", + "SMART": "smart", + "SUPERFAMILY": "ssf", + "TIGRFAMS": "tigrfams", + "NCBIfam": "ncbifam", +} + + +@dataclass(frozen=True) +class DomainHit: + database: str + domain_code: str + start: int + end: int + + +@dataclass +class ProteinMatch: + protein_id: str + hits: list[DomainHit] = field(default_factory=list) + + +@dataclass(frozen=True) +class ResolvedHit: + database: str + name: str + start: int + end: int + + +def _md5_of(sequence: str) -> str: + return hashlib.md5(sequence.encode()).hexdigest().upper() + + +def _batched(items: list, size: int): + for i in range(0, len(items), size): + yield items[i : i + size] + + +def _query_md5_batch(md5_list: Sequence[str], max_retries: int = 10, retry_delay: float = 5.0) -> ET.Element: + data = urlencode({"md5": list(md5_list)}, doseq=True).encode("ascii") + + attempts = 0 + while True: + try: + with urlopen(MATCH_LOOKUP_URL, data) as response: + payload = response.read().decode("utf-8") + return ET.fromstring(payload) + except urllib.error.HTTPError: + attempts += 1 + if attempts >= max_retries: + raise + logger.info("InterPro match-lookup busy, retrying in %.0fs.", retry_delay) + time.sleep(retry_delay) + + +def _parse_hits(hit_text: str) -> DomainHit: + columns = hit_text.split(",") + database = columns[0] + domain_code = columns[2] + # Discontinuous domains report multiple ';'-separated fragments; only the + # first fragment's boundaries are used (matching the published tool). + first_fragment = columns[6].split(";")[0] + start, end, _ = first_fragment.split("-") + return DomainHit(database, domain_code, int(start), int(end)) + + +def match_lookup(proteins: Mapping[str, str]) -> tuple[dict[str, ProteinMatch], list[str]]: + """Query InterPro's MD5 match-lookup for every sequence in ``proteins``. + + ``proteins`` maps sequence ID -> amino-acid sequence. Sequences are + queried by MD5 hash, so IDs sharing an identical sequence share a single + lookup and all receive the same result. Returns ``(matches, unmatched_ids)``. + """ + sequence_to_ids: dict[str, list[str]] = {} + for identifier, sequence in proteins.items(): + sequence_to_ids.setdefault(sequence, []).append(identifier) + + md5_to_sequence = {_md5_of(sequence): sequence for sequence in sequence_to_ids} + matches: dict[str, ProteinMatch] = {} + matched_md5s: set[str] = set() + + for batch in _batched(list(md5_to_sequence.keys()), _MD5_BATCH_SIZE): + root = _query_md5_batch(batch) + for match_el in root.findall(".//match"): + md5 = match_el.find("proteinMD5").text + sequence = md5_to_sequence[md5] + hits = [_parse_hits(hit_el.text) for hit_el in match_el.findall("hit")] + + for identifier in sequence_to_ids[sequence]: + matches[identifier] = ProteinMatch(identifier, hits) + matched_md5s.add(md5) + + unmatched_ids = [ + identifier + for md5, sequence in md5_to_sequence.items() + if md5 not in matched_md5s + for identifier in sequence_to_ids[sequence] + ] + + return matches, unmatched_ids + + +def detect_domains(proteins: Mapping[str, str]) -> tuple[dict[str, ProteinMatch], list[str]]: + """Query InterPro for domain hits on every sequence in ``proteins``.""" + matches, unmatched = match_lookup(proteins) + for protein_id in unmatched: + logger.warning("No domains found for %s; try running InterProScan directly.", protein_id) + return matches, unmatched + + +def _fetch_entry_name( + database: str, domain_code: str, session: requests.Session, max_retries: int = 3 +) -> str | None: + entry_path = _DATABASE_TO_ENTRY_PATH.get(database) + if entry_path is None: + return None + + url = f"{ENTRY_API_URL}{entry_path}/{domain_code}" + attempts = 0 + while True: + try: + response = session.get(url, headers={"Accept": "application/json"}, timeout=30) + except requests.RequestException: + attempts += 1 + if attempts > max_retries: + raise + time.sleep(61) + continue + + if response.status_code == 408: + time.sleep(61) + continue + if response.status_code in (204, 404): + logger.warning("No InterPro entry data for %s/%s.", database, domain_code) + return None + if not response.ok: + attempts += 1 + if attempts > max_retries: + response.raise_for_status() + time.sleep(61) + continue + + return response.json()["metadata"]["name"]["name"] + + +def resolve_domain_names( + matches: Mapping[str, ProteinMatch], session: requests.Session | None = None +) -> dict[str, list[ResolvedHit]]: + """Resolve each hit's InterPro domain code into a human-readable name. + + Names are cached per (database, code) pair for the whole call, since the + same domain accession is typically hit across many proteins. + """ + session = session or requests.Session() + name_cache: dict[tuple[str, str], str] = {} + resolved: dict[str, list[ResolvedHit]] = {} + + for protein_id, match in tqdm(matches.items(), total=len(matches), desc="Resolving domain names"): + resolved_hits = [] + for hit in match.hits: + cache_key = (hit.database, hit.domain_code) + if cache_key not in name_cache: + name = _fetch_entry_name(hit.database, hit.domain_code, session) + name_cache[cache_key] = name if name is not None else hit.domain_code + time.sleep(1) # stay polite to the API, same pacing as the original notebook + resolved_hits.append(ResolvedHit(hit.database, name_cache[cache_key], hit.start, hit.end)) + resolved[protein_id] = resolved_hits + + return resolved + + +def build_domain_tables( + resolved: Mapping[str, Sequence[ResolvedHit]], +) -> tuple[dict[str, list[list[int]]], dict[str, list[list[str]]], dict[str, dict[str, str]]]: + """Convert resolved hits into the (position, name, database) tables :mod:`candy.domains` expects.""" + position_dict: dict[str, list[list[int]]] = {} + domain_name_list: dict[str, list[list[str]]] = {} + name_to_database: dict[str, dict[str, str]] = {} + + for protein_id, hits in resolved.items(): + position_dict[protein_id] = [[hit.start, hit.end] for hit in hits] + domain_name_list[protein_id] = [[hit.database, hit.name] for hit in hits] + name_to_database[protein_id] = {hit.name: hit.database for hit in hits} + + return position_dict, domain_name_list, name_to_database diff --git a/src/candy/itol.py b/src/candy/itol.py new file mode 100644 index 0000000..ca98fa3 --- /dev/null +++ b/src/candy/itol.py @@ -0,0 +1,168 @@ +"""Generate iTOL annotation files: domain architectures and characterized-enzyme styling. + +The notebook used ``ete3`` to load the phylogenetic tree in both annotation +functions. In ``domain_vis`` (here: :func:`build_domain_annotation`) the tree +was loaded into a ``Tree`` object that was then never referenced again -- +dead code, dropped entirely, so this function no longer needs a tree at all. +In ``characterized_labeling`` (here: :func:`build_characterized_style_annotation`) +the tree genuinely is used, to find which leaf names correspond to a given +characterized enzyme's accession; that's done here with ``Bio.Phylo`` +instead, avoiding the largely-unmaintained ``ete3`` dependency entirely. +""" + +from __future__ import annotations + +import logging +from collections.abc import Mapping, Sequence +from pathlib import Path + +from Bio import Phylo + +logger = logging.getLogger(__name__) + +_DOMAIN_SHAPES = ["RE", "EL", "HH", "HV", "DI", "TR", "TL", "PL", "PR", "PU", "PD", "OC"] +_DOMAIN_COLORS = [ + "#BF0F0F", "#217D3B", "#23A5D9", "#F2C641", "#A276DB", "#D9560B", "#F294B6", "#1DF2DD", + "#E80C7A", "#A6A6A6", "#000000", "#352310", "#4F0B79", "#CE07F3", "#12FC4F", "#79D626", + "#4C45FE", "#DDFBE6", "#0842AF", "#83634C", "#82D4FF", +] +_CHARACTERIZED_COLORS = [ + "#BF0F0F", "#217D3B", "#23A5D9", "#F2C641", "#A276DB", "#D9560B", "#F294B6", "#1DF2DD", + "#E80C7A", "#A6A6A6", "#352310", "#4F0B79", "#CE07F3", "#12FC4F", "#79D626", "#4C45FE", + "#DDFBE6", "#0842AF", "#83634C", "#82D4FF", "#FF5733", "#33FF57", "#5733FF", "#FFD700", + "#00BFFF", "#FF69B4", "#8A2BE2", "#7FFF00", "#DC143C", "#20B2AA", "#008080", "#FF4500", + "#DA70D6", "#4682B4", "#D2691E", "#9ACD32", "#00CED1", "#FF1493", "#1E90FF", "#B22222", + "#FF6347", "#00FF7F", "#8B4513", "#2E8B57", "#6A5ACD", "#FFDAB9", "#7B68EE", "#32CD32", + "#FF00FF", "#FFDEAD", +] + + +def _umbrella_name(raw_name: str, curated_domains: Mapping[str, Sequence[str]]) -> str: + for umbrella, raw_names in curated_domains.items(): + if raw_name in raw_names: + return umbrella + return raw_name + + +def _protein_id_for_matching(record_id: str) -> str: + """Recover the bare protein ID from a '{accession}_{organism}_{taxcode}'-style header. + + Ported as-is from the notebook's ``domain_vis``: IDs whose first ``_`` + occurs at index < 2 are used verbatim (covers custom-FASTA-mode headers, + which don't follow the accession/organism/taxcode convention at all), + otherwise the leading accession segment (or, when the accession itself + contains a very short prefix before an underscore, its first two + underscore-separated segments) is used. + """ + first_underscore = record_id.find("_") + if first_underscore > 2: + return record_id.split("_")[0] + if first_underscore < 2: + return record_id + parts = record_id.split("_") + return f"{parts[0]}_{parts[1]}" + + +def build_domain_annotation( + jobname: str, + domain_architecture: Mapping[str, Sequence[tuple[list[int], str]]], + curated_domains: Mapping[str, Sequence[str]], + protein_records: Sequence, +) -> str: + """Build an iTOL DATASET_DOMAINS annotation file mapping each leaf to its domain motifs.""" + shape_dict: dict[str, str] = {} + color_dict: dict[str, str] = {} + data = "" + + for protein_id, domains in domain_architecture.items(): + motif = "" + for position, raw_name in domains: + name = _umbrella_name(raw_name, curated_domains) + if name not in shape_dict: + shape_dict[name] = _DOMAIN_SHAPES[len(shape_dict) % len(_DOMAIN_SHAPES)] + if name not in color_dict: + color_dict[name] = _DOMAIN_COLORS[len(color_dict) % len(_DOMAIN_COLORS)] + start, end = position + motif += f"\t{shape_dict[name]}|{start}|{end}|{color_dict[name]}|{name}" + + for record in protein_records: + if _protein_id_for_matching(record.id) == protein_id: + data += f"{record.id}\t{len(record.seq)}{motif}\n" + + header = ( + "DATASET_DOMAINS\n" + "SEPARATOR TAB\n" + "DATASET_LABEL\tSMART architecture export\n" + "COLOR\t#0000ff\n" + "BORDER_WIDTH\t1\n" + "GRADIENT_FILL\t1\n" + "SHOW_DOMAIN_LABELS\t1\n" + f"LEGEND_TITLE\t {jobname} domains\n" + ) + legend = ( + "LEGEND_LABELS\t\t" + "\t".join(shape_dict.keys()) + "\n" + "LEGEND_SHAPES\t\t" + "\t".join(shape_dict.values()) + "\n" + "LEGEND_COLORS\t\t" + "\t".join(color_dict[name] for name in shape_dict) + "\n" + ) + return header + legend + "DATA\n" + data + + +def build_characterized_style_annotation( + tree_path: str | Path, + characterized_ids: Sequence[str], + ec_numbers: Mapping[str, str], + legend_title: str, +) -> str: + """Build an iTOL DATASET_STYLE file coloring characterized-enzyme leaf labels by EC number. + + .. note:: + The notebook only assigned an ``ecnumber``/color for IDs present in + ``ec_numbers``, but then unconditionally used that (possibly stale, + possibly unset) ``ecnumber`` variable to style *every* ID's matching + leaf -- so a characterized enzyme with no known EC number got styled + with whichever EC color happened to be left over from the previous + one in the loop (or crashed with an ``UnboundLocalError`` if it was + first). Here, only IDs with a resolved EC number are styled. + """ + tree = Phylo.read(str(tree_path), "newick") + leaf_names = [leaf.name for leaf in tree.get_terminals()] + + color_dict: dict[str, str] = {} + data = "" + + for genbank_id in characterized_ids: + if genbank_id not in ec_numbers: + continue + + ec_values = ec_numbers[genbank_id].split(" ") + ec_number = ec_values[0] + if len(ec_values) > 1: + logger.info( + "%s has multiple EC numbers: %s. EC number %s will be used for annotation.", + genbank_id, ec_numbers[genbank_id], ec_number, + ) + if ec_number not in color_dict: + color_dict[ec_number] = _CHARACTERIZED_COLORS[len(color_dict) % len(_CHARACTERIZED_COLORS)] + + for leaf_name in leaf_names: + if genbank_id in leaf_name: + data += f"{leaf_name} label node {color_dict[ec_number]} 1 bold\n" + + header = ( + "DATASET_STYLE\n" + "SEPARATOR SPACE\n" + "DATASET_LABEL Characterized enzymes\n" + "COLOR #0000ff\n" + "BORDER_WIDTH 1\n" + "GRADIENT_FILL 1\n" + "SHOW_DOMAIN_LABELS 1\n" + f"LEGEND_TITLE {legend_title} characterized enzymes\n" + ) + # All legend swatches intentionally share shape code "1" (a single + # uniform marker), same as the published notebook. + legend = ( + "LEGEND_LABELS " + " ".join(color_dict.keys()) + "\n" + "LEGEND_SHAPES " + " ".join("1" for _ in color_dict) + "\n" + "LEGEND_COLORS " + " ".join(color_dict.values()) + "\n" + ) + return header + legend + "DATA\n" + data diff --git a/src/candy/merge.py b/src/candy/merge.py new file mode 100644 index 0000000..2c27acc --- /dev/null +++ b/src/candy/merge.py @@ -0,0 +1,81 @@ +"""Merge clustered sequences with characterized enzymes into one FASTA file. + +Characterized (experimentally studied) enzymes must always be present in +the final analysis, even if clustering discarded them as redundant. This +drops any clustered copy of a characterized sequence and re-adds the +characterized version instead (using its BLAST-homolog stand-in sequence +when the original had no UniParc entry, see :mod:`candy.blast`), so each +characterized enzyme appears exactly once, correctly annotated. +""" + +from __future__ import annotations + +from collections.abc import Mapping, Sequence +from pathlib import Path + +from Bio import SeqIO + + +def _protein_id_for_merge(record_id: str) -> str: + """Recover the bare protein ID from a FASTA record id for merge-matching purposes. + + Ported as-is from the notebook's characterized-enzyme merge cell, which + used a subtly different rule than the one in ``itol.py``: any ``_`` found + at position 0-2 joins the first two underscore-separated segments, and + only a literal absence of ``_`` falls back to the id verbatim. + """ + first_underscore = record_id.find("_") + if first_underscore > 2: + return record_id.split("_")[0] + if first_underscore != -1: + parts = record_id.split("_") + return f"{parts[0]}_{parts[1]}" + return record_id + + +def merge_characterized_sequences( + clustered_fasta: str | Path, + characterized_fasta: str | Path, + output_fasta: str | Path, + verified_ids: Sequence[str], + blast_choice: Mapping[str, str], + blast_hit_sequences: Mapping[str, str], +) -> Path: + """Write ``output_fasta`` = clustered sequences (minus characterized duplicates) + characterized sequences.""" + verified_id_set = set(verified_ids) + + with open(output_fasta, "w") as out: + with open(clustered_fasta) as handle: + for record in SeqIO.parse(handle, "fasta"): + if _protein_id_for_merge(record.id) not in verified_id_set: + SeqIO.write(record, out, "fasta") + + characterized_output = "" + with open(characterized_fasta) as handle: + for record in SeqIO.parse(handle, "fasta"): + protein_id = _protein_id_for_merge(record.id) + if protein_id not in verified_id_set: + continue + + accession = protein_id.replace("_", "") + if record.id.find("_") > 2: + description = "_".join(record.id.split("_", 1)[1:]) + else: + description = "_".join(record.id.split("_", 2)[2:]) + + if protein_id in blast_choice: + sequence = blast_hit_sequences[blast_choice[protein_id]] + else: + sequence = str(record.seq) + + # `description` (everything after the accession) already ends + # in "_{taxcode}", since candy.cazy.fetch_characterized_sequences + # writes characterized-enzyme headers as accession_organism_taxcode. + # The original notebook re-appended `taxonomydict[proteinidentifier]` + # here regardless, duplicating the code (e.g. "..._coli_B_B") and + # corrupting the organism name derived from it downstream. + characterized_output += f">{accession}_{description}\n{sequence}\n" + + out.write(characterized_output) + + return Path(output_fasta) diff --git a/src/candy/network.py b/src/candy/network.py new file mode 100644 index 0000000..324d718 --- /dev/null +++ b/src/candy/network.py @@ -0,0 +1,80 @@ +"""Protein domain co-occurrence network. + +Builds a graph where nodes are domain names and edges connect domains that +appear adjacently within the same protein's architecture, so recurring +architectural patterns across a whole CAZy family become visible at a +glance. Exported as GraphML for Cytoscape, same as the original notebook. +""" + +from __future__ import annotations + +from collections import Counter +from collections.abc import Mapping +from pathlib import Path + +import networkx as nx + +ARCHITECTURE_SEPARATOR = "--" + + +def build_cooccurrence_network(domain_architectures: Mapping[str, str]) -> nx.Graph: + """Build a domain co-occurrence graph from {protein_id: '--'-joined domain names}. + + Each node's ``size`` attribute is the number of times that domain occurs + across all architectures; each edge's ``width`` is the minimum occurrence + count of the two domains it connects. + + .. note:: + The original notebook only added a node via ``add_edge``, so a domain + that only ever occurred alone (no co-occurring partner in any + architecture) never appeared in the network at all. Here every + observed domain gets a node, including isolated ones with no edges, + since a domain's absence of co-occurrence is itself informative. + """ + graph = nx.Graph() + domain_counts: Counter[str] = Counter() + + for architecture in domain_architectures.values(): + domains = [d for d in architecture.split(ARCHITECTURE_SEPARATOR) if d] + domain_counts.update(domains) + graph.add_nodes_from(domains) + for a, b in zip(domains, domains[1:]): + graph.add_edge(a, b) + + for node in graph.nodes(): + graph.nodes[node]["size"] = domain_counts[node] + + for source, target in graph.edges(): + graph.edges[source, target]["width"] = min(domain_counts[source], domain_counts[target]) + + return graph + + +def write_graphml(graph: nx.Graph, output_path: str | Path) -> None: + nx.write_graphml(graph, output_path) + + +def plot_network(graph: nx.Graph, output_path: str | Path | None = None, *, title: str = "Domain Co-Occurrence Network"): + """Render a simple spring-layout visualisation; returns the Matplotlib figure. + + Saves to ``output_path`` if given, otherwise leaves the figure for the + caller to show/save (kept out of interactive display so this works + headlessly in scripts and CI). + """ + import matplotlib.pyplot as plt + + node_sizes = [graph.nodes[n].get("size", 1) * 100 for n in graph.nodes()] + edge_widths = [graph.edges[e].get("width", 1) for e in graph.edges()] + + fig, ax = plt.subplots() + pos = nx.spring_layout(graph, seed=42) + nx.draw_networkx_nodes(graph, pos, node_size=node_sizes, node_color="skyblue", ax=ax) + nx.draw_networkx_edges(graph, pos, width=edge_widths, edge_color="gray", alpha=0.7, ax=ax) + nx.draw_networkx_labels(graph, pos, font_size=10, font_color="black", ax=ax) + ax.set_title(title) + ax.axis("off") + + if output_path is not None: + fig.savefig(output_path, bbox_inches="tight") + + return fig diff --git a/src/candy/phylogenetics/__init__.py b/src/candy/phylogenetics/__init__.py new file mode 100644 index 0000000..2d689a9 --- /dev/null +++ b/src/candy/phylogenetics/__init__.py @@ -0,0 +1,22 @@ +"""Pluggable phylogenetic tree-building backends.""" + +from __future__ import annotations + +from pathlib import Path +from typing import Protocol + + +class TreeBuilder(Protocol): + name: str + + def build_tree(self, alignment_fasta: Path, output_newick: Path) -> Path: + """Build a tree from an aligned FASTA file, writing Newick to ``output_newick``.""" + ... + + +def get_tree_builder(name: str) -> TreeBuilder: + if name == "fasttree": + from candy.phylogenetics.fasttree import FastTreeBuilder + + return FastTreeBuilder() + raise ValueError(f"Unknown tree builder: {name}") diff --git a/src/candy/phylogenetics/fasttree.py b/src/candy/phylogenetics/fasttree.py new file mode 100644 index 0000000..c3642be --- /dev/null +++ b/src/candy/phylogenetics/fasttree.py @@ -0,0 +1,33 @@ +from __future__ import annotations + +import logging +from pathlib import Path + +from candy.external_tools import MissingDependencyError, find_binary, run_tool + +logger = logging.getLogger(__name__) + +# Different packagers ship this binary under different casings +# (bioconda historically used 'FastTree', some distros use 'fasttree'). +_CANDIDATE_NAMES = ["FastTree", "fasttree"] + + +class FastTreeBuilder: + name = "fasttree" + + def build_tree(self, alignment_fasta: Path, output_newick: Path) -> Path: + binary = self._locate_binary() + logger.info("Building phylogenetic tree with FastTree.") + run_tool([binary, str(alignment_fasta)], stdout_path=output_newick) + return output_newick + + def _locate_binary(self) -> str: + for candidate in _CANDIDATE_NAMES: + path = find_binary(candidate) + if path is not None: + return path + raise MissingDependencyError( + f"Required external tool 'FastTree' was not found on PATH " + f"(tried: {', '.join(_CANDIDATE_NAMES)}). Install it via the bundled conda " + "environment: `conda env create -f environment.yml && conda activate candy`." + ) diff --git a/src/candy/pipeline.py b/src/candy/pipeline.py new file mode 100644 index 0000000..c431c7f --- /dev/null +++ b/src/candy/pipeline.py @@ -0,0 +1,250 @@ +"""End-to-end orchestration of the CANDy pipeline. + +Replaces the notebook's implicit, cell-execution-order state (a `go` boolean +gating which cells ran, module-level globals like `family`/`blastchar` +threaded silently between cells) with one explicit function per input mode, +sharing the same downstream stages (domain detection onward). + +Digging through the notebook to write this orchestrator surfaced that the +two input modes are *far* less symmetric than they first appear: for a +custom FASTA file, CAZy extraction, UniParc filtering, formatting, +clustering, and characterized-enzyme merging are all skipped entirely -- +the uploaded file goes straight into InterPro domain detection. That's +preserved here as two genuinely different code paths that reconverge at +domain detection, rather than threading a boolean through one path. +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass +from pathlib import Path + +from Bio import SeqIO + +from candy import cazy, database, domains, fasta, interpro, itol, merge, network, uniparc +from candy.alignment import get_alignment_tool +from candy.blast import resolve_characterized_via_blast +from candy.clustering import get_clusterer +from candy.config import CAZyFamilyInput, CustomFastaInput, PipelineConfig, Taxonomy +from candy.curation import get_curation_backend +from candy.phylogenetics import get_tree_builder + +logger = logging.getLogger(__name__) + + +@dataclass +class PipelineResult: + jobname_dir: Path + database_path: Path + network_graphml_path: Path + domain_annotation_path: Path | None + characterized_annotation_path: Path | None + alignment_path: Path | None + tree_path: Path | None + sequence_count: int + + +def _accession_with_version(record_id: str) -> str: + """Recover 'ACCESSION.VERSION' from a header, e.g. 'P12345.1' from 'P12345.1_Org_B'.""" + dot = record_id.find(".") + return record_id[: dot + 2] if dot != -1 else record_id + + +def _prepare_domain_detection_input( + fasta_path: Path, + is_cazy_query: bool, + blast_choice: dict[str, str], + blast_hit_sequences: dict[str, str], +) -> dict[str, str]: + proteins: dict[str, str] = {} + with open(fasta_path) as handle: + for record in SeqIO.parse(handle, "fasta"): + seq_id = record.id.split("_")[0] if is_cazy_query else record.id + proteins[seq_id] = ( + blast_hit_sequences[blast_choice[seq_id]] if seq_id in blast_choice else str(record.seq) + ) + return proteins + + +def run_pipeline(config: PipelineConfig) -> PipelineResult: + jobname_dir = config.output_dir / config.jobname + jobname_dir.mkdir(parents=True, exist_ok=True) + + if isinstance(config.input, CAZyFamilyInput): + sequences_fasta, taxonomy_dict, characterized_ids, blast_choice, blast_hit_sequences, ec_numbers = ( + _prepare_cazy_query(config, jobname_dir) + ) + is_cazy_query = True + elif isinstance(config.input, CustomFastaInput): + sequences_fasta = config.input.fasta_path + taxonomy_dict, characterized_ids, blast_choice, blast_hit_sequences, ec_numbers = {}, set(), {}, {}, {} + is_cazy_query = False + else: + raise TypeError(f"Unsupported input type: {type(config.input)}") + + # --- Stage 6: InterPro domain detection --- + proteins = _prepare_domain_detection_input(sequences_fasta, is_cazy_query, blast_choice, blast_hit_sequences) + matches, _unmatched = interpro.detect_domains(proteins) + resolved = interpro.resolve_domain_names(matches) + position_dict, domain_name_list, name_to_database = interpro.build_domain_tables(resolved) + + # --- Stage 7: resolve overlapping domain hits into a clean architecture --- + position_to_name = domains.map_positions_to_names(position_dict, domain_name_list) + domain_architecture = domains.clean_domains( + position_dict, + position_to_name, + name_to_database, + config.domain_cleaning.database_preference, + config.domain_cleaning.overlap_percentage / 100, + config.domain_cleaning.max_domain_length, + ) + + # --- Stage 8: curate synonymous domain names --- + raw_domain_names = sorted({name for hits in domain_architecture.values() for _, name in hits}) + curation_kwargs = {} if config.curation.backend == "manual" else {"api_key": config.curation.api_key} + curation_backend = get_curation_backend(config.curation.backend, **curation_kwargs) + family_label = config.input.family if is_cazy_query else None + curated_domains = curation_backend.curate(raw_domain_names, family=family_label) + + # --- Stage 9: build the result database --- + with open(sequences_fasta) as handle: + protein_records = list(SeqIO.parse(handle, "fasta")) + + db_path = jobname_dir / f"{config.jobname}_db.db" + engine = database.create_database(db_path) + database.populate_database( + engine, + protein_records, + characterized_ids=characterized_ids, + domain_architecture=domain_architecture, + domain_to_database=name_to_database, + curated_domains=curated_domains, + is_cazy_query=is_cazy_query, + ) + + # --- Stage 10: domain co-occurrence network --- + db_rows = database.read_protein_sequences(engine) + architectures = {row.protein_sequence_id: row.domain_architecture for row in db_rows} + graph = network.build_cooccurrence_network(architectures) + graphml_path = jobname_dir / f"{config.jobname}_domain_cooccurence_network.graphml" + network.write_graphml(graph, graphml_path) + + alignment_path: Path | None = None + tree_path: Path | None = None + domain_annotation_path: Path | None = None + characterized_annotation_path: Path | None = None + + if config.build_tree: + # --- Stage 10.5: drop sequences with no detected domain before MSA/PTI --- + selected_fasta_text, ids_without_domains = fasta.select_sequences_with_domains( + sequences_fasta, domain_architecture, is_cazy_query + ) + if ids_without_domains: + logger.info("Excluding %d sequences with no detected domain from MSA/PTI.", len(ids_without_domains)) + selected_fasta = jobname_dir / f"CAZy_{config.jobname}_inclchar_selected.fasta" + selected_fasta.write_text(selected_fasta_text) + + # --- Stage 11: MSA --- + alignment_tool = get_alignment_tool(config.alignment_tool) + alignment_path = jobname_dir / f"CAZy_{config.jobname}_aligned.fasta" + alignment_tool.align(selected_fasta, alignment_path) + + # --- Stage 12: phylogenetic tree --- + tree_builder = get_tree_builder(config.tree_tool) + tree_path = jobname_dir / f"CAZy_{config.jobname}_phyltree.nwk" + tree_builder.build_tree(alignment_path, tree_path) + + # --- Stage 13: iTOL annotation files --- + with open(selected_fasta) as handle: + selected_records = list(SeqIO.parse(handle, "fasta")) + + domain_annotation_path = jobname_dir / f"iTOL_annotation_CAZy_{config.jobname}.txt" + domain_annotation_path.write_text( + itol.build_domain_annotation(config.jobname, domain_architecture, curated_domains, selected_records) + ) + + if is_cazy_query: + characterized_annotation_path = jobname_dir / f"iTOL_annotation_CAZy_{config.jobname}_characterized.txt" + characterized_annotation_path.write_text( + itol.build_characterized_style_annotation( + tree_path, sorted(characterized_ids), ec_numbers, family_label + ) + ) + + return PipelineResult( + jobname_dir=jobname_dir, + database_path=db_path, + network_graphml_path=graphml_path, + domain_annotation_path=domain_annotation_path, + characterized_annotation_path=characterized_annotation_path, + alignment_path=alignment_path, + tree_path=tree_path, + sequence_count=len(protein_records), + ) + + +def _prepare_cazy_query( + config: PipelineConfig, jobname_dir: Path +) -> tuple[Path, dict[str, str], set[str], dict[str, str], dict[str, str], dict[str, str]]: + cazy_input: CAZyFamilyInput = config.input + family = cazy_input.family + taxonomy = cazy_input.taxonomy + + # --- Stage 1: extract CAZy family sequences --- + fasta_text, taxonomy_dict = cazy.extract_family_sequences(family, taxonomy, cazy_input.email) + raw_fasta = jobname_dir / f"CAZy_{family}_{taxonomy.value}_FASTA_sequences.fasta" + raw_fasta.write_text(fasta_text) + + # --- Stage 2: exclude sequences without a UniParc entry --- + proteins = fasta.parse_fasta_to_dict(raw_fasta) + available_ids, _unavailable_ids = uniparc.filter_available_in_uniparc(proteins) + verified_fasta = jobname_dir / f"CAZy_{family}_{taxonomy.value}_Verified_FASTA_sequences.fasta" + uniparc.write_verified_fasta(raw_fasta, verified_fasta, set(available_ids)) + + # --- Stage 3: format FASTA (dedupe, sanitize headers, encode taxonomy) --- + formatted_fasta = jobname_dir / f"CAZy_{family}_{taxonomy.value}_FASTA_sequences_formatted.fasta" + fasta.format_fasta_file(verified_fasta, formatted_fasta, taxonomy_dict) + + # --- Stage 4: clustering --- + clusterer = get_clusterer(config.clustering.software) + cutoff_label = f"{config.clustering.identity_cutoff}pct" + clustered_fasta = jobname_dir / f"CAZy_{family}_{taxonomy.value}_formatted_{cutoff_label}.fasta" + clusterer.cluster(formatted_fasta, clustered_fasta, config.clustering) + + # --- Stage 5: include characterized sequences --- + characterized_tables = cazy.fetch_characterized_page(family) + char_ids, char_taxonomy_dict, ec_numbers = cazy.extract_characterized_ids(characterized_tables) + characterized_fasta_text = cazy.fetch_characterized_sequences( + char_ids, char_taxonomy_dict, taxonomy, cazy_input.email + ) + characterized_fasta = jobname_dir / f"Characterized_{family}_{taxonomy.value}_FASTA_sequences.fasta" + characterized_fasta.write_text(characterized_fasta_text) + + available_id_set = set(available_ids) + verified_ids: list[str] = [] + ids_without_uniparc: list[str] = [] + with open(characterized_fasta) as handle: + for record in SeqIO.parse(handle, "fasta"): + accession = _accession_with_version(record.id) + if accession in available_id_set: + verified_ids.append(accession) + else: + ids_without_uniparc.append(accession) + + blast_choice, blast_hit_sequences = resolve_characterized_via_blast( + characterized_fasta, ids_without_uniparc, config.blast_identity_threshold + ) + verified_ids.extend(blast_choice.keys()) + + merged_fasta = jobname_dir / f"CAZy_{family}_{taxonomy.value}_{config.clustering.identity_cutoff}pct_inclchar.fasta" + merge.merge_characterized_sequences( + clustered_fasta, + characterized_fasta, + merged_fasta, + verified_ids, + blast_choice, + blast_hit_sequences, + ) + + return merged_fasta, taxonomy_dict, set(verified_ids), blast_choice, blast_hit_sequences, ec_numbers diff --git a/src/candy/uniparc.py b/src/candy/uniparc.py new file mode 100644 index 0000000..a7e13f0 --- /dev/null +++ b/src/candy/uniparc.py @@ -0,0 +1,58 @@ +"""Filter sequences down to those with a UniParc entry. + +InterPro's domain databases key off UniParc, so a sequence without a UniParc +entry can never get a domain annotation. Filtering these out early avoids +wasting the (slower, rate-limited) domain-detection step on sequences that +are guaranteed to come back empty. +""" + +from __future__ import annotations + +import logging +from collections.abc import Mapping +from pathlib import Path + +from Bio import SeqIO + +from candy.interpro import match_lookup + +logger = logging.getLogger(__name__) + + +def filter_available_in_uniparc(proteins: Mapping[str, str]) -> tuple[list[str], list[str]]: + """Return (ids with a UniParc entry, ids without one) for ``{id: sequence}``.""" + matches, unmatched = match_lookup(proteins) + available_ids = list(matches.keys()) + + logger.info( + "%d of %d sequences have a UniParc entry and will be included.", + len(available_ids), + len(available_ids) + len(unmatched), + ) + if unmatched: + logger.info("%d sequences have no UniParc entry and will be excluded.", len(unmatched)) + + return available_ids, unmatched + + +def write_verified_fasta(input_path: str | Path, output_path: str | Path, available_ids: set[str]) -> int: + """Write only records whose id is in ``available_ids`` and aren't marked 'partial'. + + Returns the number of sequences written. + """ + written = 0 + partial_count = 0 + + with open(input_path) as handle, open(output_path, "w") as out: + for record in SeqIO.parse(handle, "fasta"): + if "partial" in record.description: + partial_count += 1 + continue + if record.id in available_ids: + out.write(f">{record.description}\n{record.seq}\n") + written += 1 + + if partial_count: + logger.info("%d partial sequences excluded from the analysis.", partial_count) + + return written diff --git a/tests/test_alignment_and_phylogenetics.py b/tests/test_alignment_and_phylogenetics.py new file mode 100644 index 0000000..dc192f9 --- /dev/null +++ b/tests/test_alignment_and_phylogenetics.py @@ -0,0 +1,35 @@ +from unittest.mock import patch + +from candy.alignment.mafft import MafftAligner +from candy.phylogenetics.fasttree import FastTreeBuilder + + +def test_mafft_aligner_invokes_binary_with_stdout_redirect(tmp_path): + input_fasta = tmp_path / "in.fasta" + output_fasta = tmp_path / "out.fasta" + + with patch("candy.alignment.mafft.require_binary", return_value="/usr/bin/mafft"), patch( + "candy.alignment.mafft.run_tool" + ) as mock_run: + MafftAligner().align(input_fasta, output_fasta) + + args, kwargs = mock_run.call_args + assert args[0] == ["/usr/bin/mafft", str(input_fasta)] + assert kwargs["stdout_path"] == output_fasta + + +def test_fasttree_builder_falls_back_to_lowercase_binary_name(tmp_path): + alignment = tmp_path / "aligned.fasta" + output = tmp_path / "tree.nwk" + + def fake_find_binary(name): + return "/usr/bin/fasttree" if name == "fasttree" else None + + with patch("candy.phylogenetics.fasttree.find_binary", side_effect=fake_find_binary), patch( + "candy.phylogenetics.fasttree.run_tool" + ) as mock_run: + FastTreeBuilder().build_tree(alignment, output) + + args, kwargs = mock_run.call_args + assert args[0] == ["/usr/bin/fasttree", str(alignment)] + assert kwargs["stdout_path"] == output diff --git a/tests/test_blast.py b/tests/test_blast.py new file mode 100644 index 0000000..34f1070 --- /dev/null +++ b/tests/test_blast.py @@ -0,0 +1,84 @@ +from types import SimpleNamespace +from unittest.mock import patch + +from candy.blast import resolve_characterized_via_blast + + +def make_blast_record(query_length, alignments): + return SimpleNamespace(query_length=query_length, alignments=alignments) + + +def make_alignment(accession, hsps): + return SimpleNamespace(accession=accession, hsps=hsps) + + +def make_hsp(identities, sbjct): + return SimpleNamespace(identities=identities, sbjct=sbjct) + + +def write_fasta(path, records): + with open(path, "w") as f: + for header, seq in records: + f.write(f">{header}\n{seq}\n") + + +def test_resolve_characterized_via_blast_picks_best_uniparc_hit(tmp_path): + fasta_path = tmp_path / "characterized.fasta" + write_fasta(fasta_path, [("P1_Escherichia_coli_B", "MKVLA")]) + + blast_records = [ + make_blast_record( + query_length=100, + alignments=[ + make_alignment("ACC1", [make_hsp(identities=96, sbjct="MKV")]), # 96% + make_alignment("ACC2", [make_hsp(identities=98, sbjct="MKL")]), # 98%, higher + ], + ) + ] + + with patch("candy.blast.NCBIWWW.qblast", return_value=object()), patch( + "candy.blast.NCBIXML.parse", return_value=iter(blast_records) + ), patch("candy.blast.match_lookup", return_value=({"ACC1": object(), "ACC2": object()}, [])): + chosen, sequences = resolve_characterized_via_blast(fasta_path, ["P1"], identity_threshold=95) + + assert chosen == {"P1": "ACC2"} + assert sequences == {"ACC2": "MKL"} + + +def test_resolve_characterized_via_blast_skips_id_with_no_hits_above_threshold(tmp_path): + fasta_path = tmp_path / "characterized.fasta" + write_fasta(fasta_path, [("P1_Org_B", "MKVLA")]) + + blast_records = [make_blast_record(query_length=100, alignments=[])] + + with patch("candy.blast.NCBIWWW.qblast", return_value=object()), patch( + "candy.blast.NCBIXML.parse", return_value=iter(blast_records) + ): + chosen, sequences = resolve_characterized_via_blast(fasta_path, ["P1"], identity_threshold=95) + + assert chosen == {} + assert sequences == {} + + +def test_resolve_characterized_via_blast_does_not_leak_state_between_ids(tmp_path): + # Regression test for the original notebook's bug: a later sequence with + # zero hits of its own must not pick up an earlier sequence's UniParc + # match and crash (or silently succeed) via shared accumulated state. + fasta_path = tmp_path / "characterized.fasta" + write_fasta(fasta_path, [("P1_Org_B", "MKVLA"), ("P2_Org_B", "AAAAA")]) + + call_count = {"n": 0} + + def fake_parse(handle): + call_count["n"] += 1 + if call_count["n"] == 1: + return iter([make_blast_record(100, [make_alignment("ACC1", [make_hsp(99, "MKV")])])]) + return iter([make_blast_record(100, [])]) # second id: no hits at all + + with patch("candy.blast.NCBIWWW.qblast", return_value=object()), patch( + "candy.blast.NCBIXML.parse", side_effect=fake_parse + ), patch("candy.blast.match_lookup", return_value=({"ACC1": object()}, [])): + chosen, sequences = resolve_characterized_via_blast(fasta_path, ["P1", "P2"], identity_threshold=95) + + assert chosen == {"P1": "ACC1"} + assert "P2" not in chosen diff --git a/tests/test_cazy.py b/tests/test_cazy.py new file mode 100644 index 0000000..7e8761e --- /dev/null +++ b/tests/test_cazy.py @@ -0,0 +1,42 @@ +import pandas as pd + +from candy.cazy import extract_characterized_ids, extract_ncbi_ids, parse_family_table +from candy.config import Taxonomy + +SAMPLE_FAMILY_TABLE = ( + "GH5\tBacteria\tEcoli\tP12345\tncbi\n" + "GH5\tArchaea\tSacc\tP67890\tjgi\n" + "GH5\tEukaryota\tYeast\tP11111\tncbi\n" +) + + +def test_parse_family_table_builds_taxonomy_dict(): + df, taxonomy_dict = parse_family_table(SAMPLE_FAMILY_TABLE) + assert taxonomy_dict == {"P12345": "B", "P67890": "A", "P11111": "E"} + assert len(df) == 3 + + +def test_extract_ncbi_ids_excludes_jgi_sourced(): + df, taxonomy_dict = parse_family_table(SAMPLE_FAMILY_TABLE) + ids = extract_ncbi_ids(df, taxonomy_dict, Taxonomy.ALL) + assert ids == ["P12345", "P11111"] + + +def test_extract_ncbi_ids_filters_by_taxonomy(): + df, taxonomy_dict = parse_family_table(SAMPLE_FAMILY_TABLE) + ids = extract_ncbi_ids(df, taxonomy_dict, Taxonomy.BACTERIA) + assert ids == ["P12345"] + + +def test_extract_characterized_ids_groups_by_taxonomy_header(): + # column 4 mimics CAZy's characterized table layout: taxonomy header rows + # interspersed with GenBank-ID rows ("ID.version"). + col4 = pd.Series(["Bacteria", "P12345.1", "P12345.1", "Archaea", "Q99999.2"]) + col1 = pd.Series(["", "1.2.4.-", "1.2.4.-", "", "3.2.1.-"]) + table = pd.DataFrame({1: col1, 4: col4}) + + ids, taxonomy_dict, activity_dict = extract_characterized_ids([None, table]) + + assert ids == ["P12345.1", "Q99999.2"] + assert taxonomy_dict == {"P12345.1": "B", "Q99999.2": "A"} + assert activity_dict == {"P12345.1": "1.2.4.-", "Q99999.2": "3.2.1.-"} diff --git a/tests/test_cli.py b/tests/test_cli.py new file mode 100644 index 0000000..6fece99 --- /dev/null +++ b/tests/test_cli.py @@ -0,0 +1,55 @@ +from typer.testing import CliRunner + +from candy.cli import _parse_family, app + +runner = CliRunner() + + +def test_parse_family_no_subfamily(): + assert _parse_family("GH5") == ("GH", 5, None) + + +def test_parse_family_with_subfamily(): + assert _parse_family("GH5_1") == ("GH", 5, "1") + + +def test_parse_family_invalid_raises(): + import typer + + try: + _parse_family("not-a-family") + assert False, "expected BadParameter" + except typer.BadParameter: + pass + + +def test_cli_help(): + result = runner.invoke(app, ["--help"]) + assert result.exit_code == 0 + assert "run" in result.output.lower() + + +def test_cli_run_help(): + result = runner.invoke(app, ["run", "--help"]) + assert result.exit_code == 0 + assert "--family" in result.output + assert "--fasta" in result.output + + +def test_cli_run_requires_family_or_fasta(): + result = runner.invoke(app, ["run", "--jobname", "test"]) + assert result.exit_code != 0 + + +def test_cli_run_rejects_both_family_and_fasta(tmp_path): + fasta_path = tmp_path / "in.fasta" + fasta_path.write_text(">a\nMKV\n") + result = runner.invoke( + app, ["run", "--jobname", "test", "--family", "GH5", "--fasta", str(fasta_path), "--email", "a@b.com"] + ) + assert result.exit_code != 0 + + +def test_cli_run_family_without_email_fails(): + result = runner.invoke(app, ["run", "--jobname", "test", "--family", "GH5"]) + assert result.exit_code != 0 diff --git a/tests/test_clustering.py b/tests/test_clustering.py new file mode 100644 index 0000000..5a4577a --- /dev/null +++ b/tests/test_clustering.py @@ -0,0 +1,55 @@ +from unittest.mock import patch + +from candy.config import ClusteringConfig, ClusteringSoftware +from candy.clustering import get_clusterer + + +def test_get_clusterer_returns_expected_implementation(): + assert get_clusterer(ClusteringSoftware.CDHIT).name == "cd-hit" + assert get_clusterer(ClusteringSoftware.MMSEQS2).name == "mmseqs2" + + +def test_cdhit_builds_expected_command(tmp_path): + from candy.clustering.cdhit import CdHitClusterer + + config = ClusteringConfig(identity_cutoff=85, cdhit_min_short_coverage=90, cdhit_min_long_coverage=90) + input_fasta = tmp_path / "in.fasta" + output_fasta = tmp_path / "out.fasta" + + with patch("candy.clustering.cdhit.require_binary", return_value="/usr/bin/cd-hit"), patch( + "candy.clustering.cdhit.run_tool" + ) as mock_run: + CdHitClusterer().cluster(input_fasta, output_fasta, config) + + args = mock_run.call_args[0][0] + assert args[0] == "/usr/bin/cd-hit" + assert "-c" in args and args[args.index("-c") + 1] == "0.85" + assert "-aL" in args and args[args.index("-aL") + 1] == "0.9" + assert "-aS" in args and args[args.index("-aS") + 1] == "0.9" + assert "-i" in args and args[args.index("-i") + 1] == str(input_fasta) + assert "-o" in args and args[args.index("-o") + 1] == str(output_fasta) + + +def test_mmseqs2_builds_expected_command_and_copies_result(tmp_path): + from candy.clustering.mmseqs2 import Mmseqs2Clusterer + + config = ClusteringConfig(identity_cutoff=85, mmseqs_min_coverage=90, mmseqs_cov_mode=0) + input_fasta = tmp_path / "in.fasta" + output_fasta = tmp_path / "out.fasta" + + def fake_run_tool(args, **kwargs): + # mimic mmseqs2 easy-cluster writing its representative-sequence output + prefix = args[3] + with open(f"{prefix}_rep_seq.fasta", "w") as f: + f.write(">rep\nMKV\n") + + with patch("candy.clustering.mmseqs2.require_binary", return_value="/usr/bin/mmseqs"), patch( + "candy.clustering.mmseqs2.run_tool", side_effect=fake_run_tool + ) as mock_run: + Mmseqs2Clusterer().cluster(input_fasta, output_fasta, config) + + args = mock_run.call_args[0][0] + assert args[0] == "/usr/bin/mmseqs" + assert args[1] == "easy-cluster" + assert "--min-seq-id" in args and args[args.index("--min-seq-id") + 1] == "0.85" + assert output_fasta.read_text() == ">rep\nMKV\n" diff --git a/tests/test_curation.py b/tests/test_curation.py new file mode 100644 index 0000000..533e866 --- /dev/null +++ b/tests/test_curation.py @@ -0,0 +1,49 @@ +from candy.curation import get_curation_backend +from candy.curation.manual import ManualCurationBackend + + +def test_manual_curation_groups_by_index(): + domains = ["Catalytic domain A", "Catalytic domain B", "CBM"] + responses = iter(["Catalytic domain", "0,1", "STOP"]) + printed = [] + + backend = ManualCurationBackend(input_fn=lambda _: next(responses), print_fn=printed.append) + result = backend.curate(domains) + + assert result == { + "Catalytic domain": ["Catalytic domain A", "Catalytic domain B"], + "CBM": ["CBM"], + } + + +def test_manual_curation_stop_immediately_keeps_all_names_as_is(): + domains = ["A", "B"] + responses = iter(["stop"]) + + backend = ManualCurationBackend(input_fn=lambda _: next(responses), print_fn=lambda _: None) + result = backend.curate(domains) + + assert result == {"A": ["A"], "B": ["B"]} + + +def test_manual_curation_ends_when_all_domains_grouped(): + domains = ["A", "B"] + responses = iter(["Grouped", "0,1"]) + + backend = ManualCurationBackend(input_fn=lambda _: next(responses), print_fn=lambda _: None) + result = backend.curate(domains) + + assert result == {"Grouped": ["A", "B"]} + + +def test_get_curation_backend_returns_manual(): + backend = get_curation_backend("manual", input_fn=lambda _: "STOP", print_fn=lambda _: None) + assert backend.name == "manual" + assert backend.curate(["A"]) == {"A": ["A"]} + + +def test_get_curation_backend_unknown_raises(): + import pytest + + with pytest.raises(ValueError): + get_curation_backend("unknown-backend") diff --git a/tests/test_database.py b/tests/test_database.py new file mode 100644 index 0000000..bb26a58 --- /dev/null +++ b/tests/test_database.py @@ -0,0 +1,111 @@ +from Bio.Seq import Seq +from Bio.SeqRecord import SeqRecord + +from candy.database import ( + DomainAssemblyRow, + DomainCurationRow, + create_database, + populate_database, + read_protein_sequences, +) + + +def make_record(id_, seq="MKV"): + return SeqRecord(Seq(seq), id=id_) + + +def test_populate_and_read_cazy_mode(tmp_path): + engine = create_database(tmp_path / "test.db") + + records = [make_record("P1_Escherichia_coli_B")] + domain_architecture = {"P1": [([1, 100], "Catalytic domain raw name")]} + domain_to_database = {"P1": {"Catalytic domain raw name": "SMART"}} + curated_domains = {"Catalytic domain": ["Catalytic domain raw name"]} + + populate_database( + engine, + records, + characterized_ids={"P1"}, + domain_architecture=domain_architecture, + domain_to_database=domain_to_database, + curated_domains=curated_domains, + is_cazy_query=True, + ) + + rows = read_protein_sequences(engine) + assert len(rows) == 1 + row = rows[0] + assert row.protein_sequence_id == "P1" + assert row.taxonomy == "Bacteria" # not the organism name -- this was the swapped-column bug + assert row.organism_name == "Escherichia coli" + assert row.characterized == "C" + assert row.domain_architecture == "Catalytic domain" + assert row.domain_locations == "[1, 100]" + assert row.domain_databases == "SMART" + + +def test_populate_database_custom_fasta_mode_uses_raw_id(tmp_path): + engine = create_database(tmp_path / "test2.db") + records = [make_record("my_custom_header")] + + populate_database( + engine, + records, + characterized_ids={}, + domain_architecture={}, + domain_to_database={}, + curated_domains={}, + is_cazy_query=False, + ) + + rows = read_protein_sequences(engine) + assert rows[0].protein_sequence_id == "my_custom_header" + assert rows[0].taxonomy == "" + assert rows[0].organism_name == "" + + +def test_populate_database_skips_duplicate_ids(tmp_path, caplog): + engine = create_database(tmp_path / "test3.db") + records = [make_record("dup"), make_record("dup")] + + populate_database( + engine, + records, + characterized_ids={}, + domain_architecture={}, + domain_to_database={}, + curated_domains={}, + is_cazy_query=False, + ) + + rows = read_protein_sequences(engine) + assert len(rows) == 1 + + +def test_populate_database_writes_assembly_and_curation_tables(tmp_path): + from sqlalchemy.orm import Session + + engine = create_database(tmp_path / "test4.db") + records = [make_record("P1_Homo_sapiens_E")] + domain_architecture = {"P1": [([1, 50], "raw")]} + domain_to_database = {"P1": {"raw": "PFAM"}} + curated_domains = {"Umbrella": ["raw"]} + + populate_database( + engine, + records, + characterized_ids={}, + domain_architecture=domain_architecture, + domain_to_database=domain_to_database, + curated_domains=curated_domains, + is_cazy_query=True, + ) + + with Session(engine) as session: + assemblies = session.query(DomainAssemblyRow).all() + curation = session.query(DomainCurationRow).all() + + assert assemblies[0].protein_domain == "Umbrella" + assert assemblies[0].protein_sequence_ids == "P1" + assert curation[0].domain_name == "Umbrella" + assert curation[0].synonymous_domain_names == "raw" diff --git a/tests/test_domains.py b/tests/test_domains.py new file mode 100644 index 0000000..2228a13 --- /dev/null +++ b/tests/test_domains.py @@ -0,0 +1,120 @@ +from candy.domains import ( + clean_domains, + group_overlapping, + map_positions_to_names, + positions_overlap, +) + +DEFAULT_PREFERENCE = ["SMART", "CDD", "PFAM", "SUPERFAMILY"] + + +def test_positions_overlap_no_overlap(): + assert positions_overlap([1, 50], [200, 250], 0.2) is False + + +def test_positions_overlap_full_containment(): + assert positions_overlap([1, 500], [100, 200], 0.2) is True + + +def test_positions_overlap_partial_above_threshold(): + # overlap region 80-100 = 20 residues; length1=100(0-100)->100 len; length2=100(80-180) + assert positions_overlap([0, 100], [80, 180], 0.1) is True + + +def test_positions_overlap_partial_below_threshold(): + assert positions_overlap([0, 100], [95, 195], 0.5) is False + + +def test_positions_overlap_zero_length_domain_is_full_identity(): + assert positions_overlap([50, 50], [40, 60], 0.99) is True + + +def test_group_overlapping_separates_distinct_domains(): + positions = [[1, 100], [150, 300]] + groups = group_overlapping(positions, 0.2) + assert groups == [[1, 100], [150, 300]] + + +def test_group_overlapping_merges_overlapping_domains(): + positions = [[1, 100], [10, 110]] + groups = group_overlapping(positions, 0.2) + assert groups == [[[1, 100], [10, 110]]] + + +def test_map_positions_to_names_duplicate_position_keeps_last_hit(): + # Two databases report the identical span [1, 100] with different names. + # Only one can survive the str(position) dict key collision; pairing by + # loop index (not value.index()) means the *current* hit's name wins, + # i.e. normal last-write-wins dict semantics rather than being stuck on + # whichever hit happened to come first in the raw API response. + position_dict = {"P1": [[1, 100], [1, 100]]} + domain_name_list = {"P1": [["SMART", "Catalytic domain"], ["PFAM", "CBM"]]} + + result = map_positions_to_names(position_dict, domain_name_list) + + assert result["P1"] == {"[1, 100]": "CBM"} + + +def test_map_positions_to_names_distinct_positions(): + position_dict = {"P1": [[1, 100], [150, 300]]} + domain_name_list = {"P1": [["SMART", "Catalytic domain"], ["PFAM", "CBM"]]} + + result = map_positions_to_names(position_dict, domain_name_list) + + assert result["P1"] == {"[1, 100]": "Catalytic domain", "[150, 300]": "CBM"} + + +def test_clean_domains_keeps_single_non_overlapping_domain(): + position_dict = {"P1": [[1, 100]]} + position_to_name = {"P1": {"[1, 100]": "Catalytic domain"}} + name_to_database = {"P1": {"Catalytic domain": "SMART"}} + + result = clean_domains( + position_dict, position_to_name, name_to_database, DEFAULT_PREFERENCE, 0.2, 800 + ) + + assert result["P1"] == [([1, 100], "Catalytic domain")] + + +def test_clean_domains_excludes_signal_peptide_and_transmembrane(): + position_dict = {"P1": [[1, 20], [30, 100]]} + position_to_name = {"P1": {"[1, 20]": "Signal peptide", "[30, 100]": "Catalytic domain"}} + name_to_database = {"P1": {"Signal peptide": "PHOBIUS", "Catalytic domain": "SMART"}} + + result = clean_domains( + position_dict, position_to_name, name_to_database, DEFAULT_PREFERENCE, 0.2, 800 + ) + + assert result["P1"] == [([30, 100], "Catalytic domain")] + + +def test_clean_domains_picks_highest_priority_database_for_overlapping_hits(): + # Same region called by both PFAM and SMART; SMART should win (rank 0). + position_dict = {"P1": [[1, 100], [5, 105]]} + position_to_name = {"P1": {"[1, 100]": "PFAM catalytic", "[5, 105]": "SMART catalytic"}} + name_to_database = {"P1": {"PFAM catalytic": "PFAM", "SMART catalytic": "SMART"}} + + result = clean_domains( + position_dict, position_to_name, name_to_database, DEFAULT_PREFERENCE, 0.2, 800 + ) + + assert len(result["P1"]) == 1 + assert result["P1"][0][1] == "SMART catalytic" + + +def test_clean_domains_filters_domains_outside_length_bounds_in_grouped_case(): + # Overlapping group where one candidate is far too long (> max_domain_length). + position_dict = {"P1": [[1, 1000], [5, 105]]} + position_to_name = {"P1": {"[1, 1000]": "Bogus long hit", "[5, 105]": "Real domain"}} + name_to_database = {"P1": {"Bogus long hit": "PFAM", "Real domain": "SMART"}} + + result = clean_domains( + position_dict, position_to_name, name_to_database, DEFAULT_PREFERENCE, 0.2, max_domain_length=200 + ) + + assert result["P1"] == [([5, 105], "Real domain")] + + +def test_clean_domains_no_domains_gives_empty_list(): + result = clean_domains({"P1": []}, {"P1": {}}, {"P1": {}}, DEFAULT_PREFERENCE, 0.2, 800) + assert result["P1"] == [] diff --git a/tests/test_external_tools.py b/tests/test_external_tools.py new file mode 100644 index 0000000..93b6775 --- /dev/null +++ b/tests/test_external_tools.py @@ -0,0 +1,38 @@ +from unittest.mock import patch + +import pytest + +from candy.external_tools import ExternalToolError, MissingDependencyError, require_binary, run_tool + + +def test_require_binary_raises_when_missing(): + with patch("candy.external_tools.shutil.which", return_value=None): + with pytest.raises(MissingDependencyError): + require_binary("nonexistent-tool") + + +def test_require_binary_returns_path_when_found(): + with patch("candy.external_tools.shutil.which", return_value="/usr/bin/mafft"): + assert require_binary("mafft") == "/usr/bin/mafft" + + +def test_run_tool_raises_with_stderr_message(tmp_path): + import sys + + script = tmp_path / "fail.py" + script.write_text("import sys; sys.stderr.write('boom'); sys.exit(2)") + + with pytest.raises(ExternalToolError, match="boom"): + run_tool([sys.executable, str(script)]) + + +def test_run_tool_succeeds_and_redirects_stdout(tmp_path): + import sys + + script = tmp_path / "ok.py" + script.write_text("print('hello')") + out_path = tmp_path / "out.txt" + + run_tool([sys.executable, str(script)], stdout_path=out_path) + + assert out_path.read_text().strip() == "hello" diff --git a/tests/test_fasta.py b/tests/test_fasta.py new file mode 100644 index 0000000..2ba6c72 --- /dev/null +++ b/tests/test_fasta.py @@ -0,0 +1,99 @@ +from Bio.SeqRecord import SeqRecord +from Bio.Seq import Seq + +from candy.fasta import ( + FormattedRecord, + count_sequences, + format_records, + parse_fasta_to_dict, + select_sequences_with_domains, + write_fasta, +) + + +def make_record(id_, description, seq): + record = SeqRecord(Seq(seq), id=id_) + record.description = description + return record + + +def test_count_sequences(tmp_path): + path = tmp_path / "seqs.fasta" + write_fasta([("a", "MKV"), ("b", "MKL")], path) + assert count_sequences(path) == 2 + + +def test_parse_fasta_to_dict_roundtrip(tmp_path): + path = tmp_path / "seqs.fasta" + write_fasta([("a", "MKV"), ("b", "MKL")], path) + assert parse_fasta_to_dict(path) == {"a": "MKV", "b": "MKL"} + + +def test_format_records_extracts_organism_and_taxonomy(): + records = [make_record("P1", "P1 some enzyme [Escherichia coli]", "MKV")] + formatted = format_records(records, {"P1": "B"}) + + assert formatted == [ + FormattedRecord(accession="P1", organism="Escherichia_coli", taxonomy_code="B", sequence="MKV") + ] + assert formatted[0].header == "P1_Escherichia_coli_B" + + +def test_format_records_handles_nested_brackets(): + records = [make_record("P2", "P2 enzyme [[Candida] auris]", "MKL")] + formatted = format_records(records, {"P2": "E"}) + assert formatted[0].organism == "Candida_auris" + + +def test_format_records_drops_sequences_without_organism_tag(): + records = [make_record("P3", "P3 enzyme, no bracket here", "MKL")] + assert format_records(records, {"P3": "E"}) == [] + + +def test_format_records_dedupes_by_sequence_and_id(): + records = [ + make_record("P4", "P4 enzyme [Homo sapiens]", "MKV"), + make_record("P5", "P5 duplicate sequence [Mus musculus]", "MKV"), + make_record("P4", "P4 duplicate id [Homo sapiens]", "MKV"), + ] + formatted = format_records(records, {"P4": "E", "P5": "E"}) + assert len(formatted) == 1 + assert formatted[0].accession == "P4" + + +def test_format_records_strips_underscore_from_accession(): + records = [make_record("AW_1234", "AW_1234 enzyme [Escherichia coli]", "MKV")] + formatted = format_records(records, {"AW_1234": "B"}) + assert formatted[0].accession == "AW1234" + + +def test_format_records_truncates_long_organism_name(): + long_name = "X" * 150 + records = [make_record("P6", f"P6 enzyme [{long_name}]", "MKV")] + formatted = format_records(records, {"P6": "B"}) + assert len(formatted[0].organism) == 100 + + +def test_select_sequences_with_domains_keeps_only_domained_cazy_records(tmp_path): + path = tmp_path / "in.fasta" + write_fasta([("P1_Org_B", "MKV"), ("P2_Org_B", "MKL")], path) + + output, no_domains = select_sequences_with_domains( + path, {"P1": [([1, 3], "Domain A")], "P2": []}, is_cazy_query=True + ) + + assert "P1_Org_B" in output + assert "P2_Org_B" not in output + assert no_domains == ["P2"] + + +def test_select_sequences_with_domains_custom_fasta_mode_uses_full_id(tmp_path): + path = tmp_path / "in.fasta" + write_fasta([("myheader", "MKV")], path) + + output, no_domains = select_sequences_with_domains( + path, {"myheader": [([1, 3], "Domain A")]}, is_cazy_query=False + ) + + assert "myheader" in output + assert no_domains == [] diff --git a/tests/test_interpro.py b/tests/test_interpro.py new file mode 100644 index 0000000..f228a60 --- /dev/null +++ b/tests/test_interpro.py @@ -0,0 +1,68 @@ +from unittest.mock import patch + +from candy.interpro import ( + DomainHit, + ProteinMatch, + ResolvedHit, + _parse_hits, + build_domain_tables, + match_lookup, +) + +SAMPLE_MATCH_XML = """ + + + + {md5} + PFAM,domain,PF00128,,,,1-100-Y + + +""" + + +def test_parse_hits_uses_first_fragment_only(): + # database, ?, domain_code, ?, ?, ?, fragments + hit_text = "PFAM,domain,PF00128,X,Y,Z,1-100-Y;150-200-Y" + hit = _parse_hits(hit_text) + assert hit == DomainHit(database="PFAM", domain_code="PF00128", start=1, end=100) + + +def test_build_domain_tables_shapes_match_domains_module_expectations(): + resolved = { + "P1": [ + ResolvedHit(database="SMART", name="Catalytic domain", start=1, end=100), + ResolvedHit(database="PFAM", name="CBM", start=150, end=300), + ] + } + + position_dict, domain_name_list, name_to_database = build_domain_tables(resolved) + + assert position_dict == {"P1": [[1, 100], [150, 300]]} + assert domain_name_list == {"P1": [["SMART", "Catalytic domain"], ["PFAM", "CBM"]]} + assert name_to_database == {"P1": {"Catalytic domain": "SMART", "CBM": "PFAM"}} + + +def test_match_lookup_fans_out_identical_sequences_to_all_ids(): + import hashlib + + sequence = "MKVLA" + md5 = hashlib.md5(sequence.encode()).hexdigest().upper() + xml_response = SAMPLE_MATCH_XML.format(md5=md5).encode() + + class FakeResponse: + def __enter__(self): + return self + + def __exit__(self, *exc): + return False + + def read(self): + return xml_response + + with patch("candy.interpro.urlopen", return_value=FakeResponse()): + matches, unmatched = match_lookup({"P1": sequence, "P2": sequence, "P3": "DIFFERENT"}) + + assert set(matches.keys()) == {"P1", "P2"} + assert matches["P1"].hits == matches["P2"].hits + assert matches["P1"].hits[0].domain_code == "PF00128" + assert unmatched == ["P3"] diff --git a/tests/test_itol.py b/tests/test_itol.py new file mode 100644 index 0000000..4386ee2 --- /dev/null +++ b/tests/test_itol.py @@ -0,0 +1,78 @@ +from Bio.Seq import Seq +from Bio.SeqRecord import SeqRecord + +from candy.itol import ( + _protein_id_for_matching, + build_characterized_style_annotation, + build_domain_annotation, +) + + +def make_record(id_, seq="MKVLAMKVLAMKVLA"): + return SeqRecord(Seq(seq), id=id_) + + +def test_protein_id_for_matching_cazy_style_header(): + assert _protein_id_for_matching("P12345_Escherichia_coli_B") == "P12345" + + +def test_protein_id_for_matching_short_prefix_before_underscore(): + assert _protein_id_for_matching("AW_12345") == "AW_12345" + + +def test_protein_id_for_matching_no_underscore_used_verbatim(): + assert _protein_id_for_matching("customheader") == "customheader" + + +def test_build_domain_annotation_assigns_shapes_colors_and_data_line(): + domain_architecture = {"P12345": [([1, 10], "Catalytic domain")]} + curated_domains = {"Catalytic domain": ["Catalytic domain"]} + records = [make_record("P12345_Escherichia_coli_B")] + + result = build_domain_annotation("myjob", domain_architecture, curated_domains, records) + + assert "DATASET_DOMAINS" in result + assert "LEGEND_LABELS\t\tCatalytic domain" in result + assert "P12345_Escherichia_coli_B\t15\tRE|1|10|#BF0F0F|Catalytic domain" in result + + +def test_build_domain_annotation_uses_umbrella_name(): + domain_architecture = {"P12345": [([1, 10], "raw name")]} + curated_domains = {"Umbrella": ["raw name"]} + records = [make_record("P12345_Org_B")] + + result = build_domain_annotation("myjob", domain_architecture, curated_domains, records) + data_section = result.split("\nDATA\n")[1] + assert "Umbrella" in data_section + assert "raw name" not in data_section + + +def test_build_characterized_style_annotation_skips_ids_without_ec_number(tmp_path): + tree_path = tmp_path / "tree.nwk" + tree_path.write_text("(P12345_Org_B:0.1,P67890_Org_B:0.1);") + + result = build_characterized_style_annotation( + tree_path, + characterized_ids=["P12345", "P67890"], + ec_numbers={"P12345": "3.2.1.4"}, + legend_title="GH5", + ) + + data_section = result.split("\nDATA\n")[1] + assert "P12345_Org_B label node" in data_section + assert "P67890_Org_B" not in data_section + + +def test_build_characterized_style_annotation_multiple_ec_numbers_uses_first(tmp_path): + tree_path = tmp_path / "tree.nwk" + tree_path.write_text("(P12345_Org_B:0.1,P67890_Org_B:0.1);") + + result = build_characterized_style_annotation( + tree_path, + characterized_ids=["P12345"], + ec_numbers={"P12345": "3.2.1.4 3.2.1.20"}, + legend_title="GH5", + ) + header_and_legend = result.split("\nDATA\n")[0] + assert "3.2.1.4" in header_and_legend + assert "3.2.1.20" not in header_and_legend diff --git a/tests/test_merge.py b/tests/test_merge.py new file mode 100644 index 0000000..685c811 --- /dev/null +++ b/tests/test_merge.py @@ -0,0 +1,62 @@ +from candy.merge import _protein_id_for_merge, merge_characterized_sequences + + +def test_protein_id_for_merge_long_accession(): + assert _protein_id_for_merge("P12345_Org_B") == "P12345" + + +def test_protein_id_for_merge_short_prefix_joins_first_two_segments(): + assert _protein_id_for_merge("AW_12345_Org_B") == "AW_12345" + + +def test_protein_id_for_merge_no_underscore(): + assert _protein_id_for_merge("plainid") == "plainid" + + +def test_merge_excludes_clustered_duplicate_and_appends_characterized(tmp_path): + clustered = tmp_path / "clustered.fasta" + characterized = tmp_path / "characterized.fasta" + output = tmp_path / "merged.fasta" + + clustered.write_text(">P12345_Org_B\nMKVAAA\n>P99999_Org2_E\nMKLBBB\n") + # Headers from candy.cazy.fetch_characterized_sequences already end in + # "_{taxcode}" -- the taxonomy code must not be duplicated on merge. + characterized.write_text(">P12345_Escherichia_coli_B\nMKVREPLACED\n") + + merge_characterized_sequences( + clustered, + characterized, + output, + verified_ids=["P12345"], + blast_choice={}, + blast_hit_sequences={}, + ) + + content = output.read_text() + assert "P99999" in content # untouched, not characterized + assert content.count("MKVAAA") == 0 # clustered copy of P12345 dropped + assert "MKVREPLACED" in content # characterized version appended + assert ">P12345_Escherichia_coli_B" in content + assert ">P12345_Escherichia_coli_B_B" not in content # no duplicated taxcode + + +def test_merge_uses_blast_hit_sequence_when_available(tmp_path): + clustered = tmp_path / "clustered.fasta" + characterized = tmp_path / "characterized.fasta" + output = tmp_path / "merged.fasta" + + clustered.write_text("") + characterized.write_text(">P12345_Org_B\nORIGINAL_SEQ\n") + + merge_characterized_sequences( + clustered, + characterized, + output, + verified_ids=["P12345"], + blast_choice={"P12345": "ACC1"}, + blast_hit_sequences={"ACC1": "BLAST_HOMOLOG_SEQ"}, + ) + + content = output.read_text() + assert "BLAST_HOMOLOG_SEQ" in content + assert "ORIGINAL_SEQ" not in content diff --git a/tests/test_network.py b/tests/test_network.py new file mode 100644 index 0000000..b3f98aa --- /dev/null +++ b/tests/test_network.py @@ -0,0 +1,28 @@ +from candy.network import build_cooccurrence_network + + +def test_build_cooccurrence_network_counts_and_edges(): + architectures = { + "P1": "CBM--Catalytic domain", + "P2": "CBM--Catalytic domain", + "P3": "Catalytic domain", + } + + graph = build_cooccurrence_network(architectures) + + assert set(graph.nodes()) == {"CBM", "Catalytic domain"} + assert graph.nodes["Catalytic domain"]["size"] == 3 + assert graph.nodes["CBM"]["size"] == 2 + assert graph.has_edge("CBM", "Catalytic domain") + assert graph.edges["CBM", "Catalytic domain"]["width"] == 2 + + +def test_build_cooccurrence_network_ignores_empty_architecture(): + graph = build_cooccurrence_network({"P1": ""}) + assert len(graph.nodes()) == 0 + + +def test_build_cooccurrence_network_single_domain_no_edges(): + graph = build_cooccurrence_network({"P1": "Catalytic domain"}) + assert set(graph.nodes()) == {"Catalytic domain"} + assert len(graph.edges()) == 0 diff --git a/tests/test_pipeline.py b/tests/test_pipeline.py new file mode 100644 index 0000000..7d8dc8f --- /dev/null +++ b/tests/test_pipeline.py @@ -0,0 +1,155 @@ +"""Integration-style tests for the pipeline orchestrator. + +External services (InterPro, NCBI/CAZy) and external CLI tools (MAFFT, +FastTree, CD-HIT/MMseqs2) are mocked at their lowest-level boundary so these +tests exercise the *wiring* between candy's modules -- correct data flow and +argument order between stages -- without needing network access or any +bioinformatics tools installed. +""" + +import xml.etree.ElementTree as ET +from unittest.mock import patch + +import pandas as pd + +from candy.config import CAZyFamilyInput, ClusteringConfig, CustomFastaInput, PipelineConfig, Taxonomy +from candy.pipeline import run_pipeline + + +class StubCurationBackend: + name = "stub" + + def curate(self, domain_names, *, family=None): + return {name: [name] for name in domain_names} + + +def _fake_match_xml(md5_list): + root = ET.Element("response") + results = ET.SubElement(root, "results") + for md5 in md5_list: + match = ET.SubElement(results, "match") + ET.SubElement(match, "proteinMD5").text = md5 + hit = ET.SubElement(match, "hit") + hit.text = "SMART,domain,SM00000,x,x,x,1-50-Y" + return root + + +def test_run_pipeline_custom_fasta_mode_without_tree(tmp_path): + fasta_path = tmp_path / "input.fasta" + fasta_path.write_text(">Protein1\nMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLA\n") + + config = PipelineConfig( + input=CustomFastaInput(fasta_path=fasta_path), + jobname="testjob", + output_dir=tmp_path / "out", + build_tree=False, + ) + + with patch("candy.interpro._query_md5_batch", side_effect=lambda batch: _fake_match_xml(batch)), patch( + "candy.interpro._fetch_entry_name", return_value="Catalytic domain" + ), patch("candy.pipeline.get_curation_backend", return_value=StubCurationBackend()): + result = run_pipeline(config) + + assert result.database_path.exists() + assert result.network_graphml_path.exists() + assert result.domain_annotation_path is None # build_tree=False + assert result.tree_path is None + assert result.sequence_count == 1 + + +def test_run_pipeline_custom_fasta_mode_with_tree(tmp_path): + fasta_path = tmp_path / "input.fasta" + fasta_path.write_text(">Protein1\nMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLA\n") + + config = PipelineConfig( + input=CustomFastaInput(fasta_path=fasta_path), + jobname="testjob2", + output_dir=tmp_path / "out", + build_tree=True, + ) + + def fake_align(self, input_fasta, output_fasta): + output_fasta.write_text(input_fasta.read_text()) + return output_fasta + + def fake_build_tree(self, alignment_fasta, output_newick): + output_newick.write_text("(Protein1:0.1);") + return output_newick + + with patch("candy.interpro._query_md5_batch", side_effect=lambda batch: _fake_match_xml(batch)), patch( + "candy.interpro._fetch_entry_name", return_value="Catalytic domain" + ), patch("candy.pipeline.get_curation_backend", return_value=StubCurationBackend()), patch( + "candy.alignment.mafft.MafftAligner.align", fake_align + ), patch("candy.phylogenetics.fasttree.FastTreeBuilder.build_tree", fake_build_tree): + result = run_pipeline(config) + + assert result.tree_path.exists() + assert result.domain_annotation_path.exists() + assert result.characterized_annotation_path is None # custom FASTA mode has no characterized enzymes + assert "DATASET_DOMAINS" in result.domain_annotation_path.read_text() + + +class FakeHandle: + """Stand-in for the object Bio.Entrez.efetch returns (used as a context manager).""" + + def __init__(self, text): + self._text = text + + def __enter__(self): + return self + + def __exit__(self, *exc): + return False + + def read(self): + return self._text + + +class IdentityClusterer: + name = "identity" + + def cluster(self, input_fasta, output_fasta, config): + output_fasta.write_text(input_fasta.read_text()) + return output_fasta + + +def test_run_pipeline_cazy_query_mode_end_to_end(tmp_path): + fasta_by_id = { + "P12345.1": ">P12345.1 alpha-glucosidase [Escherichia coli]\nMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLAMKVLA\n", + } + + def fake_efetch(db, id, rettype, retmode): + return FakeHandle("".join(fasta_by_id[i] for i in id.split(","))) + + family_page_text = "GH5\tBacteria\tEcoli\tP12345.1\tncbi\n" + characterized_table = pd.DataFrame( + {1: ["", "3.2.1.4"], 4: ["Bacteria", "P12345.1"]} + ) + + config = PipelineConfig( + input=CAZyFamilyInput(enzyme_class="GH", family_number=5, email="test@example.com", taxonomy=Taxonomy.ALL), + jobname="cazytestjob", + output_dir=tmp_path / "out", + clustering=ClusteringConfig(identity_cutoff=85), + build_tree=False, + ) + + with patch("candy.cazy.fetch_family_page", return_value=family_page_text), patch( + "candy.cazy.fetch_characterized_page", return_value=[None, characterized_table] + ), patch("candy.cazy.Entrez.efetch", side_effect=fake_efetch), patch( + "candy.interpro._query_md5_batch", side_effect=lambda batch: _fake_match_xml(batch) + ), patch("candy.interpro._fetch_entry_name", return_value="Catalytic domain"), patch( + "candy.pipeline.get_curation_backend", return_value=StubCurationBackend() + ), patch("candy.pipeline.get_clusterer", return_value=IdentityClusterer()): + result = run_pipeline(config) + + assert result.database_path.exists() + assert result.sequence_count == 1 + + import sqlite3 + + conn = sqlite3.connect(result.database_path) + rows = conn.execute("SELECT protein_sequence_id, taxonomy, organism_name, characterized FROM protein_sequences").fetchall() + conn.close() + + assert rows == [("P12345.1", "Bacteria", "Escherichia coli", "C")] diff --git a/tests/test_uniparc.py b/tests/test_uniparc.py new file mode 100644 index 0000000..ebd8b82 --- /dev/null +++ b/tests/test_uniparc.py @@ -0,0 +1,19 @@ +from candy.uniparc import write_verified_fasta + + +def test_write_verified_fasta_keeps_only_available_non_partial(tmp_path): + input_path = tmp_path / "in.fasta" + output_path = tmp_path / "out.fasta" + input_path.write_text( + ">P1 complete sequence\nMKV\n" + ">P2 partial sequence\nMKL\n" + ">P3 complete sequence\nMKA\n" + ) + + written = write_verified_fasta(input_path, output_path, available_ids={"P1", "P2"}) + + assert written == 1 + content = output_path.read_text() + assert "P1" in content + assert "P2" not in content + assert "P3" not in content From eec917be4d294df814af63e0d26bc04de63917b9 Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Thu, 6 Aug 2026 20:47:47 +0200 Subject: [PATCH 2/7] Bundle clustering/MSA/phylogenetics tools for a zero-dependency pip install Replaces the conda-only external-tool requirement with a "pip install and run" experience: FAMSA (via pyfamsa) and VeryFastTree (via veryfasttree) are now real pip dependencies with bundled binaries, and MMseqs2 is auto-downloaded and cached from the official GitHub release on first clustering run. CD-HIT/MAFFT/FastTree remain available as opt-in alternatives via the (now much smaller) conda environment for anyone who wants the original notebook's exact toolchain. Testing this for real against live binaries (not mocks) surfaced two upstream bugs, both fixed with regression tests: veryfasttree's default ext="AUTO" crashes with an access violation on Windows/AVX2 machines (now always passes an explicit ext matching the bundled binary), and MMseqs2's Windows build needs its mmseqs.bat wrapper (not mmseqs.exe directly) since its clustering workflows are internally POSIX shell scripts -- documented clearly since first use may prompt for administrator permission once, mirroring MMseqs2's own upstream guidance. Also simplifies the CLI to a single positional target (`candy GH173` / `candy sequences.fasta`, auto-detected, no `run` subcommand needed), adds CANDY_EMAIL env var support, restores jobname-collision avoidance that had been dropped from the original notebook port, and adds a --db-preference flag to reprioritize InterPro member databases without needing to spell out the full default ranking. Co-Authored-By: Claude Sonnet 5 --- README.md | 48 ++-- environment.yml | 10 +- pyproject.toml | 4 + src/candy/alignment/__init__.py | 4 + src/candy/alignment/famsa.py | 41 ++++ src/candy/cli.py | 112 ++++++--- src/candy/clustering/mmseqs2.py | 21 +- src/candy/clustering/mmseqs2_download.py | 173 ++++++++++++++ src/candy/config.py | 37 ++- src/candy/phylogenetics/__init__.py | 4 + .../phylogenetics/veryfasttree_builder.py | 39 ++++ src/candy/pipeline.py | 16 ++ tests/test_cli.py | 221 +++++++++++++++--- tests/test_clustering.py | 2 +- tests/test_config.py | 37 +++ tests/test_famsa_aligner.py | 27 +++ tests/test_mmseqs2_download.py | 96 ++++++++ tests/test_veryfasttree_builder.py | 41 ++++ 18 files changed, 851 insertions(+), 82 deletions(-) create mode 100644 src/candy/alignment/famsa.py create mode 100644 src/candy/clustering/mmseqs2_download.py create mode 100644 src/candy/phylogenetics/veryfasttree_builder.py create mode 100644 tests/test_config.py create mode 100644 tests/test_famsa_aligner.py create mode 100644 tests/test_mmseqs2_download.py create mode 100644 tests/test_veryfasttree_builder.py diff --git a/README.md b/README.md index 5a21e59..040161c 100644 --- a/README.md +++ b/README.md @@ -6,14 +6,28 @@ CANDy is a fast, FAIR and seamless protein domain analysis tool for any [CAZy](h ## Installation -CANDy depends on a handful of external bioinformatics tools (CD-HIT or MMseqs2 for clustering, MAFFT for alignment, FastTree for phylogenetics) that aren't distributed via PyPI. The supported install path is a conda environment that provides them, with CANDy itself installed via pip into that environment: +Not yet published to PyPI -- install directly from this repository for now: ```bash -conda env create -f environment.yml -conda activate candy +git clone https://github.com/PyEED/CANDy.git +cd CANDy +git checkout package-upgrade +pip install -e . ``` -This installs the external tools via bioconda and CANDy itself (editable) via pip. If you already have those tools on your PATH through some other means, you can also just `pip install -e .` directly. +(Once published, this will just be `pip install candy-cazyme`.) That's it for most users -- CANDy's default toolchain is fully bundled: + +- **Clustering**: [MMseqs2](https://github.com/soedinglab/MMseqs2) -- auto-downloaded and cached on first use (no conda needed). On Linux/macOS this just works. On **Windows**, MMseqs2's clustering workflows internally need a POSIX shell; the official Windows build handles this itself by installing a small helper (`busybox`) the first time it runs, which may ask for administrator permission **once** -- never again after that. (This mirrors upstream: MMseqs2's own docs list WSL as the recommended Windows path and this static build as the fallback for anyone who can't use WSL.) +- **MSA**: [FAMSA](https://github.com/refresh-bio/FAMSA) via [`pyfamsa`](https://github.com/althonos/pyfamsa) -- a real pip dependency, runs in-process, no download needed. +- **Phylogenetics**: [VeryFastTree](https://github.com/citiususc/veryfasttree) via [`veryfasttree`](https://github.com/citiususc/veryfasttree-python) -- also a real pip dependency, no download needed. + +If you'd rather use the original CD-HIT/MAFFT/FastTree tools instead (e.g. to reproduce results bit-for-bit against the published notebook), a conda environment with those is still provided: + +```bash +conda env create -f environment.yml # only needed for CD-HIT as a clustering alternative to MMseqs2 +conda activate candy +``` +and pass `--clustering-software cd-hit` / build a `PipelineConfig` with `alignment_tool="mafft"`, `tree_tool="fasttree"`. To also enable automated Gemini-based domain-name curation: @@ -21,19 +35,28 @@ To also enable automated Gemini-based domain-name curation: pip install -e ".[gemini]" ``` +If you'd rather not have CANDy download anything automatically (e.g. air-gapped environments), set `CANDY_NO_AUTO_DOWNLOAD=1` -- clustering will then require `mmseqs`/`cd-hit` already on PATH. + ## Usage ### Command line ```bash -# Query a CAZy family directly -candy run --jobname my_gh5_run --family GH5 --email you@example.com --tree +# Query a CAZy family directly -- TARGET is auto-detected as a family code or a file path +candy GH173 --email you@example.com --tree + +# Restrict to a taxonomic subset, and use a stricter clustering cutoff +candy GH173 --email you@example.com --taxonomy Bacteria --cluster-identity 90 + +# Reprioritize which InterPro database wins when two disagree on a domain boundary +# (only the databases you name move; everything else keeps its default order) +candy GH173 --email you@example.com --db-preference PFAM,SMART # Analyse your own FASTA file instead -candy run --jobname my_custom_run --fasta my_sequences.fasta --tree +candy my_sequences.fasta --tree ``` -Run `candy run --help` for the full list of options (taxonomy subset, clustering software/cutoff, BLAST identity threshold, max domain length, domain-overlap threshold, curation backend, ...). +`--email` falls back to the `CANDY_EMAIL` environment variable, then an interactive prompt, so `export CANDY_EMAIL=you@example.com` once and just run `candy GH173` from then on. Run `candy --help` for the full list of options. ### Python API @@ -58,7 +81,7 @@ Results are written to `{output_dir}/{jobname}/`: - FASTA files for each processing stage - A SQLite database (`{jobname}_db.db`) containing the domain annotations -- open it with [DB Browser for SQLite](https://sqlitebrowser.org/) - A protein domain co-occurrence network (`{jobname}_domain_cooccurence_network.graphml`) -- open it in [Cytoscape](https://cytoscape.org/) (yFiles Organic Layout recommended) -- If `--tree`/`build_tree=True`: a MAFFT alignment, a FastTree phylogenetic tree (Newick), and [iTOL](https://itol.embl.de/) annotation files for the domain architecture and (for CAZy family queries) characterized-enzyme activity +- If `--tree`/`build_tree=True`: a FAMSA alignment, a VeryFastTree phylogenetic tree (Newick), and [iTOL](https://itol.embl.de/) annotation files for the domain architecture and (for CAZy family queries) characterized-enzyme activity ## Acknowledgements @@ -70,10 +93,9 @@ CANDy communicates with and/or references the following separate libraries, pack - [sqlitebrowser](https://sqlitebrowser.org/) - [SQLAlchemy](https://www.sqlalchemy.org/) - [requests](https://requests.readthedocs.io/en/latest/) -- [CD-HIT](https://academic.oup.com/bioinformatics/article/22/13/1658/194225?login=true) -- [MMseqs2](https://www.nature.com/articles/nbt.3988) -- [MAFFT](https://academic.oup.com/nar/article/30/14/3059/2904316?login=true) -- [FastTree](http://www.microbesonline.org/fasttree/) +- [MMseqs2](https://www.nature.com/articles/nbt.3988) / [CD-HIT](https://academic.oup.com/bioinformatics/article/22/13/1658/194225?login=true) (clustering) +- [FAMSA](https://academic.oup.com/nar/article/44/16/e121/2468101) via [pyfamsa](https://github.com/althonos/pyfamsa) / [MAFFT](https://academic.oup.com/nar/article/30/14/3059/2904316?login=true) (alignment) +- [VeryFastTree](https://academic.oup.com/bioinformatics/article/36/17/4658/5850991) via [veryfasttree](https://github.com/citiususc/veryfasttree-python) / [FastTree](http://www.microbesonline.org/fasttree/) (phylogenetics) - [NetworkX](https://networkx.org/) - [Matplotlib](https://matplotlib.org/) diff --git a/environment.yml b/environment.yml index 5084977..337341c 100644 --- a/environment.yml +++ b/environment.yml @@ -1,3 +1,10 @@ +# Optional. CANDy's default toolchain (MMseqs2 for clustering, FAMSA for MSA, +# VeryFastTree for phylogenetics) ships via `pip install candy-cazyme` alone -- +# MSA/phylogenetics are bundled Python packages, and MMseqs2 is auto-downloaded +# on first use. This environment is only needed if you specifically want +# CD-HIT as the clustering backend instead of MMseqs2: CD-HIT has no official +# Windows build and no pip-installable bindings, so (unlike MMseqs2) it can't +# be auto-downloaded and must come from conda/bioconda. name: candy channels: - bioconda @@ -5,9 +12,6 @@ channels: dependencies: - python>=3.10 - pip - - mmseqs2 - cd-hit - - mafft - - fasttree - pip: - -e . diff --git a/pyproject.toml b/pyproject.toml index 1b9b3ba..873df54 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,6 +29,9 @@ dependencies = [ "numpy>=1.26", "typer>=0.12", "lxml>=5.0", + "pyfamsa>=0.5", + "veryfasttree>=4.0", + "platformdirs>=4.0", ] [project.optional-dependencies] @@ -46,6 +49,7 @@ where = ["src"] [tool.pytest.ini_options] testpaths = ["tests"] +addopts = "-m 'not integration'" markers = [ "integration: requires network access or external CLI tools (mafft, fasttree, mmseqs2, cd-hit)", ] diff --git a/src/candy/alignment/__init__.py b/src/candy/alignment/__init__.py index 5312b5a..e6b459f 100644 --- a/src/candy/alignment/__init__.py +++ b/src/candy/alignment/__init__.py @@ -15,6 +15,10 @@ def align(self, input_fasta: Path, output_fasta: Path) -> Path: def get_alignment_tool(name: str) -> AlignmentTool: + if name == "famsa": + from candy.alignment.famsa import FamsaAligner + + return FamsaAligner() if name == "mafft": from candy.alignment.mafft import MafftAligner diff --git a/src/candy/alignment/famsa.py b/src/candy/alignment/famsa.py new file mode 100644 index 0000000..610d6ee --- /dev/null +++ b/src/candy/alignment/famsa.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +import logging +from pathlib import Path + +from Bio import SeqIO +from pyfamsa import Aligner, Sequence + +logger = logging.getLogger(__name__) + + +class FamsaAligner: + """MSA via FAMSA (through the ``pyfamsa`` bindings). + + Runs in-process (no subprocess, no external binary) -- ``pyfamsa`` ships + prebuilt wheels for Linux/macOS/Windows, making this the default, + zero-install alignment backend. :class:`candy.alignment.mafft.MafftAligner` + remains available for anyone who specifically wants MAFFT. + """ + + name = "famsa" + + def __init__(self, threads: int = 0) -> None: + # threads=0 lets FAMSA pick a sensible default (all available cores). + self.threads = threads + + def align(self, input_fasta: Path, output_fasta: Path) -> Path: + with open(input_fasta) as handle: + records = list(SeqIO.parse(handle, "fasta")) + + sequences = [Sequence(record.id.encode(), str(record.seq).encode()) for record in records] + + logger.info("Aligning %d sequences with FAMSA.", len(sequences)) + aligner = Aligner(threads=self.threads) + alignment = aligner.align(sequences) + + with open(output_fasta, "w") as out: + for gapped in alignment: + out.write(f">{gapped.id.decode()}\n{gapped.sequence.decode()}\n") + + return output_fasta diff --git a/src/candy/cli.py b/src/candy/cli.py index 0976710..b1d3e2a 100644 --- a/src/candy/cli.py +++ b/src/candy/cli.py @@ -1,15 +1,17 @@ -"""Command-line entry point: `candy run ...`.""" +"""Command-line entry point: `candy GH173` / `candy my_sequences.fasta`.""" from __future__ import annotations import logging import re +import sys from pathlib import Path from typing import Optional import typer from candy.config import ( + DEFAULT_DATABASE_PREFERENCE, CAZyFamilyInput, ClusteringConfig, ClusteringSoftware, @@ -18,46 +20,46 @@ DomainCleaningConfig, PipelineConfig, Taxonomy, + reorder_database_preference, ) from candy.pipeline import run_pipeline -app = typer.Typer(help="CANDy: automated analysis of domain architectures in carbohydrate-active enzymes.") +app = typer.Typer( + help="CANDy: automated analysis of domain architectures in carbohydrate-active enzymes.", + add_completion=False, +) _FAMILY_RE = re.compile(r"^([A-Za-z]+)(\d+)(?:_(\d+))?$") +_SANITIZE_RE = re.compile(r"\W+") -@app.callback() -def _main() -> None: - """CANDy: automated analysis of domain architectures in carbohydrate-active enzymes. - - A no-op callback: Typer collapses a Typer() app with exactly one - @app.command() into a single top-level command (dropping the command - name entirely), which would silently break `candy run ...`. Registering - this callback keeps `run` as a real, required subcommand. - """ - - -def _parse_family(family: str) -> tuple[str, int, str | None]: - match = _FAMILY_RE.match(family) +def _try_parse_family(value: str) -> tuple[str, int, str | None] | None: + match = _FAMILY_RE.match(value) if not match: - raise typer.BadParameter( - f"Could not parse '{family}' as a CAZy family (expected e.g. 'GH5' or 'GH5_1')." - ) + return None enzyme_class, number, subfamily = match.groups() return enzyme_class, int(number), subfamily +def _sanitize_jobname(value: str) -> str: + return _SANITIZE_RE.sub("", "".join(value.split())) + + @app.command() -def run( - jobname: str = typer.Option(..., help="Job name; results are written to output-dir/jobname."), - family: Optional[str] = typer.Option( - None, help="CAZy family to query, e.g. 'GH5' or 'GH5_1'. Mutually exclusive with --fasta." +def main( + target: str = typer.Argument( + ..., help="CAZy family to query (e.g. 'GH173' or 'GH5_1'), or a path to a custom FASTA file." + ), + jobname: Optional[str] = typer.Option( + None, help="Job name; results are written to output-dir/jobname. Defaults to TARGET." ), - fasta: Optional[Path] = typer.Option( - None, exists=True, help="Custom FASTA file to analyse instead of querying CAZy." + email: Optional[str] = typer.Option( + None, + envvar="CANDY_EMAIL", + help="Email for NCBI Entrez (required for CAZy family queries, NCBI policy). " + "Falls back to the CANDY_EMAIL environment variable, then an interactive prompt.", ), - email: Optional[str] = typer.Option(None, help="Email for NCBI Entrez. Required with --family."), - taxonomy: Taxonomy = typer.Option(Taxonomy.ALL, help="Taxonomic subset to restrict a --family query to."), + taxonomy: Taxonomy = typer.Option(Taxonomy.ALL, help="Taxonomic subset to restrict a family query to."), output_dir: Path = typer.Option(Path("."), help="Directory results are written under."), clustering_software: ClusteringSoftware = typer.Option( ClusteringSoftware.MMSEQS2, help="Sequence clustering backend." @@ -70,6 +72,14 @@ def run( overlap_percentage: int = typer.Option( 20, help="Overlap threshold above which two hits are considered the same domain, percent." ), + db_preference: Optional[str] = typer.Option( + None, + "--db-preference", + help="Comma-separated InterPro member-database names to prioritize when two databases " + "disagree on a domain boundary, e.g. 'PFAM,SMART'. Named databases move to the front in " + "the given order; unlisted ones keep their default relative order after. Valid names: " + + ", ".join(DEFAULT_DATABASE_PREFERENCE) + ".", + ), build_tree: bool = typer.Option(False, "--tree/--no-tree", help="Run MSA + phylogenetics + iTOL export."), curation_backend: str = typer.Option( "manual", help="Domain-name curation backend: 'manual' or 'gemini'." @@ -79,16 +89,39 @@ def run( ), verbose: bool = typer.Option(False, "--verbose", "-v", help="Enable debug logging."), ) -> None: - """Run the CANDy pipeline end-to-end.""" + """Run the CANDy pipeline end-to-end. + + TARGET is either a CAZy family code ('GH173', 'GH5_1', ...) or a path to + a FASTA file -- whichever it looks like determines the input mode, e.g.: + + candy GH173 --email you@example.com --tree + + candy my_sequences.fasta --tree + """ logging.basicConfig(level=logging.DEBUG if verbose else logging.INFO, format="%(message)s") - if (family is None) == (fasta is None): - raise typer.BadParameter("Specify exactly one of --family or --fasta.") + target_path = Path(target) + if target_path.is_file(): + pipeline_input: CAZyFamilyInput | CustomFastaInput = CustomFastaInput(fasta_path=target_path) + default_jobname = target_path.stem + else: + parsed = _try_parse_family(target) + if parsed is None: + raise typer.BadParameter( + f"'{target}' is neither an existing FASTA file nor a valid CAZy family code " + "(expected e.g. 'GH5' or 'GH5_1')." + ) + enzyme_class, family_number, subfamily = parsed - if family is not None: if not email: - raise typer.BadParameter("--email is required when using --family (NCBI Entrez requires it).") - enzyme_class, family_number, subfamily = _parse_family(family) + if sys.stdin.isatty(): + email = typer.prompt("NCBI Entrez requires an email address") + else: + raise typer.BadParameter( + "An email is required for CAZy family queries (NCBI Entrez policy). " + "Pass --email, or set the CANDY_EMAIL environment variable." + ) + pipeline_input = CAZyFamilyInput( enzyme_class=enzyme_class, family_number=family_number, @@ -96,16 +129,25 @@ def run( email=email, taxonomy=taxonomy, ) + default_jobname = target + + if db_preference: + try: + database_preference = reorder_database_preference(db_preference.split(",")) + except ValueError as exc: + raise typer.BadParameter(str(exc)) from exc else: - pipeline_input = CustomFastaInput(fasta_path=fasta) + database_preference = list(DEFAULT_DATABASE_PREFERENCE) config = PipelineConfig( input=pipeline_input, - jobname=jobname, + jobname=_sanitize_jobname(jobname or default_jobname), output_dir=output_dir, clustering=ClusteringConfig(software=clustering_software, identity_cutoff=cluster_identity), domain_cleaning=DomainCleaningConfig( - max_domain_length=max_domain_length, overlap_percentage=overlap_percentage + max_domain_length=max_domain_length, + overlap_percentage=overlap_percentage, + database_preference=database_preference, ), curation=CurationConfig(backend=curation_backend, api_key=curation_api_key), blast_identity_threshold=blast_identity, diff --git a/src/candy/clustering/mmseqs2.py b/src/candy/clustering/mmseqs2.py index bca212f..1b60470 100644 --- a/src/candy/clustering/mmseqs2.py +++ b/src/candy/clustering/mmseqs2.py @@ -1,21 +1,38 @@ from __future__ import annotations import logging +import platform import shutil import tempfile from pathlib import Path +from candy.clustering.mmseqs2_download import resolve_mmseqs2_binary from candy.config import ClusteringConfig -from candy.external_tools import require_binary, run_tool +from candy.external_tools import run_tool logger = logging.getLogger(__name__) class Mmseqs2Clusterer: + """Sequence clustering via MMseqs2. + + Prefers an ``mmseqs`` already on PATH; otherwise transparently downloads + and caches a static build for the current platform on first use (see + :mod:`candy.clustering.mmseqs2_download`). This is the default clustering + backend precisely because that auto-download is possible for MMseqs2 but + not for CD-HIT (source-only, no Windows build). + """ + name = "mmseqs2" def cluster(self, input_fasta: Path, output_fasta: Path, config: ClusteringConfig) -> Path: - binary = require_binary("mmseqs") + binary = resolve_mmseqs2_binary() + if platform.system() == "Windows" and binary.lower().endswith(".bat"): + logger.info( + "Clustering on Windows uses MMseqs2's mmseqs.bat wrapper, which may ask for " + "administrator permission once (to install a small POSIX-shell helper it " + "needs internally). This only happens on the very first clustering run." + ) min_seq_id = config.identity_cutoff / 100 min_coverage = config.mmseqs_min_coverage / 100 diff --git a/src/candy/clustering/mmseqs2_download.py b/src/candy/clustering/mmseqs2_download.py new file mode 100644 index 0000000..25871a9 --- /dev/null +++ b/src/candy/clustering/mmseqs2_download.py @@ -0,0 +1,173 @@ +"""Locate or auto-download a static MMseqs2 binary. + +MMseqs2 has no pip-native Python bindings, but -- unlike CD-HIT, which is +source-only and has no Windows build -- upstream publishes genuine static +binaries for every major platform (Linux, macOS, Windows, ARM64, PPC64LE). +That makes it possible to give CANDy a "just pip install and run" clustering +backend: on first use, if ``mmseqs`` isn't already on PATH, download and +cache the correct static build from the official GitHub release, and reuse +it on every later run. + +This reaches out to the network and runs a downloaded binary automatically, +which is a real trust decision, not just a convenience one: the exact +release version is pinned below, downloads only ever come from the official +soedinglab/MMseqs2 GitHub release over HTTPS, and the whole thing can be +disabled with the CANDY_NO_AUTO_DOWNLOAD environment variable (falling back +to requiring `mmseqs` on PATH, e.g. via the bundled conda environment -- +except on Windows, see below, where bioconda publishes no `mmseqs2` build at +all). + +Windows note: MMseqs2's clustering workflows (``cluster``, ``linclust``, and +their ``easy-*`` wrappers) are implemented as shell scripts even inside the +compiled binary, not pure compiled code. On Linux/macOS that's invisible +since a POSIX shell is always present. On Windows there usually isn't one, +so the official release instead ships a `mmseqs.bat` wrapper plus a bundled +`busybox.exe`; on first invocation, `mmseqs.bat` tries to silently install +busybox's POSIX utilities via symlinks, and Windows only allows unprivileged +symlink creation with Developer Mode enabled -- otherwise `mmseqs.bat` falls +back to popping a UAC administrator-elevation prompt. This is documented +upstream (the MMseqs2 user guide lists WSL as the *recommended* Windows +install path, with this static build as a fallback for anyone who can't use +WSL). CANDy resolves to `mmseqs.bat` (not `mmseqs.exe` directly) on Windows +so that self-install happens automatically; expect a one-time UAC prompt on +the first clustering run on a Windows machine without Developer Mode, never +again after that (the installed helpers are cached). +""" + +from __future__ import annotations + +import logging +import os +import platform +import shutil +import stat +import tarfile +import zipfile +from pathlib import Path + +import requests +from platformdirs import user_cache_dir + +from candy.external_tools import MissingDependencyError, find_binary + +logger = logging.getLogger(__name__) + +MMSEQS2_VERSION = "18-8cc5c" +_RELEASE_BASE_URL = f"https://github.com/soedinglab/MMseqs2/releases/download/{MMSEQS2_VERSION}" + +_ASSET_BY_PLATFORM = { + ("Linux", "x86_64"): "mmseqs-linux-sse41.tar.gz", + ("Linux", "aarch64"): "mmseqs-linux-arm64.tar.gz", + ("Linux", "arm64"): "mmseqs-linux-arm64.tar.gz", + ("Linux", "ppc64le"): "mmseqs-linux-ppc64le-power9.tar.gz", + ("Darwin", "x86_64"): "mmseqs-osx-universal.tar.gz", + ("Darwin", "arm64"): "mmseqs-osx-universal.tar.gz", + ("Windows", "x86_64"): "mmseqs-win64.zip", + ("Windows", "AMD64"): "mmseqs-win64.zip", +} + +NO_AUTO_DOWNLOAD_ENV = "CANDY_NO_AUTO_DOWNLOAD" + + +def _cache_dir() -> Path: + return Path(user_cache_dir("candy")) / "mmseqs2" / MMSEQS2_VERSION + + +def _binary_name() -> str: + # On Windows we deliberately resolve to the .bat wrapper, not mmseqs.exe + # directly -- see the module docstring for why (POSIX-shell workflow + # scripts + one-time busybox self-install). + return "mmseqs.bat" if platform.system() == "Windows" else "mmseqs" + + +def _find_extracted_binary(root: Path) -> Path | None: + if not root.is_dir(): + return None + target = _binary_name() + for path in root.rglob(target): + if path.is_file(): + return path + return None + + +def _download_and_extract(asset: str, cache_dir: Path) -> Path: + url = f"{_RELEASE_BASE_URL}/{asset}" + logger.info( + "Downloading MMseqs2 %s (%s) for clustering -- this happens once and is cached at %s", + MMSEQS2_VERSION, asset, cache_dir, + ) + + cache_dir.mkdir(parents=True, exist_ok=True) + archive_path = cache_dir / asset + extract_dir = cache_dir / "extracted" + if extract_dir.exists(): + shutil.rmtree(extract_dir) + + with requests.get(url, stream=True, timeout=120) as response: + response.raise_for_status() + with open(archive_path, "wb") as f: + for chunk in response.iter_content(chunk_size=1 << 20): + f.write(chunk) + + extract_dir.mkdir(exist_ok=True) + if asset.endswith(".zip"): + with zipfile.ZipFile(archive_path) as zf: + zf.extractall(extract_dir) + else: + with tarfile.open(archive_path) as tf: + tf.extractall(extract_dir) + archive_path.unlink() + + binary = _find_extracted_binary(extract_dir) + if binary is None: + raise MissingDependencyError( + f"Downloaded MMseqs2 archive '{asset}' did not contain a '{_binary_name()}' executable." + ) + + if platform.system() != "Windows": + binary.chmod(binary.stat().st_mode | stat.S_IXUSR | stat.S_IXGRP | stat.S_IXOTH) + + return binary + + +def resolve_mmseqs2_binary() -> str: + """Return a path to a usable ``mmseqs`` binary, downloading one if needed. + + Resolution order: PATH, then the local cache, then (unless disabled via + the CANDY_NO_AUTO_DOWNLOAD env var) a fresh download of the pinned + MMseqs2 release for the current platform. + """ + on_path = find_binary("mmseqs") + if on_path: + return on_path + + cache_dir = _cache_dir() + cached = _find_extracted_binary(cache_dir) + if cached: + return str(cached) + + if os.environ.get(NO_AUTO_DOWNLOAD_ENV): + raise MissingDependencyError( + "mmseqs was not found on PATH and automatic download is disabled " + f"({NO_AUTO_DOWNLOAD_ENV} is set). Install MMseqs2 yourself (e.g. via " + "`conda env create -f environment.yml`) or unset that variable." + ) + + platform_key = (platform.system(), platform.machine()) + asset = _ASSET_BY_PLATFORM.get(platform_key) + if asset is None: + raise MissingDependencyError( + f"No known MMseqs2 static build for platform {platform_key}. Install " + "MMseqs2 yourself and ensure it is on PATH." + ) + + try: + binary = _download_and_extract(asset, cache_dir) + except requests.RequestException as exc: + raise MissingDependencyError( + f"Failed to download MMseqs2 automatically: {exc}. Install it yourself " + "(e.g. via `conda env create -f environment.yml`) or check your network connection. " + f"You can also set {NO_AUTO_DOWNLOAD_ENV}=1 to disable this download attempt." + ) from exc + + return str(binary) diff --git a/src/candy/config.py b/src/candy/config.py index 18d5f13..b0cdfae 100644 --- a/src/candy/config.py +++ b/src/candy/config.py @@ -10,6 +10,7 @@ from __future__ import annotations +from collections.abc import Sequence from dataclasses import dataclass, field from enum import Enum from pathlib import Path @@ -61,6 +62,38 @@ class ClusteringSoftware(str, Enum): ] +def reorder_database_preference( + priority: Sequence[str], base: Sequence[str] = DEFAULT_DATABASE_PREFERENCE +) -> list[str]: + """Move the named databases (in the given order) to the front of ``base``. + + Databases not named in ``priority`` keep their existing relative order, + appended after the ones that were moved -- so callers only need to name + the databases they actually want to reprioritize, not repeat the full + 13-item list every time. + + Raises: + ValueError: if ``priority`` contains a name not present in ``base``. + """ + normalized_priority = [name.strip().upper() for name in priority] + valid = set(base) + unknown = [name for name in normalized_priority if name not in valid] + if unknown: + raise ValueError( + f"Unknown database name(s): {', '.join(unknown)}. Valid names: {', '.join(base)}." + ) + + front: list[str] = [] + seen: set[str] = set() + for name in normalized_priority: + if name not in seen: + front.append(name) + seen.add(name) + + rest = [name for name in base if name not in seen] + return front + rest + + @dataclass class CAZyFamilyInput: """Query a CAZy family/subfamily directly (input option 1 in the notebook).""" @@ -132,8 +165,8 @@ class PipelineConfig: # BLAST fallback when a characterized sequence has no UniParc entry. build_tree: bool = False - alignment_tool: str = "mafft" - tree_tool: str = "fasttree" + alignment_tool: str = "famsa" + tree_tool: str = "veryfasttree" def __post_init__(self) -> None: self.output_dir = Path(self.output_dir) diff --git a/src/candy/phylogenetics/__init__.py b/src/candy/phylogenetics/__init__.py index 2d689a9..902b385 100644 --- a/src/candy/phylogenetics/__init__.py +++ b/src/candy/phylogenetics/__init__.py @@ -15,6 +15,10 @@ def build_tree(self, alignment_fasta: Path, output_newick: Path) -> Path: def get_tree_builder(name: str) -> TreeBuilder: + if name == "veryfasttree": + from candy.phylogenetics.veryfasttree_builder import VeryFastTreeBuilder + + return VeryFastTreeBuilder() if name == "fasttree": from candy.phylogenetics.fasttree import FastTreeBuilder diff --git a/src/candy/phylogenetics/veryfasttree_builder.py b/src/candy/phylogenetics/veryfasttree_builder.py new file mode 100644 index 0000000..b87a9c8 --- /dev/null +++ b/src/candy/phylogenetics/veryfasttree_builder.py @@ -0,0 +1,39 @@ +from __future__ import annotations + +import logging +from pathlib import Path + +import veryfasttree + +logger = logging.getLogger(__name__) + +_SIMD_TO_EXT = {"avx512f": "AVX512", "avx2": "AVX2", "sse2": "SSE"} + + +class VeryFastTreeBuilder: + """Phylogenetics via VeryFastTree (through the ``veryfasttree`` bindings). + + Runs a bundled, precompiled binary shipped inside the ``veryfasttree`` + wheel -- no external FastTree/VeryFastTree installation needed. This is + the default, zero-install phylogenetics backend; + :class:`candy.phylogenetics.fasttree.FastTreeBuilder` remains available + for anyone who specifically wants the original FastTree binary. + + .. note:: + ``veryfasttree.run()`` defaults to ``ext="AUTO"`` (let the binary + self-detect vector extensions). During testing, both ``AUTO`` and the + explicit scalar path ``ext="NONE"`` reliably crashed (access + violation) on a Windows/AVX2 machine, while explicitly requesting the + same SIMD tier the package itself already selected for the bundled + binary (``SSE``/``AVX2``/``AVX512``, matched via ``veryfasttree._simd``) + ran correctly. This builder always passes that explicit ``ext`` + rather than relying on the library's ``AUTO`` default. + """ + + name = "veryfasttree" + + def build_tree(self, alignment_fasta: Path, output_newick: Path) -> Path: + ext = _SIMD_TO_EXT.get(getattr(veryfasttree, "_simd", None), "SSE") + logger.info("Building phylogenetic tree with VeryFastTree (ext=%s).", ext) + veryfasttree.run(str(alignment_fasta), out=str(output_newick), quiet=True, nopr=True, ext=ext) + return output_newick diff --git a/src/candy/pipeline.py b/src/candy/pipeline.py index c431c7f..8b7da00 100644 --- a/src/candy/pipeline.py +++ b/src/candy/pipeline.py @@ -67,7 +67,23 @@ def _prepare_domain_detection_input( return proteins +def _unique_jobname(output_dir: Path, jobname: str) -> str: + """Suffix jobname with an incrementing number if it's already taken. + + Avoids silently mixing a new run's output files into a previous run's + directory (the notebook did the same, appending _0, _1, ... on collision). + """ + if not (output_dir / jobname).exists(): + return jobname + + n = 0 + while (output_dir / f"{jobname}_{n}").exists(): + n += 1 + return f"{jobname}_{n}" + + def run_pipeline(config: PipelineConfig) -> PipelineResult: + config.jobname = _unique_jobname(config.output_dir, config.jobname) jobname_dir = config.output_dir / config.jobname jobname_dir.mkdir(parents=True, exist_ok=True) diff --git a/tests/test_cli.py b/tests/test_cli.py index 6fece99..c00cafa 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,55 +1,224 @@ from typer.testing import CliRunner -from candy.cli import _parse_family, app +from candy.cli import _sanitize_jobname, _try_parse_family, app runner = CliRunner() -def test_parse_family_no_subfamily(): - assert _parse_family("GH5") == ("GH", 5, None) +def test_try_parse_family_no_subfamily(): + assert _try_parse_family("GH5") == ("GH", 5, None) -def test_parse_family_with_subfamily(): - assert _parse_family("GH5_1") == ("GH", 5, "1") +def test_try_parse_family_with_subfamily(): + assert _try_parse_family("GH5_1") == ("GH", 5, "1") -def test_parse_family_invalid_raises(): - import typer +def test_try_parse_family_invalid_returns_none(): + assert _try_parse_family("not-a-family") is None - try: - _parse_family("not-a-family") - assert False, "expected BadParameter" - except typer.BadParameter: - pass + +def test_sanitize_jobname_strips_whitespace_and_special_chars(): + assert _sanitize_jobname(" my job! ") == "myjob" def test_cli_help(): result = runner.invoke(app, ["--help"]) assert result.exit_code == 0 - assert "run" in result.output.lower() + assert "TARGET" in result.output or "target" in result.output.lower() -def test_cli_run_help(): - result = runner.invoke(app, ["run", "--help"]) - assert result.exit_code == 0 - assert "--family" in result.output - assert "--fasta" in result.output +def test_cli_bare_invocation_is_candy_target_no_subcommand_needed(): + # `candy GH5` should work directly -- no `run` subcommand required. + result = runner.invoke(app, ["not-a-real-family-and-not-a-file"]) + result_output = result.output + for border_char in "─│┌┐└┘": + result_output = result_output.replace(border_char, " ") + normalized = " ".join(result_output.split()) + assert result.exit_code != 0 + assert "neither an existing FASTA file" in normalized + assert "valid CAZy family code" in normalized -def test_cli_run_requires_family_or_fasta(): - result = runner.invoke(app, ["run", "--jobname", "test"]) +def test_cli_family_query_without_email_and_non_interactive_fails(): + result = runner.invoke(app, ["GH5"], input="") assert result.exit_code != 0 + assert "email" in result.output.lower() -def test_cli_run_rejects_both_family_and_fasta(tmp_path): - fasta_path = tmp_path / "in.fasta" +def test_cli_fasta_target_autodetected(tmp_path, monkeypatch): + fasta_path = tmp_path / "my_seqs.fasta" fasta_path.write_text(">a\nMKV\n") + + captured = {} + + def fake_run_pipeline(config): + captured["config"] = config + from candy.pipeline import PipelineResult + + return PipelineResult( + jobname_dir=tmp_path / "out", + database_path=tmp_path / "out" / "db.db", + network_graphml_path=tmp_path / "out" / "net.graphml", + domain_annotation_path=None, + characterized_annotation_path=None, + alignment_path=None, + tree_path=None, + sequence_count=1, + ) + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + + result = runner.invoke(app, [str(fasta_path)]) + + assert result.exit_code == 0, result.output + from candy.config import CustomFastaInput + + assert isinstance(captured["config"].input, CustomFastaInput) + assert captured["config"].jobname == _sanitize_jobname(fasta_path.stem) + + +def test_cli_family_target_autodetected_with_email_flag(tmp_path, monkeypatch): + captured = {} + + def fake_run_pipeline(config): + captured["config"] = config + from candy.pipeline import PipelineResult + + return PipelineResult( + jobname_dir=tmp_path / "out", + database_path=tmp_path / "out" / "db.db", + network_graphml_path=tmp_path / "out" / "net.graphml", + domain_annotation_path=None, + characterized_annotation_path=None, + alignment_path=None, + tree_path=None, + sequence_count=1, + ) + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + + result = runner.invoke(app, ["GH173", "--email", "you@example.com"]) + + assert result.exit_code == 0, result.output + from candy.config import CAZyFamilyInput + + config_input = captured["config"].input + assert isinstance(config_input, CAZyFamilyInput) + assert config_input.family == "GH173" + assert config_input.email == "you@example.com" + assert captured["config"].jobname == "GH173" + + +def test_cli_email_falls_back_to_env_var(tmp_path, monkeypatch): + captured = {} + + def fake_run_pipeline(config): + captured["config"] = config + from candy.pipeline import PipelineResult + + return PipelineResult( + jobname_dir=tmp_path / "out", + database_path=tmp_path / "out" / "db.db", + network_graphml_path=tmp_path / "out" / "net.graphml", + domain_annotation_path=None, + characterized_annotation_path=None, + alignment_path=None, + tree_path=None, + sequence_count=1, + ) + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + monkeypatch.setenv("CANDY_EMAIL", "env@example.com") + + result = runner.invoke(app, ["GH173"]) + + assert result.exit_code == 0, result.output + assert captured["config"].input.email == "env@example.com" + + +def test_cli_explicit_jobname_overrides_default(tmp_path, monkeypatch): + captured = {} + + def fake_run_pipeline(config): + captured["config"] = config + from candy.pipeline import PipelineResult + + return PipelineResult( + jobname_dir=tmp_path / "out", + database_path=tmp_path / "out" / "db.db", + network_graphml_path=tmp_path / "out" / "net.graphml", + domain_annotation_path=None, + characterized_annotation_path=None, + alignment_path=None, + tree_path=None, + sequence_count=1, + ) + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + + result = runner.invoke(app, ["GH173", "--email", "you@example.com", "--jobname", "custom_name"]) + + assert result.exit_code == 0, result.output + assert captured["config"].jobname == "custom_name" + + +def _fake_pipeline_result(tmp_path): + from candy.pipeline import PipelineResult + + return PipelineResult( + jobname_dir=tmp_path / "out", + database_path=tmp_path / "out" / "db.db", + network_graphml_path=tmp_path / "out" / "net.graphml", + domain_annotation_path=None, + characterized_annotation_path=None, + alignment_path=None, + tree_path=None, + sequence_count=1, + ) + + +def test_cli_db_preference_reorders_default(tmp_path, monkeypatch): + captured = {} + + def fake_run_pipeline(config): + captured["config"] = config + return _fake_pipeline_result(tmp_path) + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + result = runner.invoke( - app, ["run", "--jobname", "test", "--family", "GH5", "--fasta", str(fasta_path), "--email", "a@b.com"] + app, ["GH173", "--email", "you@example.com", "--db-preference", "PFAM,SMART"] ) - assert result.exit_code != 0 + + assert result.exit_code == 0, result.output + preference = captured["config"].domain_cleaning.database_preference + assert preference[:2] == ["PFAM", "SMART"] + assert "CDD" in preference # untouched entries still present -def test_cli_run_family_without_email_fails(): - result = runner.invoke(app, ["run", "--jobname", "test", "--family", "GH5"]) +def test_cli_db_preference_unknown_name_is_a_clean_error(tmp_path, monkeypatch): + monkeypatch.setattr("candy.cli.run_pipeline", lambda config: _fake_pipeline_result(tmp_path)) + + result = runner.invoke( + app, ["GH173", "--email", "you@example.com", "--db-preference", "NOTAREALDB"] + ) + assert result.exit_code != 0 + assert "Unknown database name" in result.output + + +def test_cli_no_db_preference_uses_default_order(tmp_path, monkeypatch): + from candy.config import DEFAULT_DATABASE_PREFERENCE + + captured = {} + + def fake_run_pipeline(config): + captured["config"] = config + return _fake_pipeline_result(tmp_path) + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + + result = runner.invoke(app, ["GH173", "--email", "you@example.com"]) + + assert result.exit_code == 0, result.output + assert captured["config"].domain_cleaning.database_preference == DEFAULT_DATABASE_PREFERENCE diff --git a/tests/test_clustering.py b/tests/test_clustering.py index 5a4577a..7449204 100644 --- a/tests/test_clustering.py +++ b/tests/test_clustering.py @@ -43,7 +43,7 @@ def fake_run_tool(args, **kwargs): with open(f"{prefix}_rep_seq.fasta", "w") as f: f.write(">rep\nMKV\n") - with patch("candy.clustering.mmseqs2.require_binary", return_value="/usr/bin/mmseqs"), patch( + with patch("candy.clustering.mmseqs2.resolve_mmseqs2_binary", return_value="/usr/bin/mmseqs"), patch( "candy.clustering.mmseqs2.run_tool", side_effect=fake_run_tool ) as mock_run: Mmseqs2Clusterer().cluster(input_fasta, output_fasta, config) diff --git a/tests/test_config.py b/tests/test_config.py new file mode 100644 index 0000000..fedcee6 --- /dev/null +++ b/tests/test_config.py @@ -0,0 +1,37 @@ +import pytest + +from candy.config import DEFAULT_DATABASE_PREFERENCE, reorder_database_preference + + +def test_reorder_moves_named_databases_to_front_in_given_order(): + result = reorder_database_preference(["PFAM", "SMART"]) + assert result[0] == "PFAM" + assert result[1] == "SMART" + assert set(result) == set(DEFAULT_DATABASE_PREFERENCE) + assert len(result) == len(DEFAULT_DATABASE_PREFERENCE) + + +def test_reorder_keeps_relative_order_of_unlisted_databases(): + result = reorder_database_preference(["PROSITE"]) + unlisted = [name for name in result if name != "PROSITE"] + assert unlisted == [name for name in DEFAULT_DATABASE_PREFERENCE if name != "PROSITE"] + + +def test_reorder_is_case_insensitive_and_strips_whitespace(): + result = reorder_database_preference([" pfam ", "Smart"]) + assert result[:2] == ["PFAM", "SMART"] + + +def test_reorder_empty_priority_returns_default_order(): + assert reorder_database_preference([]) == DEFAULT_DATABASE_PREFERENCE + + +def test_reorder_deduplicates_repeated_names(): + result = reorder_database_preference(["PFAM", "PFAM"]) + assert result.count("PFAM") == 1 + assert len(result) == len(DEFAULT_DATABASE_PREFERENCE) + + +def test_reorder_raises_on_unknown_database_name(): + with pytest.raises(ValueError, match="Unknown database name"): + reorder_database_preference(["NOTAREALDB"]) diff --git a/tests/test_famsa_aligner.py b/tests/test_famsa_aligner.py new file mode 100644 index 0000000..e2a4737 --- /dev/null +++ b/tests/test_famsa_aligner.py @@ -0,0 +1,27 @@ +from candy.alignment import get_alignment_tool +from candy.alignment.famsa import FamsaAligner + + +def test_get_alignment_tool_returns_famsa_by_name(): + assert get_alignment_tool("famsa").name == "famsa" + + +def test_famsa_aligner_produces_equal_length_aligned_records(tmp_path): + input_fasta = tmp_path / "in.fasta" + output_fasta = tmp_path / "out.fasta" + input_fasta.write_text( + ">a\nMKVLAMKVLA\n" + ">b\nMKVLPMKVLAEXTRA\n" + ">c\nMKVAAMKVLA\n" + ) + + FamsaAligner().align(input_fasta, output_fasta) + + from Bio import SeqIO + + records = list(SeqIO.parse(output_fasta, "fasta")) + ids = {r.id for r in records} + lengths = {len(r.seq) for r in records} + + assert ids == {"a", "b", "c"} + assert len(lengths) == 1 # all aligned records share one length diff --git a/tests/test_mmseqs2_download.py b/tests/test_mmseqs2_download.py new file mode 100644 index 0000000..b4440e0 --- /dev/null +++ b/tests/test_mmseqs2_download.py @@ -0,0 +1,96 @@ +import os +from pathlib import Path +from unittest.mock import patch + +import pytest + +from candy.external_tools import MissingDependencyError +from candy.clustering import mmseqs2_download as m + + +def test_prefers_binary_already_on_path(): + with patch("candy.clustering.mmseqs2_download.find_binary", return_value="/usr/bin/mmseqs"): + assert m.resolve_mmseqs2_binary() == "/usr/bin/mmseqs" + + +def test_uses_cached_binary_if_present(tmp_path): + cached = tmp_path / "extracted" / "mmseqs" / "bin" / m._binary_name() + cached.parent.mkdir(parents=True) + cached.write_text("fake binary") + + with patch("candy.clustering.mmseqs2_download.find_binary", return_value=None), patch( + "candy.clustering.mmseqs2_download._cache_dir", return_value=tmp_path + ): + assert m.resolve_mmseqs2_binary() == str(cached) + + +def test_raises_when_no_auto_download_env_set_and_nothing_available(tmp_path): + with patch("candy.clustering.mmseqs2_download.find_binary", return_value=None), patch( + "candy.clustering.mmseqs2_download._cache_dir", return_value=tmp_path + ), patch.dict(os.environ, {m.NO_AUTO_DOWNLOAD_ENV: "1"}): + with pytest.raises(MissingDependencyError, match="automatic download is disabled"): + m.resolve_mmseqs2_binary() + + +def test_raises_clear_error_for_unknown_platform(tmp_path): + with patch("candy.clustering.mmseqs2_download.find_binary", return_value=None), patch( + "candy.clustering.mmseqs2_download._cache_dir", return_value=tmp_path + ), patch("candy.clustering.mmseqs2_download.platform.system", return_value="Plan9"), patch( + "candy.clustering.mmseqs2_download.platform.machine", return_value="mips" + ), patch.dict(os.environ, {}, clear=False): + os.environ.pop(m.NO_AUTO_DOWNLOAD_ENV, None) + with pytest.raises(MissingDependencyError, match="No known MMseqs2 static build"): + m.resolve_mmseqs2_binary() + + +def test_download_failure_raises_missing_dependency_error(tmp_path): + import requests + + with patch("candy.clustering.mmseqs2_download.find_binary", return_value=None), patch( + "candy.clustering.mmseqs2_download._cache_dir", return_value=tmp_path + ), patch("candy.clustering.mmseqs2_download.platform.system", return_value="Linux"), patch( + "candy.clustering.mmseqs2_download.platform.machine", return_value="x86_64" + ), patch("candy.clustering.mmseqs2_download.requests.get", side_effect=requests.ConnectionError("no network")): + os.environ.pop(m.NO_AUTO_DOWNLOAD_ENV, None) + with pytest.raises(MissingDependencyError, match="Failed to download MMseqs2"): + m.resolve_mmseqs2_binary() + + +def test_asset_selection_covers_major_platforms(): + assert m._ASSET_BY_PLATFORM[("Linux", "x86_64")] == "mmseqs-linux-sse41.tar.gz" + assert m._ASSET_BY_PLATFORM[("Darwin", "arm64")] == "mmseqs-osx-universal.tar.gz" + assert m._ASSET_BY_PLATFORM[("Windows", "AMD64")] == "mmseqs-win64.zip" + + +@pytest.mark.integration +def test_real_download_and_run_mmseqs2(tmp_path): + """Exercises the real network download + extraction against the pinned release. + + Marked integration since it needs network access; skipped by default runs + that filter out that marker, but was run manually during development to + confirm the auto-download path genuinely works end-to-end. + + On Windows this deliberately does NOT invoke the resolved mmseqs.bat + directly: .bat always attempts its one-time busybox self-install (even + for a plain "version" call), which can pop a UAC prompt -- unsuitable + for an unattended/CI test run. Instead it locates the underlying + mmseqs.exe inside the same extracted tree and runs that directly, which + validates the download+extraction worked without touching the + self-install path. On Linux/macOS, resolve_mmseqs2_binary() already + returns the plain binary, so it's exercised directly either way. + """ + import platform + import subprocess + + with patch("candy.clustering.mmseqs2_download.find_binary", return_value=None), patch( + "candy.clustering.mmseqs2_download._cache_dir", return_value=tmp_path + ): + os.environ.pop(m.NO_AUTO_DOWNLOAD_ENV, None) + binary = m.resolve_mmseqs2_binary() + + if platform.system() == "Windows": + binary = str(next(Path(binary).parent.glob("bin/mmseqs.exe"))) + + result = subprocess.run([binary, "version"], capture_output=True, text=True) + assert result.returncode == 0 + assert result.stdout.strip() diff --git a/tests/test_veryfasttree_builder.py b/tests/test_veryfasttree_builder.py new file mode 100644 index 0000000..ba71555 --- /dev/null +++ b/tests/test_veryfasttree_builder.py @@ -0,0 +1,41 @@ +import random + +from candy.phylogenetics import get_tree_builder +from candy.phylogenetics.veryfasttree_builder import VeryFastTreeBuilder + + +def test_get_tree_builder_returns_veryfasttree_by_name(): + assert get_tree_builder("veryfasttree").name == "veryfasttree" + + +def _random_alignment(n_sequences=12, length=200, seed=42): + random.seed(seed) + amino_acids = "ACDEFGHIKLMNPQRSTVWY" + base = "".join(random.choice(amino_acids) for _ in range(length)) + records = [] + for i in range(n_sequences): + seq = list(base) + for _ in range(10): + pos = random.randrange(len(seq)) + seq[pos] = random.choice(amino_acids) + records.append((f"seq{i}", "".join(seq))) + return "".join(f">{name}\n{seq}\n" for name, seq in records) + + +def test_veryfasttree_builder_produces_a_tree_file(tmp_path): + # Regression test: veryfasttree's default ext="AUTO" was found to crash + # (access violation) with a realistically-sized alignment on at least one + # Windows/AVX2 machine; VeryFastTreeBuilder must always pass an explicit + # ext= to avoid that path. A tiny (~4 sequence) alignment reproduced the + # same crash, so this uses a more realistic size to keep the regression + # meaningful. + alignment_fasta = tmp_path / "aligned.fasta" + alignment_fasta.write_text(_random_alignment()) + output_newick = tmp_path / "tree.nwk" + + VeryFastTreeBuilder().build_tree(alignment_fasta, output_newick) + + tree_text = output_newick.read_text() + assert tree_text.strip().endswith(";") + for i in range(12): + assert f"seq{i}" in tree_text From 7e57477851ea845d00fca37f63ecf10cb908600b Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Thu, 6 Aug 2026 21:34:00 +0200 Subject: [PATCH 3/7] Fix real-world MMseqs2 failures on Windows found during manual testing Manually running CANDy end-to-end against a live CAZy family (GH173) surfaced two distinct MMseqs2 failures on Windows, both fixed with regression tests and verified against a real 155-sequence clustering run: - Clustering silently produced no output when the project path contained a space (e.g. under "...\OneDrive - Org\..."): MMseqs2's Windows build routes its clustering workflow through a bundled Cygwin/busybox POSIX shell layer, which mishandles such paths. Clustering now always runs from a space-free temp directory with relative filenames. - easy-cluster's internal FASTA-conversion step (result2flat) segfaulted through that same Cygwin process-spawning layer on real data, even though every preceding clustering step succeeded. Replaced the easy-cluster convenience workflow with its constituent steps run individually (createdb -> cluster -> result2repseq -> convert2fasta), using the shell-capable mmseqs.bat only for the one step that actually needs a POSIX shell (cluster's cascaded-clustering algorithm) and the plain compiled binary directly for the rest. - mmseqs.bat always exits 0 regardless of internal failure, so success is now verified by checking the expected output file actually exists rather than trusting the return code. Also broadens NCBI/CAZy fetch retry logic beyond urllib.error.HTTPError to cover the wider class of transient connection failures (a real run hit http.client.IncompleteRead mid-fetch on a large batched Entrez request, which the narrower except-clause didn't catch). Documents the interactive domain-name curation flow and how to switch to Gemini-based curation with your own API key, and prints a citation reminder (Windels et al., PLoS One 2024) at the start of every pipeline run. Co-Authored-By: Claude Sonnet 5 --- README.md | 15 ++++++ src/candy/cazy.py | 62 +++++++++++++++++----- src/candy/clustering/mmseqs2.py | 81 ++++++++++++++++++++++++----- src/candy/pipeline.py | 8 +++ tests/test_cazy.py | 73 +++++++++++++++++++++++++- tests/test_clustering.py | 91 +++++++++++++++++++++++++++++---- 6 files changed, 292 insertions(+), 38 deletions(-) diff --git a/README.md b/README.md index 040161c..6cf84d2 100644 --- a/README.md +++ b/README.md @@ -58,6 +58,21 @@ candy my_sequences.fasta --tree `--email` falls back to the `CANDY_EMAIL` environment variable, then an interactive prompt, so `export CANDY_EMAIL=you@example.com` once and just run `candy GH173` from then on. Run `candy --help` for the full list of options. +### Domain-name curation + +Partway through a run, CANDy needs to decide which raw InterPro domain names (often several near-duplicates from different member databases) should be grouped under one umbrella name. By default (`--curation-backend manual`) it asks you interactively: it prints a numbered list of the domain names still to curate, then prompts twice -- + +1. `Domain name:` -- type the umbrella name you want to use (e.g. `Catalytic domain`) +2. `Includes:` -- type the comma-separated numbers of the domains that belong under it (e.g. `0,2`) + +It repeats this until every domain is grouped; type `STOP` at the `Domain name:` prompt at any point to leave all remaining domains as their own individual groups. + +To skip this entirely, use Gemini to curate automatically instead: + +1. Install the extra: `pip install -e ".[gemini]"` +2. Get a free API key at [aistudio.google.com/app/api-keys](https://aistudio.google.com/app/api-keys) +3. Run with `--curation-backend gemini --curation-api-key YOUR_KEY`, or set it once via `$env:GOOGLE_API_KEY="YOUR_KEY"` (PowerShell) / `export GOOGLE_API_KEY=YOUR_KEY` (bash) and just pass `--curation-backend gemini` + ### Python API ```python diff --git a/src/candy/cazy.py b/src/candy/cazy.py index 54e4058..a96fde0 100644 --- a/src/candy/cazy.py +++ b/src/candy/cazy.py @@ -10,11 +10,12 @@ from __future__ import annotations +import http.client import io import logging import time from collections.abc import Iterable, Sequence -from urllib.error import HTTPError +from urllib.error import URLError import pandas as pd import requests @@ -39,19 +40,36 @@ _ENTREZ_MAX_RETRIES = 10 _ENTREZ_RETRY_DELAY = 3.0 +# urllib.error.HTTPError is a URLError subclass, so this also covers plain +# HTTP error statuses. http.client.HTTPException covers things like +# IncompleteRead (the connection dropping mid-response, which NCBI's Entrez +# endpoints hit occasionally on large batched fetches) and BadStatusLine; +# ConnectionError/TimeoutError cover socket-level drops and timeouts. +_RETRYABLE_NETWORK_ERRORS = (URLError, http.client.HTTPException, ConnectionError, TimeoutError) + def _batched(items: Sequence[str], size: int) -> Iterable[Sequence[str]]: for i in range(0, len(items), size): yield items[i : i + size] -def fetch_family_page(family: str) -> str: +def fetch_family_page(family: str, max_retries: int = 5, retry_delay: float = 5.0) -> str: """Download the raw CAZy family data page (tab-separated GenBank ID listing).""" url = CAZY_FAMILY_URL.format(family=family) logger.info("Retrieving data from %s", url) - response = requests.get(url, timeout=60) - response.raise_for_status() - return response.text + + attempts = 0 + while True: + try: + response = requests.get(url, timeout=60) + response.raise_for_status() + return response.text + except requests.RequestException as exc: + attempts += 1 + if attempts >= max_retries: + raise + logger.info("Network error retrieving CAZy data (%s); retrying in %.0fs.", exc, retry_delay) + time.sleep(retry_delay) def parse_family_table(raw_text: str) -> tuple[pd.DataFrame, dict[str, str]]: @@ -90,11 +108,14 @@ def fetch_sequences_fasta(ids: Sequence[str], email: str, batch_size: int = _ENT with Entrez.efetch(db="protein", id=",".join(batch), rettype="fasta", retmode="text") as handle: chunks.append(handle.read()) break - except HTTPError: + except _RETRYABLE_NETWORK_ERRORS as exc: attempts += 1 if attempts >= _ENTREZ_MAX_RETRIES: - logger.warning("Giving up on a batch of %d IDs after repeated HTTP errors.", len(batch)) + logger.warning( + "Giving up on a batch of %d IDs after repeated network errors (%s).", len(batch), exc + ) break + logger.info("Network error fetching a sequence batch (%s); retrying in %.0fs.", exc, _ENTREZ_RETRY_DELAY) time.sleep(_ENTREZ_RETRY_DELAY) logger.info("Fetched %d/%d sequence batches.", len(chunks), total_batches) @@ -134,11 +155,11 @@ def fetch_characterized_page(family: str, max_retries: int = 10, retry_delay: fl while True: try: return pd.read_html(url) - except HTTPError: + except _RETRYABLE_NETWORK_ERRORS as exc: attempts += 1 if attempts >= max_retries: raise - logger.info("Too many requests to CAZy; retrying in %.0fs.", retry_delay) + logger.info("Network error retrieving CAZy data (%s); retrying in %.0fs.", exc, retry_delay) time.sleep(retry_delay) @@ -191,11 +212,24 @@ def fetch_characterized_sequences( output = "" for genbank_id in filtered_ids: - try: - with Entrez.efetch(db="protein", id=genbank_id, rettype="fasta", retmode="text") as handle: - text = handle.read() - except HTTPError: - logger.warning("HTTP error fetching characterized sequence %s.", genbank_id) + text = None + attempts = 0 + while attempts < _ENTREZ_MAX_RETRIES: + try: + with Entrez.efetch(db="protein", id=genbank_id, rettype="fasta", retmode="text") as handle: + text = handle.read() + break + except _RETRYABLE_NETWORK_ERRORS as exc: + attempts += 1 + if attempts >= _ENTREZ_MAX_RETRIES: + logger.warning( + "Giving up on characterized sequence %s after repeated network errors (%s).", + genbank_id, exc, + ) + break + time.sleep(_ENTREZ_RETRY_DELAY) + + if text is None: continue from Bio import SeqIO diff --git a/src/candy/clustering/mmseqs2.py b/src/candy/clustering/mmseqs2.py index 1b60470..38b4271 100644 --- a/src/candy/clustering/mmseqs2.py +++ b/src/candy/clustering/mmseqs2.py @@ -8,7 +8,7 @@ from candy.clustering.mmseqs2_download import resolve_mmseqs2_binary from candy.config import ClusteringConfig -from candy.external_tools import run_tool +from candy.external_tools import ExternalToolError, run_tool logger = logging.getLogger(__name__) @@ -21,13 +21,42 @@ class Mmseqs2Clusterer: :mod:`candy.clustering.mmseqs2_download`). This is the default clustering backend precisely because that auto-download is possible for MMseqs2 but not for CD-HIT (source-only, no Windows build). + + .. note:: + This deliberately does *not* use MMseqs2's ``easy-cluster`` convenience + workflow, and instead runs its constituent steps individually + (``createdb`` -> ``cluster`` -> ``result2repseq`` -> ``convert2fasta``). + Three Windows-specific issues were found running the real thing: + + 1. ``mmseqs.bat`` always exits 0 regardless of whether the underlying + command actually succeeded (its final line is a hardcoded + ``exit /b 0``), so a non-zero return code can't be used to detect + failure there -- worked around by always checking that the + expected output file was actually produced. + 2. MMseqs2's clustering workflow runs through a bundled Cygwin/busybox + POSIX-shell layer internally (see ``mmseqs2_download``), which + mishandles Windows paths containing spaces (e.g. a project + directory under "...\\OneDrive - Org\\..."). Avoided by always + running from a space-free temporary directory with relative + filenames, regardless of where the real input/output paths live. + 3. ``easy-cluster``'s internal script runs its final DB-to-FASTA + conversion step (``result2flat``) through that same Cygwin process- + spawning layer, which was observed to segfault on a real (155 + sequence) clustering run despite every preceding step completing + successfully. The equivalent public step, ``convert2fasta``, runs + reliably when invoked directly (bypassing the shell entirely) -- + so only the one step that genuinely needs a POSIX shell (the + cascaded-clustering algorithm inside ``cluster`` itself) goes + through ``mmseqs.bat``; everything else uses the plain compiled + binary directly. """ name = "mmseqs2" def cluster(self, input_fasta: Path, output_fasta: Path, config: ClusteringConfig) -> Path: - binary = resolve_mmseqs2_binary() - if platform.system() == "Windows" and binary.lower().endswith(".bat"): + shell_binary = resolve_mmseqs2_binary() + direct_binary = self._direct_binary(shell_binary) + if platform.system() == "Windows" and shell_binary.lower().endswith(".bat"): logger.info( "Clustering on Windows uses MMseqs2's mmseqs.bat wrapper, which may ask for " "administrator permission once (to install a small POSIX-shell helper it " @@ -35,24 +64,50 @@ def cluster(self, input_fasta: Path, output_fasta: Path, config: ClusteringConfi ) min_seq_id = config.identity_cutoff / 100 min_coverage = config.mmseqs_min_coverage / 100 + cov_mode = str(config.mmseqs_cov_mode) logger.info("Clustering with MMseqs2 at %.0f%% identity.", config.identity_cutoff) with tempfile.TemporaryDirectory(prefix="candy_mmseqs2_") as tmp_dir: tmp_dir = Path(tmp_dir) - output_prefix = tmp_dir / "cluster" - run_tool( + shutil.copyfile(input_fasta, tmp_dir / "input.fasta") + + run_tool([direct_binary, "createdb", "input.fasta", "db"], cwd=tmp_dir) + result = run_tool( [ - binary, - "easy-cluster", - str(input_fasta), - str(output_prefix), - str(tmp_dir / "tmp"), + shell_binary, + "cluster", "db", "db_clu", "tmp", "--min-seq-id", f"{min_seq_id}", "-c", f"{min_coverage}", - "--cov-mode", str(config.mmseqs_cov_mode), - ] + "--cov-mode", cov_mode, + ], + cwd=tmp_dir, ) - representative_sequences = Path(f"{output_prefix}_rep_seq.fasta") + run_tool([direct_binary, "result2repseq", "db", "db_clu", "db_clu_rep"], cwd=tmp_dir) + run_tool([direct_binary, "convert2fasta", "db_clu_rep", "cluster_rep_seq.fasta"], cwd=tmp_dir) + + representative_sequences = tmp_dir / "cluster_rep_seq.fasta" + if not representative_sequences.exists(): + raise ExternalToolError( + "MMseqs2 clustering reported success but did not produce the expected " + f"output ({representative_sequences.name}). Captured output from the " + f"clustering step:\nstdout:\n{result.stdout}\nstderr:\n{result.stderr}" + ) shutil.copyfile(representative_sequences, output_fasta) return output_fasta + + @staticmethod + def _direct_binary(shell_binary: str) -> str: + """Resolve the plain compiled ``mmseqs`` binary, bypassing any shell wrapper. + + On Windows, ``resolve_mmseqs2_binary()`` returns ``mmseqs.bat`` (one + directory above a ``bin/`` folder containing ``mmseqs.exe``); steps + that don't need a POSIX shell run against that ``.exe`` directly + instead. Everywhere else, the shell binary already *is* the plain + binary, so this is a no-op. + """ + if platform.system() == "Windows" and shell_binary.lower().endswith(".bat"): + candidate = Path(shell_binary).parent / "bin" / "mmseqs.exe" + if candidate.is_file(): + return str(candidate) + return shell_binary diff --git a/src/candy/pipeline.py b/src/candy/pipeline.py index 8b7da00..78f8bb8 100644 --- a/src/candy/pipeline.py +++ b/src/candy/pipeline.py @@ -32,6 +32,13 @@ logger = logging.getLogger(__name__) +CITATION = ( + "If you use CANDy and publish the results, please cite:\n" + "Windels A, Franceus J, Pleiss J, Desmet T. CANDy: Automated analysis of domain architectures " + "in carbohydrate-active enzymes. PLoS One. 2024 Jul 11;19(7):e0306410. " + "doi: 10.1371/journal.pone.0306410. PMID: 38990885; PMCID: PMC11238990." +) + @dataclass class PipelineResult: @@ -83,6 +90,7 @@ def _unique_jobname(output_dir: Path, jobname: str) -> str: def run_pipeline(config: PipelineConfig) -> PipelineResult: + print(f"\n{CITATION}\n") config.jobname = _unique_jobname(config.output_dir, config.jobname) jobname_dir = config.output_dir / config.jobname jobname_dir.mkdir(parents=True, exist_ok=True) diff --git a/tests/test_cazy.py b/tests/test_cazy.py index 7e8761e..b278dea 100644 --- a/tests/test_cazy.py +++ b/tests/test_cazy.py @@ -1,6 +1,15 @@ +import http.client +from unittest.mock import patch + import pandas as pd -from candy.cazy import extract_characterized_ids, extract_ncbi_ids, parse_family_table +from candy.cazy import ( + extract_characterized_ids, + extract_ncbi_ids, + fetch_characterized_sequences, + fetch_sequences_fasta, + parse_family_table, +) from candy.config import Taxonomy SAMPLE_FAMILY_TABLE = ( @@ -40,3 +49,65 @@ def test_extract_characterized_ids_groups_by_taxonomy_header(): assert ids == ["P12345.1", "Q99999.2"] assert taxonomy_dict == {"P12345.1": "B", "Q99999.2": "A"} assert activity_dict == {"P12345.1": "1.2.4.-", "Q99999.2": "3.2.1.-"} + + +class FakeHandle: + """Stand-in for the object Bio.Entrez.efetch returns (used as a context manager).""" + + def __init__(self, text): + self._text = text + + def __enter__(self): + return self + + def __exit__(self, *exc): + return False + + def read(self): + return self._text + + +def test_fetch_sequences_fasta_retries_past_incomplete_read(): + # Regression test: a real run hit http.client.IncompleteRead mid-fetch + # (NCBI's connection dropped during a large chunked response) and the + # original except-HTTPError-only clause let it crash the whole pipeline + # instead of retrying like other transient network errors. + calls = {"n": 0} + + def fake_efetch(db, id, rettype, retmode): + calls["n"] += 1 + if calls["n"] == 1: + raise http.client.IncompleteRead(b"partial") + return FakeHandle(">P1\nMKV\n") + + with patch("candy.cazy.Entrez.efetch", side_effect=fake_efetch), patch("candy.cazy.time.sleep"): + result = fetch_sequences_fasta(["P1"], "test@example.com") + + assert result == ">P1\nMKV\n" + assert calls["n"] == 2 + + +def test_fetch_sequences_fasta_gives_up_after_max_retries(): + def fake_efetch(db, id, rettype, retmode): + raise http.client.IncompleteRead(b"") + + with patch("candy.cazy.Entrez.efetch", side_effect=fake_efetch), patch("candy.cazy.time.sleep"): + result = fetch_sequences_fasta(["P1"], "test@example.com") + + assert result == "" # gives up gracefully rather than raising + + +def test_fetch_characterized_sequences_retries_past_connection_error(): + calls = {"n": 0} + + def fake_efetch(db, id, rettype, retmode): + calls["n"] += 1 + if calls["n"] == 1: + raise ConnectionError("connection reset") + return FakeHandle(">P1.1 enzyme [Escherichia coli]\nMKV\n") + + with patch("candy.cazy.Entrez.efetch", side_effect=fake_efetch), patch("candy.cazy.time.sleep"): + result = fetch_characterized_sequences(["P1.1"], {"P1.1": "B"}, Taxonomy.ALL, "test@example.com") + + assert "MKV" in result + assert calls["n"] == 2 diff --git a/tests/test_clustering.py b/tests/test_clustering.py index 7449204..d520600 100644 --- a/tests/test_clustering.py +++ b/tests/test_clustering.py @@ -1,3 +1,5 @@ +from pathlib import Path +from types import SimpleNamespace from unittest.mock import patch from candy.config import ClusteringConfig, ClusteringSoftware @@ -30,26 +32,95 @@ def test_cdhit_builds_expected_command(tmp_path): assert "-o" in args and args[args.index("-o") + 1] == str(output_fasta) -def test_mmseqs2_builds_expected_command_and_copies_result(tmp_path): +def test_mmseqs2_runs_individual_steps_and_copies_result(tmp_path): from candy.clustering.mmseqs2 import Mmseqs2Clusterer config = ClusteringConfig(identity_cutoff=85, mmseqs_min_coverage=90, mmseqs_cov_mode=0) input_fasta = tmp_path / "in.fasta" + input_fasta.write_text(">seq\nMKV\n") output_fasta = tmp_path / "out.fasta" + calls = [] + def fake_run_tool(args, **kwargs): - # mimic mmseqs2 easy-cluster writing its representative-sequence output - prefix = args[3] - with open(f"{prefix}_rep_seq.fasta", "w") as f: - f.write(">rep\nMKV\n") + calls.append(args) + cwd = Path(kwargs["cwd"]) + # mimic convert2fasta writing the final representative-sequence FASTA + if args[1] == "convert2fasta": + (cwd / "cluster_rep_seq.fasta").write_text(">rep\nMKV\n") + return SimpleNamespace(stdout="", stderr="", returncode=0) with patch("candy.clustering.mmseqs2.resolve_mmseqs2_binary", return_value="/usr/bin/mmseqs"), patch( "candy.clustering.mmseqs2.run_tool", side_effect=fake_run_tool - ) as mock_run: + ): Mmseqs2Clusterer().cluster(input_fasta, output_fasta, config) - args = mock_run.call_args[0][0] - assert args[0] == "/usr/bin/mmseqs" - assert args[1] == "easy-cluster" - assert "--min-seq-id" in args and args[args.index("--min-seq-id") + 1] == "0.85" + subcommands = [call[1] for call in calls] + assert subcommands == ["createdb", "cluster", "result2repseq", "convert2fasta"] + # only the shell-dependent `cluster` step uses the shell-capable binary; + # everything else runs the plain binary directly (no shell involved) + assert calls[0][0] == "/usr/bin/mmseqs" # createdb + assert calls[1][0] == "/usr/bin/mmseqs" # cluster (shell binary, same path here since no .bat) + assert calls[0][2] == "input.fasta" # relative filename, not an absolute path with spaces + + cluster_call = calls[1] + assert "--min-seq-id" in cluster_call and cluster_call[cluster_call.index("--min-seq-id") + 1] == "0.85" assert output_fasta.read_text() == ">rep\nMKV\n" + + +def test_mmseqs2_uses_plain_exe_directly_for_non_shell_steps_on_windows(tmp_path): + """The cluster step goes through mmseqs.bat; createdb/result2repseq/convert2fasta bypass it.""" + from candy.clustering.mmseqs2 import Mmseqs2Clusterer + + bat_path = tmp_path / "mmseqs" / "mmseqs.bat" + exe_path = tmp_path / "mmseqs" / "bin" / "mmseqs.exe" + exe_path.parent.mkdir(parents=True) + exe_path.write_text("fake exe") + bat_path.write_text("fake bat") + + config = ClusteringConfig(identity_cutoff=85) + input_fasta = tmp_path / "in.fasta" + input_fasta.write_text(">seq\nMKV\n") + output_fasta = tmp_path / "out.fasta" + + calls = [] + + def fake_run_tool(args, **kwargs): + calls.append(args) + cwd = Path(kwargs["cwd"]) + if args[1] == "convert2fasta": + (cwd / "cluster_rep_seq.fasta").write_text(">rep\nMKV\n") + return SimpleNamespace(stdout="", stderr="", returncode=0) + + with patch("candy.clustering.mmseqs2.resolve_mmseqs2_binary", return_value=str(bat_path)), patch( + "candy.clustering.mmseqs2.platform.system", return_value="Windows" + ), patch("candy.clustering.mmseqs2.run_tool", side_effect=fake_run_tool): + Mmseqs2Clusterer().cluster(input_fasta, output_fasta, config) + + binaries_used = [call[0] for call in calls] + assert binaries_used == [str(exe_path), str(bat_path), str(exe_path), str(exe_path)] + + +def test_mmseqs2_raises_clear_error_when_output_missing_despite_success_exit_code(tmp_path): + """Regression test: mmseqs.bat always exits 0 on Windows even on internal failure.""" + from candy.clustering.mmseqs2 import Mmseqs2Clusterer + from candy.external_tools import ExternalToolError + + config = ClusteringConfig(identity_cutoff=85) + input_fasta = tmp_path / "in.fasta" + input_fasta.write_text(">seq\nMKV\n") + output_fasta = tmp_path / "out.fasta" + + def fake_run_tool(args, **kwargs): + # everything "succeeded" (no exception) but never actually produced output + return SimpleNamespace(stdout="some diagnostic output", stderr="", returncode=0) + + with patch("candy.clustering.mmseqs2.resolve_mmseqs2_binary", return_value="/usr/bin/mmseqs"), patch( + "candy.clustering.mmseqs2.run_tool", side_effect=fake_run_tool + ): + try: + Mmseqs2Clusterer().cluster(input_fasta, output_fasta, config) + assert False, "expected ExternalToolError" + except ExternalToolError as exc: + assert "did not produce the expected output" in str(exc) + assert "some diagnostic output" in str(exc) From e44b836989062b9354f38a47768dbb3111b5b4ca Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Fri, 7 Aug 2026 09:59:22 +0200 Subject: [PATCH 4/7] Fix Gemini curation default model and add --curation-model override The hardcoded gemini-2.0-flash default lost its free-tier quota, and the next hardcoded guess (gemini-2.5-flash) turned out to be restricted for new users too -- discovered via a real run, not just docs. Switch to Google's gemini-flash-latest alias, which tracks whatever their current free-tier flash model is, and add a --curation-model flag so a future rename/deprecation doesn't need a code change. Co-Authored-By: Claude Sonnet 5 --- src/candy/cli.py | 8 +++++++- src/candy/config.py | 1 + src/candy/curation/gemini.py | 2 +- src/candy/pipeline.py | 7 ++++++- 4 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/candy/cli.py b/src/candy/cli.py index b1d3e2a..078218e 100644 --- a/src/candy/cli.py +++ b/src/candy/cli.py @@ -87,6 +87,12 @@ def main( curation_api_key: Optional[str] = typer.Option( None, envvar="GOOGLE_API_KEY", help="API key for the curation backend, if it needs one." ), + curation_model: Optional[str] = typer.Option( + None, + help="Model name for the curation backend, if it needs one (e.g. a Gemini model name). " + "Defaults to the backend's own default. Useful to override without a code change if a " + "provider deprecates/renames their default model.", + ), verbose: bool = typer.Option(False, "--verbose", "-v", help="Enable debug logging."), ) -> None: """Run the CANDy pipeline end-to-end. @@ -149,7 +155,7 @@ def main( overlap_percentage=overlap_percentage, database_preference=database_preference, ), - curation=CurationConfig(backend=curation_backend, api_key=curation_api_key), + curation=CurationConfig(backend=curation_backend, api_key=curation_api_key, model=curation_model), blast_identity_threshold=blast_identity, build_tree=build_tree, ) diff --git a/src/candy/config.py b/src/candy/config.py index b0cdfae..6dae5b3 100644 --- a/src/candy/config.py +++ b/src/candy/config.py @@ -149,6 +149,7 @@ class CurationConfig: backend: str = "manual" # "manual" | "gemini" | any registered CurationBackend name api_key: str | None = None + model: str | None = None # backend-specific model name; None uses the backend's default @dataclass diff --git a/src/candy/curation/gemini.py b/src/candy/curation/gemini.py index f9a2592..2865e15 100644 --- a/src/candy/curation/gemini.py +++ b/src/candy/curation/gemini.py @@ -28,7 +28,7 @@ class GeminiCurationBackend: name = "gemini" - def __init__(self, api_key: str | None = None, model: str = "gemini-2.0-flash") -> None: + def __init__(self, api_key: str | None = None, model: str = "gemini-flash-latest") -> None: self.api_key = api_key or os.environ.get("GOOGLE_API_KEY") or os.environ.get("GEMINI_API_KEY") if not self.api_key: raise ValueError( diff --git a/src/candy/pipeline.py b/src/candy/pipeline.py index 78f8bb8..a4e94b4 100644 --- a/src/candy/pipeline.py +++ b/src/candy/pipeline.py @@ -126,7 +126,12 @@ def run_pipeline(config: PipelineConfig) -> PipelineResult: # --- Stage 8: curate synonymous domain names --- raw_domain_names = sorted({name for hits in domain_architecture.values() for _, name in hits}) - curation_kwargs = {} if config.curation.backend == "manual" else {"api_key": config.curation.api_key} + if config.curation.backend == "manual": + curation_kwargs = {} + else: + curation_kwargs = {"api_key": config.curation.api_key} + if config.curation.model: + curation_kwargs["model"] = config.curation.model curation_backend = get_curation_backend(config.curation.backend, **curation_kwargs) family_label = config.input.family if is_cazy_query else None curated_domains = curation_backend.curate(raw_domain_names, family=family_label) From c8dccab5db1d8b729902ccb7b56966c8fe3c61f1 Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Fri, 7 Aug 2026 11:14:54 +0200 Subject: [PATCH 5/7] Prepare v3.0.0 for PyPI release - Add top-level LICENSE (previously only embedded in SECURITY.md, which GitHub/PyPI license detection doesn't look at) - Add project.urls and bump version off the 3.0.0.dev0 pre-release marker - Add CI: run pytest on every push/PR, and publish to PyPI on GitHub Release via trusted publishing (OIDC, no stored API token) - Archive the original notebook under archive/, update README to reflect the package as the primary distribution rather than a WIP branch Co-Authored-By: Claude Sonnet 5 --- .github/workflows/publish.yml | 33 ++++++++++++++++++++ .github/workflows/tests.yml | 22 +++++++++++++ LICENSE | 21 +++++++++++++ README.md | 11 ++----- CANDy v2.0.ipynb => archive/CANDy v2.0.ipynb | 0 pyproject.toml | 7 ++++- 6 files changed, 85 insertions(+), 9 deletions(-) create mode 100644 .github/workflows/publish.yml create mode 100644 .github/workflows/tests.yml create mode 100644 LICENSE rename CANDy v2.0.ipynb => archive/CANDy v2.0.ipynb (100%) diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml new file mode 100644 index 0000000..4ac49d2 --- /dev/null +++ b/.github/workflows/publish.yml @@ -0,0 +1,33 @@ +name: Publish to PyPI + +on: + release: + types: [published] + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: "3.12" + - run: pip install build + - run: python -m build + - uses: actions/upload-artifact@v4 + with: + name: dist + path: dist/ + + publish: + needs: build + runs-on: ubuntu-latest + environment: pypi + permissions: + id-token: write # required for PyPI trusted publishing (OIDC) + steps: + - uses: actions/download-artifact@v4 + with: + name: dist + path: dist/ + - uses: pypa/gh-action-pypi-publish@release/v1 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..7be8e14 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,22 @@ +name: Tests + +on: + push: + branches: [main] + pull_request: + +jobs: + test: + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, windows-latest] + python-version: ["3.10", "3.11", "3.12"] + runs-on: ${{ matrix.os }} + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + - run: pip install -e ".[dev]" + - run: pytest -q diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..c9fb56c --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2023 Alex Windels + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md index 6cf84d2..318013a 100644 --- a/README.md +++ b/README.md @@ -2,20 +2,15 @@ CANDy is a fast, FAIR and seamless protein domain analysis tool for any [CAZy](http://www.cazy.org/) family. -**This branch is a work-in-progress rewrite of CANDy as an installable Python package** (CLI + Python API), replacing the original Google Colab / Jupyter Notebook implementation. The previous notebook (`CANDy v2.0.ipynb`) is still in this repository for reference, and the original is on [Google Colab](https://colab.research.google.com/drive/1ipRAwMFMDRGUinPDk2bwu1cg8fE2WY8Q?usp=sharing). +CANDy is available as an installable Python package (CLI + Python API), replacing the original Google Colab / Jupyter Notebook implementation. The previous notebook (`archive/CANDy v2.0.ipynb`) is kept in this repository for reference, and the original is on [Google Colab](https://colab.research.google.com/drive/1ipRAwMFMDRGUinPDk2bwu1cg8fE2WY8Q?usp=sharing). ## Installation -Not yet published to PyPI -- install directly from this repository for now: - ```bash -git clone https://github.com/PyEED/CANDy.git -cd CANDy -git checkout package-upgrade -pip install -e . +pip install candy-cazyme ``` -(Once published, this will just be `pip install candy-cazyme`.) That's it for most users -- CANDy's default toolchain is fully bundled: +That's it for most users -- CANDy's default toolchain is fully bundled: - **Clustering**: [MMseqs2](https://github.com/soedinglab/MMseqs2) -- auto-downloaded and cached on first use (no conda needed). On Linux/macOS this just works. On **Windows**, MMseqs2's clustering workflows internally need a POSIX shell; the official Windows build handles this itself by installing a small helper (`busybox`) the first time it runs, which may ask for administrator permission **once** -- never again after that. (This mirrors upstream: MMseqs2's own docs list WSL as the recommended Windows path and this static build as the fallback for anyone who can't use WSL.) - **MSA**: [FAMSA](https://github.com/refresh-bio/FAMSA) via [`pyfamsa`](https://github.com/althonos/pyfamsa) -- a real pip dependency, runs in-process, no download needed. diff --git a/CANDy v2.0.ipynb b/archive/CANDy v2.0.ipynb similarity index 100% rename from CANDy v2.0.ipynb rename to archive/CANDy v2.0.ipynb diff --git a/pyproject.toml b/pyproject.toml index 873df54..8890d1a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "candy-cazyme" -version = "3.0.0.dev0" +version = "3.0.0" description = "Automated analysis of domain architectures in carbohydrate-active enzymes (CAZymes)" readme = "README.md" license = { text = "MIT" } @@ -41,6 +41,11 @@ dev = [ "pytest-mock>=3.14", ] +[project.urls] +Homepage = "https://github.com/PyEED/CANDy" +Repository = "https://github.com/PyEED/CANDy" +Issues = "https://github.com/PyEED/CANDy/issues" + [project.scripts] candy = "candy.cli:app" From 7cdea36c75145d8015a9101e61b73038221771c2 Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Fri, 7 Aug 2026 11:24:05 +0200 Subject: [PATCH 6/7] Fail fast with a clear message on invalid CAZy family codes fetch_family_page previously retried a 404 (nonexistent family) up to 5 times against the same permanent failure before surfacing a raw requests.HTTPError traceback. Now it raises a clear ValueError immediately, and the CLI catches ValueError from run_pipeline to print a one-line "Error: ..." message and exit non-zero instead of a full traceback. Co-Authored-By: Claude Sonnet 5 --- src/candy/cazy.py | 8 ++++++++ src/candy/cli.py | 6 +++++- tests/test_cazy.py | 28 ++++++++++++++++++++++++++++ tests/test_cli.py | 16 ++++++++++++++++ 4 files changed, 57 insertions(+), 1 deletion(-) diff --git a/src/candy/cazy.py b/src/candy/cazy.py index a96fde0..21bb14d 100644 --- a/src/candy/cazy.py +++ b/src/candy/cazy.py @@ -65,6 +65,14 @@ def fetch_family_page(family: str, max_retries: int = 5, retry_delay: float = 5. response.raise_for_status() return response.text except requests.RequestException as exc: + if isinstance(exc, requests.HTTPError) and exc.response is not None and exc.response.status_code == 404: + # A 404 means the family doesn't exist on CAZy -- permanent, not + # worth retrying, and deserves a clear message instead of a + # raw traceback surfacing from deep inside the retry loop. + raise ValueError( + f"CAZy family '{family}' not found ({url} returned 404). " + "Check the family code, e.g. 'GH173' or 'GH5_1'." + ) from exc attempts += 1 if attempts >= max_retries: raise diff --git a/src/candy/cli.py b/src/candy/cli.py index 078218e..b5a124c 100644 --- a/src/candy/cli.py +++ b/src/candy/cli.py @@ -160,7 +160,11 @@ def main( build_tree=build_tree, ) - result = run_pipeline(config) + try: + result = run_pipeline(config) + except ValueError as exc: + typer.echo(f"Error: {exc}", err=True) + raise typer.Exit(code=1) from None typer.echo(f"\nDone. Results written to {result.jobname_dir}") typer.echo(f" Database: {result.database_path}") diff --git a/tests/test_cazy.py b/tests/test_cazy.py index b278dea..039c2d1 100644 --- a/tests/test_cazy.py +++ b/tests/test_cazy.py @@ -2,11 +2,14 @@ from unittest.mock import patch import pandas as pd +import pytest +import requests from candy.cazy import ( extract_characterized_ids, extract_ncbi_ids, fetch_characterized_sequences, + fetch_family_page, fetch_sequences_fasta, parse_family_table, ) @@ -51,6 +54,31 @@ def test_extract_characterized_ids_groups_by_taxonomy_header(): assert activity_dict == {"P12345.1": "1.2.4.-", "Q99999.2": "3.2.1.-"} +def test_fetch_family_page_raises_clear_error_on_404_without_retrying(): + # Regression test: a nonexistent family code (e.g. a typo) used to be + # retried 5 times against the same permanent 404 before surfacing a raw + # requests.HTTPError traceback. It should now fail fast with a clear + # message instead. + calls = {"n": 0} + + class FakeResponse: + status_code = 404 + + def raise_for_status(self): + raise requests.HTTPError("404 Client Error: Not Found", response=self) + + def fake_get(url, timeout): + calls["n"] += 1 + return FakeResponse() + + with patch("candy.cazy.requests.get", side_effect=fake_get), patch("candy.cazy.time.sleep") as mock_sleep: + with pytest.raises(ValueError, match="GH999"): + fetch_family_page("GH999") + + assert calls["n"] == 1 + mock_sleep.assert_not_called() + + class FakeHandle: """Stand-in for the object Bio.Entrez.efetch returns (used as a context manager).""" diff --git a/tests/test_cli.py b/tests/test_cli.py index c00cafa..78581ef 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -207,6 +207,22 @@ def test_cli_db_preference_unknown_name_is_a_clean_error(tmp_path, monkeypatch): assert "Unknown database name" in result.output +def test_cli_pipeline_value_error_is_reported_cleanly_without_a_traceback(monkeypatch): + # Regression test: a ValueError from run_pipeline (e.g. an unknown CAZy + # family) used to propagate as a raw traceback. It should now print a + # one-line "Error: ..." message and exit non-zero, no traceback. + def fake_run_pipeline(config): + raise ValueError("CAZy family 'GH999' not found (... returned 404). Check the family code.") + + monkeypatch.setattr("candy.cli.run_pipeline", fake_run_pipeline) + + result = runner.invoke(app, ["GH999", "--email", "you@example.com"]) + + assert result.exit_code == 1 + assert "Error: CAZy family 'GH999' not found" in result.output + assert result.exception is None or isinstance(result.exception, SystemExit) + + def test_cli_no_db_preference_uses_default_order(tmp_path, monkeypatch): from candy.config import DEFAULT_DATABASE_PREFERENCE From ca6367c051fbca3f2b180a1486103f5547002619 Mon Sep 17 00:00:00 2001 From: Alex Windels Date: Fri, 7 Aug 2026 11:30:54 +0200 Subject: [PATCH 7/7] Fix CI-only test failure from ANSI codes in Typer's error panel test_cli_bare_invocation_is_candy_target_no_subcommand_needed failed on every CI matrix job (but passed locally) because Typer's Rich-rendered error panel is colorized there, and the ANSI escape codes happened to land mid-phrase at a wrapped line boundary, breaking a substring check that only stripped box-drawing characters, not color codes. Add an ANSI-stripping helper and use it for all Rich-panel-derived assertions; verified locally by forcing the same colorized rendering with FORCE_COLOR=1. Co-Authored-By: Claude Sonnet 5 --- tests/test_cli.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/tests/test_cli.py b/tests/test_cli.py index 78581ef..0995a94 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,9 +1,29 @@ +import re + from typer.testing import CliRunner from candy.cli import _sanitize_jobname, _try_parse_family, app runner = CliRunner() +_ANSI_RE = re.compile(r"\x1b\[[0-9;]*m") + + +def _plain_text(output: str) -> str: + """Strip ANSI color codes and box-drawing characters from Typer's Rich-rendered output. + + Whether Rich colorizes its error panels depends on how it detects the + terminal (observed to differ between a local shell and a CI runner, even + on identical typer/rich versions) -- when it does, the color codes can + land mid-word at a wrapped line boundary, splitting a phrase that should + otherwise be a contiguous substring. Stripping them makes assertions + robust to that regardless of environment. + """ + text = _ANSI_RE.sub("", output) + for border_char in "─│┌┐└┘": + text = text.replace(border_char, " ") + return " ".join(text.split()) + def test_try_parse_family_no_subfamily(): assert _try_parse_family("GH5") == ("GH", 5, None) @@ -30,10 +50,7 @@ def test_cli_help(): def test_cli_bare_invocation_is_candy_target_no_subcommand_needed(): # `candy GH5` should work directly -- no `run` subcommand required. result = runner.invoke(app, ["not-a-real-family-and-not-a-file"]) - result_output = result.output - for border_char in "─│┌┐└┘": - result_output = result_output.replace(border_char, " ") - normalized = " ".join(result_output.split()) + normalized = _plain_text(result.output) assert result.exit_code != 0 assert "neither an existing FASTA file" in normalized assert "valid CAZy family code" in normalized @@ -204,7 +221,7 @@ def test_cli_db_preference_unknown_name_is_a_clean_error(tmp_path, monkeypatch): ) assert result.exit_code != 0 - assert "Unknown database name" in result.output + assert "Unknown database name" in _plain_text(result.output) def test_cli_pipeline_value_error_is_reported_cleanly_without_a_traceback(monkeypatch):