From 9c40b288058a7cbe0cc448135c43129ef1dafdbb Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 00:38:55 +0200 Subject: [PATCH 1/8] Recommend repairing the structure with PDB2PQR before running PROPKA and Prodes Prodes describes the structure it is given and does not correct one. A residue whose side chain was not modelled carries no charge, and nothing in the output records that this happened, so every charge-derived feature is wrong without any signal to the user. PDB2PQR repairs the structure: it rebuilds missing heavy atoms, adds hydrogens, optimises the hydrogen bond network, resolves the resulting steric clashes, adds a missing OXT, and reports backbone chain breaks and gaps it cannot repair. It is BSD-3-Clause, pure Python, and available from conda-forge. Measured on PDB entry 4HKZ, a Fab with two lysines modelled only as far as CB: rebuilding those two side chains alone moves 31 of the 54 default features, takes the formal charge at pH 7.4 from -2 to 0, moves the isoelectric point from 6.46 to 7.61, and flips the sign of the mean surface electrostatic potential. For choosing between an anion and a cation exchanger at a given buffer pH, those are different recommendations. Changes: - A "Preparing your structure" section in the README covering why the step is needed, what it fixes, what Prodes keeps doing itself, and the traps. - The quick start and the worked QSPR loop now repair the structure first and run every later step, PROPKA included, on the repaired file. Running PROPKA before the repair loses the titratable groups that the repair restores: on 4HKZ, 145 groups against 149. - docs/structure_preparation.md holds the detail, including what Prodes discards from the repaired file and why. - pdb2pqr added to environment.yml, pinned, alongside propka and for the same reason. The repaired file is written into a subdirectory under its original name rather than being renamed, because the ID column of the output is taken from the file name. PROPKA writes its .pka file into the working directory rather than beside its input, which the documented commands account for. PDB2PQR can run PROPKA itself, but that does not replace the propka3 step: it writes no .pka file for the converter to read, and it fixes a single pH where Prodes deliberately separates prediction from calculation. --- README.rst | 128 +++++++++++++++++++++++++++++----- docs/structure_preparation.md | 107 ++++++++++++++++++++++++++++ environment.yml | 5 ++ 3 files changed, 224 insertions(+), 16 deletions(-) create mode 100644 docs/structure_preparation.md diff --git a/README.rst b/README.rst index 843ccba..a5c3568 100644 --- a/README.rst +++ b/README.rst @@ -7,7 +7,7 @@ Those columns are designed to be used directly as the input features of a machin Prodes is **completely free, including for commercial use**, under the MIT licence. It needs no external electrostatics solver, no licence server and no web upload. A protein under 1000 residues should only take a few minutes on a normal desktop computer, and seconds on a dedicated linux server. -**Jump to:** `Installation`_ | `Quick start`_ | `Viewing the surface`_ | `The output bundle`_ | `Ionic strength and screening`_ | `Using the features in a model`_ | `pKa values and protonation states`_ | `Using Prodes from Python`_ | `Speed`_ | `How to cite`_ +**Jump to:** `Installation`_ | `Quick start`_ | `Viewing the surface`_ | `The output bundle`_ | `Ionic strength and screening`_ | `Using the features in a model`_ | `Preparing your structure`_ | `pKa values and protonation states`_ | `Using Prodes from Python`_ | `Speed`_ | `How to cite`_ This is a fork of `tneijenhuis/prodes `_, a package written by `Tim Neijenhuis `_ during his Ph.D. at the `Marcel Ottens group `_ at the `Delft University of Technology (TU Delft) `_. Currently this fork preserves the original algorithm. The changes are performance (a 170x speedup for many proteins, see `Speed`_) and a reduced, non-redundant default feature set (see `The reduced feature set`_). @@ -45,11 +45,17 @@ Installing as a user conda create -n prodes python=3.13 conda activate prodes pip install git+https://github.com/datacatalysis/prodes.git - conda install conda-forge::propka + conda install conda-forge::propka conda-forge::pdb2pqr Prodes depends on NumPy, pandas and `Biopython `_, which ``pip`` installs with it. Biopython reads the coordinate records of a PDB file; the rules that decide which alternate conformation of a residue to keep, and which cysteines are bonded, remain prodes' own. See `how prodes reads a PDB file `_. -The fourth line adds `PROPKA `_, which is a separate program by the Jensen group that predicts the pKa of each individual residue in your structure. Prodes does not need it to run, but you should use it: it takes seconds, it works on Windows, macOS and Linux, and it makes the charge-related features considerably more realistic. See `pKa values and protonation states`_. If you prefer, ``pip install propka`` does the same job. +The fourth line adds two separate programs. Prodes needs neither to run, and you should use both. + +`PROPKA `_, from the Jensen group, predicts the pKa of each individual residue in your structure. It takes seconds, it works on Windows, macOS and Linux, and it makes the charge-related features considerably more realistic. See `pKa values and protonation states`_. + +`PDB2PQR `_, from the Electrostatics Consortium, repairs the structure before anything measures it: it rebuilds side chains that were not modelled, adds hydrogens, and reports the damage it cannot repair. 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`_. + +Install PDB2PQR from conda-forge rather than with ``pip``: the PyPI package pins ``docutils<0.18``, which collides with Sphinx and several other common packages, and the conda-forge build does not. Installing as a developer ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -63,7 +69,7 @@ Installing as a developer conda activate prodes pre-commit install -Note that ``environment.yml`` already installs Prodes itself, in editable mode and includes PROPKA. +Note that ``environment.yml`` already installs Prodes itself, in editable mode, and includes both PROPKA and PDB2PQR. Check that it worked: @@ -95,18 +101,25 @@ This is deliberately not in ``.pre-commit-config.yaml``, which everyone runs: gg Quick start ------------ -Activate the environment, then run these three commands on your structure: +Activate the environment, then run these four commands on your structure: .. code-block:: text conda activate prodes + mkdir prepared - propka3 1GDW.pdb # predict per-residue pKa values + pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr # repair the structure + propka3 prepared/1GDW.pdb # predict per-residue pKa values python -m prodes.io.pka_converter 1GDW.pka propka -o 1GDW_pka.json # convert them for Prodes - python -m prodes 1GDW.pdb 1GDW.zip --ph 7.4 --pka 1GDW_pka.json # calculate the features + 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. +Two things about that sequence are worth noticing, and both are explained under `Preparing your structure`_: + +* **Everything after the** ``pdb2pqr`` **line runs on the repaired structure**, ``prepared/1GDW.pdb``, not on the file you downloaded. That includes PROPKA. +* **The repaired file keeps the original name, in a new directory.** Prodes takes the ``ID`` column from the file name, so writing ``1GDW_prep.pdb`` instead would label that row ``1GDW_prep``. PROPKA writes its ``.pka`` into the directory you are standing in rather than next to its input, which is why the third line reads ``1GDW.pka`` and not ``prepared/1GDW.pka``. + See `Viewing the surface`_ to look at the result. PROPKA works out the pKa of every titratable residue *in its actual structural context*, because a buried aspartate and an exposed one titrate at quite different pH values. Without it (leave off ``--pka``), Prodes falls back to one textbook pKa per residue type. See `pKa values and protonation states`_. @@ -310,7 +323,7 @@ Using the features in a model This is the point of the whole exercise, so here is the shape of a complete QSPR workflow. -**Step 1, calculate features for every structure you have measured.** Loop over a folder, predict pKa values for each structure, and let the feature rows accumulate in one CSV: +**Step 1, calculate features for every structure you have measured.** Loop over a folder, repair each structure, predict pKa values for it, and let the feature rows accumulate in one CSV: .. code-block:: python @@ -320,16 +333,31 @@ This is the point of the whole exercise, so here is the shape of a complete QSPR import prodes from prodes.io.pka_converter import convert_propka, write_json + prepared = Path("prepared") bundles = Path("bundles") + prepared.mkdir(exist_ok=True) + bundles.mkdir(exist_ok=True) for pdb in sorted(Path("structures").glob("*.pdb")): - subprocess.run(["propka3", pdb.name], cwd=pdb.parent, check=True) - pka_json = pdb.with_suffix(".pka.json") - write_json(convert_propka(str(pdb.with_suffix(".pka"))), str(pka_json)) + repaired = prepared / pdb.name + subprocess.run( + ["pdb2pqr", "--ff=PARSE", f"--pdb-output={repaired}", + str(pdb), str(repaired.with_suffix(".pqr"))], + check=True, + ) + subprocess.run(["propka3", repaired.name], cwd=prepared, check=True) + + pka_json = prepared / f"{pdb.stem}_pka.json" + write_json(convert_propka(str(prepared / f"{pdb.stem}.pka")), str(pka_json)) + + prodes.run_prodes(str(repaired), str(bundles / f"{pdb.stem}.zip"), + ph=7.4, pkas_file=str(pka_json)) - prodes.run_prodes(str(pdb), str(bundles / f"{pdb.stem}.zip"), ph=7.4, pkas_file=str(pka_json)) +The repaired file keeps the original name inside ``prepared/``, so the ``ID`` column still reads ``1GDW`` rather than ``1GDW_prep``. See `Preparing your structure`_. -``check=True`` matters here: without it a structure PROPKA chokes on would fail silently, and the loop would carry on and calculate that protein with default pKa values instead. You would end up with one row in the table quietly computed on a different basis from all the others. +``check=True`` matters in both calls: without it a structure PROPKA chokes on would fail silently, and the loop would carry on and calculate that protein with default pKa values instead. You would end up with one row in the table quietly computed on a different basis from all the others. + +PDB2PQR refuses some structures outright rather than repairing them, and with ``check=True`` that stops the loop, which is what you want. Decide what to do with a refused structure explicitly: repair it another way, run Prodes on the unrepaired file and record that you did, or drop it. What you must not do is let half your dataset be repaired and the other half not. **Step 2, join the features to your measurements** on the ``ID`` column. Your measurement file needs an ``ID`` column whose values match the structure file names (see `Input files`_): @@ -399,6 +427,68 @@ The values are not quite identical to the released 4.x numbers, for a separate r **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``. +Preparing your structure +------------------------- + +**Prodes describes the structure you give it. It does not correct one.** If a surface lysine has no side chain in the file, Prodes gives that residue no charge and says nothing about it, and every charge-derived feature comes out wrong. + +`PDB2PQR `_ repairs the structure first. It is free, BSD-3-Clause, pure Python, and one command. + +.. code-block:: text + + mkdir prepared + pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr + +``--pdb-output`` is the flag that matters. The ``.pqr`` file is PDB2PQR's normal output and Prodes cannot read it; ``--pdb-output`` writes the repaired structure as a PDB, which Prodes can. Run everything after this point, PROPKA included, on the repaired file. + +**Keep the original file name and change the directory, not the other way round.** Prodes takes the ``ID`` column of the output from the file name it was given (see `Input files`_), so a repaired file called ``1GDW_prep.pdb`` labels that row ``1GDW_prep``. Over a dataset that is tedious to undo. + +Prodes will also print ``Ignoring unrecognized record 'TER'`` on a PDB2PQR output file. That is harmless: ``TER`` is a chain separator and carries no coordinates. + +How much does it matter? +~~~~~~~~~~~~~~~~~~~~~~~~~ + +PDB entry 4HKZ is a Fab. Two of its lysines are modelled only as far as CB, which is ordinary for surface residues at that resolution. The only difference between these two columns is whether those two side chains are present: + +===================================== ======================= ===================== +Feature side chains unmodelled side chains rebuilt +===================================== ======================= ===================== +Formal charge at pH 7.4 -2 0 +Isoelectric point 6.46 7.61 +Mean surface electrostatic potential -0.013 +0.011 +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. + +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. + +What it fixes +~~~~~~~~~~~~~~ + +PDB2PQR rebuilds missing heavy atoms, adds every hydrogen, optimises the hydrogen-bond network (including flipping Asn, Gln and His where the deposited assignment is the worse one), resolves the steric clashes that creates, adds a missing ``OXT``, detects backbone chain breaks and names them, and refuses to continue on a gap too large to rebuild rather than describing a hole in silence. Prodes does none of that. + +Prodes keeps two jobs itself and does them well: it recognises disulfide bonds, so a cystine is not titrated as a free thiol, and it resolves alternate conformations by keeping the best occupied one per residue. + +**Prodes uses only the geometry from the prepared file.** It reads ``ATOM`` records, filters out hydrogens, and re-titrates every residue from its own pKa table at the pH you ask for, so the hydrogens, protonation states and partial charges PDB2PQR writes are all discarded. That is deliberate: it is what lets you run Prodes at several pH values against one prepared structure. + +You still need PROPKA +~~~~~~~~~~~~~~~~~~~~~~ + +PDB2PQR can run PROPKA itself, with ``--titration-state-method=propka``. **That does not replace the** ``propka3`` **step.** PDB2PQR does not write a ``.pka`` file: it uses the values internally to choose protonation states and then discards them, so ``prodes.io.pka_converter`` has nothing to read. It also bakes in a single pH, where Prodes deliberately separates prediction from calculation so you can predict once and run at many pH values. + +What to take from PDB2PQR is the **ordering**: it runs PROPKA after repairing the structure, and that is the right way round. On 4HKZ, PROPKA finds 145 titratable groups in the deposited file and 149 in the prepared one, the four extra being the two rebuilt lysine side chains and the two C-termini that only exist once ``OXT`` has been added. Six of the 145 shared groups shift by more than 0.5 pKa units. + +Traps +~~~~~~ + +* **Do not use** ``--ffout``. It renames residues into the force field's own scheme, putting ``ASH``, ``GLH``, ``LYN`` and ``HID`` into the file, and Prodes raises ``KeyError`` on any of those. Without it, both PARSE and AMBER keep canonical residue names. +* **PDB2PQR can fail on sequence microheterogeneity.** Crambin (1CBN) has residue 22 modelled as both PRO and SER, and PDB2PQR 3.6.1 exits with ``Unable to debump biomolecule``. Prodes handles that case correctly on its own, so running it on the original file is a reasonable fallback when PDB2PQR refuses a structure. +* **mmCIF input is unreliable.** PDB2PQR 3.6.1 accepts a ``.cif`` but produced an empty ``.pqr`` from a valid one in testing. Convert to PDB first. +* **Be consistent within a dataset**, exactly as for PROPKA. Prepare all of your structures or none of them. + +Full detail, including what Prodes discards and why: `preparing a structure `_. + pKa values and protonation states ---------------------------------- @@ -426,14 +516,20 @@ It goes into the same environment as Prodes, and is already included if you buil or equivalently ``pip install propka``. Check it with ``propka3 --version``. +PDB2PQR installs the same way, ``conda install conda-forge::pdb2pqr``, and is checked with ``pdb2pqr --version``. Prefer conda-forge to ``pip`` for that one; see `Installing as a user`_. + The three steps ~~~~~~~~~~~~~~~~ -**Step 1, predict the pKa values.** PROPKA reads your PDB file and writes a ``.pka`` file beside it: +**Step 1, predict the pKa values.** PROPKA reads your PDB file and writes a ``.pka`` file: .. code-block:: text - propka3 1GDW.pdb # writes 1GDW.pka + propka3 prepared/1GDW.pdb # writes 1GDW.pka in the current directory + +Point it at the **prepared** structure, not at the file you downloaded. A residue whose side chain is missing has no pKa to predict, so running PROPKA before the repair silently drops those groups. See `Preparing your structure`_. + +Note that PROPKA writes the ``.pka`` file into the directory you are standing in, named after the input's base name, rather than beside the input file. **Step 2, convert that file to the Prodes pKa JSON format.** Prodes ships converters for three predictors, usable from the command line @@ -688,7 +784,7 @@ What Prodes offers, then, is: 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.** 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. It does not form bonds, rebuild side chains, cap chain breaks or touch any coordinate. Preparing the structure is otherwise your responsibility. +* **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. 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. diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md new file mode 100644 index 0000000..528b443 --- /dev/null +++ b/docs/structure_preparation.md @@ -0,0 +1,107 @@ +# Preparing a structure before running Prodes + +Prodes describes the structure you give it. It does not correct one. If a surface lysine has no side chain in the file, Prodes gives that residue no charge and says nothing about it, and every charge-derived feature comes out wrong. + +[PDB2PQR](https://github.com/Electrostatics/pdb2pqr) fixes most of that. It is free, BSD-3-Clause, pure Python, and a single command. This page is why you should run it, what it changes, and the traps. + +## Why it matters + +PDB entry **4HKZ** is a Fab. Two of its lysines, A188 and B221, are modelled only as far as CB, which is ordinary for surface residues at that resolution. Prodes reads the file without comment. + +Running PDB2PQR first rebuilds those two side chains. The only difference between these two columns is whether those two side chains are present: + +| | side chains unmodelled | side chains rebuilt | +| --- | --- | --- | +| Formal charge at pH 7.4 | **-2** | **0** | +| Isoelectric point | **6.46** | **7.61** | +| Mean surface electrostatic potential | -0.013 V | **+0.011 V** | +| 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, an isoelectric point of 6.46 and one of 7.61 are different recommendations. + +This is not universal. Many deposited structures are complete, and 1AO6, 1GPB, 4NZU, 1GDW, 1CBN, 1IGT and 1BEY all are. 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. + +## What PDB2PQR does that Prodes does not + +| PDB2PQR | Prodes on its own | +| --- | --- | +| **Rebuilds missing heavy atoms** from the residue's reference geometry | nothing; the atom stays absent and everything computed from it is wrong | +| **Detects backbone chain breaks** and names them | nothing; a nine residue gap in 1BEY passes without a message | +| **Refuses to continue** when a gap is too large to rebuild, naming the residue | nothing | +| **Deletes atoms that do not belong to a residue**, logging each one | keeps them, and assigns a radius from the first letter of the atom name | +| **Adds every hydrogen**, including polar hydrogens whose position depends on the hydrogen bond network | reads hydrogens if present, then filters them out; every calculation is heavy-atom only | +| **Optimises the hydrogen bond network**, including flipping Asn, Gln and His where the deposited assignment is worse than the flip | nothing | +| **Resolves steric clashes** created by rebuilt atoms, by rotating side chain dihedrals | nothing | +| **Adds a missing `OXT`** and treats termini as proper patches, with neutral variants available | assigns the terminus but does not add `OXT`, and places the C-terminal charge on the backbone carbonyl carbon | +| **Assigns a partial charge and a radius to every atom** from a named force field | formal charges on ionisable groups only, and one radius per element from a five-entry table | +| **Reads 57 PDB record types**, including `MODEL`, `MODRES` and `LINK` | reads coordinates and `SSBOND` | +| **Fetches a structure by PDB ID** | reads a local `.pdb` or `.pdb.zip` | + +Prodes keeps hold of two of these itself, and does them well: it **recognises disulfide bonds**, so a cystine is not titrated as a free thiol, and it **resolves alternate conformations** by keeping the best occupied one per residue. + +## What Prodes actually uses from the prepared file + +Only the **geometry**. Prodes reads `ATOM` records, filters out hydrogens, and re-titrates every residue from its own pKa table at the pH you ask for. + +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. + +## The commands + +``` +conda activate prodes +conda install conda-forge::pdb2pqr # once + +mkdir prepared +pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr +propka3 prepared/1GDW.pdb # the prepared file, not the download +python -m prodes.io.pka_converter 1GDW.pka propka -o 1GDW_pka.json +python -m prodes prepared/1GDW.pdb 1GDW.zip --ph 7.4 --pka 1GDW_pka.json +``` + +`--pdb-output` is the flag that matters. The `.pqr` file is PDB2PQR's normal output and Prodes cannot read it; `--pdb-output` writes the prepared structure as a PDB, which Prodes can. + +**Keep the original file name and change the directory.** Prodes takes the `ID` column from the file name, so a prepared file called `1GDW_prep.pdb` labels that row `1GDW_prep`. Over a dataset that is tedious to undo. + +**PROPKA writes its `.pka` into the directory you are standing in**, named after the input's base name, not next to the input file. That is why the third line above reads `1GDW.pka` and not `prepared/1GDW.pka`. + +Prodes will print `Ignoring unrecognized record 'TER'` on a PDB2PQR output file. That is harmless: `TER` is a chain separator and carries no coordinates. + +Install from **conda-forge**, not pip. The PyPI package pins `docutils<0.18`, which collides with Sphinx and with several other common packages. The conda-forge build does not carry that pin. + +## Run PROPKA on the prepared file, not on the download + +PDB2PQR can run PROPKA itself, with `--titration-state-method=propka`. **That does not replace the `propka3` step**, for two reasons: + +* PDB2PQR does not write a `.pka` file. It uses the values internally to pick protonation states and then discards them, so `prodes.io.pka_converter` has nothing to read. +* PDB2PQR bakes in one pH. Prodes separates prediction from calculation so that you can predict once and then run at as many pH values as you like. + +What you should take from PDB2PQR is the **ordering**. It runs PROPKA after repairing the structure, and that is the right way round. On 4HKZ: + +| | titratable groups predicted | +| --- | --- | +| PROPKA on the deposited file | 145 | +| PROPKA on the prepared file | 149 | + +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. + +## Traps + +**Do not use `--ffout`.** It renames residues into the force field's own scheme, putting `ASH`, `GLH`, `LYN` and `HID` into the file. Prodes raises `KeyError` on any of those. Without `--ffout` both PARSE and AMBER keep canonical residue names, verified at pH 4.0 and pH 7.4. + +**PDB2PQR applies terminal patches at chain breaks.** It puts a real +1 and -1 at a position where the chain is actually continuous and merely unmodelled. Prodes assigns termini only at the first and last residue of a chain, which is the better answer, and it keeps doing so on the prepared file because it re-titrates from scratch. The charges PDB2PQR invented at the gap are discarded along with the rest of its charge assignment. + +**PDB2PQR can fail on sequence microheterogeneity.** Crambin (1CBN) has residue 22 modelled as both PRO and SER. PDB2PQR 3.6.1 picks the highest occupancy per *atom* without keeping the residue identity consistent, ends up with a proline carrying a serine hydroxyl, and exits with `Unable to debump biomolecule`. Prodes handles this case correctly on its own, so if PDB2PQR refuses a structure, running Prodes on the original file is a reasonable fallback. + +**mmCIF input does not work reliably.** PDB2PQR 3.6.1 accepts a `.cif` but produced an empty `.pqr` from a valid one during testing. Convert to PDB first. + +**Be consistent within a dataset.** The same rule as for PROPKA: prepare all of your structures or none of them. Mixing prepared and unprepared structures puts two different kinds of number in the same feature column. + +## If you would rather not add the dependency + +The single most valuable thing PDB2PQR gives you here is knowing that a residue is incomplete. Until Prodes reports that itself, you can check without installing anything by comparing each residue's atom names against the standard set for its type, and treating any structure with missing side-chain atoms on charged residues as suspect. + +## Reference + +Jurrus E, et al. Improvements to the APBS biomolecular solvation software suite. *Protein Sci* 27, 112-128 (2018). https://doi.org/10.1002/pro.3280 diff --git a/environment.yml b/environment.yml index 63e44bc..4ba80a6 100644 --- a/environment.yml +++ b/environment.yml @@ -14,6 +14,11 @@ dependencies: # per-residue pKa prediction is the recommended default workflow (see README), # and asking every user to install it separately just means most of them won't. - propka=3.5.1 + # Recommended companion tool, not an import-time dependency: prodes never + # imports pdb2pqr, it reads the repaired PDB that pdb2pqr writes. Included + # here because running it first is the recommended workflow (see README), and + # because an unrepaired structure changes feature values silently. + - pdb2pqr=3.6.1 - pip=26.0.1 - setuptools=82.0.1 - wheel=0.47.0 From f4d24b094e391eddfafc94d33660a643df4ef0d9 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 00:39:15 +0200 Subject: [PATCH 2/8] Warn about the multi-chain pKa warnings and about selecting features outside cross-validation Two notes in the README for users building a dataset, both about failure modes that are quiet rather than loud. The pKa step prints hundreds of warnings on any structure with more than one chain, because pKa values are keyed by residue number with no chain, so every chain's predictions are offered to every chain's residue with that number and most are rejected. A user meeting this for the first time on an antibody cannot tell whether their file is wrong. The note says the warnings are expected, says what the underlying defect can actually do to a value, and points at the issue that tracks it. The second note covers feature selection. Screening the 54 columns against a measurement and then cross-validating a model built from the survivors reports a score that is not out of sample, because every held-out protein helped choose the features. At the sample sizes typical of this field, tens of proteins rather than thousands, that produces impressive numbers from data with no signal in it. The note gives the correction, which is to refit the selection on each training fold. Neither note changes any behaviour. --- README.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.rst b/README.rst index a5c3568..b362670 100644 --- a/README.rst +++ b/README.rst @@ -393,6 +393,8 @@ PLS is a reasonable first choice for this kind of data because the features are **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. +**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. + Ionic strength and screening ----------------------------- @@ -564,6 +566,10 @@ Steps 1 and 2 are done **once per structure**. Step 3 can then be repeated as of Residues that appear in the file get the predicted value; every other residue keeps its default. So a prediction covering only the titratable residues, which is what these tools produce, is complete as far as Prodes is concerned. +.. note:: + + On a structure with more than one chain, such as an antibody, this step prints a long list of warnings of the form ``the pKa file gives THR 31 a ASP pKa, which that residue does not have; ignoring it``. Prodes currently keys pKa values by residue number alone, so every chain's predictions are offered to every chain's residue with that number and most of them are rejected. The warnings are noise rather than a problem with your file, but a genuine warning is easy to miss among them, and a residue that collides with the same residue type in another chain can take the wrong chain's value. Tracked as `issue #12 `_. + Cysteines and disulfide bonds ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From 508e234c2f3d8d2a253c348c68f7ff4921ce3eb5 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 02:19:53 +0200 Subject: [PATCH 3/8] Install PDB2PQR in its own environment instead of alongside Prodes Prodes never imports PDB2PQR. It only reads the repaired PDB file that PDB2PQR writes, so the two programs never have to be importable at the same time and installing them together buys nothing. What it costs is that each one's pins constrain the other's for as long as both are installed. They are independently maintained projects on separate release cycles, and the first time one moves a shared dependency the other has not caught up with, a single environment stops solving and takes a working Prodes install down with it. pdb2pqr is therefore removed from environment.yml, and the documented workflow creates a separate environment for it: conda create -n pdb2pqr -c conda-forge pdb2pqr The quick start now activates that environment for the repair step and the Prodes environment for everything after it. The worked loop under "Using the features in a model" calls conda run -n pdb2pqr instead, so a script runs start to finish without switching environments by hand. Installing into one environment with pip is documented as the alternative for anyone who prefers it, with the note that the PyPI package pins docutils below 0.18 while the conda-forge build does not. Also adds the second reference that PDB2PQR asks for and prints on every run, the 2007 Nucleic Acids Research paper describing PDB2PQR itself, alongside the 2018 Protein Science paper on the wider software suite. Both were checked against CrossRef and PubMed: the Protein Science article is volume 27, issue 1, pages 112-128, in the January 2018 issue, published online in October 2017. --- README.rst | 47 ++++++++++++++++++++++++----------- docs/structure_preparation.md | 27 +++++++++++++++++--- environment.yml | 5 ---- 3 files changed, 56 insertions(+), 23 deletions(-) diff --git a/README.rst b/README.rst index b362670..07f16e7 100644 --- a/README.rst +++ b/README.rst @@ -45,17 +45,28 @@ Installing as a user conda create -n prodes python=3.13 conda activate prodes pip install git+https://github.com/datacatalysis/prodes.git - conda install conda-forge::propka conda-forge::pdb2pqr + conda install conda-forge::propka Prodes depends on NumPy, pandas and `Biopython `_, which ``pip`` installs with it. Biopython reads the coordinate records of a PDB file; the rules that decide which alternate conformation of a residue to keep, and which cysteines are bonded, remain prodes' own. See `how prodes reads a PDB file `_. -The fourth line adds two separate programs. Prodes needs neither to run, and you should use both. +The fourth line adds `PROPKA `_, which is a separate program by the Jensen group that predicts the pKa of each individual residue in your structure. Prodes does not need it to run, but you should use it: it takes seconds, it works on Windows, macOS and Linux, and it makes the charge-related features considerably more realistic. See `pKa values and protonation states`_. If you prefer, ``pip install propka`` does the same job. + +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`_. + +**Put it in an environment of its own:** + +.. code-block:: text + + conda create -n pdb2pqr -c conda-forge pdb2pqr -`PROPKA `_, from the Jensen group, predicts the pKa of each individual residue in your structure. It takes seconds, it works on Windows, macOS and Linux, and it makes the charge-related features considerably more realistic. See `pKa values and protonation states`_. +Prodes never imports PDB2PQR. It only reads the repaired PDB file that PDB2PQR writes, so the two programs never have to be importable at the same time, and installing them together buys nothing. What it costs is that each one's pins constrain the other's for as long as both are installed. These are two independently maintained projects on separate release cycles, and the first time one of them moves a shared dependency the other has not caught up with, a single environment stops solving and takes your working Prodes install down with it. -`PDB2PQR `_, from the Electrostatics Consortium, repairs the structure before anything measures it: it rebuilds side chains that were not modelled, adds hydrogens, and reports the damage it cannot repair. 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`_. +Keeping them apart means an upgrade to either is someone else's problem. -Install PDB2PQR from conda-forge rather than with ``pip``: the PyPI package pins ``docutils<0.18``, which collides with Sphinx and several other common packages, and the conda-forge build does not. +If you would rather have one environment anyway, ``pip install pdb2pqr`` into the Prodes environment usually works. Be aware that the PyPI package pins ``docutils<0.18``, which collides with Sphinx and several other common packages; the conda-forge build does not carry that pin, which is the other reason the separate conda environment is the cleaner route. Installing as a developer ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -69,7 +80,7 @@ Installing as a developer conda activate prodes pre-commit install -Note that ``environment.yml`` already installs Prodes itself, in editable mode, and includes both PROPKA and PDB2PQR. +Note that ``environment.yml`` already installs Prodes itself, in editable mode and includes PROPKA. It deliberately does not include PDB2PQR; see `PDB2PQR, in its own environment`_. Check that it worked: @@ -101,14 +112,16 @@ This is deliberately not in ``.pre-commit-config.yaml``, which everyone runs: gg Quick start ------------ -Activate the environment, then run these four commands on your structure: +Run these four commands on your structure. PDB2PQR lives in its own environment, so the first one is run with that environment active and the rest with the Prodes environment active: .. code-block:: text - conda activate prodes mkdir prepared + conda activate pdb2pqr pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr # repair the structure + + conda activate prodes propka3 prepared/1GDW.pdb # predict per-residue pKa values 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 @@ -117,7 +130,7 @@ That writes ``1GDW.zip``, a bundle holding the 54 features, the surface points t Two things about that sequence are worth noticing, and both are explained under `Preparing your structure`_: -* **Everything after the** ``pdb2pqr`` **line runs on the repaired structure**, ``prepared/1GDW.pdb``, not on the file you downloaded. That includes PROPKA. +* **Everything after the** ``pdb2pqr`` **line runs on the repaired structure**, ``prepared/1GDW.pdb``, not on the file you downloaded. That includes PROPKA. In a script, ``conda run -n pdb2pqr pdb2pqr ...`` avoids switching environments by hand. * **The repaired file keeps the original name, in a new directory.** Prodes takes the ``ID`` column from the file name, so writing ``1GDW_prep.pdb`` instead would label that row ``1GDW_prep``. PROPKA writes its ``.pka`` into the directory you are standing in rather than next to its input, which is why the third line reads ``1GDW.pka`` and not ``prepared/1GDW.pka``. See `Viewing the surface`_ to look at the result. @@ -341,8 +354,8 @@ This is the point of the whole exercise, so here is the shape of a complete QSPR for pdb in sorted(Path("structures").glob("*.pdb")): repaired = prepared / pdb.name subprocess.run( - ["pdb2pqr", "--ff=PARSE", f"--pdb-output={repaired}", - str(pdb), str(repaired.with_suffix(".pqr"))], + ["conda", "run", "-n", "pdb2pqr", "pdb2pqr", "--ff=PARSE", + f"--pdb-output={repaired}", str(pdb), str(repaired.with_suffix(".pqr"))], check=True, ) subprocess.run(["propka3", repaired.name], cwd=prepared, check=True) @@ -353,7 +366,7 @@ This is the point of the whole exercise, so here is the shape of a complete QSPR prodes.run_prodes(str(repaired), str(bundles / f"{pdb.stem}.zip"), ph=7.4, pkas_file=str(pka_json)) -The repaired file keeps the original name inside ``prepared/``, so the ``ID`` column still reads ``1GDW`` rather than ``1GDW_prep``. See `Preparing your structure`_. +``conda run -n pdb2pqr`` calls PDB2PQR in its own environment, so this loop runs start to finish with the Prodes environment active. The repaired file keeps the original name inside ``prepared/``, so the ``ID`` column still reads ``1GDW`` rather than ``1GDW_prep``. See `Preparing your structure`_. ``check=True`` matters in both calls: without it a structure PROPKA chokes on would fail silently, and the loop would carry on and calculate that protein with default pKa values instead. You would end up with one row in the table quietly computed on a different basis from all the others. @@ -434,12 +447,18 @@ Preparing your structure **Prodes describes the structure you give it. It does not correct one.** If a surface lysine has no side chain in the file, Prodes gives that residue no charge and says nothing about it, and every charge-derived feature comes out wrong. -`PDB2PQR `_ repairs the structure first. It is free, BSD-3-Clause, pure Python, and one command. +`PDB2PQR `_ repairs the structure first. It is free, BSD-3-Clause, pure Python, and one command. Install it into an environment of its own rather than alongside Prodes; see `PDB2PQR, in its own environment`_ for why. .. code-block:: text + conda create -n pdb2pqr -c conda-forge pdb2pqr # once + mkdir prepared + conda activate pdb2pqr pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr + conda activate prodes + +Inside a script, ``conda run -n pdb2pqr pdb2pqr ...`` does the same thing without switching environments by hand. ``--pdb-output`` is the flag that matters. The ``.pqr`` file is PDB2PQR's normal output and Prodes cannot read it; ``--pdb-output`` writes the repaired structure as a PDB, which Prodes can. Run everything after this point, PROPKA included, on the repaired file. @@ -518,7 +537,7 @@ It goes into the same environment as Prodes, and is already included if you buil or equivalently ``pip install propka``. Check it with ``propka3 --version``. -PDB2PQR installs the same way, ``conda install conda-forge::pdb2pqr``, and is checked with ``pdb2pqr --version``. Prefer conda-forge to ``pip`` for that one; see `Installing as a user`_. +PDB2PQR goes into a separate environment rather than this one, and is checked with ``conda run -n pdb2pqr pdb2pqr --version``. See `PDB2PQR, in its own environment`_. The three steps ~~~~~~~~~~~~~~~~ diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index 528b443..658e4d4 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -49,17 +49,32 @@ So the hydrogens PDB2PQR adds, the protonation states it assigns and the partial ## The commands +Install PDB2PQR into an environment of its own, not alongside Prodes: + ``` -conda activate prodes -conda install conda-forge::pdb2pqr # once +conda create -n pdb2pqr -c conda-forge pdb2pqr # once +``` + +Prodes never imports PDB2PQR. It only reads the repaired PDB file that PDB2PQR writes, so the two never have to be importable at the same time, and installing them together buys nothing. What it costs is that each one's pins constrain the other's for as long as both are installed. They are independently maintained projects on separate release cycles, and the first time one of them moves a shared dependency the other has not caught up with, a single environment stops solving and takes a working Prodes install down with it. +A `pip install pdb2pqr` into the Prodes environment usually works if you would rather have one environment. Note that the PyPI package pins `docutils<0.18`, which collides with Sphinx and several other common packages; the conda-forge build does not carry that pin. + +Then: + +``` mkdir prepared + +conda activate pdb2pqr pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr + +conda activate prodes propka3 prepared/1GDW.pdb # the prepared file, not the download python -m prodes.io.pka_converter 1GDW.pka propka -o 1GDW_pka.json python -m prodes prepared/1GDW.pdb 1GDW.zip --ph 7.4 --pka 1GDW_pka.json ``` +Inside a script, `conda run -n pdb2pqr pdb2pqr ...` calls PDB2PQR in its own environment without switching by hand, so a loop can run start to finish with the Prodes environment active. + `--pdb-output` is the flag that matters. The `.pqr` file is PDB2PQR's normal output and Prodes cannot read it; `--pdb-output` writes the prepared structure as a PDB, which Prodes can. **Keep the original file name and change the directory.** Prodes takes the `ID` column from the file name, so a prepared file called `1GDW_prep.pdb` labels that row `1GDW_prep`. Over a dataset that is tedious to undo. @@ -102,6 +117,10 @@ The four extra groups are the two rebuilt lysine side chains and the two C-termi The single most valuable thing PDB2PQR gives you here is knowing that a residue is incomplete. Until Prodes reports that itself, you can check without installing anything by comparing each residue's atom names against the standard set for its type, and treating any structure with missing side-chain atoms on charged residues as suspect. -## Reference +## References + +PDB2PQR asks to be cited as both of these, and prints both on every run: + +Dolinsky TJ, Czodrowski P, Li H, Nielsen JE, Jensen JH, Klebe G, Baker NA. PDB2PQR: expanding and upgrading automated preparation of biomolecular structures for molecular simulations. *Nucleic Acids Research* **35**(Web Server issue), W522-W525 (2007). https://doi.org/10.1093/nar/gkm276 -Jurrus E, et al. Improvements to the APBS biomolecular solvation software suite. *Protein Sci* 27, 112-128 (2018). https://doi.org/10.1002/pro.3280 +Jurrus E, Engel D, Star K, Monson K, Brandi J, Felberg LE, et al. Improvements to the APBS biomolecular solvation software suite. *Protein Science* **27**(1), 112-128 (2018). https://doi.org/10.1002/pro.3280 diff --git a/environment.yml b/environment.yml index 4ba80a6..63e44bc 100644 --- a/environment.yml +++ b/environment.yml @@ -14,11 +14,6 @@ dependencies: # per-residue pKa prediction is the recommended default workflow (see README), # and asking every user to install it separately just means most of them won't. - propka=3.5.1 - # Recommended companion tool, not an import-time dependency: prodes never - # imports pdb2pqr, it reads the repaired PDB that pdb2pqr writes. Included - # here because running it first is the recommended workflow (see README), and - # because an unrepaired structure changes feature values silently. - - pdb2pqr=3.6.1 - pip=26.0.1 - setuptools=82.0.1 - wheel=0.47.0 From 7efea15c8f0f60d08cad1498af6ccfda706dfe1e Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 08:40:47 +0200 Subject: [PATCH 4/8] Ignore the per-repo direnv file and single-investigation Python scripts .envrc holds the per-repo GitHub token routing and is machine-local, so it was already untracked in practice but not ignored. *_notrack.py matches the existing *_notrack.md convention: a script written to answer one question, superseded as soon as the answer is recorded, and not something a reader in five years should have to work out the purpose of. --- .gitignore | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 6b178f3..8bbb240 100644 --- a/.gitignore +++ b/.gitignore @@ -21,8 +21,10 @@ docs/similar_tools/* # wrong. Anything worth keeping belongs in a docstring, the README, or an issue. # Name them so the intent is visible in a file listing, for example # plan_thing_notrack.md or plan_thing_REVIEW1_notrack.md. Use git add -f to -# commit one deliberately. +# commit one deliberately. The same applies to a script written to answer one +# question, which is superseded as soon as its answer is recorded. *_notrack.md +*_notrack.py # APBS comparison outputs and the capture file APBS writes into its working directory apbs_results From fb54c05d7000e1c7e62086717a69a90e805241d5 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 08:41:43 +0200 Subject: [PATCH 5/8] Record what repairing a predicted structure does, and how the terminal oxygen is scored Running PDB2PQR over 855 structures, 826 Boltz2 predicted monomers and 29 crystal structures, repaired 854 and refused 1, at a median of about one second per predicted structure. The refusal was correct: 5CHA models a residue with its backbone nitrogen and no other atom, which cannot be rebuilt from its neighbours. Prodes reads that same file without a message. Predicted structures turn out to need exactly one atom each, the C terminal OXT, which Boltz2 does not write. That is worth documenting for two reasons. Structure sources disagree about whether to write OXT at all. Over the structures measured, every AlphaFold database model has one, no Boltz2 monomer has one, and 21 of 29 crystal structures have one. So the same protein arrives with a different atom count depending on where its structure came from. Prodes assigns hydrophobicity per residue and gives every heavy atom the value of its residue, with one exception: OXT is given a fixed 1.0, which on the default mj_scaled scale is the maximum, the same as phenylalanine. A carboxylate oxygen is strongly hydrophilic, so the value is inverted. The constant has been present since the first commit and this fork preserves it along with the rest of the original algorithm. Measured on 40 Boltz2 models, adding the OXT moves a median of 9 of the 54 default features. Most of that is the ordinary effect of one more atom on the surface. Two hydrophobicity features move in one direction: NSurfPosMhp rises on 33 of the 39 structures where it moves, and SurfMhpMean rises on 12 of 13. The size of the shift is small. The problem is that it currently applies to every AlphaFold structure, no Boltz2 structure and most crystal structures, so a model trained across sources partly learns which predictor produced each file. Repairing every structure removes that, because every structure then carries exactly one OXT per chain and the artefact becomes a constant rather than a per-source difference. No code changes. The constant itself is left alone here. --- README.rst | 6 +++++- docs/structure_preparation.md | 26 ++++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 07f16e7..c56405a 100644 --- a/README.rst +++ b/README.rst @@ -484,6 +484,10 @@ Positive surface points 12,976 1 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. +Predicted structures are a different case. They are complete by construction, and running PDB2PQR over 826 Boltz2 models added exactly one atom to each and nothing else: ``OXT``, the second oxygen of the C-terminal carboxylate, which Boltz2 does not write. That one atom still matters, because structure sources disagree about whether to write it at all and Prodes treats it specially. See `the terminal oxygen `_. + +**How often does the repair succeed?** Over 855 structures, 826 Boltz2 monomers and 29 crystal structures, PDB2PQR repaired 854 and refused 1, at a median of about one second per predicted structure. The refusal was correct: in ``5CHA`` a residue is modelled with its backbone nitrogen and no other atom, and a residue with one atom cannot be rebuilt. Prodes reads that same file without a word. + What it fixes ~~~~~~~~~~~~~~ @@ -506,7 +510,7 @@ Traps * **Do not use** ``--ffout``. It renames residues into the force field's own scheme, putting ``ASH``, ``GLH``, ``LYN`` and ``HID`` into the file, and Prodes raises ``KeyError`` on any of those. Without it, both PARSE and AMBER keep canonical residue names. * **PDB2PQR can fail on sequence microheterogeneity.** Crambin (1CBN) has residue 22 modelled as both PRO and SER, and PDB2PQR 3.6.1 exits with ``Unable to debump biomolecule``. Prodes handles that case correctly on its own, so running it on the original file is a reasonable fallback when PDB2PQR refuses a structure. * **mmCIF input is unreliable.** PDB2PQR 3.6.1 accepts a ``.cif`` but produced an empty ``.pqr`` from a valid one in testing. Convert to PDB first. -* **Be consistent within a dataset**, exactly as for PROPKA. Prepare all of your structures or none of them. +* **Be consistent within a dataset**, exactly as for PROPKA. Prepare all of your structures or none of them. This matters more than it looks: because structure predictors disagree about whether to write the C-terminal ``OXT``, a half-prepared dataset carries a systematic difference between its sources. See `the terminal oxygen `_. Full detail, including what Prodes discards and why: `preparing a structure `_. diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index 658e4d4..f063f9d 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -47,6 +47,32 @@ 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. +## 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. + +**What `OXT` is.** The C terminus of a protein is a carboxylate, `-C(=O)-O`. The PDB format names the first oxygen `O`, which every residue has, and the second one `OXT`, which only the last residue of a chain has. It is an ordinary heavy atom of the real molecule, not an artefact. + +**Structure sources disagree about whether to write it.** Measured over the structures at hand: + +| source | structures with `OXT` | +| --- | --- | +| AlphaFold database | 200 of 200 sampled | +| Boltz2 predicted monomers | 0 of 200 sampled | +| crystal structures | 21 of 29 | + +Boltz2 leaves it off, AlphaFold always writes it, and a crystal structure has it or not depending on whether the terminal residue was resolved. So the same protein arrives with a different atom count depending on where its structure came from. + +**Why that matters here.** Prodes assigns hydrophobicity per residue, giving every heavy atom the value of the residue it belongs to. `OXT` is the one exception: `Property_point.set_lipo` gives it a fixed value of 1.0, which on the default `mj_scaled` scale is the maximum, the same as phenylalanine. A carboxylate oxygen is one of the most hydrophilic groups in a protein, so the value is the wrong way round. + +This is an oversight in the original package rather than a decision, and an easy one to make: Prodes was written single-handed during a PhD and released under the MIT licence, and this is a single constant in one list comprehension. It has been in the code since the first commit and this fork preserves it, along with the rest of the original algorithm. It is recorded here because it is now measurable, not as a criticism. + +**The practical consequence is inconsistency, and that is what the repair fixes.** Adding one `OXT` to a Boltz2 model moves a median of 9 of the 54 default features. Most of that is the ordinary effect of one more atom on the surface, but two hydrophobicity features move in one direction: `NSurfPosMhp` rises on 33 of the 39 structures where it moves at all, and `SurfMhpMean` rises on 12 of the 13. + +The size of that shift is small. The problem is that today it is applied to every AlphaFold structure, no Boltz2 structure, and most but not all crystal structures. A model trained across sources is then partly learning which predictor produced each file. + +Repairing every structure with PDB2PQR removes that, because every structure ends up with exactly one `OXT` per chain. The artefact becomes a constant rather than a per-source bias, which is the thing that actually damages a model. Until the constant itself is changed, consistency is what matters, and this is the cheapest way to get it. + ## The commands Install PDB2PQR into an environment of its own, not alongside Prodes: From 07c0a3a657a8f620cd601a3a93b02a256388bd44 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 09:02:08 +0200 Subject: [PATCH 6/8] Shorten the note on how the terminal oxygen is scored The paragraph explaining where the fixed hydrophobicity value came from was longer than the fact it conveyed. One sentence naming the upstream repository says the same thing. --- docs/structure_preparation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index f063f9d..a8173ca 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -65,7 +65,7 @@ Boltz2 leaves it off, AlphaFold always writes it, and a crystal structure has it **Why that matters here.** Prodes assigns hydrophobicity per residue, giving every heavy atom the value of the residue it belongs to. `OXT` is the one exception: `Property_point.set_lipo` gives it a fixed value of 1.0, which on the default `mj_scaled` scale is the maximum, the same as phenylalanine. A carboxylate oxygen is one of the most hydrophilic groups in a protein, so the value is the wrong way round. -This is an oversight in the original package rather than a decision, and an easy one to make: Prodes was written single-handed during a PhD and released under the MIT licence, and this is a single constant in one list comprehension. It has been in the code since the first commit and this fork preserves it, along with the rest of the original algorithm. It is recorded here because it is now measurable, not as a criticism. +This behaviour is inherited from the original [tneijenhuis/prodes](https://github.com/tneijenhuis/prodes) repository, which this fork preserves. **The practical consequence is inconsistency, and that is what the repair fixes.** Adding one `OXT` to a Boltz2 model moves a median of 9 of the 54 default features. Most of that is the ordinary effect of one more atom on the surface, but two hydrophobicity features move in one direction: `NSurfPosMhp` rises on 33 of the 39 structures where it moves at all, and `SurfMhpMean` rises on 12 of the 13. From e572afe5883f2879d0812a4274d7b684d9f12808 Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 09:09:54 +0200 Subject: [PATCH 7/8] Pass --keep-chain when repairing a structure, and record which output it affects 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 applies per-residue pKa values by residue number. PDB2PQR writes two files and they behave differently. The file written by --pdb-output, which is the one Prodes reads, always carries the chain identifier. The .pqr file, which is PDB2PQR's normal output, leaves the chain column blank unless --keep-chain is given. The cause is in io.print_biomolecule_atoms: the PDB branch calls atom.get_pdb_string(), which takes no chain flag, while the PQR branch calls atom.get_pqr_string(chainflag=chainflag). So the flag is not required for this pipeline, and that was verified rather than assumed. On structures of two, four, six and eight chains the --pdb-output file carries every chain identifier, and running with and without --keep-chain produces byte-for-byte identical PDB output. On the six chain structure 1F6R, Prodes reads the same six chains, the same 727 residues and the same 24 disulfide bonds from the repaired file as from the original. The flag is added to the documented commands anyway. The .pqr file is what APBS and the tools built on it consume, so a reader who uses it needs the chain column, and nobody should have to remember which of two output files preserves what. A new section in the README and in docs/structure_preparation.md records the difference between the two outputs, so this does not have to be rediscovered. --- README.rst | 17 ++++++++++++++--- docs/structure_preparation.md | 21 ++++++++++++++++++++- 2 files changed, 34 insertions(+), 4 deletions(-) diff --git a/README.rst b/README.rst index c56405a..d7c3960 100644 --- a/README.rst +++ b/README.rst @@ -119,7 +119,7 @@ Run these four commands on your structure. PDB2PQR lives in its own environment, mkdir prepared conda activate pdb2pqr - pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr # repair the structure + pdb2pqr --ff=PARSE --keep-chain --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr # repair conda activate prodes propka3 prepared/1GDW.pdb # predict per-residue pKa values @@ -354,7 +354,7 @@ This is the point of the whole exercise, so here is the shape of a complete QSPR for pdb in sorted(Path("structures").glob("*.pdb")): repaired = prepared / pdb.name subprocess.run( - ["conda", "run", "-n", "pdb2pqr", "pdb2pqr", "--ff=PARSE", + ["conda", "run", "-n", "pdb2pqr", "pdb2pqr", "--ff=PARSE", "--keep-chain", f"--pdb-output={repaired}", str(pdb), str(repaired.with_suffix(".pqr"))], check=True, ) @@ -455,13 +455,15 @@ Preparing your structure mkdir prepared conda activate pdb2pqr - pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr + pdb2pqr --ff=PARSE --keep-chain --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr conda activate prodes Inside a script, ``conda run -n pdb2pqr pdb2pqr ...`` does the same thing without switching environments by hand. ``--pdb-output`` is the flag that matters. The ``.pqr`` file is PDB2PQR's normal output and Prodes cannot read it; ``--pdb-output`` writes the repaired structure as a PDB, which Prodes can. Run everything after this point, PROPKA included, on the repaired file. +``--keep-chain`` keeps the chain identifier in the ``.pqr`` file, which PDB2PQR otherwise leaves blank. It makes no difference to the file Prodes reads, because ``--pdb-output`` always writes chains, but chain identifiers matter to Prodes and the flag costs nothing, so use it. See `chain identifiers`_. + **Keep the original file name and change the directory, not the other way round.** Prodes takes the ``ID`` column of the output from the file name it was given (see `Input files`_), so a repaired file called ``1GDW_prep.pdb`` labels that row ``1GDW_prep``. Over a dataset that is tedious to undo. Prodes will also print ``Ignoring unrecognized record 'TER'`` on a PDB2PQR output file. That is harmless: ``TER`` is a chain separator and carries no coordinates. @@ -512,6 +514,15 @@ Traps * **mmCIF input is unreliable.** PDB2PQR 3.6.1 accepts a ``.cif`` but produced an empty ``.pqr`` from a valid one in testing. Convert to PDB first. * **Be consistent within a dataset**, exactly as for PROPKA. Prepare all of your structures or none of them. This matters more than it looks: because structure predictors disagree about whether to write the C-terminal ``OXT``, a half-prepared dataset carries a systematic difference between its sources. See `the terminal oxygen `_. +Chain identifiers +~~~~~~~~~~~~~~~~~~ + +Prodes needs chain identifiers. It groups residues by chain, decides which cysteines are bonded into a disulfide using the chain and residue number together, and applies per-residue pKa values by residue number. A file whose chains have been flattened into one would change all three. + +**The file Prodes reads is safe.** ``--pdb-output`` always writes the chain identifier, whether or not ``--keep-chain`` is given: the two files are byte for byte identical. Checked on structures of two, four, six and eight chains, and on a six-chain structure the repaired file gives Prodes the same six chains, the same 727 residues and the same 24 disulfide bonds as the original. + +**PDB2PQR's own ``.pqr`` output is not.** There the chain column is blank unless ``--keep-chain`` is given. That file is what APBS and the tools built on it consume, so if you use the ``.pqr`` for anything, pass the flag. The documented commands pass it always, because remembering which of two output files preserves what is not a good use of anyone's attention. + Full detail, including what Prodes discards and why: `preparing a structure `_. pKa values and protonation states diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index a8173ca..3ae12be 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -91,7 +91,7 @@ Then: mkdir prepared conda activate pdb2pqr -pdb2pqr --ff=PARSE --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr +pdb2pqr --ff=PARSE --keep-chain --pdb-output=prepared/1GDW.pdb 1GDW.pdb prepared/1GDW.pqr conda activate prodes propka3 prepared/1GDW.pdb # the prepared file, not the download @@ -105,6 +105,8 @@ Inside a script, `conda run -n pdb2pqr pdb2pqr ...` calls PDB2PQR in its own env **Keep the original file name and change the directory.** Prodes takes the `ID` column from the file name, so a prepared file called `1GDW_prep.pdb` labels that row `1GDW_prep`. Over a dataset that is tedious to undo. +**`--keep-chain` is in the command on purpose.** It does not change the file Prodes reads, but it does change PDB2PQR's other output. See the section below. + **PROPKA writes its `.pka` into the directory you are standing in**, named after the input's base name, not next to the input file. That is why the third line above reads `1GDW.pka` and not `prepared/1GDW.pka`. Prodes will print `Ignoring unrecognized record 'TER'` on a PDB2PQR output file. That is harmless: `TER` is a chain separator and carries no coordinates. @@ -127,6 +129,23 @@ 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. +## 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 applies per-residue pKa values by residue number. Flattening the chains would change all three. + +PDB2PQR writes two files, and they behave differently. + +| output | chain identifier | +| --- | --- | +| `--pdb-output`, the file Prodes reads | **always written**, with or without `--keep-chain` | +| the `.pqr`, PDB2PQR's normal output | **blank** unless `--keep-chain` is given | + +The reason is in `io.print_biomolecule_atoms`: the PDB branch calls `atom.get_pdb_string()`, which takes no chain flag, while the PQR branch calls `atom.get_pqr_string(chainflag=chainflag)`. + +Verified rather than assumed. On structures of two, four, six and eight chains, the `--pdb-output` file carries every chain identifier, and running PDB2PQR with and without `--keep-chain` produces byte-for-byte identical PDB output. On the six-chain structure 1F6R, Prodes reads the same six chains, the same 727 residues and the same 24 disulfide bonds from the repaired file as from the original. + +So the flag is not required for this pipeline. It is in the documented command anyway, for two reasons: the `.pqr` is what APBS and the tools built on it consume, and nobody should have to remember which of two output files preserves what. + ## Traps **Do not use `--ffout`.** It renames residues into the force field's own scheme, putting `ASH`, `GLH`, `LYN` and `HID` into the file. Prodes raises `KeyError` on any of those. Without `--ffout` both PARSE and AMBER keep canonical residue names, verified at pH 4.0 and pH 7.4. From 4fe8fed4ed83ea4c3b5a6993bbf51f801b19b07c Mon Sep 17 00:00:00 2001 From: Mark Teese Date: Sat, 12 Sep 2026 21:56:46 +0200 Subject: [PATCH 8/8] Reconcile the structure preparation documentation with version 9.0 This branch was written against 7.1. Version 9.0 keys a predicted pKa by chain as well as residue number, which makes one note added here wrong and another understated. Removes the note warning that the pKa step prints hundreds of warnings on a multi-chain structure because predictions are keyed by residue number alone. That is no longer true, and the note sat a few lines above the section that now describes the fixed behaviour. The section on chain identifiers said Prodes applies pKa values by residue number. It looks them up by chain and residue number from 9.0, which is a stronger reason to keep chain identifiers through a structure repair, not a weaker one: a file whose chains were flattened would silently fall back to the behaviour 9.0 exists to remove. No change to the advice itself. The survey behind it still holds, because the line it measured in src/prodes/core/point.py is unchanged in 9.0. --- README.rst | 6 +----- docs/structure_preparation.md | 2 +- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/README.rst b/README.rst index d7c3960..1c78c53 100644 --- a/README.rst +++ b/README.rst @@ -517,7 +517,7 @@ Traps Chain identifiers ~~~~~~~~~~~~~~~~~~ -Prodes needs chain identifiers. It groups residues by chain, decides which cysteines are bonded into a disulfide using the chain and residue number together, and applies per-residue pKa values by residue number. A file whose chains have been flattened into one would change all three. +Prodes needs chain identifiers. It groups residues by chain, decides which cysteines are bonded into a disulfide using the chain and residue number together, and from version 9.0 looks up a predicted pKa by chain as well as residue number. A file whose chains have been flattened into one would change all three, and the pKa lookup would silently fall back to the behaviour that version 9.0 exists to remove. **The file Prodes reads is safe.** ``--pdb-output`` always writes the chain identifier, whether or not ``--keep-chain`` is given: the two files are byte for byte identical. Checked on structures of two, four, six and eight chains, and on a six-chain structure the repaired file gives Prodes the same six chains, the same 727 residues and the same 24 disulfide bonds as the original. @@ -600,10 +600,6 @@ Steps 1 and 2 are done **once per structure**. Step 3 can then be repeated as of Residues that appear in the file get the predicted value; every other residue keeps its default. So a prediction covering only the titratable residues, which is what these tools produce, is complete as far as Prodes is concerned. -.. note:: - - On a structure with more than one chain, such as an antibody, this step prints a long list of warnings of the form ``the pKa file gives THR 31 a ASP pKa, which that residue does not have; ignoring it``. Prodes currently keys pKa values by residue number alone, so every chain's predictions are offered to every chain's residue with that number and most of them are rejected. The warnings are noise rather than a problem with your file, but a genuine warning is easy to miss among them, and a residue that collides with the same residue type in another chain can take the wrong chain's value. Tracked as `issue #12 `_. - Cysteines and disulfide bonds ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/structure_preparation.md b/docs/structure_preparation.md index 3ae12be..9aac784 100644 --- a/docs/structure_preparation.md +++ b/docs/structure_preparation.md @@ -131,7 +131,7 @@ The four extra groups are the two rebuilt lysine side chains and the two C-termi ## 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 applies per-residue pKa values by residue number. Flattening the chains would change all three. +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. PDB2PQR writes two files, and they behave differently.