TMol How-To Guide — Metals and Nucleic Acids#
3. Metals and Nucleic Acids#
Metal ions and nucleic acids both score correctly with no special setup – this notebook is mostly worked examples confirming that. It assumes the basics from tmol_how_to_guide.ipynb and scoring/relax from 02_scoring_and_optimization.ipynb.
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
Metals#
Handled by default#
Metal ions (ZN, CA, MG, FE, …) are ordinary residue types in the default chemical database, not ligands – detecting a metal ion and building its coordinating bonds needs no prepare_ligands flag. A MetalCoordination energy term scoring that geometry is included in beta2016_score_function by default, at weight 1.0.
We use a small curated fixture from TMol’s own test suite – a tetrahedral zinc site from 3KS3 (carbonic anhydrase II) – rather than fetching the full deposited structure, since it’s already cropped to the site of interest.
[2]:
metal_pose_stack = pose_stack_from_file(data_path("metal_fixtures", "zn_tetrahedral_3ks3.cif.zst"), device)
print(f"n_poses={metal_pose_stack.n_poses}, max_n_blocks={metal_pose_stack.max_n_blocks}")
viewer = tmol.view(metal_pose_stack)
viewer.show()
HIS is seen only at a chain terminus and N, ND1, NE2 could each carry its other connection; it is read as the conventional backbone through N. A copy of the residue in a chain would settle it.
Unrecognized 3lc GOL
n_poses=1, max_n_blocks=258
/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.
Determined Coordination - TODO#
TODO - I dont really understand what determined coordination is, but this is clearly too complicated. Claude feedback below:
There is no single “report coordination” function; the real pattern is to find the metal’s block type (it declares coordination sites via tmol.database.chemical.site_connections), then walk pose_stack.inter_residue_connections64 to see which residue fills each site.
[3]:
def metal_sites(pose_stack, pose, block):
"""(metal atom, geometry, donor labels) for each metal in a block; None marks an open site."""
bt = pose_stack.block_type(pose, block)
result = []
for site in bt.metal_sites:
donors = []
for conn_name in site.site_connections:
conn = bt.connection_to_cidx[conn_name]
partner, partner_conn = pose_stack.inter_residue_connections64[pose, block, conn].tolist()
if partner < 0:
donors.append(None)
continue
partner_bt = pose_stack.block_type(pose, partner)
atom = partner_bt.atoms[partner_bt.ordered_connection_atoms[partner_conn]].name
donors.append(f"{residue_label(pose_stack, pose, partner)} {atom}")
result.append((site.metal_atom, site.geometry, donors))
return result
cif_text = fetch_cif("3F7L")
pose_stack = pose_stack_from_file(io.StringIO(cif_text), device)
metals = metal_blocks(pose_stack, 0)
# Coordination found at construction
for block in metals:
for metal_atom, geometry, donors in metal_sites(pose_stack, 0, block):
print(f"{residue_label(pose_stack, 0, block)} {metal_atom} ({geometry}):")
for site, donor in enumerate(donors):
print(f" site {site}: {donor or 'open'}")
# Relax, then compare metal_coordination energy per metal
score_function = beta2016_score_function(device)
before = block_energies(score_function, pose_stack, ScoreType.metal_coordination, 0, metals)
torch.manual_seed(SEED)
relaxed = tmol.fast_relax(
pose_stack.clone(),
score_function,
PackerPalette(),
tmol.CartesianMoveMap(),
tmol.FoldForest.reasonable_fold_forest(pose_stack),
num_repeats=1,
)
after = block_energies(score_function, relaxed, ScoreType.metal_coordination, 0, metals)
for block, e0, e1 in zip(metals, before, after):
print(f"{residue_label(pose_stack, 0, block)}: metal_coordination {e0:.3f} -> {e1:.3f}")
/home/runner/work/tmol/tmol/tmol/io/_atomworks_reader.py:893: UserWarning: Keeping one residue per site: dropping C:202 CU
array = one_residue_per_site(array)
Unrecognized 3lc ACY
Unrecognized 3lc SO4
A 201 CU1 CU (trigonal):
site 0: A 44 HIS ND1
site 1: A 46 HIS NE2
site 2: A 118 HIS NE2
A 203 ZN ZN (tetrahedral):
site 0: A 61 HIS ND1
site 1: A 69 HIS ND1
site 2: A 81 ASP OD1
site 3: A 78 HIS ND1
A 204 NA NA (irregular):
site 0: A 90 ASP OD1
site 1: open
site 2: open
site 3: open
site 4: open
site 5: open
site 6: open
site 7: open
A 201 CU1: metal_coordination 1.784 -> 0.400
A 203 ZN: metal_coordination 1.082 -> 0.466
A 204 NA: metal_coordination 21.387 -> 0.018
Relax it (using the ligand-extended parameter database, since the glycerol needs it) and read off the metal_coordination term.
[4]:
metal_score_function = beta2016_score_function(device, param_db=param_db)
metal_fold_forest = tmol.FoldForest.reasonable_fold_forest(metal_pose_stack)
torch.manual_seed(SEED)
metal_relaxed = tmol.fast_relax(
metal_pose_stack.clone(), metal_score_function, PackerPalette(), tmol.CartesianMoveMap(), metal_fold_forest,
num_repeats=1,
)
metal_scorer = metal_score_function.render_whole_pose_scoring_module(metal_relaxed)
weighted_terms = metal_scorer(metal_relaxed.coords, sum_terms=False, apply_weights=True)
metal_coordination_index = metal_score_function.all_score_types().index(ScoreType.metal_coordination)
print(f"metal_coordination energy: {float(weighted_terms[metal_coordination_index, 0].detach()):.2f}")
display(tmol.overlay_view(
{"original": metal_pose_stack, "relaxed": metal_relaxed},
notes={
"original": "pre-relax",
"relaxed": f"metal_coordination energy: {float(weighted_terms[metal_coordination_index, 0].detach()):.2f}",
},
))
metal_coordination energy: 0.28
/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"))
Metal cluster example#
Multi-ion cofactors (iron-sulfur clusters, and similar) are supported the same way – still no prepare_ligands flag, no cluster-specific API. Here a [2Fe-2S] cluster from a ferredoxin (6E6R), each iron coordinated by two cysteines.
[5]:
fes_pose_stack = pose_stack_from_file(data_path("metal_fixtures", "fes_ferredoxin_6e6r.cif.zst"), device)
fes_score_function = beta2016_score_function(device, param_db=param_db)
fes_scorer = fes_score_function.render_whole_pose_scoring_module(fes_pose_stack)
fes_weighted_terms = fes_scorer(fes_pose_stack.coords, sum_terms=False, apply_weights=True)
print(f"n_poses={fes_pose_stack.n_poses}, max_n_blocks={fes_pose_stack.max_n_blocks}")
print(f"metal_coordination energy: {float(fes_weighted_terms[metal_coordination_index, 0]):.2f}")
viewer = tmol.view(fes_pose_stack)
viewer.show()
n_poses=1, max_n_blocks=57
metal_coordination energy: 37.32
/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.
Nucleic acids#
Handled by default#
DNA and RNA residues are likewise ordinary chemistry in the default database. An NaTorsion term scoring their backbone, sugar-pucker, and glycosidic torsions is included in beta2016_score_function by default; it’s an exact zero on protein-only poses, so mixed protein/NA structures score correctly with no special handling either.
Example#
1HDD, an engrailed homeodomain bound to its DNA operator, gives us both an NA torsion energy and a protein-DNA interface in one structure.
[6]:
dna_cif_text = fetch_cif("1HDD")
dna_pose_stack = pose_stack_from_file(io.StringIO(dna_cif_text), device)
print(f"n_poses={dna_pose_stack.n_poses}, max_n_blocks={dna_pose_stack.max_n_blocks}")
viewer = tmol.view(dna_pose_stack)
viewer.show()
na_score_function = beta2016_score_function(device, param_db=param_db)
na_fold_forest = tmol.FoldForest.reasonable_fold_forest(dna_pose_stack)
torch.manual_seed(SEED)
dna_relaxed = tmol.fast_relax(
dna_pose_stack.clone(), na_score_function, PackerPalette(), tmol.CartesianMoveMap(), na_fold_forest,
num_repeats=1,
)
dna_scorer = na_score_function.render_whole_pose_scoring_module(dna_relaxed)
dna_weighted_terms = dna_scorer(dna_relaxed.coords, sum_terms=False, apply_weights=True)
na_torsion_index = na_score_function.all_score_types().index(ScoreType.na_torsion)
print(f"na_torsion energy: {float(dna_weighted_terms[na_torsion_index, 0].detach()):.2f}")
display(tmol.overlay_view(
{"original": dna_pose_stack, "relaxed": dna_relaxed},
notes={
"original": "pre-relax",
"relaxed": f"na_torsion energy: {float(dna_weighted_terms[na_torsion_index, 0].detach()):.2f}",
},
))
ARG is seen only at a chain terminus and N, NE, NH1, NH2 could each carry its other connection; it is read as the conventional backbone through N. A copy of the residue in a chain would settle it.
n_poses=1, max_n_blocks=156
/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.
na_torsion energy: 34.96
/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"))
The protein-DNA interface energy reuses the block-pair scoring technique from 02_scoring_and_optimization.ipynb, with a polymer_type-based mask standing in for a chain mask – selections like this are covered in general in 05_selections.ipynb.
[7]:
pbt = dna_relaxed.packed_block_types
protein_mask = torch.zeros((1, dna_relaxed.max_n_blocks), dtype=torch.bool)
dna_mask = torch.zeros((1, dna_relaxed.max_n_blocks), dtype=torch.bool)
for block in range(dna_relaxed.max_n_blocks):
bt_ind = int(dna_relaxed.block_type_ind64[0, block])
if bt_ind < 0:
continue
polymer_type = pbt.active_block_types[bt_ind].properties.polymer.polymer_type
protein_mask[0, block] = polymer_type == "amino_acid"
dna_mask[0, block] = polymer_type == "nucleic_acid"
dna_block_scorer = na_score_function.render_block_pair_scoring_module(dna_relaxed)
dna_block_pair = dna_block_scorer(dna_relaxed.coords)[0].detach()
interface_score = (
dna_block_pair[protein_mask[0]][:, dna_mask[0]].sum()
+ dna_block_pair[dna_mask[0]][:, protein_mask[0]].sum()
)
print(f"protein blocks: {int(protein_mask.sum())}, DNA blocks: {int(dna_mask.sum())}")
print(f"protein-DNA interface score: {float(interface_score):.2f}")
protein blocks: 114, DNA blocks: 42
protein-DNA interface score: -205.33
The same recipe works unchanged for RNA. 1EHT is a theophylline-binding RNA aptamer; here the interesting “other side” is its bound ligand rather than a second polymer chain, covered in 04_ligands_and_ncaas.ipynb, so we report just the NA torsion energy.
[8]:
rna_cif_text = fetch_cif("1EHT")
rna_pose_stack = pose_stack_from_file(io.StringIO(rna_cif_text), device)
print(f"n_poses={rna_pose_stack.n_poses}, max_n_blocks={rna_pose_stack.max_n_blocks}")
viewer = tmol.view(rna_pose_stack)
viewer.show()
torch.manual_seed(SEED)
rna_relaxed = tmol.fast_relax(
rna_pose_stack.clone(), na_score_function, PackerPalette(), tmol.CartesianMoveMap(),
tmol.FoldForest.reasonable_fold_forest(rna_pose_stack), num_repeats=1,
)
rna_scorer = na_score_function.render_whole_pose_scoring_module(rna_relaxed)
rna_weighted_terms = rna_scorer(rna_relaxed.coords, sum_terms=False, apply_weights=True)
print(f"na_torsion energy: {float(rna_weighted_terms[na_torsion_index, 0].detach()):.2f}")
display(tmol.overlay_view(
{"original": rna_pose_stack, "relaxed": rna_relaxed},
notes={
"original": "pre-relax",
"relaxed": f"na_torsion energy: {float(rna_weighted_terms[na_torsion_index, 0].detach()):.2f}",
},
))
G is seen only at a chain terminus and N1, N2 could each carry its other connection; none of them is a conventional backbone, so it cannot be told apart. A copy of the residue in a chain would settle it.
Unrecognized 3lc TEP
n_poses=1, max_n_blocks=33
/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.
na_torsion energy: 62.66
/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"))