Prodes: Protein Descriptors

Prodes turns a 3D protein structure into numbers describing the surface properties. You give it a PDB file, and it gives you a row in a CSV file with 54 columns describing the surface of that protein: how much of it is charged, how much of it is hydrophobic, how those properties are spread out across the surface, and how the picture changes with pH.

Those columns are designed to be used directly as the input features of a machine learning or QSPR model. If you have measured something about a set of proteins (a retention time, an aggregation onset, a viscosity, a titre) and you want to predict it for proteins you have not measured, Prodes gives you the X side of that problem from structure alone.

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

What you can use it for

Prodes was originally built and validated for predicting retention times in anion-exchange and cation-exchange chromatography, and it is well tested there. But surface charge and surface hydrophobicity drive a great many things in protein science, so the same features have been applied to:

• Hydrophobic interaction chromatography — surface hydrophobicity profiling
• Aggregation propensity — identifying aggregation-prone surface regions
• Protein–surface interactions — non-specific binding to chromatography resins, filtration membranes and container surfaces
• Binding affinity — the electrostatic contribution to protein–ligand and protein–protein interfaces
• Developability of biologics — screening monoclonal antibodies and other therapeutic proteins for surface liabilities such as high surface hydrophobicity and charge asymmetry
• Protein stability — correlating surface properties with aggregation propensity and shelf life in liquid formulations
• Formulation development — predicting colloidal stability and viscosity behaviour from the surface charge distribution
• Biologics manufacturing — surface property screening for process development and purification design

Installation

Python Environment

We suggest setting up a specific python environment for Prodes, because Prodes pins specific versions of NumPy and pandas for increased reproducibility and stability. Installing it into the same environment as your other analysis scripts could change the NumPy version underneath them and break something that worked yesterday.

We suggest using conda to manage your environments. You can install via Miniforge, which is the free, no-strings conda distribution. mamba is a drop-in faster replacement; everywhere below you can type mamba instead of conda if you have it.

Prodes requires Python 3.13.

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

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.

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:

conda create -n pdb2pqr -c conda-forge pdb2pqr

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.

Keeping them apart means an upgrade to either is someone else's problem.

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

git clone https://github.com/datacatalysis/prodes.git
cd prodes
conda env create -n prodes -f environment.yml
conda env update -n prodes -f environment_dev.yml
conda activate prodes
pre-commit install

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:

pytest                      # run the test suite
pre-commit run --all-files  # lint and type-check the whole repository

Maintainers, strongly recommended. Secrets are scanned in CI by GitGuardian on every push, and blocked at the push itself by GitHub push protection. Anyone with push access should also catch one locally, before it ever leaves the machine. Install the opt-in hook once per clone, alongside the ordinary one:

pre-commit install
pre-commit install -c .pre-commit-config-ggshield.yaml --hook-type pre-push

Those are two different git hooks, .git/hooks/pre-commit and .git/hooks/pre-push, so they coexist. The second scans the commits being pushed and refuses the push if it finds a credential.

It needs a GitGuardian account. If you do not have one, set one up, then either export GITGUARDIAN_API_KEY or run ggshield auth login once:

pipx install ggshield
ggshield auth login           # on a server, see below

ggshield auth login opens a browser and waits for the OAuth callback on localhost, so it cannot work over SSH. On a server, create a personal access token in the GitGuardian dashboard with the scan scope and hand it to ggshield auth login --method token, which prompts for it and stores it. Exporting GITGUARDIAN_API_KEY skips the login entirely, which is what CI does. ggshield install -m global is an alternative that covers every repository you clone rather than this one.

This is deliberately not in .pre-commit-config.yaml, which everyone runs: ggshield needs an account, so a hook there would break the first commit of anyone who clones the repository. Outside contributors need only pre-commit install and are unaffected. See SECURITY.md.

Quick start

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:

mkdir prepared

conda activate pdb2pqr
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
python -m prodes.io.pka_converter 1GDW.pka propka -o 1GDW_pka.json  # convert them for Prodes
python -m prodes prepared/1GDW.pdb 1GDW.zip --ph 7.4 --pka 1GDW_pka.json  # calculate the features

That writes 1GDW.zip, a bundle holding the 54 features, the surface points they were calculated from, and ready-to-open viewer scripts. On a small protein the whole thing takes seconds.

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

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.

To see every option:

python -m prodes --help

The ones you are most likely to want:

• --ph — the pH at which to compute protonation states (default 7). Charge features change substantially with this, so set it to the pH of your buffer.
• --pka — supply per-residue pKa values from PROPKA or a similar tool. Recommended, as above; see pKa values and protonation states.
• --probe — the radius of the solvent probe used for the surface calculation (default 1.4 Å, i.e. water).
• --full-features — write the original 105-feature set instead of the reduced 54. See The reduced feature set.
• --n-workers, --chunksize, --mem-limit — CPU and memory tuning. You can ignore all three; the defaults are sensible. See Speed.

Viewing the surface

Prodes does not only return numbers. Every run also produces two ready-made views of the protein surface it measured, so you can see where the charge and the hydrophobic patches actually are.

Electrostatic potential. Red is negative, blue positive, white near zero.	Hydrophobicity. Grey is hydrophilic, pale green hydrophobic, dark green strongly so.

Unpack the bundle and open one of the PyMOL scripts from inside the directory:

unzip 1GDW.zip
cd 1GDW
pymol 1GDW_ep.pml                 # electrostatic potential
pymol 1GDW_hydrophobicity.pml     # hydrophobicity

The .pml files are PyMOL scripts: plain text files of PyMOL commands that load the structure and colour its surface in one step, so there is nothing to set up by hand. The .cxc files do the same for ChimeraX. The paths inside them are relative, so they only work when opened from inside the unpacked directory. If PyMOL is already open, use File, Run Script instead.

Reading the colours

