Skip to content

Commit 76ffd16

Browse files
j-wagsj-wags
authored
Raise an error when loading a PDB file would lead to mangled atom order (#2107)
1 parent 5b4941c commit 76ffd16

File tree

6 files changed

+207
-0
lines changed

6 files changed

+207
-0
lines changed

‎docs/releasehistory.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@ Releases follow the `major.minor.micro` scheme recommended by [PEP440](https://w
1313
### Behavior changes
1414

1515
### Bugfixes
16+
- [PR #2107](https://github.com/openforcefield/openff-toolkit/pull/2107): Fixes [Issue #2093](https://github.com/openforcefield/openff-toolkit/issues/2093), where in rare cases `Topology.from_pdb` would mangle the order/coordinates of loaded atoms, by making it raise a `PDBMoleculeHasNoncontiguousAtomIndicesError` instead.
17+
1618

1719
### New features
1820

‎openff/toolkit/_tests/test_topology.py‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,7 @@
6666
MissingUniqueMoleculesError,
6767
MoleculeNotInTopologyError,
6868
NonUniqueSubstructureName,
69+
PDBMoleculeHasNoncontiguousAtomIndicesError,
6970
SubstructureAtomSmartsInvalid,
7071
SubstructureBondSmartsInvalid,
7172
UnassignedChemistryInPDBError,
@@ -906,6 +907,27 @@ def test_from_pdb_additional_substructures(self):
906907
expected_mol = Molecule.from_file(get_data_file_path("proteins/ace-ZZZ-gly-nme.sdf"))
907908
assert top.molecule(0).is_isomorphic_with(expected_mol, atom_stereochemistry_matching=False)
908909

910+
@requires_rdkit
911+
def test_from_pdb_molecule_atom_ordering_noncontiguous(self):
912+
"""
913+
Test that from_pdb correctly raises an error when the PDB atom ordering has atoms
914+
in the same molecule on noncontiguous lines.
915+
See https://github.com/openforcefield/openff-toolkit/issues/2093
916+
"""
917+
with pytest.raises(
918+
PDBMoleculeHasNoncontiguousAtomIndicesError, match="Atom indices 23 and 48 are in molecule 0"
919+
):
920+
Topology.from_pdb(get_data_file_path("proteins/split_chain.pdb"))
921+
922+
# Ensure that following the prior error message's advice and naively reordering the PDB lines works
923+
top = Topology.from_pdb(get_data_file_path("proteins/split_chain_reordered.pdb"))
924+
925+
# Check for sanity by ensuring that no bonds have unreasonable geometry
926+
for mol in top.molecules:
927+
for bond in mol.bonds:
928+
length = np.linalg.norm(mol.conformers[0][bond.atom1_index] - mol.conformers[0][bond.atom2_index])
929+
assert 0.4 < length.m_as(unit.angstrom) < 2.5
930+
909931
@requires_pkg("mdtraj")
910932
def test_from_mdtraj(self):
911933
"""Test construction of an OpenFF Topology from an MDTraj Topology object"""

‎openff/toolkit/data/proteins/split_chain.pdb‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
ATOM 1 C ACE 1 13.754 22.500 -5.565 1.00 0.00 C
2+
ATOM 2 O ACE 1 13.328 22.545 -4.403 1.00 0.00 O
3+
ATOM 3 CH3 ACE 1 13.637 23.697 -6.484 1.00 0.00 C
4+
ATOM 4 1HH3 ACE 1 14.165 24.540 -6.063 1.00 0.00 H
5+
ATOM 5 2HH3 ACE 1 14.063 23.465 -7.449 1.00 0.00 H
6+
ATOM 6 3HH3 ACE 1 12.599 23.963 -6.615 1.00 0.00 H
7+
ATOM 7 N CYS 2 14.305 21.418 -5.998 1.00 0.00 N
8+
ATOM 8 CA CYS 2 14.412 20.269 -5.103 1.00 0.00 C
9+
ATOM 9 C CYS 2 15.753 19.590 -5.252 1.00 0.00 C
10+
ATOM 10 O CYS 2 16.612 19.998 -6.047 1.00 0.00 O
11+
ATOM 11 CB CYS 2 13.235 19.317 -5.385 1.00 0.00 C
12+
ATOM 12 SG CYS 2 13.206 17.963 -4.187 1.00 0.00 S
13+
ATOM 13 H CYS 2 14.657 21.365 -6.943 1.00 0.00 H
14+
ATOM 14 HA CYS 2 14.343 20.634 -4.061 1.00 0.00 H
15+
ATOM 15 2HB CYS 2 13.288 18.898 -6.410 1.00 0.00 H
16+
ATOM 16 3HB CYS 2 12.267 19.846 -5.321 1.00 0.00 H
17+
ATOM 17 N NME 3 16.013 18.553 -4.529 1.00 0.00 N
18+
ATOM 18 CH3 NME 3 17.295 17.883 -4.655 1.00 0.00 C
19+
ATOM 19 H NME 3 15.324 18.208 -3.876 1.00 0.00 H
20+
ATOM 20 1HH3 NME 3 17.275 17.201 -5.492 1.00 0.00 H
21+
ATOM 21 2HH3 NME 3 18.079 18.610 -4.815 1.00 0.00 H
22+
ATOM 22 3HH3 NME 3 17.511 17.327 -3.755 1.00 0.00 H
23+
TER
24+
ATOM 23 C ACE A 1 14.483 12.299 -1.006 1.00 0.00 A C
25+
ATOM 24 O ACE A 1 15.094 11.306 -1.423 1.00 0.00 A O
26+
ATOM 25 CH3 ACE A 1 13.026 12.207 -0.607 1.00 0.00 A C
27+
ATOM 26 1HH3 ACE A 1 12.906 11.492 0.194 1.00 0.00 A H
28+
ATOM 27 2HH3 ACE A 1 12.673 13.171 -0.271 1.00 0.00 A H
29+
ATOM 28 3HH3 ACE A 1 12.431 11.891 -1.451 1.00 0.00 A H
30+
ATOM 29 N THR A 2 15.097 13.429 -0.910 1.00 0.00 A N
31+
ATOM 30 CA THR A 2 16.502 13.505 -1.298 1.00 0.00 A C
32+
ATOM 31 C THR A 2 17.286 14.359 -0.330 1.00 0.00 A C
33+
ATOM 32 O THR A 2 16.870 15.429 0.065 1.00 0.00 A O
34+
ATOM 33 CB THR A 2 16.641 14.055 -2.758 1.00 0.00 A C
35+
ATOM 34 CG2 THR A 2 18.078 14.199 -3.302 1.00 0.00 A C
36+
ATOM 35 OG1 THR A 2 15.996 13.187 -3.682 1.00 0.00 A O
37+
ATOM 36 H THR A 2 14.614 14.249 -0.572 1.00 0.00 A H
38+
ATOM 37 HA THR A 2 16.929 12.486 -1.261 1.00 0.00 A H
39+
ATOM 38 HB THR A 2 16.148 15.050 -2.805 1.00 0.00 A H
40+
ATOM 39 1HG THR A 2 16.057 13.617 -4.539 1.00 0.00 A H
41+
ATOM 40 1HG2 THR A 2 18.095 14.570 -4.343 1.00 0.00 A H
42+
ATOM 41 2HG2 THR A 2 18.681 14.919 -2.718 1.00 0.00 A H
43+
ATOM 42 3HG2 THR A 2 18.621 13.235 -3.288 1.00 0.00 A H
44+
ATOM 43 N NME A 3 18.431 13.938 0.089 1.00 0.00 A N
45+
ATOM 44 CH3 NME A 3 19.203 14.742 1.021 1.00 0.00 A C
46+
ATOM 45 H NME A 3 18.791 13.050 -0.230 1.00 0.00 A H
47+
ATOM 46 1HH3 NME A 3 19.761 15.495 0.485 1.00 0.00 A H
48+
ATOM 47 2HH3 NME A 3 18.543 15.229 1.724 1.00 0.00 A H
49+
ATOM 48 3HH3 NME A 3 19.894 14.115 1.565 1.00 0.00 A H
50+
TER
51+
ATOM 49 C ACE B 1 8.816 15.357 -0.581 1.00 0.00 B C
52+
ATOM 50 O ACE B 1 9.427 14.364 -0.998 1.00 0.00 B O
53+
ATOM 51 CH3 ACE B 1 7.359 15.265 -0.182 1.00 0.00 B C
54+
ATOM 52 1HH3 ACE B 1 7.239 14.550 0.619 1.00 0.00 B H
55+
ATOM 53 2HH3 ACE B 1 7.006 16.229 0.154 1.00 0.00 B H
56+
ATOM 54 3HH3 ACE B 1 6.764 14.949 -1.026 1.00 0.00 B H
57+
ATOM 55 N CYS B 2 9.431 16.487 -0.485 1.00 0.00 B N
58+
ATOM 56 CA CYS B 2 10.836 16.564 -0.873 1.00 0.00 B C
59+
ATOM 57 C CYS B 2 11.619 17.418 0.097 1.00 0.00 B C
60+
ATOM 58 O CYS B 2 11.194 18.510 0.501 1.00 0.00 B O
61+
ATOM 59 CB CYS B 2 10.917 17.089 -2.319 1.00 0.00 B C
62+
ATOM 60 SG CYS B 2 12.620 17.038 -2.924 1.00 0.00 B S
63+
ATOM 61 H CYS B 2 8.947 17.306 -0.147 1.00 0.00 B H
64+
ATOM 62 HA CYS B 2 11.264 15.545 -0.832 1.00 0.00 B H
65+
ATOM 63 2HB CYS B 2 10.532 18.126 -2.397 1.00 0.00 B H
66+
ATOM 64 3HB CYS B 2 10.297 16.481 -3.001 1.00 0.00 B H
67+
ATOM 65 N NME B 3 12.764 16.998 0.516 1.00 0.00 B N
68+
ATOM 66 CH3 NME B 3 13.535 17.802 1.449 1.00 0.00 B C
69+
ATOM 67 H NME B 3 13.125 16.111 0.197 1.00 0.00 B H
70+
ATOM 68 1HH3 NME B 3 14.093 18.556 0.913 1.00 0.00 B H
71+
ATOM 69 2HH3 NME B 3 12.874 18.289 2.152 1.00 0.00 B H
72+
ATOM 70 3HH3 NME B 3 14.226 17.176 1.993 1.00 0.00 B H
73+
TER
74+
END

‎openff/toolkit/data/proteins/split_chain_reordered.pdb‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
ATOM 1 C ACE 1 13.754 22.500 -5.565 1.00 0.00 C
2+
ATOM 2 O ACE 1 13.328 22.545 -4.403 1.00 0.00 O
3+
ATOM 3 CH3 ACE 1 13.637 23.697 -6.484 1.00 0.00 C
4+
ATOM 4 1HH3 ACE 1 14.165 24.540 -6.063 1.00 0.00 H
5+
ATOM 5 2HH3 ACE 1 14.063 23.465 -7.449 1.00 0.00 H
6+
ATOM 6 3HH3 ACE 1 12.599 23.963 -6.615 1.00 0.00 H
7+
ATOM 7 N CYS 2 14.305 21.418 -5.998 1.00 0.00 N
8+
ATOM 8 CA CYS 2 14.412 20.269 -5.103 1.00 0.00 C
9+
ATOM 9 C CYS 2 15.753 19.590 -5.252 1.00 0.00 C
10+
ATOM 10 O CYS 2 16.612 19.998 -6.047 1.00 0.00 O
11+
ATOM 11 CB CYS 2 13.235 19.317 -5.385 1.00 0.00 C
12+
ATOM 12 SG CYS 2 13.206 17.963 -4.187 1.00 0.00 S
13+
ATOM 13 H CYS 2 14.657 21.365 -6.943 1.00 0.00 H
14+
ATOM 14 HA CYS 2 14.343 20.634 -4.061 1.00 0.00 H
15+
ATOM 15 2HB CYS 2 13.288 18.898 -6.410 1.00 0.00 H
16+
ATOM 16 3HB CYS 2 12.267 19.846 -5.321 1.00 0.00 H
17+
ATOM 17 N NME 3 16.013 18.553 -4.529 1.00 0.00 N
18+
ATOM 18 CH3 NME 3 17.295 17.883 -4.655 1.00 0.00 C
19+
ATOM 19 H NME 3 15.324 18.208 -3.876 1.00 0.00 H
20+
ATOM 20 1HH3 NME 3 17.275 17.201 -5.492 1.00 0.00 H
21+
ATOM 21 2HH3 NME 3 18.079 18.610 -4.815 1.00 0.00 H
22+
ATOM 22 3HH3 NME 3 17.511 17.327 -3.755 1.00 0.00 H
23+
TER
24+
ATOM 49 C ACE B 1 8.816 15.357 -0.581 1.00 0.00 B C
25+
ATOM 50 O ACE B 1 9.427 14.364 -0.998 1.00 0.00 B O
26+
ATOM 51 CH3 ACE B 1 7.359 15.265 -0.182 1.00 0.00 B C
27+
ATOM 52 1HH3 ACE B 1 7.239 14.550 0.619 1.00 0.00 B H
28+
ATOM 53 2HH3 ACE B 1 7.006 16.229 0.154 1.00 0.00 B H
29+
ATOM 54 3HH3 ACE B 1 6.764 14.949 -1.026 1.00 0.00 B H
30+
ATOM 55 N CYS B 2 9.431 16.487 -0.485 1.00 0.00 B N
31+
ATOM 56 CA CYS B 2 10.836 16.564 -0.873 1.00 0.00 B C
32+
ATOM 57 C CYS B 2 11.619 17.418 0.097 1.00 0.00 B C
33+
ATOM 58 O CYS B 2 11.194 18.510 0.501 1.00 0.00 B O
34+
ATOM 59 CB CYS B 2 10.917 17.089 -2.319 1.00 0.00 B C
35+
ATOM 60 SG CYS B 2 12.620 17.038 -2.924 1.00 0.00 B S
36+
ATOM 61 H CYS B 2 8.947 17.306 -0.147 1.00 0.00 B H
37+
ATOM 62 HA CYS B 2 11.264 15.545 -0.832 1.00 0.00 B H
38+
ATOM 63 2HB CYS B 2 10.532 18.126 -2.397 1.00 0.00 B H
39+
ATOM 64 3HB CYS B 2 10.297 16.481 -3.001 1.00 0.00 B H
40+
ATOM 65 N NME B 3 12.764 16.998 0.516 1.00 0.00 B N
41+
ATOM 66 CH3 NME B 3 13.535 17.802 1.449 1.00 0.00 B C
42+
ATOM 67 H NME B 3 13.125 16.111 0.197 1.00 0.00 B H
43+
ATOM 68 1HH3 NME B 3 14.093 18.556 0.913 1.00 0.00 B H
44+
ATOM 69 2HH3 NME B 3 12.874 18.289 2.152 1.00 0.00 B H
45+
ATOM 70 3HH3 NME B 3 14.226 17.176 1.993 1.00 0.00 B H
46+
TER
47+
ATOM 23 C ACE A 1 14.483 12.299 -1.006 1.00 0.00 A C
48+
ATOM 24 O ACE A 1 15.094 11.306 -1.423 1.00 0.00 A O
49+
ATOM 25 CH3 ACE A 1 13.026 12.207 -0.607 1.00 0.00 A C
50+
ATOM 26 1HH3 ACE A 1 12.906 11.492 0.194 1.00 0.00 A H
51+
ATOM 27 2HH3 ACE A 1 12.673 13.171 -0.271 1.00 0.00 A H
52+
ATOM 28 3HH3 ACE A 1 12.431 11.891 -1.451 1.00 0.00 A H
53+
ATOM 29 N THR A 2 15.097 13.429 -0.910 1.00 0.00 A N
54+
ATOM 30 CA THR A 2 16.502 13.505 -1.298 1.00 0.00 A C
55+
ATOM 31 C THR A 2 17.286 14.359 -0.330 1.00 0.00 A C
56+
ATOM 32 O THR A 2 16.870 15.429 0.065 1.00 0.00 A O
57+
ATOM 33 CB THR A 2 16.641 14.055 -2.758 1.00 0.00 A C
58+
ATOM 34 CG2 THR A 2 18.078 14.199 -3.302 1.00 0.00 A C
59+
ATOM 35 OG1 THR A 2 15.996 13.187 -3.682 1.00 0.00 A O
60+
ATOM 36 H THR A 2 14.614 14.249 -0.572 1.00 0.00 A H
61+
ATOM 37 HA THR A 2 16.929 12.486 -1.261 1.00 0.00 A H
62+
ATOM 38 HB THR A 2 16.148 15.050 -2.805 1.00 0.00 A H
63+
ATOM 39 1HG THR A 2 16.057 13.617 -4.539 1.00 0.00 A H
64+
ATOM 40 1HG2 THR A 2 18.095 14.570 -4.343 1.00 0.00 A H
65+
ATOM 41 2HG2 THR A 2 18.681 14.919 -2.718 1.00 0.00 A H
66+
ATOM 42 3HG2 THR A 2 18.621 13.235 -3.288 1.00 0.00 A H
67+
ATOM 43 N NME A 3 18.431 13.938 0.089 1.00 0.00 A N
68+
ATOM 44 CH3 NME A 3 19.203 14.742 1.021 1.00 0.00 A C
69+
ATOM 45 H NME A 3 18.791 13.050 -0.230 1.00 0.00 A H
70+
ATOM 46 1HH3 NME A 3 19.761 15.495 0.485 1.00 0.00 A H
71+
ATOM 47 2HH3 NME A 3 18.543 15.229 1.724 1.00 0.00 A H
72+
ATOM 48 3HH3 NME A 3 19.894 14.115 1.565 1.00 0.00 A H
73+
TER
74+
END

‎openff/toolkit/utils/exceptions.py‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -406,6 +406,12 @@ class MultipleMoleculesInPDBError(OpenFFToolkitException):
406406
"""Error raised when a multiple molecules are found when one was expected"""
407407

408408

409+
class PDBMoleculeHasNoncontiguousAtomIndicesError(OpenFFToolkitException):
410+
"""Error raised when attempting to load a PDB file where at least one molecule has a gap in its atom
411+
indices (this can not be represented in an OpenFF Topology as atoms in a molecule must be
412+
contiguously indexed)"""
413+
414+
409415
class MultipleComponentsInMoleculeWarning(UserWarning):
410416
"""Warning emitted when user attempts to make an OpenFF Molecule with multiple disconnected components"""
411417

‎openff/toolkit/utils/rdkit_wrapper.py‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@
3939
MultipleComponentsInMoleculeWarning,
4040
NonUniqueSubstructureName,
4141
NotAttachedToMoleculeError,
42+
PDBMoleculeHasNoncontiguousAtomIndicesError,
4243
RadicalsNotSupportedError,
4344
SMILESParseError,
4445
SubstructureAtomSmartsInvalid,
@@ -315,6 +316,33 @@ def _polymer_openmm_pdbfile_to_offtop(
315316

316317
rdkit_mol = self._polymer_openmm_topology_to_rdmol(omm_top, substructure_dictionary)
317318

319+
# If the PDB file contains a single molecule that doesn't have contiguous atom indices,
320+
# we would need to rearrange the atoms away from the order that the user expects. We could
321+
# do this with some additional logic, but until a user makes a compelling case for this to
322+
# be supported, we'll just raise an error.
323+
sorted_mol_frags = [tuple(sorted(i)) for i in Chem.GetMolFrags(rdkit_mol)]
324+
325+
for mol_idx, frag_idxs in enumerate(sorted_mol_frags):
326+
# frag_idxs is sorted
327+
min_frag_idx = frag_idxs[0]
328+
max_frag_idx = frag_idxs[-1]
329+
if len(frag_idxs) != (max_frag_idx - min_frag_idx) + 1:
330+
expected_idxs = set(range(min_frag_idx, max_frag_idx + 1))
331+
missing_idxs = expected_idxs.difference(set(frag_idxs))
332+
raise PDBMoleculeHasNoncontiguousAtomIndicesError(
333+
"At least one molecule in the PDB file has noncontiguous atom indices. "
334+
"This is not currently supported by Topology.from_pdb.\n\n"
335+
f"Atom indices {min(missing_idxs) + 1} and {max(missing_idxs) + 1} "
336+
f"are in molecule {mol_idx}, but some or all of the intervening indices are not.\n\n"
337+
"This can happen when a crosslink is introduced between two molecules, but other "
338+
"molecules are present between them in the PDB atom ordering. "
339+
"You may be able to get around this error by rearranging the atom order in your PDB file to "
340+
"ensure that all the atoms in a single covalently bonded molecule are listed on a contiguous "
341+
"series of lines (don't change atom/residue/chain naming or numbering or any "
342+
"CONECT records, just the order of the ATOM/HETATM records). "
343+
"For more information see https://github.com/openforcefield/openff-toolkit/issues/2093"
344+
)
345+
318346
rdmol_conformer = Chem.Conformer()
319347
for atom_idx in range(rdkit_mol.GetNumAtoms()):
320348
x, y, z = coords_angstrom[atom_idx, :]
@@ -602,6 +630,7 @@ def _polymer_openmm_topology_to_rdmol(self, omm_top, substructure_library):
602630
# Get a tuple of tuples of atom indices belonging to separate molecules in this RDMol
603631
# (note that this rdmol may actually be a solvated protein-ligand system)
604632
sorted_mol_frags = [tuple(sorted(i)) for i in Chem.GetMolFrags(mol)]
633+
605634
query_number = 0
606635
for res_idx, res_name in enumerate(substructure_library):
607636
# TODO: This is a hack for the moment since we don't have a more sophisticated way to resolve clashes

0 commit comments

Comments
 (0)