TMol How-To Guide — Ligands and NCAAs#

Open In Colab

4. Ligands and NCAAs#

Small-molecule ligands and non-canonical amino acids both go through the same preparation pipeline (prepare_ligands=True) seen briefly in section 1. This notebook goes further: reusing prepared chemistry, saving/loading it, splitting a ligand into independently-scored fragments, and the handful of NCAAs that don’t need the pipeline at all.

Setup#

Wheel install, imports, seeding, and device selection live in notebook_setup.py.

[1]:

try:
    import google.colab  # noqa: F401
except ImportError:
    IN_COLAB = False
else:
    IN_COLAB = True

if IN_COLAB:
    from urllib.request import urlretrieve

    urlretrieve(
        "https://raw.githubusercontent.com/RosettaCommons/tmol/"
        "main/docs/tutorial/notebook_setup.py",
        "notebook_setup.py",
    )

from notebook_setup import *  # noqa: F401,F403

device, param_db = setup()

Environment variable CCD_MIRROR_PATH not set. Will not be able to use function requiring this variable. To set it you may:
  (1) add the line 'export VAR_NAME=path/to/variable' to your .bashrc or .zshrc file
  (2) set it in your current shell with 'export VAR_NAME=path/to/variable'
  (3) write it to a .env file in the root of the atomworks.io repository
Environment variable PDB_MIRROR_PATH not set. Will not be able to use function requiring this variable. To set it you may:
  (1) add the line 'export VAR_NAME=path/to/variable' to your .bashrc or .zshrc file
  (2) set it in your current shell with 'export VAR_NAME=path/to/variable'
  (3) write it to a .env file in the root of the atomworks.io repository

Ligands#

Automatic pipeline#

We return to 3PTB (trypsin + the small-molecule inhibitor benzamidine, BEN) from section 1 and go further: relax the whole complex, then read off a ligand-protein interaction energy and the ligand’s own internal strain.

[2]:

from tmol.io import atom_array_from_file, pose_stack_from_atom_array, build_context_from_atom_array

ligand_cif_text = fetch_cif("3PTB")
ligand_structure = atom_array_from_file(io.StringIO(ligand_cif_text))
pose_stack, build_context = pose_stack_from_atom_array(
    ligand_structure, device, prepare_ligands=True, param_db=param_db, return_context=True
)
score_function = beta2016_score_function(device, param_db=build_context.parameter_database)
fold_forest = tmol.FoldForest.reasonable_fold_forest(pose_stack)

torch.manual_seed(SEED)
relaxed = tmol.fast_relax(
    pose_stack.clone(), score_function, PackerPalette(), tmol.CartesianMoveMap(), fold_forest, num_repeats=1
)
print(f"n_poses={relaxed.n_poses}, max_n_blocks={relaxed.max_n_blocks}")

ligand_atom_mask = res_mask_to_coord_mask(pose_stack, tmol.get_ligand_mask(pose_stack))
viewer = tmol.view(pose_stack, highlighted=ligand_atom_mask[pose_stack.real_atoms])
viewer.show()

display(tmol.overlay_view(
    {"original": pose_stack, "relaxed": relaxed},
    notes={
        "original": "pre-relax (3PTB, trypsin + BEN)",
        "relaxed": f"post fast_relax: n_poses={relaxed.n_poses}, max_n_blocks={relaxed.max_n_blocks}",
    },
))

==============================
*** Open Babel Warning  in WriteMolecule
  No 2D or 3D coordinates exist. Stereochemical information will be stored using an Open Babel extension. To generate 2D or 3D coordinates instead use --gen2D or --gen3D.

n_poses=1, max_n_blocks=225

