diff --git a/README.rst b/README.rst index b8e3fef..f25ee47 100644 --- a/README.rst +++ b/README.rst @@ -218,6 +218,18 @@ Name your files the way you want your rows labelled. **One structure, one model, one conformation.** Prodes describes a single molecule, and a PDB file can hold more than one. A file with several ``MODEL`` records, which is what an NMR structure is, is described by its **first** model; the count reaches ``prodes_run.json`` as ``models_in_file``. A residue modelled in several alternate conformations is described by its best occupied one. Both are logged when they happen, and both are decisions the file did not make for you: see `what Prodes calls one residue `_ and `alternate conformations `_. If you want the ensemble averaged rather than one member of it described, run Prodes over the models separately and average the CSVs yourself. +**Only the 20 standard amino acids.** Prodes reads each residue's radii, pKa and hydrophobicity from a table that holds the 20 standard amino acids and nothing else. A residue in the ``ATOM`` records with any other name is **left out, with a warning** that names each one by chain and number, for example ``left out 2 residues that are not one of the 20 standard amino acids: UNK A1, UNK B1``. The commonest case is a structure predicted from a sequence holding an ``X``: Boltz-2 and other predictors build it as ``UNK``. Modified residues (``PCA``, ``MSE``), capping groups (``ACE``, ``NH2``) and force-field residue names (``HIE``, ``HSD``, ``CYX``) are left out the same way when they are written as ``ATOM`` records, so a file renamed by a force-field tool can lose every histidine: read the warning. ``HETATM`` records are not read for the features and are not checked. A structure in which no residue is one of the 20 is refused. + +The run finishes rather than stopping, because one or two such residues move the features of a protein of hundreds of residues little, and a run that stops on them gives nothing. If you would rather it stopped, ``--no-drop-unknown-residues`` (``drop_unknown_residues=False`` in Python) raises a ``ValueError`` naming every such residue as soon as the file is read. What leaving them out means for the numbers: + +* The features describe the structure **with those atoms deleted**, not with the residue replaced by anything. Its surface is gone, the surface of its neighbours that it covered is exposed and counts as theirs, and molecular weight and every area fraction change with it. If the residue would have carried a side-chain charge, that charge is gone too, and with it the net charge and the isoelectric point move. +* A chain end that was removed takes its terminal charge with it: the residue now at that end gets no ``N+`` or ``C-`` group, because in the molecule it is bonded to the removed residue. That is right for an N-terminal pyroglutamate, which has no free amine and is a common reason an antibody heavy chain starts with ``X``. It is wrong for an ordinary residue that merely has an unknown name, whose amine or carboxylate is still charged in the real protein. +* A pKa from your ``--pka`` file for a removed residue, or a terminal pKa for the residue next to a removed end, is not applied. +* A residue set naming a removed residue loses it, with a warning; a set naming nothing else is refused. +* The removed residues are logged, and listed in ``prodes_run.json`` as ``removed_residues``, for example ``["UNK A1", "UNK B1"]``. Use that field to flag or exclude the rows of a dataset whose structures lost residues, since they are not described the same way as the rest. + +If you know what the residue is, put the standard residue into the structure or into the sequence it was predicted from, rather than dropping it: a germline residue in an antibody constant domain, for instance. For a pyroglutamate, ``GLN`` restores the side chain's surface but also gives the chain a free N-terminal amine and its charge, which the pyroglutamate does not have. + The output bundle ----------------- @@ -854,6 +866,7 @@ Improvements in this fork * **A blank element column is inferred rather than fatal** — from version 8.1 an atom whose element column is blank has its element guessed from its name, the same convention ``Bio.PDB.Atom`` uses. A blank element used to reach the van der Waals radius lookup and raise ``KeyError``, so a PDB file with no element column produced no features at all; no shipped structure has one, so no feature value moves. The count of atoms inferred is logged once per file and recorded in ``prodes_run.json`` as ``inferred_elements``. See `how Prodes reads a PDB file `_. * **pKa predictions are keyed by chain** — from version 9.0 a predicted pKa is looked up by chain as well as residue number, rather than by number alone. A structure with more than one chain, an antibody's heavy and light chain included, could have a value predicted for one chain applied to another chain's same-numbered residue instead, silently whenever the two residues shared a type. ``convert_propka`` and ``convert_pypka`` now read the chain out of their source files; ``convert_hpp`` cannot, because H++'s own output carries no chain, and marks its predictions accordingly. A ``--pka`` file written before this version has no chain in it at all and is still read, applied to every chain exactly as before, with a warning recommending it be reconverted. See `pKa values and protonation states`_. * **Surface patches** — from version 10.0, every run groups neighbourhood-averaged hydrophobicity and potential into contiguous patches and reports their areas, as 11 more default features, over the whole surface and over any named residue set, such as antibody CDRs. Every surface point in the bundle also carries its area, the residue owning it and its patch. See `surface patches `_. +* **Unknown residues left out with a warning** — from version 11.0 a residue that is not one of the 20 standard amino acids, such as the ``UNK`` a structure predictor builds for an ``X``, is left out of the calculation, named in a warning by chain and number, and listed in ``prodes_run.json``. Such a structure used to fail with a ``KeyError`` naming only the residue type, once the surface calculation had started. ``--no-drop-unknown-residues`` raises a ``ValueError`` naming them instead. See `Input files`_. * **Bug fixes** — trimean edge case for small arrays, ``surface_exit`` ``None`` handling in shell potential mapping, and a ``read_propka`` argument bug in the PDB parser. * **Test suite** — unit and regression tests added, including a committed reference output file (``tests/data/ARH96693_prodes_orig_output.csv``) generated by the original unrefactored code. The regression test verifies that every feature column and every feature value produced by this fork matches the original output within tolerance. @@ -864,6 +877,8 @@ With ``--full-features``, the output has the same 105 columns in the same order The **values** are identical only with ``--ionic-strength 0``, and then only for a structure with no disulfide bonds and no alternate conformations. From version 9.0 a predicted pKa is applied by chain as well as residue number; this only changes a structure's features if its ``--pka`` file is reconverted with the current ``convert_propka`` or ``convert_pypka``, and then only if the structure has more than one chain and a chain-specific value lands on a residue that previously took another chain's value or none at all: on the 4NZU Fab, reconverting the same PROPKA run moves 13 of the 105 columns, the isoelectric point among them, from 6.796 to 6.808. A ``--pka`` file that predates this version, or one produced by the current ``convert_hpp``, is still applied to every chain exactly as before and changes nothing. See `pKa values and protonation states`_. From version 8.0 a residue that carries an insertion code is its own residue rather than part of the one before it, which moves the features of any structure that has one, antibodies numbered by the Kabat or Chothia schemes included: on the 4NZU Fab 69 of the 105 columns move and 36 do not. Nothing geometric moves, because no atom moves and no atom changes its own residue name: ``Area``, ``NSurfPoints``, both ``Shape`` columns and the whole hydrophobicity block are unchanged, while everything that reads a residue's type, mass or charge does move. A file holding more than one model is described by its first, which changes every feature of an NMR ensemble, all of them from nonsense to a number. Structures with neither an insertion code nor a second model are unaffected. See `what Prodes calls one residue `_. From version 7.0 only one conformation of a disordered residue is described, which changes the features of any structure that has one; the size of the change tracks how many *titratable* residues are disordered rather than how many residues are, so a structure whose alternates sit on uncharged residues barely moves. AlphaFold models contain no alternate conformations and are unaffected. See `alternate conformations `_. From version 6.0 a cysteine bonded into a disulfide is not titrated, which changes the charge-derived features of any structure that has one. See `cysteines and disulfide bonds `_. From version 5.0 the electrostatic potential is screened by default, which changes the 38 ``SurfEp`` columns; the other 67, including every ``ShellEp`` column, are unaffected either way. The regression test in ``tests/test_sasa.py`` therefore runs at zero ionic strength when comparing against the committed reference file ``tests/data/ARH96693_prodes_orig_output.csv``. See `Ionic strength and screening`_. +From version 11.0 ``prodes.read`` and ``PDBparser.parse`` leave out, with a warning, every residue in the ``ATOM`` records outside the 20 standard amino acids, such as ``UNK``, which they used to return as part of the structure. ``calculate`` used to fail on such a structure and now describes it without those residues. No feature value of any other structure changes. See `Input files`_. + From version 10.0 every run appends the 11 surface patch columns to the feature table, after all the others, and the surface point table in the bundle has five more columns: ``area_A2``, ``residue``, ``local_mhp``, ``local_ep_volts`` and ``mhp_patch_id``. A model trained on the columns of an earlier version has to select those columns by name. No earlier column changes value. The default reduced output is a strict subset of those columns, in the same order and with identical values. No column is renamed and no column is recomputed, so a reduced run and a full run of the same structure agree exactly on the 65 columns they share. diff --git a/docs/pdb_reading.md b/docs/pdb_reading.md index 2884491..9332272 100644 --- a/docs/pdb_reading.md +++ b/docs/pdb_reading.md @@ -18,6 +18,8 @@ file text │ ├─ prodes.io.conformers.elect_conformers one conformation per residue, the 7.0 rules │ + ├─ prodes.io.parser.remove_unknown_residues residues outside the 20 standard amino acids: left out with a warning, from 11.0 + │ └─ prodes.io.parser.build_structure Structure, Chain, Residue, Atom ``` diff --git a/pyproject.toml b/pyproject.toml index b31946f..f6b1c5f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "prodes" -version = "10.0.0" +version = "11.0.0" description = "A package to calculate protein surface descriptors" readme = "README.rst" license = {text = "MIT"} diff --git a/src/prodes/__init__.py b/src/prodes/__init__.py index 9a8b5f5..1048aab 100644 --- a/src/prodes/__init__.py +++ b/src/prodes/__init__.py @@ -14,6 +14,7 @@ def run_prodes( mem_limit_mb=None, ionic_strength_molar=run.DEFAULT_IONIC_STRENGTH_MOLAR, patches=None, + drop_unknown_residues=True, ): """Runs Prodes similar as the commandline tool. @@ -40,6 +41,9 @@ def run_prodes( patches: a prodes.calculations.patches.PatchSettings for the surface patch features, such as residue sets to restrict them to. None means the default settings and no residue sets. + drop_unknown_residues: leave out the residues that are not one of the + 20 standard amino acids, such as UNK, with a warning (the + default), or raise when False. See prodes.run.calculate. """ calculate( @@ -53,19 +57,26 @@ def run_prodes( mem_limit_mb=mem_limit_mb, ionic_strength_molar=ionic_strength_molar, patches=patches, + drop_unknown_residues=drop_unknown_residues, ) -def read(file): +def read(file, drop_unknown_residues=True): """Loads files List of supported extentions: .pdb .pdb.zip .pka: a pKa file in the json layout written by prodes.io.pka_converter, - rather than the raw output of PROPKA or pypka""" + rather than the raw output of PROPKA or pypka + + A residue that is not one of the 20 standard amino acids, such as UNK, is + left out with a warning and listed in structure.removed_residues. With + drop_unknown_residues=False it raises a ValueError naming it instead. See + prodes.io.parser.remove_unknown_residues. + """ if file.endswith((".pdb", ".pdb.zip")): - structure = parser.PDBparser().parse(file) + structure = parser.PDBparser().parse(file, drop_unknown_residues=drop_unknown_residues) return structure elif file[-4:] == ".pka": diff --git a/src/prodes/core/chain.py b/src/prodes/core/chain.py index 2070c73..59b659b 100644 --- a/src/prodes/core/chain.py +++ b/src/prodes/core/chain.py @@ -8,6 +8,12 @@ def __init__(self, name, structure): self.atoms = np.empty([0]) self.name = name self.structure = structure + # Set by prodes.io.parser.clear_removed_termini when the residue at that + # end was left out as unknown. The residue now at that end is bonded to + # the removed one in the molecule, so Residue.pkas gives it no terminal + # group. + self.n_terminus_removed = False + self.c_terminus_removed = False self._x = None self._y = None self._z = None diff --git a/src/prodes/core/residue.py b/src/prodes/core/residue.py index fd2db2b..973e96d 100644 --- a/src/prodes/core/residue.py +++ b/src/prodes/core/residue.py @@ -78,9 +78,12 @@ def pkas(self): # every consumer of pkas agrees about it. if pka is not None and self.disulfide_partner is None: pkas.append({self.name: pka}) - if self == self.chain.residues[0]: + # A chain whose end residue was left out as unknown has no free + # terminus at that end: the residue there now is bonded to the + # removed one in the molecule. + if self == self.chain.residues[0] and not self.chain.n_terminus_removed: pkas.append({"N+": data.n_term_pka()}) - elif self == self.chain.residues[-1]: + elif self == self.chain.residues[-1] and not self.chain.c_terminus_removed: pkas.append({"C-": data.c_term_pka()}) if len(pkas) > 0: self._pka = pkas diff --git a/src/prodes/core/structure.py b/src/prodes/core/structure.py index 511352d..0c3ffbc 100644 --- a/src/prodes/core/structure.py +++ b/src/prodes/core/structure.py @@ -41,6 +41,10 @@ def __init__(self, name=None): # for a single model file read for one record type with no alternate # conformations. self.inferred_elements = 0 + # The residues prodes.io.parser.remove_unknown_residues left out because + # they are not one of the 20 standard amino acids, as "UNK A1" labels. + # Empty unless the structure was read with drop_unknown_residues. + self.removed_residues = [] def _compute_centroid(self): coords = np.array([[a.x, a.y, a.z] for a in self.heavy_atoms]) diff --git a/src/prodes/io/parser.py b/src/prodes/io/parser.py index dcac98f..0ce7085 100644 --- a/src/prodes/io/parser.py +++ b/src/prodes/io/parser.py @@ -10,7 +10,7 @@ from prodes.core.chain import Chain from prodes.core.residue import Residue from prodes.core.structure import Structure -from prodes.io.conformers import elect_conformers, report_alternate_conformers +from prodes.io.conformers import elect_conformers, report_alternate_conformers, residue_key, residue_label from prodes.io.pdb_reader import COORDINATE_DECIMALS, pdb_lines, read_atom_records logger = logging.getLogger(__name__) @@ -377,7 +377,7 @@ def build_structure(name, records): return structure -def parse_pdb_text(text, name, identifier="ATOM"): +def parse_pdb_text(text, name, identifier="ATOM", drop_unknown_residues=True): """Reads the text of a PDB file into a Structure. The records are read in full before anything is built, because which @@ -401,10 +401,15 @@ def parse_pdb_text(text, name, identifier="ATOM"): text: the text of a PDB file. name: the name to give the structure, which becomes the output ID. identifier: the record type to build atoms from, ATOM or HETATM. + drop_unknown_residues: leave out the protein residues that are not + one of the 20 standard amino acids with a warning, the default, or + raise when False. See remove_unknown_residues. Raises: - ValueError: if the file holds no records of the requested type, or if - any coordinate record could not be read. + ValueError: if the file holds no records of the requested type, if + any coordinate record could not be read, or if a protein residue + is not one of the 20 standard amino acids and + drop_unknown_residues is False. """ lines = pdb_lines(text) @@ -428,21 +433,152 @@ def parse_pdb_text(text, name, identifier="ATOM"): kept, report = elect_conformers(records) + # After the election, because a file can model one position as two + # different amino acids and only the elected one is the residue. Only for + # protein atoms: a HETATM parse reads ligands and waters, none of which is + # in the residue table. + removed = [] + ends_removed = {} + if identifier == "ATOM": + kept, removed, ends_removed = remove_unknown_residues(kept, name, drop_unknown_residues) + # A residue the election renamed to an unknown name and that was then + # removed is described by the removal, not by a rename the features + # never saw, so the two lists in the run record do not contradict. + removed_labels = {label.split(" ", 1)[1] for label in removed} + report.renamed = [entry for entry in report.renamed if residue_label(entry[0]) not in removed_labels] + structure = build_structure(name, kept) structure.models = models structure.alternate_conformers = report structure.inferred_elements = sum(1 for record in kept if record.element_inferred) + structure.removed_residues = removed report_alternate_conformers(report, name) structure.residues[-1].terminus = "C" + clear_removed_termini(structure, ends_removed) assign_disulfides(structure, ssbond_records(lines) if identifier == "ATOM" else []) return structure +def remove_unknown_residues(records, name, drop): + """Removes the records of residues that are not in the residue table, or raises. + + The residue table holds the 20 standard amino acids. Anything else, such + as the UNK a structure predictor builds for an X in its input sequence, + has no radius, pKa or hydrophobicity to read, and used to stop the run + with a KeyError deep in the surface calculation. + + With drop set, the default of every caller, the residues are removed with + a warning naming each one, and the features describe the structure + without their atoms. One or two such residues in an antibody of 1300 move + the features little, and a run that stops on them gives nothing at all, + which for an unattended caller is the worse outcome. Without drop, a + ValueError names every such residue, by chain and number, before any slow + phase has started. Which chain ends were among them is + returned as well: a residue removed from the end of a chain takes its free + terminal charge with it rather than handing it to its neighbour, which is + bonded to something in the molecule and has no free amine or carboxylate. + + Returns: + (kept records, removed residues as "UNK A1" labels in file order, + {chain name: set of "N" and/or "C" for each end that was removed}). + + Raises: + ValueError: if a residue is not in the residue table and drop is not + set, or if every residue is. + """ + + from prodes import data + + known = data.all_residues() + # A residue is unknown if any of its records names an unknown residue, + # since every atom reads its own residue name for its radius. The ends of + # a chain are its first and last residue in order of first appearance, + # which is the order build_structure gives chain.residues. + unknown = {} + chain_keys = {} + for record in records: + key = residue_key(record) + chain_keys.setdefault(record.chain_name, {}).setdefault(key, None) + if record.residue_name not in known: + unknown.setdefault(key, record.residue_name) + + if not unknown: + return records, [], {} + + labels = [f"{residue_name} {residue_label(key)}" for key, residue_name in unknown.items()] + if not drop: + raise ValueError( + f"{name} holds {count_residues(labels)} that {'is' if len(labels) == 1 else 'are'} not one of the 20 standard amino acids, so Prodes has no data for {'it' if len(labels) == 1 else 'them'}: " + f"{', '.join(labels)}. Replace them in the structure, or leave drop_unknown_residues at its default of True (do not pass --no-drop-unknown-residues) to leave them out with a warning." + ) + + kept = [record for record in records if residue_key(record) not in unknown] + if not kept: + raise ValueError(f"{name} holds no residue that is one of the 20 standard amino acids") + + ends_removed = {} + for chain_name, keys in chain_keys.items(): + keys = list(keys) + ends = set() + if keys[0] in unknown: + ends.add("N") + if keys[-1] in unknown: + ends.add("C") + if ends: + ends_removed[chain_name] = ends + + logger.warning( + "%s: left out %s that %s not one of the 20 standard amino acids: %s. The features describe the structure without their atoms.", + name, + count_residues(labels), + "is" if len(labels) == 1 else "are", + ", ".join(labels), + ) + return kept, labels, ends_removed + + +def count_residues(labels): + """Returns "1 residue" or "3 residues", for a message.""" + + return f"{len(labels)} residue{'' if len(labels) == 1 else 's'}" + + +def clear_removed_termini(structure, ends_removed): + """Takes the terminus away from a chain whose end residue was removed. + + The residue now at that end is bonded to the removed one in the molecule, + so it carries no terminal charge. Two things decide a terminal charge and + both are set here, by the same rule: the flags on the chain, which + Residue.pkas reads to decide which terminal group a residue has, and + residue.terminus, which Atom.charge reads to decide which atom carries it. + Runs after parse_pdb_text has written the structure's last C terminus, so + that nothing writes one back afterwards. + + Only the chains that lost an end are touched, so every other chain keeps + the termini build_structure gave it. + """ + + for chain in structure.chains: + ends = ends_removed.get(chain.name) + if not ends: + continue + chain.n_terminus_removed = "N" in ends + chain.c_terminus_removed = "C" in ends + for residue in chain.residues: + residue.terminus = None + # The same if/elif as Residue.pkas, so a chain left with one residue + # gets the terminus pkas gives it a group for. + if not chain.n_terminus_removed: + chain.residues[0].terminus = "N" + elif not chain.c_terminus_removed: + chain.residues[-1].terminus = "C" + + class PDBparser: - def parse(self, file, identifier="ATOM"): + def parse(self, file, identifier="ATOM", drop_unknown_residues=True): """Parses pdb files, either plain or held in a zip archive A zipped structure is read straight out of the archive, so that callers @@ -469,7 +605,7 @@ def parse(self, file, identifier="ATOM"): else: raise ValueError("File extension not recognized by parser") - return parse_pdb_text(read_pdb_text(path), name, identifier) + return parse_pdb_text(read_pdb_text(path), name, identifier, drop_unknown_residues) def read_pka(file): diff --git a/src/prodes/output.py b/src/prodes/output.py index db45d50..765d88d 100644 --- a/src/prodes/output.py +++ b/src/prodes/output.py @@ -129,7 +129,7 @@ def prodes_version(): return "unknown" -def run_metadata(pdb_file, settings, n_points, ep, disulfides=0, alternate_conformers=None, models=1, inferred_elements=0): +def run_metadata(pdb_file, settings, n_points, ep, disulfides=0, alternate_conformers=None, models=1, inferred_elements=0, removed_residues=None): """Returns the record of what was run, written into every bundle. The disulfide count is here rather than in the feature table because it @@ -149,6 +149,11 @@ def run_metadata(pdb_file, settings, n_points, ep, disulfides=0, alternate_confo The inferred element count says how many atoms of the structure had their radius, and so their contribution to the surface, decided by a guess at their element rather than by the file's own columns. See docs/pdb_reading.md. + + The removed residues are the ones left out with drop_unknown_residues + because they are not one of the 20 standard amino acids, as "UNK A1" + labels. The shipped file still holds them, so this is the only place that + says the features describe the structure without their atoms. """ return { @@ -160,6 +165,7 @@ def run_metadata(pdb_file, settings, n_points, ep, disulfides=0, alternate_confo "alternate_conformers": alternate_conformers or {}, "models_in_file": int(models), "inferred_elements": int(inferred_elements), + "removed_residues": list(removed_residues or []), "surface_points": int(n_points), "ep_min_volts": round(float(ep.min()), 3) if n_points else None, "ep_max_volts": round(float(ep.max()), 3) if n_points else None, diff --git a/src/prodes/run.py b/src/prodes/run.py index 3015e2f..a6c3e59 100644 --- a/src/prodes/run.py +++ b/src/prodes/run.py @@ -1,5 +1,7 @@ import argparse +import dataclasses import faulthandler +import logging import math import os @@ -24,6 +26,8 @@ from prodes.parallel import run_tasks, validate_settings, worker_mem_limit_mb from prodes.viz import surface_point_arrays +logger = logging.getLogger(__name__) + load_dotenv() # Shared state for the worker functions below. The parent fills these in before @@ -111,6 +115,12 @@ def parse_arguments(): default=None, help="Number of worker processes, Linux only (default: from PRODES_N_WORKERS env var, or half the logical CPUs)", ) + parser.add_argument( + "--drop-unknown-residues", + action=argparse.BooleanOptionalAction, + default=True, + help="Leave out residues that are not one of the 20 standard amino acids, such as UNK, with a warning, and list them in prodes_run.json (default). --no-drop-unknown-residues stops with an error naming them instead", + ) parser.add_argument( "--chunksize", type=int, @@ -154,7 +164,36 @@ def parse_arguments(): residue_sets=read_residue_sets(arg.residue_set), ) - return pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar, patches + return pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar, patches, arg.drop_unknown_residues + + +def without_removed_residues(patches, removed_residues): + """Returns the patch settings with the removed residues taken out of every residue set. + + A residue left out with drop_unknown_residues has no surface to restrict a + patch to, and a set naming it would otherwise fail the check that every + label is in the structure, which is the check for a typo. Taking it out is + said in the log, and the run record lists the removed residues already. + + Raises: + ValueError: if a residue set names nothing but removed residues. + """ + + if not removed_residues or not patches.residue_sets: + return patches + + # Removed residues are written "UNK A1", and residue sets name residues by + # the label alone. + removed = {label.split(" ", 1)[1] for label in removed_residues} + residue_sets = {} + for name, labels in patches.residue_sets.items(): + kept = [label for label in labels if label not in removed] + if not kept: + raise ValueError(f"residue set {name} names only residues that were left out as unknown: {', '.join(labels)}") + if len(kept) < len(labels): + logger.warning("residue set %s loses the residues that were left out as unknown: %s", name, ", ".join(sorted(set(labels) - set(kept)))) + residue_sets[name] = kept + return dataclasses.replace(patches, residue_sets=residue_sets) def read_residue_sets(specifications): @@ -513,10 +552,10 @@ def calculate_structure_features(structure, ph, r_probe, features: dict, full_fe return features -def prepare_structure(pdb_file, pkas_file): +def prepare_structure(pdb_file, pkas_file, drop_unknown_residues=True): """loads and prepares the structure for the calculations""" - structure = ps.PDBparser().parse(pdb_file) + structure = ps.PDBparser().parse(pdb_file, drop_unknown_residues=drop_unknown_residues) if pkas_file: pkas = ps.read_pka(pkas_file) structure.redo_pkas(pkas) @@ -554,6 +593,7 @@ def calculate( mem_limit_mb=None, ionic_strength_molar=DEFAULT_IONIC_STRENGTH_MOLAR, patches=None, + drop_unknown_residues=True, ): """Calculates a list of supported features and returns a csv file @@ -580,6 +620,13 @@ def calculate( patch features, such as residue sets to restrict them to. None, the default, means PatchSettings(): the default settings and no residue sets. The patch features are part of every run. + drop_unknown_residues: True, the default, leaves out the residues + that are not one of the 20 standard amino acids, such as the UNK a + structure predictor builds for an X in its sequence, with a + warning. They are listed in prodes_run.json as removed_residues, a + residue set naming one of them loses it, and a chain end that was + removed carries no terminal charge. False raises a ValueError + naming them instead. See README, "Input files". Thread safety Safe to run in several processes, one structure each, which is how the @@ -598,7 +645,8 @@ def calculate( # Both checks are cheap, and the surface phase that follows is not. patches.resolved_hydrophobic_threshold(hydro_scale) - structure = prepare_structure(pdb_file, pkas_file) + structure = prepare_structure(pdb_file, pkas_file, drop_unknown_residues) + patches = without_removed_residues(patches, structure.removed_residues) check_residue_sets(patches.residue_sets, [residue.label for residue in structure.residues]) # Reported before the slow phases, so that a structure whose disulfides were # not found can be spotted without waiting for the run to finish. @@ -648,6 +696,7 @@ def calculate( "full_features": full_features, "ionic_strength_molar": ionic_strength_molar, "patches": patches.record(hydro_scale), + "drop_unknown_residues": drop_unknown_residues, } metadata = run_metadata( pdb_file, @@ -658,6 +707,7 @@ def calculate( alternate_conformers=structure.alternate_conformers.summary() if structure.alternate_conformers else None, models=structure.models, inferred_elements=structure.inferred_elements, + removed_residues=structure.removed_residues, ) return write_bundle(out_file, structure.name, calculated_features, coords, ep, lipo, pdb_file, metadata, hydro_scale, point_columns=point_columns) diff --git a/tests/test_cli.py b/tests/test_cli.py index 72a3b8b..32da5d5 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -27,7 +27,9 @@ def test_defaults_are_the_documented_ones(monkeypatch): """A bare invocation gives the reduced feature set and no explicit memory budget.""" monkeypatch.delenv("PRODES_FULL_FEATURES", raising=False) - pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar, patches = parse(monkeypatch) + pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar, patches, drop_unknown_residues = parse( + monkeypatch + ) assert (pdb_file, out_file, pkas_file) == ("in.pdb", "out.zip", None) assert (ph, r_probe, hydro_scale) == (7, 1.4, "mj_scaled") @@ -35,6 +37,7 @@ def test_defaults_are_the_documented_ones(monkeypatch): assert mem_limit_mb is None assert ionic_strength_molar == DEFAULT_IONIC_STRENGTH_MOLAR assert patches == PatchSettings() + assert drop_unknown_residues is True def test_n_workers_flag_sets_the_environment_variable(monkeypatch): @@ -58,7 +61,7 @@ def test_chunksize_flag_sets_the_environment_variable(monkeypatch): def test_mem_limit_flag_is_returned_rather_than_exported(monkeypatch): """--mem-limit is passed to calculate directly, so it does not need the environment.""" - *_, mem_limit_mb, _ionic_strength, _patches = parse(monkeypatch, "--mem-limit", "256") + *_, mem_limit_mb, _ionic_strength, _patches, _drop = parse(monkeypatch, "--mem-limit", "256") assert mem_limit_mb == 256 @@ -93,7 +96,9 @@ def fake_calculate(*args, **kwargs): captured["kwargs"] = kwargs monkeypatch.setattr("prodes.run.calculate", fake_calculate) - monkeypatch.setattr("sys.argv", ["prodes", str(tmp_path / "in.pdb"), str(tmp_path / "out.zip"), "--ph", "5.5", "--full-features"]) + monkeypatch.setattr( + "sys.argv", ["prodes", str(tmp_path / "in.pdb"), str(tmp_path / "out.zip"), "--ph", "5.5", "--full-features", "--no-drop-unknown-residues"] + ) main() @@ -111,6 +116,7 @@ def fake_calculate(*args, **kwargs): None, DEFAULT_IONIC_STRENGTH_MOLAR, PatchSettings(), + False, ) assert captured["kwargs"] == {} @@ -129,6 +135,7 @@ def test_calculate_signature_matches_the_parsed_order(): "mem_limit_mb", "ionic_strength_molar", "patches", + "drop_unknown_residues", ] @@ -137,7 +144,7 @@ def test_patch_options_build_the_settings(monkeypatch, tmp_path): labels = tmp_path / "cdr.txt" labels.write_text("# CDR H3\nH95\nH100A # an insertion\n\nL91\n") - *_, patches = parse( + *_, patches, _drop = parse( monkeypatch, "--patch-radius", "4", diff --git a/tests/test_output_bundle.py b/tests/test_output_bundle.py index abb4c7f..9e684c6 100644 --- a/tests/test_output_bundle.py +++ b/tests/test_output_bundle.py @@ -107,10 +107,12 @@ def test_the_run_record_says_what_was_run(bundle_dir): "full_features": False, "ionic_strength_molar": 0.15, "patches": {"radius_A": 5.0, "min_area_A2": 100.0, "hydrophobic_threshold": 0.0, "ep_threshold_volts": 0.5, "residue_sets": {}}, + "drop_unknown_residues": True, } # ARH96693 is an AlphaFold model with three cysteines and no disulfide, so # the key has to be present and zero rather than absent. assert record["disulfides"] == 0 + assert record["removed_residues"] == [] assert record["surface_points"] > 0 assert record["ep_min_volts"] <= record["ep_max_volts"] assert record["input_file"].endswith("ARH96693.pdb.zip") diff --git a/tests/test_run_prodes.py b/tests/test_run_prodes.py index 8f5bf22..06d0bcf 100644 --- a/tests/test_run_prodes.py +++ b/tests/test_run_prodes.py @@ -61,6 +61,7 @@ def fake_calculate(pdb_file, out_file, **kwargs): full_features=True, mem_limit_mb=512, patches=PatchSettings(radius=4.0), + drop_unknown_residues=False, ) assert captured == { @@ -74,6 +75,7 @@ def fake_calculate(pdb_file, out_file, **kwargs): "mem_limit_mb": 512, "ionic_strength_molar": 0.15, "patches": PatchSettings(radius=4.0), + "drop_unknown_residues": False, } diff --git a/tests/test_unknown_residues.py b/tests/test_unknown_residues.py new file mode 100644 index 0000000..699a283 --- /dev/null +++ b/tests/test_unknown_residues.py @@ -0,0 +1,307 @@ +"""Residues that are not one of the 20 standard amino acids, issue #37. + +A structure predicted from a sequence holding an X has a residue named UNK, +which the residue table has no radius, pKa or hydrophobicity for. By default +it is left out with a warning and recorded, and a chain end it occupied +carries no terminal charge; with drop_unknown_residues=False the structure is +refused when it is read, naming every such residue. +""" + +import json +import logging +import zipfile + +import pytest + +import prodes +from prodes.calculations.patches import PatchSettings +from prodes.io import parser +from prodes.io.pka_converter import ANY_CHAIN +from prodes.output import run_metadata +from prodes.run import calculate, without_removed_residues +from tests.pdb_records import atom_line, write_structure + +# Backbone and CB, the atoms Boltz-2 builds for an UNK. GLY has no CB, but an +# extra atom changes nothing the parser decides, so every residue gets one. +BACKBONE = [("N", 0.0, 0.0, 0.0, "N"), ("CA", 1.2, 0.5, 0.6, "C"), ("C", 2.4, 0.0, 0.0, "C"), ("O", 2.4, -1.2, -0.6, "O"), ("CB", 1.2, 2.0, 1.4, "C")] + + +def residue_lines(serial, residue_name, chain, number, x, y, altloc="", occupancy=1.00): + """Returns the five heavy atoms of one residue, placed at x along its chain.""" + + return [ + atom_line(serial + offset, name, residue_name, chain, number, x + dx, y + dy, dz, element, altloc=altloc, occupancy=occupancy) + for offset, (name, dx, dy, dz, element) in enumerate(BACKBONE) + ] + + +def chain_lines(chain, names, y=0.0, serial=1): + """Returns the records of one chain, residue n named names[n - 1], 3.8 A apart.""" + + lines = [] + for number, residue_name in enumerate(names, start=1): + lines += residue_lines(serial + len(lines), residue_name, chain, number, 3.8 * number, y) + return lines + + +def pdb_text(*chains): + """Returns PDB text for the given chain records, chains 20 A apart.""" + + lines = [] + for offset, (chain, names) in enumerate(chains): + lines += chain_lines(chain, names, y=20.0 * offset, serial=len(lines) + 1) + return "\n".join(lines) + "\nEND\n" + + +# Unknown at the start, in the middle and at the end of chain A; chain B has none. +MIXED = pdb_text(("A", ["UNK", "VAL", "ALA", "UNK", "GLY", "UNK"]), ("B", ["SER", "LYS", "ASP"])) + + +def parse(text, drop): + return parser.parse_pdb_text(text, "test", drop_unknown_residues=drop) + + +def residue(structure, chain, number): + return next(r for r in structure.residues if r.chain.name == chain and r.number == number) + + +def terminal_groups(residue): + return [key for group in (residue.pkas or []) for key in group if key in ("N+", "C-")] + + +def test_by_default_an_unknown_residue_is_left_out_with_a_warning(caplog): + """An unattended run finishes: the structure is read, and the log names what was left out.""" + + with caplog.at_level(logging.WARNING, logger="prodes.io.parser"): + structure = parser.parse_pdb_text(MIXED, "test") + + assert structure.removed_residues == ["UNK A1", "UNK A4", "UNK A6"] + assert "left out 3 residues that are not one of the 20 standard amino acids: UNK A1, UNK A4, UNK A6" in caplog.text + + +def test_without_the_option_an_unknown_residue_is_refused_when_the_structure_is_read(): + """Every unknown residue is named by chain and number, and the setting that removes them is named.""" + + with pytest.raises(ValueError) as raised: + parse(MIXED, drop=False) + + message = str(raised.value) + assert "UNK A1, UNK A4, UNK A6" in message + assert "drop_unknown_residues" in message + assert "--no-drop-unknown-residues" in message + + +def test_one_unknown_residue_is_counted_in_the_singular(): + """One X in a sequence is the common case, so its message has to read right.""" + + with pytest.raises(ValueError, match="holds 1 residue that is not one of the 20 standard amino acids, so Prodes has no data for it: UNK A1"): + parse(pdb_text(("A", ["UNK", "VAL", "SER"])), drop=False) + + +def test_prodes_read_has_the_same_default_and_the_same_option(tmp_path): + """The library entry point reads through the same parser, so it needs the same option.""" + + path = tmp_path / "mixed.pdb" + path.write_text(MIXED) + + assert [r.label for r in prodes.read(str(path)).residues] == ["A2", "A3", "A5", "B1", "B2", "B3"] + with pytest.raises(ValueError, match="UNK A1"): + prodes.read(str(path), drop_unknown_residues=False) + + +def test_the_option_removes_them_and_lists_them_in_file_order(caplog): + """The removed residues are kept on the structure and named in the log, in the order the file holds them.""" + + structure = parse(pdb_text(("B", ["UNK", "VAL"]), ("A", ["SER", "UNK", "LYS"])), drop=True) + assert structure.removed_residues == ["UNK B1", "UNK A2"] + + with caplog.at_level(logging.WARNING, logger="prodes.io.parser"): + structure = parse(MIXED, drop=True) + + assert [r.label for r in structure.residues] == ["A2", "A3", "A5", "B1", "B2", "B3"] + assert all(atom.residue_name != "UNK" for atom in structure.atoms) + assert structure.removed_residues == ["UNK A1", "UNK A4", "UNK A6"] + assert "UNK A1, UNK A4, UNK A6" in caplog.text + + +def test_a_removed_first_residue_takes_the_n_terminus_with_it(): + """The new first residue is bonded to the removed one, so it has no free amine.""" + + structure = parse(MIXED, drop=True) + first = residue(structure, "A", 2) + + assert first.chain.n_terminus_removed is True + assert first.terminus is None + assert "N+" not in terminal_groups(first) + assert next(atom for atom in first.atoms if atom.name == "N").charge(ph=7) == 0 + + +def test_a_removed_last_residue_takes_the_c_terminus_with_it(): + """The same for the carboxylate at the other end.""" + + structure = parse(MIXED, drop=True) + last = residue(structure, "A", 5) + + assert last.chain.c_terminus_removed is True + assert last.terminus is None + assert "C-" not in terminal_groups(last) + assert next(atom for atom in last.atoms if atom.name == "C").charge(ph=7) == 0 + + +def test_a_chain_with_nothing_removed_keeps_both_termini(): + """Chain B is untouched, termini, groups and charges alike.""" + + structure = parse(MIXED, drop=True) + first, last = residue(structure, "B", 1), residue(structure, "B", 3) + + assert (first.terminus, last.terminus) == ("N", "C") + assert terminal_groups(first) == ["N+"] + assert terminal_groups(last) == ["C-"] + assert next(atom for atom in first.atoms if atom.name == "N").charge(ph=7) == 1 + assert next(atom for atom in last.atoms if atom.name == "C").charge(ph=7) == -1 + + +def test_a_removed_middle_residue_changes_no_terminus(): + """Removing A4 leaves A3 and A5 as ordinary residues.""" + + structure = parse(pdb_text(("A", ["SER", "VAL", "ALA", "UNK", "GLY", "LYS"])), drop=True) + + assert structure.removed_residues == ["UNK A4"] + assert [r.terminus for r in structure.residues] == ["N", None, None, None, "C"] + assert terminal_groups(residue(structure, "A", 1)) == ["N+"] + assert terminal_groups(residue(structure, "A", 6)) == ["C-"] + + +@pytest.mark.parametrize( + ("names", "kept", "terminus", "group", "atom_name", "charge"), + [ + (["UNK", "VAL"], 2, "C", "C-", "C", -1), + (["VAL", "UNK"], 1, "N", "N+", "N", 1), + ], +) +def test_a_chain_left_with_one_residue_keeps_the_end_that_was_not_removed(names, kept, terminus, group, atom_name, charge): + """The terminus attribute, the pKa group and the atom charge agree on which end survived.""" + + structure = parse(pdb_text(("A", names)), drop=True) + remaining = residue(structure, "A", kept) + + assert remaining.terminus == terminus + assert terminal_groups(remaining) == [group] + assert next(atom for atom in remaining.atoms if atom.name == atom_name).charge(ph=7) == charge + + +def test_a_structure_of_nothing_but_unknown_residues_is_refused_either_way(): + """There is nothing left to describe, with or without the option.""" + + text = pdb_text(("A", ["UNK", "UNK"])) + + with pytest.raises(ValueError, match="UNK A1, UNK A2"): + parse(text, drop=False) + with pytest.raises(ValueError, match="no residue that is one of the 20 standard amino acids"): + parse(text, drop=True) + + +def test_a_structure_without_unknown_residues_is_the_same_with_and_without_the_option(): + """The option changes nothing for a structure it has nothing to do for.""" + + text = pdb_text(("A", ["SER", "VAL", "LYS"]), ("B", ["ASP", "GLY"])) + plain, dropped = parse(text, drop=False), parse(text, drop=True) + + def summary(structure): + return [(r.label, r.name, r.terminus, r.pkas, len(r.atoms)) for r in structure.residues] + + assert summary(plain) == summary(dropped) + for ph in (4, 7, 10): + assert [atom.charge(ph=ph, formal=False) for atom in plain.atoms] == [atom.charge(ph=ph, formal=False) for atom in dropped.atoms] + assert dropped.removed_residues == [] + assert not any(chain.n_terminus_removed or chain.c_terminus_removed for chain in dropped.chains) + + +@pytest.mark.parametrize(("winner", "occupancies"), [("ALA", (0.6, 0.4)), ("UNK", (0.4, 0.6))]) +def test_the_elected_conformation_decides_whether_a_residue_is_unknown(winner, occupancies): + """A position modelled as ALA and as UNK is judged by the amino acid the election keeps.""" + + lines = chain_lines("A", ["SER"]) + lines += residue_lines(len(lines) + 1, "ALA", "A", 2, 7.6, 0.0, altloc="A", occupancy=occupancies[0]) + lines += residue_lines(len(lines) + 1, "UNK", "A", 2, 7.6, 0.0, altloc="B", occupancy=occupancies[1]) + lines += residue_lines(len(lines) + 1, "LYS", "A", 3, 11.4, 0.0) + text = "\n".join(lines) + "\nEND\n" + + if winner == "ALA": + assert residue(parse(text, drop=False), "A", 2).name == "ALA" + else: + with pytest.raises(ValueError, match="UNK A2"): + parse(text, drop=False) + dropped = parse(text, drop=True) + assert dropped.removed_residues == ["UNK A2"] + # Removed, so not also reported as renamed in the run record. + assert dropped.alternate_conformers.summary()["residues_renamed"] == [] + + +def test_a_hetatm_parse_is_not_checked(): + """Ligands and waters are never in the residue table, and are read separately.""" + + text = MIXED.replace("END\n", "") + atom_line(999, "O", "HOH", "W", 1, 0.0, 50.0, 0.0, "O").replace("ATOM ", "HETATM", 1) + "\nEND\n" + + structure = parser.parse_pdb_text(text, "test", identifier="HETATM") + + assert [r.name for r in structure.residues] == ["HOH"] + + +def test_a_predicted_terminal_pka_for_the_new_end_residue_is_ignored_and_reported(caplog): + """A pKa file can give the residue after a removed one an N+ value; it has no such group.""" + + structure = parse(MIXED, drop=True) + + with caplog.at_level(logging.WARNING, logger="prodes.core.structure"): + structure.redo_pkas({"A": {2: [{"N+": 7.5}]}, ANY_CHAIN: {}}) + + assert "N+" not in terminal_groups(residue(structure, "A", 2)) + assert "VAL A2 a N+ pKa, which that residue does not have" in caplog.text + + +def test_the_run_record_lists_the_removed_residues(): + """Empty when none were removed, so the key is always there to read.""" + + import numpy as np + + assert run_metadata("x.pdb", {}, 0, np.array([]))["removed_residues"] == [] + assert run_metadata("x.pdb", {}, 0, np.array([]), removed_residues=["UNK A1"])["removed_residues"] == ["UNK A1"] + + +def test_a_residue_set_loses_the_removed_residues(caplog): + """A residue with no atoms has no surface to restrict a patch to; the rest of the set is kept.""" + + patches = PatchSettings(residue_sets={"Fv": ["A1", "A2", "A3"]}) + + with caplog.at_level(logging.WARNING, logger="prodes.run"): + kept = without_removed_residues(patches, ["UNK A1"]) + + assert kept.residue_sets == {"Fv": ["A2", "A3"]} + assert "residue set Fv loses the residues that were left out as unknown: A1" in caplog.text + assert without_removed_residues(patches, []) is patches + + +def test_a_residue_set_of_nothing_but_removed_residues_is_refused(): + """An empty set would describe nothing, which is not what whoever named it asked for.""" + + with pytest.raises(ValueError, match="residue set CDR names only residues that were left out as unknown"): + without_removed_residues(PatchSettings(residue_sets={"CDR": ["A1"]}), ["UNK A1"]) + + +def test_calculate_runs_by_default_and_records_what_it_left_out(tmp_path): + """End to end: the bundle is written and its run record says what was left out.""" + + pdb = write_structure(tmp_path / "mixed.pdb", MIXED.removesuffix("END\n").splitlines()) + out = tmp_path / "mixed.zip" + + with pytest.raises(ValueError, match="UNK A1"): + calculate(pdb, str(out), drop_unknown_residues=False) + + calculate(pdb, str(out), patches=PatchSettings(residue_sets={"A": ["A1", "A2", "A3"]})) + + with zipfile.ZipFile(out) as bundle: + record = json.loads(bundle.read(next(name for name in bundle.namelist() if name.endswith("prodes_run.json")))) + assert record["removed_residues"] == ["UNK A1", "UNK A4", "UNK A6"] + assert record["settings"]["drop_unknown_residues"] is True + assert record["settings"]["patches"]["residue_sets"] == {"A": ["A2", "A3"]}