From 854c909278aeb944f3a6c63ad8b803bdc8f29bf4 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Thu, 15 Oct 2026 20:25:57 +0200 Subject: [PATCH 1/4] #35 Add contiguous surface patch features with residue-set restriction The default features summarise the surface by the sign of each point's value, so every hydrophobic point is pooled whether the points form one region or many, and every feature covers the whole molecule. Hydrophobic interaction chromatography retention of antibodies has been reported to correlate with the area of contiguous hydrophobic surface, in particular on the CDRs, which no existing feature measures. A new, opt-in patch step (--patches, or calculate(patches=PatchSettings(...))) averages the projected hydrophobicity and the formal-charge potential of each surface point over a neighbourhood (default 5 A, at most 15 A), marks the points past a threshold, joins marked points that touch on the 1 A surface lattice into patches, and discards patches under a minimum area (default 100 A^2, about half the solvent accessible area of one exposed phenylalanine). It reports the summed area, largest area and count of the hydrophobic, positive and negative patches, plus the maximum and the area-weighted positive sum of the local hydrophobicity. --residue-set NAME=FILE reports every one of these again over the residues listed, with patches found on the whole surface and only their in-set share counted. Each surface point now carries the solvent accessible area of its Shrake-Rupley dots and the residue owning most of it, computed while the surface cell still holds those dots. The point table always writes both, and read_surface_points reads the residue as text. The default hydrophobic threshold is refused on a scale with no negative values, residue labels and the threshold are checked before the surface is built, patch options given without --patches are refused, and the patch settings and residue sets are recorded in prodes_run.json. The patch columns are appended after the default ones, so without the flag the feature table is unchanged. The feature dictionary documents the patch columns in a separate file, features_patches.yaml, and resolves a residue-set suffix on patch codes only. The README passages that said Prodes finds no discrete patches and that its hydrophobicity potential has the length scale of SAP are corrected: the exp(-d) weight has a 1 A decay length, which is why the neighbourhood average is needed. Version 9.1.0. --- README.rst | 37 ++- docs/surface_patches.md | 101 ++++++ pyproject.toml | 2 +- src/prodes/__init__.py | 5 + src/prodes/calculations/grid_wizard.py | 52 ++- src/prodes/calculations/patches.py | 357 +++++++++++++++++++++ src/prodes/core/point.py | 7 +- src/prodes/data/features_patches.yaml | 99 ++++++ src/prodes/feature_dictionary.py | 81 ++++- src/prodes/output.py | 33 +- src/prodes/run.py | 118 ++++++- tests/test_cli.py | 65 +++- tests/test_output_bundle.py | 2 +- tests/test_patches.py | 424 +++++++++++++++++++++++++ tests/test_run_prodes.py | 3 + 15 files changed, 1345 insertions(+), 41 deletions(-) create mode 100644 docs/surface_patches.md create mode 100644 src/prodes/calculations/patches.py create mode 100644 src/prodes/data/features_patches.yaml create mode 100644 tests/test_patches.py diff --git a/README.rst b/README.rst index 1c78c53..1660218 100644 --- a/README.rst +++ b/README.rst @@ -228,7 +228,8 @@ Unpacked, the bundle holds: 1GDW/ 1GDW_features.csv the 54 features, one row, first column ID - 1GDW_surface_points.csv every surface point: x, y, z, potential, hydrophobicity + 1GDW_surface_points.csv every surface point: x, y, z, potential, hydrophobicity, + area, and the residue owning it 1GDW_ep.pml PyMOL, coloured by electrostatic potential 1GDW_hydrophobicity.pml PyMOL, coloured by hydrophobicity 1GDW_ep.cxc ChimeraX, electrostatic potential @@ -256,7 +257,7 @@ To read the features back: from prodes.output import read_features, read_surface_points features = read_features("1GDW.zip") - points = read_surface_points("1GDW.zip") # x, y, z, ep_volts, hydrophobicity + points = read_surface_points("1GDW.zip") # x, y, z, ep_volts, hydrophobicity, area_A2, residue What the features actually are ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -264,7 +265,7 @@ What the features actually are Prodes builds a dotted surface over the protein (a Shrake-Rupley solvent-accessible surface), assigns each surface point an electrostatic potential and a hydrophobicity, and then summarises those distributions. The 54 default columns are, broadly: * **Whole-molecule properties** — molecular weight, total surface area, formal charge at your chosen pH, isoelectric point, dipole moment, and shape descriptors. -* **Surface electrostatic potential** — the maximum, minimum, mean and standard deviation of the potential over the whole surface, then the same statistics computed separately over just the positive and just the negative regions, plus a count of how many surface points are positive. A protein with a strongly positive area of surface and a strongly negative one can have a near-zero net charge but very distinctive values here, which is exactly the kind of thing a net-charge calculation misses. Note that these are statistics over all positive or all negative surface points; Prodes does not group them into individual patches (see `Similar tools`_). +* **Surface electrostatic potential** — the maximum, minimum, mean and standard deviation of the potential over the whole surface, then the same statistics computed separately over just the positive and just the negative regions, plus a count of how many surface points are positive. A protein with a strongly positive area of surface and a strongly negative one can have a near-zero net charge but very distinctive values here, which is exactly the kind of thing a net-charge calculation misses. Note that these are statistics over all positive or all negative surface points; grouping them into individual patches is the optional patch step (see `Surface patches`_). * **Surface hydrophobicity** — the same treatment applied to a molecular hydrophobicity potential mapped onto the surface. * **Far-field shell electrostatics** — the potential projected out onto shells around the protein, which captures how the charge distribution looks to another molecule approaching from a distance rather than at contact. * **Per-residue surface fractions** — how much of the accessible surface each amino acid type contributes. @@ -274,6 +275,19 @@ Every feature has a plot label, a one-line description, a full explanation of ho * `features_reduced.yaml `_ — the **54 features calculated by default** * `features_full_only.yaml `_ — the **other 51**, written only under ``--full-features``, each with the reason it is not in the default set +A third file, `features_patches.yaml `_, documents the surface patch features, which are written only with ``--patches``. + +Surface patches +~~~~~~~~~~~~~~~ + +``--patches`` adds a block of features that group the surface into contiguous regions. Each point's hydrophobicity and potential are averaged over a 5 Å neighbourhood, points above a threshold that touch are joined into patches, and the run reports the summed area, the largest area and the number of hydrophobic, positive and negative patches, plus the highest and the summed positive local hydrophobicity. ``--residue-set NAME=FILE`` reports every one of them again restricted to the residues listed in the file, such as the CDRs of an antibody, under columns ending in ``_NAME``: + +.. code-block:: bash + + prodes Fv.pdb Fv.zip --patches --residue-set CDR=cdr_residues.txt + +The default 54 columns are unchanged by the flag; the patch columns are appended after them. The thresholds and the minimum patch area are settings with starting values rather than calibrated choices. Method, settings, columns and limitations: `surface patches `_. + Looking features up from Python ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -633,7 +647,7 @@ The command line and the Python interface do the same work. To reproduce a comma import prodes prodes.run_prodes("./tests/data/1GDW.pdb.zip", "example.zip") -The full signature is ``run_prodes(pdb_file, out_file, pkas_file=None, ph=7, r_probe=1.4, hydro_scale="mj_scaled", full_features=False, mem_limit_mb=None)``. +The full signature is ``run_prodes(pdb_file, out_file, pkas_file=None, ph=7, r_probe=1.4, hydro_scale="mj_scaled", full_features=False, mem_limit_mb=None, ionic_strength_molar=0.15, patches=None)``, where ``patches`` takes a ``prodes.calculations.patches.PatchSettings``; see `Surface patches`_. The lower-level pieces are importable too, if you want a single property rather than the whole feature set. Calculating just the surface area, for example: @@ -784,7 +798,7 @@ In industry this category is dominated by three licensed packages, and the antib * **Discovery Studio** (`BIOVIA, Dassault Systèmes `_) — the ``Calculate Protein Features`` protocol (structure-based) and ``Calculate Sequence Descriptors`` (sequence-based, the only genuinely sequence-only protocol of the three), plus named developability predictors: the Developability Index built on **SAP**, the spatial charge map **SCM** for viscosity, and solubility and pI calculators. -Two of the methods these suites ship are the conceptual ancestors of most of the open tools above. **SAP** (Chennamsetty et al., *PNAS* 2009, `DOI: 10.1073/pnas.0904191106 `_) established averaging a hydrophobicity scale over a sphere around each atom as the standard way to score a protein surface, and **AggScore** (Sankar et al., *Proteins* 2018, `DOI: 10.1002/prot.25594 `_) extended it to the distribution of hydrophobic and charged patches. Neither has a canonical free implementation. Prodes' molecular hydrophobicity potential — an ``exp(-d)``-weighted sum over non-hydrogen atoms within 10 Å — belongs to the same family. +Two of the methods these suites ship are the conceptual ancestors of most of the open tools above. **SAP** (Chennamsetty et al., *PNAS* 2009, `DOI: 10.1073/pnas.0904191106 `_) established averaging a hydrophobicity scale over a sphere around each atom as the standard way to score a protein surface, and **AggScore** (Sankar et al., *Proteins* 2018, `DOI: 10.1002/prot.25594 `_) extended it to the distribution of hydrophobic and charged patches. Neither has a canonical free implementation. Prodes' molecular hydrophobicity potential — an ``exp(-d)``-weighted sum over non-hydrogen atoms within 10 Å — belongs to the same family at a much shorter length scale; the 5 Å neighbourhood average of the patch step (see `Surface patches`_) is the closer match. Related, but not alternatives ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -795,7 +809,7 @@ These two come up constantly in this literature, but neither competes with Prode * `NanoShaper `_ / `NanoShaperWeb `_ (Decherchi & Rocchia, IIT) — **GPL-3.0**, so fine to run as a standalone binary but copyleft if you link it into a product you distribute. Two distinct roles, and it is worth keeping them apart: - * **As a surface engine**, it triangulates the molecular surface analytically by ray casting rather than on a grid, and this is the mesh layer beneath both PROPERMAB and Genentech's descriptors. When those papers report a patch area in Ų, the number rests on a NanoShaper triangulation. Prodes does not need it, because it uses a dotted Shrake-Rupley solvent-accessible surface — a point cloud with no per-point area, which is why Prodes reports point counts where mesh-based tools report areas. + * **As a surface engine**, it triangulates the molecular surface analytically by ray casting rather than on a grid, and this is the mesh layer beneath both PROPERMAB and Genentech's descriptors. When those papers report a patch area in Ų, the number rests on a NanoShaper triangulation. Prodes does not need it, because it uses a dotted Shrake-Rupley solvent-accessible surface. Its default features report point counts where mesh-based tools report areas; the patch features give each point the solvent-accessible area of its dots, which is larger than the molecular-surface area a mesh tool reports for the same region. * **As a descriptor generator**, it is aimed at **pockets and cavities rather than the outer surface**. Pockets are found from the volumetric difference between two solvent-excluded surfaces built with different probe radii, and NanoShaperWeb then characterises each pocket with the DrugPred druggability descriptor set — volume, entrance area, compactness, donor/acceptor/hydrophobic surface fractions. That is a small-molecule binding-site question — is this cavity druggable? — not the whole-surface charge and hydrophobicity profiling that Prodes and every other tool listed above is built for. The two are not substitutes in either direction. Decherchi & Rocchia, *PLoS ONE* 2013, `DOI: 10.1371/journal.pone.0059744 `_; Abate et al., *J. Chem. Inf. Model.* 2025, `DOI: 10.1021/acs.jcim.5c00821 `_. @@ -803,11 +817,11 @@ These two come up constantly in this literature, but neither competes with Prode Where Prodes sits ~~~~~~~~~~~~~~~~~~ -**Currently, Prodes does not identify discrete patches.** It builds a dotted surface, gives every point an electrostatic potential and a hydrophobicity, and then summarises those values over the whole surface and over the positive and negative subsets. Those subsets are defined by the sign of the value at a point, not by spatial contiguity: a positive-subset statistic pools every positive surface point on the protein, whether they form one large region or twenty scattered ones. +**By default, Prodes does not identify discrete patches.** It builds a dotted surface, gives every point an electrostatic potential and a hydrophobicity, and then summarises those values over the whole surface and over the positive and negative subsets. Those subsets are defined by the sign of the value at a point, not by spatial contiguity: a positive-subset statistic pools every positive surface point on the protein, whether they form one large region or twenty scattered ones. -So if you need the area of the single largest hydrophobic patch, or the list of residues making up a particular patch so you can map it back onto the structure, use PEP-Patch. Prodes cannot give you that. +With ``--patches`` it does: neighbourhood-averaged values are thresholded and touching points joined into patches, whose areas are reported over the whole surface and over any residue set you name, and the point table records which hydrophobic patch each point belongs to and which residue owns it. See `Surface patches`_. The patches are found on a solvent-accessible point cloud rather than a triangulated molecular surface, so their areas are not directly comparable with those of mesh-based tools such as PEP-Patch. -It is worth being clear about what this does *not* mean, though, because the word "patch" is used very loosely in this field. **Prodes is not blind to where things sit on the surface.** Its molecular hydrophobicity potential gives each surface point an ``exp(-d)``-weighted sum over the non-hydrogen atoms within 10 Å, which is the same construction as SAP, as Genentech's HPATCH (a 10 Å average) and as Protein-Sol's NPP ratio (a 13 Å ratio). Only three of the tools above actually segment the surface into discrete regions: PEP-Patch (connected components on the surface graph), PROPERMAB (DBSCAN on mesh triangles) and MOE's Protein Patch Analyzer. The rest, TAP and TNP included, compute neighbourhood-weighted quantities and then aggregate them, which is what Prodes does too. The difference is in the last step: they take a maximum or a region sum, Prodes takes the distribution statistics. +It is worth being clear about what this does *not* mean, though, because the word "patch" is used very loosely in this field. **Prodes is not blind to where things sit on the surface.** Its molecular hydrophobicity potential gives each surface point an ``exp(-d)``-weighted sum over the non-hydrogen atoms within 10 Å. The weight has a decay length of 1 Å, so the 10 Å cutoff hardly matters and the value mostly reflects the one or two nearest atoms, a much shorter length scale than SAP, Genentech's HPATCH (a 10 Å average) or Protein-Sol's NPP ratio (a 13 Å ratio). The patch step's 5 Å neighbourhood average is what brings it to the scale of a surface region. Only three of the tools above actually segment the surface into discrete regions: PEP-Patch (connected components on the surface graph), PROPERMAB (DBSCAN on mesh triangles) and MOE's Protein Patch Analyzer. The rest, TAP and TNP included, compute neighbourhood-weighted quantities and then aggregate them, which is what Prodes does too. The difference is in the last step: they take a maximum or a region sum, the default Prodes features take the distribution statistics, and the patch features take a maximum and a sum, over the whole surface or a residue set. What Prodes offers, then, is: @@ -821,7 +835,7 @@ Three things Prodes does **not** do, which the better tools here do: * **No conformational averaging.** Genentech's main finding is that descriptors taken from a single static structure are unstable, which is why they average over a 5 ns accelerated-MD ensemble; MOE averages over a LowModeMD ensemble, and Aggrescan3D offers a coarse-grained dynamic mode. Prodes computes from whatever one structure you give it, and inherits that instability. You can approximate the ensemble average by running Prodes over several structures and averaging yourself, but nothing in the package does it for you. * **Almost no structure preparation of its own.** The commercial suites correct the input before they describe it — resolve alternate conformations, rebuild missing side chains, cap chain breaks, form disulfide bonds, optimise the hydrogen-bond network. Prodes does two of those and none of the rest. It *recognises* existing disulfide bonds, so that a cystine is not titrated as a free thiol, and it *resolves alternate conformations* by keeping the best occupied one per residue, which it does more reliably than PDB2PQR. It does not form bonds, rebuild side chains, cap chain breaks or touch any coordinate. The free fix is to run PDB2PQR first, which is why that is now the recommended workflow; see `Preparing your structure`_. -* **No region restriction, and no pockets.** Every Prodes feature covers the whole molecule. There is no way to compute over a CDR, a domain or an interface, and nothing equivalent to NanoShaper's cavity and pocket detection. +* **Little region restriction, and no pockets.** Every default feature covers the whole molecule. Only the patch features can be restricted to a residue set, such as a CDR, a domain or an interface, and there is nothing equivalent to NanoShaper's cavity and pocket detection. A caveat that applies to every tool in this list, Prodes included, and is made forcefully in both the Genentech and PROPERMAB papers: descriptor values are sensitive to the structure model, the protonation assignment and the software version, and different packages computing nominally the same quantity often disagree. Treat any of these numbers as reproducible only within one pipeline held fixed. @@ -839,6 +853,7 @@ Improvements in this fork * **Insertion codes and NMR models** — from version 8.0 a residue is its chain, its number **and** its insertion code, and a file holding several models is described by its first. Kabat and Chothia numbering write every antibody CDR insertion as an insertion code, and reading ``H100`` and ``H100A`` as one residue put both side chains in one residue: on the 4NZU Fab that lost 12 of its 434 residues, moved the molecular weight by 1.3 kDa, titrated a lysine and a cysteine against aspartate's pKa, and hid a cysteine well enough that one of its five disulfide bonds went undetected. An NMR ensemble read whole was worse: 1PIT came out as 58 residues holding 17 780 atoms, 16 901 of them in a single residue, at a formal charge of -1096 for a protein whose real charge at pH 7 is about +6. See `what Prodes calls one residue `_. * **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 9.1, ``--patches`` groups neighbourhood-averaged hydrophobicity and potential into contiguous patches and reports their areas over the whole surface and over named residue sets, such as antibody CDRs. Off by default; the default columns do not change. Every surface point in the bundle now also carries its area and the residue owning it. See `surface patches `_. * **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. @@ -849,6 +864,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 9.1 the surface point table in the bundle has two more columns, ``area_A2`` and ``residue``, on every run. The feature table is unchanged unless ``--patches`` is given, which appends the patch columns after the others. + 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 54 columns they share. How to cite diff --git a/docs/surface_patches.md b/docs/surface_patches.md new file mode 100644 index 0000000..c6e2d66 --- /dev/null +++ b/docs/surface_patches.md @@ -0,0 +1,101 @@ +# Surface patches + +The default Prodes features summarise the surface by the sign of each point's value: every hydrophobic point on the protein is pooled, whether the points form one region or twenty scattered ones. The patch step adds the spatial grouping. It reports how much of the surface is covered by contiguous hydrophobic, positive and negative regions, how large the largest one is, and the same quantities restricted to a named set of residues, such as the CDRs of an antibody. + +The step is off by default. Switching it on appends columns to the feature table and leaves every default column unchanged. + +## Running it + +Command line: + +```bash +prodes Fv.pdb Fv.zip --patches --residue-set CDR=cdr_residues.txt --residue-set Fv=fv_residues.txt +``` + +Python: + +```python +from prodes.run import calculate +from prodes.calculations.patches import PatchSettings + +calculate("Fv.pdb", "Fv.zip", patches=PatchSettings(residue_sets={"CDR": ["H27", "H28", "H100A", "L91"]})) +``` + +A residue set file holds one residue label per line. A label is the chain, the residue number and the insertion code, written as `H100A` or `L91`; anything after a `#` is ignored. A label the structure does not contain stops the run with an error naming it, before the surface is built. Which residues form the CDRs depends on the numbering scheme of the structure and is decided by the caller, not by Prodes. + +## Method + +| Step | What happens | +| --- | --- | +| Per-point value | Each point of the 1 A surface grid carries the projected hydrophobicity and the formal-charge electrostatic potential that the default features are calculated from. | +| Per-point area and residue | Each point is one occupied 1 A cell of the Shrake-Rupley solvent accessible surface. Its area is the summed area of the surface dots in the cell, weighted as `Atom.surface_area` weights them, so the point areas add up to `Area` apart from internal cavities (see Limitations). Its residue is the residue whose dots carry the most area in the cell. | +| Local value | Each point's value is replaced by the area-weighted mean over every surface point within the radius, the point itself included. | +| Threshold | A point is hydrophobic when its local hydrophobicity is above the hydrophobic threshold, positive when its local potential is above the electrostatic threshold, and negative when it is below the negative of the electrostatic threshold. | +| Patch | Marked points that touch on the surface lattice (any of the 26 neighbours, the connectivity the surface itself is built with) form one patch. Patches smaller than the minimum area are discarded. | +| Restriction | Patches are found on the whole surface. A restricted column counts only the area of patch points owned by a residue in the set. | + +The local average is needed because the per-point projected hydrophobicity is short-ranged. Each atom contributes with a weight of `exp(-d)`, a decay length of 1 A: 0.05 at 3 A and 4.5e-5 at 10 A. A surface point sits about 3 A from the nearest atom centres, so its value mostly reflects the one or two atoms closest to it. Averaging over 5 A gives the value the length scale of a surface region. Waibl et al. (Front. Mol. Biosci. 2022) compared averaging radii for hydrophobicity descriptors of antibodies and found 5 A better than 10 A in almost every case. + +Finding patches on the whole surface and then restricting them means a patch that crosses the border of a residue set is one patch, and the minimum area is judged on the whole of it. Restricting first would cut patches at the border and make the minimum area depend on where the border falls. + +## Settings + +| Setting | CLI | Default | Notes | +| --- | --- | --- | --- | +| `radius` | `--patch-radius` | 5.0 A | Radius of the local average, at most 15 A. The cost of the average grows with the cube of the radius. | +| `min_area` | `--patch-min-area` | 100.0 A^2 | Smallest patch kept. A starting value; see below. | +| `hydrophobic_threshold` | `--patch-hydrophobic-threshold` | 0.0 | The sign of the local hydrophobicity. Refused for a scale with no negative values, where 0 would mark every point hydrophobic; pass a threshold explicitly there. | +| `ep_threshold` | `--patch-ep-threshold` | 0.5 V | Applied symmetrically to positive and negative patches. | +| `residue_sets` | `--residue-set NAME=FILE` | none | Repeatable. `NAME` is letters and digits only. | + +The `--patch-*` and `--residue-set` options are refused without `--patches`, since they would otherwise do nothing. + +The hydrophobicity scale is the run's `--hydro` setting. The thresholds are in the units of the local value, so a threshold chosen for one scale does not carry over to another. + +The defaults are starting values to be tuned against measured data, not calibrated choices: + +- **Minimum area.** Areas here are solvent accessible surface areas with a 1.4 A probe, which are larger than molecular surface areas of the same region, and by an amount that depends on curvature, so a minimum area quoted by a tool that measures on a molecular surface does not transfer to this one. 100 A^2 is about half the solvent accessible area of one fully exposed phenylalanine residue, 210 to 240 A^2 (Tien et al., PLoS ONE 2013): a patch has to be larger than the exposed half of one large hydrophobic side chain. +- **Electrostatic threshold.** At 0 V, the surface with the same sign as the net charge becomes one or two patches covering most of the protein, and the patch area then says little beyond `NSurfPosEpFormal` or the total area. 0.5 V was chosen from the local potential at the default ionic strength of 0.15 M on the eight structures in `tests/data` (1AO6, 1CBN, 1GDW, 1GPB, 1PIT, 4NZU, ARH96693 and ARH98503): over those structures the median of the 75th percentile of the absolute local potential was 0.39 V and the median of the 90th percentile 0.61 V. At 0.5 V the marked points covered between 13 and 26 per cent of the surface area on six of the eight structures, 2 per cent on 1CBN and 32 per cent on 1PIT. The potential is a screened Coulomb sum at a relative permittivity of 4, so thresholds from Poisson-Boltzmann tools, usually given in kT/e, do not carry over. +- **Hydrophobic threshold.** 0 on `mj_scaled` separates hydrophobic from hydrophilic residues. It is lower than the 0.03 the hydrophobicity viewer uses to colour points, because the viewer shows the per-point value and the patches use the local value, whose range is narrower. The local value on `mj_scaled` is negative over most of a typical surface, so at 0 only the few most hydrophobic regions pass, and the hydrophobic patch area changes steeply with the threshold. Choose it against measured data rather than relying on the default. + +All settings, with the hydrophobic threshold resolved and every residue set's labels, are recorded in `prodes_run.json` under `settings.patches`. + +## Columns + +The whole-surface columns are, in output order: + +| Code | Meaning | +| --- | --- | +| `SurfMhpPatchArea` | summed area of the hydrophobic patches, A^2 | +| `SurfMhpPatchMaxArea` | area of the largest hydrophobic patch, A^2 | +| `NSurfMhpPatch` | number of hydrophobic patches | +| `SurfMhpLocalMax` | highest local hydrophobicity | +| `SurfMhpLocalPosSum` | sum over points of area times the positive part of the local hydrophobicity, A^2 | +| `SurfEpPosPatchArea`, `SurfEpPosPatchMaxArea`, `NSurfEpPosPatch` | the same three for positive patches | +| `SurfEpNegPatchArea`, `SurfEpNegPatchMaxArea`, `NSurfEpNegPatch` | the same three for negative patches | + +For each residue set `NAME` the run then writes `SurfArea_NAME`, the surface area owned by the set, followed by every code above with `_NAME` appended. Within a set, `NSurf*Patch_NAME` counts the patches with at least one point in the set, and `*PatchMaxArea_NAME` is the largest share of one patch that lies in the set, not the largest patch lying entirely inside it. A framework patch that only touches the CDRs contributes its small CDR share. + +`SurfMhpLocalMax` and `SurfMhpLocalPosSum` need no threshold and no minimum area, so they are not affected by either choice. They are similar in spirit to spatial aggregation propensity (Chennamsetty et al., PNAS 2009), which also averages a hydrophobicity scale over a neighbourhood and sums the positive part, although that method works per atom and weights each residue by its exposed side-chain fraction. `SurfMhpLocalPosSum` is in square angstrom times the units of the hydrophobicity scale, so its value depends on the scale. + +Every code is documented in `src/prodes/data/features_patches.yaml`, and `FeatureDictionary` looks up a code with or without a residue set suffix. + +The point table, `*_surface_points.csv`, always carries `area_A2` and `residue` for each point. With the patch step it also carries `local_mhp`, `local_ep_volts` and `mhp_patch_id`, the hydrophobic patch the point belongs to (-1 for none), so a patch can be traced back to its residues. + +## Limitations + +- **Truncated structures.** On an Fv or Fab cut out of a larger molecule, the surface that faces the removed domains is exposed only because of the cut, and the terminal oxygen (`OXT`) contributes 1.0, the top of the default `mj_scaled` scale; see [the terminal oxygen](structure_preparation.md#what-happens-to-a-predicted-structure-the-terminal-oxygen). Both can form or enlarge whole-surface patches. Columns restricted to a residue set away from the cut, such as the CDRs, are not affected. +- **Clefts.** The neighbourhood is a Euclidean sphere, so at 5 A it can reach across a narrow cleft to the facing surface. +- **Internal cavities and separate surfaces.** The surface grid is built by a flood fill from one point, so a surface component that does not touch the rest, such as an internal cavity or a second molecule that does not touch the first, is not part of it. This applies to every surface feature, not only the patches, and is why the point areas can add up to less than `Area`. +- **Residues with little exposed surface.** A surface point belongs to the residue whose dots carry the most area in its 1 A cell, so a residue with only a sliver of exposed surface can own no point. A residue set none of whose residues owns a point stops the run with an error rather than reporting zeros. +- **Memory.** The local average places the points on a dense array the size of the structure's bounding box, two arrays of 8 bytes per cubic angstrom. This is a few MB for an Fv and a few tens of MB for a whole antibody, and it is not counted against `--mem-limit`. +- **Electrostatic patches and chromatography.** In hydrophobic interaction chromatography the salt concentration is high and the Debye length is a few angstrom, against the 0.15 M Prodes uses by default. The electrostatic patches describe charge distribution at the run's ionic strength and are not intended as hydrophobic interaction descriptors. +- **One structure.** As with every Prodes feature, the values come from one static structure. Patch areas are known to vary between structure models of the same molecule. + +## References + +- Chennamsetty et al., *PNAS* 2009, DOI: 10.1073/pnas.0904191106. Spatial aggregation propensity. +- Waibl et al., *Front. Mol. Biosci.* 2022. Hydrophobicity scales and averaging radii for antibody descriptors. +- Hoerschinger et al., *J. Chem. Inf. Model.* 2023, DOI: 10.1021/acs.jcim.3c01490. PEP-Patch, connected-component surface patches. +- Li et al., *mAbs* 2025, 17:2474521. PROPERMAB, patch descriptors on the Fv with CDR-localised versions. +- Tien et al., *PLoS ONE* 2013, DOI: 10.1371/journal.pone.0080635. Maximum solvent accessible surface areas of the amino acids. diff --git a/pyproject.toml b/pyproject.toml index b45687c..0addc1c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "prodes" -version = "9.0.0" +version = "9.1.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 4d65fa2..4b217e6 100644 --- a/src/prodes/__init__.py +++ b/src/prodes/__init__.py @@ -13,6 +13,7 @@ def run_prodes( full_features=False, mem_limit_mb=None, ionic_strength_molar=run.DEFAULT_IONIC_STRENGTH_MOLAR, + patches=None, ): """Runs Prodes similar as the commandline tool. @@ -36,6 +37,9 @@ def run_prodes( mem_limit_mb: memory budget in MB for the intermediate NumPy arrays of this run, divided among the workers. Pass explicitly in worker processes rather than relying on the PRODES_MEM_LIMIT_MB environment variable. + patches: a prodes.calculations.patches.PatchSettings to also report + contiguous surface patches, optionally over named residue sets. + None leaves them out. """ calculate( @@ -48,6 +52,7 @@ def run_prodes( full_features=full_features, mem_limit_mb=mem_limit_mb, ionic_strength_molar=ionic_strength_molar, + patches=patches, ) diff --git a/src/prodes/calculations/grid_wizard.py b/src/prodes/calculations/grid_wizard.py index 690dba9..ce5f2c9 100644 --- a/src/prodes/calculations/grid_wizard.py +++ b/src/prodes/calculations/grid_wizard.py @@ -138,8 +138,52 @@ def grid_content(self, *content_filter, cells=None): return content -def property_points_on_surface(grid, to_consider="Surface_point"): - """Finds the property points on a surface""" +def surface_cell_area_and_residue(cell, probe_r, dot_weights): + """Returns the area a cell's Shrake-Rupley dots stand for, and the residue owning most of it. + + Each dot weighs what Atom.surface_area gives it, the atom's maximal + accessible area divided by its maximal number of dots, so the areas of all + the cells add up to Structure.surface_area over the same dots. The weight + is cached per atom in dot_weights, since a large protein has some 200,000 + dots and a few thousand atoms. + + The residue is the one whose dots carry the most area in the cell. On an + exact tie the alphabetically first label wins, so the assignment does not + depend on the order the dots were filed in. + + Returns: + area: summed dot area in square angstrom. + residue: the winning residue's label, such as H100A, or "" when the + dots carry no residue. + """ + + area = 0.0 + by_residue: dict = {} + for dot in cell.filtered_content("Surface_point"): + atom = dot.atom + weight = dot_weights.get(id(atom)) + if weight is None: + weight = atom.max_area(probe_r=probe_r) / atom.max_n_points(probe_r=probe_r) + dot_weights[id(atom)] = weight + area += weight + + label = atom.residue.label if atom.residue is not None else "" + by_residue[label] = by_residue.get(label, 0.0) + weight + + residue = min(by_residue.items(), key=lambda item: (-item[1], item[0]))[0] if by_residue else "" + + return area, residue + + +def property_points_on_surface(grid, to_consider="Surface_point", probe_r=1.4): + """Finds the property points on a surface. + + Each property point is one occupied cell, reached by a flood fill from the + first occupied cell through occupied neighbours. A surface component that + never touches that first one, such as an internal cavity, is not reached. + Each point carries the area its cell's dots stand for and the residue that + owns most of it; see surface_cell_area_and_residue. + """ from prodes.core.point import Property_point @@ -155,12 +199,14 @@ def property_points_on_surface(grid, to_consider="Surface_point"): surface_points = [] done_ids = set() todo_cells = [start] + dot_weights: dict = {} while todo_cells: next_todo = [] next_todo_ids = set() for cell in todo_cells: - surface_points.append(Property_point(cell.x, cell.y, cell.z)) + area, residue = surface_cell_area_and_residue(cell, probe_r, dot_weights) + surface_points.append(Property_point(cell.x, cell.y, cell.z, area=area, residue=residue)) done_ids.add(id(cell)) surrounding = grid.find_surrounding_cells(cell) diff --git a/src/prodes/calculations/patches.py b/src/prodes/calculations/patches.py new file mode 100644 index 0000000..5ab7c09 --- /dev/null +++ b/src/prodes/calculations/patches.py @@ -0,0 +1,357 @@ +"""Finds contiguous patches of hydrophobic, positive and negative surface, and reports them over named residue sets. + +The per-point values the surface grid phase computes (projected hydrophobicity +and the formal-charge electrostatic potential) are summarised elsewhere over the +whole surface by sign alone. This module adds the spatial step: it averages each +value over a neighbourhood of the point, marks the points whose average passes a +threshold, links marked points that touch on the 1 A surface lattice into +patches, and reports the area of those patches. Every descriptor can also be +reported restricted to a named set of residues, such as the CDRs of an antibody. + +The method and the reasons for each setting are in docs/surface_patches.md. + +Patches are found on the whole surface and only then restricted. A patch that +crosses the border of a residue set is one patch, and the minimum area applies +to the whole of it; the restricted columns count only its share inside the set. +""" + +import re +from collections import deque +from dataclasses import dataclass, field +from itertools import product + +import numpy as np + +#: Neighbourhood radius in angstrom over which each point's value is averaged. +DEFAULT_PATCH_RADIUS = 5.0 + +#: Patches smaller than this, in square angstrom of solvent accessible surface, +#: are discarded: about half the solvent accessible area of one fully exposed +#: phenylalanine. A starting value; see docs/surface_patches.md. +DEFAULT_PATCH_MIN_AREA = 100.0 + +#: The largest radius accepted. The neighbourhood average reads the whole point +#: set once per lattice offset within the radius, and the offsets grow with its +#: cube: 515 at 5 A, 14,147 at 15 A. +MAX_PATCH_RADIUS = 15.0 + +#: A point is hydrophobic when its neighbourhood-averaged hydrophobicity is +#: above this. 0 is the sign of the projected hydrophobicity, which only means +#: something on a scale that has negative values. +DEFAULT_HYDROPHOBIC_THRESHOLD = 0.0 + +#: A point is positive when its neighbourhood-averaged potential is above this +#: many volts, and negative when below minus this. See docs/surface_patches.md +#: for how the value was chosen. +DEFAULT_EP_THRESHOLD = 0.5 + +RESIDUE_SET_NAME = re.compile(r"^[A-Za-z0-9]+$") + + +@dataclass(frozen=True) +class PatchSettings: + """The settings of the patch step. Passing one to calculate switches the step on. + + Attributes: + radius: neighbourhood radius in angstrom for the average. + min_area: smallest patch kept, in square angstrom. + hydrophobic_threshold: neighbourhood-averaged hydrophobicity above which + a point is hydrophobic. None means DEFAULT_HYDROPHOBIC_THRESHOLD, + which is refused for a scale with no negative values. + ep_threshold: neighbourhood-averaged potential in volts above which a + point is positive and below the negative of which it is negative. + residue_sets: name to residue labels (such as H100A), each reported as + its own set of restricted columns. + """ + + radius: float = DEFAULT_PATCH_RADIUS + min_area: float = DEFAULT_PATCH_MIN_AREA + hydrophobic_threshold: float | None = None + ep_threshold: float = DEFAULT_EP_THRESHOLD + # Left out of the hash, since a dict cannot be hashed; equal settings still + # hash equal, which is all a hash has to guarantee. + residue_sets: dict = field(default_factory=dict, hash=False) + + def __post_init__(self): + for name, value in (("radius", self.radius), ("min_area", self.min_area), ("ep_threshold", self.ep_threshold)): + if isinstance(value, bool) or not isinstance(value, int | float) or not np.isfinite(value) or value < 0: + raise ValueError(f"patch {name} must be a finite number and not negative, got {value!r}") + if not 1 <= self.radius <= MAX_PATCH_RADIUS: + raise ValueError(f"patch radius must be between the 1 A lattice spacing and {MAX_PATCH_RADIUS} A, got {self.radius}") + threshold = self.hydrophobic_threshold + if threshold is not None and (isinstance(threshold, bool) or not isinstance(threshold, int | float) or not np.isfinite(threshold)): + raise ValueError(f"patch hydrophobic_threshold must be a finite number or None, got {threshold!r}") + if not isinstance(self.residue_sets, dict): + raise ValueError(f"residue_sets must be a dict of name to residue labels, got {type(self.residue_sets).__name__}") + for name, labels in self.residue_sets.items(): + if not isinstance(name, str) or not RESIDUE_SET_NAME.match(name): + raise ValueError(f"residue set name {name!r} must be letters and digits only, since it becomes part of a column name") + if isinstance(labels, str) or not all(isinstance(label, str) for label in labels): + raise ValueError(f"residue set {name} must be a list of residue labels such as ['H95', 'H100A'], got {labels!r}") + if len(labels) == 0: + raise ValueError(f"residue set {name} is empty") + + def resolved_hydrophobic_threshold(self, hydro_scale): + """Returns the hydrophobic threshold to use with this scale. + + Raises: + ValueError: when no threshold was given and the scale has no + negative values, so that 0 would mark every point hydrophobic. + """ + + if self.hydrophobic_threshold is not None: + return self.hydrophobic_threshold + + from prodes.data import hydrophobic_scale + + if min(hydrophobic_scale(hydro_scale).values()) >= 0: + raise ValueError( + f"the hydrophobicity scale {hydro_scale} has no negative values, so the default patch threshold of " + f"{DEFAULT_HYDROPHOBIC_THRESHOLD} would mark every point hydrophobic; pass a hydrophobic threshold explicitly" + ) + + return DEFAULT_HYDROPHOBIC_THRESHOLD + + def record(self, hydro_scale): + """Returns the settings as written into prodes_run.json, with the threshold resolved.""" + + return { + "radius_A": self.radius, + "min_area_A2": self.min_area, + "hydrophobic_threshold": self.resolved_hydrophobic_threshold(hydro_scale), + "ep_threshold_volts": self.ep_threshold, + "residue_sets": {name: list(labels) for name, labels in self.residue_sets.items()}, + } + + +def lattice_indices(coords): + """Returns the integer lattice position of each point of the 1 A surface grid. + + The property points are cell centres of a grid of 1 A cells, so their + offsets from the lowest corner are whole numbers up to rounding. Working on + the integers makes every neighbour and radius test exact. + + Raises: + ValueError: if the points are not on a 1 A lattice, or two share a position. + """ + + offsets = coords - coords.min(axis=0) + indices = np.rint(offsets).astype(np.int64) + if not np.allclose(offsets, indices, atol=1e-6): + raise ValueError("the surface points are not on a 1 A lattice") + if len(np.unique(indices, axis=0)) != len(indices): + raise ValueError("two surface points share a lattice position") + + return indices + + +def neighbourhood_offsets(radius): + """Returns the integer offsets whose length is at most radius, the origin included.""" + + reach = int(np.floor(radius)) + steps = np.arange(-reach, reach + 1) + offsets = np.array(list(product(steps, steps, steps))) + + return offsets[(offsets**2).sum(axis=1) <= radius**2] + + +def local_average(indices, values, areas, radius): + """Returns each point's area-weighted mean value over the points within radius of it. + + The point itself is included. Built by placing area and area x value on a + dense array and reading it back at every offset within the radius, which is + a few hundred gathers of the whole point set rather than a pairwise distance + matrix. + """ + + offsets = neighbourhood_offsets(radius) + pad = int(np.floor(radius)) + shape = tuple(indices.max(axis=0) + 1 + 2 * pad) + + area_grid = np.zeros(shape) + weighted_grid = np.zeros(shape) + x, y, z = (indices + pad).T + area_grid[x, y, z] = areas + weighted_grid[x, y, z] = areas * values + + area_sum = np.zeros(len(indices)) + weighted_sum = np.zeros(len(indices)) + for dx, dy, dz in offsets: + area_sum += area_grid[x + dx, y + dy, z + dz] + weighted_sum += weighted_grid[x + dx, y + dy, z + dz] + + return weighted_sum / area_sum + + +def label_patches(indices, mask): + """Returns the patch index of each point, or -1 for a point outside the mask. + + Two marked points belong to the same patch when they touch on the lattice, + which is any of the 26 neighbours, the same connectivity the flood fill that + builds the surface uses. Patches are numbered in the order their first point + appears. + """ + + position = {tuple(index): row for row, index in enumerate(indices) if mask[row]} + steps = [step for step in product((-1, 0, 1), repeat=3) if step != (0, 0, 0)] + + labels = np.full(len(indices), -1, dtype=np.int64) + patch = 0 + for row in range(len(indices)): + if not mask[row] or labels[row] >= 0: + continue + + labels[row] = patch + queue = deque([row]) + while queue: + ix, iy, iz = indices[queue.popleft()] + for dx, dy, dz in steps: + neighbour = position.get((ix + dx, iy + dy, iz + dz)) + if neighbour is not None and labels[neighbour] < 0: + labels[neighbour] = patch + queue.append(neighbour) + patch += 1 + + return labels + + +def kept_patches(labels, areas, min_area): + """Returns labels with every patch smaller than min_area set to -1.""" + + if not np.any(labels >= 0): + return labels + + patch_areas = np.bincount(labels[labels >= 0], weights=areas[labels >= 0]) + small = patch_areas < min_area + kept = labels.copy() + kept[(labels >= 0) & small[np.maximum(labels, 0)]] = -1 + + return kept + + +def patch_statistics(name, labels, areas, in_set): + """Returns the summed area, largest area and count of the kept patches inside a set. + + Args: + name: the population prefix, such as SurfMhp. + labels: kept patch index per point, -1 for none. + areas: area per point. + in_set: boolean per point, True for the points counted. + """ + + counted = (labels >= 0) & in_set + if not np.any(counted): + return {f"{name}PatchArea": 0.0, f"{name}PatchMaxArea": 0.0, f"N{name}Patch": 0} + + per_patch = np.bincount(labels[counted], weights=areas[counted]) + present = per_patch > 0 + + return { + f"{name}PatchArea": round(float(per_patch.sum()), 3), + f"{name}PatchMaxArea": round(float(per_patch.max()), 3), + f"N{name}Patch": int(np.count_nonzero(present)), + } + + +def hydrophobic_strength(local_mhp, areas, in_set): + """Returns the maximum and the area-weighted sum of the positive part of the local hydrophobicity.""" + + return { + "SurfMhpLocalMax": round(float(local_mhp[in_set].max()), 3), + "SurfMhpLocalPosSum": round(float(np.sum(np.maximum(local_mhp[in_set], 0) * areas[in_set])), 3), + } + + +def check_residue_sets(residue_sets, structure_labels): + """Raises ValueError if a residue set names a residue the structure does not have. + + Cheap, so calculate runs it as soon as the structure is read rather than + after the surface has been built. + """ + + known = set(structure_labels) + for name, labels in residue_sets.items(): + missing = sorted(set(labels) - known) + if missing: + raise ValueError(f"residue set {name} names residues the structure does not have: {', '.join(missing)}") + + +def residue_set_masks(residue_sets, point_residues, structure_labels): + """Returns a boolean point mask per residue set. + + A surface point belongs to the residue whose dots carry the most area in + its 1 A cell, so a residue with only a sliver of exposed surface can own + no point at all. + + Raises: + ValueError: if a set names a residue the structure does not have, or + none of its residues owns a point of the surface. + """ + + check_residue_sets(residue_sets, structure_labels) + + masks = {} + for name, labels in residue_sets.items(): + mask = np.isin(point_residues, list(labels)) + if not np.any(mask): + raise ValueError( + f"none of the residues of set {name} owns a surface point: each is buried, exposed too little to own one " + "of the 1 A surface cells, or on a part of the surface that does not touch the rest; see docs/surface_patches.md" + ) + masks[name] = mask + + return masks + + +def calculate_patch_features(coords, areas, residues, lipo, ep, settings, hydro_scale, structure_labels): + """Returns the patch features and the per-point columns for the bundle. + + Args: + coords: (N, 3) property point coordinates on the 1 A lattice. + areas: (N,) area per point in square angstrom. + residues: (N,) residue label per point. + lipo: (N,) projected hydrophobicity per point. + ep: (N,) formal-charge electrostatic potential per point, in volts. + settings: a PatchSettings. + hydro_scale: the hydrophobicity scale lipo was projected with. + structure_labels: every residue label in the structure, to validate the sets. + + Returns: + features: whole-surface columns, then per residue set SurfArea_ + and the same columns with _ appended. + point_columns: local_mhp, local_ep and mhp_patch_id per point. + """ + + hydrophobic_threshold = settings.resolved_hydrophobic_threshold(hydro_scale) + masks = residue_set_masks(settings.residue_sets, residues, structure_labels) + if np.any(areas <= 0): + raise ValueError("every surface point needs a positive area, or its neighbourhood average is undefined") + + indices = lattice_indices(coords) + local_mhp = local_average(indices, lipo, areas, settings.radius) + local_ep = local_average(indices, ep, areas, settings.radius) + + labels = { + "SurfMhp": kept_patches(label_patches(indices, local_mhp > hydrophobic_threshold), areas, settings.min_area), + "SurfEpPos": kept_patches(label_patches(indices, local_ep > settings.ep_threshold), areas, settings.min_area), + "SurfEpNeg": kept_patches(label_patches(indices, local_ep < -settings.ep_threshold), areas, settings.min_area), + } + + # Column order is the order of src/prodes/data/features_patches.yaml, which + # documents each one; tests/test_patches.py checks the two agree. + def block(in_set): + return { + **patch_statistics("SurfMhp", labels["SurfMhp"], areas, in_set), + **hydrophobic_strength(local_mhp, areas, in_set), + **patch_statistics("SurfEpPos", labels["SurfEpPos"], areas, in_set), + **patch_statistics("SurfEpNeg", labels["SurfEpNeg"], areas, in_set), + } + + features = block(np.ones(len(coords), dtype=bool)) + for name, mask in masks.items(): + features[f"SurfArea_{name}"] = round(float(areas[mask].sum()), 3) + features.update({f"{code}_{name}": value for code, value in block(mask).items()}) + + point_columns = {"local_mhp": np.round(local_mhp, 4), "local_ep_volts": np.round(local_ep, 4), "mhp_patch_id": labels["SurfMhp"]} + + return features, point_columns diff --git a/src/prodes/core/point.py b/src/prodes/core/point.py index e15d0f3..e362b26 100644 --- a/src/prodes/core/point.py +++ b/src/prodes/core/point.py @@ -40,12 +40,17 @@ def __init__(self, x, y, z, atom): class Property_point: """Class which represents the points onto which properties are projected""" - def __init__(self, x, y, z, ep=None, lipo=None): + def __init__(self, x, y, z, ep=None, lipo=None, area=0.0, residue=""): self.x = x self.y = y self.z = z self.__ep = ep self.__lipo = lipo + # Square angstrom of solvent accessible surface this point stands for, + # and the label of the residue owning most of it. Set by + # prodes.calculations.grid_wizard.property_points_on_surface. + self.area = area + self.residue = residue @property def ep(self): diff --git a/src/prodes/data/features_patches.yaml b/src/prodes/data/features_patches.yaml new file mode 100644 index 0000000..35854e5 --- /dev/null +++ b/src/prodes/data/features_patches.yaml @@ -0,0 +1,99 @@ +# This file is the single source of truth for which patch features Prodes writes; +# prodes.feature_dictionary reads its feature codes and order from here. +# Maintained by hand: edit this file to add, drop or re-describe a patch feature. +meta: + title: Prodes surface patch features + description: The features written when the patch step is switched on, in output order. With a residue set NAME, every code is written a second time with _NAME appended, restricted to that set, preceded by SurfArea_NAME. + feature_count: 11 + written_when: --patches / patches=PatchSettings(...) + method: docs/surface_patches.md + default_settings: + radius_A: 5.0 + min_area_A2: 100.0 + hydrophobic_threshold: 0.0 + ep_threshold_volts: 0.5 + conventions: + local_value: The local value of a surface point is the area-weighted mean of the per-point value over every surface point within the radius, the point itself included. + patch: A patch is a set of surface points whose local value passes the threshold and which touch on the 1 A surface lattice. Patches smaller than the minimum area are discarded, judged on the whole patch. + area: Area is solvent accessible surface area with the 1.4 A probe, in A^2, which is larger than a molecular surface area of the same region. + residue_set: Patches are found on the whole surface. A restricted column counts only the area of patch points owned by a residue in the set. NSurf*Patch_NAME counts the kept patches with at least one point in the set; *PatchMaxArea_NAME is the largest in-set share of any one patch. + fields: + plot_name: Short label for an axis or legend. + short_explanation: One-line description of what the feature is. + long_explanation: How the feature is calculated and how to read it. + original_explanation: Always null; these features were added after Neijenhuis et al. 2024. + unit: Unit of measurement. Null when the feature is a count or dimensionless. +features: + SurfMhpPatchArea: + plot_name: Hydrophobic patch area (A^2) + short_explanation: Summed area of the hydrophobic surface patches. + long_explanation: Surface points whose local projected hydrophobicity is above the hydrophobic threshold are linked into patches where they touch; patches under the minimum area are discarded and the area of the rest is summed. Measures how much of the surface is covered by contiguous hydrophobic regions rather than by scattered hydrophobic points. + unit: Ų + original_explanation: null + SurfMhpPatchMaxArea: + plot_name: Largest hydrophobic patch (A^2) + short_explanation: Area of the largest hydrophobic surface patch. + long_explanation: The area of the single largest patch counted in SurfMhpPatchArea. 0 when no patch reaches the minimum area. + unit: Ų + original_explanation: null + NSurfMhpPatch: + plot_name: Number of hydrophobic patches + short_explanation: Number of hydrophobic surface patches at or above the minimum area. + long_explanation: The number of patches counted in SurfMhpPatchArea. + unit: null + original_explanation: null + SurfMhpLocalMax: + plot_name: Max local MHP + short_explanation: Highest local projected hydrophobicity on the surface. + long_explanation: The maximum over surface points of the projected hydrophobicity averaged over the neighbourhood radius, in the units of the hydrophobicity scale. The per-point value is dominated by the nearest one or two atoms, so the average is what gives it the length scale of a surface region; this is the most hydrophobic such region. + unit: null + original_explanation: null + SurfMhpLocalPosSum: + plot_name: Summed positive local MHP + short_explanation: Area-weighted sum of the positive part of the local projected hydrophobicity. + long_explanation: The sum over surface points of the point area times the local projected hydrophobicity where that is positive, zero elsewhere. Combines how hydrophobic and how large the hydrophobic regions are, without a threshold or a minimum patch area. The unit is square angstrom times the units of the hydrophobicity scale, so values from different scales are not comparable. On a scale with no negative values every point contributes. + unit: Ų + original_explanation: null + SurfEpPosPatchArea: + plot_name: Positive patch area (A^2) + short_explanation: Summed area of the positive electrostatic surface patches. + long_explanation: As SurfMhpPatchArea, for points whose local formal-charge electrostatic potential is above the electrostatic threshold. + unit: Ų + original_explanation: null + SurfEpPosPatchMaxArea: + plot_name: Largest positive patch (A^2) + short_explanation: Area of the largest positive electrostatic surface patch. + long_explanation: The area of the single largest patch counted in SurfEpPosPatchArea. + unit: Ų + original_explanation: null + NSurfEpPosPatch: + plot_name: Number of positive patches + short_explanation: Number of positive electrostatic surface patches at or above the minimum area. + long_explanation: The number of patches counted in SurfEpPosPatchArea. + unit: null + original_explanation: null + SurfEpNegPatchArea: + plot_name: Negative patch area (A^2) + short_explanation: Summed area of the negative electrostatic surface patches. + long_explanation: As SurfMhpPatchArea, for points whose local formal-charge electrostatic potential is below the negative of the electrostatic threshold. + unit: Ų + original_explanation: null + SurfEpNegPatchMaxArea: + plot_name: Largest negative patch (A^2) + short_explanation: Area of the largest negative electrostatic surface patch. + long_explanation: The area of the single largest patch counted in SurfEpNegPatchArea. + unit: Ų + original_explanation: null + NSurfEpNegPatch: + plot_name: Number of negative patches + short_explanation: Number of negative electrostatic surface patches at or above the minimum area. + long_explanation: The number of patches counted in SurfEpNegPatchArea. + unit: null + original_explanation: null +residue_set_only_features: + SurfArea: + plot_name: Surface area of the set (A^2) + short_explanation: Solvent accessible surface area owned by the residues of the set. + long_explanation: The summed area of the surface points whose residue is in the set. Written only as SurfArea_NAME, first among the columns of each set, so that the restricted patch areas can be normalised by the size of the region. + unit: Ų + original_explanation: null diff --git a/src/prodes/feature_dictionary.py b/src/prodes/feature_dictionary.py index 34142f7..0d09ff5 100644 --- a/src/prodes/feature_dictionary.py +++ b/src/prodes/feature_dictionary.py @@ -54,6 +54,12 @@ Why the other 51 are dropped by default is set out in ``docs/redundant_feature_analysis.md``. + +A third file, ``features_patches.yaml``, documents the surface patch features, +which are written only when the patch step is switched on (``--patches``). They +are in neither list above; ``get_patch_feature_codes`` returns them, and every +getter also accepts a patch code with a residue set suffix, such as +``SurfMhpPatchArea_CDR``. The method is in ``docs/surface_patches.md``. """ import importlib.resources @@ -65,13 +71,14 @@ REDUCED_DICTIONARY_FILE = "features_reduced.yaml" FULL_ONLY_DICTIONARY_FILE = "features_full_only.yaml" +PATCH_DICTIONARY_FILE = "features_patches.yaml" -def _read_features(filename: str) -> dict: - """Returns the ``features`` mapping of one YAML dictionary file.""" +def _read_features(filename: str, key: str = "features") -> dict: + """Returns one mapping, ``features`` by default, of one YAML dictionary file.""" with importlib.resources.open_text("prodes.data", filename, encoding="utf-8") as handle: - return yaml.safe_load(handle)["features"] + return yaml.safe_load(handle)[key] # Each file is read exactly once, at import. The two are kept as separate dicts @@ -94,6 +101,39 @@ def _read_features(filename: str) -> dict: #: need that. FULL_FEATURE_CODES = REDUCED_FEATURE_CODES + DROPPED_FEATURE_CODES +# The patch features are a separate block, written only when the patch step is +# switched on, so they are kept out of the reduced and full sets above. +_PATCH_ENTRIES = _read_features(PATCH_DICTIONARY_FILE) +_RESIDUE_SET_ONLY_ENTRIES = _read_features(PATCH_DICTIONARY_FILE, "residue_set_only_features") + +#: The whole-surface patch features, in the order Prodes writes them. With a +#: residue set NAME, each is written again as ``_NAME``. +PATCH_FEATURE_CODES = tuple(_PATCH_ENTRIES) + + +def residue_set_feature_codes(name: str) -> list[str]: + """Returns the columns Prodes writes for the residue set NAME, in output order.""" + + return [f"{code}_{name}" for code in (*_RESIDUE_SET_ONLY_ENTRIES, *PATCH_FEATURE_CODES)] + + +def _patch_entry(feature_code: str) -> dict | None: + """Returns the entry of a patch code, either whole-surface or ``_NAME``, or None. + + Only patch codes take a residue set suffix, so ``Area_CDR`` is not resolved + to ``Area``. + """ + + if feature_code in _PATCH_ENTRIES: + return _PATCH_ENTRIES[feature_code] + + code, separator, name = feature_code.rpartition("_") + # The rule PatchSettings applies to a set name: ASCII letters and digits. + if separator and name.isascii() and name.isalnum(): + return _PATCH_ENTRIES.get(code) or _RESIDUE_SET_ONLY_ENTRIES.get(code) + + return None + class FeatureDictionary: """Translates a Prodes feature code into a plot label or an explanation. @@ -117,7 +157,7 @@ def __repr__(self) -> str: return f"{type(self).__name__}({len(REDUCED_FEATURE_CODES)} reduced, {len(FULL_FEATURE_CODES)} full feature codes)" def __contains__(self, feature_code: str) -> bool: - return feature_code in self._entries + return feature_code in self._entries or _patch_entry(feature_code) is not None # --- the feature code lists ------------------------------------------------- @@ -138,6 +178,20 @@ def get_feature_codes(self, full: bool = False) -> list[str]: return list(FULL_FEATURE_CODES if full else REDUCED_FEATURE_CODES) + def get_patch_feature_codes(self, residue_sets=()) -> list[str]: + """Returns the patch feature codes a run with the patch step writes, in output order. + + Args: + residue_sets: the names of the residue sets of the run, in the order + they were given. + """ + + codes = list(PATCH_FEATURE_CODES) + for name in residue_sets: + codes.extend(residue_set_feature_codes(name)) + + return codes + def get_dropped_feature_codes(self) -> list[str]: """Returns the feature codes in the full set that the reduced set leaves out.""" @@ -211,15 +265,22 @@ def get_entry(self, feature_code: str) -> dict: (``dropped_because``, ``redundant_with``, ``drop_note``). """ - try: + if feature_code in self._entries: return dict(self._entries[feature_code]) - except KeyError as err: - raise KeyError( - f"{feature_code!r} is not a Prodes feature code. Use get_feature_codes(full=True) for the {len(FULL_FEATURE_CODES)} valid codes." - ) from err + + entry = _patch_entry(feature_code) + if entry is not None: + return dict(entry) + + raise KeyError( + f"{feature_code!r} is not a Prodes feature code. Use get_feature_codes(full=True) for the {len(FULL_FEATURE_CODES)} " + "valid codes, and get_patch_feature_codes() for the patch features." + ) def get_dictionary(self) -> dict: - """Returns every entry, keyed by feature code. + """Returns every entry of the reduced and full sets, keyed by feature code. + + The patch features are not included; get_entry looks them up. Returns: A new dict of new dicts, safe to mutate. diff --git a/src/prodes/output.py b/src/prodes/output.py index 42e6d8e..17debc7 100644 --- a/src/prodes/output.py +++ b/src/prodes/output.py @@ -75,7 +75,11 @@ {name}_features.csv the calculated features, one row {name}_surface_points.csv every surface point: x, y, z, potential, - hydrophobicity + hydrophobicity, the area it stands for, and + the residue owning most of that area; with + --patches also the neighbourhood averages + and the hydrophobic patch it belongs to + (-1 for none) {name}_ep.pml PyMOL, coloured by electrostatic potential {name}_hydrophobicity.pml PyMOL, coloured by hydrophobicity {name}_ep.cxc ChimeraX, electrostatic potential @@ -162,13 +166,17 @@ def run_metadata(pdb_file, settings, n_points, ep, disulfides=0, alternate_confo } -def bundle_files(directory, name, features, coords, ep, lipo, pdb_file, metadata, hydro_scale): - """Writes every member of the bundle into an existing directory.""" +def bundle_files(directory, name, features, coords, ep, lipo, pdb_file, metadata, hydro_scale, point_columns=None): + """Writes every member of the bundle into an existing directory. + + point_columns holds further per-point columns for the point table, in + order, such as the area and residue of each point. + """ structure_file, point_file = f"{name}.pdb", f"{name}_ep.pdb" features.to_csv(directory / f"{name}_features.csv", index=False) - pd.DataFrame({"x": coords[:, 0], "y": coords[:, 1], "z": coords[:, 2], "ep_volts": ep, "hydrophobicity": lipo}).to_csv( + pd.DataFrame({"x": coords[:, 0], "y": coords[:, 1], "z": coords[:, 2], "ep_volts": ep, "hydrophobicity": lipo, **(point_columns or {})}).to_csv( directory / f"{name}_surface_points.csv", index=False ) @@ -206,7 +214,7 @@ def bundle_files(directory, name, features, coords, ep, lipo, pdb_file, metadata (directory / "README.txt").write_text(README.format(name=name, underline="=" * (len(name) + 26))) -def read_member(bundle, suffix): +def read_member(bundle, suffix, **read_csv_options): """Reads the one bundle member whose name ends in suffix, as a dataframe.""" with zipfile.ZipFile(bundle) as archive: @@ -216,7 +224,7 @@ def read_member(bundle, suffix): raise ValueError(f"{bundle} holds no member ending in {suffix}") from None with archive.open(name) as handle: - return pd.read_csv(handle) + return pd.read_csv(handle, **read_csv_options) def read_features(bundle): @@ -230,12 +238,17 @@ def read_features(bundle): def read_surface_points(bundle): - """Reads the surface point table out of a bundle: coordinates, potential, hydrophobicity.""" + """Reads the surface point table out of a bundle: coordinates, potential, hydrophobicity, area, residue. + + The residue column is read as text, so that a label such as 100 from a + structure with no chain identifier stays the string "100", and a blank one + stays "" rather than becoming NaN. + """ - return read_member(bundle, "_surface_points.csv") + return read_member(bundle, "_surface_points.csv", converters={"residue": str}) -def write_bundle(out_file, name, features, coords, ep, lipo, pdb_file, metadata, hydro_scale): +def write_bundle(out_file, name, features, coords, ep, lipo, pdb_file, metadata, hydro_scale, point_columns=None): """Writes the bundle to out_file, which must be a .zip path. One structure in, one archive out. Supporting a directory as well would mean @@ -254,7 +267,7 @@ def write_bundle(out_file, name, features, coords, ep, lipo, pdb_file, metadata, with TemporaryDirectory() as scratch: staged = Path(scratch) / name staged.mkdir() - bundle_files(staged, name, features, coords, ep, lipo, pdb_file, metadata, hydro_scale) + bundle_files(staged, name, features, coords, ep, lipo, pdb_file, metadata, hydro_scale, point_columns=point_columns) with zipfile.ZipFile(out_file, "w", zipfile.ZIP_DEFLATED) as archive: for member in sorted(staged.iterdir()): diff --git a/src/prodes/run.py b/src/prodes/run.py index 46c2347..fc455d3 100644 --- a/src/prodes/run.py +++ b/src/prodes/run.py @@ -8,6 +8,14 @@ from prodes.calculations import geometry, grid_wizard from prodes.calculations.distance_functions import debye_length +from prodes.calculations.patches import ( + DEFAULT_EP_THRESHOLD, + DEFAULT_PATCH_MIN_AREA, + DEFAULT_PATCH_RADIUS, + PatchSettings, + calculate_patch_features, + check_residue_sets, +) from prodes.calculations.sasa import shape, shrake_rupley from prodes.calculations.standard_equations import trimean from prodes.io import parser as ps @@ -60,6 +68,33 @@ def parse_arguments(): default=DEFAULT_IONIC_STRENGTH_MOLAR, ) parser.add_argument("--hydro", help="Abbreviation of the hydrophobicity scale to be used", type=str, default="mj_scaled") + parser.add_argument( + "--patches", + action="store_true", + help="Also report contiguous hydrophobic, positive and negative surface patches; see docs/surface_patches.md", + ) + parser.add_argument( + "--residue-set", + action="append", + default=[], + metavar="NAME=FILE", + help="With --patches, also report every patch feature restricted to the residues listed in FILE, one label " + "such as H100A per line, under columns ending in _NAME. Repeatable", + ) + parser.add_argument("--patch-radius", type=float, default=None, help=f"Neighbourhood radius in A for the patch average (default: {DEFAULT_PATCH_RADIUS})") + parser.add_argument("--patch-min-area", type=float, default=None, help=f"Smallest patch kept, in A^2 (default: {DEFAULT_PATCH_MIN_AREA})") + parser.add_argument( + "--patch-hydrophobic-threshold", + type=float, + default=None, + help="Neighbourhood-averaged hydrophobicity above which a point is hydrophobic (default: 0, refused for a scale with no negative values)", + ) + parser.add_argument( + "--patch-ep-threshold", + type=float, + default=None, + help=f"Neighbourhood-averaged potential in volts beyond which a point is positive or negative (default: {DEFAULT_EP_THRESHOLD})", + ) parser.add_argument( "--full-features", action=argparse.BooleanOptionalAction, @@ -113,7 +148,50 @@ def parse_arguments(): if not math.isfinite(ionic_strength_molar) or ionic_strength_molar < 0: raise ValueError(f"--ionic-strength must be a finite number and not negative, got {ionic_strength_molar}") - return pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar + patch_options = { + "--patch-radius": arg.patch_radius, + "--patch-min-area": arg.patch_min_area, + "--patch-hydrophobic-threshold": arg.patch_hydrophobic_threshold, + "--patch-ep-threshold": arg.patch_ep_threshold, + } + patches = None + if arg.patches: + patches = PatchSettings( + radius=DEFAULT_PATCH_RADIUS if arg.patch_radius is None else arg.patch_radius, + min_area=DEFAULT_PATCH_MIN_AREA if arg.patch_min_area is None else arg.patch_min_area, + hydrophobic_threshold=arg.patch_hydrophobic_threshold, + ep_threshold=DEFAULT_EP_THRESHOLD if arg.patch_ep_threshold is None else arg.patch_ep_threshold, + residue_sets=read_residue_sets(arg.residue_set), + ) + else: + # Without --patches these would be silently ignored. + given = [option for option, value in patch_options.items() if value is not None] + (["--residue-set"] if arg.residue_set else []) + if given: + raise ValueError(f"{', '.join(given)} only apply to the patch features, so they need --patches") + + return pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar, patches + + +def read_residue_sets(specifications): + """Reads the --residue-set NAME=FILE options into a name to labels mapping. + + Each file holds one residue label per line, such as H100A. Blank lines and + anything after a # are ignored. + """ + + residue_sets = {} + for specification in specifications: + name, separator, path = specification.partition("=") + if not separator or not name or not path: + raise ValueError(f"--residue-set takes NAME=FILE, got {specification!r}") + if name in residue_sets: + raise ValueError(f"--residue-set {name} is given twice") + + with open(path, encoding="utf-8") as handle: + labels = [line.split("#", 1)[0].strip() for line in handle] + residue_sets[name] = [label for label in labels if label] + + return residue_sets def standard_features(values, name="", reduced=False): @@ -475,7 +553,7 @@ def construct_surface_grid(structure, r_probe): surface_grid.construct_cells(surface) surface_grid.fill_cells(surface) - property_points = grid_wizard.property_points_on_surface(surface_grid) + property_points = grid_wizard.property_points_on_surface(surface_grid, probe_r=r_probe) return property_points @@ -490,6 +568,7 @@ def calculate( full_features=False, mem_limit_mb=None, ionic_strength_molar=DEFAULT_IONIC_STRENGTH_MOLAR, + patches=None, ): """Calculates a list of supported features and returns a csv file @@ -512,6 +591,10 @@ def calculate( run, divided among the workers. If None, falls back to the PRODES_MEM_LIMIT_MB env var (default 2048). Pass explicitly when calling as a library to avoid relying on environment variables. + patches: a prodes.calculations.patches.PatchSettings to also report + contiguous surface patches, optionally restricted to named residue + sets; None, the default, leaves them out and the feature table is + the same as without this option. Thread safety Safe to run in several processes, one structure each, which is how the @@ -524,6 +607,10 @@ def calculate( validate_settings(mem_limit_mb) structure = prepare_structure(pdb_file, pkas_file) + # Both checks are cheap, and the surface phase that follows is not. + if patches is not None: + patches.resolved_hydrophobic_threshold(hydro_scale) + 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. print(f"{structure.name}: {len(structure.disulfides)} disulfide bonds detected") @@ -540,10 +627,31 @@ def calculate( # formal charges, so a snapshot taken afterwards would hold a different # physical quantity under the same name, and only when --full-features is set. coords, ep, lipo = surface_point_arrays(surface_points) + areas = np.array([point.area for point in surface_points]) + residues = np.array([point.residue for point in surface_points]) + point_columns = {"area_A2": areas, "residue": residues} + + # Computed from the snapshot above, so the average charge phase overwriting + # the potential on the points cannot reach it, and appended after the shell + # features so that the default columns keep their positions. + patch_features: dict = {} + if patches is not None: + patch_features, patch_columns = calculate_patch_features( + coords, + areas, + residues, + lipo, + ep, + patches, + hydro_scale, + [residue.label for residue in structure.residues], + ) + point_columns.update(patch_columns) if full_features: features = calculate_average_chargesurface_grid_features(structure, surface_points, ph, features, ionic_strength_molar=ionic_strength_molar) features = calculate_shell_features(structure, surface_points, ph, features, full_features=full_features, mem_limit_mb=mem_limit_mb) + features.update(patch_features) calculated_features = pd.Series(features).to_frame().transpose() settings = { @@ -553,6 +661,8 @@ def calculate( "full_features": full_features, "ionic_strength_molar": ionic_strength_molar, } + if patches is not None: + settings["patches"] = patches.record(hydro_scale) metadata = run_metadata( pdb_file, settings, @@ -564,11 +674,13 @@ def calculate( inferred_elements=structure.inferred_elements, ) - return write_bundle(out_file, structure.name, calculated_features, coords, ep, lipo, pdb_file, metadata, hydro_scale) + return write_bundle(out_file, structure.name, calculated_features, coords, ep, lipo, pdb_file, metadata, hydro_scale, point_columns=point_columns) def main(): parsed_arguments = parse_arguments() + # Positional, so a new option of calculate has to be appended to both + # signatures in the same place; tests/test_cli.py pins the order. calculate(*parsed_arguments) diff --git a/tests/test_cli.py b/tests/test_cli.py index 33acedb..7ab1d17 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -11,6 +11,7 @@ import pytest +from prodes.calculations.patches import PatchSettings from prodes.run import DEFAULT_IONIC_STRENGTH_MOLAR, calculate, main, parse_arguments @@ -26,13 +27,14 @@ 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 = parse(monkeypatch) + pdb_file, out_file, pkas_file, ph, r_probe, hydro_scale, full_features, mem_limit_mb, ionic_strength_molar, patches = 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") assert full_features is False assert mem_limit_mb is None assert ionic_strength_molar == DEFAULT_IONIC_STRENGTH_MOLAR + assert patches is None def test_n_workers_flag_sets_the_environment_variable(monkeypatch): @@ -56,7 +58,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 = parse(monkeypatch, "--mem-limit", "256") + *_, mem_limit_mb, _ionic_strength, _patches = parse(monkeypatch, "--mem-limit", "256") assert mem_limit_mb == 256 @@ -108,6 +110,7 @@ def fake_calculate(*args, **kwargs): True, None, DEFAULT_IONIC_STRENGTH_MOLAR, + None, ) assert captured["kwargs"] == {} @@ -125,4 +128,62 @@ def test_calculate_signature_matches_the_parsed_order(): "full_features", "mem_limit_mb", "ionic_strength_molar", + "patches", ] + + +def test_patches_flag_builds_the_settings(monkeypatch, tmp_path): + """--patches and its options reach calculate as one PatchSettings, residue sets read from file.""" + + labels = tmp_path / "cdr.txt" + labels.write_text("# CDR H3\nH95\nH100A # an insertion\n\nL91\n") + *_, patches = parse( + monkeypatch, + "--patches", + "--patch-radius", + "4", + "--patch-min-area", + "30", + "--patch-hydrophobic-threshold", + "0.02", + "--patch-ep-threshold", + "0.3", + "--residue-set", + f"CDR={labels}", + ) + + assert patches == PatchSettings(radius=4, min_area=30, hydrophobic_threshold=0.02, ep_threshold=0.3, residue_sets={"CDR": ["H95", "H100A", "L91"]}) + + +def test_residue_set_without_patches_is_refused(monkeypatch, tmp_path): + """A residue set restricts the patch features, so on its own it would silently do nothing.""" + + labels = tmp_path / "cdr.txt" + labels.write_text("H95\n") + + with pytest.raises(ValueError, match="need --patches"): + parse(monkeypatch, "--residue-set", f"CDR={labels}") + + +@pytest.mark.parametrize("specification", ["CDR", "=cdr.txt", "CDR="]) +def test_a_malformed_residue_set_is_refused(monkeypatch, specification): + """NAME=FILE is the only accepted form.""" + + with pytest.raises(ValueError, match="NAME=FILE"): + parse(monkeypatch, "--patches", "--residue-set", specification) + + +@pytest.mark.parametrize( + "option", [["--patch-radius", "3"], ["--patch-min-area", "80"], ["--patch-hydrophobic-threshold", "0.1"], ["--patch-ep-threshold", "0.2"]] +) +def test_a_patch_option_without_patches_is_refused(monkeypatch, option): + """Otherwise the option would be accepted and silently do nothing.""" + + with pytest.raises(ValueError, match=f"{option[0]} only apply to the patch features"): + parse(monkeypatch, *option) + + +def test_patches_alone_gives_the_documented_defaults(monkeypatch): + *_, patches = parse(monkeypatch, "--patches") + + assert patches == PatchSettings() diff --git a/tests/test_output_bundle.py b/tests/test_output_bundle.py index 9c29bcb..5d763c1 100644 --- a/tests/test_output_bundle.py +++ b/tests/test_output_bundle.py @@ -80,7 +80,7 @@ def test_the_tables_read_back_out_of_the_archive(bundle_zip): assert len(features) == 1 assert features["ID"].iloc[0] == "ARH96693" - assert list(points.columns) == ["x", "y", "z", "ep_volts", "hydrophobicity"] + assert list(points.columns) == ["x", "y", "z", "ep_volts", "hydrophobicity", "area_A2", "residue"] assert len(points) > 0 diff --git a/tests/test_patches.py b/tests/test_patches.py new file mode 100644 index 0000000..0023c82 --- /dev/null +++ b/tests/test_patches.py @@ -0,0 +1,424 @@ +"""Tests for the surface patch step, prodes.calculations.patches. + +The synthetic tests build points directly on the 1 A lattice, where the right +answer can be counted by hand. The tests on a real structure check what the +synthetic ones cannot: that the step is wired into calculate without changing +any default feature, that the point areas add up to the surface area, and that +a residue set selects the points it names. +""" + +import tempfile +from functools import cache +from pathlib import Path +from types import SimpleNamespace + +import numpy as np +import pytest + +from prodes import run as run_module +from prodes.calculations import grid_wizard +from prodes.calculations.patches import ( + DEFAULT_HYDROPHOBIC_THRESHOLD, + PatchSettings, + calculate_patch_features, + kept_patches, + label_patches, + lattice_indices, + local_average, + neighbourhood_offsets, + patch_statistics, + residue_set_masks, +) +from prodes.feature_dictionary import REDUCED_FEATURE_CODES, FeatureDictionary +from prodes.output import read_features, read_surface_points, write_bundle +from prodes.run import calculate, construct_surface_grid, prepare_structure + +STRUCTURE = "tests/data/ARH96693.pdb.zip" + +# The first twenty residues of ARH96693, a single chain named A. +N_TERMINUS = [f"A{number}" for number in range(1, 21)] + + +def block(origin, size): + """Returns the lattice points of a cube of side size starting at origin.""" + + steps = range(size) + return [(origin[0] + x, origin[1] + y, origin[2] + z) for x in steps for y in steps for z in steps] + + +# --- the lattice -------------------------------------------------------------- + + +def test_lattice_indices_recover_whole_number_offsets(): + """Cell centres carry float rounding; the indices must not.""" + + coords = np.array([[10.1, 5.3, -2.7], [11.1, 5.3, -2.7], [10.1 + 1e-9, 7.3, -1.7]]) + + assert lattice_indices(coords).tolist() == [[0, 0, 0], [1, 0, 0], [0, 2, 1]] + + +def test_points_off_the_lattice_are_refused(): + with pytest.raises(ValueError, match="1 A lattice"): + lattice_indices(np.array([[0.0, 0.0, 0.0], [0.5, 0.0, 0.0]])) + + +def test_points_sharing_a_position_are_refused(): + with pytest.raises(ValueError, match="share a lattice position"): + lattice_indices(np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])) + + +def test_the_neighbourhood_boundary_is_exact(): + """(3, 4, 0) lies exactly 5 A away and is inside; (3, 4, 1) is not.""" + + offsets = {tuple(offset) for offset in neighbourhood_offsets(5.0)} + + assert (3, 4, 0) in offsets + assert (0, 0, 5) in offsets + assert (3, 4, 1) not in offsets + assert (0, 0, 0) in offsets + + +# --- the neighbourhood average ------------------------------------------------ + + +def test_local_average_on_a_line(): + """One hydrophobic point at the end of a line of eleven, radius 2.""" + + coords = np.array([[x, 0.0, 0.0] for x in range(11)]) + values = np.zeros(11) + values[0] = 1.0 + + local = local_average(lattice_indices(coords), values, np.ones(11), 2.0) + + # Point 0 sees points 0, 1, 2; point 2 sees 0 to 4; point 3 no longer sees 0. + assert local[0] == pytest.approx(1 / 3) + assert local[2] == pytest.approx(1 / 5) + assert local[3] == 0 + + +def test_local_average_is_weighted_by_area(): + coords = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]]) + + local = local_average(lattice_indices(coords), np.array([1.0, 0.0]), np.array([3.0, 1.0]), 1.0) + + assert local == pytest.approx([0.75, 0.75]) + + +# --- patches ------------------------------------------------------------------ + + +def synthetic_surface(): + """Three separated cubes: 27 points, 8 points and 1 point, every point 1 A^2.""" + + points = block((0, 0, 0), 3) + block((10, 0, 0), 2) + [(20, 0, 0)] + coords = np.array(points, dtype=float) + return lattice_indices(coords), np.ones(len(points)) + + +def test_separated_regions_are_separate_patches(): + indices, _ = synthetic_surface() + + labels = label_patches(indices, np.ones(len(indices), dtype=bool)) + + assert sorted(np.bincount(labels).tolist()) == [1, 8, 27] + + +def test_diagonal_neighbours_are_linked(): + """The surface is built on the 26-neighbourhood, so a corner contact joins two points.""" + + indices = lattice_indices(np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0], [3.0, 3.0, 3.0]])) + + assert label_patches(indices, np.ones(3, dtype=bool)).tolist() == [0, 0, 1] + + +def test_points_outside_the_mask_are_in_no_patch(): + indices, _ = synthetic_surface() + mask = np.zeros(len(indices), dtype=bool) + mask[:27] = True + + labels = label_patches(indices, mask) + + assert set(labels[27:]) == {-1} + assert set(labels[:27]) == {0} + + +def test_patches_under_the_minimum_area_are_discarded(): + indices, areas = synthetic_surface() + labels = kept_patches(label_patches(indices, np.ones(len(indices), dtype=bool)), areas, 8.0) + + statistics = patch_statistics("SurfMhp", labels, areas, np.ones(len(indices), dtype=bool)) + + assert statistics == {"SurfMhpPatchArea": 35.0, "SurfMhpPatchMaxArea": 27.0, "NSurfMhpPatch": 2} + + +def test_restriction_counts_the_in_set_share_of_whole_patches(): + """A set holding 9 points of the 27-point patch and nothing of the 1-point one.""" + + indices, areas = synthetic_surface() + labels = kept_patches(label_patches(indices, np.ones(len(indices), dtype=bool)), areas, 8.0) + in_set = np.zeros(len(indices), dtype=bool) + in_set[:9] = True + in_set[-1] = True + + statistics = patch_statistics("SurfMhp", labels, areas, in_set) + + # The minimum area was judged on the whole patch, so the 9-point share + # counts although it is itself above 8, and the 1-point patch, discarded + # on the whole surface, does not come back. + assert statistics == {"SurfMhpPatchArea": 9.0, "SurfMhpPatchMaxArea": 9.0, "NSurfMhpPatch": 1} + + +def test_no_patch_gives_zeros(): + indices, areas = synthetic_surface() + labels = np.full(len(indices), -1) + + assert patch_statistics("SurfEpPos", labels, areas, np.ones(len(indices), dtype=bool)) == { + "SurfEpPosPatchArea": 0.0, + "SurfEpPosPatchMaxArea": 0.0, + "NSurfEpPosPatch": 0, + } + + +# --- settings ----------------------------------------------------------------- + + +def test_the_default_threshold_is_refused_on_a_scale_with_no_negative_values(): + with pytest.raises(ValueError, match="no negative values"): + PatchSettings().resolved_hydrophobic_threshold("mj") + + +def test_an_explicit_threshold_is_accepted_on_any_scale(): + assert PatchSettings(hydrophobic_threshold=6.0).resolved_hydrophobic_threshold("mj") == 6.0 + assert PatchSettings().resolved_hydrophobic_threshold("mj_scaled") == DEFAULT_HYDROPHOBIC_THRESHOLD + + +@pytest.mark.parametrize( + ("arguments", "message"), + [ + ({"radius": -1.0}, "radius"), + ({"radius": 0.5}, "lattice spacing"), + ({"radius": 50.0}, "between"), + ({"radius": "5"}, "radius"), + ({"hydrophobic_threshold": "0.1"}, "hydrophobic_threshold"), + ({"residue_sets": {"CDR": "H95"}}, "list of residue labels"), + ({"residue_sets": None}, "dict"), + ({"min_area": float("nan")}, "min_area"), + ({"ep_threshold": -0.1}, "ep_threshold"), + ({"residue_sets": {"CDR-H3": ["H95"]}}, "letters and digits"), + ({"residue_sets": {"CDR": []}}, "empty"), + ], +) +def test_bad_settings_are_refused(arguments, message): + with pytest.raises(ValueError, match=message): + PatchSettings(**arguments) + + +def test_a_residue_the_structure_does_not_have_is_refused(): + with pytest.raises(ValueError, match="H999"): + residue_set_masks({"CDR": ["H95", "H999"]}, np.array(["H95"]), ["H95", "H96"]) + + +def test_a_residue_set_with_no_surface_points_is_refused(): + with pytest.raises(ValueError, match="owns a surface point"): + residue_set_masks({"Core": ["H96"]}, np.array(["H95"]), ["H95", "H96"]) + + +# --- the residue of a surface cell -------------------------------------------- + + +def dot(residue_label, radius): + """A Shrake-Rupley dot on an atom of the given radius, in the given residue.""" + + atom = SimpleNamespace( + residue=SimpleNamespace(label=residue_label), + max_area=lambda probe_r: 4 * np.pi * (radius + probe_r) ** 2, + max_n_points=lambda probe_r: int(4 * np.pi * (radius + probe_r) ** 2 * 2), + ) + return SimpleNamespace(atom=atom) + + +def cell_of(*dots): + return SimpleNamespace(filtered_content=lambda *_: list(dots)) + + +def test_a_cell_belongs_to_the_residue_with_most_area(): + """Two dots of A5 outweigh one of A4.""" + + area, residue = grid_wizard.surface_cell_area_and_residue(cell_of(dot("A4", 1.7), dot("A5", 1.7), dot("A5", 1.7)), 1.4, {}) + + assert residue == "A5" + assert area == pytest.approx(3 * 4 * np.pi * 3.1**2 / int(4 * np.pi * 3.1**2 * 2)) + + +def test_a_tie_goes_to_the_first_label_whatever_the_order(): + first = grid_wizard.surface_cell_area_and_residue(cell_of(dot("A5", 1.7), dot("A4", 1.7)), 1.4, {}) + second = grid_wizard.surface_cell_area_and_residue(cell_of(dot("A4", 1.7), dot("A5", 1.7)), 1.4, {}) + + assert first[1] == second[1] == "A4" + + +# --- on a real structure ------------------------------------------------------ + + +@cache +def _run(full_features, patched): + with tempfile.TemporaryDirectory(prefix="prodes_patch_tests_") as directory: + out_file = Path(directory) / "bundle.zip" + patches = PatchSettings(residue_sets={"Nterm": N_TERMINUS}) if patched else None + calculate(STRUCTURE, str(out_file), full_features=full_features, patches=patches) + return read_features(out_file), read_surface_points(out_file) + + +def run(full_features=False, patched=False): + """Returns copies of the feature and point tables of one run per setting. + + With patched, the patch step runs with one residue set, Nterm. + """ + + features, points = _run(full_features, patched) + return features.copy(), points.copy() + + +def test_point_areas_add_up_to_the_surface_area(): + """ARH96693 has no internal cavity, so the flood fill reaches every dot.""" + + structure = prepare_structure(STRUCTURE, None) + points = construct_surface_grid(structure, 1.4) + + assert sum(point.area for point in points) == pytest.approx(structure.surface_area(1.4), rel=1e-9) + assert all(point.residue for point in points) + + +def test_the_default_features_are_unchanged_by_the_patch_step(): + """The anion exchange mode depends on these values; the patch step only appends.""" + + default, _ = run() + patched, _ = run(patched=True) + + assert list(default.columns) == ["ID", *REDUCED_FEATURE_CODES] + assert list(patched.columns[: len(default.columns)]) == list(default.columns) + assert patched[default.columns].equals(default) + + +def test_the_columns_are_the_ones_the_dictionary_documents(): + patched, _ = run(patched=True) + + assert list(patched.columns[1 + len(REDUCED_FEATURE_CODES) :]) == FeatureDictionary().get_patch_feature_codes(["Nterm"]) + + +def test_the_patch_step_reads_the_formal_potential_under_full_features(): + """The partial-charge phase overwrites the potential on the points; the patch step must run before it.""" + + reduced, _ = run(patched=True) + full, _ = run(full_features=True, patched=True) + codes = FeatureDictionary().get_patch_feature_codes(["Nterm"]) + + assert full[codes].equals(reduced[codes]) + + +def test_the_residue_set_selects_the_points_it_names(): + features, points = run(patched=True) + in_set = points["residue"].isin(N_TERMINUS) + + assert features["SurfArea_Nterm"].iloc[0] == pytest.approx(points.loc[in_set, "area_A2"].sum(), abs=1e-3) + assert set(points.loc[in_set, "residue"]) <= set(N_TERMINUS) + assert 0 < in_set.sum() < len(points) + + +def test_the_point_table_agrees_with_the_hydrophobic_patch_area(): + features, points = run(patched=True) + in_patch = points["mhp_patch_id"] >= 0 + + assert features["SurfMhpPatchArea"].iloc[0] == pytest.approx(points.loc[in_patch, "area_A2"].sum(), abs=1e-3) + assert features["NSurfMhpPatch"].iloc[0] == points.loc[in_patch, "mhp_patch_id"].nunique() + + +def test_restricted_values_do_not_exceed_whole_surface_ones(): + features = run(patched=True)[0].iloc[0] + + for code in ("SurfMhpPatchArea", "SurfEpPosPatchArea", "SurfEpNegPatchArea", "SurfMhpLocalPosSum"): + assert features[f"{code}_Nterm"] <= features[code] + assert features["SurfMhpLocalMax_Nterm"] <= features["SurfMhpLocalMax"] + + +def test_calculate_patch_features_refuses_an_unknown_residue_before_any_work(): + with pytest.raises(ValueError, match="H999"): + calculate_patch_features( + np.zeros((1, 3)), + np.ones(1), + np.array(["A1"]), + np.zeros(1), + np.zeros(1), + PatchSettings(residue_sets={"CDR": ["H999"]}), + "mj_scaled", + ["A1"], + ) + + +def test_settings_with_residue_sets_can_be_hashed(): + """Frozen settings are hashable even though the residue sets are a dict.""" + + assert hash(PatchSettings(residue_sets={"CDR": ["H95"]})) == hash(PatchSettings(residue_sets={"CDR": ["H96"]})) + + +def test_a_point_with_no_area_is_refused(): + with pytest.raises(ValueError, match="positive area"): + calculate_patch_features(np.zeros((1, 3)), np.zeros(1), np.array(["A1"]), np.zeros(1), np.zeros(1), PatchSettings(), "mj_scaled", ["A1"]) + + +def test_calculate_refuses_bad_patch_settings_before_building_the_surface(monkeypatch, tmp_path): + """An unknown residue or a refused threshold fails in the first second, not after the surface phase.""" + + def surface_was_built(*args, **kwargs): + raise AssertionError("the surface was built before the patch settings were checked") + + monkeypatch.setattr(run_module, "construct_surface_grid", surface_was_built) + + with pytest.raises(ValueError, match="A999"): + calculate(STRUCTURE, str(tmp_path / "out.zip"), patches=PatchSettings(residue_sets={"X": ["A999"]})) + with pytest.raises(ValueError, match="no negative values"): + calculate(STRUCTURE, str(tmp_path / "out.zip"), hydro_scale="mj", patches=PatchSettings()) + + +def test_residue_labels_read_back_as_text(tmp_path): + """A chainless label such as 100 must not come back as the integer 100, nor a blank one as NaN.""" + + import pandas as pd + + bundle = write_bundle( + tmp_path / "out.zip", + "x", + pd.DataFrame([{"ID": "x", "Area": 1.0}]), + np.zeros((2, 3)), + np.zeros(2), + np.zeros(2), + STRUCTURE, + {"prodes_version": "test"}, + "mj_scaled", + point_columns={"area_A2": np.ones(2), "residue": np.array(["100", ""])}, + ) + + assert read_surface_points(bundle)["residue"].tolist() == ["100", ""] + + +# --- the feature dictionary --------------------------------------------------- + + +def test_the_dictionary_resolves_a_residue_set_suffix_on_patch_codes_only(): + fd = FeatureDictionary() + + assert fd.get_entry("SurfMhpPatchArea_CDR") == fd.get_entry("SurfMhpPatchArea") + assert "SurfArea_CDR" in fd + assert "SurfArea" not in fd + assert "Area_CDR" not in fd + assert "SurfMhpPatchArea_CDR_H3" not in fd + assert "SurfMhpPatchArea_CDR²" not in fd + with pytest.raises(KeyError, match="get_patch_feature_codes"): + fd.get_entry("Area_CDR") + + +def test_the_patch_codes_are_kept_out_of_the_reduced_and_full_sets(): + fd = FeatureDictionary() + + assert set(fd.get_patch_feature_codes()).isdisjoint(fd.get_feature_codes(full=True)) + assert fd.get_patch_feature_codes(["CDR"])[11] == "SurfArea_CDR" diff --git a/tests/test_run_prodes.py b/tests/test_run_prodes.py index 3e86ed5..8f5bf22 100644 --- a/tests/test_run_prodes.py +++ b/tests/test_run_prodes.py @@ -11,6 +11,7 @@ import prodes from prodes import run as run_module +from prodes.calculations.patches import PatchSettings from prodes.feature_dictionary import FULL_FEATURE_CODES, REDUCED_FEATURE_CODES from prodes.output import read_features @@ -59,6 +60,7 @@ def fake_calculate(pdb_file, out_file, **kwargs): hydro_scale="kd", full_features=True, mem_limit_mb=512, + patches=PatchSettings(radius=4.0), ) assert captured == { @@ -71,6 +73,7 @@ def fake_calculate(pdb_file, out_file, **kwargs): "full_features": True, "mem_limit_mb": 512, "ionic_strength_molar": 0.15, + "patches": PatchSettings(radius=4.0), } From d91d3aaca8e401eb77f9e6628866aee957cdc5c5 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Fri, 16 Oct 2026 21:38:38 +0200 Subject: [PATCH 2/4] #35 Print the Python stack on a crash, and document missing pKa groups and silent exits A Prodes run that ends with a non-zero exit code, no traceback and no bundle gives its caller nothing to act on. The command line now enables faulthandler, so a crash in compiled code prints the Python stack of every thread to stderr before the process exits; forked workers inherit it. A SIGKILL, typically from the out-of-memory killer, cannot be caught, so docs/parallelism_and_memory.md now gives the return codes of both cases under subprocess and where to look for an out-of-memory kill. docs/structure_preparation.md now states that hydrogens do not change any feature value (stripping all hydrogens from an ABodyBuilder3 Fv model left all 77 columns of a --patches run unchanged), and names the two causes of the warning about titratable groups missing from the pKa file. Both are gaps in PROPKA's own output rather than values dropped by convert_propka: a C-terminus without OXT gets no C- group, and PROPKA 3.5.1 gives no side-chain pKa for an aspartate that is the first residue of a chain, as on most antibody light chains. --- docs/parallelism_and_memory.md | 11 +++++++++++ docs/structure_preparation.md | 7 +++++++ src/prodes/run.py | 6 ++++++ 3 files changed, 24 insertions(+) diff --git a/docs/parallelism_and_memory.md b/docs/parallelism_and_memory.md index c064d60..7e5f729 100644 --- a/docs/parallelism_and_memory.md +++ b/docs/parallelism_and_memory.md @@ -111,6 +111,17 @@ The budget covers the chunked arrays only. The parsed structure, the grids and t When running several proteins side by side, one process each, RAM is **cumulative** across those processes: each one gets its own budget. Divide the machine's memory by the number of concurrent runs, and leave headroom for the operating system. +### A run that stops with no traceback + +A Python error in a worker is reported, and the phase is repeated serially with a warning naming it. A run that ends with a non-zero exit code, no traceback and no bundle was instead stopped by a signal. Under `subprocess`, the return code is the negative signal number: + +| Return code | Signal | Usual cause | +| --- | --- | --- | +| -9 (137 from a shell) | SIGKILL | the kernel's out-of-memory killer, or a scheduler enforcing a limit | +| -11 (139 from a shell) | SIGSEGV | a crash in compiled code | + +The command line enables Python's `faulthandler`, so a crash in compiled code prints the Python stack to stderr before the process exits. A SIGKILL cannot be caught by any process, so it leaves nothing behind; record the return code, and check `dmesg` or `journalctl -k` for `Out of memory` at that time. Several Prodes runs side by side each have their own memory budget, and the structure, grids and surface points come on top of it. + ## BLAS threads Workers force single-threaded BLAS on themselves (`OPENBLAS_NUM_THREADS=1` and friends, applied via `threadpoolctl`), and the parent is held to one thread for as long as the workers exist. This is deliberate: parallelism should happen across processes rather than within them, and letting BLAS spawn its own threads inside each worker would oversubscribe the machine badly. diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index 9aac784..0ce7494 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -47,6 +47,8 @@ Only the **geometry**. Prodes reads `ATOM` records, filters out hydrogens, and r So the hydrogens PDB2PQR adds, the protonation states it assigns and the partial charges it writes are all discarded. That is deliberate: it is what lets you run Prodes at several pH values against one prepared structure. What you gain from the preparation step is that the atoms are there and in sensible positions. +A structure without hydrogens is therefore a valid input. Stripping all 1,625 hydrogens from an ABodyBuilder3 Fv model, with the same pKa file, left every one of the 77 feature columns of a `--patches` run unchanged. What does change the values is a missing heavy atom, and on a predicted structure that is almost always `OXT`; see the next section. + ## What happens to a predicted structure: the terminal oxygen Run PDB2PQR over a set of predicted structures and you will usually see it add exactly one atom to each, and nothing else. That atom is `OXT`, and it is worth understanding because of what Prodes currently does with it. @@ -129,6 +131,11 @@ What you should take from PDB2PQR is the **ordering**. It runs PROPKA after repa The four extra groups are the two rebuilt lysine side chains and the two C-termini that only exist once `OXT` has been added. Among the 145 groups both runs share, six shift by more than 0.5 pKa units. +A group PROPKA does not report keeps the textbook pKa of its residue type, and Prodes names every such group in one warning per run, such as `6 titratable groups are not in the pKa file and keep the textbook value for their residue type: C- A451, ...`. Two causes account for most of them, and neither is `convert_propka` dropping a value: the group is absent from PROPKA's own output. + +* **A C-terminus without `OXT`.** PROPKA reports no `C-` group for a chain whose last residue has no `OXT`. Running PDB2PQR first adds it. +* **An aspartate as the first residue of a chain.** PROPKA 3.5.1 reports the N-terminal `N+` group of such a chain but no side-chain pKa for the aspartate, while it does report one for an N-terminal glutamate. Seen on the light chains of whole antibodies, which commonly start with Asp. Structure preparation does not change this. + ## Chain identifiers Prodes needs chain identifiers. It groups residues by chain, decides which cysteines are bonded into a disulfide from the chain and residue number together, and from version 9.0 looks up a predicted pKa by chain as well as residue number. Flattening the chains would change all three. diff --git a/src/prodes/run.py b/src/prodes/run.py index fc455d3..326c9a1 100644 --- a/src/prodes/run.py +++ b/src/prodes/run.py @@ -1,4 +1,5 @@ import argparse +import faulthandler import math import os @@ -678,6 +679,11 @@ def calculate( def main(): + # A crash in compiled code would otherwise end the run with no message at + # all; this prints the Python stack of every thread to stderr first. Forked + # workers inherit it. A SIGKILL, such as from the out-of-memory killer, + # cannot be caught by anything. + faulthandler.enable() parsed_arguments = parse_arguments() # Positional, so a new option of calculate has to be appended to both # signatures in the same place; tests/test_cli.py pins the order. From 5267011cd4161826561b71a2095071d161e613c1 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sun, 18 Oct 2026 20:25:06 +0200 Subject: [PATCH 3/4] #35 Record that the patch defaults were checked against antibody HIC data and kept The patch settings were tuned against HIC retention time on the GDPa1 antibody panel (236 antibodies) over radius 3 to 12 A, minimum area 0 to 500 A^2, hydrophobic threshold -0.08 to 0.06 and electrostatic threshold 0.1 to 1.0 V. The best setting on GDPa1 raised the Spearman R of SurfMhpPatchArea_CDR there from 0.34 to 0.46, but lowered it on the held-out GDPa3 panel (79 antibodies) from 0.355 to 0.328, and an HIC model fitted on GDPa1 scored 0.589 on GDPa3 with the tuned setting against 0.615 with the defaults. The defaults therefore stay as they are. docs/surface_patches.md, the README and the comment on DEFAULT_PATCH_MIN_AREA now say so, in place of calling the defaults uncalibrated starting values, and state that any re-tuning needs an independent panel to check it on. --- README.rst | 2 +- docs/surface_patches.md | 10 +++++++--- src/prodes/calculations/patches.py | 2 +- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/README.rst b/README.rst index 1660218..c5cdb79 100644 --- a/README.rst +++ b/README.rst @@ -286,7 +286,7 @@ Surface patches prodes Fv.pdb Fv.zip --patches --residue-set CDR=cdr_residues.txt -The default 54 columns are unchanged by the flag; the patch columns are appended after them. The thresholds and the minimum patch area are settings with starting values rather than calibrated choices. Method, settings, columns and limitations: `surface patches `_. +The default 54 columns are unchanged by the flag; the patch columns are appended after them. The thresholds and the minimum patch area are settings; their defaults were checked against antibody HIC data and kept, because settings tuned on one antibody panel did worse on a held-out one. Method, settings, columns and limitations: `surface patches `_. Looking features up from Python ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/surface_patches.md b/docs/surface_patches.md index c6e2d66..9ecaa3c 100644 --- a/docs/surface_patches.md +++ b/docs/surface_patches.md @@ -43,7 +43,7 @@ Finding patches on the whole surface and then restricting them means a patch tha | Setting | CLI | Default | Notes | | --- | --- | --- | --- | | `radius` | `--patch-radius` | 5.0 A | Radius of the local average, at most 15 A. The cost of the average grows with the cube of the radius. | -| `min_area` | `--patch-min-area` | 100.0 A^2 | Smallest patch kept. A starting value; see below. | +| `min_area` | `--patch-min-area` | 100.0 A^2 | Smallest patch kept; see below. | | `hydrophobic_threshold` | `--patch-hydrophobic-threshold` | 0.0 | The sign of the local hydrophobicity. Refused for a scale with no negative values, where 0 would mark every point hydrophobic; pass a threshold explicitly there. | | `ep_threshold` | `--patch-ep-threshold` | 0.5 V | Applied symmetrically to positive and negative patches. | | `residue_sets` | `--residue-set NAME=FILE` | none | Repeatable. `NAME` is letters and digits only. | @@ -52,11 +52,15 @@ The `--patch-*` and `--residue-set` options are refused without `--patches`, sin The hydrophobicity scale is the run's `--hydro` setting. The thresholds are in the units of the local value, so a threshold chosen for one scale does not carry over to another. -The defaults are starting values to be tuned against measured data, not calibrated choices: +The defaults were checked against antibody HIC retention time on the GDPa1 panel (236 antibodies). Settings tuned on that panel raised the correlation there but not on the held-out GDPa3 panel (79 antibodies), so the defaults were kept. Any re-tuning needs an independent panel to check it on. + +The tuning searched radius 3 to 12 A, minimum area 0 to 500 A^2, hydrophobic threshold -0.08 to 0.06 and electrostatic threshold 0.1 to 1.0 V, with patches on ABodyBuilder3 Fv models and residue sets for the CDRs and the Fv. The best setting on GDPa1 (radius 3 A, minimum area 0, hydrophobic threshold -0.06, electrostatic threshold 0.1 V) raised the Spearman R of `SurfMhpPatchArea_CDR` against retention time from 0.34 to 0.46 there, but lowered it from 0.355 to 0.328 on GDPa3. A random-forest HIC model on sequence features plus the Prodes columns, fitted on GDPa1, scored a Spearman R of 0.589 on GDPa3 with the tuned setting against 0.615 with the defaults, lower on all five random-forest seeds. Choosing the setting inside each training fold of GDPa1 instead gave 0.579 against 0.581 on its own five published folds. + +How each default was arrived at: - **Minimum area.** Areas here are solvent accessible surface areas with a 1.4 A probe, which are larger than molecular surface areas of the same region, and by an amount that depends on curvature, so a minimum area quoted by a tool that measures on a molecular surface does not transfer to this one. 100 A^2 is about half the solvent accessible area of one fully exposed phenylalanine residue, 210 to 240 A^2 (Tien et al., PLoS ONE 2013): a patch has to be larger than the exposed half of one large hydrophobic side chain. - **Electrostatic threshold.** At 0 V, the surface with the same sign as the net charge becomes one or two patches covering most of the protein, and the patch area then says little beyond `NSurfPosEpFormal` or the total area. 0.5 V was chosen from the local potential at the default ionic strength of 0.15 M on the eight structures in `tests/data` (1AO6, 1CBN, 1GDW, 1GPB, 1PIT, 4NZU, ARH96693 and ARH98503): over those structures the median of the 75th percentile of the absolute local potential was 0.39 V and the median of the 90th percentile 0.61 V. At 0.5 V the marked points covered between 13 and 26 per cent of the surface area on six of the eight structures, 2 per cent on 1CBN and 32 per cent on 1PIT. The potential is a screened Coulomb sum at a relative permittivity of 4, so thresholds from Poisson-Boltzmann tools, usually given in kT/e, do not carry over. -- **Hydrophobic threshold.** 0 on `mj_scaled` separates hydrophobic from hydrophilic residues. It is lower than the 0.03 the hydrophobicity viewer uses to colour points, because the viewer shows the per-point value and the patches use the local value, whose range is narrower. The local value on `mj_scaled` is negative over most of a typical surface, so at 0 only the few most hydrophobic regions pass, and the hydrophobic patch area changes steeply with the threshold. Choose it against measured data rather than relying on the default. +- **Hydrophobic threshold.** 0 on `mj_scaled` separates hydrophobic from hydrophilic residues. It is lower than the 0.03 the hydrophobicity viewer uses to colour points, because the viewer shows the per-point value and the patches use the local value, whose range is narrower. The local value on `mj_scaled` is negative over most of a typical surface, so at 0 only the few most hydrophobic regions pass, and the hydrophobic patch area changes steeply with the threshold. All settings, with the hydrophobic threshold resolved and every residue set's labels, are recorded in `prodes_run.json` under `settings.patches`. diff --git a/src/prodes/calculations/patches.py b/src/prodes/calculations/patches.py index 5ab7c09..f7eecf4 100644 --- a/src/prodes/calculations/patches.py +++ b/src/prodes/calculations/patches.py @@ -27,7 +27,7 @@ #: Patches smaller than this, in square angstrom of solvent accessible surface, #: are discarded: about half the solvent accessible area of one fully exposed -#: phenylalanine. A starting value; see docs/surface_patches.md. +#: phenylalanine. Checked against antibody HIC data and kept; see docs/surface_patches.md. DEFAULT_PATCH_MIN_AREA = 100.0 #: The largest radius accepted. The neighbourhood average reads the whole point From 31fe8942166c3d56300a35a88165627fed2a06ba Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Mon, 19 Oct 2026 10:55:25 +0200 Subject: [PATCH 4/4] #35 Calculate the surface patch features on every run and add them to the default set Added to a random-forest model of antibody hydrophobic interaction chromatography retention built on sequence features, the 11 surface patch features raised the Spearman R on the held-out GDPa3 panel from 0.527 to 0.620 and improved every one of the five published GDPa1 folds, and they cost about a second per structure. They were nevertheless only written under --patches, so a caller had to know to ask for them. The patch step now runs on every structure and --patches is removed. The 11 patch entries move from features_patches.yaml, which is deleted, to the end of features_reduced.yaml, so the default feature set is 65 columns and --full-features writes 116; each is marked per_residue_set, and SurfArea moves with them under residue_set_only_features. The patch options (--patch-radius, --patch-min-area, --patch-hydrophobic-threshold, --patch-ep-threshold, --residue-set) now always apply, and calculate(patches=None) means PatchSettings() rather than no patches. prodes_run.json always records the patch settings, and the surface point table always carries local_mhp, local_ep_volts and mhp_patch_id. FeatureDictionary.get_feature_codes takes residue_sets to list the columns of a run with residue sets, and get_patch_feature_codes returns the 11 patch codes. This is a breaking change, so the version is 10.0.0: the default feature table has 11 more columns, a --full-features table no longer matches the original 105 columns, drop_redundant_features requires the patch columns, and a run with --hydro mj, whose values are all positive, is refused unless --patch-hydrophobic-threshold is given. No earlier column changes value. The tests that pinned the original 105 columns now name the patch codes as the one deliberate addition, and the README, docs/surface_patches.md and docs/redundant_feature_analysis.md describe the new default set, noting that the patch features were not part of the redundancy analysis. A docstring in tests/test_pka.py no longer names a particular downstream caller. --- README.rst | 56 +++++++-------- docs/redundant_feature_analysis.md | 2 + docs/structure_preparation.md | 2 +- docs/surface_patches.md | 16 +++-- pyproject.toml | 2 +- src/prodes/__init__.py | 6 +- src/prodes/calculations/patches.py | 7 +- src/prodes/data/features_patches.yaml | 99 --------------------------- src/prodes/data/features_reduced.yaml | 95 ++++++++++++++++++++++++- src/prodes/feature_dictionary.py | 67 +++++++++--------- src/prodes/output.py | 8 +-- src/prodes/run.py | 95 +++++++++++-------------- tests/test_cli.py | 38 +++------- tests/test_feature_dictionary.py | 20 ++++-- tests/test_output_bundle.py | 3 +- tests/test_patches.py | 51 +++++++------- tests/test_pka.py | 6 +- tests/test_sasa.py | 8 ++- 18 files changed, 279 insertions(+), 302 deletions(-) delete mode 100644 src/prodes/data/features_patches.yaml diff --git a/README.rst b/README.rst index c5cdb79..b8e3fef 100644 --- a/README.rst +++ b/README.rst @@ -1,7 +1,7 @@ Prodes: Protein Descriptors ============================ -**Prodes turns a 3D protein structure into numbers describing the surface properties.** You give it a PDB file, and it gives you a row in a CSV file with 54 columns describing the surface of that protein: how much of it is charged, how much of it is hydrophobic, how those properties are spread out across the surface, and how the picture changes with pH. +**Prodes turns a 3D protein structure into numbers describing the surface properties.** You give it a PDB file, and it gives you a row in a CSV file with 65 columns describing the surface of that protein: how much of it is charged, how much of it is hydrophobic, how those properties are spread out across the surface, how large the contiguous hydrophobic and charged patches are, and how the picture changes with pH. Those columns are designed to be used directly as the input features of a machine learning or QSPR model. If you have measured something about a set of proteins (a retention time, an aggregation onset, a viscosity, a titre) and you want to predict it for proteins you have not measured, Prodes gives you the X side of that problem from structure alone. @@ -54,7 +54,7 @@ The fourth line adds `PROPKA `_, which is PDB2PQR, in its own environment ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -You should also use `PDB2PQR `_, which repairs the structure before anything measures it. On one antibody structure, rebuilding two unmodelled lysine side chains moved 31 of the 54 features and shifted the isoelectric point by more than a pH unit. See `Preparing your structure`_. +You should also use `PDB2PQR `_, which repairs the structure before anything measures it. On one antibody structure, rebuilding two unmodelled lysine side chains moved 31 of the 54 non-patch default features and shifted the isoelectric point by more than a pH unit. See `Preparing your structure`_. **Put it in an environment of its own:** @@ -126,7 +126,7 @@ Run these four commands on your structure. PDB2PQR lives in its own environment, python -m prodes.io.pka_converter 1GDW.pka propka -o 1GDW_pka.json # convert them for Prodes python -m prodes prepared/1GDW.pdb 1GDW.zip --ph 7.4 --pka 1GDW_pka.json # calculate the features -That writes ``1GDW.zip``, a bundle holding the 54 features, the surface points they were calculated from, and ready-to-open viewer scripts. On a small protein the whole thing takes seconds. +That writes ``1GDW.zip``, a bundle holding the 65 features, the surface points they were calculated from, and ready-to-open viewer scripts. On a small protein the whole thing takes seconds. Two things about that sequence are worth noticing, and both are explained under `Preparing your structure`_: @@ -148,7 +148,8 @@ The ones you are most likely to want: * ``--ph`` — the pH at which to compute protonation states (default 7). Charge features change substantially with this, so set it to the pH of your buffer. * ``--pka`` — supply per-residue pKa values from PROPKA or a similar tool. Recommended, as above; see `pKa values and protonation states`_. * ``--probe`` — the radius of the solvent probe used for the surface calculation (default 1.4 Å, i.e. water). -* ``--full-features`` — write the original 105-feature set instead of the reduced 54. See `The reduced feature set`_. +* ``--full-features`` — write the original 105 features plus the 11 patch features, 116 in all, instead of the reduced 65. See `The reduced feature set`_. +* ``--residue-set NAME=FILE`` — also report every surface patch feature restricted to the residues listed in the file, such as the CDRs of an antibody. See `Surface patches`_. * ``--n-workers``, ``--chunksize``, ``--mem-limit`` — CPU and memory tuning. You can ignore all three; the defaults are sensible. See `Speed`_. Viewing the surface @@ -227,7 +228,7 @@ Unpacked, the bundle holds: .. code-block:: text 1GDW/ - 1GDW_features.csv the 54 features, one row, first column ID + 1GDW_features.csv the 65 features, one row, first column ID 1GDW_surface_points.csv every surface point: x, y, z, potential, hydrophobicity, area, and the residue owning it 1GDW_ep.pml PyMOL, coloured by electrostatic potential @@ -262,7 +263,7 @@ To read the features back: What the features actually are ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Prodes builds a dotted surface over the protein (a Shrake-Rupley solvent-accessible surface), assigns each surface point an electrostatic potential and a hydrophobicity, and then summarises those distributions. The 54 default columns are, broadly: +Prodes builds a dotted surface over the protein (a Shrake-Rupley solvent-accessible surface), assigns each surface point an electrostatic potential and a hydrophobicity, and then summarises those distributions. The 65 default columns are, broadly: * **Whole-molecule properties** — molecular weight, total surface area, formal charge at your chosen pH, isoelectric point, dipole moment, and shape descriptors. * **Surface electrostatic potential** — the maximum, minimum, mean and standard deviation of the potential over the whole surface, then the same statistics computed separately over just the positive and just the negative regions, plus a count of how many surface points are positive. A protein with a strongly positive area of surface and a strongly negative one can have a near-zero net charge but very distinctive values here, which is exactly the kind of thing a net-charge calculation misses. Note that these are statistics over all positive or all negative surface points; grouping them into individual patches is the optional patch step (see `Surface patches`_). @@ -272,21 +273,20 @@ Prodes builds a dotted surface over the protein (a Shrake-Rupley solvent-accessi Every feature has a plot label, a one-line description, a full explanation of how it is calculated, its unit, and the wording used for it in the original publication. These live in two human-readable YAML files, which are also the single source of truth for which features Prodes writes, so you can read the whole list without running anything: -* `features_reduced.yaml `_ — the **54 features calculated by default** +* `features_reduced.yaml `_ — the **65 features calculated by default**, the last 11 of them the surface patch features * `features_full_only.yaml `_ — the **other 51**, written only under ``--full-features``, each with the reason it is not in the default set -A third file, `features_patches.yaml `_, documents the surface patch features, which are written only with ``--patches``. Surface patches ~~~~~~~~~~~~~~~ -``--patches`` adds a block of features that group the surface into contiguous regions. Each point's hydrophobicity and potential are averaged over a 5 Å neighbourhood, points above a threshold that touch are joined into patches, and the run reports the summed area, the largest area and the number of hydrophobic, positive and negative patches, plus the highest and the summed positive local hydrophobicity. ``--residue-set NAME=FILE`` reports every one of them again restricted to the residues listed in the file, such as the CDRs of an antibody, under columns ending in ``_NAME``: +The last 11 default columns group the surface into contiguous regions. Each point's hydrophobicity and potential are averaged over a 5 Å neighbourhood, points above a threshold that touch are joined into patches, and the run reports the summed area, the largest area and the number of hydrophobic, positive and negative patches, plus the highest and the summed positive local hydrophobicity. ``--residue-set NAME=FILE`` reports every one of them again restricted to the residues listed in the file, such as the CDRs of an antibody, under columns ending in ``_NAME``: .. code-block:: bash - prodes Fv.pdb Fv.zip --patches --residue-set CDR=cdr_residues.txt + prodes Fv.pdb Fv.zip --residue-set CDR=cdr_residues.txt -The default 54 columns are unchanged by the flag; the patch columns are appended after them. The thresholds and the minimum patch area are settings; their defaults were checked against antibody HIC data and kept, because settings tuned on one antibody panel did worse on a held-out one. Method, settings, columns and limitations: `surface patches `_. +The patch columns come after all the others, and the residue-set columns after those. Added to a random-forest model of antibody hydrophobic interaction chromatography retention built on sequence features, the patch columns raised its Spearman R on a held-out panel from 0.527 to 0.620. The thresholds and the minimum patch area are settings; their defaults were checked against antibody HIC data and kept, because settings tuned on one antibody panel did worse on a held-out one. Method, settings, columns and limitations: `surface patches `_. Looking features up from Python ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -308,7 +308,7 @@ The same information is available programmatically, which is what you want when fd.get_reason_dropped("SurfEpMeanAverage") # why it is not in the default set fd.get_reason_dropped("SurfEpMeanFormal") # None: this one is kept -The lookups are one-directional by design: the feature code is the key, and there is no route back from a label or description to a code. ``fd.get_entry(code)`` returns everything known about one feature and ``fd.get_dictionary()`` returns all 105 entries at once. +The lookups are one-directional by design: the feature code is the key, and there is no route back from a label or description to a code. ``fd.get_entry(code)`` returns everything known about one feature and ``fd.get_dictionary()`` returns all 116 entries at once. Getting the feature codes ^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -321,8 +321,8 @@ Both lists come from the YAML dictionaries shipped in ``prodes/data``, so a pipe fd = FeatureDictionary() - fd.get_feature_codes() # the 54 written by default - fd.get_feature_codes(full=True) # all 105 legacy columns + fd.get_feature_codes() # the 65 written by default + fd.get_feature_codes(full=True) # the 116 written with --full-features fd.get_dropped_feature_codes() # the 51 the default leaves out Each returns a fresh ``list`` of feature codes, in the order Prodes writes the columns. ``ID`` is a row label rather than a feature and is in neither list; it is available as ``prodes.feature_dictionary.ID_COLUMN``. The underlying tuples are exposed as ``FULL_FEATURE_CODES``, ``REDUCED_FEATURE_CODES`` and ``DROPPED_FEATURE_CODES`` if you would rather not construct the class. @@ -418,9 +418,9 @@ Using ``validate="one_to_one"`` is worth the extra keystrokes: it raises if a st PLS is a reasonable first choice for this kind of data because the features are correlated with one another by construction, and PLS handles that gracefully. Random forests and gradient boosting also work well and need no scaling. -**A note on sample size.** 54 features is a lot if you have measured 20 proteins. With fewer observations than features, almost any model will fit the training data perfectly and predict nothing. The reduced feature set exists precisely to push that ratio in your favour; do not reach for ``--full-features`` to get more columns unless you have the observations to support them (see `The reduced feature set`_). Always report a score against a holdout blind test set, never a training-set score. +**A note on sample size.** 65 features is a lot if you have measured 20 proteins. With fewer observations than features, almost any model will fit the training data perfectly and predict nothing. The reduced feature set exists precisely to push that ratio in your favour; do not reach for ``--full-features`` to get more columns unless you have the observations to support them (see `The reduced feature set`_). Always report a score against a holdout blind test set, never a training-set score. -**If you select features, select them inside the cross-validation loop.** Screening all 54 columns for the ones that correlate with your measurement, and then cross-validating a model built from the survivors, reports a score that is not out of sample: every held-out protein helped choose the features. On small datasets this reliably manufactures impressive numbers from nothing. ``sklearn.feature_selection.SelectKBest`` inside a ``Pipeline`` does it correctly, because the selection is refitted on each training fold. +**If you select features, select them inside the cross-validation loop.** Screening all 65 columns for the ones that correlate with your measurement, and then cross-validating a model built from the survivors, reports a score that is not out of sample: every held-out protein helped choose the features. On small datasets this reliably manufactures impressive numbers from nothing. ``sklearn.feature_selection.SelectKBest`` inside a ``Pipeline`` does it correctly, because the selection is refitted on each training fold. Ionic strength and screening ----------------------------- @@ -454,7 +454,7 @@ The values are not quite identical to the released 4.x numbers, for a separate r **What this is and is not.** The screening length is the Debye length in water, where the relative permittivity is about 78.5. Prodes evaluates its sum at a uniform relative permittivity of 4, where the self consistent value would be about 1.8 Angstrom. Using the water value inside a kernel with a protein dielectric is an empirical correction that mimics screening. It is used because it is what was measured to agree best with a Poisson-Boltzmann reference, not because the two are consistent, and the potential is still not comparable to an APBS calculation despite both being reported in volts. -**Scope.** Screening applies to the surface features, the ``SurfEp`` family: 9 of the 54 default features, and 38 of the 105 with ``--full-features``, the extra 19 being the ``SurfEp*Average`` columns computed from partial rather than formal charges. The ``ShellEp`` features, 9 by default and 19 with the full set, are unchanged, and that is a measured decision rather than an omission: they are computed by a different route which divides the path at the molecular surface and weights the solvent leg with a permittivity of 80 against 4 for the protein, so a distant charge is already damped about twentyfold and the offset that affected the surface never built up. Checked against an APBS equivalent, the unscreened shell agrees at a Spearman of 0.877, and adding screening moved that to 0.855. See ``docs/screening_validation.md``. +**Scope.** Screening applies to the surface features, the ``SurfEp`` family: 9 of the 65 default features, and 38 of the 116 with ``--full-features``, the extra 19 being the ``SurfEp*Average`` columns computed from partial rather than formal charges. The ``ShellEp`` features, 9 by default and 19 with the full set, are unchanged, and that is a measured decision rather than an omission: they are computed by a different route which divides the path at the molecular surface and weights the solvent leg with a permittivity of 80 against 4 for the protein, so a distant charge is already damped about twentyfold and the offset that affected the surface never built up. Checked against an APBS equivalent, the unscreened shell agrees at a Spearman of 0.877, and adding screening moved that to 0.855. See ``docs/screening_validation.md``. Preparing your structure ------------------------- @@ -496,7 +496,7 @@ Mean surface electrostatic potential -0.013 + Positive surface points 12,976 13,684 ===================================== ======================= ===================== -**31 of the 54 default features move.** The net charge is wrong by 2, the isoelectric point is wrong by 1.14 units, and the mean surface potential has the wrong sign. For choosing between an anion and a cation exchanger at a given buffer pH, a pI of 6.46 and a pI of 7.61 are different recommendations. +**31 of the 54 non-patch default features move.** The net charge is wrong by 2, the isoelectric point is wrong by 1.14 units, and the mean surface potential has the wrong sign. For choosing between an anion and a cation exchanger at a given buffer pH, a pI of 6.46 and a pI of 7.61 are different recommendations. This is not universal. Many deposited structures are complete. It is a property of resolution and of surface disorder, and the residues most often left unmodelled are the long charged ones, which are exactly the residues that carry the charge. @@ -548,7 +548,7 @@ Left to itself, Prodes uses **one textbook pKa per residue type**: every asparta **PROPKA predicts a pKa for each individual residue from its structural environment**, and Prodes reads those predictions. This is the recommended default. It is a small, free, MIT-licensed program from the Jensen group, it installs on Windows, macOS and Linux, and on a normal protein it finishes in seconds. There is very little reason not to use it. -How much does it matter? On 1GDW, feeding in PROPKA values changes **19 of the 54 features**, moving the isoelectric point from 10.34 to 10.57 and the formal charge from +7 to +8. In a larger test across 819 AlphaFold structures it changed 22 features, the most affected being ``SurfEpMinFormal`` at only R² = 0.82 between the two runs (`full report `_). +How much does it matter? On 1GDW, feeding in PROPKA values changes **19 of the 54 non-patch default features**, moving the isoelectric point from 10.34 to 10.57 and the formal charge from +7 to +8. In a larger test across 819 AlphaFold structures it changed 22 features, the most affected being ``SurfEpMinFormal`` at only R² = 0.82 between the two runs (`full report `_). .. important:: @@ -647,7 +647,7 @@ The command line and the Python interface do the same work. To reproduce a comma import prodes prodes.run_prodes("./tests/data/1GDW.pdb.zip", "example.zip") -The full signature is ``run_prodes(pdb_file, out_file, pkas_file=None, ph=7, r_probe=1.4, hydro_scale="mj_scaled", full_features=False, mem_limit_mb=None, ionic_strength_molar=0.15, patches=None)``, where ``patches`` takes a ``prodes.calculations.patches.PatchSettings``; see `Surface patches`_. +The full signature is ``run_prodes(pdb_file, out_file, pkas_file=None, ph=7, r_probe=1.4, hydro_scale="mj_scaled", full_features=False, mem_limit_mb=None, ionic_strength_molar=0.15, patches=None)``, where ``patches`` takes a ``prodes.calculations.patches.PatchSettings`` to change the patch settings or name residue sets, and None means the defaults; see `Surface patches`_. The lower-level pieces are importable too, if you want a single property rather than the whole feature set. Calculating just the surface area, for example: @@ -672,7 +672,7 @@ The lower-level pieces are importable too, if you want a single property rather The reduced feature set ------------------------ -The 105 numeric features in the original Prodes output were highly redundant. They have been reduced to **54 features** to lower the risk of overfitting and to cut calculation time. Nothing new is calculated and no kept feature changes value: the default output is a strict column subset of the original, in the original column order. +The 105 numeric features in the original Prodes output were highly redundant. They have been reduced to **54 features** to lower the risk of overfitting and to cut calculation time, and the default output adds the 11 surface patch features after them, **65** in all. No kept feature changes value: apart from the patch columns at the end, the default output is a strict column subset of the original, in the original column order. The patch features were not part of the redundancy analysis. To calculate the full set anyway, use the ``--full-features`` flag, set ``PRODES_FULL_FEATURES=true`` in the environment, or pass ``full_features=True`` to ``prodes.run.calculate()``. @@ -817,9 +817,9 @@ These two come up constantly in this literature, but neither competes with Prode Where Prodes sits ~~~~~~~~~~~~~~~~~~ -**By default, Prodes does not identify discrete patches.** It builds a dotted surface, gives every point an electrostatic potential and a hydrophobicity, and then summarises those values over the whole surface and over the positive and negative subsets. Those subsets are defined by the sign of the value at a point, not by spatial contiguity: a positive-subset statistic pools every positive surface point on the protein, whether they form one large region or twenty scattered ones. +**Prodes identifies discrete patches, in its last 11 default features.** It builds a dotted surface and gives every point an electrostatic potential and a hydrophobicity. Most of its features then summarise those values over the whole surface and over the positive and negative subsets, which are defined by the sign of the value at a point, not by spatial contiguity: a positive-subset statistic pools every positive surface point on the protein, whether they form one large region or twenty scattered ones. -With ``--patches`` it does: neighbourhood-averaged values are thresholded and touching points joined into patches, whose areas are reported over the whole surface and over any residue set you name, and the point table records which hydrophobic patch each point belongs to and which residue owns it. See `Surface patches`_. The patches are found on a solvent-accessible point cloud rather than a triangulated molecular surface, so their areas are not directly comparable with those of mesh-based tools such as PEP-Patch. +The patch features add the contiguity: neighbourhood-averaged values are thresholded and touching points joined into patches, whose areas are reported over the whole surface and over any residue set you name, and the point table records which hydrophobic patch each point belongs to and which residue owns it. See `Surface patches`_. The patches are found on a solvent-accessible point cloud rather than a triangulated molecular surface, so their areas are not directly comparable with those of mesh-based tools such as PEP-Patch. It is worth being clear about what this does *not* mean, though, because the word "patch" is used very loosely in this field. **Prodes is not blind to where things sit on the surface.** Its molecular hydrophobicity potential gives each surface point an ``exp(-d)``-weighted sum over the non-hydrogen atoms within 10 Å. The weight has a decay length of 1 Å, so the 10 Å cutoff hardly matters and the value mostly reflects the one or two nearest atoms, a much shorter length scale than SAP, Genentech's HPATCH (a 10 Å average) or Protein-Sol's NPP ratio (a 13 Å ratio). The patch step's 5 Å neighbourhood average is what brings it to the scale of a surface region. Only three of the tools above actually segment the surface into discrete regions: PEP-Patch (connected components on the surface graph), PROPERMAB (DBSCAN on mesh triangles) and MOE's Protein Patch Analyzer. The rest, TAP and TNP included, compute neighbourhood-weighted quantities and then aggregate them, which is what Prodes does too. The difference is in the last step: they take a maximum or a region sum, the default Prodes features take the distribution statistics, and the patch features take a maximum and a sum, over the whole surface or a residue set. @@ -846,27 +846,27 @@ Improvements in this fork * **Vectorised shell feature computation** — ``find_exit``, ``project_point``, and ``map_ep_to_plane`` have batch counterparts (``_batch``) that process all charged atoms simultaneously with ``np.einsum`` and broadcasting (~10x speedup on the shell phase, previously the slowest part of the pipeline). * **Vectorised surface grid construction** — grid construction and cell filling refactored to use NumPy arrays throughout. * **Multi-core parallelism** — the SASA, surface grid and shell phases are spread across worker processes on Linux. -* **Reduced feature set** — 54 non-redundant features by default, with the original 105 available via ``--full-features``. +* **Reduced feature set** — 54 non-redundant features of the original 105 by default, with all 105 available via ``--full-features``. * **Disulfide bonds recognised** — a cysteine bonded into a disulfide is no longer titrated as a free thiol, which corrected the charge-derived features of every structure with a disulfide. See `cysteines and disulfide bonds `_. * **Alternate conformations resolved** — the ``altLoc`` column used to be read as part of the atom name, so a disordered residue kept every conformation in the surface and, because the package looks up charge by atom name, silently carried no charge at all. From version 7.0 one conformation per residue is kept and the atom name is read on its own. See `alternate conformations `_. * **PDB records read by Biopython** — from version 7.1 the coordinate records of a PDB file are read by ``Bio.PDB.PDBParser``, driven with a record-collecting builder of prodes' own rather than through Biopython's entity tree, which represents disorder rather than resolving it. No feature value moves on any well-formed structure. The insertion code and the model number now arrive as their own fields, which is what the two parser defects still open at that release needed; both were fixed in 8.0. See `how prodes reads a PDB file `_. * **Insertion codes and NMR models** — from version 8.0 a residue is its chain, its number **and** its insertion code, and a file holding several models is described by its first. Kabat and Chothia numbering write every antibody CDR insertion as an insertion code, and reading ``H100`` and ``H100A`` as one residue put both side chains in one residue: on the 4NZU Fab that lost 12 of its 434 residues, moved the molecular weight by 1.3 kDa, titrated a lysine and a cysteine against aspartate's pKa, and hid a cysteine well enough that one of its five disulfide bonds went undetected. An NMR ensemble read whole was worse: 1PIT came out as 58 residues holding 17 780 atoms, 16 901 of them in a single residue, at a formal charge of -1096 for a protein whose real charge at pH 7 is about +6. See `what Prodes calls one residue `_. * **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 9.1, ``--patches`` groups neighbourhood-averaged hydrophobicity and potential into contiguous patches and reports their areas over the whole surface and over named residue sets, such as antibody CDRs. Off by default; the default columns do not change. Every surface point in the bundle now also carries its area and the residue owning it. See `surface patches `_. +* **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 `_. * **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. Output compatibility --------------------- -With ``--full-features``, the output has the same 105 columns in the same order as the original Prodes. +With ``--full-features``, the output has the same 105 columns in the same order as the original Prodes, followed from version 10.0 by the 11 surface patch columns. 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 9.1 the surface point table in the bundle has two more columns, ``area_A2`` and ``residue``, on every run. The feature table is unchanged unless ``--patches`` is given, which appends the patch columns after the others. +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 54 columns they share. +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. How to cite ------------ diff --git a/docs/redundant_feature_analysis.md b/docs/redundant_feature_analysis.md index bf7bbf2..442cb42 100644 --- a/docs/redundant_feature_analysis.md +++ b/docs/redundant_feature_analysis.md @@ -138,6 +138,8 @@ the README for how to enable them and when that is advisable. Every distinct physical property the original calculates is still represented. What is gone is the second, third and fourth summary statistic of the same distribution, and the second charge model applied to the same surface. +From version 10.0 the default set also holds the 11 surface patch features after these 54, 65 in all; see `docs/surface_patches.md`. They were added because they improved a model of antibody hydrophobic interaction chromatography retention on a held-out panel, and were not part of the R² analysis on this page, so whether any of them is redundant with the 54 has not been measured. + ### Why the reduced set is the default The original package calculates all 105 features and leaves feature selection to the user. Tim Neijenhuis, the original author, makes the fair point that anyone building a model can simply pick the subset they need later, and for an experienced modeller that is exactly right. diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index 0ce7494..f2fc594 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -47,7 +47,7 @@ Only the **geometry**. Prodes reads `ATOM` records, filters out hydrogens, and r So the hydrogens PDB2PQR adds, the protonation states it assigns and the partial charges it writes are all discarded. That is deliberate: it is what lets you run Prodes at several pH values against one prepared structure. What you gain from the preparation step is that the atoms are there and in sensible positions. -A structure without hydrogens is therefore a valid input. Stripping all 1,625 hydrogens from an ABodyBuilder3 Fv model, with the same pKa file, left every one of the 77 feature columns of a `--patches` run unchanged. What does change the values is a missing heavy atom, and on a predicted structure that is almost always `OXT`; see the next section. +A structure without hydrogens is therefore a valid input. Stripping all 1,625 hydrogens from an ABodyBuilder3 Fv model, with the same pKa file, left every one of the 77 feature columns of a run with a CDR residue set unchanged. What does change the values is a missing heavy atom, and on a predicted structure that is almost always `OXT`; see the next section. ## What happens to a predicted structure: the terminal oxygen diff --git a/docs/surface_patches.md b/docs/surface_patches.md index 9ecaa3c..c2625a3 100644 --- a/docs/surface_patches.md +++ b/docs/surface_patches.md @@ -1,18 +1,20 @@ # Surface patches -The default Prodes features summarise the surface by the sign of each point's value: every hydrophobic point on the protein is pooled, whether the points form one region or twenty scattered ones. The patch step adds the spatial grouping. It reports how much of the surface is covered by contiguous hydrophobic, positive and negative regions, how large the largest one is, and the same quantities restricted to a named set of residues, such as the CDRs of an antibody. +Most Prodes features summarise the surface by the sign of each point's value: every hydrophobic point on the protein is pooled, whether the points form one region or twenty scattered ones. The patch step adds the spatial grouping. It reports how much of the surface is covered by contiguous hydrophobic, positive and negative regions, how large the largest one is, and the same quantities restricted to a named set of residues, such as the CDRs of an antibody. -The step is off by default. Switching it on appends columns to the feature table and leaves every default column unchanged. +The patch step runs on every structure. Its 11 columns are the last of the 65 default features, and the residue-set columns, when a run is given residue sets, come after them. It adds about 1 s to a run of 57,000 surface points. + +Added to a random-forest model of antibody hydrophobic interaction chromatography (HIC) retention built on sequence features, the patch columns raised the mean per-fold Spearman R on the five published folds of the GDPa1 panel (236 antibodies) from 0.533 to 0.573 with ABodyBuilder3 Fv models and from 0.529 to 0.589 with Boltz-2 models of the whole antibody, improving every fold, and raised it on the held-out GDPa3 panel (79 antibodies) from 0.527 to 0.620. Added to models of anion and cation exchange retention, they changed the Spearman R by less than 0.01. ## Running it Command line: ```bash -prodes Fv.pdb Fv.zip --patches --residue-set CDR=cdr_residues.txt --residue-set Fv=fv_residues.txt +prodes Fv.pdb Fv.zip --residue-set CDR=cdr_residues.txt --residue-set Fv=fv_residues.txt ``` -Python: +Python, where `patches` is optional and only needed to name residue sets or change a setting: ```python from prodes.run import calculate @@ -48,7 +50,7 @@ Finding patches on the whole surface and then restricting them means a patch tha | `ep_threshold` | `--patch-ep-threshold` | 0.5 V | Applied symmetrically to positive and negative patches. | | `residue_sets` | `--residue-set NAME=FILE` | none | Repeatable. `NAME` is letters and digits only. | -The `--patch-*` and `--residue-set` options are refused without `--patches`, since they would otherwise do nothing. +With `--hydro mj`, whose values are all positive, the default hydrophobic threshold of 0 would mark every point hydrophobic, so a run with that scale is refused unless `--patch-hydrophobic-threshold` is given. The hydrophobicity scale is the run's `--hydro` setting. The thresholds are in the units of the local value, so a threshold chosen for one scale does not carry over to another. @@ -82,9 +84,9 @@ For each residue set `NAME` the run then writes `SurfArea_NAME`, the surface are `SurfMhpLocalMax` and `SurfMhpLocalPosSum` need no threshold and no minimum area, so they are not affected by either choice. They are similar in spirit to spatial aggregation propensity (Chennamsetty et al., PNAS 2009), which also averages a hydrophobicity scale over a neighbourhood and sums the positive part, although that method works per atom and weights each residue by its exposed side-chain fraction. `SurfMhpLocalPosSum` is in square angstrom times the units of the hydrophobicity scale, so its value depends on the scale. -Every code is documented in `src/prodes/data/features_patches.yaml`, and `FeatureDictionary` looks up a code with or without a residue set suffix. +Every code is documented in `src/prodes/data/features_reduced.yaml`, where the patch features are marked `per_residue_set`. `FeatureDictionary` looks up a code with or without a residue set suffix, and `get_feature_codes(residue_sets=[...])` lists the columns of a run with residue sets. -The point table, `*_surface_points.csv`, always carries `area_A2` and `residue` for each point. With the patch step it also carries `local_mhp`, `local_ep_volts` and `mhp_patch_id`, the hydrophobic patch the point belongs to (-1 for none), so a patch can be traced back to its residues. +The point table, `*_surface_points.csv`, carries `area_A2`, `residue`, `local_mhp`, `local_ep_volts` and `mhp_patch_id`, the hydrophobic patch the point belongs to (-1 for none), for each point, so a patch can be traced back to its residues. ## Limitations diff --git a/pyproject.toml b/pyproject.toml index 0addc1c..b31946f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "prodes" -version = "9.1.0" +version = "10.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 4b217e6..9a8b5f5 100644 --- a/src/prodes/__init__.py +++ b/src/prodes/__init__.py @@ -37,9 +37,9 @@ def run_prodes( mem_limit_mb: memory budget in MB for the intermediate NumPy arrays of this run, divided among the workers. Pass explicitly in worker processes rather than relying on the PRODES_MEM_LIMIT_MB environment variable. - patches: a prodes.calculations.patches.PatchSettings to also report - contiguous surface patches, optionally over named residue sets. - None leaves them out. + 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. """ calculate( diff --git a/src/prodes/calculations/patches.py b/src/prodes/calculations/patches.py index f7eecf4..f407e61 100644 --- a/src/prodes/calculations/patches.py +++ b/src/prodes/calculations/patches.py @@ -50,7 +50,7 @@ @dataclass(frozen=True) class PatchSettings: - """The settings of the patch step. Passing one to calculate switches the step on. + """The settings of the patch step, which runs on every structure. calculate uses PatchSettings() when given none. Attributes: radius: neighbourhood radius in angstrom for the average. @@ -107,7 +107,8 @@ def resolved_hydrophobic_threshold(self, hydro_scale): if min(hydrophobic_scale(hydro_scale).values()) >= 0: raise ValueError( f"the hydrophobicity scale {hydro_scale} has no negative values, so the default patch threshold of " - f"{DEFAULT_HYDROPHOBIC_THRESHOLD} would mark every point hydrophobic; pass a hydrophobic threshold explicitly" + f"{DEFAULT_HYDROPHOBIC_THRESHOLD} would mark every point hydrophobic; pass one explicitly, with " + "--patch-hydrophobic-threshold or PatchSettings(hydrophobic_threshold=...)" ) return DEFAULT_HYDROPHOBIC_THRESHOLD @@ -337,7 +338,7 @@ def calculate_patch_features(coords, areas, residues, lipo, ep, settings, hydro_ "SurfEpNeg": kept_patches(label_patches(indices, local_ep < -settings.ep_threshold), areas, settings.min_area), } - # Column order is the order of src/prodes/data/features_patches.yaml, which + # Column order is the order of src/prodes/data/features_reduced.yaml, which # documents each one; tests/test_patches.py checks the two agree. def block(in_set): return { diff --git a/src/prodes/data/features_patches.yaml b/src/prodes/data/features_patches.yaml deleted file mode 100644 index 35854e5..0000000 --- a/src/prodes/data/features_patches.yaml +++ /dev/null @@ -1,99 +0,0 @@ -# This file is the single source of truth for which patch features Prodes writes; -# prodes.feature_dictionary reads its feature codes and order from here. -# Maintained by hand: edit this file to add, drop or re-describe a patch feature. -meta: - title: Prodes surface patch features - description: The features written when the patch step is switched on, in output order. With a residue set NAME, every code is written a second time with _NAME appended, restricted to that set, preceded by SurfArea_NAME. - feature_count: 11 - written_when: --patches / patches=PatchSettings(...) - method: docs/surface_patches.md - default_settings: - radius_A: 5.0 - min_area_A2: 100.0 - hydrophobic_threshold: 0.0 - ep_threshold_volts: 0.5 - conventions: - local_value: The local value of a surface point is the area-weighted mean of the per-point value over every surface point within the radius, the point itself included. - patch: A patch is a set of surface points whose local value passes the threshold and which touch on the 1 A surface lattice. Patches smaller than the minimum area are discarded, judged on the whole patch. - area: Area is solvent accessible surface area with the 1.4 A probe, in A^2, which is larger than a molecular surface area of the same region. - residue_set: Patches are found on the whole surface. A restricted column counts only the area of patch points owned by a residue in the set. NSurf*Patch_NAME counts the kept patches with at least one point in the set; *PatchMaxArea_NAME is the largest in-set share of any one patch. - fields: - plot_name: Short label for an axis or legend. - short_explanation: One-line description of what the feature is. - long_explanation: How the feature is calculated and how to read it. - original_explanation: Always null; these features were added after Neijenhuis et al. 2024. - unit: Unit of measurement. Null when the feature is a count or dimensionless. -features: - SurfMhpPatchArea: - plot_name: Hydrophobic patch area (A^2) - short_explanation: Summed area of the hydrophobic surface patches. - long_explanation: Surface points whose local projected hydrophobicity is above the hydrophobic threshold are linked into patches where they touch; patches under the minimum area are discarded and the area of the rest is summed. Measures how much of the surface is covered by contiguous hydrophobic regions rather than by scattered hydrophobic points. - unit: Ų - original_explanation: null - SurfMhpPatchMaxArea: - plot_name: Largest hydrophobic patch (A^2) - short_explanation: Area of the largest hydrophobic surface patch. - long_explanation: The area of the single largest patch counted in SurfMhpPatchArea. 0 when no patch reaches the minimum area. - unit: Ų - original_explanation: null - NSurfMhpPatch: - plot_name: Number of hydrophobic patches - short_explanation: Number of hydrophobic surface patches at or above the minimum area. - long_explanation: The number of patches counted in SurfMhpPatchArea. - unit: null - original_explanation: null - SurfMhpLocalMax: - plot_name: Max local MHP - short_explanation: Highest local projected hydrophobicity on the surface. - long_explanation: The maximum over surface points of the projected hydrophobicity averaged over the neighbourhood radius, in the units of the hydrophobicity scale. The per-point value is dominated by the nearest one or two atoms, so the average is what gives it the length scale of a surface region; this is the most hydrophobic such region. - unit: null - original_explanation: null - SurfMhpLocalPosSum: - plot_name: Summed positive local MHP - short_explanation: Area-weighted sum of the positive part of the local projected hydrophobicity. - long_explanation: The sum over surface points of the point area times the local projected hydrophobicity where that is positive, zero elsewhere. Combines how hydrophobic and how large the hydrophobic regions are, without a threshold or a minimum patch area. The unit is square angstrom times the units of the hydrophobicity scale, so values from different scales are not comparable. On a scale with no negative values every point contributes. - unit: Ų - original_explanation: null - SurfEpPosPatchArea: - plot_name: Positive patch area (A^2) - short_explanation: Summed area of the positive electrostatic surface patches. - long_explanation: As SurfMhpPatchArea, for points whose local formal-charge electrostatic potential is above the electrostatic threshold. - unit: Ų - original_explanation: null - SurfEpPosPatchMaxArea: - plot_name: Largest positive patch (A^2) - short_explanation: Area of the largest positive electrostatic surface patch. - long_explanation: The area of the single largest patch counted in SurfEpPosPatchArea. - unit: Ų - original_explanation: null - NSurfEpPosPatch: - plot_name: Number of positive patches - short_explanation: Number of positive electrostatic surface patches at or above the minimum area. - long_explanation: The number of patches counted in SurfEpPosPatchArea. - unit: null - original_explanation: null - SurfEpNegPatchArea: - plot_name: Negative patch area (A^2) - short_explanation: Summed area of the negative electrostatic surface patches. - long_explanation: As SurfMhpPatchArea, for points whose local formal-charge electrostatic potential is below the negative of the electrostatic threshold. - unit: Ų - original_explanation: null - SurfEpNegPatchMaxArea: - plot_name: Largest negative patch (A^2) - short_explanation: Area of the largest negative electrostatic surface patch. - long_explanation: The area of the single largest patch counted in SurfEpNegPatchArea. - unit: Ų - original_explanation: null - NSurfEpNegPatch: - plot_name: Number of negative patches - short_explanation: Number of negative electrostatic surface patches at or above the minimum area. - long_explanation: The number of patches counted in SurfEpNegPatchArea. - unit: null - original_explanation: null -residue_set_only_features: - SurfArea: - plot_name: Surface area of the set (A^2) - short_explanation: Solvent accessible surface area owned by the residues of the set. - long_explanation: The summed area of the surface points whose residue is in the set. Written only as SurfArea_NAME, first among the columns of each set, so that the restricted patch areas can be normalised by the size of the region. - unit: Ų - original_explanation: null diff --git a/src/prodes/data/features_reduced.yaml b/src/prodes/data/features_reduced.yaml index 80d2cf9..13bc923 100644 --- a/src/prodes/data/features_reduced.yaml +++ b/src/prodes/data/features_reduced.yaml @@ -4,23 +4,32 @@ meta: title: Prodes reduced feature set description: The features Prodes calculates by default. Keys are the feature codes written as CSV column headings, in output order. - feature_count: 54 + feature_count: 65 written_when: always source: prodes (Neijenhuis et al. 2024), https://github.com/Tneijenhuis/prodes default_settings: ph: 7 probe_radius_A: 1.4 hydrophobicity_scale: mj_scaled (Miyazawa-Jernigan, scaled) + patch_radius_A: 5.0 + patch_min_area_A2: 100.0 + patch_hydrophobic_threshold: 0.0 + patch_ep_threshold_volts: 0.5 + patch_method: docs/surface_patches.md conventions: stat_suite: For any property X, the suite is Mean, Trimean, Median, Sum, Std. Trimean = (Q1 + 2*Q2 + Q3)/4. 'Max' and 'Min' are also reported separately. 'NPos*' / 'NNeg*' are the counts of surface points with strictly positive / negative values, and the 'Pos' / 'Neg' stats are computed over those subsets only. formal_vs_average: Formal = electrostatic potential computed using integer formal charges of titratable residues at the given pH (each residue is fully protonated or fully deprotonated). Average = computed using partial/pKa-weighted (Henderson-Hasselbalch) charges, so each titratable site contributes a fractional charge. units: Lengths in angstrom (A), area in A^2, mass in Da, charge in elementary units e. Electrostatic potentials are in arbitrary Coulombic units (charge_simple with eps_r=4); use only as relative descriptors. MHP is unitless. + local_value: For the patch features, the local value of a surface point is the area-weighted mean of the per-point value over every surface point within the patch radius, the point itself included. + patch: A patch is a set of surface points whose local value passes the threshold and which touch on the 1 A surface lattice. Patches smaller than the minimum area are discarded, judged on the whole patch. Patch areas are solvent accessible surface areas with the 1.4 A probe, which are larger than molecular surface areas of the same region. + residue_set: With a residue set NAME, every feature marked per_residue_set is written a second time as _NAME, restricted to that set, preceded by SurfArea_NAME. Patches are found on the whole surface; a restricted column counts only the area of patch points owned by a residue in the set. NSurf*Patch_NAME counts the kept patches with at least one point in the set; *PatchMaxArea_NAME is the largest in-set share of any one patch. fields: plot_name: Short label for an axis or legend. short_explanation: One-line description of what the feature is. long_explanation: How the feature is calculated and how to read it. original_explanation: The description from Supplemental Table 1 of Neijenhuis et al. 2024. Null for the few features added to the code after publication. unit: Unit of measurement as given by that paper. Null when the feature is dimensionless, such as a fraction or a ratio, or when the paper does not list it. + per_residue_set: True for the patch features, which are repeated for every residue set the run is given. Absent otherwise. features: Molecular weight: plot_name: Molecular weight (Da) @@ -346,3 +355,87 @@ features: long_explanation: Spread of negative shell-EP values across directions. unit: v original_explanation: The standard deviation of the negative electrostatic potentials calculated using binary charges of +1, 0 or -1 + SurfMhpPatchArea: + plot_name: Hydrophobic patch area (A^2) + short_explanation: Summed area of the hydrophobic surface patches. + long_explanation: Surface points whose local projected hydrophobicity is above the hydrophobic threshold are linked into patches where they touch; patches under the minimum area are discarded and the area of the rest is summed. Measures how much of the surface is covered by contiguous hydrophobic regions rather than by scattered hydrophobic points. + unit: Ų + original_explanation: null + per_residue_set: true + SurfMhpPatchMaxArea: + plot_name: Largest hydrophobic patch (A^2) + short_explanation: Area of the largest hydrophobic surface patch. + long_explanation: The area of the single largest patch counted in SurfMhpPatchArea. 0 when no patch reaches the minimum area. + unit: Ų + original_explanation: null + per_residue_set: true + NSurfMhpPatch: + plot_name: Number of hydrophobic patches + short_explanation: Number of hydrophobic surface patches at or above the minimum area. + long_explanation: The number of patches counted in SurfMhpPatchArea. + unit: null + original_explanation: null + per_residue_set: true + SurfMhpLocalMax: + plot_name: Max local MHP + short_explanation: Highest local projected hydrophobicity on the surface. + long_explanation: The maximum over surface points of the projected hydrophobicity averaged over the neighbourhood radius, in the units of the hydrophobicity scale. The per-point value is dominated by the nearest one or two atoms, so the average is what gives it the length scale of a surface region; this is the most hydrophobic such region. + unit: null + original_explanation: null + per_residue_set: true + SurfMhpLocalPosSum: + plot_name: Summed positive local MHP + short_explanation: Area-weighted sum of the positive part of the local projected hydrophobicity. + long_explanation: The sum over surface points of the point area times the local projected hydrophobicity where that is positive, zero elsewhere. Combines how hydrophobic and how large the hydrophobic regions are, without a threshold or a minimum patch area. The unit is square angstrom times the units of the hydrophobicity scale, so values from different scales are not comparable. On a scale with no negative values every point contributes. + unit: Ų + original_explanation: null + per_residue_set: true + SurfEpPosPatchArea: + plot_name: Positive patch area (A^2) + short_explanation: Summed area of the positive electrostatic surface patches. + long_explanation: As SurfMhpPatchArea, for points whose local formal-charge electrostatic potential is above the electrostatic threshold. + unit: Ų + original_explanation: null + per_residue_set: true + SurfEpPosPatchMaxArea: + plot_name: Largest positive patch (A^2) + short_explanation: Area of the largest positive electrostatic surface patch. + long_explanation: The area of the single largest patch counted in SurfEpPosPatchArea. + unit: Ų + original_explanation: null + per_residue_set: true + NSurfEpPosPatch: + plot_name: Number of positive patches + short_explanation: Number of positive electrostatic surface patches at or above the minimum area. + long_explanation: The number of patches counted in SurfEpPosPatchArea. + unit: null + original_explanation: null + per_residue_set: true + SurfEpNegPatchArea: + plot_name: Negative patch area (A^2) + short_explanation: Summed area of the negative electrostatic surface patches. + long_explanation: As SurfMhpPatchArea, for points whose local formal-charge electrostatic potential is below the negative of the electrostatic threshold. + unit: Ų + original_explanation: null + per_residue_set: true + SurfEpNegPatchMaxArea: + plot_name: Largest negative patch (A^2) + short_explanation: Area of the largest negative electrostatic surface patch. + long_explanation: The area of the single largest patch counted in SurfEpNegPatchArea. + unit: Ų + original_explanation: null + per_residue_set: true + NSurfEpNegPatch: + plot_name: Number of negative patches + short_explanation: Number of negative electrostatic surface patches at or above the minimum area. + long_explanation: The number of patches counted in SurfEpNegPatchArea. + unit: null + original_explanation: null + per_residue_set: true +residue_set_only_features: + SurfArea: + plot_name: Surface area of the set (A^2) + short_explanation: Solvent accessible surface area owned by the residues of the set. + long_explanation: The summed area of the surface points whose residue is in the set. Written only as SurfArea_NAME, first among the columns of each set, so that the restricted patch areas can be normalised by the size of the region. + unit: Ų + original_explanation: null diff --git a/src/prodes/feature_dictionary.py b/src/prodes/feature_dictionary.py index 0d09ff5..145c505 100644 --- a/src/prodes/feature_dictionary.py +++ b/src/prodes/feature_dictionary.py @@ -25,9 +25,10 @@ line up datasets calculated under either setting without re-running Prodes or diffing CSV headers by hand: - fd.get_feature_codes() # the 54 written by default - fd.get_feature_codes(full=True) # the 105 legacy columns + fd.get_feature_codes() # the 65 written by default + fd.get_feature_codes(full=True) # the 116: the 105 legacy columns and the 11 patch ones fd.get_dropped_feature_codes() # the 51 the default leaves out + fd.get_feature_codes(residue_sets=["CDR"]) # plus the columns of a residue set # cut a full-feature CSV down so it lines up with a reduced-feature one df = fd.drop_redundant_features(pd.read_csv("full_output.csv")) @@ -41,7 +42,7 @@ :data:`ID_COLUMN` when you need it. Everything, including which feature codes exist, comes from two human-readable -YAML files in ``prodes/data``: ``features_reduced.yaml`` for the 54 calculated by +YAML files in ``prodes/data``: ``features_reduced.yaml`` for the 65 calculated by default and ``features_full_only.yaml`` for the other 51, the latter carrying the reason each one was dropped. Those files are the single source of truth, so the list of features is never written down twice and cannot disagree with itself. Each @@ -55,11 +56,13 @@ Why the other 51 are dropped by default is set out in ``docs/redundant_feature_analysis.md``. -A third file, ``features_patches.yaml``, documents the surface patch features, -which are written only when the patch step is switched on (``--patches``). They -are in neither list above; ``get_patch_feature_codes`` returns them, and every -getter also accepts a patch code with a residue set suffix, such as -``SurfMhpPatchArea_CDR``. The method is in ``docs/surface_patches.md``. +The last 11 of the default features are the surface patch features, marked +``per_residue_set`` in ``features_reduced.yaml``. A run given a residue set NAME +writes each of them a second time as ``_NAME``, preceded by +``SurfArea_NAME``. ``get_patch_feature_codes`` returns the 11, every getter also +accepts a patch code with a residue set suffix, such as ``SurfMhpPatchArea_CDR``, +and ``get_feature_codes(residue_sets=...)`` includes the suffixed columns. The +method is in ``docs/surface_patches.md``. """ import importlib.resources @@ -71,7 +74,6 @@ REDUCED_DICTIONARY_FILE = "features_reduced.yaml" FULL_ONLY_DICTIONARY_FILE = "features_full_only.yaml" -PATCH_DICTIONARY_FILE = "features_patches.yaml" def _read_features(filename: str, key: str = "features") -> dict: @@ -101,13 +103,13 @@ def _read_features(filename: str, key: str = "features") -> dict: #: need that. FULL_FEATURE_CODES = REDUCED_FEATURE_CODES + DROPPED_FEATURE_CODES -# The patch features are a separate block, written only when the patch step is -# switched on, so they are kept out of the reduced and full sets above. -_PATCH_ENTRIES = _read_features(PATCH_DICTIONARY_FILE) -_RESIDUE_SET_ONLY_ENTRIES = _read_features(PATCH_DICTIONARY_FILE, "residue_set_only_features") +# The patch features are part of the reduced set; they are the entries marked +# per_residue_set, which a run repeats for every residue set it is given. +_PATCH_ENTRIES = {code: entry for code, entry in _REDUCED_ENTRIES.items() if entry.get("per_residue_set")} +_RESIDUE_SET_ONLY_ENTRIES = _read_features(REDUCED_DICTIONARY_FILE, "residue_set_only_features") -#: The whole-surface patch features, in the order Prodes writes them. With a -#: residue set NAME, each is written again as ``_NAME``. +#: The surface patch features, in the order Prodes writes them. With a residue +#: set NAME, each is written again as ``_NAME``. PATCH_FEATURE_CODES = tuple(_PATCH_ENTRIES) @@ -124,9 +126,6 @@ def _patch_entry(feature_code: str) -> dict | None: to ``Area``. """ - if feature_code in _PATCH_ENTRIES: - return _PATCH_ENTRIES[feature_code] - code, separator, name = feature_code.rpartition("_") # The rule PatchSettings applies to a set name: ASCII letters and digits. if separator and name.isascii() and name.isalnum(): @@ -161,37 +160,35 @@ def __contains__(self, feature_code: str) -> bool: # --- the feature code lists ------------------------------------------------- - def get_feature_codes(self, full: bool = False) -> list[str]: + def get_feature_codes(self, full: bool = False, residue_sets=()) -> list[str]: """Returns the feature codes, excluding ``ID``. The reduced set comes back in the order Prodes writes those columns. The full set is the reduced codes followed by the full-only ones, which is not the column order of a ``--full-features`` run: that interleaves the two. - Read the CSV header if you need the exact column order. + Read the CSV header if you need the exact column order. The columns of + the residue sets come last, in the order Prodes writes them. Args: - full: return the full legacy set instead of the reduced default. + full: return the full set instead of the reduced default. + residue_sets: the names of the run's residue sets, in the order + they were given. Returns: A new list of feature codes. """ - return list(FULL_FEATURE_CODES if full else REDUCED_FEATURE_CODES) - - def get_patch_feature_codes(self, residue_sets=()) -> list[str]: - """Returns the patch feature codes a run with the patch step writes, in output order. - - Args: - residue_sets: the names of the residue sets of the run, in the order - they were given. - """ - - codes = list(PATCH_FEATURE_CODES) + codes = list(FULL_FEATURE_CODES if full else REDUCED_FEATURE_CODES) for name in residue_sets: codes.extend(residue_set_feature_codes(name)) return codes + def get_patch_feature_codes(self) -> list[str]: + """Returns the 11 surface patch feature codes, which a run repeats for every residue set.""" + + return list(PATCH_FEATURE_CODES) + def get_dropped_feature_codes(self) -> list[str]: """Returns the feature codes in the full set that the reduced set leaves out.""" @@ -274,13 +271,15 @@ def get_entry(self, feature_code: str) -> dict: raise KeyError( f"{feature_code!r} is not a Prodes feature code. Use get_feature_codes(full=True) for the {len(FULL_FEATURE_CODES)} " - "valid codes, and get_patch_feature_codes() for the patch features." + "valid codes; a patch code may also carry a residue set suffix, such as SurfMhpPatchArea_CDR." ) def get_dictionary(self) -> dict: """Returns every entry of the reduced and full sets, keyed by feature code. - The patch features are not included; get_entry looks them up. + The patch features are included under their plain codes; the residue + set columns, such as ``SurfMhpPatchArea_CDR``, are not, but get_entry + looks them up. Returns: A new dict of new dicts, safe to mutate. diff --git a/src/prodes/output.py b/src/prodes/output.py index 17debc7..db45d50 100644 --- a/src/prodes/output.py +++ b/src/prodes/output.py @@ -76,10 +76,10 @@ {name}_features.csv the calculated features, one row {name}_surface_points.csv every surface point: x, y, z, potential, hydrophobicity, the area it stands for, and - the residue owning most of that area; with - --patches also the neighbourhood averages - and the hydrophobic patch it belongs to - (-1 for none) + the residue owning most of that area, the + neighbourhood averages of hydrophobicity and + potential, and the hydrophobic patch it + belongs to (-1 for none) {name}_ep.pml PyMOL, coloured by electrostatic potential {name}_hydrophobicity.pml PyMOL, coloured by hydrophobicity {name}_ep.cxc ChimeraX, electrostatic potential diff --git a/src/prodes/run.py b/src/prodes/run.py index 326c9a1..3015e2f 100644 --- a/src/prodes/run.py +++ b/src/prodes/run.py @@ -69,31 +69,28 @@ def parse_arguments(): default=DEFAULT_IONIC_STRENGTH_MOLAR, ) parser.add_argument("--hydro", help="Abbreviation of the hydrophobicity scale to be used", type=str, default="mj_scaled") - parser.add_argument( - "--patches", - action="store_true", - help="Also report contiguous hydrophobic, positive and negative surface patches; see docs/surface_patches.md", - ) parser.add_argument( "--residue-set", action="append", default=[], metavar="NAME=FILE", - help="With --patches, also report every patch feature restricted to the residues listed in FILE, one label " - "such as H100A per line, under columns ending in _NAME. Repeatable", + help="Also report every surface patch feature restricted to the residues listed in FILE, one label " + "such as H100A per line, under columns ending in _NAME. Repeatable; see docs/surface_patches.md", + ) + parser.add_argument( + "--patch-radius", type=float, default=DEFAULT_PATCH_RADIUS, help=f"Neighbourhood radius in A for the patch average (default: {DEFAULT_PATCH_RADIUS})" ) - parser.add_argument("--patch-radius", type=float, default=None, help=f"Neighbourhood radius in A for the patch average (default: {DEFAULT_PATCH_RADIUS})") - parser.add_argument("--patch-min-area", type=float, default=None, help=f"Smallest patch kept, in A^2 (default: {DEFAULT_PATCH_MIN_AREA})") + parser.add_argument("--patch-min-area", type=float, default=DEFAULT_PATCH_MIN_AREA, help=f"Smallest patch kept, in A^2 (default: {DEFAULT_PATCH_MIN_AREA})") parser.add_argument( "--patch-hydrophobic-threshold", type=float, default=None, - help="Neighbourhood-averaged hydrophobicity above which a point is hydrophobic (default: 0, refused for a scale with no negative values)", + help="Neighbourhood-averaged hydrophobicity above which a point is hydrophobic (default: 0; required with a scale that has no negative values, such as mj)", ) parser.add_argument( "--patch-ep-threshold", type=float, - default=None, + default=DEFAULT_EP_THRESHOLD, help=f"Neighbourhood-averaged potential in volts beyond which a point is positive or negative (default: {DEFAULT_EP_THRESHOLD})", ) parser.add_argument( @@ -149,26 +146,13 @@ def parse_arguments(): if not math.isfinite(ionic_strength_molar) or ionic_strength_molar < 0: raise ValueError(f"--ionic-strength must be a finite number and not negative, got {ionic_strength_molar}") - patch_options = { - "--patch-radius": arg.patch_radius, - "--patch-min-area": arg.patch_min_area, - "--patch-hydrophobic-threshold": arg.patch_hydrophobic_threshold, - "--patch-ep-threshold": arg.patch_ep_threshold, - } - patches = None - if arg.patches: - patches = PatchSettings( - radius=DEFAULT_PATCH_RADIUS if arg.patch_radius is None else arg.patch_radius, - min_area=DEFAULT_PATCH_MIN_AREA if arg.patch_min_area is None else arg.patch_min_area, - hydrophobic_threshold=arg.patch_hydrophobic_threshold, - ep_threshold=DEFAULT_EP_THRESHOLD if arg.patch_ep_threshold is None else arg.patch_ep_threshold, - residue_sets=read_residue_sets(arg.residue_set), - ) - else: - # Without --patches these would be silently ignored. - given = [option for option, value in patch_options.items() if value is not None] + (["--residue-set"] if arg.residue_set else []) - if given: - raise ValueError(f"{', '.join(given)} only apply to the patch features, so they need --patches") + patches = PatchSettings( + radius=arg.patch_radius, + min_area=arg.patch_min_area, + hydrophobic_threshold=arg.patch_hydrophobic_threshold, + ep_threshold=arg.patch_ep_threshold, + 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 @@ -592,10 +576,10 @@ def calculate( run, divided among the workers. If None, falls back to the PRODES_MEM_LIMIT_MB env var (default 2048). Pass explicitly when calling as a library to avoid relying on environment variables. - patches: a prodes.calculations.patches.PatchSettings to also report - contiguous surface patches, optionally restricted to named residue - sets; None, the default, leaves them out and the feature table is - the same as without this option. + patches: a prodes.calculations.patches.PatchSettings for the surface + 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. Thread safety Safe to run in several processes, one structure each, which is how the @@ -607,11 +591,15 @@ def calculate( validate_settings(mem_limit_mb) - structure = prepare_structure(pdb_file, pkas_file) + if patches is None: + patches = PatchSettings() + # Checked before the structure is read, since the threshold depends only on + # the scale; the residue labels are checked as soon as it has been read. # Both checks are cheap, and the surface phase that follows is not. - if patches is not None: - patches.resolved_hydrophobic_threshold(hydro_scale) - check_residue_sets(patches.residue_sets, [residue.label for residue in structure.residues]) + patches.resolved_hydrophobic_threshold(hydro_scale) + + structure = prepare_structure(pdb_file, pkas_file) + 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. print(f"{structure.name}: {len(structure.disulfides)} disulfide bonds detected") @@ -634,20 +622,18 @@ def calculate( # Computed from the snapshot above, so the average charge phase overwriting # the potential on the points cannot reach it, and appended after the shell - # features so that the default columns keep their positions. - patch_features: dict = {} - if patches is not None: - patch_features, patch_columns = calculate_patch_features( - coords, - areas, - residues, - lipo, - ep, - patches, - hydro_scale, - [residue.label for residue in structure.residues], - ) - point_columns.update(patch_columns) + # features, which is their place in features_reduced.yaml. + patch_features, patch_columns = calculate_patch_features( + coords, + areas, + residues, + lipo, + ep, + patches, + hydro_scale, + [residue.label for residue in structure.residues], + ) + point_columns.update(patch_columns) if full_features: features = calculate_average_chargesurface_grid_features(structure, surface_points, ph, features, ionic_strength_molar=ionic_strength_molar) @@ -661,9 +647,8 @@ def calculate( "hydro_scale": hydro_scale, "full_features": full_features, "ionic_strength_molar": ionic_strength_molar, + "patches": patches.record(hydro_scale), } - if patches is not None: - settings["patches"] = patches.record(hydro_scale) metadata = run_metadata( pdb_file, settings, diff --git a/tests/test_cli.py b/tests/test_cli.py index 7ab1d17..72a3b8b 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -34,7 +34,7 @@ def test_defaults_are_the_documented_ones(monkeypatch): assert full_features is False assert mem_limit_mb is None assert ionic_strength_molar == DEFAULT_IONIC_STRENGTH_MOLAR - assert patches is None + assert patches == PatchSettings() def test_n_workers_flag_sets_the_environment_variable(monkeypatch): @@ -110,7 +110,7 @@ def fake_calculate(*args, **kwargs): True, None, DEFAULT_IONIC_STRENGTH_MOLAR, - None, + PatchSettings(), ) assert captured["kwargs"] == {} @@ -132,14 +132,13 @@ def test_calculate_signature_matches_the_parsed_order(): ] -def test_patches_flag_builds_the_settings(monkeypatch, tmp_path): - """--patches and its options reach calculate as one PatchSettings, residue sets read from file.""" +def test_patch_options_build_the_settings(monkeypatch, tmp_path): + """The patch options reach calculate as one PatchSettings, residue sets read from file.""" labels = tmp_path / "cdr.txt" labels.write_text("# CDR H3\nH95\nH100A # an insertion\n\nL91\n") *_, patches = parse( monkeypatch, - "--patches", "--patch-radius", "4", "--patch-min-area", @@ -155,35 +154,16 @@ def test_patches_flag_builds_the_settings(monkeypatch, tmp_path): assert patches == PatchSettings(radius=4, min_area=30, hydrophobic_threshold=0.02, ep_threshold=0.3, residue_sets={"CDR": ["H95", "H100A", "L91"]}) -def test_residue_set_without_patches_is_refused(monkeypatch, tmp_path): - """A residue set restricts the patch features, so on its own it would silently do nothing.""" - - labels = tmp_path / "cdr.txt" - labels.write_text("H95\n") - - with pytest.raises(ValueError, match="need --patches"): - parse(monkeypatch, "--residue-set", f"CDR={labels}") - - @pytest.mark.parametrize("specification", ["CDR", "=cdr.txt", "CDR="]) def test_a_malformed_residue_set_is_refused(monkeypatch, specification): """NAME=FILE is the only accepted form.""" with pytest.raises(ValueError, match="NAME=FILE"): - parse(monkeypatch, "--patches", "--residue-set", specification) + parse(monkeypatch, "--residue-set", specification) -@pytest.mark.parametrize( - "option", [["--patch-radius", "3"], ["--patch-min-area", "80"], ["--patch-hydrophobic-threshold", "0.1"], ["--patch-ep-threshold", "0.2"]] -) -def test_a_patch_option_without_patches_is_refused(monkeypatch, option): - """Otherwise the option would be accepted and silently do nothing.""" - - with pytest.raises(ValueError, match=f"{option[0]} only apply to the patch features"): - parse(monkeypatch, *option) +def test_the_patches_option_no_longer_exists(monkeypatch): + """The patch features are part of every run, so there is nothing to switch on.""" - -def test_patches_alone_gives_the_documented_defaults(monkeypatch): - *_, patches = parse(monkeypatch, "--patches") - - assert patches == PatchSettings() + with pytest.raises(SystemExit): + parse(monkeypatch, "--patches") diff --git a/tests/test_feature_dictionary.py b/tests/test_feature_dictionary.py index 882e55a..2848a96 100644 --- a/tests/test_feature_dictionary.py +++ b/tests/test_feature_dictionary.py @@ -23,6 +23,7 @@ FULL_FEATURE_CODES, FULL_ONLY_DICTIONARY_FILE, ID_COLUMN, + PATCH_FEATURE_CODES, REDUCED_DICTIONARY_FILE, REDUCED_FEATURE_CODES, FeatureDictionary, @@ -39,7 +40,7 @@ def test_get_the_reduced_feature_codes(): fd = FeatureDictionary() codes = fd.get_feature_codes() - assert len(codes) == 54 + assert len(codes) == 65 assert codes[0] == "Molecular weight" assert "SurfEpMeanFormal" in codes assert "SurfEpMeanAverage" not in codes @@ -50,7 +51,7 @@ def test_get_the_full_feature_codes(): fd = FeatureDictionary() codes = fd.get_feature_codes(full=True) - assert len(codes) == 105 + assert len(codes) == 116 assert "SurfEpMeanAverage" in codes @@ -91,7 +92,7 @@ def test_label_a_whole_dataframe_of_results(): fd = FeatureDictionary() labels = {code: fd.get_plot_name(code) for code in fd.get_feature_codes()} - assert len(labels) == 54 + assert len(labels) == 65 assert all(labels.values()) @@ -158,9 +159,15 @@ def test_full_codes_match_the_original_reference_output(): reference = {column for column in pd.read_csv(ORIG_OUTPUT).columns if column != ID_COLUMN} ours = set(FULL_FEATURE_CODES) - added = sorted(ours - reference) + # The 11 surface patch features are the one deliberate addition. They were + # put in the default set in version 10.0, after a model of antibody HIC + # retention improved with them on a held-out panel; see + # docs/surface_patches.md. Their redundancy against the others has not been + # analysed the way docs/redundant_feature_analysis.md analysed the 105. + added = sorted(ours - reference - set(PATCH_FEATURE_CODES)) missing = sorted(reference - ours) + assert len(PATCH_FEATURE_CODES) == 11 assert not added, f"in the YAML dictionaries but not in the original Prodes output: {added}. If this is a genuinely new feature, see this test's docstring." assert ( not missing @@ -379,11 +386,12 @@ def test_original_explanation_is_none_for_features_the_paper_omits(): assert fd.get_original_explanation("SurfNegMhpMedian") is None -def test_the_paper_documents_all_but_two_features(): +def test_the_paper_documents_every_original_feature_but_two(): + """The patch features postdate the paper, so none of them has its wording.""" fd = FeatureDictionary() undocumented = [code for code in FULL_FEATURE_CODES if fd.get_original_explanation(code) is None] - assert undocumented == ["SurfPosMhpMedian", "SurfNegMhpMedian"] + assert set(undocumented) == {"SurfPosMhpMedian", "SurfNegMhpMedian", *PATCH_FEATURE_CODES} def test_our_description_and_the_paper_agree_on_the_statistic(): diff --git a/tests/test_output_bundle.py b/tests/test_output_bundle.py index 5d763c1..abb4c7f 100644 --- a/tests/test_output_bundle.py +++ b/tests/test_output_bundle.py @@ -80,7 +80,7 @@ def test_the_tables_read_back_out_of_the_archive(bundle_zip): assert len(features) == 1 assert features["ID"].iloc[0] == "ARH96693" - assert list(points.columns) == ["x", "y", "z", "ep_volts", "hydrophobicity", "area_A2", "residue"] + assert list(points.columns) == ["x", "y", "z", "ep_volts", "hydrophobicity", "area_A2", "residue", "local_mhp", "local_ep_volts", "mhp_patch_id"] assert len(points) > 0 @@ -106,6 +106,7 @@ def test_the_run_record_says_what_was_run(bundle_dir): "hydro_scale": "mj_scaled", "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": {}}, } # ARH96693 is an AlphaFold model with three cysteines and no disulfide, so # the key has to be present and zero rather than absent. diff --git a/tests/test_patches.py b/tests/test_patches.py index 0023c82..187af88 100644 --- a/tests/test_patches.py +++ b/tests/test_patches.py @@ -29,7 +29,7 @@ patch_statistics, residue_set_masks, ) -from prodes.feature_dictionary import REDUCED_FEATURE_CODES, FeatureDictionary +from prodes.feature_dictionary import PATCH_FEATURE_CODES, REDUCED_FEATURE_CODES, FeatureDictionary from prodes.output import read_features, read_surface_points, write_bundle from prodes.run import calculate, construct_surface_grid, prepare_structure @@ -261,21 +261,21 @@ def test_a_tie_goes_to_the_first_label_whatever_the_order(): @cache -def _run(full_features, patched): +def _run(full_features, with_set): with tempfile.TemporaryDirectory(prefix="prodes_patch_tests_") as directory: out_file = Path(directory) / "bundle.zip" - patches = PatchSettings(residue_sets={"Nterm": N_TERMINUS}) if patched else None + patches = PatchSettings(residue_sets={"Nterm": N_TERMINUS}) if with_set else None calculate(STRUCTURE, str(out_file), full_features=full_features, patches=patches) return read_features(out_file), read_surface_points(out_file) -def run(full_features=False, patched=False): +def run(full_features=False, with_set=False): """Returns copies of the feature and point tables of one run per setting. - With patched, the patch step runs with one residue set, Nterm. + With with_set, the run is given one residue set, Nterm. """ - features, points = _run(full_features, patched) + features, points = _run(full_features, with_set) return features.copy(), points.copy() @@ -289,35 +289,35 @@ def test_point_areas_add_up_to_the_surface_area(): assert all(point.residue for point in points) -def test_the_default_features_are_unchanged_by_the_patch_step(): - """The anion exchange mode depends on these values; the patch step only appends.""" +def test_a_residue_set_only_appends_columns(): + """Patches are found on the whole surface, so a residue set cannot change a whole-surface value.""" - default, _ = run() - patched, _ = run(patched=True) + plain, _ = run() + restricted, _ = run(with_set=True) - assert list(default.columns) == ["ID", *REDUCED_FEATURE_CODES] - assert list(patched.columns[: len(default.columns)]) == list(default.columns) - assert patched[default.columns].equals(default) + assert list(plain.columns) == ["ID", *REDUCED_FEATURE_CODES] + assert list(restricted.columns[: len(plain.columns)]) == list(plain.columns) + assert restricted[plain.columns].equals(plain) def test_the_columns_are_the_ones_the_dictionary_documents(): - patched, _ = run(patched=True) + restricted, _ = run(with_set=True) - assert list(patched.columns[1 + len(REDUCED_FEATURE_CODES) :]) == FeatureDictionary().get_patch_feature_codes(["Nterm"]) + assert list(restricted.columns[1:]) == FeatureDictionary().get_feature_codes(residue_sets=["Nterm"]) def test_the_patch_step_reads_the_formal_potential_under_full_features(): """The partial-charge phase overwrites the potential on the points; the patch step must run before it.""" - reduced, _ = run(patched=True) - full, _ = run(full_features=True, patched=True) - codes = FeatureDictionary().get_patch_feature_codes(["Nterm"]) + reduced, _ = run(with_set=True) + full, _ = run(full_features=True, with_set=True) + codes = [code for code in FeatureDictionary().get_feature_codes(residue_sets=["Nterm"]) if code in PATCH_FEATURE_CODES or code.endswith("_Nterm")] assert full[codes].equals(reduced[codes]) def test_the_residue_set_selects_the_points_it_names(): - features, points = run(patched=True) + features, points = run(with_set=True) in_set = points["residue"].isin(N_TERMINUS) assert features["SurfArea_Nterm"].iloc[0] == pytest.approx(points.loc[in_set, "area_A2"].sum(), abs=1e-3) @@ -326,7 +326,7 @@ def test_the_residue_set_selects_the_points_it_names(): def test_the_point_table_agrees_with_the_hydrophobic_patch_area(): - features, points = run(patched=True) + features, points = run(with_set=True) in_patch = points["mhp_patch_id"] >= 0 assert features["SurfMhpPatchArea"].iloc[0] == pytest.approx(points.loc[in_patch, "area_A2"].sum(), abs=1e-3) @@ -334,7 +334,7 @@ def test_the_point_table_agrees_with_the_hydrophobic_patch_area(): def test_restricted_values_do_not_exceed_whole_surface_ones(): - features = run(patched=True)[0].iloc[0] + features = run(with_set=True)[0].iloc[0] for code in ("SurfMhpPatchArea", "SurfEpPosPatchArea", "SurfEpNegPatchArea", "SurfMhpLocalPosSum"): assert features[f"{code}_Nterm"] <= features[code] @@ -413,12 +413,13 @@ def test_the_dictionary_resolves_a_residue_set_suffix_on_patch_codes_only(): assert "Area_CDR" not in fd assert "SurfMhpPatchArea_CDR_H3" not in fd assert "SurfMhpPatchArea_CDR²" not in fd - with pytest.raises(KeyError, match="get_patch_feature_codes"): + with pytest.raises(KeyError, match="residue set suffix"): fd.get_entry("Area_CDR") -def test_the_patch_codes_are_kept_out_of_the_reduced_and_full_sets(): +def test_the_patch_codes_are_the_last_of_the_default_set(): fd = FeatureDictionary() - assert set(fd.get_patch_feature_codes()).isdisjoint(fd.get_feature_codes(full=True)) - assert fd.get_patch_feature_codes(["CDR"])[11] == "SurfArea_CDR" + assert fd.get_feature_codes()[-11:] == fd.get_patch_feature_codes() == list(PATCH_FEATURE_CODES) + assert set(PATCH_FEATURE_CODES) <= set(fd.get_feature_codes(full=True)) + assert fd.get_feature_codes(residue_sets=["CDR"])[65:67] == ["SurfArea_CDR", "SurfMhpPatchArea_CDR"] diff --git a/tests/test_pka.py b/tests/test_pka.py index dc10e73..59086b1 100644 --- a/tests/test_pka.py +++ b/tests/test_pka.py @@ -11,9 +11,9 @@ 1GDW a real PROPKA prediction moves 19 of the 54 reduced features, formal charge included, so a refactor that quietly stopped applying the values would change more than a third of the output while every other test still passed. -Downstream code already depends on it: biochai runs PROPKA and imports -convert_propka to read the result, so this is a public API with an outside -caller. +Downstream code already depends on it: callers typically run PROPKA and +import convert_propka to read the result, so this is a public API with outside +callers. The chain has four links and one test each below:: diff --git a/tests/test_sasa.py b/tests/test_sasa.py index 812c3d8..deab5db 100644 --- a/tests/test_sasa.py +++ b/tests/test_sasa.py @@ -24,6 +24,7 @@ from prodes.calculations import sasa from prodes.calculations.geometry import Sunflower_sphere from prodes.calculations.grid_wizard import Grid +from prodes.feature_dictionary import PATCH_FEATURE_CODES from prodes.io import parser as ps from tests.pipeline_output import pipeline_output @@ -59,10 +60,13 @@ def test_all_columns_present(): - """Every column in the original output should appear in the calculated output.""" + """Every column in the original output should appear in the calculated output. + + The only columns the original did not write are the surface patch features. + """ calc_cols = set(FULL_OUTPUT.columns) orig_cols = set(ORIGINAL_OUTPUT.columns) - assert calc_cols == orig_cols, f"Column mismatch: missing={orig_cols - calc_cols}, extra={calc_cols - orig_cols}" + assert calc_cols - set(PATCH_FEATURE_CODES) == orig_cols, f"Column mismatch: missing={orig_cols - calc_cols}, extra={calc_cols - orig_cols}" def test_feature_values_match_original():