• Electrostatic potential. The scale is taken from each protein's own range, so no patch is ever clipped. The limit used is written into the script, and comparing two proteins directly means setting the same limit on both.
• Hydrophobicity. The two green cutoffs are fixed rather than per-protein, so the same green means the same hydrophobicity on any structure and two proteins can be compared as they are.

Both views are drawn from the same surface points, so a patch in one lines up exactly with the same place in the other.

Showing the residues underneath

Both scripts already load the structure behind the surface, as a grey cartoon with the relevant side chains picked out: red and blue for acidic and basic in the potential view, forest green for the hydrophobic residues in the hydrophobicity view. It is hidden only because the cloud is opaque. To see which residues produce a patch, make the cloud see-through:

set sphere_transparency, 0.4, surface_ep        # potential view
set sphere_transparency, 0.4, hydrophobicity    # hydrophobicity view

Raise or lower the number to taste. Much above 0.5 and the surface colours start to wash out. Toggling the surface object off in the PyMOL object panel works too, and leaves the residues on their own.

The residues shown in green are those the selected hydrophobicity scale scores as hydrophobic, so they follow --hydro rather than being a fixed list. Note that histidine is only partly charged at pH 7, and that Prodes also places a charge at each chain terminus, so a charged patch may have no coloured side chain beneath it.

ChimeraX

The .cxc scripts are provided for ChimeraX users and are written to mirror the PyMOL ones, but they have not been tested. If one does not behave, the PyMOL scripts are the reference; please open an issue.

Input files

Prodes reads .pdb files and .pdb.zip archives holding exactly one structure; an archive is unpacked to a temporary directory and read from there.

The ID column of the output is the file name with its last extension removed, taken from the file you named rather than from anything inside it. So 1GDW.pdb and 1GDW.pdb.zip both give 1GDW, an archive called bar.pdb.zip gives bar whatever the member inside it is called, and a file called 1abc.ent.pdb gives 1abc.ent rather than 1abc.

Name your files the way you want your rows labelled.

One structure, one model, one conformation. Prodes describes a single molecule, and a PDB file can hold more than one. A file with several MODEL records, which is what an NMR structure is, is described by its first model; the count reaches prodes_run.json as models_in_file. A residue modelled in several alternate conformations is described by its best occupied one. Both are logged when they happen, and both are decisions the file did not make for you: see what Prodes calls one residue and alternate conformations. If you want the ensemble averaged rather than one member of it described, run Prodes over the models separately and average the CSVs yourself.

The output bundle

One run takes one structure and writes one zip bundle. Nothing is appended to anything, so runs are independent and safe to parallelise. The output path must end in .zip.

Unpacked, the bundle holds:

1GDW/
  1GDW_features.csv              the 54 features, one row, first column ID
  1GDW_surface_points.csv        every surface point: x, y, z, potential, hydrophobicity
  1GDW_ep.pml                    PyMOL, coloured by electrostatic potential
  1GDW_hydrophobicity.pml        PyMOL, coloured by hydrophobicity
  1GDW_ep.cxc                    ChimeraX, electrostatic potential
  1GDW_hydrophobicity.cxc        ChimeraX, hydrophobicity
  1GDW_ep.pdb                    the points, potential in the B-factor column
  1GDW_hydrophobicity.pdb        the points, hydrophobicity in the B-factor column
  1GDW.pdb                       the structure the run was given
  prodes_run.json                version, settings, time of the run, and what
                                 the structure was read as: disulfide bonds
                                 found, alternate conformations collapsed,
                                 models the file held, and elements guessed
                                 from a blank column
  README.txt                     the same explanation, inside the bundle

Both point clouds hold the same coordinates and differ only in the value carried in the B-factor column, so the two views describe exactly the same surface.

The point cloud and the features come out of the same calculation, so a figure can never disagree with a feature value.

To read the features back:

from prodes.output import read_features

from prodes.output import read_features, read_surface_points

features = read_features("1GDW.zip")
points = read_surface_points("1GDW.zip")     # x, y, z, ep_volts, hydrophobicity

What the features actually are

Prodes builds a dotted surface over the protein (a Shrake-Rupley solvent-accessible surface), assigns each surface point an electrostatic potential and a hydrophobicity, and then summarises those distributions. The 54 default columns are, broadly:

• Whole-molecule properties — molecular weight, total surface area, formal charge at your chosen pH, isoelectric point, dipole moment, and shape descriptors.
• Surface electrostatic potential — the maximum, minimum, mean and standard deviation of the potential over the whole surface, then the same statistics computed separately over just the positive and just the negative regions, plus a count of how many surface points are positive. A protein with a strongly positive area of surface and a strongly negative one can have a near-zero net charge but very distinctive values here, which is exactly the kind of thing a net-charge calculation misses. Note that these are statistics over all positive or all negative surface points; Prodes does not group them into individual patches (see Similar tools).
• Surface hydrophobicity — the same treatment applied to a molecular hydrophobicity potential mapped onto the surface.
• Far-field shell electrostatics — the potential projected out onto shells around the protein, which captures how the charge distribution looks to another molecule approaching from a distance rather than at contact.
• Per-residue surface fractions — how much of the accessible surface each amino acid type contributes.

Every feature has a plot label, a one-line description, a full explanation of how it is calculated, its unit, and the wording used for it in the original publication. These live in two human-readable YAML files, which are also the single source of truth for which features Prodes writes, so you can read the whole list without running anything:

• features_reduced.yaml — the 54 features calculated by default
• features_full_only.yaml — the other 51, written only under --full-features, each with the reason it is not in the default set

Looking features up from Python

The same information is available programmatically, which is what you want when labelling a plot or a feature-importance table. A feature code is the short identifier Prodes writes as the CSV column heading, such as SurfEpMeanFormal, and every lookup takes one and returns something more readable:

from prodes.feature_dictionary import FeatureDictionary

fd = FeatureDictionary()

fd.get_plot_name("SurfEpMeanFormal")        # 'Surface EP mean (formal)'
fd.get_description("SurfEpMeanFormal")      # one line, suitable for a table
fd.get_long_description("SurfEpMeanFormal") # how it is calculated
fd.get_unit("Molecular weight")             # 'Da'; None if dimensionless
fd.get_original_explanation("SurfEpMeanFormal")   # wording from the 2024 paper

