Tsimech2000/virulent_molecular_structures

1"""Input validation, PDB inspection, and sequence-aware alignment."""2 3from __future__ import annotations4 5import io6from dataclasses import dataclass7from pathlib import Path8 9from Bio.PDB import PDBParser, PPBuilder, Superimposer10from Bio.Align import PairwiseAligner11 12 13MAX_UPLOAD_BYTES = 25 * 1024 * 102414 15 16class StructureError(ValueError):17    """Raised for invalid or scientifically unusable structure input."""18 19 20@dataclass(frozen=True)21class PDBReport:22    models: int23    chains: int24    residues: int25    atoms: int26    hetero_residues: int27    waters: int28    alternate_locations: int29    warnings: tuple[str, ...]30 31    def as_dict(self) -> dict[str, int]:32        return {33            "models": self.models,34            "chains": self.chains,35            "residues": self.residues,36            "atoms": self.atoms,37            "hetero residues": self.hetero_residues,38            "waters": self.waters,39            "alternate locations": self.alternate_locations,40        }41 42 43def validate_upload(data: bytes, filename: str, allowed_suffixes: set[str]) -> bytes:44    suffix = Path(filename).suffix.lower()45    if suffix not in allowed_suffixes:46        raise StructureError(f"Unsupported file extension: {suffix or 'none'}")47    if not data:48        raise StructureError("The uploaded file is empty.")49    if len(data) > MAX_UPLOAD_BYTES:50        raise StructureError("The uploaded file exceeds the 25 MB research-app limit.")51    if b"\x00" in data[:4096]:52        raise StructureError("Binary content was detected where a text molecular file was expected.")53    return data54 55 56def _parse_pdb_text(pdb_text: str, structure_id: str = "structure"):57    try:58        structure = PDBParser(QUIET=True).get_structure(structure_id, io.StringIO(pdb_text))59    except Exception as exc:60        raise StructureError(f"Invalid PDB structure: {exc}") from exc61    if not list(structure.get_atoms()):62        raise StructureError("The PDB file contains no readable atoms.")63    return structure64 65 66def inspect_pdb(pdb_text: str) -> PDBReport:67    structure = _parse_pdb_text(pdb_text)68    models = list(structure)69    chains = [chain for model in models for chain in model]70    residues = [residue for chain in chains for residue in chain]71    atoms = [atom for residue in residues for atom in residue]72    hetero = [r for r in residues if r.id[0].strip() and r.id[0] != "W"]73    waters = [r for r in residues if r.id[0] == "W" or r.resname in {"HOH", "WAT"}]74    altlocs = sum(1 for atom in atoms if atom.is_disordered())75    warnings: list[str] = []76    if len(models) > 1:77        warnings.append("Multiple models detected; preparation currently uses the first Open Babel-readable model.")78    if hetero:79        warnings.append("Hetero residues are present. Decide deliberately which cofactors, metals, and additives to retain.")80    if waters:81        warnings.append("Waters are present. Retain only waters supported by the intended docking protocol.")82    if altlocs:83        warnings.append("Alternate atom locations are present and require deliberate conformer selection.")84    return PDBReport(len(models), len(chains), len(residues), len(atoms), len(hetero), len(waters), altlocs, tuple(warnings))85 86 87def sequence_aware_superposition(reference_path: Path, mobile_path: Path):88    reference = _parse_pdb_text(reference_path.read_text(errors="replace"), "reference")89    mobile = _parse_pdb_text(mobile_path.read_text(errors="replace"), "mobile")90    ref_peptides = PPBuilder().build_peptides(reference)91    mob_peptides = PPBuilder().build_peptides(mobile)92    if not ref_peptides or not mob_peptides:93        raise StructureError("Both structures must contain a recognizable polypeptide chain.")94    ref_peptide = max(ref_peptides, key=len)95    mob_peptide = max(mob_peptides, key=len)96    aligner = PairwiseAligner()97    aligner.mode = "global"98    aligner.match_score = 299    aligner.mismatch_score = -1100    aligner.open_gap_score = -10101    aligner.extend_gap_score = -0.5102    alignment = aligner.align(str(ref_peptide.get_sequence()), str(mob_peptide.get_sequence()))[0]103    ref_residues, mob_residues = list(ref_peptide), list(mob_peptide)104    ref_atoms, mob_atoms = [], []105    for (ref_start, ref_end), (mob_start, mob_end) in zip(alignment.aligned[0], alignment.aligned[1]):106        for ref_index, mob_index in zip(range(ref_start, ref_end), range(mob_start, mob_end)):107            if "CA" in ref_residues[ref_index] and "CA" in mob_residues[mob_index]:108                ref_atoms.append(ref_residues[ref_index]["CA"])109                mob_atoms.append(mob_residues[mob_index]["CA"])110    if len(ref_atoms) < 3:111        raise StructureError("Fewer than three aligned Cα pairs were found; the structures are not comparable.")112    superimposer = Superimposer()113    superimposer.set_atoms(ref_atoms, mob_atoms)114    superimposer.apply(mobile.get_atoms())115    return reference, mobile, len(ref_atoms), float(superimposer.rms)116