Evaluate ligand poses and protein interfaces

Compare predictions with PLINDER references using OpenStructure, and check ligand plausibility with PoseBusters.

Evaluate two Boltz predictions for the 8c3u protein–ligand complex from Runs N’ Poses and a GeoDock prediction for the 2e31 protein interface from PINDER’s evaluation tutorial.

Set up

Install OpenStructure 2.12.0 or newer and the evaluation extra in your notebook environment:

conda install -c conda-forge -c bioconda 'openstructure>=2.12.0'
pip install 'plinder[eval]'

The ost command must be on your PATH. Set PLINDER_DATA_DIR to use a local release; otherwise PLINDER uses the configured release. The evaluation reference describes the options and output files.

import os
import tempfile
from pathlib import Path
from urllib.request import urlretrieve

import yaml
import pandas as pd

from plinder.core import PlinderRelease
from plinder.eval import evaluate_predictions

data_dir = os.environ.get("PLINDER_DATA_DIR")
release = PlinderRelease(data_dir=Path(data_dir)) if data_dir else PlinderRelease()
work_dir = Path(
    os.environ.get("PLINDER_EVAL_WORK_DIR")
    or tempfile.mkdtemp(prefix="plinder-evaluation-")
).resolve()
predictions = work_dir / "predictions"
work_dir

/home/runner/micromamba/envs/plinder/lib/python3.12/site-packages/vcs_versioning/_backends/_git.py:709: UserWarning: "/home/runner/work/plinder/plinder" is shallow and may cause errors
  pre_parse(wd)

PosixPath('/tmp/plinder-evaluation-o9tu1crr')

Choose prediction files

Search and evaluation share an input table with input_id, structure_path, and optional ligand_path. Evaluation adds a reference_id: a ligand-system ID, an interface ID, or a PDB ID. The evaluation reference describes the columns and folder alternative.

predictions/
├── <ligand-system ID>/
│   ├── boltz_seed_1372115236.cif
│   └── boltz_seed_555476360.cif
└── <interface ID>/
    └── geodock_model_1.pdb

For your own models, use the same folder layout and skip this example-data setup. Interface evaluation accepts PDB or mmCIF files. For ligand evaluation, pair a receptor model.pdb or model.cif with model.sdf, or a model.ligands/ folder containing one SDF per ligand. Each prediction then carries its own ligand chemistry and coordinates.

Complete-mmCIF predictions are also supported, as shown with the two Boltz models below. Their input YAML provides the ligand SMILES.

system_id = "8c3u__1__1.A__1.C"
revision = "78bfe19e15b086e39de4e03f30142e0dcb38beb0"
example_url = f"https://raw.githubusercontent.com/plinder-org/runs-n-poses/{revision}/examples"
ligand_folder = predictions / system_id
ligand_folder.mkdir(parents=True, exist_ok=True)

for seed in (1372115236, 555476360):
    model = ligand_folder / f"boltz_seed_{seed}.cif"
    if not model.is_file():
        urlretrieve(
            f"{example_url}/outputs/boltz/{system_id}/{seed}/"
            "boltz_results_input/predictions/input/input_model_0.cif",
            model,
        )

input_yaml = work_dir / "boltz_input.yaml"
if not input_yaml.is_file():
    urlretrieve(f"{example_url}/inputs/boltz/{system_id}/input.yaml", input_yaml)
boltz_input = yaml.safe_load(input_yaml.read_text())
ligand_smiles = next(
    item["ligand"]["smiles"] for item in boltz_input["sequences"] if "ligand" in item
)

interface_id = "2e31__1__1.A--1.B"
interface_folder = predictions / interface_id
interface_folder.mkdir(parents=True, exist_ok=True)

pinder_revision = "20487ac403516f453d8be1b0bdc07a9f7a77d9df"
pinder_id = "2e31__A1_Q80UW2--2e31__B1_P63208"
geodock_pdb = interface_folder / "geodock_model_1.pdb"
if not geodock_pdb.is_file():
    urlretrieve(
        f"https://raw.githubusercontent.com/pinder-org/pinder/{pinder_revision}/"
        f"tests/test_data/method_eval/geodock/{pinder_id}/holo_decoys/model_1.pdb",
        geodock_pdb,
    )

input_table = pd.DataFrame([
    {
        "input_id": path.stem,
        "structure_path": str(path.resolve()),
        "reference_id": path.parent.name,
    }
    for path in sorted(predictions.glob("*/*")) if path.is_file()
])
input_table

	input_id	structure_path	reference_id
0	geodock_model_1	/tmp/plinder-evaluation-o9tu1crr/predictions/2...	2e31__1__1.A--1.B
1	boltz_seed_1372115236	/tmp/plinder-evaluation-o9tu1crr/predictions/8...	8c3u__1__1.A__1.C
2	boltz_seed_555476360	/tmp/plinder-evaluation-o9tu1crr/predictions/8...	8c3u__1__1.A__1.C

Evaluate the predictions