fd.get_reason_dropped("SurfEpMeanAverage")  # why it is not in the default set
fd.get_reason_dropped("SurfEpMeanFormal")   # None: this one is kept

The lookups are one-directional by design: the feature code is the key, and there is no route back from a label or description to a code. fd.get_entry(code) returns everything known about one feature and fd.get_dictionary() returns all 105 entries at once.

Getting the feature codes

Both lists come from the YAML dictionaries shipped in prodes/data, so a pipeline can line up datasets calculated under either setting without re-running Prodes or comparing CSV headers by hand:

from prodes.feature_dictionary import FeatureDictionary

fd = FeatureDictionary()

fd.get_feature_codes()             # the 54 written by default
fd.get_feature_codes(full=True)    # all 105 legacy columns
fd.get_dropped_feature_codes()     # the 51 the default leaves out

Each returns a fresh list of feature codes, in the order Prodes writes the columns. ID is a row label rather than a feature and is in neither list; it is available as prodes.feature_dictionary.ID_COLUMN. The underlying tuples are exposed as FULL_FEATURE_CODES, REDUCED_FEATURE_CODES and DROPPED_FEATURE_CODES if you would rather not construct the class.

If you have proteins measured with the full set and others with the reduced set, drop_redundant_features cuts the full ones down so the two are directly comparable:

import pandas as pd
from prodes.feature_dictionary import FeatureDictionary

fd = FeatureDictionary()

combined = pd.concat([
    fd.drop_redundant_features(pd.read_csv("measured_with_full_features.csv")),
    pd.read_csv("measured_with_reduced_features.csv"),
])

It keeps the ID column in front by default (pass keep_id=False to drop it) and raises KeyError naming the absent columns if the input is missing any reduced feature, rather than silently returning a narrower frame.

The lists are checked against the real output of calculate() on every test run (tests/test_feature_dictionary.py), so they cannot drift from what the code produces.

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, repair each structure, predict pKa values for it, and let the feature rows accumulate in one CSV:

import subprocess
from pathlib import Path

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")):
    repaired = prepared / pdb.name
    subprocess.run(
        ["conda", "run", "-n", "pdb2pqr", "pdb2pqr", "--ff=PARSE", "--keep-chain",
         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))

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.

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

import pandas as pd

from prodes.output import read_features

X = pd.concat([read_features(bundle) for bundle in sorted(Path("bundles").glob("*.zip"))])
y = pd.read_csv("measurements.csv")       # columns: ID, retention_time

data = X.merge(y, on="ID", validate="one_to_one")

Using validate="one_to_one" is worth the extra keystrokes: it raises if a structure was calculated twice or a measurement is duplicated, which is the failure mode that quietly halves your effective sample size.

Step 3, fit a model. The features are on wildly different scales (a molecular weight in the tens of thousands next to a surface fraction between 0 and 1), so scale them:

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_score

feature_cols = [c for c in data.columns if c not in ("ID", "retention_time")]

model = make_pipeline(StandardScaler(), PLSRegression(n_components=5))
scores = cross_val_score(model, data[feature_cols], data["retention_time"],
                         cv=5, scoring="r2")

PLS is a reasonable first choice for this kind of data because the features are correlated with one another by construction, and PLS handles that gracefully. Random forests and gradient boosting also work well and need no scaling.

A note on sample size. 54 features is a lot if you have measured 20 proteins. With fewer observations than features, almost any model will fit the training data perfectly and predict nothing. The reduced feature set exists precisely to push that ratio in your favour; do not reach for --full-features to get more columns unless you have the observations to support them (see The reduced feature set). Always report a score against a holdout blind test set, never a training-set score.

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

The projected electrostatic potential is damped with distance to account for the mobile ions in a buffer:

V(point) = sum over charged atoms of   q / (4 pi eps0 eps_r d)  *  exp(-d / lambda)

lambda = 3.04 / sqrt(I)   Angstrom, the Debye screening length

I is the ionic strength in mol/L, set with --ionic-strength. The default is 0.15, roughly physiological, which gives a screening length of 7.9 Angstrom.

Without this damping every charged atom contributes in full to every surface point, including atoms 60 Angstrom away on the far side of the protein. On a net negative protein, and most soluble proteins are, that adds a large smooth negative offset to the whole surface at once. The local pattern survives underneath it, so rank based comparisons looked fine, but the whole distribution is pushed below zero and genuinely positive patches are reported as negative. With screening, a charge one screening length away keeps about a third of its contribution and one four lengths away keeps under two per cent, so a surface point describes its own neighbourhood.

python -m prodes 1GDW.pdb 1GDW.zip --ionic-strength 0.035  # a 20 mM sodium phosphate buffer, pH 7
python -m prodes 1GDW.pdb 1GDW.zip --ionic-strength 0      # no screening, the pre-5.0 physics

Ionic strength is not the same as the molarity printed on the bottle. For a 1:1 salt such as NaCl they coincide, but a multivalent buffer contributes more: 20 mM sodium phosphate at pH 7 has an ionic strength of about 0.035 mol/L, not 0.02. Ionic strength is 0.5 * sum(c_i * z_i^2) over every ion present.

Range. Agreement with a Poisson-Boltzmann reference was measured to be flat between screening lengths of 5 and 12 Angstrom, which is an ionic strength of roughly 0.06 to 0.37 mol/L. Values well outside that, and a 20 mM buffer is outside it, are extrapolation rather than something that has been checked.

Which value to use. The ionic strength of the buffer the protein is binding in, not the one it elutes at. In a gradient the protein binds at the starting buffer and elutes when the salt rises enough to compete, so the binding condition is a known constant of the experiment rather than the quantity being predicted.

Setting it to 0 gives the unscreened potential of versions before 5.0. That is pinned point by point in tests/test_screening.py against a Coulomb sum written out independently of the code under test.

