#35: surface patch features in every run and in the default feature set (10.0.0, breaking)

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.

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.

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

docs/surface_patches.md

@@ -0,0 +1,107 @@
+# Surface patches
+
+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 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 --residue-set CDR=cdr_residues.txt --residue-set Fv=fv_residues.txt
+```
+
+Python, where `patches` is optional and only needed to name residue sets or change a setting:
+
+```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; 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. |
+
+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.
+
+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.
+
+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_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`, 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
+
+- **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.

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "prodes"
-version = "9.0.0"
+version = "10.0.0"
 description = "A package to calculate protein surface descriptors"
 readme = "README.rst"
 license = {text = "MIT"}

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 for the surface
+            patch features, such as residue sets to restrict them to. None means
+            the default settings and no residue sets.
     """
 
     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,
     )