[teese]

Problem

Prodes computes a table of protein surface features from a structure, for use in machine learning models of chromatography. It gives each point of the protein surface two values, a projected hydrophobicity (the molecular hydrophobicity potential, Mhp in column names) and an electrostatic potential (Ep), and summarises each over the whole surface by its sign. Every hydrophobic point is pooled whether the points form one region or twenty scattered ones, and no feature can be restricted to part of the molecule, such as the CDRs of an antibody. Hydrophobic interaction chromatography retention of antibodies has been reported to correlate with the area of contiguous hydrophobic surface on the CDRs, which no Prodes feature measures.

Closes #35

What changed

• New step in every run, in src/prodes/calculations/patches.py. Its settings and residue sets are passed as CLI options or calculate(..., patches=PatchSettings(...)).
• Each surface point's hydrophobicity and potential are averaged over a 5 Å neighbourhood; points past a threshold that touch on the 1 Å surface grid form a patch; patches under 100 Ų are dropped.
• 11 new default columns, at the end of src/prodes/data/features_reduced.yaml: summed area, largest area and count of the hydrophobic patches (SurfMhpPatchArea, SurfMhpPatchMaxArea, NSurfMhpPatch), the same three for positive and for negative potential (SurfEpPos*, SurfEpNeg*), and the highest and the area-weighted positive sum of the averaged hydrophobicity (SurfMhpLocalMax, SurfMhpLocalPosSum).
• --residue-set NAME=FILE repeats every column as <column>_NAME, restricted to the residues listed in the file.
• The surface point table in the output zip carries each point's area (area_A2), owning residue (residue), neighbourhood averages (local_mhp, local_ep_volts) and hydrophobic patch (mhp_patch_id).
• Breaking, version 10.0.0: the default feature table has 65 columns instead of 54, --full-features writes 116 instead of 105, --patches is gone, and --hydro mj needs --patch-hydrophobic-threshold. No earlier column changes value.
• Method, settings and limitations in docs/surface_patches.md.
• The command line enables faulthandler, so a crash in compiled code prints a Python stack instead of exiting silently. docs/parallelism_and_memory.md explains the return codes of a run killed by a signal.
• docs/surface_patches.md records that the default patch settings were checked against antibody HIC data and kept: settings tuned on GDPa1 did worse on the held-out GDPa3 panel.
• docs/structure_preparation.md states that hydrogens do not change any feature value, and why PROPKA leaves some groups without a pKa (a C-terminus without OXT; an N-terminal aspartate).

Results

Measured outside this repository, on version 9.1.0 of this branch, which wrote the same patch columns under a --patches option, with a random-forest model of HIC retention built on sequence-derived features, to which the patch columns were added. Patches were computed with two residue sets per antibody: the CDRs, and the whole variable region (Fv), both from AHo numbering. Scores are Spearman R between predicted and measured retention on the public GDPa1 panel (about 236 antibodies), using the five cross-validation folds its authors publish, the same folds in every row; and on the public GDPa3 panel (79 antibodies), predicted by a model fitted on all of GDPa1. Each figure is the mean over five random-forest seeds.

Scored on	Structure source	Sequence model	Sequence model + patch columns	Improved in
GDPa1, mean per fold	ABodyBuilder3 Fv	0.533	0.573	5 of 5 folds
GDPa1, mean per fold	Boltz-2 whole antibody	0.529	0.589	5 of 5 folds
GDPa3, held out	ABodyBuilder3 Fv	0.527	0.620	5 of 5 seeds

Adding the 54 default Prodes columns to the sequence model instead changed the GDPa1 mean per-fold Spearman R by between -0.011 and -0.001, depending on the structure source. The patch columns changed anion exchange (pooled out-of-fold Spearman R 0.763 to 0.768) and cation exchange (0.822 to 0.814) by less than 0.01. On 423 structures (Fv models, whole antibodies, anion and cation exchange proteins) the 54 earlier default columns were identical with and without the patch step.

The patch step takes about 1 s for 57,000 surface points; the whole Prodes run on the largest test structure, PDB entry 1AO6, takes 67 s including it.

Checks on tests/data/ARH96693.pdb.zip, the single-chain structure the regression tests use:

Check	Result
A run with a residue set against one without: every whole-surface column	identical
Sum of the point areas against the Area feature (total solvent accessible surface area)	equal
Patch features with --full-features (the full 105-feature legacy set) against without	identical

CI not run yet. 488/488 tests pass locally.

Acceptance criteria still open

All acceptance criteria met.

Not verified

• The default settings (radius 5 Å, minimum patch area 100 Ų, hydrophobic threshold 0, electrostatic threshold 0.5 V) were checked on antibody HIC data only. Tuned on GDPa1, the best setting raised the correlation of SurfMhpPatchArea_CDR with retention there from 0.34 to 0.46 but lowered it on GDPa3 from 0.355 to 0.328, and the GDPa3 score of the HIC model fell from 0.615 to 0.589, so the defaults were kept. They have not been checked for other proteins or other chromatography modes.
• One of 26 runs of one whole antibody, at 12 workers, gave different ShellEp*Formal columns; every other column agreed. These columns summarise the electrostatic potential projected onto 120 planes placed around the protein, and in that run the summed potential of one plane changed sign. 25 further runs of the same input agreed: 8 complete runs, and 17 repeats of only the step that computes the 120 planes, at 1, 4 and 12 workers, which agreed to within 1e-15. The cause was not found.
• The 11 patch columns were not part of the redundancy analysis behind the 54 earlier default columns (docs/redundant_feature_analysis.md).
• The dense arrays of the neighbourhood average are not counted against --mem-limit: a few MB for an Fv, a few tens of MB for a whole antibody (estimated, not measured).

Detail

• Patches are found on the whole surface, then restricted: a patch crossing the border of a residue set is one patch, the minimum area applies to the whole patch, and the restricted columns count only its share inside the set.
• Each point's area is the summed area of the Shrake-Rupley dots in its 1 Å cell, weighted as Atom.surface_area weights them; its residue is the residue whose dots carry the most area in the cell, ties going to the alphabetically first label.
• The 0.5 V electrostatic default lies between the medians of the 75th (0.39 V) and 90th (0.61 V) percentiles of the absolute local potential over the eight structures in tests/data, at 0.15 M.
• Refused with an error: the default hydrophobic threshold of 0 on a hydrophobicity scale with no negative values (mj, the unscaled Miyazawa-Jernigan scale); a residue label the structure does not have, checked before the surface is built; a residue set that owns no surface point; a radius above 15 Å.
• FeatureDictionary resolves <code>_NAME for patch codes only, so Area_CDR is not a feature code.