#37: leave out residues outside the 20 standard amino acids with a warning, and name them

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 <docs/residue_identity.md>`_ and `alternate conformations <docs/alternate_conformations.md>`_. 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 <docs/pdb_reading.md>`_.
 * **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 <docs/surface_patches.md>`_.
+* **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 <docs/residue_identity.md>`_. 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 <docs/alternate_conformations.md>`_. 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 <docs/disulfide_bond_detection.md>`_. 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.

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
 ```
 

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"}

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":

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

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

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])