omic/propka
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

De-lint protonate.py.

Browse Source
This commit is contained in:
Nathan Baker
2020-05-25 13:05:16 -07:00
parent b3e2685e94
commit cd0e9e5c3d
1 changed files with 236 additions and 318 deletions
Download Patch File Download Diff File Expand all files Collapse all files

554
propka/protonate.py
View File

@@ -1,441 +1,359 @@
#!/usr/bin/python
"""Protonate a structure."""
import math
import propka.bonds
import propka.atom
from propka.vector_algebra import rotate_vector_around_an_axis, vector
from propka.lib import warning, debug
from __future__ import division
from __future__ import print_function
from propka.vector_algebra import *
import propka.bonds, propka.pdb, propka.atom
from propka.lib import info, warning, debug
class Protonate:
""" Protonates atoms using VSEPR theory """
def __init__(self, verbose=False):
self.verbose=verbose
self.valence_electrons = {'H': 1,
'He':2,
'Li':1,
'Be':2,
'B': 3,
'C': 4,
'N': 5,
'O': 6,
'F': 7,
'Ne':8,
'Na':1,
'Mg':2,
'Al':3,
'Si':4,
'P': 5,
'S': 6,
'Cl':7,
'Ar':8,
'K': 1,
'Ca':2,
'Sc':2,
'Ti':2,
'Va':2,
'Cr':1,
'Mn':2,
'Fe':2,
'Co':2,
'Ni':2,
'Cu':1,
'Zn':2,
'Ga':3,
'Ge':4,
'As':5,
'Se':6,
'Br':7,
'Kr':8,
'I':7,
}
self.standard_charges= {'ARG-NH1':1.0,
'ASP-OD2':-1.0,
'GLU-OE2':-1.0,
'HIS-ND1':1.0,
'LYS-NZ':1.0,
'N+':1.0,
'C-':-1.0}
self.sybyl_charges = {'N.pl3':+1,
'N.3':+1,
'N.4':+1,
'N.ar':+1,
'O.co2-':-1}
self.bond_lengths = {'C':1.09,
'N':1.01,
'O':0.96,
'F':0.92,
'Cl':1.27,
'Br':1.41,
'I':1.61,
'S':1.35}
# protonation_methods[steric_number] = method
self.protonation_methods = {4:self.tetrahedral,
3:self.trigonal}
return
"""Initialize with flag for verbosity
Args:
verbose: True for verbose output
"""
self.verbose = verbose
self.valence_electrons = {'H': 1, 'He': 2, 'Li': 1, 'Be': 2, 'B': 3,
'C': 4, 'N': 5, 'O': 6, 'F': 7, 'Ne': 8,
'Na': 1, 'Mg': 2, 'Al': 3, 'Si': 4, 'P': 5,
'S': 6, 'Cl': 7, 'Ar': 8, 'K': 1, 'Ca': 2,
'Sc': 2, 'Ti': 2, 'Va': 2, 'Cr': 1, 'Mn': 2,
'Fe': 2, 'Co': 2, 'Ni': 2, 'Cu': 1, 'Zn': 2,
'Ga': 3, 'Ge': 4, 'As': 5, 'Se': 6, 'Br': 7,
'Kr': 8, 'I': 7}
# TODO - consider putting charges in a configuration file
self.standard_charges = {'ARG-NH1': 1.0, 'ASP-OD2': -1.0,
'GLU-OE2': -1.0, 'HIS-ND1': 1.0,
'LYS-NZ': 1.0, 'N+': 1.0, 'C-': -1.0}
self.sybyl_charges = {'N.pl3': 1, 'N.3': 1, 'N.4': 1, 'N.ar': 1,
'O.co2-': 1}
# TODO - consider putting bond lengths in a configuration file
self.bond_lengths = {'C': 1.09, 'N': 1.01, 'O': 0.96, 'F': 0.92,
'Cl': 1.27, 'Br': 1.41, 'I': 1.61, 'S': 1.35}
self.protonation_methods = {4: self.tetrahedral, 3: self.trigonal}
def protonate(self, molecules):
""" Will protonate all atoms in the molecular container """
"""Protonate all atoms in the molecular container.
Args:
molecules: molecular containers
"""
debug('----- Protonation started -----')
# Remove all currently present hydrogen atoms
self.remove_all_hydrogen_atoms(molecules)
# protonate all atoms
for name in molecules.conformation_names:
non_H_atoms = molecules.conformations[name].get_non_hydrogen_atoms()
for atom in non_H_atoms:
non_h_atoms = (molecules.conformations[name]
.get_non_hydrogen_atoms())
for atom in non_h_atoms:
self.protonate_atom(atom)
# fix hydrogen names
#self.set_proton_names(non_H_atoms)
@staticmethod
def remove_all_hydrogen_atoms(molecular_container):
"""Remove all hydrogen atoms from molecule.
return
def remove_all_hydrogen_atoms(self, molecular_container):
Args:
molecular_container: molecule to remove hydrogens from
"""
for name in molecular_container.conformation_names:
molecular_container.conformations[name].atoms = molecular_container.conformations[name].get_non_hydrogen_atoms()
return
molecular_container.conformations[name].atoms = (
molecular_container.conformations[name]
.get_non_hydrogen_atoms())
def set_charge(self, atom):
"""Set charge for atom.
Args:
atom: atom to be charged
"""
# atom is a protein atom
if atom.type=='atom':
key = '%3s-%s'%(atom.res_name, atom.name)
if atom.type == 'atom':
key = '%3s-%s' % (atom.res_name, atom.name)
if atom.terminal:
debug(atom.terminal)
key=atom.terminal
if key in list(self.standard_charges.keys()):
key = atom.terminal
if key in self.standard_charges:
atom.charge = self.standard_charges[key]
debug('Charge', atom, atom.charge)
atom.charge_set = True
# atom is a ligand atom
elif atom.type=='hetatm':
if atom.sybyl_type in list(self.sybyl_charges.keys()):
elif atom.type == 'hetatm':
if atom.sybyl_type in self.sybyl_charges:
atom.charge = self.sybyl_charges[atom.sybyl_type]
atom.sybyl_type = atom.sybyl_type.replace('-','')
atom.sybyl_type = atom.sybyl_type.replace('-', '')
atom.charge_set = True
return
def protonate_atom(self, atom):
if atom.is_protonated: return
if atom.element == 'H': return
"""Protonate an atom.
Args:
atom: atom to be protonated
"""
if atom.is_protonated:
return
if atom.element == 'H':
return
self.set_charge(atom)
self.set_number_of_protons_to_add(atom)
self.set_steric_number_and_lone_pairs(atom)
self.add_protons(atom)
atom.is_protonated = True
return
def set_proton_names(self, heavy_atoms):
@staticmethod
def set_proton_names(heavy_atoms):
"""Set names for protons.
Args:
heavy_atoms: list of heavy atoms with protons to be named
"""
for heavy_atom in heavy_atoms:
i = 1
for bonded in heavy_atom.bonded_atoms:
if bonded.element == 'H':
bonded.name+='%d'%i
i+=1
return
bonded.name += '%d' % i
i += 1
def set_number_of_protons_to_add(self, atom):
debug('*'*10)
debug('Setting number of protons to add for',atom)
atom.number_of_protons_to_add = 8
debug(' %4d'%8)
atom.number_of_protons_to_add -= self.valence_electrons[atom.element]
debug('Valence eletrons: %4d'%-self.valence_electrons[atom.element])
atom.number_of_protons_to_add -= len(atom.bonded_atoms)
debug('Number of bonds: %4d'%- len(atom.bonded_atoms))
atom.number_of_protons_to_add -= atom.num_pi_elec_2_3_bonds
debug('Pi electrons: %4d'%-atom.num_pi_elec_2_3_bonds)
atom.number_of_protons_to_add += int(atom.charge)
debug('Charge: %4.1f'%atom.charge)
"""Set the number of protons to add to this atom.
Args:
atom: atom for calculation
"""
debug('*'*10)
debug('Setting number of protons to add for', atom)
atom.number_of_protons_to_add = 8
debug(' %4d' % 8)
atom.number_of_protons_to_add -= self.valence_electrons[atom.element]
debug('Valence eletrons: %4d' % -self.valence_electrons[atom.element])
atom.number_of_protons_to_add -= len(atom.bonded_atoms)
debug('Number of bonds: %4d' % -len(atom.bonded_atoms))
atom.number_of_protons_to_add -= atom.num_pi_elec_2_3_bonds
debug('Pi electrons: %4d' % -atom.num_pi_elec_2_3_bonds)
atom.number_of_protons_to_add += int(atom.charge)
debug('Charge: %4.1f' % atom.charge)
debug('-'*10)
debug(atom.number_of_protons_to_add)
return
def set_steric_number_and_lone_pairs(self, atom):
"""Set steric number and lone pairs for atom.
Args:
atom: atom for calculation
"""
# If we already did this, there is no reason to do it again
if atom.steric_num_lone_pairs_set:
return
debug('='*10)
debug('Setting steric number and lone pairs for',atom)
# costumly set the N backbone atoms up for peptide bond trigonal planer shape
#if atom.name == 'N' and len(atom.bonded_atoms) == 2:
# atom.steric_number = 3
# atom.number_of_lone_pairs = 0
# self.display 'Peptide bond: steric number is %d and number of lone pairs is %s'%(atom.steric_number,
# atom.number_of_lone_pairs)
# return
debug('Setting steric number and lone pairs for', atom)
atom.steric_number = 0
debug('%65s: %4d'%('Valence electrons',self.valence_electrons[atom.element]))
debug('%65s: %4d' % ('Valence electrons',
self.valence_electrons[atom.element]))
atom.steric_number += self.valence_electrons[atom.element]
debug('%65s: %4d'%('Number of bonds',len(atom.bonded_atoms)))
debug('%65s: %4d' % ('Number of bonds',
len(atom.bonded_atoms)))
atom.steric_number += len(atom.bonded_atoms)
debug('%65s: %4d'%('Number of hydrogen atoms to add',atom.number_of_protons_to_add))
debug('%65s: %4d' % ('Number of hydrogen atoms to add',
atom.number_of_protons_to_add))
atom.steric_number += atom.number_of_protons_to_add
debug('%65s: %4d'%('Number of pi-electrons in double and triple bonds(-)',atom.num_pi_elec_2_3_bonds))
debug('%65s: %4d' % ('Number of pi-electrons in double '
'and triple bonds(-)',
atom.num_pi_elec_2_3_bonds))
atom.steric_number -= atom.num_pi_elec_2_3_bonds
debug('%65s: %4d'%('Number of pi-electrons in conjugated double and triple bonds(-)',atom.num_pi_elec_conj_2_3_bonds))
debug('%65s: %4d' % ('Number of pi-electrons in conjugated double and '
'triple bonds(-)',
atom.num_pi_elec_conj_2_3_bonds))
atom.steric_number -= atom.num_pi_elec_conj_2_3_bonds
debug('%65s: %4d'%('Number of donated co-ordinated bonds',0))
debug('%65s: %4d' % ('Number of donated co-ordinated bonds', 0))
atom.steric_number += 0
debug('%65s: %4.1f'%('Charge(-)',atom.charge))
debug('%65s: %4.1f' % ('Charge(-)', atom.charge))
atom.steric_number -= atom.charge
atom.steric_number = math.floor(atom.steric_number/2.0)
atom.number_of_lone_pairs = atom.steric_number - len(atom.bonded_atoms) - atom.number_of_protons_to_add
atom.number_of_lone_pairs = (atom.steric_number
- len(atom.bonded_atoms)
- atom.number_of_protons_to_add)
debug('-'*70)
debug('%65s: %4d'%('Steric number',atom.steric_number))
debug('%65s: %4d'%('Number of lone pairs',atom.number_of_lone_pairs))
debug('%65s: %4d' % ('Steric number', atom.steric_number))
debug('%65s: %4d' % ('Number of lone pairs',
atom.number_of_lone_pairs))
atom.steric_num_lone_pairs_set = True
return
def add_protons(self, atom):
"""Add protons to atom.
Args:
atom: atom for calculation
"""
# decide which method to use
debug('PROTONATING',atom)
debug('PROTONATING', atom)
if atom.steric_number in list(self.protonation_methods.keys()):
self.protonation_methods[atom.steric_number](atom)
else:
warning('Do not have a method for protonating', atom, '(steric number: %d)' % atom.steric_number)
return
warning('Do not have a method for protonating',
atom, '(steric number: %d)' % atom.steric_number)
def trigonal(self, atom):
debug('TRIGONAL - %d bonded atoms'%(len(atom.bonded_atoms)))
"""Add hydrogens in trigonal geometry.
Args:
atom: atom to protonate
"""
debug('TRIGONAL - %d bonded atoms' % len(atom.bonded_atoms))
rot_angle = math.radians(120.0)
c = vector(atom1 = atom)
cvec = vector(atom1=atom)
# 0 bonds
if len(atom.bonded_atoms) == 0:
pass
# 1 bond
if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first one
a = vector(atom1 = atom, atom2 = atom.bonded_atoms[0])
avec = vector(atom1=atom, atom2=atom.bonded_atoms[0])
# use plane of bonded trigonal atom - e.g. arg
self.set_steric_number_and_lone_pairs(atom.bonded_atoms[0])
if atom.bonded_atoms[0].steric_number == 3 and len(atom.bonded_atoms[0].bonded_atoms)>1:
# use other atoms bonded to the neighbour to establish the plane, if possible
if (atom.bonded_atoms[0].steric_number == 3
and len(atom.bonded_atoms[0].bonded_atoms) > 1):
# use other atoms bonded to the neighbour to establish the
# plane, if possible
other_atom_indices = []
for i in range(len(atom.bonded_atoms[0].bonded_atoms)):
if atom.bonded_atoms[0].bonded_atoms[i] != atom:
for i, bonded_atom in enumerate(atom.bonded_atoms[0].bonded_atoms):
if bonded_atom != atom:
other_atom_indices.append(i)
v1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0])
v2 = vector(atom1 = atom.bonded_atoms[0],
atom2 = atom.bonded_atoms[0].bonded_atoms[other_atom_indices[0]])
axis = v1**v2
# this is a trick to make sure that the order of atoms doesn't influence
# the final postions of added protons
if len(other_atom_indices)>1:
v3 = vector(atom1 = atom.bonded_atoms[0],
atom2 = atom.bonded_atoms[0].bonded_atoms[other_atom_indices[1]])
axis2 = v1**v3
if axis * axis2>0:
vec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0])
vec2 = vector(atom1=atom.bonded_atoms[0],
atom2=atom.bonded_atoms[0]
.bonded_atoms[other_atom_indices[0]])
axis = vec1**vec2
# this is a trick to make sure that the order of atoms doesn't
# influence the final postions of added protons
if len(other_atom_indices) > 1:
vec3 = vector(atom1=atom.bonded_atoms[0],
atom2=atom.bonded_atoms[0]
.bonded_atoms[other_atom_indices[1]])
axis2 = vec1**vec3
if axis*axis2 > 0:
axis = axis+axis2
else:
axis = axis-axis2
else:
axis = a.orthogonal()
a = rotate_vector_around_an_axis(rot_angle, axis, a)
a = self.set_bond_distance(a, atom.element)
self.add_proton(atom, c+a)
axis = avec.orthogonal()
avec = rotate_vector_around_an_axis(rot_angle, axis, avec)
avec = self.set_bond_distance(avec, atom.element)
self.add_proton(atom, cvec+avec)
# 2 bonds
if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first two
a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0)
a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0)
avec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
new_a = -a1 - a2
new_a = -avec1 - avec2
new_a = self.set_bond_distance(new_a, atom.element)
self.add_proton(atom, c+new_a)
return
self.add_proton(atom, cvec+new_a)
def tetrahedral(self, atom):
debug('TETRAHEDRAL - %d bonded atoms'%(len(atom.bonded_atoms)))
"""Protonate atom in tetrahedral geometry.
Args:
atom: atom to protonate.
"""
debug('TETRAHEDRAL - %d bonded atoms' % len(atom.bonded_atoms))
# TODO - might be good to move tetrahedral angle to constant
rot_angle = math.radians(109.5)
# sanity check
# if atom.number_of_protons_to_add + len(atom.bonded_atoms) != 4:
# self.display 'Error: Attempting tetrahedral structure with %d bonds'%(atom.number_of_protons_to_add +
# len(atom.bonded_atoms))
c = vector(atom1 = atom)
cvec = vector(atom1=atom)
# 0 bonds
if len(atom.bonded_atoms) == 0:
pass
# 1 bond
if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first one
a = vector(atom1 = atom, atom2 = atom.bonded_atoms[0])
axis = a.orthogonal()
a = rotate_vector_around_an_axis(rot_angle, axis, a)
a = self.set_bond_distance(a, atom.element)
self.add_proton(atom, c+a)
avec = vector(atom1=atom, atom2=atom.bonded_atoms[0])
axis = avec.orthogonal()
avec = rotate_vector_around_an_axis(rot_angle, axis, avec)
avec = self.set_bond_distance(avec, atom.element)
self.add_proton(atom, cvec+avec)
# 2 bonds
if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first two
a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0)
a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0)
axis = a1 + a2
new_a = rotate_vector_around_an_axis(math.radians(90), axis, -a1)
avec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
axis = avec1 + avec2
new_a = rotate_vector_around_an_axis(math.radians(90), axis,
-avec1)
new_a = self.set_bond_distance(new_a, atom.element)
self.add_proton(atom, c+new_a)
self.add_proton(atom, cvec+new_a)
# 3 bonds
if len(atom.bonded_atoms) == 3 and atom.number_of_protons_to_add > 0:
a1 = vector(atom1 = atom, atom2 = atom.bonded_atoms[0]).rescale(1.0)
a2 = vector(atom1 = atom, atom2 = atom.bonded_atoms[1]).rescale(1.0)
a3 = vector(atom1 = atom, atom2 = atom.bonded_atoms[2]).rescale(1.0)
new_a = -a1-a2-a3
avec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
avec3 = vector(atom1=atom, atom2=atom.bonded_atoms[2]).rescale(1.0)
new_a = -avec1-avec2-avec3
new_a = self.set_bond_distance(new_a, atom.element)
self.add_proton(atom, c+new_a)
self.add_proton(atom, cvec+new_a)
return
@staticmethod
def add_proton(atom, position):
"""Add a proton to an atom at a specific position.
def add_proton(self, atom, position):
Args:
atom: atom to protonate
position: position for proton
"""
# Create the new proton
new_H = propka.atom.Atom()
new_H.set_property(numb = None,
name = 'H%s'%atom.name[1:],
res_name = atom.res_name,
chain_id = atom.chain_id,
res_num = atom.res_num,
x = round(position.x,3), # round of to three digimal points
y = round(position.y,3), # to avoid round-off differences
z = round(position.z,3), # when input file
occ = None,
beta = None)
new_H.element = 'H'
new_H.type = atom.type
new_H.bonded_atoms = [atom]
new_H.charge = 0
new_H.steric_number = 0
new_H.number_of_lone_pairs = 0
new_H.number_of_protons_to_add = 0
new_H.num_pi_elec_2_3_bonds = 0
new_H.is_protonates = True
atom.bonded_atoms.append(new_H)
atom.number_of_protons_to_add -=1
atom.conformation_container.add_atom(new_H)
new_h = propka.atom.Atom()
new_h.set_property(
numb=None,
name='H%s' % atom.name[1:],
res_name=atom.res_name,
chain_id=atom.chain_id,
res_num=atom.res_num,
x=round(position.x, 3), # round of to three decimal points to
# avoid round-off differences in input file
y=round(position.y, 3),
z=round(position.z, 3),
occ=None,
beta=None)
new_h.element = 'H'
new_h.type = atom.type
new_h.bonded_atoms = [atom]
new_h.charge = 0
new_h.steric_number = 0
new_h.number_of_lone_pairs = 0
new_h.number_of_protons_to_add = 0
new_h.num_pi_elec_2_3_bonds = 0
new_h.is_protonates = True
atom.bonded_atoms.append(new_h)
atom.number_of_protons_to_add -= 1
atom.conformation_container.add_atom(new_h)
# update names of all protons on this atom
new_H.residue_label = "%-3s%4d%2s" % (new_H.name,new_H.res_num, new_H.chain_id)
new_h.residue_label = "%-3s%4d%2s" % (new_h.name, new_h.res_num,
new_h.chain_id)
no_protons = atom.count_bonded_elements('H')
if no_protons > 1:
i = 1
for proton in atom.get_bonded_elements('H'):
proton.name = 'H%s%d'%(atom.name[1:],i)
proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.res_num, proton.chain_id)
i+=1
proton.name = 'H%s%d' % (atom.name[1:], i)
proton.residue_label = "%-3s%4d%2s" % (proton.name,
proton.res_num,
proton.chain_id)
i += 1
debug('added', new_h, 'to', atom)
def set_bond_distance(self, bvec, element):
"""Set bond distance between atom and element.
debug('added',new_H, 'to',atom)
return
def set_bond_distance(self, a, element):
d = 1.0
Args:
bvec: bond vector
element: bonded element
Returns:
scaled bond vector
"""
dist = 1.0
if element in list(self.bond_lengths.keys()):
d = self.bond_lengths[element]
dist = self.bond_lengths[element]
else:
warning('Bond length for %s not found, using the standard value of %f' % (element, d))
a = a.rescale(d)
return a
if __name__ == '__main__':
import protein, pdb, sys,os
arguments = sys.argv
if len(arguments) != 2:
info('Usage: protonate.py <pdb_file>')
sys.exit(0)
filename = arguments[1]
if not os.path.isfile(filename):
info('Error: Could not find \"%s\"' % filename)
sys.exit(1)
p = Protonate()
pdblist = pdb.readPDB(filename)
my_protein = protein.Protein(pdblist,'test.pdb')
p.remove_all_hydrogen_atoms_from_protein(my_protein)
my_protein.write_pdb('before_protonation.pdb')
p.protonate_protein(my_protein)
## write out protonated file
my_protein.write_pdb('protonated.pdb')
str_ = ('Bond length for %s not found, using the standard value '
'of %f' % (element, dist))
warning(str_)
bvec = bvec.rescale(dist)
return bvec