Rewrite CANDy as an installable Python package (CLI + API)

.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

.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

.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

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.

environment.yml

@@ -0,0 +1,17 @@
+# 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
+  - conda-forge
+dependencies:
+  - python>=3.10
+  - pip
+  - cd-hit
+  - pip:
+      - -e .

pyproject.toml

@@ -0,0 +1,60 @@
+[build-system]
+requires = ["setuptools>=68", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "candy-cazyme"
+version = "3.0.0"
+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",
+    "pyfamsa>=0.5",
+    "veryfasttree>=4.0",
+    "platformdirs>=4.0",
+]
+
+[project.optional-dependencies]
+gemini = ["google-genai>=0.3"]
+dev = [
+    "pytest>=8.0",
+    "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"
+
+[tool.setuptools.packages.find]
+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)",
+]

src/candy/__init__.py

@@ -0,0 +1,3 @@
+"""CANDy: automated analysis of domain architectures in carbohydrate-active enzymes."""
+
+__version__ = "3.0.0.dev0"

src/candy/alignment/__init__.py

@@ -0,0 +1,26 @@
+"""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 == "famsa":
+        from candy.alignment.famsa import FamsaAligner
+
+        return FamsaAligner()
+    if name == "mafft":
+        from candy.alignment.mafft import MafftAligner
+
+        return MafftAligner()
+    raise ValueError(f"Unknown alignment tool: {name}")

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

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

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