One call evaluates both kinds of prediction. Choose mode="ligands" or mode="interfaces" for just one kind. num_workers limits the number of predictions processed at once.

By default, ligand results cover proper reference ligands. Set include_all_ligands=True to include reference ions and artifacts too.

The two Boltz predictions share the same ligand, so pass its SMILES as ligand_smiles={"LIG": ligand_smiles}. A SMILES must follow the CIF heavy-atom order. With complete-mmCIF inputs, ligand_ccd_codes supplies CCD templates and ligand_chains=["B"] selects a peptide ligand by its label asym ID. For receptor + SDF inputs, the SDF files supply the chemistry.

results = evaluate_predictions(
    input_table,
    output_dir=work_dir / "results",
    release=release,
    mode="both",
    num_workers=2,
    ligand_smiles={"LIG": ligand_smiles},
)

display(results["failures"])
assert results["failures"].empty, "Inspect the errors and tool logs before interpreting results."

/home/runner/micromamba/envs/plinder/lib/python3.12/site-packages/vcs_versioning/_backends/_git.py:709: UserWarning: "/home/runner/work/plinder/plinder" is shallow and may cause errors
  pre_parse(wd)
/home/runner/micromamba/envs/plinder/lib/python3.12/site-packages/vcs_versioning/_backends/_git.py:709: UserWarning: "/home/runner/work/plinder/plinder" is shallow and may cause errors
  pre_parse(wd)

	prediction	system_id	stage	error

Read ligand results

Each row represents one reference ligand and one prediction. rmsd is the binding-site-superposed, symmetry-corrected ligand RMSD; bb_rmsd describes the backbone fit, lddt_lp the local binding-site structure, and lddt_pli the protein–ligand contacts. Lower RMSD is better; lDDT values range from zero to one.

status="success" means the comparison completed. Both Boltz examples here have a large ligand RMSD. A missing or unassigned reference ligand has a row with unavailable metrics and an unassigned reason.

ligands = results["ligands"]
ligands[
    [
        "prediction", "ligand_id", "status",
        "rmsd", "bb_rmsd", "lddt_lp", "lddt_pli",
        "rmsd_unassigned", "lddt_pli_unassigned",
    ]
]

	prediction	ligand_id	status	rmsd	bb_rmsd	lddt_lp	lddt_pli	rmsd_unassigned	lddt_pli_unassigned
0	boltz_seed_1372115236	8c3u__1__1.C	success	21.360781	2.621234	0.639289	0.000949	NaN	NaN
1	boltz_seed_555476360	8c3u__1__1.C	success	21.320765	2.660444	0.638495	0.000000	NaN	NaN

Attach PoseBusters checks

PoseBusters checks all candidate predicted ligands. To attach those checks to a reference-ligand result, use the predicted ligand selected by the metric you are analysing. RMSD and lDDT-PLI can choose different assignments.

Below, the join follows rmsd_model_ligand. For lDDT-PLI, use lddt_pli_model_ligand instead. PoseBusters uses each prediction’s own receptor coordinates. Failed checks appear as False; PoseBusters fields for missing ligands are empty.

checks = results["posebusters"][
    [
        "prediction", "model_ligand", "sanitization",
        "internal_steric_clash", "minimum_distance_to_protein",
    ]
]
ligands_with_checks = ligands.merge(
    checks,
    left_on=["prediction", "rmsd_model_ligand"],
    right_on=["prediction", "model_ligand"],
    how="left",
    validate="many_to_one",
)
ligands_with_checks[
    [
        "prediction", "ligand_id", "rmsd_model_ligand", "rmsd",
        "sanitization", "internal_steric_clash", "minimum_distance_to_protein",
    ]
]

	prediction	ligand_id	rmsd_model_ligand	rmsd	sanitization	internal_steric_clash	minimum_distance_to_protein
0	boltz_seed_1372115236	8c3u__1__1.C	B	21.360781	True	True	False
1	boltz_seed_555476360	8c3u__1__1.C	B	21.320765	True	True	False

Read interface results

OST maps the GeoDock model’s chains R and L to reference chains 1.A and 1.B.

ilddt describes local agreement at the interface; QS-score describes assembly agreement. Use dockq_ave_full or dockq_wave_full when comparing predictions that may have missing chains: these include the penalty for unmapped reference interfaces.

results["interfaces"][
    [
        "prediction", "system_id", "status", "lddt", "ilddt",
        "qs_global", "qs_best", "dockq_ave_full", "dockq_wave_full",
    ]
]

	prediction	system_id	status	lddt	ilddt	qs_global	qs_best	dockq_ave_full	dockq_wave_full
0	geodock_model_1	2e31__1__1.A--1.B	success	0.216	0.029	0.187	0.187	0.066	0.066

Access reference-ligand assignments in results["ligands"] using rmsd_model_ligand and lddt_pli_model_ligand, with reasons for unmatched ligands in rmsd_unassigned and lddt_pli_unassigned.

Reload saved tables with pandas, for example pd.read_parquet(work_dir / "results" / "ligands.parquet").