The values are not quite identical to the released 4.x numbers, for a separate reason. The potential is now stored to three decimals rather than two. At two, a point whose potential fell between 0 and 0.005 rounded to zero, and since the positive count tests for greater than zero, that point silently left the positive population. Nothing ever rounded into it, so the error only ever subtracted, undercounting NSurfPosEp by up to about six per cent on the screened potential and around one per cent without it. Three decimals brings that under one per cent. The extra digit is not physically meaningful, but it keeps the threshold from eating real points.

What this is and is not. The screening length is the Debye length in water, where the relative permittivity is about 78.5. Prodes evaluates its sum at a uniform relative permittivity of 4, where the self consistent value would be about 1.8 Angstrom. Using the water value inside a kernel with a protein dielectric is an empirical correction that mimics screening. It is used because it is what was measured to agree best with a Poisson-Boltzmann reference, not because the two are consistent, and the potential is still not comparable to an APBS calculation despite both being reported in volts.

Scope. Screening applies to the surface features, the SurfEp family: 9 of the 54 default features, and 38 of the 105 with --full-features, the extra 19 being the SurfEp*Average columns computed from partial rather than formal charges. The ShellEp features, 9 by default and 19 with the full set, are unchanged, and that is a measured decision rather than an omission: they are computed by a different route which divides the path at the molecular surface and weights the solvent leg with a permittivity of 80 against 4 for the protein, so a distant charge is already damped about twentyfold and the offset that affected the surface never built up. Checked against an APBS equivalent, the unscreened shell agrees at a Spearman of 0.877, and adding screening moved that to 0.855. See docs/screening_validation.md.

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. Install it into an environment of its own rather than alongside Prodes; see PDB2PQR, in its own environment for why.

conda create -n pdb2pqr -c conda-forge pdb2pqr      # once

mkdir prepared
conda activate pdb2pqr
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.

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.

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

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

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

Every charge-related feature Prodes calculates depends on which residues are protonated at your chosen pH, and that depends on their pKa values.

Left to itself, Prodes uses one textbook pKa per residue type: every aspartate in the structure is assumed to titrate at the same pH as every other aspartate. That is not true in a real protein. A buried aspartate next to another carboxylate can be shifted by several pH units from an exposed one on the far side of the molecule, and at your working pH the two may not even carry the same charge.

PROPKA predicts a pKa for each individual residue from its structural environment, and Prodes reads those predictions. This is the recommended default. It is a small, free, MIT-licensed program from the Jensen group, it installs on Windows, macOS and Linux, and on a normal protein it finishes in seconds. There is very little reason not to use it.

How much does it matter? On 1GDW, feeding in PROPKA values changes 19 of the 54 features, moving the isoelectric point from 10.34 to 10.57 and the formal charge from +7 to +8. In a larger test across 819 AlphaFold structures it changed 22 features, the most affected being SurfEpMinFormal at only R² = 0.82 between the two runs (full report).

Important

Be consistent within a dataset. Use PROPKA values for all of your structures or for none of them. Mixing the two puts two different kinds of number in the same feature column, and any model you fit will partly be learning which structures you happened to run PROPKA on.

Installing PROPKA

It goes into the same environment as Prodes, and is already included if you built your environment from environment.yml. Otherwise:

conda activate prodes
conda install conda-forge::propka

or equivalently pip install propka. Check it with propka3 --version.

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

Step 1, predict the pKa values. PROPKA reads your PDB file and writes a .pka file:

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

python -m prodes.io.pka_converter 1GDW.pka propka -o 1GDW_pka.json

or from Python, which is what you want inside a pipeline

from prodes.io.pka_converter import convert_propka, write_json

write_json(convert_propka("1GDW.pka"), "1GDW_pka.json")

Step 3, pass the converted file to Prodes

python -m prodes 1GDW.pdb 1GDW.zip --ph 7.4 --pka 1GDW_pka.json

import prodes
prodes.run_prodes("1GDW.pdb", "1GDW.zip", ph=7.4, pkas_file="1GDW_pka.json")

Steps 1 and 2 are done once per structure. Step 3 can then be repeated as often as you like at different pH values against the same JSON file, which is the reason the prediction and the calculation are separate commands rather than one.

Note

--pka takes the converted JSON, not PROPKA's own output. Passing a raw .pka file straight to --pka fails with a JSONDecodeError.

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.

Cysteines and disulfide bonds

A cysteine bonded into a disulfide has no thiol proton and does not titrate. From version 6.0 Prodes detects those bonds, from SSBOND records where the file has them and from the SG-SG distance otherwise, and gives the two cysteines no pKa. Before that every CYS was titrated as a free thiol, which on lysozyme put the formal charge at pH 8.5 at -1 when it is +7.

The number of bonds found is printed at the start of the run and recorded in prodes_run.json. For the cutoff, how records and geometry are combined, and the cases that are still wrong, see cysteines and disulfide bonds.

Other pKa predictors

Prodes does not run any predictor itself and does not import any of them; it only reads their output. Besides PROPKA, converters are shipped for H++ and pypka.

The second positional argument on the command line selects the converter and is one of propka, hpp or pypka; the Python equivalents are convert_propka, convert_hpp and convert_pypka. convert_propka and convert_pypka return a {chain: {residue_number: [{identifier: pka}]}} mapping, where the identifier is the three-letter residue name, or N+ and C- for the termini, and chain matches a residue's own chain identifier. convert_hpp returns the same shape, but keyed under the single reserved chain "*" (prodes.io.pka_converter.ANY_CHAIN) rather than a real chain identifier, because H++'s own output does not say which chain a residue belongs to; its predictions are applied to every chain that has a residue of the given number and type. Anything else that can produce that mapping can be fed in the same way.

Chains, and files from before version 9.0

A predicted pKa is keyed by chain as well as residue number, so it is only offered to the residue it was predicted for. This matters most for an antibody: a heavy and a light chain both number from residue 1, and before version 9.0 a value predicted for one was offered to the other's same-numbered residue too, applied whenever the two residues happened to share a type.

