🆕 Print the Summary of Protein-Ligand Interactions for a Set of Molecules
Problem
You want to perceive protein-ligand interactions and print them to the console.
Ingredients
|
Difficulty Level
🌶️ 🌶️
Download
Source Code
interactions2summary
#!/usr/bin/env python3
# (C) 2026 Cadence Design Systems, Inc. (Cadence)
# All rights reserved.
# TERMS FOR USE OF SAMPLE CODE The software below ("Sample Code") is
# provided to current licensees or subscribers of Cadence products or
# SaaS offerings (each a "Customer").
# Customer is hereby permitted to use, copy, and modify the Sample Code,
# subject to these terms. Cadence claims no rights to Customer's
# modifications. Modification of Sample Code is at Customer's sole and
# exclusive risk. Sample Code may require Customer to have a then
# current license or subscription to the applicable Cadence offering.
# THE SAMPLE CODE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED. OPENEYE DISCLAIMS ALL WARRANTIES, INCLUDING, BUT
# NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
# PARTICULAR PURPOSE AND NONINFRINGEMENT. In no event shall Cadence be
# liable for any damages or liability in connection with the Sample Code
# or its use.
"""Print the summary of the interaction of a protein and a series of molecules to console."""
import argparse
import os
import pathlib
import sys
from collections import Counter
from collections.abc import Iterator
import rich
import rich.console
import rich.table
import rich.text
from openeye import oechem
from rich.progress import BarColumn, Progress, TextColumn, TimeElapsedColumn
from rich_argparse import HelpPreviewAction, RichHelpFormatter
__SCRIPT_NAME__ = pathlib.Path(__file__).absolute().stem
__SCRIPT_DESC__ = (
"Print the summary of the interaction of a protein and a series of molecules."
)
__SCRIPT_TOOLKITS__ = ["oechem"]
__SCRIPT_KEYWORDS__ = [
"perception",
"active-site",
"protein-ligand",
"interactions",
"display",
]
__SCRIPT_CATEGORIES__ = ["protein-ligand interactions"]
def parse_options() -> argparse.Namespace:
"""Set up command line options."""
parser = argparse.ArgumentParser(
add_help=True,
formatter_class=RichHelpFormatter,
description="[yellow]" + __SCRIPT_DESC__ + "[/yellow]",
)
# input options
input_group = parser.add_argument_group("Input")
input_group.add_argument(
"--protein",
type=str,
required=False,
metavar="PDB/DU-FILE",
help="design unit or apo protein structure (.pdb, .oedu)",
)
input_group.add_argument(
"--mol",
type=str,
required=False,
metavar="MOL-FILE",
help="input of molecules at active site (.sdf, .oeb)",
)
input_group.add_argument(
"--serialized",
type=str,
required=False,
metavar="OEB/JSON-FILE",
help="input file with serialized protein-ligand interactions (.oeb or .json)",
)
parser.add_argument("--help-image", action=HelpPreviewAction)
parser.add_argument(
"--save-console-svg",
default=False,
action="store_true",
help=f"run command and capture console output in {__SCRIPT_NAME__}.svg file",
)
args = parser.parse_args()
# validate that either (--protein and --mol) or --serialized is provided
if not any([args.serialized, (args.protein and args.mol)]):
parser.error("Either (--protein and --mol) or --serialized must be provided")
if args.serialized and (args.protein or args.mol):
parser.error("Cannot use --serialized together with --protein or --mol")
if (args.protein and not args.mol) or (args.mol and not args.protein):
parser.error("Both --protein and --mol must be provided together")
return args
def main() -> int:
"""Print interactions to console."""
args = parse_options()
console = rich.console.Console(record=args.save_console_svg)
mols: list[oechem.OEMolBase] = []
if args.protein and args.mol:
protein, mols = get_molecules(args.protein, args.mol)
else:
protein, mols = get_molecules_from_serialized(args.serialized)
interactions: list[tuple[str, oechem.OEInteractionHintTypeBase]] = list(
retrieve_interactions(
protein, mols, console, serialized=args.serialized is not None
)
)
unique_residues: list[str] = list({r for r, _ in interactions})
unique_residues.sort()
only_contact_residues: list[str] = []
accumulated_interactions: list[tuple[str, Counter]] = []
for residue_str in unique_residues:
residue_interactions = [
inter_type for r, inter_type in interactions if r == residue_str
]
non_contact_residue_interactions = [
inter_type
for inter_type in residue_interactions
if (inter_type != oechem.OEContactInteractionHint())
]
if len(non_contact_residue_interactions) == 0:
only_contact_residues.append(residue_str)
else:
interactions_count = Counter(
i.GetName().removeprefix("bio:active-site:")
for i in non_contact_residue_interactions
)
accumulated_interactions.append((residue_str, interactions_count))
console.print(
f"[green]Residue(s) with only contact interactions:[/green] {'; '.join(only_contact_residues)} \n",
highlight=False,
)
max_bar_width = 40
max_interaction_label_length = max(
len(inter_name)
for _, interactions_count in accumulated_interactions
for inter_name in interactions_count
)
max_interaction_count = max(
count
for _, interactions_count in accumulated_interactions
for count in interactions_count.values()
)
main_table = rich.table.Table(
title=rich.text.Text(
f"Interactions Summary for {len(mols)} molecules", style="bold"
),
box=rich.box.DOUBLE_EDGE,
)
main_table.add_column("Residues")
main_table.add_column("Interactions", justify="center")
for residue_str, interactions_count in accumulated_interactions:
sub_table = rich.table.Table(show_header=False, show_edge=False)
for inter_name, count in interactions_count.items():
inter_name_text = rich.text.Text(
f"{inter_name:>{max_interaction_label_length}} | ", style="bold"
)
bar = rich.text.Text(
"█" * max(1, int(count * max_bar_width / max_interaction_count)),
style=f"bold {color_by_interaction_name(inter_name)}",
)
sub_table.add_row(
inter_name_text + bar + rich.text.Text(f" {count}", style="bold white")
)
main_table.add_row(residue_str, sub_table)
main_table.add_section()
console.print(main_table)
if args.save_console_svg:
console.save_svg(f"{__SCRIPT_NAME__}.svg", title="output")
return os.EX_OK
def retrieve_interactions(
receptor: oechem.OEMolBase,
mols: list[oechem.OEMolBase],
console: rich.console.Console,
serialized: bool,
) -> Iterator[tuple[str, oechem.OEInteractionHintTypeBase]]:
"""Accumulate interactions for all molecules."""
task_name = (
"Retrieving interactions from serialized data"
if serialized
else "Perceiving interactions"
)
with Progress(
TextColumn("{task.description}"),
BarColumn(bar_width=60),
TextColumn("{task.percentage:3.1f}%"),
TimeElapsedColumn(),
transient=False,
disable=len(mols) < 10, # noqa: PLR2004
console=console,
) as progress:
perception_task = progress.add_task(f"[blue]{task_name}", total=len(mols))
for idx, mol in enumerate(mols, start=1):
if serialized:
active_site = oechem.OEInteractionHintContainer()
if not oechem.OEConstructInteractionHintContainer(
active_site, receptor, mol
):
console.print(
f"[red]Error: Cannot construct active site for molecule at index {idx}![/red]"
)
continue
else:
active_site = oechem.OEInteractionHintContainer(receptor, mol)
if not oechem.OEIsValidActiveSite(active_site):
console.print(
f"[red]Error: Cannot initialize active site for molecule at index {idx}![/red]"
)
continue
if not oechem.OEPerceiveInteractionHints(active_site):
console.print(
f"[red]Error: Cannot perceive interactions for molecule at index {idx}![/red]"
)
continue
for inter in active_site.GetInteractions(oechem.OEIsInterInteractionHint()):
residue_str = _get_interaction_residue_str(inter)
if residue_str != "":
yield (residue_str, inter.GetInteractionType().CreateCopy())
progress.update(perception_task, advance=1)
def _get_interaction_residue_str(inter: oechem.OEInteractionHint) -> str:
prot_frag = inter.GetFragment(oechem.OEProteinInteractionHintComponent())
if prot_frag is None:
return ""
# in very rare cases atoms interactions can belong to multiple residues
residues: set[oechem.OEResidue] = set()
for prot_atom in prot_frag.GetAtoms():
residues.add(oechem.OEAtomGetResidue(prot_atom))
return ";".join(_get_residue_str(r) for r in residues)
def _get_residue_str(residue: oechem.OEResidue) -> str:
return (
f"{residue.GetName():3s} {residue.GetResidueNumber():4d} {residue.GetChainID()}"
)
def get_molecules(
protein_filename: str, mol_filename: str
) -> tuple[oechem.OEMolBase, list[oechem.OEMolBase]]:
"""Load the protein molecule from a file."""
protein = oechem.OEGraphMol()
# read protein
du = oechem.OEDesignUnit()
if oechem.OEIsReadableDesignUnit(protein_filename) and oechem.OEReadDesignUnit(
protein_filename, du
):
if not du.GetComponents(
protein, oechem.OEDesignUnitComponents_TargetComplexNoSolvent
):
oechem.OEThrow.Fatal("Could not extract protein from the design unit.")
else:
# fall back to reading protein as a molecule
ifs = oechem.oemolistream()
if not ifs.open(protein_filename):
oechem.OEThrow.Fatal(f"Unable to open {protein_filename} for reading")
if not oechem.OEReadMolecule(ifs, protein):
oechem.OEThrow.Fatal(f"Unable to read protein from {protein_filename}")
if oechem.OEGetDimensionFromCoords(protein) != 3: # noqa: PLR2004
oechem.OEThrow.Fatal(f"Protein molecule in {protein_filename} is not 3D")
protein.SetDimension(3)
# read molecules at active site
ifs = oechem.oemolistream()
if not ifs.open(mol_filename):
oechem.OEThrow.Fatal(f"Unable to open {mol_filename} for reading!")
mols: list[oechem.OEMolBase] = [oechem.OEGraphMol(m) for m in ifs.GetOEGraphMols()]
if len(mols) == 0:
oechem.OEThrow.Fatal(f"No molecules read from {mol_filename}!")
return protein, mols
def get_molecules_from_serialized(
filename: str,
) -> tuple[oechem.OEMolBase, list[oechem.OEMolBase]]:
"""Load the protein molecule and molecules at active site from a file with serialized interactions."""
ifs = oechem.oemolistream()
if not ifs.open(filename):
oechem.OEThrow.Fatal(f"Unable to open {filename} for reading!")
if ifs.GetFormat() not in [oechem.OEFormat_OEB, oechem.OEFormat_JSON]:
oechem.OEThrow.Fatal(
"Unsupported file format for serialization! Only .oeb and .json are supported."
)
mols: list[oechem.OEMolBase] = [oechem.OEGraphMol(m) for m in ifs.GetOEGraphMols()]
if not all(oechem.OEHasInteractionsHintSerializationData(m) for m in mols):
oechem.OEThrow.Fatal(
f"All molecules in {filename} should have serialized interactions!"
)
if len(mols) < 2: # noqa: PLR2004
oechem.OEThrow.Fatal(
f"File should have protein and at least one other molecule {filename}!"
)
# protein is expected to be first
protein: oechem.OEMolBase = mols.pop(0)
if not oechem.OEHasInteractionHintSerializationIds(protein):
oechem.OEThrow.Fatal(
f"The first molecule in {filename} should be the protein with serialization ids!"
)
if protein.NumAtoms() < max(m.NumAtoms() for m in mols):
oechem.OEThrow.Fatal(f"Protein should be the first in {filename}!")
return protein, mols
def color_by_interaction_name(inter_name: str) -> str: # noqa: C901, PLR0911, PLR0912
"""Return a color string for a given interaction type."""
if inter_name.startswith("bio:active-site"):
inter_name = inter_name.removeprefix("bio:active-site:")
match inter_name:
case "covalent":
return "#bfbfbf"
case "clash":
return "#872924"
case _ if "clash" in inter_name:
return "red"
case _ if inter_name.startswith("halogen"):
return "#f5af91"
case _ if inter_name.startswith("stacking"):
return "#5faf5f"
case _ if inter_name.startswith("cationpi"):
return "#aedae1"
case _ if inter_name.startswith("chelator"):
return "#ffdb7f"
case _ if inter_name.startswith("salt-bridge:ligand-protein+"):
return "#5f5fff"
case _ if inter_name.startswith("salt-bridge:ligand+protein-"):
return "#ff8c8c"
case _ if inter_name.startswith("hbond:ligand2protein"):
return "#ffc8c8"
case _ if inter_name.startswith("hbond:non-ideal-ligand2protein"):
return "#ffe6e6"
case _ if inter_name.startswith("hbond:protein2ligand"):
return "#c8c8ff"
case _ if inter_name.startswith("hbond:non-ideal-protein2ligand"):
return "#e6e6ff"
case _:
return "white"
setattr(main, "__SCRIPT_NAME__", __SCRIPT_NAME__)
setattr(main, "__SCRIPT_DESC__", __SCRIPT_DESC__)
setattr(main, "__SCRIPT_TOOLKITS__", __SCRIPT_TOOLKITS__)
setattr(main, "__SCRIPT_KEYWORDS__", __SCRIPT_KEYWORDS__)
setattr(main, "__SCRIPT_CATEGORIES__", __SCRIPT_CATEGORIES__)
if __name__ == "__main__":
sys.exit(main())
Discussion
Usage
See Download section to download the script.
> interactions2summary --help
- --protein
- --mol
Printing a summary of interactions for the protein of CDK5.oedu
and molecules of CDK5-hits.sdf:
> interactions2summary --protein CDK5.oedu --mol CDK5-hits.sdf
- --serialized
The same summary can be printed from the serialized interactions file
CDK5-protein-hits.oeb file
generated with 🆕 Serialize Protein-Ligand Interactions script.
> interactions2summary --serialized CDK5-protein-hits.oeb
See also in OEChem TK manual
API
OEDesignUnit class
OEHasInteractionHintType predicate
OEInteractionHint class
OEPerceiveInteractionHints function
OEResidue class