Skip to content

P
pdbfixer
Commit 816500da by peastman
Options

Merge pull request #60 from peastman/master 

Initial positions of added residues are more intelligently chosen
parents 6d75bb33 feff9522
Showing with 40 additions and 8 deletions
pdbfixer/pdbfixer.py
...@@ -115,6 +115,20 @@ def _overlayPoints(points1, points2): ...@@ -115,6 +115,20 @@ def _overlayPoints(points1, points2):
(u, s, v) = lin.svd(R) (u, s, v) = lin.svd(R)
return (-1*center2, np.dot(u, v).transpose(), center1) return (-1*center2, np.dot(u, v).transpose(), center1)
def _findUnoccupiedDirection(point, positions):
"""Given a point in space and a list of atom positions, find the direction in which the local density of atoms is lowest."""
point = point.value_in_unit(unit.nanometers)
direction = mm.Vec3(0, 0, 0)
for pos in positions.value_in_unit(unit.nanometers):
delta = pos-point
distance = unit.norm(delta)
if distance > 0.1:
distance2 = distance*distance
direction -= delta/(distance2*distance2)
direction /= unit.norm(direction)
return direction
class PDBFixer(object): class PDBFixer(object):
"""PDBFixer implements many tools for fixing problems in PDB files. """PDBFixer implements many tools for fixing problems in PDB files.
""" """
...@@ -250,6 +264,7 @@ class PDBFixer(object): ...@@ -250,6 +264,7 @@ class PDBFixer(object):
existingAtomMap = {} existingAtomMap = {}
addedAtomMap = {} addedAtomMap = {}
addedOXT = [] addedOXT = []
residueCenters = [self._computeResidueCenter(res).value_in_unit(unit.nanometers) for res in self.topology.residues()]*unit.nanometers
for chain in self.topology.chains(): for chain in self.topology.chains():
if omitUnknownMolecules and not any(residue.name in self.templates for residue in chain.residues()): if omitUnknownMolecules and not any(residue.name in self.templates for residue in chain.residues()):
continue continue
...@@ -264,13 +279,15 @@ class PDBFixer(object): ...@@ -264,13 +279,15 @@ class PDBFixer(object):
endPosition = self._computeResidueCenter(residue) endPosition = self._computeResidueCenter(residue)
if indexInChain > 0: if indexInChain > 0:
startPosition = self._computeResidueCenter(chainResidues[indexInChain-1]) startPosition = self._computeResidueCenter(chainResidues[indexInChain-1])
loopDirection = _findUnoccupiedDirection((startPosition+endPosition)/2, residueCenters)
else: else:
outward = endPosition-self.centroid outward = _findUnoccupiedDirection(endPosition, residueCenters)*unit.nanometers
norm = unit.norm(outward) norm = unit.norm(outward)
if norm > 0*unit.nanometer: if norm > 0*unit.nanometer:
outward *= len(insertHere)*0.5*unit.nanometer/norm outward *= len(insertHere)*0.5*unit.nanometer/norm
startPosition = endPosition+outward startPosition = endPosition+outward
self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, residue, newAtoms, newPositions) loopDirection = None
self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, loopDirection, residue, newAtoms, newPositions)
# Create the new residue and add existing heavy atoms. # Create the new residue and add existing heavy atoms.
...@@ -320,12 +337,12 @@ class PDBFixer(object): ...@@ -320,12 +337,12 @@ class PDBFixer(object):
insertHere = self.missingResidues[(chain.index, indexInChain+1)] insertHere = self.missingResidues[(chain.index, indexInChain+1)]
if len(insertHere) > 0: if len(insertHere) > 0:
startPosition = self._computeResidueCenter(residue) startPosition = self._computeResidueCenter(residue)
outward = startPosition-self.centroid outward = _findUnoccupiedDirection(startPosition, residueCenters)*unit.nanometers
norm = unit.norm(outward) norm = unit.norm(outward)
if norm > 0*unit.nanometer: if norm > 0*unit.nanometer:
outward *= len(insertHere)*0.5*unit.nanometer/norm outward *= len(insertHere)*0.5*unit.nanometer/norm
endPosition = startPosition+outward endPosition = startPosition+outward
self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, residue, newAtoms, newPositions) self._addMissingResiduesToChain(newChain, insertHere, startPosition, endPosition, None, residue, newAtoms, newPositions)
newResidue = list(newChain.residues())[-1] newResidue = list(newChain.residues())[-1]
if newResidue.name in proteinResidues: if newResidue.name in proteinResidues:
terminalsToAdd = ['OXT'] terminalsToAdd = ['OXT']
...@@ -357,9 +374,23 @@ class PDBFixer(object): ...@@ -357,9 +374,23 @@ class PDBFixer(object):
"""Compute the centroid of a residue.""" """Compute the centroid of a residue."""
return unit.sum([self.pdb.positions[atom.index] for atom in residue.atoms()])/len(list(residue.atoms())) return unit.sum([self.pdb.positions[atom.index] for atom in residue.atoms()])/len(list(residue.atoms()))
def _addMissingResiduesToChain(self, chain, residueNames, startPosition, endPosition, orientTo, newAtoms, newPositions): def _addMissingResiduesToChain(self, chain, residueNames, startPosition, endPosition, loopDirection, orientTo, newAtoms, newPositions):
"""Add a series of residues to a chain.""" """Add a series of residues to a chain."""
orientToPositions = dict((atom.name, self.positions[atom.index]) for atom in orientTo.atoms()) orientToPositions = dict((atom.name, self.positions[atom.index]) for atom in orientTo.atoms())
if loopDirection is None:
loopDirection = mm.Vec3(0, 0, 0)
# We'll add the residues in an arc connecting the endpoints. Figure out the height of that arc.
length = unit.norm(endPosition-startPosition)
numResidues = len(residueNames)
if length > numResidues*0.3*unit.nanometers:
loopHeight = 0*unit.nanometers
else:
loopHeight = (numResidues*0.3*unit.nanometers-length)/2
# Add the residues.
for i, residueName in enumerate(residueNames): for i, residueName in enumerate(residueNames):
template = self.templates[residueName] template = self.templates[residueName]
...@@ -376,7 +407,8 @@ class PDBFixer(object): ...@@ -376,7 +407,8 @@ class PDBFixer(object):
# Create the new residue. # Create the new residue.
newResidue = chain.topology.addResidue(residueName, chain) newResidue = chain.topology.addResidue(residueName, chain)
translate = startPosition+(endPosition-startPosition)*(i+1.0)/(len(residueNames)+1.0) fraction = (i+1.0)/(numResidues+1.0)
translate = startPosition + (endPosition-startPosition)*fraction + loopHeight*math.sin(fraction*math.pi)*loopDirection
templateAtoms = list(template.topology.atoms()) templateAtoms = list(template.topology.atoms())
if newResidue == next(chain.residues()): if newResidue == next(chain.residues()):
templateAtoms = [atom for atom in templateAtoms if atom.name not in ('P', 'OP1', 'OP2')] templateAtoms = [atom for atom in templateAtoms if atom.name not in ('P', 'OP1', 'OP2')]
...@@ -782,7 +814,7 @@ class PDBFixer(object): ...@@ -782,7 +814,7 @@ class PDBFixer(object):
mm.LocalEnergyMinimizer.minimize(context) mm.LocalEnergyMinimizer.minimize(context)
state = context.getState(getPositions=True) state = context.getState(getPositions=True)
nearest = self._findNearestDistance(context, newTopology, newAtoms) nearest = self._findNearestDistance(context, newTopology, newAtoms)
if nearest < 0.15: if nearest < 0.13:
# Some atoms are very close together. Run some dynamics while slowly increasing the strength of the # Some atoms are very close together. Run some dynamics while slowly increasing the strength of the
# repulsive interaction to try to improve the result. # repulsive interaction to try to improve the result.
...@@ -794,7 +826,7 @@ class PDBFixer(object): ...@@ -794,7 +826,7 @@ class PDBFixer(object):
if d > nearest: if d > nearest:
nearest = d nearest = d
state = context.getState(getPositions=True) state = context.getState(getPositions=True)
if nearest >= 0.15: if nearest >= 0.13:
break break
context.setState(state) context.setState(state)
state = context.getState(getPositions=True) state = context.getState(getPositions=True)
......
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