A --pka file written before version 9.0 is a flat {residue_number: [{identifier: pka}]} mapping, with no chain in it at all. Prodes still reads one: it is applied to every chain with a matching residue, exactly as the whole file used to be applied, and a warning names the file and recommends reconverting it. Reconverting means re-running the same converter that produced it, with the current version of Prodes, over the same predictor output; the predictor itself does not need to be re-run.

Chain still does not distinguish an insertion code: two residues that differ only by an insertion code and share a chain, number and type still collide, because PROPKA's summary line has no insertion-code column. See what Prodes calls one residue.

Using Prodes from Python

The command line and the Python interface do the same work. To reproduce a command-line run:

import prodes
prodes.run_prodes("./tests/data/1GDW.pdb.zip", "example.zip")

The full signature is run_prodes(pdb_file, out_file, pkas_file=None, ph=7, r_probe=1.4, hydro_scale="mj_scaled", full_features=False, mem_limit_mb=None).

The lower-level pieces are importable too, if you want a single property rather than the whole feature set. Calculating just the surface area, for example:

from prodes.io.parser import PDBparser
from prodes.calculations import grid_wizard, sasa

structure = PDBparser().parse("./tests/data/1GDW.pdb.zip")
grid = grid_wizard.Grid(10)
grid.construct_cells(structure.heavy_atoms)
grid.fill_cells(structure.heavy_atoms)

sasa.shrake_rupley(grid)

print(structure.surface_area())

Warning

Call Prodes from one thread at a time within a process. Several processes each running one structure is fine and is the normal way to scale up; several threads in one process is not, and will return wrong values rather than raising an error. See parallelism and memory.

The reduced feature set

The 105 numeric features in the original Prodes output were highly redundant. They have been reduced to 54 features to lower the risk of overfitting and to cut calculation time. Nothing new is calculated and no kept feature changes value: the default output is a strict column subset of the original, in the original column order.

To calculate the full set anyway, use the --full-features flag, set PRODES_FULL_FEATURES=true in the environment, or pass full_features=True to prodes.run.calculate().

The full feature set is only recommended if:

• You have a large number (>200) of independent observations in your dataset
• You have a well-established feature reduction pipeline that includes removal of correlated features
• You are using PCA or another dimensionality reduction technique to reduce the feature space
• You are using algorithms that are resistant to overfitting

For the analysis behind the reduction, see docs/redundant_feature_analysis.md.

Speed

How long will my protein take?

Measured on a 16-core Linux server, full 105-feature set, no pKa file:

Structure	Residues	Default (8 cores)	One core
ARH96693	60	2.0 s	4.2 s
1GDW	130	3.0 s	9.9 s
ARH98503	410	10.2 s	51.4 s
1GPB	823	30.2 s	184.6 s

One core is what you get on Windows and macOS, and on Linux with --n-workers 1.

Cost grows faster than protein size. Over this range a power law fits t ∝ n^1.44 on one core (R² = 0.99) and t ∝ n^1.04 on eight (R² = 0.97). Extra cores do not change how the work grows, they only hide more of it behind more hardware. Four structures over a 14x size range cannot separate a power law from an exponential, so read these as a description of the measured range rather than a formula to extrapolate with.

Nothing here is measured above ~800 residues. The only evidence at larger sizes is an earlier study of 51 Boltz-2 multimers from 218 to 1788 residues, which found the same strongly convex growth. Its absolute times are obsolete by roughly a factor of 11, because they predate the charged atom fix and multiprocessing, but its shape remains the best available guide for large multimers. See docs/calculation_time_benchmark.md.

If a structure is taking too long, split it into individual chains or domains where that is biologically meaningful. Because the cost is superlinear, two halves are genuinely cheaper than one whole, which is more than any setting will buy you.

Speedup over the original prodes/tneijenhuis codebase

On the protein structure with PDB code 1GPB, the calculation drops from 5146 s (86 minutes) to 30 s: a 171x speedup. Of that, 28x comes from rewriting the three hotspots in NumPy, which needs no special hardware and applies on every platform, and a further 5.6x from spreading one protein across eight CPU cores, which happens automatically on Linux.

The full benchmark is in docs/benchmark/benchmark_summary.md.

Using more than one CPU core

On Linux, this is already switched on and uses half your logical cores. You do not have to do anything.

On Windows and macOS, Prodes always runs on one core. This is a platform limitation rather than a setting: the parallel version shares memory with its worker processes through fork, which those operating systems do not provide. Everything still works and the 28x vectorisation speedup still applies; only the extra 5.6x is unavailable.

If you want to control this, the only setting most people need is:

python -m prodes in.pdb out.zip --n-workers 4

There is one trap worth knowing about even if you read nothing else. If you are processing many proteins, the best throughput comes from running one Prodes process per protein, with PRODES_N_WORKERS=1 set so that each process stays on one core. Workers are cumulative, so ten processes at the default eight workers each would ask your machine for eighty cores.

Everything else — worker counts, chunk sizes, memory budgets, thread safety and the platform reasoning in full — is in docs/parallelism_and_memory.md.

Memory

The default budget is 2048 MB for a whole run, which is enough for structures well past 1000 residues and is divided among the workers rather than taken per worker. A run stays under 2 GB whatever the worker count.

Raise or lower it with --mem-limit, PRODES_MEM_LIMIT_MB, or the mem_limit_mb argument. Lowering it costs remarkably little time: the same 1GPB run took 127 s at 2048 MB and 117 s at 64 MB, while peak RAM fell from 923 MB to 327 MB. Details and sizing tables are in docs/parallelism_and_memory.md.

Similar tools

Prodes sits in the same space as a number of other tools that derive surface charge and surface hydrophobicity descriptors from a protein structure. Only tools that actually produce such descriptors are listed here; a final subsection covers two pieces of infrastructure that are often mistaken for alternatives. We have not benchmarked Prodes against any of them, so what follows describes what each tool is, where to get it and how it is licensed, rather than claiming an advantage over it.

Open source, any protein

• PEP-Patch / surface_analyses (Liedl Lab, Innsbruck) — Permissive MIT licence. Cuts the surface into discrete electrostatic and hydrophobic patches by finding connected components on a triangulated surface, and reports the area and main residue of each. Requires APBS and PDB2PQR. Hoerschinger et al., J. Chem. Inf. Model. 2023, DOI: 10.1021/acs.jcim.3c01490.
• Protein-Sol Patches (Warwicker lab, Manchester) — web server only; the downloadable package on that site is the sequence solubility tool, not the patch analysis. Colours the surface by FDPB electrostatic potential (fixed at pH 6.3, no pKa calculation) and by a non-polar/polar SASA ratio taken over a 13 Å sphere around each atom, and reports the most non-polar region against a benchmark Fab distribution. Hebditch & Warwicker, Sci. Rep. 2019, DOI: 10.1038/s41598-018-36950-8.
• Aggrescan3D (Ventura / Kmiecik) — MIT. Projects experimentally derived aggregation propensities onto a structure, with an optional coarse-grained dynamics mode and an automated solubilising-mutation search. Aggregation-specific rather than a general descriptor generator. Kuriata et al., Nucleic Acids Res. 2019, DOI: 10.1093/nar/gkz321.

Antibody and nanobody specific

• PROPERMAB (Regeneron) — academic use only; commercial use is explicitly prohibited by the licence. The closest tool to Prodes in intent: 9 sequence and 26 structure features fed to machine learning models, including DBSCAN-segmented patch areas and CDR-localised versions. Fv only, structures predicted internally with ABodyBuilder2. Li et al., mAbs 2025, 17:2474521.
• TAP (OPIG, Oxford) — web application, no source released. Five metrics with red/amber/green flags set by percentiles of the clinical-stage therapeutic distribution. Its PSH/PPC/PNC "patch" metrics are 1/r²-weighted sums over surface residue pairs within 7.5 Å, not segmented patches. Raybould et al., PNAS 2019, DOI: 10.1073/pnas.1810576116.
• TNP (OPIG, Oxford) — BSD-3-Clause. TAP rebuilt for nanobodies, adding CDR3 compactness. The only permissively licensed implementation of the PSH/PPC/PNC family. Gordon et al., Commun. Biol. 2026, DOI: 10.1038/s42003-026-09594-y.
• HPATCH / APBS surface descriptors (Park & Izadi, Genentech) — no code released, but the closest published analogue to what Prodes computes, and the most careful study of how sensitive these descriptors are to structure model, protonation and conformational sampling. They build a NanoShaper triangulated surface, map APBS potentials onto its vertices, push the values down to atoms and then residues, and integrate the positive and negative parts separately over the Fab, Fv or CDR to give APBS_pos, APBS_neg and APBS_sum. HPATCH is the hydrophobicity analogue: a residue scale averaged over neighbouring vertices within 10 Å, then the positive residues summed. Everything is averaged over a 5 ns accelerated-MD ensemble rather than taken from one structure. mAbs 2024, DOI: 10.1080/19420862.2024.2362788.

Commercial suites

In industry this category is dominated by three licensed packages, and the antibody developability literature benchmarks against the same three almost every time. All require a licence server, all are per-seat, and none is free for commercial use. All three also bundle a homology-modelling step, so the workflow looks sequence-in from the user's side even though the descriptors are structure-based; Prodes takes a PDB file and stops there.

• MOE (Chemical Computing Group, CCG) — the Protein Properties application, with the Protein Patch Analyzer and Protein Patch 2D Maps panels for patch segmentation and visual QC. Generates a LowModeMD conformational ensemble (with extra CDR loop sampling for antibodies) and averages the descriptors over it, which is where the ens_* and avg_cdr_* descriptor families come from. Roughly 250 descriptors, including dedicated HIC retention models. Structure preparation is handled by QuickPrep, which corrects the structure, forms disulfides and assigns protonation states with Protonate3D. Also available as the browser-based BioMOE.
• BioLuminate (Schrödinger) — the Protein Descriptors panel and the calc_protein_descriptors.py batch script, alongside the Protein Surface Analyzer for hydrophobic and charged surface patches and AggScore for aggregation-prone regions. Descriptor counts reported in the literature range from about 900 to 1600 depending on release and on how many pH values are sampled.
• Discovery Studio (BIOVIA, Dassault Systèmes) — the Calculate Protein Features protocol (structure-based) and Calculate Sequence Descriptors (sequence-based, the only genuinely sequence-only protocol of the three), plus named developability predictors: the Developability Index built on SAP, the spatial charge map SCM for viscosity, and solubility and pI calculators.

Two of the methods these suites ship are the conceptual ancestors of most of the open tools above. SAP (Chennamsetty et al., PNAS 2009, DOI: 10.1073/pnas.0904191106) established averaging a hydrophobicity scale over a sphere around each atom as the standard way to score a protein surface, and AggScore (Sankar et al., Proteins 2018, DOI: 10.1002/prot.25594) extended it to the distribution of hydrophobic and charged patches. Neither has a canonical free implementation. Prodes' molecular hydrophobicity potential — an exp(-d)-weighted sum over non-hydrogen atoms within 10 Å — belongs to the same family.

Related, but not alternatives

These two come up constantly in this literature, but neither competes with Prodes: one is a solver several of the tools above depend on, and the other answers a different question.

• APBS + PDB2PQR (Electrostatics Consortium) — BSD-3-Clause. A Poisson-Boltzmann solver, not a descriptor tool. It produces a 3D electrostatic potential grid, and PEP-Patch, PROPERMAB and the Genentech pipeline all build on it. Jurrus et al., Protein Sci. 2018, DOI: 10.1002/pro.3280.

• NanoShaper / NanoShaperWeb (Decherchi & Rocchia, IIT) — GPL-3.0, so fine to run as a standalone binary but copyleft if you link it into a product you distribute. Two distinct roles, and it is worth keeping them apart:

  • As a surface engine, it triangulates the molecular surface analytically by ray casting rather than on a grid, and this is the mesh layer beneath both PROPERMAB and Genentech's descriptors. When those papers report a patch area in Ų, the number rests on a NanoShaper triangulation. Prodes does not need it, because it uses a dotted Shrake-Rupley solvent-accessible surface — a point cloud with no per-point area, which is why Prodes reports point counts where mesh-based tools report areas.
  • As a descriptor generator, it is aimed at pockets and cavities rather than the outer surface. Pockets are found from the volumetric difference between two solvent-excluded surfaces built with different probe radii, and NanoShaperWeb then characterises each pocket with the DrugPred druggability descriptor set — volume, entrance area, compactness, donor/acceptor/hydrophobic surface fractions. That is a small-molecule binding-site question — is this cavity druggable? — not the whole-surface charge and hydrophobicity profiling that Prodes and every other tool listed above is built for. The two are not substitutes in either direction.

  Decherchi & Rocchia, PLoS ONE 2013, DOI: 10.1371/journal.pone.0059744; Abate et al., J. Chem. Inf. Model. 2025, DOI: 10.1021/acs.jcim.5c00821.

Where Prodes sits

Currently, Prodes does not identify discrete patches. It builds a dotted surface, gives every point an electrostatic potential and a hydrophobicity, and then summarises those values over the whole surface and over the positive and negative subsets. Those subsets are defined by the sign of the value at a point, not by spatial contiguity: a positive-subset statistic pools every positive surface point on the protein, whether they form one large region or twenty scattered ones.

So if you need the area of the single largest hydrophobic patch, or the list of residues making up a particular patch so you can map it back onto the structure, use PEP-Patch. Prodes cannot give you that.

It is worth being clear about what this does not mean, though, because the word "patch" is used very loosely in this field. Prodes is not blind to where things sit on the surface. Its molecular hydrophobicity potential gives each surface point an exp(-d)-weighted sum over the non-hydrogen atoms within 10 Å, which is the same construction as SAP, as Genentech's HPATCH (a 10 Å average) and as Protein-Sol's NPP ratio (a 13 Å ratio). Only three of the tools above actually segment the surface into discrete regions: PEP-Patch (connected components on the surface graph), PROPERMAB (DBSCAN on mesh triangles) and MOE's Protein Patch Analyzer. The rest, TAP and TNP included, compute neighbourhood-weighted quantities and then aggregate them, which is what Prodes does too. The difference is in the last step: they take a maximum or a region sum, Prodes takes the distribution statistics.

What Prodes offers, then, is:

• Any protein. Most of the field is antibody-specific, often Fv-specific and tied to IMGT numbering.
• A fixed-width feature table, designed as the X matrix of a QSPR or machine learning model. Of the tools above only PROPERMAB shares that intent, and it is academic-only and antibody-only.
• pH and per-residue pKa handled properly, via PROPKA at a pH you choose. Protein-Sol Patches is fixed at pH 6.3 with no titration, PEP-Patch does not titrate, PROPERMAB is fixed at pH 7.4, and TAP assigns charges by residue type.
• No external solver. Nothing else here computes its own electrostatics. The cost is a much cruder physical model — a distance-weighted Coulomb sum at fixed permittivity, damped for ionic strength, rather than a Poisson-Boltzmann solution with an explicit dielectric boundary. For ranking a set of related proteins in a regression this is usually adequate; for absolute accuracy of the potential it is not, and APBS is the right tool.
• The MIT licence, and no licence server. For commercial work this rules out PROPERMAB, SAP and AggScore outright, puts MOE, BioLuminate and Discovery Studio behind a per-seat purchase, and makes anything built on NanoShaper awkward to redistribute. PEP-Patch, TNP and Aggrescan3D are the permissively licensed peers.

Three things Prodes does not do, which the better tools here do:

• No conformational averaging. Genentech's main finding is that descriptors taken from a single static structure are unstable, which is why they average over a 5 ns accelerated-MD ensemble; MOE averages over a LowModeMD ensemble, and Aggrescan3D offers a coarse-grained dynamic mode. Prodes computes from whatever one structure you give it, and inherits that instability. You can approximate the ensemble average by running Prodes over several structures and averaging yourself, but nothing in the package does it for you.
• Almost no structure preparation of its own. The commercial suites correct the input before they describe it — resolve alternate conformations, rebuild missing side chains, cap chain breaks, form disulfide bonds, optimise the hydrogen-bond network. Prodes does two of those and none of the rest. It recognises existing disulfide bonds, so that a cystine is not titrated as a free thiol, and it resolves alternate conformations by keeping the best occupied one per residue, which it does more reliably than PDB2PQR. It does not form bonds, rebuild side chains, cap chain breaks or touch any coordinate. The free fix is to run PDB2PQR first, which is why that is now the recommended workflow; see Preparing your structure.
• No region restriction, and no pockets. Every Prodes feature covers the whole molecule. There is no way to compute over a CDR, a domain or an interface, and nothing equivalent to NanoShaper's cavity and pocket detection.

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.

Improvements in this fork