/home/runner/work/tmol/tmol/tmol/io/_pose_stack_from_biotite.py:1401: UserWarning: Renumbering chain A of the output AtomArray: biotite's CIF writer does not support numbering that decreases within a chain
  warnings.warn(
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

/home/runner/work/tmol/tmol/tmol/io/_pose_stack_from_biotite.py:1401: UserWarning: Renumbering chain A of the output AtomArray: biotite's CIF writer does not support numbering that decreases within a chain
  warnings.warn(
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))
/home/runner/work/tmol/tmol/tmol/io/_pose_stack_from_biotite.py:1401: UserWarning: Renumbering chain A of the output AtomArray: biotite's CIF writer does not support numbering that decreases within a chain
  warnings.warn(
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))
/home/runner/work/tmol/tmol/tmol/io/_pose_stack_from_biotite.py:1401: UserWarning: Renumbering chain A of the output AtomArray: biotite's CIF writer does not support numbering that decreases within a chain
  warnings.warn(
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))
/home/runner/work/tmol/tmol/tmol/io/_pose_stack_from_biotite.py:1401: UserWarning: Renumbering chain A of the output AtomArray: biotite's CIF writer does not support numbering that decreases within a chain
  warnings.warn(
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

tmol.ops.calculate_block_pair_ddg sums the cross-mask block-pair interaction between a residue mask and everything else – here, between the ligand block and the rest of the complex. (Despite its name it reports that raw interaction, not a true bound/unbound ddG.)

[3]:

ligand_mask = torch.zeros((1, relaxed.max_n_blocks), dtype=torch.bool)
for block in range(relaxed.max_n_blocks):
    bt_ind = int(relaxed.block_type_ind64[0, block])
    if bt_ind >= 0 and relaxed.packed_block_types.active_block_types[bt_ind].name3 == "BEN":
        ligand_mask[0, block] = True

interaction = calculate_block_pair_ddg(relaxed, ligand_mask, sfxn=score_function, minimize=False, pack=False)
print(f"ligand-protein interaction energy: {float(interaction[0]):.2f}")

ligand-protein interaction energy: -27.49

gen_torsions is the generic-atom-type torsional term that shows up in a per-term score breakdown (see section 2) – it fills in wherever a residue lacks dedicated per-residue cartbonded parameters (canonical amino acids score an exact zero here), which in practice means ligands and other novel chemistry. It’s a reasonable proxy for “ligand strain.”

[4]:

relaxed_scorer = score_function.render_whole_pose_scoring_module(relaxed)
weighted_terms = relaxed_scorer(relaxed.coords, sum_terms=False, apply_weights=True)
gen_torsions_index = score_function.all_score_types().index(ScoreType.gen_torsions)
print(f"gen_torsions (ligand strain): {float(weighted_terms[gen_torsions_index, 0].detach()):.2f}")

gen_torsions (ligand strain): 0.61

Reusing a build context#

Preparing a ligand’s chemistry (typing, charges, conformer) is the expensive part of loading. 1BTY is a different trypsin-benzamidine structure – passing the earlier build_context into context= skips re-preparing BEN entirely. build_context_from_atom_array extends a copy of the parameter database rather than mutating it, so the original database is untouched.

[5]:

other_cif_text = fetch_cif("1BTY")
other_structure = atom_array_from_file(io.StringIO(other_cif_text))
other_pose_stack = pose_stack_from_atom_array(other_structure, device, context=build_context, no_optH=True)
print(f"reused context to build a second pose: {other_pose_stack.max_n_blocks} blocks, no re-preparation needed")

other_ligand_mask = res_mask_to_coord_mask(other_pose_stack, tmol.get_ligand_mask(other_pose_stack))
viewer = tmol.view(other_pose_stack, highlighted=other_ligand_mask[other_pose_stack.real_atoms])
viewer.show()

reused context to build a second pose: 225 blocks, no re-preparation needed

/home/runner/work/tmol/tmol/tmol/io/_pose_stack_from_biotite.py:1401: UserWarning: Renumbering chain A of the output AtomArray: biotite's CIF writer does not support numbering that decreases within a chain
  warnings.warn(
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Saving and loading .tmol ligand parameters#

prepare_ligands (the lower-level function prepare_ligands=True calls into) can also write everything it derives for a ligand to a .tmol YAML file via params_output, given just that ligand’s own atoms – not the whole complex, since the pipeline expects a self-contained molecule.

[6]:

import tempfile
from pathlib import Path
from tmol.ligand._preparation import prepare_ligands

ligand_only_atoms = ligand_structure[ligand_structure.res_name == "BEN"]
tmol_path = Path(tempfile.gettempdir()) / "ben.tmol"
_, _ = prepare_ligands(ligand_only_atoms, param_db=param_db, params_output=str(tmol_path), seed=SEED)
print(f"wrote {tmol_path} ({tmol_path.stat().st_size} bytes)")

wrote /tmp/ben.tmol (9485 bytes)

inject_params_file reads a .tmol file directly into a ParameterDatabase – here the global default – without touching RDKit or OpenBabel again. Loading a full complex can then skip prepare_ligands=True entirely and just point ligand_params_files at the file.

[7]:

from tmol.ligand import inject_params_file

reloaded_db = inject_params_file(ParameterDatabase.get_default(), str(tmol_path))
reloaded_pose = pose_stack_from_atom_array(
    ligand_structure, device, param_db=ParameterDatabase.get_default(),
    ligand_params_files=[str(tmol_path)], no_optH=True,
)
print(f"rebuilt from the saved file: {reloaded_pose.max_n_blocks} blocks")

rebuilt from the saved file: 225 blocks

Ligand “splitting”#

A large or multi-fragment ligand can be cut into separately-tracked blocks: chemistry is perceived once for the whole molecule, and a per-atom tmol_fragment_id annotation (FRAGMENT_ID_ANNOTATION) marks which side of a chosen bond each atom falls on. There’s no automatic cut-point selection – the cut below comes from picking a bridging bond in a test fixture ligand and splitting its bond graph into the two halves that remain connected without it.

[8]:

import numpy as np
import biotite.structure as struc
import biotite.structure.io
from tmol.ligand import FRAGMENT_ID_ANNOTATION, load_params_file
split_structure = biotite.structure.io.load_structure(
    str(data_path("protein_ligand_test", "ace.tmol.nomin.cif")), model=1, include_bonds=True
)
if isinstance(split_structure, struc.AtomArrayStack):
    split_structure = split_structure[0]
split_params_path = data_path("protein_ligand_test", "ace.xtal-lig.mmff94.tmol")
ligand_restype = load_params_file(split_params_path)[0].residue_type

cut = ("C1", "C2")
adjacency = {atom.name: set() for atom in ligand_restype.atoms}
for a, b, *_ in ligand_restype.bonds:
    if {a, b} == set(cut):
        continue
    adjacency[a].add(b)
    adjacency[b].add(a)


def connected_component(start):
    seen, stack = {start}, [start]
    while stack:
        atom = stack.pop()
        for neighbor in adjacency[atom]:
            if neighbor not in seen:
                seen.add(neighbor)
                stack.append(neighbor)
    return seen


fragment_1_atoms = connected_component(cut[0])
fragment_2_atoms = connected_component(cut[1])
print(f"cutting bond {cut} splits the ligand into fragments of {len(fragment_1_atoms)} and {len(fragment_2_atoms)} atoms")

cutting bond ('C1', 'C2') splits the ligand into fragments of 17 and 42 atoms

[9]:

ligand_atom_mask = split_structure.res_name == "LG1"
fragment_ids = np.zeros(split_structure.array_length(), dtype=np.int32)
for atom_index in np.flatnonzero(ligand_atom_mask):
    name = str(split_structure.atom_name[atom_index]).strip()
    if name in fragment_1_atoms:
        fragment_ids[atom_index] = 1
    elif name in fragment_2_atoms:
        fragment_ids[atom_index] = 2
split_structure.set_annotation(FRAGMENT_ID_ANNOTATION, fragment_ids)

split_pose_stack, split_build_context = pose_stack_from_atom_array(
    split_structure, device, prepare_ligands=True, ligand_params_files=[str(split_params_path)],
    no_optH=True, param_db=ParameterDatabase.get_default(), return_context=True,
)
fragment_blocks = [entry.block_ind for entry in split_pose_stack.split_block_mapping.entries]
print(f"ligand split into blocks {fragment_blocks} of {split_pose_stack.max_n_blocks} total")

split_ligand_block_mask = torch.zeros((1, split_pose_stack.max_n_blocks), dtype=torch.bool)
split_ligand_block_mask[0, fragment_blocks] = True
split_ligand_atom_mask = res_mask_to_coord_mask(split_pose_stack, split_ligand_block_mask)
viewer = tmol.view(split_pose_stack, highlighted=split_ligand_atom_mask[split_pose_stack.real_atoms])
viewer.show()

ligand split into blocks [574, 575] of 576 total

/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Each fragment block now scores independently via the same block-pair technique seen earlier (section 2).

[10]:

split_score_function = beta2016_score_function(device, param_db=split_build_context.parameter_database)
split_block_scorer = split_score_function.render_block_pair_scoring_module(split_pose_stack)
split_block_pair = split_block_scorer(split_pose_stack.coords)[0]
for block in fragment_blocks:
    contribution = split_block_pair[block].sum() + split_block_pair[:, block].sum() - split_block_pair[block, block]
    print(f"fragment block {block}: {float(contribution):.2f}")

fragment block 574: -0.73
fragment block 575: 12.12

NCAAs#

D-amino acids and AIB, handled by default#

D-amino acids and AIB (alpha-aminoisobutyric acid) are first-class entries in the default chemical database, not ligands – they load like any other residue, no prepare_ligands flag needed. We use a small curated fixture, a helical peptide built from AIB and leucine.

[11]:

aib_structure = atom_array_from_file(data_path("ncaa_fixtures", "uic1_aib_9bzt.cif"))
aib_pose_stack = pose_stack_from_atom_array(aib_structure, device, param_db=param_db, no_optH=True)
print(f"n_poses={aib_pose_stack.n_poses}, max_n_blocks={aib_pose_stack.max_n_blocks}")
print("sequence (AIB has no one-letter code, shown as X):", tmol.get_sequence(aib_pose_stack))

aib_block_mask = torch.zeros((1, aib_pose_stack.max_n_blocks), dtype=torch.bool)
for block in range(aib_pose_stack.max_n_blocks):
    if aib_pose_stack.block_type(0, block).name3 == "AIB":
        aib_block_mask[0, block] = True
aib_atom_mask = res_mask_to_coord_mask(aib_pose_stack, aib_block_mask)
viewer = tmol.view(aib_pose_stack, highlighted=aib_atom_mask[aib_pose_stack.real_atoms])
viewer.show()

n_poses=1, max_n_blocks=9
sequence (AIB has no one-letter code, shown as X): LXAXLXQXL

/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Most NCAAs go through the ligand pipeline, like SEP#

Less common modifications – SEP (phosphoserine) here – aren’t in the default database, so they load like a free-molecule ligand: as a residue, with its polymer connections replaced by chemical stubs, run through the same prepare_ligands=True pipeline described above, then the stubs restored.

[12]:

sep_structure = atom_array_from_file(data_path("ncaa_fixtures", "phosphopeptide_5ema.cif"))
sep_pose_stack, sep_build_context = pose_stack_from_atom_array(
    sep_structure, device, prepare_ligands=True, param_db=param_db, return_context=True
)
sep_block = next(
    block for block in range(sep_pose_stack.max_n_blocks)
    if sep_pose_stack.block_type(0, block).name3 == "SEP"
)
print(f"n_poses={sep_pose_stack.n_poses}, max_n_blocks={sep_pose_stack.max_n_blocks}, SEP at block {sep_block}")

sep_block_mask = torch.zeros((1, sep_pose_stack.max_n_blocks), dtype=torch.bool)
sep_block_mask[0, sep_block] = True
sep_atom_mask = res_mask_to_coord_mask(sep_pose_stack, sep_block_mask)
viewer = tmol.view(sep_pose_stack, highlighted=sep_atom_mask[sep_pose_stack.real_atoms])
viewer.show()

n_poses=1, max_n_blocks=7, SEP at block 3

/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Packing at a SEP position#

The same samplers as any repack (the Dunbrack library, FixedAAChiSampler, IncludeCurrentSampler) apply here too. build_rotamers exposes the actual generated conformers before packing picks among them.

[13]:

from tmol.pack import SetPackerTask
from tmol.pack.rotamer import build_rotamers

sep_task = PackerTask(sep_pose_stack, PackerPalette())
sep_task.restrict_to_repacking()
sep_shell = torch.zeros_like(sep_task.is_real_block)
sep_shell[0, sep_block] = True
sep_task.disable_packing_by_block_mask(~sep_shell)
sep_task.add_conformer_sampler(create_dunbrack_sampler_from_database(sep_build_context.parameter_database, device))
sep_task.add_conformer_sampler(FixedAAChiSampler())
sep_task.add_conformer_sampler(IncludeCurrentSampler())

_, sep_rotamer_set = build_rotamers(
    sep_pose_stack, SetPackerTask.from_packer_task(sep_task), sep_pose_stack.packed_block_types.chem_db
)
print(f"rotamers generated for the SEP block: {int(sep_rotamer_set.n_rots_for_block[0, sep_block])}")

sep_score_function = beta2016_score_function(device, param_db=sep_build_context.parameter_database)
torch.manual_seed(SEED)
sep_packed = pack_rotamers(sep_pose_stack, sep_score_function, sep_task)
print("packed SEP position successfully")

display(tmol.overlay_view(
    {"original": sep_pose_stack, "packed": sep_packed},
    notes={
        "original": f"SEP at block {sep_block}, pre-packing",
        "packed": f"packed from {int(sep_rotamer_set.n_rots_for_block[0, sep_block])} sampled rotamers",
    },
))

rotamers generated for the SEP block: 244
packed SEP position successfully

/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))
/home/runner/work/tmol/tmol/.venv/lib/python3.12/site-packages/biotite/structure/io/pdb/file.py:629: DeprecationWarning: The chararray class is deprecated and will be removed in a future release. Use an ndarray with a string or bytes dtype instead.
  record = np.char.array(np.where(array.hetero, "HETATM", "ATOM"))

Download this notebook