Analysis of a network of RBFE SepTop calculations#
In this notebook we show how to analyze a network of transformations, run with the OpenFE Separated Topologies protocol. This notebook shows you how to extract
The overall difference in binding affinity between two ligands (DDG)
The MLE-derived absolute binding free energies
The contribution from the different legs (complex and solvent) of a transformation.
Downloading the example dataset#
First let’s download some example SepTop results. Please skip this section if you have already done this!
[1]:
!wget https://zenodo.org/records/17435569/files/septop_results.zip
!unzip septop_results.zip
--2025-10-14 20:00:18-- https://zenodo.org/records/17435569/files/septop_results.zip
Resolving zenodo.org (zenodo.org)... 188.185.45.92, 188.185.48.194, 188.185.43.25
Connecting to zenodo.org (zenodo.org)|188.185.45.92|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 5974197 (5.7M) [application/octet-stream]
Saving to: ‘septop_results.zip’
septop_results.zip 100%[===================>] 5.70M 559KB/s in 11s
2025-10-14 20:00:30 (529 KB/s) - ‘septop_results.zip’ saved [5974197/5974197]
Archive: septop_results.zip
creating: septop_results
creating: septop_results/results_2
inflating: septop_results/results_2/rbfe_7a_7b.json
inflating: septop_results/results_2/rbfe_1_7b.json
inflating: septop_results/results_2/rbfe_1_7a.json
creating: septop_results/results_1
inflating: septop_results/results_1/rbfe_7a_7b.json
inflating: septop_results/results_1/rbfe_1_7a.json
creating: septop_results/results_0
inflating: septop_results/results_0/rbfe_7a_7b.json
inflating: septop_results/results_0/rbfe_1_7b.json
inflating: septop_results/results_0/rbfe_1_7a.json
inflating: septop_results/results_0/rbfe_1_25.json
Imports#
Here are a bunch of imports we will need later in the notebook.
[3]:
import numpy as np
import os
import pathlib
from gufe.tokenization import JSON_HANDLER
import pandas as pd
from openff.units import unit
from openfecli.commands.gather import (
format_estimate_uncertainty,
_collect_result_jsons,
load_json,
)
from cinnabar import Measurement, FEMap
Some helper methods to load and format the SepTop results#
Over the next few cells, we define some helper methods that we will use to load and format the SepTop results.
Note: you do not need to directly interact with any of these, unless you are looking to change the behaviour of how data is being processed
[3]:
def _load_valid_result_json(
fpath: os.PathLike | str,
) -> tuple[tuple | None, dict | None]:
"""Load the data from a results JSON into a dict.
Parameters
----------
fpath : os.PathLike | str
The path to deserialized results.
Returns
-------
dict | None
A dict containing data from the results JSON,
or None if the JSON file is invalid or missing.
Raises
------
ValueError
If the JSON file contains an ``estimate`` or ``uncertainty`` key with the
value ``None``.
If
"""
# TODO: only load this once during collection, then pass namedtuple(fname, dict) into this function
# for now though, it's not the bottleneck on performance
result = load_json(fpath)
try:
names = _get_names(result)
except (ValueError, IndexError):
print(f"{fpath}: Missing ligand names. Skipping.")
return None, None
if result["estimate"] is None:
errormsg = f"{fpath}: No 'estimate' found, assuming to be a failed simulation."
raise ValueError(errormsg)
return names, result
[4]:
def _get_legs_from_result_jsons(
result_fns: list[pathlib.Path],
) -> dict[str, dict[str, list]]:
"""
Iterate over a list of result JSONs and populate a dict of dicts with all data needed
for results processing.
Parameters
----------
result_fns : list[pathlib.Path]
List of filepaths containing results formatted as JSON.
report : Literal["dg", "ddg", "raw"]
Type of report to generate.
Returns
-------
legs: dict[str, dict[str, list]]
Data extracted from the given result JSONs, organized by the leg's ligand names and simulation type.
"""
from collections import defaultdict
ddgs = defaultdict(lambda: defaultdict(list))
for result_fn in result_fns:
names, result = _load_valid_result_json(result_fn)
if names is None: # this means it couldn't find names and/or simtype
continue
ddgs[names]["overall"].append([result["estimate"], result["uncertainty"]])
proto_key = [
k
for k in result["unit_results"].keys()
if k.startswith("ProtocolUnitResult")
]
for p in proto_key:
if "unit_estimate" in result["unit_results"][p]["outputs"]:
simtype = result["unit_results"][p]["outputs"]["simtype"]
dg = result["unit_results"][p]["outputs"]["unit_estimate"]
dg_error = result["unit_results"][p]["outputs"]["unit_estimate_error"]
ddgs[names][simtype].append([dg, dg_error])
elif "standard_state_correction_A" in result["unit_results"][p]["outputs"]:
corr_A = result["unit_results"][p]["outputs"][
"standard_state_correction_A"
]
corr_B = result["unit_results"][p]["outputs"][
"standard_state_correction_B"
]
ddgs[names]["standard_state_correction_A"].append(
[corr_A, 0 * unit.kilocalorie_per_mole]
)
ddgs[names]["standard_state_correction_B"].append(
[corr_B, 0 * unit.kilocalorie_per_mole]
)
else:
continue
return ddgs
[5]:
def _get_names(result: dict) -> tuple[str, str]:
"""Get the ligand names from a unit's results data.
Parameters
----------
result : dict
A results dict.
Returns
-------
tuple[str, str]
Ligand names corresponding to the results.
"""
try:
nm = list(result["unit_results"].values())[0]["name"]
except KeyError:
raise ValueError("Failed to guess names")
toks = nm.split(",")
toks = toks[1].split()
return toks[1], toks[3]
[6]:
def _error_std(r):
"""
Calculate the error of the estimate as the std of the repeats
"""
return np.std([v[0].m for v in r["overall"]])
[7]:
def _error_mbar(r):
"""
Calculate the error of the estimate using the reported MBAR errors.
This also takes into account that repeats may have been run for this edge by using the average MBAR error
"""
complex_errors = [x[1].m for x in r["complex"]]
solvent_errors = [x[1].m for x in r["solvent"]]
return np.sqrt(np.mean(complex_errors) ** 2 + np.mean(solvent_errors) ** 2)
Methods to extract and manipulate the SepTop results#
The next four methods allow you to extract SepTop results (extract_results_dict) and then manipulate them to get different types of results.
These include:
generate_ddg: to get the ddG values.generate_dg_mle: to get the MLE-derived dG values.generate_dg_raw: to get the raw dG values for each individual legs in the SepTop transformation cycles.
[8]:
def extract_results_dict(
results_files: list[os.PathLike | str],
) -> dict[str, dict[str, list]]:
"""
Get a dictionary of SepTop results from a list of directories.
Parameters
----------
results_files : list[ps.PathLike | str]
A list of directors with SepTop result files to process.
Returns
-------
sim_results : dict[str, dict[str, list]]
Simulation results, organized by the leg's ligand names and simulation type.
"""
# find and filter result jsons
result_fns = _collect_result_jsons(results_files)
# pair legs of simulations together into dict of dicts
sim_results = _get_legs_from_result_jsons(result_fns)
return sim_results
[9]:
def generate_ddg(results_dict: dict[str, dict[str, list]]) -> pd.DataFrame:
"""Compute and write out DDG values for the given results.
Parameters
----------
results_dict : dict[str, dict[str, list]]
Dictionary of results created by ``extract_results_dict``.
Returns
-------
pd.DataFrame
A pandas DataFrame with the ddG results for each ligand pair.
"""
data = []
# check the type of error which should be used based on the number of repeats
repeats = {len(v["overall"]) for v in results_dict.values()}
error_func = _error_mbar if 1 in repeats else _error_std
for ligpair, results in sorted(results_dict.items()):
ddg = np.mean([v[0].m for v in results["overall"]])
error = error_func(results)
m, u = format_estimate_uncertainty(ddg, error, unc_prec=2)
data.append((ligpair[0], ligpair[1], m, u))
df = pd.DataFrame(
data,
columns=[
"ligand_i",
"ligand_j",
"DDG(i->j) (kcal/mol)",
"uncertainty (kcal/mol)",
],
)
return df
[10]:
def generate_dg_mle(results_dict: dict[str, dict[str, list]]) -> pd.DataFrame:
"""Compute and write out MLE-derived DG values for the given results.
Parameters
----------
results_dict : dict[str, dict[str, list]]
Dictionary of results created by ``extract_results_dict``.
Returns
-------
pd.DataFrame
A pandas DataFrame with the dG results for each ligand pair.
"""
DDGs = generate_ddg(results_dict)
fe_results = []
for inx, row in DDGs.iterrows():
ligA, ligB, DDGbind, bind_unc = row.tolist()
m = Measurement(
labelA=ligA,
labelB=ligB,
DG=DDGbind * unit.kilocalorie_per_mole,
uncertainty=bind_unc * unit.kilocalorie_per_mole,
computational=True,
)
fe_results.append(m)
# Feed into the FEMap object
femap = FEMap()
for entry in fe_results:
femap.add_measurement(entry)
femap.generate_absolute_values()
df = femap.get_absolute_dataframe()
df = df.iloc[:, :3]
df.rename({"label": "ligand"}, axis="columns", inplace=True)
return df
[11]:
def generate_dg_raw(results_dict: dict[str, dict[str, list]]) -> pd.DataFrame:
"""
Get all the transformation cycle legs found and their DG values.
Parameters
----------
results_dict : dict[str, dict[str, list]]
Dictionary of results created by ``extract_results_dict``.
Returns
-------
pd.DataFrame
A pandas DataFrame with the individual cycle leg dG results.
"""
data = []
for ligpair, results in sorted(results_dict.items()):
for simtype, repeats in sorted(results.items()):
if simtype != "overall":
for repeat in repeats:
m, u = format_estimate_uncertainty(
repeat[0].m, repeat[1].m, unc_prec=2
)
data.append((simtype, ligpair[0], ligpair[1], m, u))
df = pd.DataFrame(
data,
columns=[
"leg",
"ligand_i",
"ligand_j",
"DG(i->j) (kcal/mol)",
"uncertainty (kcal/mol)",
],
)
return df
Analyzing your results#
Now that we have defined a set of methods to help us extract results. Let’s analyze the results!
Specify result directories and gather results#
Let’s start by gathering all our simulation results. First we define all the directories where our SepTop results exist. Here we assume that our simulation repeats sit in three different results directories under septop_results, named from results_0 to results_2.
[12]:
# Specify paths to result directories
results_dir = [
pathlib.Path("septop_results/results_0"),
pathlib.Path("septop_results/results_1"),
pathlib.Path("septop_results/results_2"),
]
ddgs = extract_results_dict(results_dir)
/home/ialibay/software/mambaforge/install/envs/openfe/lib/python3.12/site-packages/openmoltools/utils.py:9: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
from pkg_resources import resource_filename
/home/ialibay/software/mambaforge/install/envs/openfe/lib/python3.12/site-packages/Bio/Application/__init__.py:39: BiopythonDeprecationWarning: The Bio.Application modules and modules relying on it have been deprecated.
Due to the on going maintenance burden of keeping command line application
wrappers up to date, we have decided to deprecate and eventually remove these
modules.
We instead now recommend building your command line and invoking it directly
with the subprocess module.
warnings.warn(
Obtain the overall difference in binding affinity for all edges in the network#
With these extracted results, we can now get the ddG prediction between each pair of ligand.
Note: if only a single repeat was run, the MBAR error is used as uncertainty estimate, while the standard deviation is used when results from more than one repeat are provided.
[13]:
df_ddg = generate_ddg(ddgs)
df_ddg.to_csv("ddg.tsv", sep="\t", lineterminator="\n", index=False)
[14]:
df_ddg
[14]:
| ligand_i | ligand_j | DDG(i->j) (kcal/mol) | uncertainty (kcal/mol) | |
|---|---|---|---|---|
| 0 | 1 | 25 | 2.0 | 1.6 |
| 1 | 1 | 7a | 0.6 | 1.5 |
| 2 | 1 | 7b | 0.1 | 1.5 |
| 3 | 7a | 7b | 1.9 | 1.5 |
Obtain the MLE-derived absolute binding affinities#
We can also get the estimated binding dG for each ligand.
This uses the maximum-likelihood method as implemented in cinnabar.
[15]:
df_dg = generate_dg_mle(ddgs)
df_dg.to_csv("dg.tsv", sep="\t", lineterminator="\n", index=False)
[16]:
df_dg
[16]:
| ligand | DG (kcal/mol) | uncertainty (kcal/mol) | |
|---|---|---|---|
| 0 | 1 | -0.675 | 0.664267 |
| 1 | 25 | 1.325 | 1.311964 |
| 2 | 7a | -0.875 | 0.903466 |
| 3 | 7b | 0.225 | 0.903466 |
Obtain the raw DGs of every leg in the thermodynamic cycle#
If needed, you can also get the individual dG results for each leg of the SepTop transformation cycles. This can be useful in various different situation, such as when trying to diagnose which part of your simulation has the highest uncertainty.
[17]:
df_raw = generate_dg_raw(ddgs)
df_raw.to_csv("ddg_raw.tsv", sep="\t", lineterminator="\n", index=False)
[18]:
df_raw
[18]:
| leg | ligand_i | ligand_j | DG(i->j) (kcal/mol) | uncertainty (kcal/mol) | |
|---|---|---|---|---|---|
| 0 | complex | 1 | 25 | 39.88 | 0.61 |
| 1 | solvent | 1 | 25 | 37.9 | 1.4 |
| 2 | standard_state_correction_A | 1 | 25 | -9.2 | 0.0 |
| 3 | standard_state_correction_B | 1 | 25 | 9.3 | 0.0 |
| 4 | complex | 1 | 7a | -0.82 | 0.71 |
| 5 | complex | 1 | 7a | -0.35 | 0.69 |
| 6 | complex | 1 | 7a | -2.10 | 0.62 |
| 7 | solvent | 1 | 7a | -1.9 | 1.2 |
| 8 | solvent | 1 | 7a | -1.5 | 1.4 |
| 9 | solvent | 1 | 7a | -1.6 | 1.4 |
| 10 | standard_state_correction_A | 1 | 7a | -8.8 | 0.0 |
| 11 | standard_state_correction_A | 1 | 7a | -9.0 | 0.0 |
| 12 | standard_state_correction_A | 1 | 7a | -9.1 | 0.0 |
| 13 | standard_state_correction_B | 1 | 7a | 9.1 | 0.0 |
| 14 | standard_state_correction_B | 1 | 7a | 9.1 | 0.0 |
| 15 | standard_state_correction_B | 1 | 7a | 9.0 | 0.0 |
| 16 | complex | 1 | 7b | 2.48 | 0.84 |
| 17 | complex | 1 | 7b | 3.63 | 0.74 |
| 18 | solvent | 1 | 7b | 2.8 | 1.3 |
| 19 | solvent | 1 | 7b | 3.1 | 1.3 |
| 20 | standard_state_correction_A | 1 | 7b | -9.0 | 0.0 |
| 21 | standard_state_correction_A | 1 | 7b | -9.3 | 0.0 |
| 22 | standard_state_correction_B | 1 | 7b | 9.1 | 0.0 |
| 23 | standard_state_correction_B | 1 | 7b | 9.3 | 0.0 |
| 24 | complex | 7a | 7b | 3.17 | 0.80 |
| 25 | complex | 7a | 7b | 1.65 | 0.68 |
| 26 | complex | 7a | 7b | -2.06 | 0.58 |
| 27 | solvent | 7a | 7b | -1.2 | 1.3 |
| 28 | solvent | 7a | 7b | -0.9 | 1.3 |
| 29 | solvent | 7a | 7b | -1.1 | 1.2 |
| 30 | standard_state_correction_A | 7a | 7b | -9.0 | 0.0 |
| 31 | standard_state_correction_A | 7a | 7b | -9.4 | 0.0 |
| 32 | standard_state_correction_A | 7a | 7b | -9.0 | 0.0 |
| 33 | standard_state_correction_B | 7a | 7b | 8.9 | 0.0 |
| 34 | standard_state_correction_B | 7a | 7b | 9.3 | 0.0 |
| 35 | standard_state_correction_B | 7a | 7b | 9.0 | 0.0 |