• Vectorised Shrake-Rupley SASA — the per-atom sphere point / neighbour distance test is now computed via NumPy broadcasting instead of Python loops (~8x speedup on the SASA phase).
• Vectorised shell feature computation — find_exit, project_point, and map_ep_to_plane have batch counterparts (_batch) that process all charged atoms simultaneously with np.einsum and broadcasting (~10x speedup on the shell phase, previously the slowest part of the pipeline).
• Vectorised surface grid construction — grid construction and cell filling refactored to use NumPy arrays throughout.
• Multi-core parallelism — the SASA, surface grid and shell phases are spread across worker processes on Linux.
• Reduced feature set — 54 non-redundant features by default, with the original 105 available via --full-features.
• Disulfide bonds recognised — a cysteine bonded into a disulfide is no longer titrated as a free thiol, which corrected the charge-derived features of every structure with a disulfide. See cysteines and disulfide bonds.
• Alternate conformations resolved — the altLoc column used to be read as part of the atom name, so a disordered residue kept every conformation in the surface and, because the package looks up charge by atom name, silently carried no charge at all. From version 7.0 one conformation per residue is kept and the atom name is read on its own. See alternate conformations.
• PDB records read by Biopython — from version 7.1 the coordinate records of a PDB file are read by Bio.PDB.PDBParser, driven with a record-collecting builder of prodes' own rather than through Biopython's entity tree, which represents disorder rather than resolving it. No feature value moves on any well-formed structure. The insertion code and the model number now arrive as their own fields, which is what the two parser defects still open at that release needed; both were fixed in 8.0. See how prodes reads a PDB file.
• Insertion codes and NMR models — from version 8.0 a residue is its chain, its number and its insertion code, and a file holding several models is described by its first. Kabat and Chothia numbering write every antibody CDR insertion as an insertion code, and reading H100 and H100A as one residue put both side chains in one residue: on the 4NZU Fab that lost 12 of its 434 residues, moved the molecular weight by 1.3 kDa, titrated a lysine and a cysteine against aspartate's pKa, and hid a cysteine well enough that one of its five disulfide bonds went undetected. An NMR ensemble read whole was worse: 1PIT came out as 58 residues holding 17 780 atoms, 16 901 of them in a single residue, at a formal charge of -1096 for a protein whose real charge at pH 7 is about +6. See what Prodes calls one residue.
• A blank element column is inferred rather than fatal — from version 8.1 an atom whose element column is blank has its element guessed from its name, the same convention Bio.PDB.Atom uses. A blank element used to reach the van der Waals radius lookup and raise KeyError, so a PDB file with no element column produced no features at all; no shipped structure has one, so no feature value moves. The count of atoms inferred is logged once per file and recorded in prodes_run.json as inferred_elements. See how Prodes reads a PDB file.
• pKa predictions are keyed by chain — from version 9.0 a predicted pKa is looked up by chain as well as residue number, rather than by number alone. A structure with more than one chain, an antibody's heavy and light chain included, could have a value predicted for one chain applied to another chain's same-numbered residue instead, silently whenever the two residues shared a type. convert_propka and convert_pypka now read the chain out of their source files; convert_hpp cannot, because H++'s own output carries no chain, and marks its predictions accordingly. A --pka file written before this version has no chain in it at all and is still read, applied to every chain exactly as before, with a warning recommending it be reconverted. See pKa values and protonation states.
• Bug fixes — trimean edge case for small arrays, surface_exit None handling in shell potential mapping, and a read_propka argument bug in the PDB parser.
• Test suite — unit and regression tests added, including a committed reference output file (tests/data/ARH96693_prodes_orig_output.csv) generated by the original unrefactored code. The regression test verifies that every feature column and every feature value produced by this fork matches the original output within tolerance.

Output compatibility

With --full-features, the output has the same 105 columns in the same order as the original Prodes.

The values are identical only with --ionic-strength 0, and then only for a structure with no disulfide bonds and no alternate conformations. From version 9.0 a predicted pKa is applied by chain as well as residue number; this only changes a structure's features if its --pka file is reconverted with the current convert_propka or convert_pypka, and then only if the structure has more than one chain and a chain-specific value lands on a residue that previously took another chain's value or none at all: on the 4NZU Fab, reconverting the same PROPKA run moves 13 of the 105 columns, the isoelectric point among them, from 6.796 to 6.808. A --pka file that predates this version, or one produced by the current convert_hpp, is still applied to every chain exactly as before and changes nothing. See pKa values and protonation states. From version 8.0 a residue that carries an insertion code is its own residue rather than part of the one before it, which moves the features of any structure that has one, antibodies numbered by the Kabat or Chothia schemes included: on the 4NZU Fab 69 of the 105 columns move and 36 do not. Nothing geometric moves, because no atom moves and no atom changes its own residue name: Area, NSurfPoints, both Shape columns and the whole hydrophobicity block are unchanged, while everything that reads a residue's type, mass or charge does move. A file holding more than one model is described by its first, which changes every feature of an NMR ensemble, all of them from nonsense to a number. Structures with neither an insertion code nor a second model are unaffected. See what Prodes calls one residue. From version 7.0 only one conformation of a disordered residue is described, which changes the features of any structure that has one; the size of the change tracks how many titratable residues are disordered rather than how many residues are, so a structure whose alternates sit on uncharged residues barely moves. AlphaFold models contain no alternate conformations and are unaffected. See alternate conformations. From version 6.0 a cysteine bonded into a disulfide is not titrated, which changes the charge-derived features of any structure that has one. See cysteines and disulfide bonds. From version 5.0 the electrostatic potential is screened by default, which changes the 38 SurfEp columns; the other 67, including every ShellEp column, are unaffected either way. The regression test in tests/test_sasa.py therefore runs at zero ionic strength when comparing against the committed reference file tests/data/ARH96693_prodes_orig_output.csv. See Ionic strength and screening.

The default reduced output is a strict subset of those columns, in the same order and with identical values. No column is renamed and no column is recomputed, so a reduced run and a full run of the same structure agree exactly on the 54 columns they share.

How to cite

If this package is useful for you, please cite the original Prodes publication:

Neijenhuis, T., Le Bussy, O., Geldhof, G., Klijn, M. E., & Ottens, M. (2024). Predicting protein retention in ion-exchange chromatography using an open source QSPR workflow. Biotechnology Journal, 19, e2300708. https://doi.org/10.1002/biot.202300708

Contact, Maintenance, and Improvements

Contributions are welcome. To report a security problem rather than a bug, see SECURITY.md, which asks you not to open a public issue for it.

Currently the code is maintained by Mark Teese of 22DataCatalysis GmbH. Please raise a GitHub issue or contact us via the contact page on our website if you encounter any problems or have suggestions for improvements.