Skip to content

P
pdbfixer
Commit 92044f3c by Peter Eastman Committed by GitHub
Options

Support nonstandard residues (#296) 

* Can add templates for new residues

* Added documentation on templates

* Bug fixes
parent 5cf78a2b
Showing with 318 additions and 49 deletions
Manual.html
......@@ -281,5 +281,34 @@ fixer.nonstandardResidues = [(residue, 'ALA') for residue, replacement in fixer.
<span style="color:grey">fixer.replaceNonstandardResidues()</span>
</pre></tt>
<h3>Nonstandard Residues</h3>
PDBFixer knows how to build all the standard amino acids and nucleotides. If you want to apply it to other residues, you need to tell it how. This is done by defining templates.
A template is a description of a residue: what atoms it contains, how they are bonded to each other, and what conformation they take on. For residues defined in the PDB Chemical Component Dictionary, you can simply tell PDBFixer to download the definition of whatever residue you need. For example, if you want to model a protein containing a selenocysteine (PDB code SEC), call
<tt><pre>
fixer.downloadTemplate('SEC')
</pre></tt>
Now you can call other methods such as <tt>findMissingAtoms()</tt>, <tt>addMissingAtoms()</tt>, and <tt>addMissingHydrogens()</tt> in the usual way. PDBFixer will fix problems in the SEC residue exactly as it does for standard ones.
<p/>
If your structure contains molecules or residues that are not present in the PDB Chemical Component Dictionary, you can still add templates for them, but you need to define them yourself. This is done by calling <tt>registerTemplate()</tt>. It takes two required arguments: an OpenMM Topology object describing the structure of the residue, and a list of positions describing a typical conformation of the residue. These can be easily obtained from, for example, a PDB or PDBx/mmCIF file containing just the single residue.
<tt><pre>
pdbx = PDBxFile('MOL.cif')
fixer.registerTemplate(pdbx.topology, pdbx.positions)
</pre></tt>
For residues that appear as part of larger molecules rather than in isolation, you need to provide one more piece of information: which atoms should only be present in terminal residues. For example, a C-terminal amino acid is capped with OXT and HXT atoms, but these atoms are omitted when the residue appears in the middle of a chain and is bonded to another residue. Likewise, the nitrogen in an N-terminal amino acid has two hydrogens attached to it (H and H2), but when it appears in the middle of a chain and the nitrogen is bonded to another residue, the H2 hydrogen is omitted. You can pass a third argument containing a list of boolean flags for each atom, specifying which ones should be omitted from non-terminal residues.
<tt><pre>
pdbx = PDBxFile('RES.cif')
terminal = [atom.name in ('H2', 'OXT', 'HXT') for atom in pdbx.topology.atoms()]
fixer.registerTemplate(pdbx.topology, pdbx.positions, terminal)
</pre></tt>
Strictly speaking, it is the bonds that matter, not the position in the chain. For example, the sulfur in a CYS residue has a hydrogen that must be omitted when it forms a disulfide bond to another residue. What matters is whether the sulfur has a bond connecting it to another residue, not the position of the CYS residue in its chain.
</body>
</html>
pdbfixer/pdbfixer.py
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2013-2023 Stanford University and the Authors.
Portions copyright (c) 2013-2024 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -54,6 +54,7 @@ import sys
import os
import os.path
import math
from collections import defaultdict
from pkg_resources import resource_filename
......@@ -98,6 +99,15 @@ class ModifiedResidue(object):
self.residueName = residueName
self.standardName = standardName
class Template:
"""Template represents a standard residue, or a nonstandard one registered with registerTemplate()."""
def __init__(self, topology, positions, terminal=None):
self.topology = topology
self.positions = positions
if terminal is None:
terminal = [False]*topology.getNumAtoms()
self.terminal = terminal
def _guessFileFormat(file, filename):
"""Guess whether a file is PDB or PDBx/mmCIF based on its filename and contents."""
filename = filename.lower()
......@@ -276,7 +286,7 @@ class PDBFixer(object):
for file in os.listdir(templatesPath):
templatePdb = app.PDBFile(os.path.join(templatesPath, file))
name = next(templatePdb.topology.residues()).name
self.templates[name] = templatePdb
self.templates[name] = Template(templatePdb.topology, templatePdb.positions)
def _initializeFromPDB(self, file):
"""Initialize this object by reading a PDB file."""
......@@ -344,6 +354,96 @@ class PDBFixer(object):
for row in modData.getRowList():
self.modifiedResidues.append(ModifiedResidue(row[asymIdCol], int(row[resNumCol]), row[resNameCol], row[standardResCol]))
def _getTemplate(self, name):
"""Return the template with a name. If none has been registered, this will return None."""
if name in self.templates:
return self.templates[name]
return None
def registerTemplate(self, topology, positions, terminal=None):
"""Register a template for a nonstandard residue. This allows PDBFixer to add missing residues of this type,
to add missing atoms to existing residues, and to mutate other residues to it.
Parameters
----------
topology: openmm.app.Topology
A Topology containing a single chain with a single residue, describing the nonstandard residue
being registered.
positions: array of shape (n_atoms, 3)
The positions of the atoms in the residue in a typical conformation. These positions are used
when adding missing atoms or residues.
terminal: optional list of bool
If this is present, it should be a list of length equal to the number of atoms in the residue.
If an element is True, that indicates the corresponding atom should only be added to terminal
residues.
"""
residues = list(topology.residues())
if len(residues) != 1:
raise ValueError('The Topology must contain a single residue')
if topology.getNumAtoms() != len(positions):
raise ValueError('The number of positions does not match the number of atoms in the Topology')
if terminal is not None and len(terminal) != topology.getNumAtoms():
raise ValueError('The number of terminal flags does not match the number of atoms in the Topology')
self.templates[residues[0].name] = Template(topology, positions, terminal)
def downloadTemplate(self, name):
"""Attempt to download a residue definition from the PDB and register a template for it.
Parameters
----------
name: str
The name of the residue, as specified in the PDB Chemical Component Dictionary.
Returns
-------
True if a template was successfully registered, false otherwise.
"""
name = name.upper()
try:
file = urlopen(f'https://files.rcsb.org/ligands/download/{name}.cif')
contents = file.read().decode('utf-8')
file.close()
except:
return False
# Load the atoms.
from openmm.app.internal.pdbx.reader.PdbxReader import PdbxReader
reader = PdbxReader(StringIO(contents))
data = []
reader.read(data)
block = data[0]
atomData = block.getObj('chem_comp_atom')
atomNameCol = atomData.getAttributeIndex('atom_id')
symbolCol = atomData.getAttributeIndex('type_symbol')
leavingCol = atomData.getAttributeIndex('pdbx_leaving_atom_flag')
xCol = atomData.getAttributeIndex('pdbx_model_Cartn_x_ideal')
yCol = atomData.getAttributeIndex('pdbx_model_Cartn_y_ideal')
zCol = atomData.getAttributeIndex('pdbx_model_Cartn_z_ideal')
topology = app.Topology()
chain = topology.addChain()
residue = topology.addResidue(name, chain)
positions = []
atomByName = {}
terminal = []
for row in atomData.getRowList():
atomName = row[atomNameCol]
atom = topology.addAtom(atomName, app.Element.getBySymbol(row[symbolCol]), residue)
atomByName[atomName] = atom
terminal.append(row[leavingCol] == 'Y')
positions.append(mm.Vec3(float(row[xCol]), float(row[yCol]), float(row[zCol]))*0.1)
positions = positions*unit.nanometers
# Load the bonds.
bondData = block.getObj('chem_comp_bond')
if bondData is not None:
atom1Col = bondData.getAttributeIndex('atom_id_1')
atom2Col = bondData.getAttributeIndex('atom_id_2')
for row in bondData.getRowList():
topology.addBond(atomByName[row[atom1Col]], atomByName[row[atom2Col]])
self.registerTemplate(topology, positions, terminal)
return True
def _addAtomsToTopology(self, heavyAtomsOnly, omitUnknownMolecules):
"""Create a new Topology in which missing atoms have been added.
......@@ -375,7 +475,7 @@ class PDBFixer(object):
addedOXT = []
residueCenters = [self._computeResidueCenter(res).value_in_unit(unit.nanometers) for res in self.topology.residues()]*unit.nanometers
for chain in self.topology.chains():
if omitUnknownMolecules and not any(residue.name in self.templates for residue in chain.residues()):
if omitUnknownMolecules and all(self._getTemplate(residue.name) is None for residue in chain.residues()):
continue
chainResidues = list(chain.residues())
newChain = newTopology.addChain(chain.id)
......@@ -413,7 +513,7 @@ class PDBFixer(object):
# Find corresponding atoms in the residue and the template.
template = self.templates[residue.name]
template = self._getTemplate(residue.name)
atomPositions = dict((atom.name, self.positions[atom.index]) for atom in residue.atoms())
points1 = []
points2 = []
......@@ -508,7 +608,7 @@ class PDBFixer(object):
# Add the residues.
for i, residueName in enumerate(residueNames):
template = self.templates[residueName]
template = self._getTemplate(residueName)
# Find a translation that best matches the adjacent residue.
......@@ -703,7 +803,7 @@ class PDBFixer(object):
replacement = modres[key]
if replacement == 'DU':
replacement = 'DT'
if replacement in self.templates:
if self._getTemplate(replacement) != None:
nonstandard[residue] = replacement
self.nonstandardResidues = [(r, nonstandard[r]) for r in sorted(nonstandard, key=lambda r: r.index)]
......@@ -731,7 +831,7 @@ class PDBFixer(object):
for residue, replaceWith in self.nonstandardResidues:
residue.name = replaceWith
template = self.templates[replaceWith]
template = self._getTemplate(replaceWith)
standardAtoms = set(atom.name for atom in template.topology.atoms())
for atom in residue.atoms():
if atom.element in (None, hydrogen) or atom.name not in standardAtoms:
......@@ -760,6 +860,10 @@ class PDBFixer(object):
Notes
-----
If a target residue is not a standard amino acid, and if no template
has been registered for it with registerTemplate(), this function
attempts to look it up from the PDB and create a new template for it.
We require three letter codes to avoid possible ambiguitities.
We can't guarantee that the resulting model is a good one; for
significant changes in sequence, you should probably be using
......@@ -800,9 +904,10 @@ class PDBFixer(object):
if residue.name != old_name:
raise(ValueError("You asked to mutate chain %s residue %d name %s, but that residue is actually %s!" % (chain_id, resSeq, old_name, residue.name)))
try:
template = self.templates[new_name]
except KeyError:
if self._getTemplate(new_name) is None:
# Try to download a template from the PDB.
self.downloadTemplate(new_name)
if self._getTemplate(new_name) is None:
raise(KeyError("Cannot find residue %s in template library!" % new_name))
# Store mutation
......@@ -816,7 +921,7 @@ class PDBFixer(object):
for residue in residue_map.keys():
replaceWith = residue_map[residue]
residue.name = replaceWith
template = self.templates[replaceWith]
template = self._getTemplate(replaceWith)
standardAtoms = set(atom.name for atom in template.topology.atoms())
for atom in residue.atoms():
if atom.element in (None, hydrogen) or atom.name not in standardAtoms:
......@@ -858,23 +963,60 @@ class PDBFixer(object):
missingAtoms = {}
missingTerminals = {}
# Determine which atoms have an external bond to another residue.
hasExternal = defaultdict(bool)
for atom1, atom2 in self.topology.bonds():
if atom1.residue != atom2.residue:
hasExternal[(atom1.residue, atom1.name)] = True
hasExternal[(atom2.residue, atom2.name)] = True
for chain in self.topology.chains():
chainResidues = list(chain.residues())
for residue in chain.residues():
atomNames = [atom.name for atom in residue.atoms()]
if all(name in atomNames for name in ['C', 'O', 'CA']):
# We'll be adding peptide bonds.
if residue != chainResidues[0]:
hasExternal[(residue, 'N')] = True
if residue != chainResidues[-1]:
hasExternal[(residue, 'C')] = True
# Loop over residues.
for chain in self.topology.chains():
nucleic = any(res.name in dnaResidues or res.name in rnaResidues for res in chain.residues())
chainResidues = list(chain.residues())
for residue in chain.residues():
if residue.name in self.templates:
template = self.templates[residue.name]
template = self._getTemplate(residue.name)
if template is not None:
# If an atom is marked as terminal only, and if it is bonded to any atom that has an external bond
# to another residue, we need to omit that atom and any other terminal-only atom bonded to it.
bondedTo = defaultdict(set)
for atom1, atom2 in template.topology.bonds():
bondedTo[atom1].add(atom2)
bondedTo[atom2].add(atom1)
skip = set()
for atom, terminal in zip(template.topology.atoms(), template.terminal):
if terminal:
for atom2 in bondedTo[atom]:
if hasExternal[(residue, atom2.name)]:
skip.add(atom)
for atom, terminal in zip(template.topology.atoms(), template.terminal):
if terminal:
for atom2 in bondedTo[atom]:
if atom2 in skip:
skip.add(atom)
atomNames = set(atom.name for atom in residue.atoms())
templateAtoms = list(template.topology.atoms())
if residue == chainResidues[0] and (chain.index, 0) not in self.missingResidues:
templateAtoms = [atom for atom in template.topology.atoms() if atom not in skip]
if nucleic and residue == chainResidues[0] and (chain.index, 0) not in self.missingResidues:
templateAtoms = [atom for atom in templateAtoms if atom.name not in ('P', 'OP1', 'OP2')]
# Add atoms from the template that are missing.
missing = []
for atom in templateAtoms:
if atom.name not in atomNames:
if atom.name not in atomNames and atom.element != app.element.hydrogen:
missing.append(atom)
if len(missing) > 0:
missingAtoms[residue] = missing
......@@ -1118,9 +1260,22 @@ class PDBFixer(object):
def _describeVariant(self, residue, definitions):
"""Build the variant description to pass to addHydrogens() for a residue."""
if residue.name not in definitions:
return None
if residue.name not in app.PDBFile._standardResidues and self._getTemplate(residue.name) is not None:
# The user has registered a template for this residue. Use the hydrogens from it.
template = self._getTemplate(residue.name)
atoms = [(atom.name, atom.element.symbol.upper(), terminal) for atom, terminal in zip(template.topology.atoms(), template.terminal)]
resAtoms = dict((atom.name, atom) for atom in residue.atoms())
bonds = []
for atom1, atom2 in template.topology.bonds():
if atom1.element == app.element.hydrogen and atom2.name in resAtoms:
bonds.append((atom1.name, atom2.name))
elif atom2.element == app.element.hydrogen and atom1.name in resAtoms:
bonds.append((atom2.name, atom1.name))
elif residue.name in definitions:
# We downloaded a definition.
atoms, bonds = definitions[residue.name]
else:
return None
# See if the heavy atoms are identical.
......
pdbfixer/tests/data/CAS.pdb 0 → 100644
ATOM 1 N CAS A 1 -1.377 -1.201 -2.417 1.00 10.00 N
ATOM 2 CA CAS A 1 -0.664 0.081 -2.349 1.00 10.00 C
ATOM 3 CB CAS A 1 0.084 0.179 -1.018 1.00 10.00 C
ATOM 4 C CAS A 1 0.319 0.169 -3.487 1.00 10.00 C
ATOM 5 O CAS A 1 0.824 -0.836 -3.928 1.00 10.00 O
ATOM 6 OXT CAS A 1 0.635 1.364 -4.010 1.00 10.00 O
ATOM 7 SG CAS A 1 -1.100 0.073 0.350 1.00 10.00 S
ATOM 8 AS CAS A 1 0.299 0.246 2.090 1.00 10.00 AS
ATOM 9 CE1 CAS A 1 -0.737 0.161 3.787 1.00 10.00 C
ATOM 10 CE2 CAS A 1 1.610 -1.249 2.038 1.00 10.00 C
ATOM 11 H CAS A 1 -0.678 -1.925 -2.343 1.00 10.00 H
ATOM 12 H2 CAS A 1 -1.945 -1.261 -1.585 1.00 10.00 H
ATOM 13 HA CAS A 1 -1.380 0.899 -2.424 1.00 10.00 H
ATOM 14 HB2 CAS A 1 0.613 1.131 -0.969 1.00 10.00 H
ATOM 15 HB3 CAS A 1 0.800 -0.638 -0.943 1.00 10.00 H
ATOM 16 HXT CAS A 1 1.267 1.420 -4.740 1.00 10.00 H
ATOM 17 HE11 CAS A 1 -0.056 0.245 4.633 1.00 10.00 H
ATOM 18 HE12 CAS A 1 -1.455 0.980 3.815 1.00 10.00 H
ATOM 19 HE13 CAS A 1 -1.268 -0.789 3.841 1.00 10.00 H
ATOM 20 HE21 CAS A 1 2.178 -1.202 1.109 1.00 10.00 H
ATOM 21 HE22 CAS A 1 2.291 -1.165 2.884 1.00 10.00 H
ATOM 22 HE23 CAS A 1 1.079 -2.200 2.092 1.00 10.00 H
CONECT 1 2 11 12
CONECT 2 1 3 4 13
CONECT 3 2 7 14 15
CONECT 4 2 5 6
CONECT 5 4
CONECT 6 4 16
CONECT 7 3 8
CONECT 8 7 9 10
CONECT 9 8 17 18 19
CONECT 10 8 20 21 22
CONECT 11 1
CONECT 12 1
CONECT 13 2
CONECT 14 3
CONECT 15 3
CONECT 16 6
CONECT 17 9
CONECT 18 9
CONECT 19 9
CONECT 20 10
CONECT 21 10
CONECT 22 10
END
pdbfixer/tests/test_add_hydrogens.py
import pdbfixer
import openmm.app as app
from pathlib import Path
from io import StringIO
def test_nonstandard():
"""Test adding hydrogens to nonstandard residues."""
content = (Path(__file__).parent / "data" / "4JSV.pdb").read_text()
fixer = pdbfixer.PDBFixer(pdbfile=StringIO(content))
fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "4JSV.pdb").as_posix())
fixer.removeChains(chainIndices=[0, 1, 2])
fixer.addMissingHydrogens()
for residue in fixer.topology.residues():
......@@ -17,10 +16,29 @@ def test_nonstandard():
def test_leaving_atoms():
"""Test adding hydrogens to a nonstandard residue with leaving atoms."""
content = (Path(__file__).parent / "data" / "1BHL.pdb").read_text()
fixer = pdbfixer.PDBFixer(pdbfile=StringIO(content))
fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "1BHL.pdb").as_posix())
fixer.addMissingHydrogens()
for residue in fixer.topology.residues():
count = sum(1 for atom in residue.atoms() if atom.element.symbol == 'H')
if residue.name == 'CAS':
assert count == 10
def test_registered_template():
"""Test adding hydrogens based on a template registered by the user."""
fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "1BHL.pdb").as_posix())
# Register a template for CAS from which a single hydrogen has been removed.
pdb = app.PDBFile((Path(__file__).parent / "data" / "CAS.pdb").as_posix())
modeller = app.Modeller(pdb.topology, pdb.positions)
modeller.delete([list(modeller.topology.atoms())[-1]])
terminal = [atom.name in ('H2', 'OXT', 'HXT') for atom in modeller.topology.atoms()]
fixer.registerTemplate(modeller.topology, modeller.positions, terminal)
# See if the correct hydrogens get added.
fixer.addMissingHydrogens()
for residue in fixer.topology.residues():
count = sum(1 for atom in residue.atoms() if atom.element.symbol == 'H')
if residue.name == 'CAS':
assert count == 9
pdbfixer/tests/test_build_and_simulate.py
from __future__ import print_function
import pdbfixer
import openmm
from pdbfixer.pdbfixer import PDBFixer
from openmm import app
import os
import os.path
......@@ -97,7 +97,7 @@ def test_build_and_simulate():
]
# DEBUG: A few small test cases.
pdbcodes_to_build = ['110D', '116D', '117D', '118D', '134D', '135D', '136D', '138D', '143D', '148D', '151D', '152D', '159D', '177D', '17RA', '183D', '184D', '186D', '187D', '188D', '189D', '1A11', '1A13', '1A1P', '1A3P', '1A51', '1A60', '1A83', '1A9L', '1AAF', '1AB1', '1ABZ', '1AC7', '1ACW', '1AD7', '1ADX', '1AFP', '1AFT', '1AFX', '1AG7', '1AGG', '1AGL', '1AGT', '1AHL', '1AIE', '1AJ1', '1AJF', '1AJJ', '1AJU', '1AKG', '1AKX', '1AL1', '1ALE', '1ALF', '1ALG', '1AM0', '1AMB', '1AMC', '1AML', '1ANP', '1ANR', '1ANS', '1AOO']
pdbcodes_to_build = ['110D', '116D', '117D', '118D', '134D', '135D', '136D', '138D', '143D', '148D', '151D', '152D', '177D', '17RA', '183D', '184D', '186D', '187D', '188D', '189D', '1A11', '1A13', '1A1P', '1A3P', '1A51', '1A60', '1A83', '1A9L', '1AAF', '1AB1', '1ABZ', '1AC7', '1ACW', '1AD7', '1ADX', '1AFP', '1AFT', '1AFX', '1AG7', '1AGG', '1AGL', '1AGT', '1AHL', '1AIE', '1AJ1', '1AJF', '1AJJ', '1AKG', '1AL1', '1ALE', '1ALF', '1ALG', '1AM0', '1AMB', '1AMC', '1AML', '1ANP', '1ANR', '1ANS', '1AOO', '1BH7', '1BX8', '1CEK']
# Don't simulate any.
pdbcodes_to_simulate = []
......@@ -105,20 +105,12 @@ def test_build_and_simulate():
# Keep track of list of failures.
failures = list()
forcefield = app.ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
for pdbcode in pdbcodes_to_build:
print("------------------------------------------------")
print(pdbcode)
output_pdb_filename = 'output.pdb'
# PDB setup parameters.
# TODO: Try several combinations?
from openmm import unit
pH = 7.0
ionic = 50.0 * unit.millimolar
box = 10.0 * unit.angstrom
positiveIon = 'Na+'
negativeIon = 'Cl-'
outfile = tempfile.NamedTemporaryFile(mode='w', delete=False)
output_pdb_filename = outfile.name
......@@ -127,10 +119,16 @@ def test_build_and_simulate():
watchdog = Watchdog(timeout_seconds)
build_successful = False
try:
from pdbfixer.pdbfixer import PDBFixer
from openmm import app
stage = "Creating PDBFixer..."
fixer = PDBFixer(pdbid=pdbcode)
stage = "Deleting hydrogens..."
if pdbcode in ['135D', '136D', '177D', '1A83', '1AGG', '1AJ1']:
# These input files include extra hydrogens that aren't supported by the force field.
# To avoid problems, delete all pre-existing hydrogens.
modeller = app.Modeller(fixer.topology, fixer.positions)
modeller.delete([a for a in fixer.topology.atoms() if a.element == app.element.hydrogen])
fixer.topology = modeller.topology
fixer.positions = modeller.positions
stage = "Finding missing residues..."
fixer.findMissingResidues()
stage = "Finding nonstandard residues..."
......@@ -147,6 +145,8 @@ def test_build_and_simulate():
fixer.addMissingHydrogens(pH)
stage = "Writing PDB file..."
app.PDBFile.writeFile(fixer.topology, fixer.positions, outfile)
stage = "Create System..."
forcefield.createSystem(fixer.topology)
stage = "Done."
outfile.close()
build_successful = True
......
pdbfixer/tests/test_mutate.py
import openmm.app as app
import pdbfixer
import tempfile
from pytest import raises
from pathlib import Path
def test_mutate_1():
fixer = pdbfixer.PDBFixer(pdbid='1VII')
......@@ -10,10 +10,6 @@ def test_mutate_1():
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
with tempfile.NamedTemporaryFile(mode='w+') as temp_pdb:
app.PDBFile.writeFile(fixer.topology, fixer.positions, temp_pdb)
temp_pdb.flush()
pdb = app.PDBFile(temp_pdb.name)
new_residue57 = list(fixer.topology.residues())[16]
assert new_residue57.name == "GLY", "Name of mutated residue did not change correctly!"
......@@ -31,7 +27,6 @@ def test_mutate_2():
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
temp_pdb = tempfile.NamedTemporaryFile(mode='w+')
new_residue57 = list(fixer.topology.residues())[16]
new_residue56 = list(fixer.topology.residues())[15]
assert new_residue57.name == "LEU", "Name of mutated residue did not change correctly!"
......@@ -56,15 +51,42 @@ def test_mutate_3_fails():
def test_mutate_4_fails():
with raises(KeyError):
fixer = pdbfixer.PDBFixer(pdbid='1VII')
fixer.applyMutations(["ALA-57-WTF", "SER-56-ALA"], "A")
def test_mutate_5_fails():
with raises(KeyError):
fixer = pdbfixer.PDBFixer(pdbid='1VII')
fixer.applyMutations(["ALA-1000-GLY", "SER-56-ALA"], "A")
def test_mutate_multiple_copies_of_chain_A():
fixer = pdbfixer.PDBFixer(pdbid='1OL5')
fixer.applyMutations(['TPO-287-THR','TPO-288-THR'], "A")
def test_mutate_to_nonstandard():
"""Test mutating to a nonstandard residue defined with registerTemplate()."""
fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "1BHL.pdb").as_posix())
pdb = app.PDBFile((Path(__file__).parent / "data" / "CAS.pdb").as_posix())
terminal = [atom.name in ('H2', 'OXT', 'HXT') for atom in pdb.topology.atoms()]
fixer.registerTemplate(pdb.topology, pdb.positions, terminal)
fixer.applyMutations(["SER-57-CAS", "ILE-60-CAS", "ASP-207-CAS"], "A")
fixer.findMissingResidues()
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
residues = list(fixer.topology.residues())
for i in (0, 3, 134):
assert residues[i].name == "CAS"
atoms = list(residues[i].atoms())
assert sum(1 for a in atoms if a.name == 'H2') == (1 if i == 0 else 0)
assert sum(1 for a in atoms if a.name == 'OXT') == (1 if i == 134 else 0)
assert sum(1 for a in atoms if a.name == 'HXT') == (1 if i == 134 else 0)
def test_download_template():
"""Test mutating to a nonstandard residue defined in the PDB."""
fixer = pdbfixer.PDBFixer(filename=(Path(__file__).parent / "data" / "1BHL.pdb").as_posix())
fixer.applyMutations(["SER-57-SEP", "ILE-60-SEP", "ASP-207-SEP"], "A")
fixer.findMissingResidues()
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(7.0)
residues = list(fixer.topology.residues())
for i in (0, 3, 134):
assert residues[i].name == "SEP"
atoms = list(residues[i].atoms())
assert sum(1 for a in atoms if a.element == app.element.phosphorus) == 1
assert sum(1 for a in atoms if a.name == 'OXT') == (1 if i == 134 else 0)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment