From 7931367dc89937688fc6c60fa79fd55bd70e0436 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Wed, 23 Sep 2026 16:26:00 +0100 Subject: [PATCH 1/5] feat(sist): add sist Python CLI package --- pyproject.toml | 44 +++- src/sist/__init__.py | 9 + src/sist/__main__.py | 6 + src/sist/_index.py | 27 +++ src/sist/argspec.py | 207 ++++++++++++++++++ src/sist/binaries.py | 64 ++++++ src/sist/cli.py | 54 +++++ src/sist/energetics.py | 260 +++++++++++++++++++++++ src/sist/ir_finder.py | 443 +++++++++++++++++++++++++++++++++++++++ src/sist/runner.py | 179 ++++++++++++++++ tests/conftest.py | 93 +++++++- tests/test_energetics.py | 188 +++++++++++++++++ tests/test_ir_finder.py | 123 +++++++++++ tests/test_python_cli.py | 104 +++++++++ 14 files changed, 1784 insertions(+), 17 deletions(-) create mode 100644 src/sist/__init__.py create mode 100644 src/sist/__main__.py create mode 100644 src/sist/_index.py create mode 100644 src/sist/argspec.py create mode 100644 src/sist/binaries.py create mode 100644 src/sist/cli.py create mode 100644 src/sist/energetics.py create mode 100644 src/sist/ir_finder.py create mode 100644 src/sist/runner.py create mode 100644 tests/test_energetics.py create mode 100644 tests/test_ir_finder.py create mode 100644 tests/test_python_cli.py diff --git a/pyproject.toml b/pyproject.toml index 83c0c35..23de95a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,11 +4,14 @@ build-backend = "setuptools.build_meta" [project] name = "SIST" -version = "0.0.1" +dynamic = ["version"] description = "Stress-Induced Structural Transitions in superhelical DNA." readme = "README.md" requires-python = ">=3.12" -dependencies = [] +dependencies = [ + "biopython>=1.85,<2.0", + "beautifulsoup4>=4.12,<5.0", +] authors = [ { name = "Craig Benham", email = "cjbenham@ucdavis.edu" }, @@ -28,15 +31,21 @@ maintainers = [ [project.urls] Repository = "https://github.com/CCPBioSim/SIST" +[project.scripts] +sist = "sist.cli:main" + [project.optional-dependencies] testing = [ "pytest>=9.0,<10.0", + "mypy>=1.14,<2.0", ] -[tool.setuptools] -# SIST is currently a Perl/C++ application. This metadata manages the Python -# test harness and its development dependencies; it does not package SIST. -packages = [] +[tool.setuptools.packages.find] +where = ["src"] +include = ["sist*"] + +[tool.setuptools.dynamic] +version = { attr = "sist.__version__" } [tool.pytest.ini_options] testpaths = ["tests"] @@ -45,3 +54,26 @@ markers = [ "regression: end-to-end regression tests against reference outputs", "slow: long-running tests excluded from the default test suite", ] + +[tool.ruff] +line-length = 88 +target-version = "py312" + +[tool.ruff.lint] +select = ["E", "F", "I", "B", "UP"] + +[tool.ruff.format] +quote-style = "double" +indent-style = "space" + +[tool.mypy] +python_version = "3.12" +packages = ["sist"] +mypy_path = "src" +disallow_untyped_defs = true +warn_unused_ignores = true +warn_redundant_casts = true + +[[tool.mypy.overrides]] +module = ["Bio.*", "bs4.*"] +ignore_missing_imports = true diff --git a/src/sist/__init__.py b/src/sist/__init__.py new file mode 100644 index 0000000..5242464 --- /dev/null +++ b/src/sist/__init__.py @@ -0,0 +1,9 @@ +""" +sist + +SIST calculates stress-induced structural transition probabilities in +superhelical DNA, including strand separation, Z-DNA formation, +cruciform extrusion, and competition between these transitions. +""" + +__version__ = "1.0.0" diff --git a/src/sist/__main__.py b/src/sist/__main__.py new file mode 100644 index 0000000..b4baf1a --- /dev/null +++ b/src/sist/__main__.py @@ -0,0 +1,6 @@ +from __future__ import annotations + +from sist.cli import main + +if __name__ == "__main__": + main() diff --git a/src/sist/_index.py b/src/sist/_index.py new file mode 100644 index 0000000..d19c55d --- /dev/null +++ b/src/sist/_index.py @@ -0,0 +1,27 @@ +"""Shared helper for indexing with negative wraparound and no bounds errors. + +Used by modules that need to look at neighbouring positions in an +alignment without special-casing the edges of the sequence. +""" + +from __future__ import annotations + +from collections.abc import Sequence + + +def index_or_none[T](sequence: Sequence[T], index: int) -> T | None: + """Index a sequence, returning ``None`` instead of raising if out of range. + + Negative indices wrap from the end, as with normal Python indexing. + + Args: + sequence: Sequence to index. + index: Index to read. + + Returns: + The element at ``index``, or ``None`` if out of range. + """ + if -len(sequence) <= index < len(sequence): + return sequence[index] + + return None diff --git a/src/sist/argspec.py b/src/sist/argspec.py new file mode 100644 index 0000000..a681f15 --- /dev/null +++ b/src/sist/argspec.py @@ -0,0 +1,207 @@ +"""Command-line argument specification for sist. + +This module provides a declarative argument specification (`ARG_SPECS`) +used to build an `argparse.ArgumentParser` for the sist CLI's short and +long-form flags. + +""" + +from __future__ import annotations + +import argparse +from dataclasses import dataclass +from typing import Any + +ALGORITHM_CHOICES = ("M", "Z", "C", "A") + + +@dataclass(frozen=True) +class ArgSpec: + """Argument specification used to build an argparse parser. + + Attributes: + flags: Option strings, e.g. ``("-f", "--file")``. + dest: Attribute name the parsed value is stored under. + help: Help text shown in CLI usage. + default: Default value if not provided via CLI. + type: Python type for parsing, such as int, float, or str. + action: Optional argparse action, such as "store_true". + choices: Optional set of allowed values. + required: Whether the argument must be provided. + """ + + dest: str + help: str + flags: tuple[str, ...] + default: Any = None + type: Any = None + action: str | None = None + choices: tuple[str, ...] | None = None + required: bool = False + + +ARG_SPECS: tuple[ArgSpec, ...] = ( + ArgSpec( + dest="file", + flags=("-f", "--file"), + type=str, + required=True, + help="Sequence file to analyse.", + ), + ArgSpec( + dest="algorithm", + flags=("-a", "--algorithm"), + type=str, + required=True, + choices=ALGORITHM_CHOICES, + help=( + "Algorithm to run: M (melting), Z (Z-DNA), C (cruciform), or A " + "(competition between all three). C and A require an Inverted " + "Repeat Finder (IRF) executable on PATH." + ), + ), + ArgSpec( + dest="temperature", + flags=("-T", "--temperature"), + type=float, + default=310.0, + help="Temperature in Kelvin. Defaults to %(default)s.", + ), + ArgSpec( + dest="superhelical_density", + flags=("-s", "--superhelical-density"), + type=float, + default=0.06, + help="Superhelical density. Defaults to %(default)s.", + ), + ArgSpec( + dest="salt", + flags=("-i", "--salt"), + type=float, + default=0.01, + help="Ionic strength (salt concentration) in mol/L. Defaults to %(default)s.", + ), + ArgSpec( + dest="energy_threshold", + flags=("-th", "--energy-threshold"), + type=float, + default=12.0, + help=( + "Energy threshold for the transition search. Defaults to " + "%(default)s; 10 is recommended for the competition algorithm." + ), + ), + ArgSpec( + dest="circular", + flags=("-c", "--circular"), + action="store_true", + help="Treat the molecule as circular. Defaults to linear.", + ), + ArgSpec( + dest="nearest_neighbor", + flags=("-n", "--nearest-neighbor"), + action="store_true", + help=( + "Use nearest-neighbor melting energetics instead of " + "copolymeric. Ignored for the Z-DNA and cruciform algorithms." + ), + ), + ArgSpec( + dest="print_base_pair", + flags=("-b", "--print-base-pair"), + action="store_true", + help="Print the base pair for each position.", + ), + ArgSpec( + dest="print_parameters", + flags=("-p", "--print-parameters"), + action="store_true", + help="Print algorithm parameters.", + ), + ArgSpec( + dest="print_ensemble_average", + flags=("-r", "--print-ensemble-average"), + action="store_true", + help="Print ensemble average results.", + ), + ArgSpec( + dest="output_file", + flags=("-o", "--output-file"), + type=str, + default=None, + help="Write output to this file instead of printing to stdout.", + ), +) + + +class SistArgumentParser: + """Builds and validates the sist CLI's argument parser.""" + + def __init__(self, arg_specs: tuple[ArgSpec, ...] | None = None) -> None: + """Initialise the parser builder. + + Args: + arg_specs: Optional override for argument specs. If omitted, + uses `ARG_SPECS`. + """ + self._arg_specs = arg_specs if arg_specs is not None else ARG_SPECS + + def build_parser(self) -> argparse.ArgumentParser: + """Build an ArgumentParser from the argument specs. + + Returns: + Configured argparse.ArgumentParser. + """ + parser = argparse.ArgumentParser( + prog="sist", + description=( + "SIST: stress-induced structural transition probabilities in " + "superhelical DNA." + ), + ) + + for spec in self._arg_specs: + kwargs: dict[str, Any] = {"dest": spec.dest, "help": spec.help} + + if spec.action is not None: + kwargs["action"] = spec.action + else: + kwargs["type"] = spec.type + kwargs["default"] = spec.default + + if spec.choices is not None: + kwargs["choices"] = spec.choices + + if spec.required: + kwargs["required"] = True + + parser.add_argument(*spec.flags, **kwargs) + + return parser + + @staticmethod + def validate(args: argparse.Namespace) -> None: + """Validate parsed arguments against sensible runtime constraints. + + Args: + args: Parsed CLI arguments. + + Raises: + ValueError: If a parameter is invalid. + """ + if args.temperature <= 0: + raise ValueError( + f"Invalid 'temperature': {args.temperature}. Temperature must " + "be greater than 0 K." + ) + + if args.salt <= 0: + raise ValueError( + f"Invalid 'salt': {args.salt}. Salt must be greater than 0 mol/L." + ) + + if args.energy_threshold <= 0: + raise ValueError( + f"Invalid 'energy_threshold': {args.energy_threshold}. It must " + "be greater than 0." + ) diff --git a/src/sist/binaries.py b/src/sist/binaries.py new file mode 100644 index 0000000..19ea429 --- /dev/null +++ b/src/sist/binaries.py @@ -0,0 +1,64 @@ +"""Locate the compiled qsidd binaries used by sist's orchestration layer.""" + +from __future__ import annotations + +import os +import sys +from pathlib import Path + + +class BinaryLocator: + """Locates the compiled qsidd binaries, for either trans_three or trans_compete.""" + + TRANS_THREE_ENV_VAR = "SIST_TRANS_THREE_BIN" + TRANS_COMPETE_ENV_VAR = "SIST_TRANS_COMPETE_BIN" + + def trans_three_binary(self) -> Path: + """Locate the melting/Z-DNA/cruciform qsidd binary. + + Returns: + Path to the trans_three qsidd executable. + """ + return self._resolve( + self.TRANS_THREE_ENV_VAR, "libexec/sist/src/trans_three/qsidd" + ) + + def trans_compete_binary(self) -> Path: + """Locate the competition qsidd binary. + + Returns: + Path to the trans_compete qsidd executable. + """ + return self._resolve( + self.TRANS_COMPETE_ENV_VAR, "libexec/sist/src/trans_compete/qsidd" + ) + + @staticmethod + def _resolve(env_var: str, installed_relative_path: str) -> Path: + """Resolve a qsidd binary path via an env var override or install layout. + + Args: + env_var: Environment variable that, if set, is used directly as + the binary path. + installed_relative_path: Path to the binary relative to + ``sys.prefix``, matching the conda package's install layout. + + Returns: + Path to an existing qsidd executable. + + Raises: + FileNotFoundError: If neither the env var nor the install + layout resolves to an existing file. + """ + override = os.environ.get(env_var) + path = ( + Path(override) if override else Path(sys.prefix) / installed_relative_path + ) + + if not path.is_file(): + raise FileNotFoundError( + f"Could not find a qsidd executable at {path}. Set {env_var} to " + "its path, or install sist so it is available there." + ) + + return path diff --git a/src/sist/cli.py b/src/sist/cli.py new file mode 100644 index 0000000..d1e4e97 --- /dev/null +++ b/src/sist/cli.py @@ -0,0 +1,54 @@ +"""Command-line entry point for sist. + +The entry point is intentionally small and only responsible for: + 1) Building the parser and parsing arguments. + 2) Validating those arguments. + 3) Constructing a SistRunner and running it. + 4) Handling fatal errors with a non-zero exit code. +""" + +from __future__ import annotations + +import logging + +from sist.argspec import SistArgumentParser +from sist.runner import SistRunner + +logger = logging.getLogger(__name__) + + +def main(argv: list[str] | None = None) -> None: + """Parse CLI arguments and run a SIST calculation. + + Args: + argv: Argument list to parse, or None to use `sys.argv[1:]`. + + Raises: + SystemExit: Exits with status code 1 on any unhandled exception. + """ + logging.basicConfig(level=logging.WARNING, format="%(levelname)s: %(message)s") + + argument_parser = SistArgumentParser() + parser = argument_parser.build_parser() + args = parser.parse_args(argv) + + try: + argument_parser.validate(args) + runner = SistRunner( + file=args.file, + algorithm=args.algorithm, + temperature=args.temperature, + superhelical_density=args.superhelical_density, + salt=args.salt, + energy_threshold=args.energy_threshold, + circular=args.circular, + nearest_neighbor=args.nearest_neighbor, + print_base_pair=args.print_base_pair, + print_parameters=args.print_parameters, + print_ensemble_average=args.print_ensemble_average, + output_file=args.output_file, + ) + runner.run() + except Exception: + logger.exception("Fatal error during SIST calculation") + raise SystemExit(1) from None diff --git a/src/sist/energetics.py b/src/sist/energetics.py new file mode 100644 index 0000000..933b1b5 --- /dev/null +++ b/src/sist/energetics.py @@ -0,0 +1,260 @@ +"""Thermodynamic energy calculations for cruciform extrusion candidates. + +Used to score candidate inverted repeats (IRs) reported by the Inverted +Repeat Finder (IRF) tool before they are handed to the ``qsidd`` +cruciform/competition binaries. + +Reference: Zhabinskaya & Benham, Nucleic Acids Res, 41(21), 9610. +""" + +from __future__ import annotations + +import math +from collections.abc import Sequence + +from sist._index import index_or_none + +BASE_PAIR_COMPLEMENTS: dict[str, str] = {"A": "T", "T": "A", "C": "G", "G": "C"} + +WATSON_CRICK_ENERGIES: dict[str, float] = { + "AA": -1.0, + "TT": -1.0, + "AT": -0.88, + "TA": -0.58, + "CA": -1.45, + "AC": -1.45, + "GT": -1.44, + "TG": -1.44, + "CT": -1.28, + "TC": -1.28, + "GA": -1.3, + "AG": -1.3, + "CG": -2.17, + "GC": -2.24, + "GG": -1.84, + "CC": -1.84, +} + +MISMATCH_ENERGIES: dict[str, float] = { + "GA.CA": 0.17, + "GA.CC": 0.81, + "GA.CG": -0.25, + "GC.CA": 0.47, + "GC.CC": 0.79, + "GC.CT": 0.62, + "GG.CA": -0.52, + "GG.CG": -1.11, + "GG.CT": 0.08, + "GT.CC": 0.98, + "GT.CG": -0.59, + "GT.CT": 0.45, + "CA.GA": 0.43, + "CA.GC": 0.75, + "CA.GG": 0.03, + "CC.GA": 0.79, + "CC.GC": 0.70, + "CC.GT": 0.62, + "CG.GA": 0.11, + "CG.GG": -0.11, + "CG.GT": -0.47, + "CT.GC": 0.40, + "CT.GG": -0.32, + "CT.GT": -0.12, + "AA.TA": 0.61, + "AA.TC": 0.88, + "AA.TG": 0.14, + "AC.TA": 0.77, + "AC.TC": 1.33, + "AC.TT": 0.64, + "AG.TA": 0.02, + "AG.TG": -0.13, + "AG.TT": 0.71, + "AT.TC": 0.73, + "AT.TG": 0.07, + "AT.TT": 0.69, + "TA.AA": 0.69, + "TA.AC": 0.92, + "TA.AG": 0.42, + "TC.AA": 1.33, + "TC.AC": 1.05, + "TC.AT": 0.97, + "TG.AA": 0.74, + "TG.AG": 0.44, + "TG.AT": 0.43, + "TT.AC": 0.75, + "TT.AG": 0.34, + "TT.AT": 0.68, +} + +IMPERFECTION_ENERGY = 2.42 + + +class CruciformEnergetics: + """Energy calculations for cruciform extrusion candidates. + + All calculations are evaluated at a fixed temperature and salt + concentration, set once when the calculator is constructed. + + Attributes: + temperature: Temperature in Kelvin. + salt: Salt (ionic strength) concentration in mol/L. + """ + + def __init__(self, *, temperature: float, salt: float) -> None: + """Initialise the calculator. + + Args: + temperature: Temperature in Kelvin. + salt: Salt (ionic strength) concentration in mol/L. + """ + self.temperature = temperature + self.salt = salt + + @property + def gas_constant_rt(self) -> float: + """RT in kcal/mol at this calculator's temperature.""" + return 1.9872 * self.temperature / 1000.0 + + def melting_energy(self, base: str) -> float: + """Compute the melting energy of a single base pair. + + Args: + base: The base at this position. ``A``/``T`` (case-insensitive) + use the AT melting energy; anything else uses the GC + melting energy. + + Returns: + Melting energy in kcal/mol. + """ + at_melting_temperature = 354.65 + 16.6 * math.log10(self.salt) + gc_melting_temperature = at_melting_temperature + 41 + + at_energy = 7.2464 * (1 - self.temperature / at_melting_temperature) + gc_energy = 9.0172 * (1 - self.temperature / gc_melting_temperature) + + if base in ("a", "t", "A", "T"): + return at_energy + + return gc_energy + + def loop_energy(self, loop_sequence: Sequence[str]) -> float: + """Compute the free energy of a cruciform loop. + + Args: + loop_sequence: Bases making up the loop. + + Returns: + Loop free energy in kcal/mol, including the loop entropy term. + """ + energy = sum(self.melting_energy(base) for base in loop_sequence) + energy += 2 * 2.44 * self.gas_constant_rt * math.log(len(loop_sequence)) + + return energy + + @staticmethod + def mismatch_energy(triplet1: str, triplet2: str) -> float: + """Compute the energy penalty of a base-pair mismatch. + + Args: + triplet1: Three consecutive bases on one strand, centred on the + mismatch. + triplet2: The corresponding three bases on the opposite strand. + + Returns: + Mismatch energy penalty in kcal/mol. + """ + watson_crick_energy = ( + WATSON_CRICK_ENERGIES[triplet1[0] + triplet1[1]] + + WATSON_CRICK_ENERGIES[triplet1[1] + triplet1[2]] + ) + + stack1 = f"{triplet1[0]}{triplet1[1]}.{triplet2[0]}{triplet2[1]}" + stack2 = f"{triplet2[2]}{triplet2[1]}.{triplet1[2]}{triplet1[1]}" + mismatch_stack_energy = MISMATCH_ENERGIES[stack1] + MISMATCH_ENERGIES[stack2] + + return mismatch_stack_energy - watson_crick_energy + + def imperfection_energies( + self, right_arm: Sequence[str], left_arm: Sequence[str] + ) -> list[float]: + """Compute per-position imperfection energies along an IR arm alignment. + + ``right_arm`` and ``left_arm`` are aligned, equal-length sequences of + single-character tokens describing one arm of a candidate inverted + repeat: a base letter for a bulge, ``"-"`` for a gap, or (on + ``right_arm``) ``"*"`` for a perfectly matched position. + + Args: + right_arm: Right-arm alignment tokens, starting with ``"*"``. + left_arm: Left-arm alignment tokens, aligned with ``right_arm``. + + Returns: + One imperfection energy (kcal/mol) per position, ``0.0`` where + there is no imperfection. + + Raises: + ValueError: If ``right_arm`` does not start with ``"*"``. + """ + if right_arm[0] != "*": + raise ValueError("Doesn't start with a star") + + energies = [0.0] * len(left_arm) + + for index, left_base in enumerate(left_arm): + if left_base == "-": + energies[index] = IMPERFECTION_ENERGY + self.melting_energy( + right_arm[index] + ) + + for index, right_base in enumerate(right_arm): + if right_base == "-": + energies[index] = IMPERFECTION_ENERGY + self.melting_energy( + left_arm[index] + ) + continue + + if right_base in ("*", "-"): + continue + + if left_arm[index] == "-": + continue + + previous_right = index_or_none(right_arm, index - 1) + next_right = index_or_none(right_arm, index + 1) + + if previous_right != "*" or next_right != "*": + energies[index] = ( + 2 * IMPERFECTION_ENERGY + + self.melting_energy(right_base) + + self.melting_energy(left_arm[index]) + ) + continue + + previous_left = index_or_none(left_arm, index - 1) + next_left = index_or_none(left_arm, index + 1) + assert previous_left is not None + assert next_left is not None + + if previous_left != "-" and next_left != "-": + triplet1 = f"{previous_left}{left_arm[index]}{next_left}" + triplet2 = ( + f"{BASE_PAIR_COMPLEMENTS[previous_left]}" + f"{right_base}" + f"{BASE_PAIR_COMPLEMENTS[next_left]}" + ) + energies[index] = self.mismatch_energy(triplet1, triplet2) + + return energies + + def cruciform_initiation_energy(self) -> float: + """Compute the cruciform initiation energy, Ecr. + + Returns: + Cruciform initiation energy in kcal/mol. + """ + return ( + 192.5 + - self.temperature * 0.565 + - 4 * self.melting_energy("A") + - 2 * 2.44 * self.gas_constant_rt * math.log(4) + ) diff --git a/src/sist/ir_finder.py b/src/sist/ir_finder.py new file mode 100644 index 0000000..0f947fa --- /dev/null +++ b/src/sist/ir_finder.py @@ -0,0 +1,443 @@ +"""Inverted repeat scoring for the cruciform and competition algorithms. + +This module runs the external Inverted Repeat Finder (``irf``) tool +against a sequence, parses its HTML report, and scores each candidate +inverted repeat (IR) using the energy model in :mod:`sist.energetics`. +The result is the ``"start,length,energy,|"`` string consumed by +``qsidd -X``. + +IRF download page: http://tandem.bu.edu/irf/irf.download.html +Reference: Zhabinskaya & Benham, Nucleic Acids Res, 41(21), 9610. +""" + +from __future__ import annotations + +import subprocess +from pathlib import Path + +from Bio import SeqIO +from bs4 import BeautifulSoup + +from sist._index import index_or_none +from sist.energetics import BASE_PAIR_COMPLEMENTS, CruciformEnergetics + +_DNA_BASES = frozenset("ATCG") + + +class IRFinder: + """Scores candidate inverted repeats for the cruciform/competition algorithms. + + Attributes: + temperature: Temperature in Kelvin. + shape: Either ``"linear"`` or ``"circular"``. + energetics: Energy calculator for this temperature. Uses a fixed + salt concentration independent of (and always overriding) any + ``-i``/salt value passed to the ``qsidd`` binaries. + """ + + IRF_EXECUTABLE = "irf" + IRF_MATCH = 2 + IRF_MISMATCH = 10 + IRF_DELTA = 10 + IRF_PM = 80 + IRF_PI = 10 + IRF_MIN_SCORE = 20 + IRF_MAX_LENGTH = 10000 + IRF_MAX_LOOP = 100 + IRF_MIN_LOOP = 3 + + CRUCIFORM_SALT = 0.01 + + def __init__(self, *, temperature: float, shape: str) -> None: + """Initialise the finder. + + Args: + temperature: Temperature in Kelvin. + shape: Either ``"linear"`` or ``"circular"``. + """ + self.temperature = temperature + self.shape = shape + self.energetics = CruciformEnergetics( + temperature=temperature, salt=self.CRUCIFORM_SALT + ) + + self._seen: dict[tuple[int, int, int], int] = {} + self._ir_string_parts: list[str] = ["1,0,10000,|"] + + def compute_cruciform_energy_string(self, file_path: str | Path) -> str: + """Score every candidate inverted repeat in a sequence. + + Converts the input sequence to IRF's required format, runs IRF + (twice, for circular sequences, to also score IRs spanning the + origin), and returns the ``"start,length,energy,|..."`` string + consumed by ``qsidd -X``. + + Args: + file_path: Path to the input FASTA sequence file. + + Returns: + The energy-scored inverted repeat string for the sequence. + """ + one_line_path = self.convert_to_one_line(Path(file_path)) + sequence = self.read_sequence(one_line_path) + + cruciform_initiation_energy = self.energetics.cruciform_initiation_energy() + + self._score_reports( + one_line_path, + circular_pass=False, + sequence=sequence, + sequence_length=len(sequence), + cruciform_initiation_energy=cruciform_initiation_energy, + ) + + if self.shape == "circular": + sequence_length, shifted_sequence = self.sequence_length_and_shift( + one_line_path, circular=True + ) + + circular_path = Path(f"circ.{one_line_path.name}") + circular_path.write_text( + f">{circular_path.name}\n{shifted_sequence}", encoding="utf-8" + ) + + self._score_reports( + circular_path, + circular_pass=True, + sequence=shifted_sequence, + sequence_length=sequence_length, + cruciform_initiation_energy=cruciform_initiation_energy, + ) + + return "".join(self._ir_string_parts) + + @staticmethod + def convert_to_one_line(sequence_path: Path) -> Path: + """Convert a FASTA file into IRF's required single-line format. + + Uses Biopython to parse well-formed FASTA input. A file with no + ``>`` header isn't valid FASTA, so Biopython rejects it; in that + case, the first line is discarded entirely (not kept as sequence + data) and a synthetic header is written instead. + + Args: + sequence_path: Path to the input sequence file. The output + file is written in the current working directory, named + after this path's basename. + + Returns: + Path to the newly written ``one_line.`` file. + """ + output_path = Path(f"one_line.{sequence_path.name}") + + try: + records = list(SeqIO.parse(sequence_path, "fasta")) + except ValueError: + records = None + + if records: + header = records[0].description + sequence = str(records[0].seq) + elif records == []: + output_path.write_text("", encoding="utf-8") + return output_path + else: + lines = sequence_path.read_text(encoding="utf-8").splitlines() + header = output_path.name + sequence = "".join(lines[1:]) + + output_path.write_text(f">{header}\n{sequence}", encoding="utf-8") + + return output_path + + @staticmethod + def read_sequence(one_line_path: Path) -> str: + """Read the sequence out of a one-line FASTA file. + + Args: + one_line_path: Path to a file in the single-line FASTA format + produced by :meth:`convert_to_one_line`. + + Returns: + The bare sequence, with no header line. + """ + return str(SeqIO.read(one_line_path, "fasta").seq) + + @classmethod + def sequence_length_and_shift( + cls, one_line_path: Path, *, circular: bool + ) -> tuple[int, str]: + """Read a one-line sequence file and optionally shift it to the middle. + + Args: + one_line_path: Path to a file in the single-line FASTA format + produced by :meth:`convert_to_one_line`. + circular: If True, also compute the sequence with its start + position shifted to the middle (for circular re-indexing). + + Returns: + A tuple of (sequence length, shifted sequence or "" if not + circular). + """ + sequence = cls.read_sequence(one_line_path) + sequence_length = len(sequence) + shifted_sequence = "" + + if circular: + midpoint = sequence_length // 2 + shifted_sequence = sequence[midpoint:] + sequence[:midpoint] + + return sequence_length, shifted_sequence + + @staticmethod + def parse_loop_header(line: str) -> tuple[int, int, int, int, int]: + """Parse an IRF report ``Loop:`` line. + + Args: + line: A line such as + ``" Indices: 2982--2992,2996--3006 Loop: 3 Score: 22"``. + + Returns: + A tuple of (loop_length, loop_start, ir_start, ir_end, ir_length). + """ + tokens = line.split() + loop_length = int(tokens[3]) + + left_range, right_range = tokens[1].split(",") + ir_start, left_end = (int(value) for value in left_range.split("--")) + _right_start, ir_end = (int(value) for value in right_range.split("--")) + + loop_start = left_end + 1 + ir_length = ir_end - ir_start + 1 + + return loop_length, loop_start, ir_start, ir_end, ir_length + + @staticmethod + def parse_arm(line: str) -> list[str]: + """Parse one alignment token out of an IRF arm line. + + Args: + line: A line such as ``" 2982 >> AAACCACCGCT >> 2992"`` or + ``" 3006 << *********** << 2996"``. + + Returns: + The arm's alignment tokens, one character each. + """ + return list(line.split()[2]) + + @staticmethod + def _read_report_lines(report_path: Path) -> list[str]: + """Read an IRF HTML report's text content, one entry per line. + + Args: + report_path: Path to an IRF ``*.txt.html`` report. + + Returns: + The report's text content (HTML tags stripped), split into + lines. + """ + soup = BeautifulSoup(report_path.read_text(encoding="utf-8"), "html.parser") + pre = soup.find("pre") + text = pre.get_text() if pre is not None else soup.get_text() + + return text.splitlines() + + @classmethod + def _extend_loop_to_minimum( + cls, + loop_bases: list[str], + left_arm: list[str], + right_arm: list[str], + *, + loop_start: int, + ) -> tuple[list[str], int, int]: + """Borrow bases from the arms until the loop reaches ``IRF_MIN_LOOP``. + + Args: + loop_bases: Bases in the loop, mutated in place. + left_arm: Left-arm alignment tokens. + right_arm: Right-arm alignment tokens, aligned with + ``left_arm``. + loop_start: Loop start position, before extension. + + Returns: + A tuple of (extended loop bases, updated loop_start, + arm bases consumed by the extension). + """ + shorten_arm = 0 + index = -1 + + while len(loop_bases) < cls.IRF_MIN_LOOP: + left_value = index_or_none(left_arm, index) + + if left_value is not None and left_value in _DNA_BASES: + right_value = index_or_none(right_arm, index) + + if right_value == "-": + loop_bases.append(left_value) + elif right_value == "*": + loop_bases.append(left_value) + loop_bases.append(BASE_PAIR_COMPLEMENTS[left_value]) + else: + assert right_value is not None + loop_bases.append(left_value) + loop_bases.append(right_value) + else: + # Both arms may be exhausted here, giving None; melting_energy() + # already treats an unrecognised/missing base as GC, so the + # None is left to flow through unchanged. + loop_bases.append(index_or_none(right_arm, index)) # type: ignore[arg-type] + shorten_arm += 1 + + loop_start -= 1 + shorten_arm += 1 + index -= 1 + + return loop_bases, loop_start, shorten_arm + + def _score_reports( + self, + sequence_file: Path, + *, + circular_pass: bool, + sequence: str, + sequence_length: int, + cruciform_initiation_energy: float, + ) -> None: + """Run IRF against ``sequence_file`` and append scored IRs to the result. + + Args: + sequence_file: One-line FASTA file to run IRF against. + circular_pass: True when scoring the circularly shifted + sequence. + sequence: The bare sequence read from ``sequence_file``. + sequence_length: Length of the (unshifted) original sequence. + cruciform_initiation_energy: Precomputed Ecr for this + temperature. + """ + subprocess.run( + [ + self.IRF_EXECUTABLE, + str(sequence_file), + str(self.IRF_MATCH), + str(self.IRF_MISMATCH), + str(self.IRF_DELTA), + str(self.IRF_PM), + str(self.IRF_PI), + str(self.IRF_MIN_SCORE), + str(self.IRF_MAX_LENGTH), + str(self.IRF_MAX_LOOP), + ], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL, + check=False, + ) + + report_prefix = ( + f"{sequence_file}.{self.IRF_MATCH}.{self.IRF_MISMATCH}.{self.IRF_DELTA}." + f"{self.IRF_PM}.{self.IRF_PI}.{self.IRF_MIN_SCORE}.{self.IRF_MAX_LENGTH}." + f"{self.IRF_MAX_LOOP}" + ) + + report_index = 1 + + while ( + report_path := Path(f"{report_prefix}.{report_index}.txt.html") + ).is_file(): + report_index += 1 + + reading = False + skip = False + left_arm: list[str] = [] + right_arm: list[str] = [] + loop_length = loop_start = ir_start = ir_end = ir_length = 0 + + for line in self._read_report_lines(report_path): + if "frequent" in line: + reading = False + + if "Done" in line: + break + + if "Loop:" in line: + skip = False + loop_length, loop_start, ir_start, ir_end, ir_length = ( + self.parse_loop_header(line) + ) + + if loop_length > self.IRF_MAX_LOOP: + continue + + self._seen[(loop_start, loop_length, ir_start)] = ir_length + reading = True + + if not reading: + continue + + if ">>" in line and "LF" not in line: + left_arm.extend(self.parse_arm(line)) + + if "<<" in line and "RF" not in line: + right_arm.extend(self.parse_arm(line)) + + if "Statistics" in line: + loop_bases = list(sequence[loop_start : loop_start + loop_length]) + + loop_bases, loop_start, shorten_arm = self._extend_loop_to_minimum( + loop_bases, left_arm, right_arm, loop_start=loop_start + ) + loop_length = len(loop_bases) + + if circular_pass: + half = sequence_length // 2 + + if loop_start > half: + offset = half + (sequence_length % 2) + loop_start -= offset + ir_start -= offset + else: + loop_start += half + ir_start += half + + seen_length = self._seen.get( + (loop_start, loop_length, ir_start) + ) + + if seen_length is not None and seen_length == ir_length: + skip = True + left_arm.clear() + right_arm.clear() + + if skip: + continue + + loop_free_energy = self.energetics.loop_energy(loop_bases) + + arm_length = len(left_arm) - shorten_arm + left_arm = left_arm[:arm_length] + right_arm = right_arm[:arm_length] + + imperfection_energies = self.energetics.imperfection_energies( + right_arm, left_arm + ) + imperfection_energies.reverse() + + total_energy = cruciform_initiation_energy + loop_free_energy + ir_position = loop_start + ir_length_running = loop_length + extension_count = 0 + + for position in range(arm_length): + total_energy += 2 * imperfection_energies[position] + + if left_arm[arm_length - 1 - position] != "-": + extension_count += 1 + ir_position = loop_start - extension_count + ir_length_running = loop_length + 2 * extension_count + + self._ir_string_parts.append( + f"{ir_position},{ir_length_running},{total_energy},|" + ) + + left_arm = [] + right_arm = [] diff --git a/src/sist/runner.py b/src/sist/runner.py new file mode 100644 index 0000000..56734b3 --- /dev/null +++ b/src/sist/runner.py @@ -0,0 +1,179 @@ +"""Orchestration that dispatches a SIST calculation to the qsidd binaries. + +Melting (M) and Z-DNA (Z) calculations run the single-transition +``qsidd`` binary directly; cruciform (C) and competition (A) calculations +first score candidate inverted repeats via :mod:`sist.ir_finder`, then run +the single-transition or competition ``qsidd`` binary respectively. +""" + +from __future__ import annotations + +import subprocess +from pathlib import Path + +from sist.binaries import BinaryLocator +from sist.ir_finder import IRFinder + + +class SistRunner: + """Coordinate a single SIST calculation: build, run, and report on it. + + Attributes: + file: Path to the input FASTA sequence file. + algorithm: One of ``"M"`` (melting), ``"Z"`` (Z-DNA), + ``"C"`` (cruciform), or ``"A"`` (competition). + temperature: Temperature in Kelvin. + superhelical_density: Superhelical density (sigma). + salt: Salt (ionic strength) concentration in mol/L. + energy_threshold: Energy threshold for the transition search. + circular: Whether the input molecule is circular (default linear). + nearest_neighbor: Use nearest-neighbor melting energetics. + Ignored for the Z-DNA and cruciform algorithms. + print_base_pair: Print the base pair for each position. + print_parameters: Print algorithm parameters. + print_ensemble_average: Print ensemble average results. + output_file: Path to write output to, or None to print to stdout. + """ + + def __init__( + self, + *, + file: str, + algorithm: str, + temperature: float, + superhelical_density: float, + salt: float, + energy_threshold: float, + circular: bool, + nearest_neighbor: bool, + print_base_pair: bool, + print_parameters: bool, + print_ensemble_average: bool, + output_file: str | None, + ) -> None: + """Initialise the runner with a single calculation's configuration.""" + self.file = file + self.algorithm = algorithm + self.temperature = temperature + self.superhelical_density = superhelical_density + self.salt = salt + self.energy_threshold = energy_threshold + self.circular = circular + self.nearest_neighbor = nearest_neighbor + self.print_base_pair = print_base_pair + self.print_parameters = print_parameters + self.print_ensemble_average = print_ensemble_average + self.output_file = output_file + + self._binaries = BinaryLocator() + + @property + def shape(self) -> str: + """Molecule shape as the qsidd binaries expect it: linear or circular.""" + return "circular" if self.circular else "linear" + + def run(self) -> None: + """Run the calculation and write or print its output. + + Raises: + ValueError: If ``algorithm`` is not one of M, Z, C, or A. + subprocess.CalledProcessError: If the qsidd binary exits + non-zero. + """ + command = self._build_command() + result = subprocess.run(command, capture_output=True, text=True, check=True) + + if self.output_file: + Path(self.output_file).write_text(result.stdout, encoding="utf-8") + else: + print(result.stdout, end="") + + def _build_command(self) -> list[str]: + """Build the qsidd command line for this calculation's algorithm.""" + shared_parameters = self._shared_parameters() + + if self.algorithm == "M": + return [ + str(self._binaries.trans_three_binary()), + *self._flags(nearest_neighbor=self.nearest_neighbor), + *shared_parameters, + "-f", + self.file, + ] + + if self.algorithm == "Z": + return [ + str(self._binaries.trans_three_binary()), + *self._flags(nearest_neighbor=False), + *shared_parameters, + "-Z", + "-f", + self.file, + ] + + if self.algorithm == "C": + ir_string = self._cruciform_energy_string() + + return [ + str(self._binaries.trans_three_binary()), + *self._flags(nearest_neighbor=False), + *shared_parameters, + "-C", + "-X", + ir_string, + "-f", + self.file, + ] + + if self.algorithm == "A": + ir_string = self._cruciform_energy_string() + + return [ + str(self._binaries.trans_compete_binary()), + *self._flags(nearest_neighbor=self.nearest_neighbor), + *shared_parameters, + "-X", + ir_string, + "-f", + self.file, + ] + + raise ValueError( + f"Unknown algorithm: {self.algorithm!r}. Expected one of M, Z, C, A." + ) + + def _cruciform_energy_string(self) -> str: + """Score candidate inverted repeats for the C/A algorithms.""" + finder = IRFinder(temperature=self.temperature, shape=self.shape) + + return finder.compute_cruciform_energy_string(self.file) + + def _flags(self, *, nearest_neighbor: bool) -> list[str]: + """Build the qsidd boolean flags.""" + flags = [] + + if self.print_base_pair: + flags.append("-b") + if self.print_parameters: + flags.append("-p") + if self.print_ensemble_average: + flags.append("-r") + if nearest_neighbor: + flags.append("-n") + if self.circular: + flags.append("-c") + + return flags + + def _shared_parameters(self) -> list[str]: + """Build the qsidd numeric parameter flags shared by every algorithm.""" + return [ + "-T", + str(self.temperature), + "-s", + str(self.superhelical_density), + "-i", + str(self.salt), + "-t", + str(self.energy_threshold), + ] diff --git a/tests/conftest.py b/tests/conftest.py index 4f1bf48..0c7ce5f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -3,12 +3,12 @@ import os import shutil import subprocess +import sys from dataclasses import dataclass from pathlib import Path import pytest - REPOSITORY_ROOT = Path(__file__).resolve().parents[1] SIST_TRANSITIONS = ( @@ -118,8 +118,17 @@ def build_sist( @pytest.fixture(scope="session") -def sist_command( +def built_sist_copy( tmp_path_factory: pytest.TempPathFactory, +) -> Path: + """Build SIST's qsidd binaries once, shared by every command fixture.""" + + return build_sist(tmp_path_factory) + + +@pytest.fixture(scope="session") +def sist_command( + built_sist_copy: Path, ) -> list[str]: """ Return the SIST command under test. @@ -139,14 +148,45 @@ def sist_command( return [executable] - built_sist = build_sist(tmp_path_factory) - return [ "perl", - str(built_sist / "master.pl"), + str(built_sist_copy / "master.pl"), ] +@pytest.fixture(scope="session") +def python_sist_command( + built_sist_copy: Path, +) -> list[str]: + """ + Return the new `sist` Python CLI command under test. + + During conda-build testing, use the installed console script. Otherwise + run `python -m sist` against the maintained source tree, with its qsidd + binaries resolved via env vars pointing at the freshly built copy. + """ + + if os.environ.get("CONDA_BUILD_STATE") == "TEST": + executable = shutil.which("sist") + + if executable is None: + pytest.fail( + "The installed sist command was not found in PATH", + pytrace=False, + ) + + return [executable] + + os.environ["SIST_TRANS_THREE_BIN"] = str( + built_sist_copy / "src" / "trans_three" / "qsidd" + ) + os.environ["SIST_TRANS_COMPETE_BIN"] = str( + built_sist_copy / "src" / "trans_compete" / "qsidd" + ) + + return [sys.executable, "-m", "sist"] + + def run_sist_calculation( sist_command: list[str], tmp_path_factory: pytest.TempPathFactory, @@ -158,12 +198,7 @@ def run_sist_calculation( runtime_directory = tmp_path_factory.mktemp(f"sist-{name}") - source_input = ( - REPOSITORY_ROOT - / "tests" - / "data" - / "pbr322.toy.fa" - ) + source_input = REPOSITORY_ROOT / "tests" / "data" / "pbr322.toy.fa" runtime_input = runtime_directory / "pbr322.toy.fa" shutil.copy2(source_input, runtime_input) @@ -228,3 +263,39 @@ def transition_run( name=name, algorithm=algorithm, ) + + +@pytest.fixture(scope="session") +def python_competition_run( + python_sist_command: list[str], + tmp_path_factory: pytest.TempPathFactory, +) -> SistRun: + """Run the new sist Python CLI's competition calculation once.""" + + return run_sist_calculation( + python_sist_command, + tmp_path_factory, + name="python-competition", + algorithm="A", + ) + + +@pytest.fixture( + scope="session", + params=SIST_TRANSITIONS, +) +def python_transition_run( + request: pytest.FixtureRequest, + python_sist_command: list[str], + tmp_path_factory: pytest.TempPathFactory, +) -> SistRun: + """Run each sist Python CLI transition calculation once.""" + + name, algorithm = request.param + + return run_sist_calculation( + python_sist_command, + tmp_path_factory, + name=f"python-{name}", + algorithm=algorithm, + ) diff --git a/tests/test_energetics.py b/tests/test_energetics.py new file mode 100644 index 0000000..ba84744 --- /dev/null +++ b/tests/test_energetics.py @@ -0,0 +1,188 @@ +"""Unit tests for sist.energetics, hand-verified against the reference paper. + +Expected values are computed independently from the published formulas +(not by re-calling the module under test), so these tests catch +transcription errors. +""" + +from __future__ import annotations + +import math + +import pytest + +from sist import energetics + +TEMPERATURE = 310.0 +SALT = 0.01 + + +@pytest.fixture +def calculator() -> energetics.CruciformEnergetics: + return energetics.CruciformEnergetics(temperature=TEMPERATURE, salt=SALT) + + +def _reference_melting_energy(base: str, *, temperature: float, salt: float) -> float: + at_melting_temperature = 354.65 + 16.6 * math.log10(salt) + gc_melting_temperature = at_melting_temperature + 41 + at_energy = 7.2464 * (1 - temperature / at_melting_temperature) + gc_energy = 9.0172 * (1 - temperature / gc_melting_temperature) + return at_energy if base in ("a", "t", "A", "T") else gc_energy + + +@pytest.mark.parametrize("base", ["A", "a", "T", "t"]) +def test_melting_energy_at_bases( + calculator: energetics.CruciformEnergetics, base: str +) -> None: + expected = _reference_melting_energy(base, temperature=TEMPERATURE, salt=SALT) + + actual = calculator.melting_energy(base) + + assert actual == pytest.approx(expected) + + +@pytest.mark.parametrize("base", ["G", "g", "C", "c"]) +def test_melting_energy_gc_bases( + calculator: energetics.CruciformEnergetics, base: str +) -> None: + expected = _reference_melting_energy(base, temperature=TEMPERATURE, salt=SALT) + + actual = calculator.melting_energy(base) + + assert actual == pytest.approx(expected) + + +def test_gas_constant_rt(calculator: energetics.CruciformEnergetics) -> None: + assert calculator.gas_constant_rt == pytest.approx(1.9872 * TEMPERATURE / 1000.0) + + +def test_loop_energy_matches_reference_formula( + calculator: energetics.CruciformEnergetics, +) -> None: + loop = ["A", "T", "G", "C", "A"] + + expected = sum( + _reference_melting_energy(base, temperature=TEMPERATURE, salt=SALT) + for base in loop + ) + expected += 2 * 2.44 * (1.9872 * TEMPERATURE / 1000.0) * math.log(len(loop)) + + actual = calculator.loop_energy(loop) + + assert actual == pytest.approx(expected) + + +def test_mismatch_energy_matches_reference_tables() -> None: + triplet1, triplet2 = "GAA", "CAT" + + e_wc = ( + energetics.WATSON_CRICK_ENERGIES["GA"] + energetics.WATSON_CRICK_ENERGIES["AA"] + ) + e_miss = ( + energetics.MISMATCH_ENERGIES["GA.CA"] + energetics.MISMATCH_ENERGIES["TA.AA"] + ) + expected = e_miss - e_wc + + actual = energetics.CruciformEnergetics.mismatch_energy(triplet1, triplet2) + + assert actual == pytest.approx(expected) + + +def test_cruciform_initiation_energy_matches_reference_formula( + calculator: energetics.CruciformEnergetics, +) -> None: + rt = 1.9872 * TEMPERATURE / 1000.0 + expected = ( + 192.5 + - TEMPERATURE * 0.565 + - 4 * _reference_melting_energy("A", temperature=TEMPERATURE, salt=SALT) + - 2 * 2.44 * rt * math.log(4) + ) + + actual = calculator.cruciform_initiation_energy() + + assert actual == pytest.approx(expected) + + +def test_imperfection_energies_requires_star_at_start( + calculator: energetics.CruciformEnergetics, +) -> None: + with pytest.raises(ValueError, match="star"): + calculator.imperfection_energies(["A", "*"], ["A", "A"]) + + +def test_imperfection_energies_perfect_match_is_all_zero( + calculator: energetics.CruciformEnergetics, +) -> None: + right_arm = ["*", "*", "*", "*"] + left_arm = ["A", "T", "G", "C"] + + energies = calculator.imperfection_energies(right_arm, left_arm) + + assert energies == [0.0, 0.0, 0.0, 0.0] + + +def test_imperfection_energies_bulge_on_left_arm( + calculator: energetics.CruciformEnergetics, +) -> None: + # Right arm has a gap at position 1 -> bulge on the left arm base "T". + right_arm = ["*", "-", "*"] + left_arm = ["A", "T", "G"] + + energies = calculator.imperfection_energies(right_arm, left_arm) + + expected_bulge = energetics.IMPERFECTION_ENERGY + _reference_melting_energy( + "T", temperature=TEMPERATURE, salt=SALT + ) + + assert energies[1] == pytest.approx(expected_bulge) + assert energies[0] == pytest.approx(0.0) + assert energies[2] == pytest.approx(0.0) + + +def test_imperfection_energies_bulge_on_right_arm( + calculator: energetics.CruciformEnergetics, +) -> None: + # Left arm has a gap at position 1 -> bulge on the right arm base "C". + right_arm = ["*", "C", "*"] + left_arm = ["A", "-", "G"] + + energies = calculator.imperfection_energies(right_arm, left_arm) + + expected_bulge = energetics.IMPERFECTION_ENERGY + _reference_melting_energy( + "C", temperature=TEMPERATURE, salt=SALT + ) + + assert energies[1] == pytest.approx(expected_bulge) + + +def test_imperfection_energies_internal_loop( + calculator: energetics.CruciformEnergetics, +) -> None: + # Two adjacent open positions (right arm not '*' on both sides of index 1). + right_arm = ["*", "-", "A", "-", "*"] + left_arm = ["A", "T", "T", "T", "G"] + + energies = calculator.imperfection_energies(right_arm, left_arm) + + expected_internal_loop = ( + 2 * energetics.IMPERFECTION_ENERGY + + _reference_melting_energy("A", temperature=TEMPERATURE, salt=SALT) + + _reference_melting_energy("T", temperature=TEMPERATURE, salt=SALT) + ) + + assert energies[2] == pytest.approx(expected_internal_loop) + + +def test_imperfection_energies_mismatch( + calculator: energetics.CruciformEnergetics, +) -> None: + # A single mismatch flanked by perfect matches on both sides. + right_arm = ["*", "A", "*"] + left_arm = ["G", "G", "T"] + + energies = calculator.imperfection_energies(right_arm, left_arm) + + expected = energetics.CruciformEnergetics.mismatch_energy("GGT", "CAA") + + assert energies[1] == pytest.approx(expected) diff --git a/tests/test_ir_finder.py b/tests/test_ir_finder.py new file mode 100644 index 0000000..1e0390f --- /dev/null +++ b/tests/test_ir_finder.py @@ -0,0 +1,123 @@ +"""Unit tests for sist.ir_finder.""" + +from __future__ import annotations + +import shutil +from pathlib import Path + +import pytest + +from sist.ir_finder import IRFinder + + +def test_convert_to_one_line_keeps_existing_header(tmp_path: Path, monkeypatch) -> None: + monkeypatch.chdir(tmp_path) + Path("seq.fa").write_text(">seq.fa\nACGT\nTTTT\n", encoding="utf-8") + + output_path = IRFinder.convert_to_one_line(Path("seq.fa")) + + assert output_path == Path("one_line.seq.fa") + assert output_path.read_text(encoding="utf-8") == ">seq.fa\nACGTTTTT" + + +def test_convert_to_one_line_discards_first_line_without_header( + tmp_path: Path, monkeypatch +) -> None: + """A headerless first line is dropped entirely, not kept as sequence data.""" + monkeypatch.chdir(tmp_path) + Path("seq.fa").write_text("DROPPED_LINE\nACGT\nTTTT\n", encoding="utf-8") + + output_path = IRFinder.convert_to_one_line(Path("seq.fa")) + + assert output_path.read_text(encoding="utf-8") == ">one_line.seq.fa\nACGTTTTT" + + +def test_convert_to_one_line_empty_file(tmp_path: Path, monkeypatch) -> None: + monkeypatch.chdir(tmp_path) + Path("seq.fa").write_text("", encoding="utf-8") + + output_path = IRFinder.convert_to_one_line(Path("seq.fa")) + + assert output_path.read_text(encoding="utf-8") == "" + + +def test_read_sequence(tmp_path: Path) -> None: + one_line = tmp_path / "one_line.seq.fa" + one_line.write_text(">one_line.seq.fa\nABCDEFGH", encoding="utf-8") + + assert IRFinder.read_sequence(one_line) == "ABCDEFGH" + + +def test_sequence_length_and_shift_linear(tmp_path: Path) -> None: + one_line = tmp_path / "one_line.seq.fa" + one_line.write_text(">one_line.seq.fa\nABCDEFGH", encoding="utf-8") + + length, shifted = IRFinder.sequence_length_and_shift(one_line, circular=False) + + assert length == 8 + assert shifted == "" + + +def test_sequence_length_and_shift_circular(tmp_path: Path) -> None: + one_line = tmp_path / "one_line.seq.fa" + one_line.write_text(">one_line.seq.fa\nABCDEFGH", encoding="utf-8") + + length, shifted = IRFinder.sequence_length_and_shift(one_line, circular=True) + + assert length == 8 + assert shifted == "EFGHABCD" + + +def test_parse_loop_header_first_example() -> None: + # From example/one_line.pbr322.toy.fa...txt.html + line = " Indices: 2982--2992,2996--3006 Loop: 3 Score: 22" + + loop_length, loop_start, ir_start, ir_end, ir_length = IRFinder.parse_loop_header( + line + ) + + assert loop_length == 3 + assert loop_start == 2993 + assert ir_start == 2982 + assert ir_end == 3006 + assert ir_length == 25 + + +def test_parse_loop_header_second_example() -> None: + line = " Indices: 3089--3110,3115--3136 Loop: 4 Score: 44" + + loop_length, loop_start, ir_start, ir_end, ir_length = IRFinder.parse_loop_header( + line + ) + + assert loop_length == 4 + assert loop_start == 3111 + assert ir_start == 3089 + assert ir_end == 3136 + assert ir_length == 48 + + +def test_parse_arm_left() -> None: + line = " 2982 >> AAACCACCGCT >> 2992" + + assert IRFinder.parse_arm(line) == list("AAACCACCGCT") + + +def test_parse_arm_right() -> None: + line = " 3006 << *********** << 2996" + + assert IRFinder.parse_arm(line) == list("***********") + + +@pytest.mark.skipif(shutil.which("irf") is None, reason="requires the irf executable") +def test_compute_cruciform_energy_string_smoke(tmp_path: Path, monkeypatch) -> None: + monkeypatch.chdir(tmp_path) + Path("toy.fa").write_text( + ">toy.fa\nAAACCACCGCTTTTTTTTTTTTTTTTTTTTTTTTTGGCGGTGGTTT\n", + encoding="utf-8", + ) + + finder = IRFinder(temperature=310.0, shape="linear") + result = finder.compute_cruciform_energy_string("toy.fa") + + assert result.startswith("1,0,10000,|") diff --git a/tests/test_python_cli.py b/tests/test_python_cli.py new file mode 100644 index 0000000..776af8b --- /dev/null +++ b/tests/test_python_cli.py @@ -0,0 +1,104 @@ +"""Parity tests for the `sist` Python CLI against the v1.0.0 reference baselines. + +These prove the Python implementation is a safe drop-in replacement: it +reuses the exact same tests/reference/v1.0.0/ baselines the existing +regression suite in test_regression.py is already checked against. +""" + +from __future__ import annotations + +import pytest +from conftest import SistRun +from test_regression import ( + COMPETITION_REFERENCE, + REFERENCE_DIRECTORY, + assert_metadata_matches, + assert_metrics_match, + parse_competition_output, + parse_transition_output, +) + +pytestmark = pytest.mark.regression + + +def test_python_competition_command_succeeds(python_competition_run: SistRun) -> None: + """The Python CLI's competition calculation should complete successfully.""" + + process = python_competition_run.process + + assert process.returncode == 0, ( + f"{python_competition_run.name} command failed with exit code " + f"{process.returncode}\n\nstdout:\n{process.stdout}\n\nstderr:\n{process.stderr}" + ) + + +def test_python_transition_command_succeeds(python_transition_run: SistRun) -> None: + """Each Python CLI transition calculation should complete successfully.""" + + process = python_transition_run.process + + assert process.returncode == 0, ( + f"{python_transition_run.name} command failed with exit code " + f"{process.returncode}\n\nstdout:\n{process.stdout}\n\nstderr:\n{process.stderr}" + ) + + +def test_python_competition_matches_baseline(python_competition_run: SistRun) -> None: + """The Python CLI's competition output should match the 1.0.0 baseline.""" + + expected = parse_competition_output(COMPETITION_REFERENCE) + actual = parse_competition_output(python_competition_run.output_path) + + assert_metadata_matches( + name="python competition", expected=expected.metadata, actual=actual.metadata + ) + assert_metrics_match( + name="python competition", expected=expected.metrics, actual=actual.metrics + ) + + assert actual.profile.keys() == expected.profile.keys(), ( + "Python competition reported sequence positions changed" + ) + + for position, expected_row in expected.profile.items(): + actual_row = actual.profile[position] + assert actual_row == expected_row, ( + f"Python competition profile changed at position {position}: " + f"expected {expected_row!r}, actual {actual_row!r}" + ) + + +def test_python_transition_matches_baseline(python_transition_run: SistRun) -> None: + """Each Python CLI transition output should match its 1.0.0 baseline.""" + + # Reference files are named after the SIST_TRANSITIONS fixture names. + name_by_algorithm = {"M": "melting", "Z": "z-dna", "C": "cruciform"} + reference_output = ( + REFERENCE_DIRECTORY + / f"{name_by_algorithm[python_transition_run.algorithm]}.txt" + ) + + expected = parse_transition_output(reference_output) + actual = parse_transition_output(python_transition_run.output_path) + + assert_metadata_matches( + name=python_transition_run.name, + expected=expected.metadata, + actual=actual.metadata, + ) + assert_metrics_match( + name=python_transition_run.name, + expected=expected.metrics, + actual=actual.metrics, + ) + + assert actual.profile.keys() == expected.profile.keys(), ( + f"{python_transition_run.name} reported sequence positions changed" + ) + + for position, expected_row in expected.profile.items(): + actual_row = actual.profile[position] + assert actual_row == expected_row, ( + f"{python_transition_run.name} profile changed at position {position}: " + f"expected {expected_row!r}, actual {actual_row!r}" + ) From 9ee34f3e4e8a9ec2fef351e42fe46ef4cf59fbf9 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Wed, 23 Sep 2026 16:41:07 +0100 Subject: [PATCH 2/5] feat(sist): cut over to sist and remove Perl --- .github/workflows/release.yaml | 6 + IR_finder.pl | 401 ------------------------------- conda-recipe/build.sh | 22 +- conda-recipe/meta.yaml | 11 +- conda-recipe/run_test.sh | 3 +- docs/source/development.rst | 8 +- docs/source/installation.rst | 25 +- docs/source/source-usage.rst | 71 +++--- example/README.txt | 4 +- master.pl | 106 -------- tests/conftest.py | 64 ----- tests/reference/v1.0.0/README.md | 8 +- tests/test_python_cli.py | 104 -------- 13 files changed, 79 insertions(+), 754 deletions(-) delete mode 100755 IR_finder.pl delete mode 100755 master.pl delete mode 100644 tests/test_python_cli.py diff --git a/.github/workflows/release.yaml b/.github/workflows/release.yaml index 54b8a03..bf09da2 100644 --- a/.github/workflows/release.yaml +++ b/.github/workflows/release.yaml @@ -69,6 +69,11 @@ jobs: "s/{% set version = \".*\" %}/{% set version = \"${{ github.event.inputs.version }}\" %}/" \ conda-recipe/meta.yaml + # Update the sist Python package version + sed -i -E \ + "s/^(__version__\s*=\s*).*/\1\"${{ github.event.inputs.version }}\"/" \ + src/sist/__init__.py + # Update CITATION.cff version and date-released if [ -f CITATION.cff ]; then sed -i -E \ @@ -90,6 +95,7 @@ jobs: body: | Update version - Update the Conda recipe with new release + - Update the sist Python package version - Update CITATION.cff version & date-released - Auto-generated by [CI] committer: version-updater diff --git a/IR_finder.pl b/IR_finder.pl deleted file mode 100755 index adbf282..0000000 --- a/IR_finder.pl +++ /dev/null @@ -1,401 +0,0 @@ -#!/usr/bin/perl -#by Dina Zhabinskaya -#analyze a sequence with IR finder and output length of sequence, and position, length, and energy string used as input for SIDD-code for each IR -#this code can analyze circular plasmids in order to find possible IR's at the junction between the start and end position of the original sequence - -use strict; use warnings; - -my $usage = "\nusage: $0 \n\n". -"The script runs Inverted Repeat Finder (IRF), which can be downloaded at: http://tandem.bu.edu/irf/irf.download.html.\n". -"See below and description in reference (Nucleic Acids Res, 41(21), 9610) for IRF parameters used.\n". -"The output is a string of all possible IRs, with their start positions, lengths, and energies of extrusion.\n". -"Example of output: 1,0,10000,|2992,5,17.2738039658653,|2991,7,17.2738039658653,|2990,9,17.2738039658653,|...\n". -"Input requires temperature, shape of DNA, and sequence file:\n". -"temperature: numerical value in units of Kelvin\n". -"shape: circular or linear\n". -"sequence_file: provide a sequenece file\n". -"Sequence file (file_name) will be converted to a format appropriate for IRF: file one_line.file_name will be created.\n". -"An additional file circ.one_line.file_name will be created for circular sequences.\n"; - -if($#ARGV < 2) { - die $usage; -} - -my $code = "irf"; - -my $temp = $ARGV[0]; #temperature -my $shape = $ARGV[1]; #linear or circular -my $file = $ARGV[2]; #sequence fasta file -my @seq_name = split("/",$file); -my $in_file = $seq_name[-1]; - -use constant PI => 4*atan2(1,1); -my %bp = ('A','T','T','A','C','G','G','C'); - -#convert fasta file to fit IR_finder format -open(my $in,$in_file) || die "error opening $in_file\n"; -my $seq_file = "one_line.$in_file"; -open(my $out,">$seq_file") || die "error creating $seq_file\n"; - -my $first = 1; -while (my $line = <$in>) { - if ($first) { - if($line =~ m/>/) { - chomp($line); - print $out "$line\n"; - } - else { - print $out ">$seq_file\n"; - } - $first = 0; - next; - } - chomp($line); - print $out "$line"; -} -close($in); -close($out); - -#IR finder parameters -my $match = 2; -my $mismatch = 10; -my $delta = 10; -my $pm = 80; -my $pi = 10; -my $minscore = 20; -my $maxlength = 10000; -my $maxloop = 100; -my $minloop = 3; - -#energy parameters -my $salt = 0.01; #salt concentration -my $RT = 1.9872*$temp/1000.0; -my $Eb_c = 2.42; #imperfection energy -my $Ecr = 192.5-$temp*0.565-4*energy_melt("A")-2*2.44*$RT*log(4); #cruciform initiation energy - -#set and initialize parameters -my @imper; -my @R = (); my @L = (); -my ($s_loop,$l_loop,$length_IR,$start_IR,$end_IR,$El); -$s_loop=$l_loop=$length_IR=$start_IR=$end_IR=$El=0; -my ($Et,$Eimp); -$Eimp=$Et=0; -my $shorten_arm; -my @IR_array; - -my ($length_seq, $shift_seq); -my $circ_seq=0; -if ($shape eq "circular") { - ($length_seq,$shift_seq)=get_info($seq_file); - #create sequence file shifting the start position to the middle - $circ_seq = "circ.$seq_file"; - open(my $out,">$circ_seq") || die "error creating $circ_seq\n"; - print $out ">$circ_seq\n"; - print $out "$shift_seq"; - close($out); -} - -my $string_IR = "1,0,10000,|"; #format of string required for SIDD-code -$string_IR = get_result($seq_file,0); -if ($shape eq "circular") { - $string_IR = get_result($circ_seq,1); -} -print "$string_IR\n"; - -#Functions: - -#obtain a string containing IR start position, length, and energy as required for SIDD input -sub get_result { - my ($in_file,$shape) = @_; - my $skip = 0; - #PARAMETERS: Match: $match Mismatch: $mismatch Delta: $delta pm: $pm pi: $pi minscore:$minscore maxlength: $maxlength maxloop: $maxloop \n"; - system("$code $in_file $match $mismatch $delta $pm $pi $minscore $maxlength $maxloop > /dev/null"); - my $num_files = `ls -l $in_file.$match.$mismatch.$delta.$pm.$pi.$minscore.$maxlength.$maxloop.*.txt.html | wc -l`; - for(my $i = 1; $i <= $num_files; $i++) { - my $out_file = "$in_file.$match.$mismatch.$delta.$pm.$pi.$minscore.$maxlength.$maxloop.$i.txt.html"; - open(my $output,$out_file) || die "error opening $out_file\n"; - my $read = 0; - my @energies; - while (my $line = <$output>) { - $read = 0 if($line =~ m/frequent/); - last if ($line =~m/Done/); #no more IRs - if ($line =~ m/Loop:/) { - $skip = 0; - my @par_IR = get_IR($line); - $l_loop = $par_IR[0]; #Loop length - next if ($l_loop > $maxloop); #ignore loops of in this range - $s_loop = $par_IR[1]; #Start position of loop - $start_IR = $par_IR[2]; #start of IR - $end_IR = $par_IR[3]; #end of IR - $length_IR = $par_IR[4]; #length of IR - $IR_array[$s_loop][$l_loop][$start_IR]=$length_IR; #make array of positions and lengths - $read = 1; - } - if ($read) { - if ($line =~m/>>/ and $line !~m/LF/) { #Left arm - @L = get_arm($line,\@L); - } - if ($line =~m/<int($length_seq/2)) { - $s_loop-=int($length_seq/2)+$length_seq%2;; - $start_IR-=int($length_seq/2)+$length_seq%2;; - - } - else { - $s_loop+=int($length_seq/2); - $start_IR+=int($length_seq/2); - } - if ($IR_array[$s_loop][$l_loop][$start_IR]) { - if ($IR_array[$s_loop][$l_loop][$start_IR]==$length_IR) { - $skip = 1; - @L = (); #initialize arm arrays - @R = (); - } - } - } - next if ($skip); #for circular sequences - $El = energy_loop(@l_seq); #energy of loop - my $length_arm = @L - $shorten_arm; - @L = splice(@L,0,$length_arm); - @R = splice(@R,0,$length_arm); - @energies = energy_imperfections(\@R,\@L); - @energies = reverse(@energies); - $Et = $Ecr + $El; #energy of imperfections - my $s_IR = $s_loop; - my $l_IR = $l_loop; - my $j=0; - for(my $i = 0; $i < $length_arm; $i++) { - $Et += 2*$energies[$i]; - if ($L[$length_arm-1-$i] ne "-") { - $s_IR = $s_loop - ($j+1); - $l_IR = $l_loop + 2*($j+1); - $j++; - } - #string used as input in SIDD - $string_IR = join("",$string_IR,$s_IR,",",$l_IR,",",$Et,",","|"); - } - @L = (); #initialize arm arrays - @R = (); - } - } - } - close($output); - } - return $string_IR; -} - -#get length of sequence and shifted sequence when circular -sub get_info { - my $seq_file = $_[0]; - open(my $in,$seq_file) || die "error opening $seq_file\n"; - my $l_seq; - my $sf=0; - while (my $line = <$in>) { - next if($line =~ m/>/); #first line of fasta file - $l_seq = length($line); - if ($shape eq "circular") { - my $s1 = substr($line,0,int($l_seq/2)); - my $s2 = substr($line,int($l_seq/2)); - #this generates a string where the start position of the sequence is shifted to the middle - $sf = join("",$s2,$s1); - } - } - return ($l_seq,$sf); - close($in); -} - -#get loop bp -sub get_sequence { - my $in = $_[0]; - my $start = $_[1]; - my $l_loop = $_[2]; - my @seq; - my @loop; - while (my $line = <$in>) { - next if($line =~ m/>/); - @seq = split("",$line); - } - for(my $i = 0; $i < $l_loop; $i++) { - $loop[$i]=$seq[$start+$i]; - } - return @loop; -} - -# convert an IR arm into an array -sub get_arm { - my ($line,$A) = @_; - my @results = split(" ",$line); - my $arm = $results[2]; - my @left = split("",$arm); - my $l_arm = @left; - @left = splice(@left,0,$l_arm); - push(@$A,@left); - return @$A; -} - -#obtain the length of arms and loop -sub get_IR { - my $line = $_[0]; - my @results = split(" ",$line); - my $l_loop = $results[3]; #lenght of loop - my $arms = $results[1]; - my @rrl = split(",",$arms); - my @la = split("--",$rrl[0]); - my @ra = split("--",$rrl[1]); - my $start = $la[0]; #start position of IR - my $end = $ra[1]; #end of IR - my $s_loop = $la[1]+1; #start position of loop - my $l_IR = $end-$start+1; #length of IR - return ($l_loop,$s_loop,$start,$end,$l_IR); -} - -#calculate loop free energy -sub energy_loop { - my @seq = @_; - my $length = @seq; - my $El = 0; - for(my $i = 0; $i < $length; $i++) { - $El += energy_melt($seq[$i]); - } - $El+=2*2.44*$RT*log($length); #loop entropy term - return $El; -} - -#calculate imperfections energy -sub energy_imperfections { - my ($RA,$LA) =@_; - my @R = @$RA; - my @L = @$LA; - my @energies; - for(my $i = 0; $i < @L; $i++) { - $energies[$i] = 0; - } - my $Eb = 0; - my $E_miss = 0; - for(my $i = 0; $i < @L; $i++) { #bulges on right arm - if ($L[$i] eq "-") { - my $melt = energy_melt($R[$i]); - $energies[$i] = $Eb_c +$melt; - } - } - for(my $i = 0; $i < @R; $i++) { - if ($i == 0 and $R[$i] ne '*') { #test - die "Doesn't start with a star\n"; - } - if ($R[$i] eq "-") { #bulges on left arm - my $melt = energy_melt($L[$i]); - $energies[$i] = $Eb_c + $melt; - } - if ($R[$i] ne '*' and $R[$i] ne '-') { - next if($L[$i] eq '-'); #case accounted for above - if ($R[$i-1] ne '*' or $R[$i+1] ne '*') { #internal loop with 2bp open - $energies[$i] = 2*$Eb_c + energy_melt($R[$i]) + energy_melt($L[$i]); - } - # missmatch - if ($R[$i-1] eq '*' and $R[$i+1] eq '*' and $L[$i-1] ne '-' and $L[$i+1] ne '-') { - my $tri1 = join("",$L[$i-1],$L[$i],$L[$i+1]); - my $tri2 = join("",$bp{$L[$i-1]},$R[$i],$bp{$L[$i+1]}); - $energies[$i]= energy_missmatch($tri1,$tri2); - } - } - } - return @energies; -} - -#energy of melting each base pair -sub energy_melt { - my $base = $_[0]; - my $Ta = 354.65 + 16.6*log($salt)/log(10); - my $Tg = $Ta + 41; - my $Ea = 7.2464*(1-$temp/$Ta); - my $Eg = 9.0172*(1-$temp/$Tg); - if ($base eq 'a' or $base eq 't'or $base eq 'A' or $base eq 'T') { - return $Ea; - } - else { - return $Eg; - } -} - -#calculate the energy penalty of a missmatch -sub energy_missmatch { - my ($tri1,$tri2) = @_; - my @s1 = split("",$tri1); - my @s2 = split("",$tri2); - my $d1a = join('',$s1[0],$s1[1]); - my $d1b = join('',$s1[1],$s1[2]); - my %WC = energy_WC(); - my $e_WC=$WC{$d1a}+$WC{$d1b}; - my $d2a = join('',$s2[0],$s2[1]); - my $d1r = join('',$s1[2],$s1[1]); - my $d2r = join('',$s2[2],$s2[1]); - my $m1 = join('.',$d1a,$d2a); - my $m2 = join('.',$d2r,$d1r); - my %miss = energy_miss(); - my $e_miss=$miss{$m1}+$miss{$m2}; - my $E_miss = $e_miss-$e_WC; - return $E_miss; -} - -#Watson-Crick base-specific energies (Table 1 in paper) -sub energy_WC { - my %WC = ("AA",-1.0,"TT",-1.0,"AT",-0.88,"TA",-0.58, - "CA",-1.45,"AC",-1.45,"GT",-1.44,"TG",-1.44, - "CT",-1.28,"TC",-1.28,"GA",-1.3,"AG",-1.3, - "CG",-2.17,"GC",-2.24,"GG",-1.84,"CC",-1.84 - ); - return %WC; -} - -#Missmatch base-specific energies (Table 2 in paper) -sub energy_miss { - my %miss = ("GA.CA",0.17,"GA.CC",0.81,"GA.CG",-0.25, - "GC.CA",0.47,"GC.CC",0.79,"GC.CT",0.62, - "GG.CA",-0.52,"GG.CG",-1.11,"GG.CT",0.08, - "GT.CC",0.98,"GT.CG",-0.59,"GT.CT",0.45, - "CA.GA",0.43,"CA.GC",0.75,"CA.GG",0.03, - "CC.GA",0.79,"CC.GC",0.70,"CC.GT",0.62, - "CG.GA",0.11,"CG.GG",-0.11,"CG.GT",-0.47, - "CT.GC",0.40,"CT.GG",-0.32,"CT.GT",-0.12, - "AA.TA",0.61,"AA.TC",0.88,"AA.TG",0.14, - "AC.TA",0.77,"AC.TC",1.33,"AC.TT",0.64, - "AG.TA",0.02,"AG.TG",-0.13,"AG.TT",0.71, - "AT.TC",0.73,"AT.TG",0.07,"AT.TT",0.69, - "TA.AA",0.69,"TA.AC",0.92,"TA.AG",0.42, - "TC.AA",1.33,"TC.AC",1.05,"TC.AT",0.97, - "TG.AA",0.74,"TG.AG",0.44,"TG.AT",0.43, - "TT.AC",0.75,"TT.AG",0.34,"TT.AT",0.68, - ); - return %miss; -} diff --git a/conda-recipe/build.sh b/conda-recipe/build.sh index 1ca8168..80b758b 100755 --- a/conda-recipe/build.sh +++ b/conda-recipe/build.sh @@ -7,19 +7,9 @@ make -C src/trans_three make -C src/trans_compete clean make -C src/trans_compete -install -d "${PREFIX}/bin" -install -d "${PREFIX}/libexec/sist" install -d "${PREFIX}/libexec/sist/src/trans_three" install -d "${PREFIX}/libexec/sist/src/trans_compete" -install -m 755 \ - master.pl \ - "${PREFIX}/libexec/sist/master.pl" - -install -m 755 \ - IR_finder.pl \ - "${PREFIX}/libexec/sist/IR_finder.pl" - install -m 755 \ src/trans_three/qsidd \ "${PREFIX}/libexec/sist/src/trans_three/qsidd" @@ -28,14 +18,4 @@ install -m 755 \ src/trans_compete/qsidd \ "${PREFIX}/libexec/sist/src/trans_compete/qsidd" -cat > "${PREFIX}/bin/sist" << 'EOF' -#!/usr/bin/env bash - -PREFIX="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" - -exec "${PREFIX}/bin/perl" \ - "${PREFIX}/libexec/sist/master.pl" \ - "$@" -EOF - -chmod 755 "${PREFIX}/bin/sist" \ No newline at end of file +"${PYTHON}" -m pip install . --no-deps --no-build-isolation -vv diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index c989395..d54831d 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -16,8 +16,15 @@ requirements: - {{ compiler('cxx') }} - make + host: + - python >=3.12 + - pip + - setuptools >=77,<84 + run: - - perl + - python >=3.12 + - biopython >=1.85,<2.0 + - beautifulsoup4 >=4.12,<5.0 - irf >=3.08,<3.09 test: @@ -40,4 +47,4 @@ about: SIST calculates stress-induced structural transition probabilities in superhelical DNA, including strand separation, Z-DNA formation, cruciform extrusion, and competition between these transitions. - dev_url: https://github.com/CCPBioSim/SIST \ No newline at end of file + dev_url: https://github.com/CCPBioSim/SIST diff --git a/conda-recipe/run_test.sh b/conda-recipe/run_test.sh index f6bace2..62d6569 100755 --- a/conda-recipe/run_test.sh +++ b/conda-recipe/run_test.sh @@ -3,8 +3,7 @@ set -euo pipefail command -v sist command -v irf -command -v perl command -v python command -v pytest -python -m pytest tests -vv \ No newline at end of file +python -m pytest tests -vv diff --git a/docs/source/development.rst b/docs/source/development.rst index 4729749..00d015e 100644 --- a/docs/source/development.rst +++ b/docs/source/development.rst @@ -30,7 +30,6 @@ also requires: * GNU Make * a C++ compiler -* Perl * IRF 3.08 on ``PATH`` Install the Python testing dependencies: @@ -86,18 +85,17 @@ Build requirements Runtime requirements ~~~~~~~~~~~~~~~~~~~~ -* Perl +* Python, Biopython, and Beautiful Soup * IRF ``>=3.08,<3.09`` * compiler runtime libraries resolved by Conda Package-test requirements ~~~~~~~~~~~~~~~~~~~~~~~~~ -* Python * pytest -Python and pytest are package-test dependencies only; they are not required for -normal use of the installed SIST package. +pytest is a package-test dependency only; it is not required for normal use +of the installed SIST package. Build and test the Conda package -------------------------------- diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 784dc59..06e3d3b 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -13,15 +13,15 @@ The recommended way to install SIST is with Conda: The package installs the ``sist`` command together with the runtime dependencies required by SIST, including: -* Perl +* Python, Biopython, and Beautiful Soup * Inverted Repeats Finder (IRF) * the required C++ runtime libraries SIST 1.0.0 is validated with IRF 3.08, and the Conda package constrains the runtime dependency to ``>=3.08,<3.09``. -Python and pytest are used for testing and are not required in a normal SIST -runtime environment. +pytest is used for testing and is not required in a normal SIST runtime +environment. Verify the installation ----------------------- @@ -31,7 +31,6 @@ Confirm that the installed command and runtime dependencies are available: .. code-block:: bash command -v sist - command -v perl command -v irf The commands should resolve inside the active Conda environment. @@ -52,7 +51,7 @@ A source build requires: * a C++ compiler * GNU Make -* Perl +* Python 3.12 or later * IRF 3.08 available as ``irf`` on ``PATH`` Build both C++ components from the repository root: @@ -62,11 +61,25 @@ Build both C++ components from the repository root: make -C src/trans_three make -C src/trans_compete +Install the ``sist`` Python package: + +.. code-block:: bash + + python -m pip install . + +``sist`` looks for the compiled ``qsidd`` binaries at a Conda install layout +by default, so point it at the binaries just built from source instead: + +.. code-block:: bash + + export SIST_TRANS_THREE_BIN="$(pwd)/src/trans_three/qsidd" + export SIST_TRANS_COMPETE_BIN="$(pwd)/src/trans_compete/qsidd" + The source-tree pipeline can then be run with: .. code-block:: bash - perl master.pl -a M -f sequence.fa + sist -a M -f sequence.fa For cruciform and competition calculations, IRF must be available on ``PATH``. diff --git a/docs/source/source-usage.rst b/docs/source/source-usage.rst index 2385cab..550423f 100644 --- a/docs/source/source-usage.rst +++ b/docs/source/source-usage.rst @@ -1,8 +1,8 @@ Source Usage ============ -SIST consists of a Perl pipeline, an IRF integration script, and two C++ -implementations of the transition calculations. +SIST consists of a Python package, ``sist``, and two C++ implementations of +the transition calculations. Normal installed use should go through the ``sist`` command. The interfaces described on this page are useful when building, inspecting, or running the @@ -11,13 +11,11 @@ source tree directly. Source components ----------------- -``master.pl`` - Pipeline used to run the supported SIST calculations. - -``IR_finder.pl`` - Processes Inverted Repeats Finder (IRF) output and produces the inverted - repeat information required for cruciform calculations, including start - positions, possible extrusion lengths, and cruciform formation energies. +``src/sist/`` + The ``sist`` Python package: the command-line interface, orchestration + that dispatches to the ``qsidd`` binaries, and the Inverted Repeats Finder + (IRF) integration and energy calculations required for cruciform + calculations. ``src/trans_three/`` C++ implementation for analysing strand separation, Z-DNA, and cruciform @@ -27,20 +25,21 @@ Source components C++ implementation for analysing competition between strand separation, Z-DNA, and cruciform extrusion. -Running ``master.pl`` ---------------------- +Running ``sist`` from source +----------------------------- -After building both C++ components, run the source-tree pipeline with Perl: +After building both C++ components and installing the ``sist`` package (see +:doc:`installation`), run the source-tree pipeline with: .. code-block:: bash - perl master.pl -f -a [options] + sist -f -a [options] For example: .. code-block:: bash - perl master.pl -a M -f sequence.fa + sist -a M -f sequence.fa The available algorithm types are: @@ -49,16 +48,16 @@ The available algorithm types are: * ``-a C``: cruciform transition only * ``-a A``: competition between melting, Z-DNA, and cruciform transitions -Running ``perl master.pl`` without the required arguments displays the -available command-line options. +Running ``sist`` without the required arguments displays the available +command-line options. IRF --- Cruciform and competition calculations require Inverted Repeats Finder. -``IR_finder.pl`` invokes ``irf`` from ``PATH``. For source builds, install a -compatible IRF 3.08 executable and ensure that: +``sist``'s IR-finder component invokes ``irf`` from ``PATH``. For source +builds, install a compatible IRF 3.08 executable and ensure that: .. code-block:: bash @@ -91,13 +90,17 @@ Cruciform and competition component workflow -------------------------------------------- When running the components directly, cruciform and competition calculations -require the output produced by ``IR_finder.pl``. - -Run ``IR_finder.pl`` first: +require the inverted-repeat energy string that ``sist`` normally computes +internally via :mod:`sist.ir_finder` before invoking ``qsidd -X``. To produce +it directly from Python: .. code-block:: bash - perl IR_finder.pl temperature shape sequence_file + python -c " + from sist.ir_finder import IRFinder + finder = IRFinder(temperature=310.0, shape='linear') + print(finder.compute_cruciform_energy_string('sequence_file')) + " For a cruciform calculation using ``src/trans_three``: @@ -111,17 +114,17 @@ For a competition calculation using ``src/trans_compete``: ./qsidd -X "string" -f sequence_file -Here, ``string`` is the output produced by ``IR_finder.pl``. +Here, ``string`` is the output produced above. -``master.pl`` coordinates this workflow automatically and is normally the +``sist`` coordinates this workflow automatically and is normally the preferred source-tree entry point. Working directory ----------------- -For cruciform and competition calculations, ``IR_finder.pl`` uses the basename -of the input sequence. The sequence file should therefore be present in the -current working directory when these calculations are run. +For cruciform and competition calculations, :mod:`sist.ir_finder` uses the +basename of the input sequence. The sequence file should therefore be present +in the current working directory when these calculations are run. Example calculation ------------------- @@ -129,19 +132,7 @@ Example calculation The repository contains an example competition calculation based on ``pbr322.toy.fa``. -The source-tree command is: - -.. code-block:: bash - - perl master.pl \ - -f pbr322.toy.fa \ - -a A \ - -o pbr322.toy.compete.txt \ - -b \ - -p \ - -r - -The equivalent installed command is: +The command is: .. code-block:: bash diff --git a/example/README.txt b/example/README.txt index 35da5f2..bd83ab1 100644 --- a/example/README.txt +++ b/example/README.txt @@ -1,14 +1,14 @@ Example of output generated when the competion algorithm analyzes the pbr322.toy.fa sequence. These results are described in the Bioinformatics citation. To execute type the following on the command line: -perl master.pl -f pbr322.toy.fa -a A -o pbr322.toy.compete.txt -b -p -r +sist -f pbr322.toy.fa -a A -o pbr322.toy.compete.txt -b -p -r Input files: 1. pbr322.toy.fa Output files: 1. one_line.pbr322.toy.fa : sequences file format required by Inverted Repeat Finder (IRF) -2. one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html: IRF output file +2. one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.html: IRF output file 3. one_line.pbr322.toy.fa.2.10.10.80.10.20.10000.100.1.txt.html: IRF output file 4. pbr322.toy.compete.txt : output file diff --git a/master.pl b/master.pl deleted file mode 100755 index 833b35b..0000000 --- a/master.pl +++ /dev/null @@ -1,106 +0,0 @@ -#!/usr/bin/perl -#by Dina Zhabinskaya -#This code is called by compete_many.thread.pl to run IR_finder and SIDD on single-transition code -#at input temperature for an input sequence file - -use strict; use warnings; -use Getopt::Long; -use FindBin qw($RealBin); -my ($file,$trans, $out_file); -my $temp = 310; -my $sig = 0.06; -my $theta = 12; -my $salt = 0.01; -my $shape = "linear"; -my $c = ""; -my $b = ""; -my $p = ""; -my $r = ""; -my $n = ""; - -my $usage = "\nusage: $0 -f -a (choose algorithm_type: M, Z, C, or A) [options]\n\n". -"Input requires a sequence file and algorithm type.\n". -"Script analyzes superhelically induced structural transition probabilities for each base pair.\n". -"Algorithm types: M (melting), Z (Z-DNA), C (cruciforms), and A (competition between all three).\n". -"Sequence file will be converted to the format required by the algorithm.\n". -"Code in directory src/trans_three/ will handle M, Z, and C algorithm types.\n". -"Code in directory src/trans_compete/ will handle A algorithm type.\n". -"For algorithm type options -a C and -a A user will need an Inverted Repeat Finder (IRF) executable compatible with user's operating system.\n". -"IRF download page: http://tandem.bu.edu/irf/irf.download.html.\n". -"Selected output will be printed to the screen.\n\n". -"Options:\n". -"-f Required: specify sequence file \n". -"-a Required: specify algorithm type: M, Z, C, or A\n". -"-t Optional: set temperature (default 310)\n". -"-s Optional: set superhelical density (default -0.06)\n". -"-i Optional: set ionic strength (default 0.01)\n". -"-th Optional: set energy threshold (default 12), -t 10 is recommended for -a A (competition algorithm)\n". -"-c Optional: flag to set molecular type to circular (default linear)\n". -"-n Optional: flag to set melting energetics to nearest neighbor (default copolymeric)\n". -"-b Optional: flag to print base pair for each position (default null)\n". -"-p Optional: flag to print algorithm parameters (default null)\n". -"-r Optional: flag to print ensemble average results (default null)\n". -"-o Optional: specify output file name\n"; - - -if($#ARGV < 2) { - die $usage; -} - -GetOptions ( -"f=s" => \$file, # sequence file -"a=s" => \$trans, # algorithm -"T=s" => \$temp, # temperature -"s=s" => \$sig, # superhelical density -"i=s" => \$salt, # ionic concentration -"th=s" => \$theta, # energy threshold -"c" => \$c, # molecular shape -"b" => \$b, # print sequence -"n" => \$n, # nearest neighbor energetics for melting -"p" => \$p, # print parameters -"r" => \$r, # print average results -"o=s" => \$out_file) # output file -or die("Error in command line arguments\n$usage\n"); - -$b="-b" if($b); -$p="-p" if($p); -$r="-r" if($r); -$n="-n" if($n); -if($c) { - $c="-c"; - $shape="circular"; -} - - -my @name = split("/",$file); -my $single_exe = "$RealBin/src/trans_three/qsidd"; -my $compete_exe = "$RealBin/src/trans_compete/qsidd"; - -my $output_IR; -if ($trans eq "M") { - #run single transition algorithm for melting - $output_IR = `$single_exe $b $p $r $n $c -T $temp -s $sig -i $salt -t $theta -f $file`; -} -if ($trans eq "Z") { - #run single transition algorithm for Z-DNA - $output_IR = `$single_exe $b $p $r $c -T $temp -s $sig -i $salt -t $theta -Z -f $file`; -} -if ($trans eq "C" or $trans eq "A") { - my $code_IR = "$RealBin/IR_finder.pl"; - my $IR_results = `perl $code_IR $temp $shape $file`; #run IR_finder - if($trans eq "C") { - #run single transition algorithm for cruciforms - $output_IR = `$single_exe $b $p $r $c -T $temp -s $sig -i $salt -t $theta -C -X "$IR_results" -f $file`; - } - if ($trans eq "A") { - #run competition algorith: melting, Z-DNA, and cruciforms competing - $output_IR = `$compete_exe $b $p $r $n $c -T $temp -s $sig -i $salt -t $theta -X "$IR_results" -f $file`; - } -} -if ($out_file) { - open(my $out,">$out_file") || die "error creating $out_file\n"; - print $out $output_IR; -} -else { - print "$output_IR"; -} diff --git a/tests/conftest.py b/tests/conftest.py index 0c7ce5f..f6f1baa 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -133,34 +133,6 @@ def sist_command( """ Return the SIST command under test. - During conda-build testing, use the installed package. Otherwise build - and test the maintained source tree. - """ - - if os.environ.get("CONDA_BUILD_STATE") == "TEST": - executable = shutil.which("sist") - - if executable is None: - pytest.fail( - "The installed sist command was not found in PATH", - pytrace=False, - ) - - return [executable] - - return [ - "perl", - str(built_sist_copy / "master.pl"), - ] - - -@pytest.fixture(scope="session") -def python_sist_command( - built_sist_copy: Path, -) -> list[str]: - """ - Return the new `sist` Python CLI command under test. - During conda-build testing, use the installed console script. Otherwise run `python -m sist` against the maintained source tree, with its qsidd binaries resolved via env vars pointing at the freshly built copy. @@ -263,39 +235,3 @@ def transition_run( name=name, algorithm=algorithm, ) - - -@pytest.fixture(scope="session") -def python_competition_run( - python_sist_command: list[str], - tmp_path_factory: pytest.TempPathFactory, -) -> SistRun: - """Run the new sist Python CLI's competition calculation once.""" - - return run_sist_calculation( - python_sist_command, - tmp_path_factory, - name="python-competition", - algorithm="A", - ) - - -@pytest.fixture( - scope="session", - params=SIST_TRANSITIONS, -) -def python_transition_run( - request: pytest.FixtureRequest, - python_sist_command: list[str], - tmp_path_factory: pytest.TempPathFactory, -) -> SistRun: - """Run each sist Python CLI transition calculation once.""" - - name, algorithm = request.param - - return run_sist_calculation( - python_sist_command, - tmp_path_factory, - name=f"python-{name}", - algorithm=algorithm, - ) diff --git a/tests/reference/v1.0.0/README.md b/tests/reference/v1.0.0/README.md index c8b7011..73c2e3b 100644 --- a/tests/reference/v1.0.0/README.md +++ b/tests/reference/v1.0.0/README.md @@ -20,6 +20,12 @@ Repeated calculations produced identical deterministic output. The Run time valu ## Commands +These were the commands used to generate the reference outputs below, run +against the `master.pl` pipeline that shipped in SIST 1.0.0 and has since +been replaced by the `sist` Python package. They are kept here as a record +of provenance; the baselines themselves remain the source of truth +regardless of which implementation produced them. + ### Melting perl master.pl \ @@ -103,4 +109,4 @@ Reference outputs must not be updated solely to make a failing regression test p Any numerical difference from these baselines should be investigated and understood before a reference is changed. -Intentional scientific changes should document why the established SIST 1.0.0 behaviour has changed. \ No newline at end of file +Intentional scientific changes should document why the established SIST 1.0.0 behaviour has changed. diff --git a/tests/test_python_cli.py b/tests/test_python_cli.py deleted file mode 100644 index 776af8b..0000000 --- a/tests/test_python_cli.py +++ /dev/null @@ -1,104 +0,0 @@ -"""Parity tests for the `sist` Python CLI against the v1.0.0 reference baselines. - -These prove the Python implementation is a safe drop-in replacement: it -reuses the exact same tests/reference/v1.0.0/ baselines the existing -regression suite in test_regression.py is already checked against. -""" - -from __future__ import annotations - -import pytest -from conftest import SistRun -from test_regression import ( - COMPETITION_REFERENCE, - REFERENCE_DIRECTORY, - assert_metadata_matches, - assert_metrics_match, - parse_competition_output, - parse_transition_output, -) - -pytestmark = pytest.mark.regression - - -def test_python_competition_command_succeeds(python_competition_run: SistRun) -> None: - """The Python CLI's competition calculation should complete successfully.""" - - process = python_competition_run.process - - assert process.returncode == 0, ( - f"{python_competition_run.name} command failed with exit code " - f"{process.returncode}\n\nstdout:\n{process.stdout}\n\nstderr:\n{process.stderr}" - ) - - -def test_python_transition_command_succeeds(python_transition_run: SistRun) -> None: - """Each Python CLI transition calculation should complete successfully.""" - - process = python_transition_run.process - - assert process.returncode == 0, ( - f"{python_transition_run.name} command failed with exit code " - f"{process.returncode}\n\nstdout:\n{process.stdout}\n\nstderr:\n{process.stderr}" - ) - - -def test_python_competition_matches_baseline(python_competition_run: SistRun) -> None: - """The Python CLI's competition output should match the 1.0.0 baseline.""" - - expected = parse_competition_output(COMPETITION_REFERENCE) - actual = parse_competition_output(python_competition_run.output_path) - - assert_metadata_matches( - name="python competition", expected=expected.metadata, actual=actual.metadata - ) - assert_metrics_match( - name="python competition", expected=expected.metrics, actual=actual.metrics - ) - - assert actual.profile.keys() == expected.profile.keys(), ( - "Python competition reported sequence positions changed" - ) - - for position, expected_row in expected.profile.items(): - actual_row = actual.profile[position] - assert actual_row == expected_row, ( - f"Python competition profile changed at position {position}: " - f"expected {expected_row!r}, actual {actual_row!r}" - ) - - -def test_python_transition_matches_baseline(python_transition_run: SistRun) -> None: - """Each Python CLI transition output should match its 1.0.0 baseline.""" - - # Reference files are named after the SIST_TRANSITIONS fixture names. - name_by_algorithm = {"M": "melting", "Z": "z-dna", "C": "cruciform"} - reference_output = ( - REFERENCE_DIRECTORY - / f"{name_by_algorithm[python_transition_run.algorithm]}.txt" - ) - - expected = parse_transition_output(reference_output) - actual = parse_transition_output(python_transition_run.output_path) - - assert_metadata_matches( - name=python_transition_run.name, - expected=expected.metadata, - actual=actual.metadata, - ) - assert_metrics_match( - name=python_transition_run.name, - expected=expected.metrics, - actual=actual.metrics, - ) - - assert actual.profile.keys() == expected.profile.keys(), ( - f"{python_transition_run.name} reported sequence positions changed" - ) - - for position, expected_row in expected.profile.items(): - actual_row = actual.profile[position] - assert actual_row == expected_row, ( - f"{python_transition_run.name} profile changed at position {position}: " - f"expected {expected_row!r}, actual {actual_row!r}" - ) From 57b32e9c19d531bb82341e1c96e719a83dd930f1 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Thu, 24 Sep 2026 10:40:49 +0100 Subject: [PATCH 3/5] fix(sist): resolve built_sist_copy lazily to skip C++ rebuild during conda-build tests --- tests/conftest.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/conftest.py b/tests/conftest.py index f6f1baa..09632e9 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -128,14 +128,14 @@ def built_sist_copy( @pytest.fixture(scope="session") def sist_command( - built_sist_copy: Path, + request: pytest.FixtureRequest, ) -> list[str]: """ Return the SIST command under test. During conda-build testing, use the installed console script. Otherwise run `python -m sist` against the maintained source tree, with its qsidd - binaries resolved via env vars pointing at the freshly built copy. + binaries resolved via env vars pointing at a freshly built copy. """ if os.environ.get("CONDA_BUILD_STATE") == "TEST": @@ -149,6 +149,8 @@ def sist_command( return [executable] + built_sist_copy: Path = request.getfixturevalue("built_sist_copy") + os.environ["SIST_TRANS_THREE_BIN"] = str( built_sist_copy / "src" / "trans_three" / "qsidd" ) From 11c1854cccd8b0d405462f63cd5b7ffa0ce56ec3 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Thu, 24 Sep 2026 12:17:10 +0100 Subject: [PATCH 4/5] chore: remove unneeded comments --- src/sist/ir_finder.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/sist/ir_finder.py b/src/sist/ir_finder.py index 0f947fa..ac919c4 100644 --- a/src/sist/ir_finder.py +++ b/src/sist/ir_finder.py @@ -283,10 +283,7 @@ def _extend_loop_to_minimum( loop_bases.append(left_value) loop_bases.append(right_value) else: - # Both arms may be exhausted here, giving None; melting_energy() - # already treats an unrecognised/missing base as GC, so the - # None is left to flow through unchanged. - loop_bases.append(index_or_none(right_arm, index)) # type: ignore[arg-type] + loop_bases.append(index_or_none(right_arm, index)) shorten_arm += 1 loop_start -= 1 From 457c27e94326c9308985f4d7087df0c2150f46cf Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Thu, 24 Sep 2026 15:19:30 +0100 Subject: [PATCH 5/5] feat(sist): add deprecated master.pl/IR_finder.pl entrypoint aliases --- docs/source/installation.rst | 10 ++++ docs/source/source-usage.rst | 29 +++++++++++- pyproject.toml | 4 ++ src/sist/cli.py | 15 ++++++ src/sist/ir_finder.py | 2 +- src/sist/ir_finder_cli.py | 54 +++++++++++++++++++++ tests/test_legacy_entrypoints.py | 81 ++++++++++++++++++++++++++++++++ 7 files changed, 192 insertions(+), 3 deletions(-) create mode 100644 src/sist/ir_finder_cli.py create mode 100644 tests/test_legacy_entrypoints.py diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 06e3d3b..c013cec 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -86,3 +86,13 @@ For cruciform and competition calculations, IRF must be available on ``PATH``. The source build uses the same calculation modes and command-line parameters as the installed ``sist`` command. The Source Usage page describes the individual source components and direct component workflow in more detail. + +Deprecated ``master.pl``/``IR_finder.pl`` aliases +-------------------------------------------------- + +For scripted workflows that still invoke the old Perl script names directly, +``master.pl`` and ``IR_finder.pl`` are also installed as console-script +aliases, backed by the same Python implementation as ``sist``. Both print a +deprecation warning and will be removed in the next release; switch scripts +over to ``sist`` (and the Python API, for direct use of the IR-finder +component) in the meantime. diff --git a/docs/source/source-usage.rst b/docs/source/source-usage.rst index 550423f..94e68ad 100644 --- a/docs/source/source-usage.rst +++ b/docs/source/source-usage.rst @@ -51,6 +51,20 @@ The available algorithm types are: Running ``sist`` without the required arguments displays the available command-line options. +Deprecated ``master.pl`` alias +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For scripted workflows that still invoke the old Perl script name directly, +``master.pl`` is also installed as a console-script alias for ``sist``, +accepting the exact same arguments: + +.. code-block:: bash + + master.pl -a M -f sequence.fa + +It prints a deprecation warning to stderr and will be removed in the next +release; switch scripts over to ``sist`` in the meantime. + IRF --- @@ -91,8 +105,19 @@ Cruciform and competition component workflow When running the components directly, cruciform and competition calculations require the inverted-repeat energy string that ``sist`` normally computes -internally via :mod:`sist.ir_finder` before invoking ``qsidd -X``. To produce -it directly from Python: +internally via :mod:`sist.ir_finder` before invoking ``qsidd -X``. + +For scripted workflows that still invoke the old Perl script name directly, +this is also available as a deprecated ``IR_finder.pl`` console-script alias, +accepting the same positional arguments as the original script and printing +the same string to stdout: + +.. code-block:: bash + + IR_finder.pl 310 linear sequence_file + +It prints a deprecation warning to stderr and will be removed in the next +release. To produce the string directly from Python instead: .. code-block:: bash diff --git a/pyproject.toml b/pyproject.toml index 23de95a..f2fab85 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,10 @@ Repository = "https://github.com/CCPBioSim/SIST" [project.scripts] sist = "sist.cli:main" +# Deprecated aliases for scripted workflows that still invoke the old Perl +# script names directly. +"master.pl" = "sist.cli:main" +"IR_finder.pl" = "sist.ir_finder_cli:main" [project.optional-dependencies] testing = [ diff --git a/src/sist/cli.py b/src/sist/cli.py index d1e4e97..6fc6f4b 100644 --- a/src/sist/cli.py +++ b/src/sist/cli.py @@ -5,17 +5,25 @@ 2) Validating those arguments. 3) Constructing a SistRunner and running it. 4) Handling fatal errors with a non-zero exit code. + +It is also installed under the deprecated console-script alias ``master.pl`` +(see ``[project.scripts]`` in ``pyproject.toml``) for scripted workflows that +still invoke it by that name. That alias will be removed in the next release. """ from __future__ import annotations import logging +import sys +from pathlib import Path from sist.argspec import SistArgumentParser from sist.runner import SistRunner logger = logging.getLogger(__name__) +DEPRECATED_ALIAS = "master.pl" + def main(argv: list[str] | None = None) -> None: """Parse CLI arguments and run a SIST calculation. @@ -28,6 +36,13 @@ def main(argv: list[str] | None = None) -> None: """ logging.basicConfig(level=logging.WARNING, format="%(levelname)s: %(message)s") + if Path(sys.argv[0]).name == DEPRECATED_ALIAS: + print( + f"WARNING: '{DEPRECATED_ALIAS}' is a deprecated alias for 'sist' " + "and will be removed in the next release. Use 'sist' instead.", + file=sys.stderr, + ) + argument_parser = SistArgumentParser() parser = argument_parser.build_parser() args = parser.parse_args(argv) diff --git a/src/sist/ir_finder.py b/src/sist/ir_finder.py index ac919c4..236add7 100644 --- a/src/sist/ir_finder.py +++ b/src/sist/ir_finder.py @@ -283,7 +283,7 @@ def _extend_loop_to_minimum( loop_bases.append(left_value) loop_bases.append(right_value) else: - loop_bases.append(index_or_none(right_arm, index)) + loop_bases.append(index_or_none(right_arm, index)) # type: ignore[arg-type] shorten_arm += 1 loop_start -= 1 diff --git a/src/sist/ir_finder_cli.py b/src/sist/ir_finder_cli.py new file mode 100644 index 0000000..a5886d6 --- /dev/null +++ b/src/sist/ir_finder_cli.py @@ -0,0 +1,54 @@ +"""Deprecated console-script entry point for the old `IR_finder.pl` interface. + +`IR_finder.pl` took positional arguments (`` +``) rather than flags, and printed the resulting IR string to +stdout. This module reproduces that exact interface on top of +:class:`sist.ir_finder.IRFinder` for scripted workflows that still invoke it +by that name (see ``[project.scripts]`` in ``pyproject.toml``). It will be +removed in the next release; use the ``sist`` command instead. +""" + +from __future__ import annotations + +import logging +import sys + +from sist.ir_finder import IRFinder + +logger = logging.getLogger(__name__) + +USAGE = "usage: IR_finder.pl " + + +def main(argv: list[str] | None = None) -> None: + """Parse `IR_finder.pl`-style positional arguments and print the IR string. + + Args: + argv: Argument list to parse, or None to use `sys.argv[1:]`. + + Raises: + SystemExit: Exits with status code 1 if the arguments are missing, + or on any unhandled exception. + """ + logging.basicConfig(level=logging.WARNING, format="%(levelname)s: %(message)s") + + print( + "WARNING: 'IR_finder.pl' is a deprecated alias for 'sist' and will " + "be removed in the next release. Use 'sist' instead.", + file=sys.stderr, + ) + + args = sys.argv[1:] if argv is None else argv + + if len(args) < 3: + print(USAGE, file=sys.stderr) + raise SystemExit(1) + + temperature, shape, sequence_file = args[0], args[1], args[2] + + try: + finder = IRFinder(temperature=float(temperature), shape=shape) + print(finder.compute_cruciform_energy_string(sequence_file)) + except Exception: + logger.exception("Fatal error during IR_finder calculation") + raise SystemExit(1) from None diff --git a/tests/test_legacy_entrypoints.py b/tests/test_legacy_entrypoints.py new file mode 100644 index 0000000..4a9f141 --- /dev/null +++ b/tests/test_legacy_entrypoints.py @@ -0,0 +1,81 @@ +"""Tests for the deprecated `master.pl`/`IR_finder.pl` console-script aliases. + +These exist for scripted workflows that still invoke the old Perl script +names directly; both aliases will be removed in the next release. +""" + +from __future__ import annotations + +import os +import shutil +import subprocess +from pathlib import Path + +import pytest +from conftest import REPOSITORY_ROOT + +pytestmark = pytest.mark.regression + + +def test_master_pl_alias_matches_sist_output( + request: pytest.FixtureRequest, + tmp_path_factory: pytest.TempPathFactory, +) -> None: + """The deprecated `master.pl` alias should behave exactly like `sist`. + + `built_sist_copy` is looked up lazily (rather than taken as a normal + fixture parameter) so conda-build testing never triggers it: that build + only has `tests/` and `pyproject.toml` available, not the C++ sources. + """ + + environment = os.environ.copy() + + if os.environ.get("CONDA_BUILD_STATE") != "TEST": + built_sist_copy: Path = request.getfixturevalue("built_sist_copy") + environment["SIST_TRANS_THREE_BIN"] = str( + built_sist_copy / "src" / "trans_three" / "qsidd" + ) + environment["SIST_TRANS_COMPETE_BIN"] = str( + built_sist_copy / "src" / "trans_compete" / "qsidd" + ) + + runtime_directory = tmp_path_factory.mktemp("master-pl-alias") + input_path = runtime_directory / "pbr322.toy.fa" + shutil.copy2(REPOSITORY_ROOT / "tests" / "data" / "pbr322.toy.fa", input_path) + + result = subprocess.run( + ["master.pl", "-f", input_path.name, "-a", "M", "-b", "-p", "-r"], + cwd=runtime_directory, + env=environment, + capture_output=True, + text=True, + check=False, + ) + + assert result.returncode == 0, result.stderr + assert "deprecated" in result.stderr.lower() + assert "Sequence Length" in result.stdout + + +@pytest.mark.skipif(shutil.which("irf") is None, reason="requires the irf executable") +def test_ir_finder_pl_alias_keeps_original_interface( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """The deprecated `IR_finder.pl` alias should keep its positional-arg interface.""" + + monkeypatch.chdir(tmp_path) + Path("toy.fa").write_text( + ">toy.fa\nAAACCACCGCTTTTTTTTTTTTTTTTTTTTTTTTTGGCGGTGGTTT\n", + encoding="utf-8", + ) + + result = subprocess.run( + ["IR_finder.pl", "310.0", "linear", "toy.fa"], + capture_output=True, + text=True, + check=False, + ) + + assert result.returncode == 0, result.stderr + assert "deprecated" in result.stderr.lower() + assert result.stdout.strip().startswith("1,0,10000,|")