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

Merge pull request #33 from Electrostatics/nathan/delint

Browse Source 
De-lint atoms.py and part of bonds.py
This commit is contained in:
Nathan Baker
2020-05-23 09:39:02 -07:00
committed by GitHub
parent 62f7dfce1a 120718fbd8
commit 3cfcc8f6d3
11 changed files with 615 additions and 699 deletions
Download Patch File Download Diff File Expand all files Collapse all files

525
propka/atom.py
View File

@@ -1,66 +1,74 @@
"""Atom class - contains all atom information found in the PDB file"""
from __future__ import division import string
from __future__ import print_function import propka.lib
import propka.group
import string, propka.lib, propka.group
from . import hybrid36 from . import hybrid36
class Atom:
""" class Atom(object):
Atom class - contains all atom information found in the pdbfile """Atom class - contains all atom information found in the PDB file"""
"""
def __init__(self, line=None, verbose=False): def __init__(self, line=None, verbose=False):
"""Initialize Atom object.
self.set_properties(line) Args:
line: Line from a PDB file to set properties of atom.
self.residue_label = "%-3s%4d%2s" % (self.name,self.resNumb, self.chainID) verbose: TODO - this does not appear to be used. Can we remove it?
"""
self.groups_extracted = 0 self.occ = None
self.numb = None
self.res_name = None
self.type = None
self.chain_id = None
self.beta = None
self.icode = None
self.res_num = None
self.name = None
self.element = None
self.x = None
self.y = None
self.z = None
self.group = None self.group = None
self.group_type = None self.group_type = None
self.number_of_bonded_elements = {} self.number_of_bonded_elements = {}
self.cysteine_bridge = False self.cysteine_bridge = False
self.bonded_atoms = [] self.bonded_atoms = []
self.residue = None self.residue = None
self.conformation_container = None self.conformation_container = None
self.molecular_container = None self.molecular_container = None
self.is_protonated = False self.is_protonated = False
self.steric_number_and_lone_pairs_set = False self.steric_num_lone_pairs_set = False
self.terminal = None self.terminal = None
self.charge = 0 self.charge = 0
self.charge_set = False self.charge_set = False
self.steric_number = 0 self.steric_number = 0
self.number_of_lone_pairs = 0 self.number_of_lone_pairs = 0
self.number_of_protons_to_add = 0 self.number_of_protons_to_add = 0
self.number_of_pi_electrons_in_double_and_triple_bonds = 0 self.num_pi_elec_2_3_bonds = 0
self.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = 0 self.num_pi_elec_conj_2_3_bonds = 0
self.groups_extracted = 0
self.set_properties(line)
self.residue_label = "%-3s%4d%2s" % (self.name, self.res_num, self.chain_id)
# ligand atom types # ligand atom types
self.sybyl_type = '' self.sybyl_type = ''
self.sybyl_assigned = False self.sybyl_assigned = False
self.marvin_pka = False self.marvin_pka = False
return
def set_properties(self, line): def set_properties(self, line):
"""Line from PDB file to set properties of atom.
Args:
line: PDB file line
"""
self.name = '' self.name = ''
self.numb = 0 self.numb = 0
self.x = 0.0 self.x = 0.0
self.y = 0.0 self.y = 0.0
self.z = 0.0 self.z = 0.0
self.resNumb = 0 self.res_num = 0
self.resName = '' self.res_name = ''
self.chainID = 'A' self.chain_id = 'A'
self.type = '' self.type = ''
self.occ = '1.0' self.occ = '1.0'
self.beta = '0.0' self.beta = '0.0'
@@ -69,19 +77,19 @@ class Atom:
if line: if line:
self.name = line[12:16].strip() self.name = line[12:16].strip()
self.numb = int( hybrid36.decode(line[ 6:11]) ) self.numb = int(hybrid36.decode(line[6:11]))
self.x = float( line[30:38].strip() ) self.x = float(line[30:38].strip())
self.y = float( line[38:46].strip() ) self.y = float(line[38:46].strip())
self.z = float( line[46:54].strip() ) self.z = float(line[46:54].strip())
self.resNumb = int( line[22:26].strip() ) self.res_num = int(line[22:26].strip())
self.resName = "%-3s" % (line[17:20].strip()) self.res_name = "%-3s" % (line[17:20].strip())
self.chainID = line[21] self.chain_id = line[21]
# Set chain id to "_" if it is just white space. # Set chain id to "_" if it is just white space.
if not self.chainID.strip(): if not self.chain_id.strip():
self.chainID = '_' self.chain_id = '_'
self.type = line[:6].strip().lower() self.type = line[:6].strip().lower()
if self.resName in ['DA ','DC ','DG ','DT ']: if self.res_name in ['DA ', 'DC ', 'DG ', 'DT ']:
self.type = 'hetatm' self.type = 'hetatm'
self.occ = line[55:60].strip() self.occ = line[55:60].strip()
@@ -92,104 +100,136 @@ class Atom:
self.element = line[12:14].strip().strip(string.digits) self.element = line[12:14].strip().strip(string.digits)
if len(self.name) == 4: if len(self.name) == 4:
self.element = self.element[0] self.element = self.element[0]
if len(self.element)==2: if len(self.element) == 2:
self.element = '%1s%1s'%(self.element[0], self.element[1].lower()) self.element = '%1s%1s' % (self.element[0], self.element[1].lower())
return def set_group_type(self, type_):
"""Set group type of atom.
def set_group_type(self, type): Args:
self.group_type = type type_: group type of atom
return """
self.group_type = type_
def count_bonded_elements(self, element): def count_bonded_elements(self, element):
res = 0 """Count number of bonded atoms with same element.
for ba in self.bonded_atoms:
if element == ba.element:
res +=1
return res
Args:
element: element type for test.
Returns:
number of bonded atoms.
"""
return len(self.get_bonded_elements(element))
def get_bonded_elements(self, element): def get_bonded_elements(self, element):
"""Get bonded atoms with same element.
Args:
element: element type for test.
Returns:
array of bonded atoms.
"""
res = [] res = []
for ba in self.bonded_atoms: for ba in self.bonded_atoms:
if ba.element == element: if ba.element == element:
res.append(ba) res.append(ba)
return res return res
def get_bonded_heavy_atoms(self): def get_bonded_heavy_atoms(self):
return [ba for ba in self.bonded_atoms if ba.element!='H'] """Get the atoms bonded to this one that aren't hydrogen.
Returns:
list of atoms.
"""
return [ba for ba in self.bonded_atoms if ba.element != 'H']
def is_atom_within_bond_distance(self, other_atom, max_bonds, cur_bond): def is_atom_within_bond_distance(self, other_atom, max_bonds, cur_bond):
""" check if <other_atom> is found within <max_bonds> bonds of self """ """Check if <other_atom> is found within <max_bonds> bonds of self.
Args:
other_atom: atom to check
max_bonds: number of bonds to check for other atom bonding to self
Returns:
Boolean for atom bond distance
"""
for ba in self.bonded_atoms: for ba in self.bonded_atoms:
if ba == other_atom: if ba == other_atom:
return True return True
if max_bonds > cur_bond: if max_bonds > cur_bond:
if ba.is_atom_within_bond_distance(other_atom, max_bonds, cur_bond+1): if ba.is_atom_within_bond_distance(other_atom, max_bonds, cur_bond+1):
return True return True
return False return False
def set_property(self, numb=None, name=None, res_name=None, chain_id=None,
res_num=None, x=None, y=None, z=None, occ=None, beta=None):
"""Set properties of the atom object.
Args:
def setProperty(self, numb: Atom number
numb = None, name: Atom name
name = None, res_name: residue name
resName = None, chain_id: chain ID
chainID = None, res_num: residue number
resNumb = None, x: atom x-coordinate
x = None, y: atom y-coordinate
y = None, z: atom z-coordinate
z = None, occ: atom occupancy
occ = None, beta: atom temperature factor
beta = None):
""" """
sets properties of the atom object if numb is not None:
self.numb = numb
if name is not None:
self.name = name
if res_name is not None:
self.res_name = res_name
if chain_id is not None:
self.chain_id = chain_id
if res_num is not None:
self.res_num = res_num
if x is not None:
self.x = x
if y is not None:
self.y = y
if z is not None:
self.z = z
if occ is not None:
self.occ = occ
if beta is not None:
self.beta = beta
def make_copy(self):
"""Make a copy of this atom.
Returns:
Another atom object copy of this one.
""" """
new_atom = Atom()
if numb != None: self.numb = numb new_atom.type = self.type
if name != None: self.name = name new_atom.numb = self.numb
if resName != None: self.resName = resName new_atom.name = self.name
if chainID != None: self.chainID = chainID new_atom.element = self.element
if resNumb != None: self.resNumb = resNumb new_atom.res_name = self.res_name
if x != None: self.x = x new_atom.res_num = self.res_num
if y != None: self.y = y new_atom.chain_id = self.chain_id
if z != None: self.z = z new_atom.x = self.x
if occ != None: self.occ = occ new_atom.y = self.y
if beta != None: self.beta = beta new_atom.z = self.z
new_atom.occ = self.occ
new_atom.beta = self.beta
new_atom.terminal = self.terminal
def makeCopy(self): new_atom.residue_label = self.residue_label
""" new_atom.icode = self.icode
making a copy of this atom return new_atom
"""
newAtom = Atom()
newAtom.type = self.type
newAtom.numb = self.numb
newAtom.name = self.name
newAtom.element = self.element
newAtom.resName = self.resName
newAtom.resNumb = self.resNumb
newAtom.chainID = self.chainID
newAtom.x = self.x
newAtom.y = self.y
newAtom.z = self.z
newAtom.occ = self.occ
newAtom.beta = self.beta
newAtom.terminal = self.terminal
newAtom.residue_label = self.residue_label
newAtom.icode = self.icode
return newAtom
def make_input_line(self): def make_input_line(self):
# making string """PDB line for this atom.
TODO - Could be @property method/attribute
TODO - figure out difference between make_pdb_line, make_input_line, and make_pdb_line2
Returns:
String with PDB-format line.
"""
group = '-' group = '-'
model_pka = '-' model_pka = '-'
if self.group: if self.group:
@@ -199,197 +239,168 @@ class Atom:
if self.group.titratable: if self.group.titratable:
model_pka = '%6.2f'%self.group.model_pka model_pka = '%6.2f'%self.group.model_pka
str_ = "%-6s%5d %s " % (self.type.upper(), self.numb,
propka.lib.makeTidyAtomLabel(self.name, self.element))
str = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s \n" % (self.type.upper(), str_ += "%s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s \n" % (self.res_name, self.chain_id,
self.numb, self.res_num, self.x, self.y,
propka.lib.makeTidyAtomLabel(self.name, self.z, group, model_pka)
self.element), return str_
self.resName,
self.chainID,
self.resNumb,
self.x,
self.y,
self.z,
group,
model_pka)
return str
def make_conect_line(self): def make_conect_line(self):
res = 'CONECT%5d'%self.numb """PDB line for bonding within this molecule.
Returns:
String with PDB line.
"""
res = 'CONECT%5d' % self.numb
# extract and sort numbers of bonded residues
bonded = [] bonded = []
for atom in self.bonded_atoms: for atom in self.bonded_atoms:
bonded.append(atom.numb) bonded.append(atom.numb)
bonded.sort() bonded.sort()
# write 'em out
for b in bonded: for b in bonded:
res += '%5d'%b res += '%5d'%b
res += '\n' res += '\n'
return res return res
def get_input_parameters(self): def get_input_parameters(self):
""" Method for getting the input parameters stored in the """Extract the input parameters stored in the occupancy and b-factor
occupancy and b-factor fields in input files""" fields in input files"""
# Set the group type # Set the group type
if self.occ != '-': if self.occ != '-':
# make sure to set the terminal # make sure to set the terminal
if self.occ in ['N+','C-']: if self.occ in ['N+', 'C-']:
self.terminal = self.occ self.terminal = self.occ
# save the ligand group charge # save the ligand group charge
if self.occ =='BLG': if self.occ == 'BLG':
self.charge=+1 self.charge = +1
elif self.occ =='ALG': elif self.occ == 'ALG':
self.charge=-1 self.charge = -1
# generic ions # generic ions
if self.occ in ['1P','2P','1N','2N']: if self.occ in ['1P', '2P', '1N', '2N']:
self.resName=self.occ self.res_name = self.occ
self.occ='Ion' self.occ = 'Ion'
# correct the group type # correct the group type
self.occ = self.occ.replace('N+','Nterm') self.occ = self.occ.replace('N+', 'Nterm')
self.occ = self.occ.replace('C-','Cterm') self.occ = self.occ.replace('C-', 'Cterm')
self.occ = self.occ.replace('ION','Ion') self.occ = self.occ.replace('ION', 'Ion')
self.occ = self.occ.replace('ALG','titratable_ligand') self.occ = self.occ.replace('ALG', 'titratable_ligand')
self.occ = self.occ.replace('BLG','titratable_ligand') self.occ = self.occ.replace('BLG', 'titratable_ligand')
self.occ = self.occ.replace('LG','non_titratable_ligand') self.occ = self.occ.replace('LG', 'non_titratable_ligand')
# try to initialise the group # try to initialise the group
try: try:
exec('self.group = propka.group.%s_group(self)'%self.occ) # TODO - get rid of this exec() statement for security reasons
exec('self.group = propka.group.%s_group(self)' % self.occ)
except: except:
raise Exception('%s in input_file is not recognized as a group'%self.occ) raise Exception('%s in input_file is not recognized as a group' % self.occ)
# set the model pKa value # set the model pKa value
if self.beta != '-': if self.beta != '-':
self.group.model_pka = float(self.beta) self.group.model_pka = float(self.beta)
self.group.model_pka_set = True self.group.model_pka_set = True
# set occ and beta to standard values # set occ and beta to standard values
self.occ = '1.00' self.occ = '1.00'
self.beta = '0.00' self.beta = '0.00'
return
def make_pdb_line(self): def make_pdb_line(self):
"""Create PDB line.
# making string TODO - this could/should be a @property method/attribute
str = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s\n" % (self.type.upper(), TODO - figure out difference between make_pdb_line, make_input_line, and make_pdb_line2
self.numb,
propka.lib.makeTidyAtomLabel(self.name,
self.element),
self.resName,
self.chainID,
self.resNumb,
self.x,
self.y,
self.z,
self.occ,
self.beta)
Returns:
return str String with PDB line.
def make_mol2_line(self,id):
#1 S1 3.6147 2.0531 1.4795 S.3 1 noname -0.1785
# making string
str = "%-4d %-4s %10.4f %10.4f %10.4f %6s %6d %10s %10.4f\n" % (id,
propka.lib.makeTidyAtomLabel(self.name,
self.element),
self.x,
self.y,
self.z,
self.sybyl_type.replace('-',''),
self.resNumb,
self.resName,
0.0)#self.charge)
return str
def makePDBLine(self,
numb = None,
name = None,
resName = None,
chainID = None,
resNumb = None,
x = None,
y = None,
z = None,
occ = None,
beta = None):
""" """
returns a pdb ATOM-line for various purposes; str_ = "%-6s%5d " % (self.type.upper(), self.numb)
specifying arguments over-writes. str_ += "%s %s" % (propka.lib.makeTidyAtomLabel(self.name, self.element),
self.res_name)
str_ += "%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s\n" % (self.chain_id, self.res_num,
self.x, self.y, self.z,
self.occ, self.beta)
return str_
def make_mol2_line(self, id_):
"""Create MOL2 line.
Format: 1 S1 3.6147 2.0531 1.4795 S.3 1 noname -0.1785
TODO - this could/should be a @property method/attribute
Returns:
String with MOL2 line.
""" """
if numb == None: numb = self.numb str_ = "%-4d %-4s " % (id_, propka.lib.makeTidyAtomLabel(self.name,
if name == None: name = self.name self.element))
if resName == None: resName = self.resName str_ += "%10.4f %10.4f %10.4f " % (self.x, self.y, self.z)
if chainID == None: chainID = self.chainID str_ += "%6s %6d %10s %10.4f\n" % (self.sybyl_type.replace('-', ''),
if resNumb == None: resNumb = self.resNumb self.res_num, self.res_name, 0.0)
if x == None: x = self.x return str_
if y == None: y = self.y
if z == None: z = self.z
if occ == None: occ = self.occ
if beta == None: beta = self.beta
# making string def make_pdb_line2(self, numb=None, name=None, res_name=None, chain_id=None,
str = "ATOM " res_num=None, x=None, y=None, z=None, occ=None, beta=None):
str += "%6d" % (numb) """Create a PDB line.
str += " %s" % (propka.lib.makeTidyAtomLabel(name,self.element))
str += " %s" % (resName)
str += "%2s" % (chainID)
str += "%4d" % (resNumb)
str += "%12.3lf" % (x)
str += "%8.3lf" % (y)
str += "%8.3lf" % (z)
str += "%6.2lf" % (occ)
str += "%6.2lf" % (beta)
str += '\n'
return str
TODO - this could/should be a @property method/attribute
TODO - figure out difference between make_pdb_line, make_input_line, and make_pdb_line2
def getTidyLabel(self): Returns:
String with PDB line.
""" """
Returns a 'tidier' atom label for printing the new pdbfile if numb is None:
""" numb = self.numb
return propka.lib.makeTidyAtomLabel(self.name,self.element) if name is None:
name = self.name
if res_name is None:
res_name = self.res_name
if chain_id is None:
chain_id = self.chain_id
if res_num is None:
res_num = self.res_num
if x is None:
x = self.x
if y is None:
y = self.y
if z is None:
z = self.z
if occ is None:
occ = self.occ
if beta is None:
beta = self.beta
str_ = "ATOM "
str_ += "%6d" % (numb)
str_ += " %s" % (propka.lib.makeTidyAtomLabel(name, self.element))
str_ += " %s" % (res_name)
str_ += "%2s" % (chain_id)
str_ += "%4d" % (res_num)
str_ += "%12.3lf" % (x)
str_ += "%8.3lf" % (y)
str_ += "%8.3lf" % (z)
str_ += "%6.2lf" % (occ)
str_ += "%6.2lf" % (beta)
str_ += '\n'
return str_
def get_tidy_label(self):
"""Returns a 'tidier' atom label for printing the new pdbfile
TODO - this could/should be a @property method/attribute
Returns:
String with label"""
return propka.lib.makeTidyAtomLabel(self.name, self.element)
def __str__(self): def __str__(self):
return '%5d-%4s %5d-%3s (%1s) [%8.3f %8.3f %8.3f] %s' %(self.numb, self.name, self.resNumb, self.resName, self.chainID, self.x, self.y, self.z,self.element) """Return an undefined-format string version of this atom."""
return '%5d-%4s %5d-%3s (%1s) [%8.3f %8.3f %8.3f] %s' % (self.numb, self.name,
self.res_num, self.res_name,
# def get_element(self): self.chain_id, self.x, self.y,
# """ try to extract element if not already done""" self.z, self.element)
# if self.element == '':
# self.element = self.name.strip(string.digits)[0]
# return self.element
def set_residue(self, residue): def set_residue(self, residue):
""" Makes a references to the parent residue""" """ Makes a reference to the parent residue
if self.residue == None:
Args:
residue: the parent residue
"""
if self.residue is None:
self.residue = residue self.residue = residue

548
propka/bonds.py
View File

@@ -1,336 +1,334 @@
"""PROPKA representation of bonds."""
from __future__ import division # TODO - is pickle still used?
from __future__ import print_function import pickle
# TODO - eliminate use of sys
import pickle,sys,os,math,propka.calculations import sys
# TODO - eliminate use of os
import os
import math
import json import json
import pkg_resources import pkg_resources
import propka.calculations
# TODO - replace the info/warning imports with logging functionality
from propka.lib import info, warning from propka.lib import info, warning
class bondmaker:
# TODO - should these constants be defined higher up in the module?
# TODO - I don't know what some of these constants mean
DISULFIDE_DISTANCE = 2.5
FLUORIDE_DISTANCE = 1.7
HYDROGEN_DISTANCE = 1.5
DEFAULT_DISTANCE = 2.0
BOX_SIZE = 2.5
POS_MAX = 1e6
class BondMaker:
"""Makes bonds?
TODO - the documentation for this class needs to be improved.
"""
def __init__(self): def __init__(self):
# predefined bonding distances # predefined bonding distances
self.distances = { self.distances = {'S-S' : DISULFIDE_DISTANCE, 'F-F' : FLUORIDE_DISTANCE}
'S-S':2.5,
'F-F':1.7}
self.distances_squared = {} self.distances_squared = {}
for k in self.distances.keys(): for k in self.distances.keys():
self.distances_squared[k]=self.distances[k]*self.distances[k] self.distances_squared[k] = self.distances[k] * self.distances[k]
self.H_dist = HYDROGEN_DISTANCE
self.H_dist = 1.5; self.default_dist = DEFAULT_DISTANCE
self.default_dist = 2.0; self.H_dist_squared = self.H_dist * self.H_dist
self.H_dist_squared = self.H_dist * self.H_dist
self.default_dist_squared = self.default_dist * self.default_dist self.default_dist_squared = self.default_dist * self.default_dist
distances = list(self.distances_squared.values()) + [self.default_dist_squared]
self.max_sq_distance = max(list(self.distances_squared.values())+[self.default_dist_squared]) self.max_sq_distance = max(distances)
# protein bonding data # protein bonding data
self.data_file_name = pkg_resources.resource_filename(__name__, 'protein_bonds.json') self.data_file_name = pkg_resources.resource_filename(__name__, 'protein_bonds.json')
with open(self.data_file_name,'rt') as json_file: with open(self.data_file_name, 'rt') as json_file:
self.protein_bonds = json.load(json_file) self.protein_bonds = json.load(json_file)
self.intra_residue_backbone_bonds = {'N': ['CA'], 'CA': ['N', 'C'],
'C': ['CA', 'O'], 'O': ['C']}
self.intra_residue_backbone_bonds = {'N': ['CA'], self.num_pi_elec_bonds_backbone = {'C': 1, 'O': 1}
'CA':['N','C'], self.num_pi_elec_conj_bonds_backbone = {'N': 1}
'C': ['CA','O'], self.num_pi_elec_bonds_sidechains = {'ARG-CZ' : 1, 'ARG-NH1': 1,
'O': ['C']} 'ASN-OD1': 1, 'ASN-CG' : 1,
'ASP-OD1': 1, 'ASP-CG' : 1,
self.number_of_pi_electrons_in_bonds_in_backbone = {'C':1, 'GLU-OE1': 1, 'GLU-CD' : 1,
'O':1} 'GLN-OE1': 1, 'GLN-CD' : 1,
'HIS-CG' : 1, 'HIS-CD2': 1,
self.number_of_pi_electrons_in_conjugate_bonds_in_backbone = {'N':1} 'HIS-ND1': 1, 'HIS-CE1': 1,
'PHE-CG' : 1, 'PHE-CD1': 1,
self.number_of_pi_electrons_in_bonds_in_sidechains = {'ARG-CZ' :1, 'PHE-CE1': 1, 'PHE-CZ' : 1,
'ARG-NH1':1, 'PHE-CE2': 1, 'PHE-CD2': 1,
'ASN-OD1':1, 'TRP-CG' : 1, 'TRP-CD1': 1,
'ASN-CG' :1, 'TRP-CE2': 1, 'TRP-CD2': 1,
'ASP-OD1':1, 'TRP-CE3': 1, 'TRP-CZ3': 1,
'ASP-CG' :1, 'TRP-CH2': 1, 'TRP-CZ2': 1,
'GLU-OE1':1, 'TYR-CG' : 1, 'TYR-CD1': 1,
'GLU-CD' :1, 'TYR-CE1': 1, 'TYR-CZ' : 1,
'GLN-OE1':1, 'TYR-CE2': 1, 'TYR-CD2': 1}
'GLN-CD' :1, self.num_pi_elec_conj_bonds_sidechains = {'ARG-NE': 1, 'ARG-NH2': 1,
'HIS-CG' :1, 'ASN-ND2': 1, 'GLN-NE2': 1,
'HIS-CD2':1, 'HIS-NE2': 1, 'TRP-NE1': 1}
'HIS-ND1':1, self.num_pi_elec_bonds_ligands = {'C.ar': 1, 'N.pl3': 0, 'C.2': 1,
'HIS-CE1':1, 'O.2': 1, 'O.co2': 1, 'N.ar': 1,
'PHE-CG' :1, 'C.1': 2, 'N.1': 2}
'PHE-CD1':1, self.num_pi_elec_conj_bonds_ligands = {'N.am': 1, 'N.pl3': 1}
'PHE-CE1':1,
'PHE-CZ' :1,
'PHE-CE2':1,
'PHE-CD2':1,
'TRP-CG' :1,
'TRP-CD1':1,
'TRP-CE2':1,
'TRP-CD2':1,
'TRP-CE3':1,
'TRP-CZ3':1,
'TRP-CH2':1,
'TRP-CZ2':1,
'TYR-CG' :1,
'TYR-CD1':1,
'TYR-CE1':1,
'TYR-CZ' :1,
'TYR-CE2':1,
'TYR-CD2':1}
self.number_of_pi_electrons_in_conjugate_bonds_in_sidechains = {'ARG-NE' :1,
'ARG-NH2':1,
'ASN-ND2':1,
'GLN-NE2':1,
'HIS-NE2':1,
'TRP-NE1':1}
self.number_of_pi_electrons_in_bonds_ligands = {'C.ar':1,
'N.pl3':0,
'C.2':1,
'O.2':1,
'O.co2':1,
'N.ar':1,
'C.1':2,
'N.1':2}
self.number_of_pi_electrons_in_conjugate_bonds_in_ligands = {'N.am':1,'N.pl3':1}
self.backbone_atoms = list(self.intra_residue_backbone_bonds.keys()) self.backbone_atoms = list(self.intra_residue_backbone_bonds.keys())
self.terminal_oxygen_names = ['OXT', 'O\'\'']
self.terminal_oxygen_names = ['OXT','O\'\'']
return
def find_bonds_for_protein(self, protein): def find_bonds_for_protein(self, protein):
""" Finds bonds proteins based on the way atoms """Finds bonds proteins based on the way atoms normally bond in proteins.
normally bond in proteins"""
Args:
protein: the protein to search for bonds
"""
info('++++ Side chains ++++') info('++++ Side chains ++++')
# side chains # side chains
for chain in protein.chains: for chain in protein.chains:
for residue in chain.residues: for residue in chain.residues:
if residue.resName.replace(' ','') not in ['N+','C-']: if residue.res_name.replace(' ', '') not in ['N+', 'C-']:
self.find_bonds_for_side_chain(residue.atoms) self.find_bonds_for_side_chain(residue.atoms)
info('++++ Backbones ++++') info('++++ Backbones ++++')
# backbone # backbone
last_residues = [] last_residues = []
for chain in protein.chains: for chain in protein.chains:
for i in range(1,len(chain.residues)): for i in range(1, len(chain.residues)):
if chain.residues[i-1].resName.replace(' ','') not in ['N+','C-']: if chain.residues[i-1].res_name.replace(' ', '') not in ['N+', 'C-']:
if chain.residues[i].resName.replace(' ','') not in ['N+','C-']: if chain.residues[i].res_name.replace(' ', '') not in ['N+', 'C-']:
self.connect_backbone(chain.residues[i-1], chain.residues[i]) self.connect_backbone(chain.residues[i-1], chain.residues[i])
last_residues.append(chain.residues[i]) last_residues.append(chain.residues[i])
info('++++ terminal oxygen ++++') info('++++ terminal oxygen ++++')
# terminal OXT # terminal OXT
for last_residue in last_residues: for last_residue in last_residues:
self.find_bonds_for_terminal_oxygen(last_residue) self.find_bonds_for_terminal_oxygen(last_residue)
info('++++ cysteines ++++') info('++++ cysteines ++++')
# Cysteines # Cysteines
for chain in protein.chains: for chain in protein.chains:
for i in range(0,len(chain.residues)): for i in range(0, len(chain.residues)):
if chain.residues[i].resName == 'CYS': if chain.residues[i].res_name == 'CYS':
for j in range(0,len(chain.residues)): for j in range(0, len(chain.residues)):
if chain.residues[j].resName == 'CYS' and j != i: if chain.residues[j].res_name == 'CYS' and j != i:
self.check_for_cysteine_bonds(chain.residues[i], self.check_for_cysteine_bonds(chain.residues[i],
chain.residues[j]) chain.residues[j])
return
def check_for_cysteine_bonds(self, cys1, cys2): def check_for_cysteine_bonds(self, cys1, cys2):
"""Looks for potential bonds between two cysteines.
Args:
cys1: one of the cysteines to check
cys1: one of the cysteines to check
"""
for atom1 in cys1.atoms: for atom1 in cys1.atoms:
if atom1.name == 'SG': if atom1.name == 'SG':
for atom2 in cys2.atoms: for atom2 in cys2.atoms:
if atom2.name == 'SG': if atom2.name == 'SG':
if propka.calculations.squared_distance(atom1,atom2) < self.SS_dist_squared: dist = propka.calculations.squared_distance(atom1, atom2)
# TODO - is SS_dist_squared an attribute of this class?
if dist < self.SS_dist_squared:
self.make_bond(atom1, atom2) self.make_bond(atom1, atom2)
return
def find_bonds_for_terminal_oxygen(self, residue): def find_bonds_for_terminal_oxygen(self, residue):
"""Look for bonds for terminal oxygen.
Args:
residue - test residue
"""
for atom1 in residue.atoms: for atom1 in residue.atoms:
if atom1.name in self.terminal_oxygen_names: if atom1.name in self.terminal_oxygen_names:
for atom2 in residue.atoms: for atom2 in residue.atoms:
if atom2.name == 'C': if atom2.name == 'C':
self.make_bond(atom1, atom2) self.make_bond(atom1, atom2)
return # TODO - stopped here.
def connect_backbone(self, residue1, residue2): def connect_backbone(self, residue1, residue2):
""" Sets up bonds in the backbone """ """Sets up bonds in the backbone
# residue 1
Args:
residue1: first residue to connect
residue2: second residue to connect
"""
self.find_bonds_for_residue_backbone(residue1) self.find_bonds_for_residue_backbone(residue1)
# residue 2
self.find_bonds_for_residue_backbone(residue2) self.find_bonds_for_residue_backbone(residue2)
# inter-residue bond
for atom1 in residue1.atoms: for atom1 in residue1.atoms:
if atom1.name == 'C': if atom1.name == 'C':
for atom2 in residue2.atoms: for atom2 in residue2.atoms:
if atom2.name == 'N': if atom2.name == 'N':
if propka.calculations.squared_distance(atom1,atom2) < self.default_dist_squared: if propka.calculations.squared_distance(atom1, atom2) \
< self.default_dist_squared:
self.make_bond(atom1, atom2) self.make_bond(atom1, atom2)
return
def find_bonds_for_residue_backbone(self, residue): def find_bonds_for_residue_backbone(self, residue):
"""Find bonds for this residue's backbone.
Args:
residue: reside to search for backbone bonds.
"""
for atom1 in residue.atoms: for atom1 in residue.atoms:
if atom1.name in list(self.number_of_pi_electrons_in_bonds_in_backbone.keys()): if atom1.name in list(self.num_pi_elec_bonds_backbone.keys()):
atom1.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_in_backbone[atom1.name] atom1.num_pi_elec_2_3_bonds \
if atom1.name in list(self.number_of_pi_electrons_in_conjugate_bonds_in_backbone.keys()) and len(atom1.bonded_atoms)>1: # last part to avoid including N-term = self.num_pi_elec_bonds_backbone[atom1.name]
atom1.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_backbone[atom1.name] if atom1.name in \
list(self.num_pi_elec_conj_bonds_backbone.keys()) \
and len(atom1.bonded_atoms) > 1: # last part to avoid including N-term
atom1.num_pi_elec_conj_2_3_bonds \
= self.num_pi_elec_conj_bonds_backbone[atom1.name]
if atom1.name in self.backbone_atoms: if atom1.name in self.backbone_atoms:
for atom2 in residue.atoms: for atom2 in residue.atoms:
if atom2.name in self.intra_residue_backbone_bonds[atom1.name]: if atom2.name in self.intra_residue_backbone_bonds[atom1.name]:
self.make_bond(atom1, atom2) self.make_bond(atom1, atom2)
return
def find_bonds_for_side_chain(self, atoms): def find_bonds_for_side_chain(self, atoms):
""" Finds bonds for a side chain """ """Finds bonds for a side chain.
for atom1 in atoms:
key = '%s-%s'%(atom1.resName,atom1.name) Args:
if key in list(self.number_of_pi_electrons_in_bonds_in_sidechains.keys()): atoms: list of atoms to check for bonds
atom1.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_in_sidechains[key] """
if key in list(self.number_of_pi_electrons_in_conjugate_bonds_in_sidechains.keys()): for atom1 in atoms:
atom1.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_sidechains[key] key = '%s-%s' % (atom1.res_name, atom1.name)
if key in list(self.num_pi_elec_bonds_sidechains.keys()):
atom1.num_pi_elec_2_3_bonds \
= self.num_pi_elec_bonds_sidechains[key]
if key in list(self.num_pi_elec_conj_bonds_sidechains.keys()):
atom1.num_pi_elec_conj_2_3_bonds \
= self.num_pi_elec_conj_bonds_sidechains[key]
if not atom1.name in self.backbone_atoms: if not atom1.name in self.backbone_atoms:
if not atom1.name in self.terminal_oxygen_names: if not atom1.name in self.terminal_oxygen_names:
for atom2 in atoms: for atom2 in atoms:
if atom2.name in self.protein_bonds[atom1.resName][atom1.name]: if atom2.name in self.protein_bonds[atom1.res_name][atom1.name]:
self.make_bond(atom1,atom2) self.make_bond(atom1, atom2)
return
def find_bonds_for_ligand(self, ligand): def find_bonds_for_ligand(self, ligand):
""" Finds bonds for all atoms in the molecule """ """Finds bonds for all atoms in the ligand molecule
Args:
ligand: ligand molecule to search for bonds
"""
# identify bonding atoms # identify bonding atoms
self.find_bonds_for_atoms(ligand.atoms) self.find_bonds_for_atoms(ligand.atoms)
self.add_pi_electron_table_info(ligand.atoms) self.add_pi_electron_table_info(ligand.atoms)
return
def add_pi_electron_table_info(self, atoms): def add_pi_electron_table_info(self, atoms):
"""Add table information on pi electrons
Args:
atoms: list of atoms for pi electron table information checking
"""
# apply table information on pi-electrons # apply table information on pi-electrons
for atom in atoms: for atom in atoms:
# for ligands # for ligands
if atom.type == 'hetatm': if atom.type == 'hetatm':
if atom.sybyl_type in self.number_of_pi_electrons_in_bonds_ligands.keys(): if atom.sybyl_type in self.num_pi_elec_bonds_ligands.keys():
atom.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_ligands[atom.sybyl_type] atom.num_pi_elec_2_3_bonds = self.num_pi_elec_bonds_ligands[atom.sybyl_type]
if atom.sybyl_type in self.number_of_pi_electrons_in_conjugate_bonds_in_ligands.keys(): if atom.sybyl_type in self.num_pi_elec_conj_bonds_ligands.keys():
atom.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_ligands[atom.sybyl_type] atom.num_pi_elec_conj_2_3_bonds \
= self.num_pi_elec_conj_bonds_ligands[atom.sybyl_type]
# for protein # for protein
if atom.type == 'atom': if atom.type == 'atom':
key = '%s-%s'%(atom.resName,atom.name) key = '%s-%s' % (atom.res_name, atom.name)
if key in list(self.number_of_pi_electrons_in_bonds_in_sidechains.keys()): if key in list(self.num_pi_elec_bonds_sidechains.keys()):
atom.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_in_sidechains[key] atom.num_pi_elec_2_3_bonds = self.num_pi_elec_bonds_sidechains[key]
if key in list(self.number_of_pi_electrons_in_conjugate_bonds_in_sidechains.keys()): if key in list(self.num_pi_elec_conj_bonds_sidechains.keys()):
atom.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_sidechains[key] atom.num_pi_elec_conj_2_3_bonds = self.num_pi_elec_conj_bonds_sidechains[key]
if atom.name in list(self.number_of_pi_electrons_in_bonds_in_backbone.keys()):
atom.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_in_backbone[atom.name]
if atom.name in list(self.number_of_pi_electrons_in_conjugate_bonds_in_backbone.keys()) and len(atom.bonded_atoms)>1: # last part to avoid including N-term
atom.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_backbone[atom.name]
return
if atom.name in list(self.num_pi_elec_bonds_backbone.keys()):
atom.num_pi_elec_2_3_bonds = self.num_pi_elec_bonds_backbone[atom.name]
if atom.name in list(self.num_pi_elec_conj_bonds_backbone.keys()) \
and len(atom.bonded_atoms) > 1:
# last part to avoid including N-term
atom.num_pi_elec_conj_2_3_bonds \
= self.num_pi_elec_conj_bonds_backbone[atom.name]
def find_bonds_for_protein_by_distance(self, molecule): def find_bonds_for_protein_by_distance(self, molecule):
""" Finds bonds for all atoms in the molecule """ """Finds bonds for all atoms in the molecule.
Args:
molecule: molecule in which to find bonds.
Returns:
list of atoms
"""
#self.find_bonds_for_protein(molecule) #self.find_bonds_for_protein(molecule)
atoms = [] atoms = []
for chain in molecule.chains: for chain in molecule.chains:
for residue in chain.residues: for residue in chain.residues:
if residue.resName.replace(' ','') not in ['N+','C-']: if residue.res_name.replace(' ', '') not in ['N+', 'C-']:
for atom in residue.atoms: for atom in residue.atoms:
atoms.append(atom) atoms.append(atom)
self.find_bonds_for_atoms_using_boxes(atoms)
self.find_bonds_for_atoms_using_boxes(atoms) #####
#self.find_bonds_for_atoms(atoms) #####
return atoms return atoms
def find_bonds_for_atoms(self, atoms): def find_bonds_for_atoms(self, atoms):
""" Finds all bonds for a list of atoms""" """Finds all bonds for a list of atoms
Args:
atoms: list of atoms in which to find bonds.
"""
no_atoms = len(atoms) no_atoms = len(atoms)
for i in range(no_atoms): for i in range(no_atoms):
for j in range(i+1,no_atoms): for j in range(i+1, no_atoms):
if atoms[i] in atoms[j].bonded_atoms: if atoms[i] in atoms[j].bonded_atoms:
continue continue
if self.check_distance(atoms[i], atoms[j]): if self.check_distance(atoms[i], atoms[j]):
self.make_bond(atoms[i],atoms[j]) self.make_bond(atoms[i], atoms[j])
# di-sulphide bonds # di-sulphide bonds
if atoms[i].element == 'S' and atoms[j].element == 'S': if atoms[i].element == 'S' and atoms[j].element == 'S':
atoms[i].cysteine_bridge = True atoms[i].cysteine_bridge = True
atoms[j].cysteine_bridge = True atoms[j].cysteine_bridge = True
return
def check_distance(self, atom1, atom2): def check_distance(self, atom1, atom2):
sq_dist = propka.calculations.squared_distance(atom1, atom2) """Check distance between two atoms
Args:
atom1: first atom for distance check
atom2: second atom for distance check
Returns:
True if within distance, False otherwise
"""
sq_dist = propka.calculations.squared_distance(atom1, atom2)
if sq_dist > self.max_sq_distance: if sq_dist > self.max_sq_distance:
return False return False
key = '%s-%s' % (atom1.element, atom2.element)
key = '%s-%s'%(atom1.element,atom2.element)
h_count = key.count('H') h_count = key.count('H')
if sq_dist < self.H_dist_squared and h_count == 1:
if sq_dist < self.H_dist_squared and h_count==1:
return True return True
if sq_dist < self.default_dist_squared and h_count==0: if sq_dist < self.default_dist_squared and h_count == 0:
return True return True
if key in self.distances_squared.keys(): if key in self.distances_squared.keys():
if sq_dist < self.distances_squared[key]: if sq_dist < self.distances_squared[key]:
return True return True
return False return False
def find_bonds_for_molecules_using_boxes(self, molecules): def find_bonds_for_molecules_using_boxes(self, molecules):
""" Finds all bonds for a molecular container""" """ Finds all bonds for a molecular container.
Args:
molecules: list of molecules for finding bonds.
"""
for name in molecules.conformation_names: for name in molecules.conformation_names:
self.find_bonds_for_atoms_using_boxes(molecules.conformations[name].atoms) self.find_bonds_for_atoms_using_boxes(molecules.conformations[name].atoms)
#for chain in molecules.conformations[name].chains:
# self.find_bonds_for_atoms_using_boxes(molecules.conformations[name].get_chain(chain))
return
def add_pi_electron_information(self, molecules): def add_pi_electron_information(self, molecules):
"""Add pi electron information to a molecule.
Args:
molecules: list of molecules for adding pi electron information.
"""
for name in molecules.conformation_names: for name in molecules.conformation_names:
self.add_pi_electron_table_info(molecules.conformations[name].atoms) self.add_pi_electron_table_info(molecules.conformations[name].atoms)
return
def find_bonds_for_atoms_using_boxes(self, atoms): def find_bonds_for_atoms_using_boxes(self, atoms):
""" Finds all bonds for a list of atoms""" """Finds all bonds for a list of atoms.
box_size = 2.5 Args:
atoms: list of atoms for finding bonds
# find min and max coordinates """
xmin = 1e6; xmax = -1e6 box_size = BOX_SIZE
ymin = 1e6; ymax = -1e6 xmin = POS_MAX
zmin = 1e6; zmax = -1e6 xmax = -1.0 * POS_MAX
ymin = POS_MAX
ymax = -1.0 * POS_MAX
zmin = POS_MAX
zmax = -1.0 * POS_MAX
for atom in atoms: for atom in atoms:
if atom.x > xmax: if atom.x > xmax:
xmax = atom.x xmax = atom.x
@@ -344,147 +342,95 @@ class bondmaker:
ymin = atom.y ymin = atom.y
if atom.z < zmin: if atom.z < zmin:
zmin = atom.z zmin = atom.z
xlen = xmax - xmin
xlen = xmax-xmin ylen = ymax - ymin
ylen = ymax-ymin zlen = zmax - zmin
zlen = zmax-zmin self.num_box_x = max(1, int(math.ceil(xlen/box_size)))
self.num_box_y = max(1, int(math.ceil(ylen/box_size)))
#info('x range: [%6.2f;%6.2f] %6.2f'%(xmin,xmax,xlen)) self.num_box_z = max(1, int(math.ceil(zlen/box_size)))
#info('y range: [%6.2f;%6.2f] %6.2f'%(ymin,ymax,ylen))
#info('z range: [%6.2f;%6.2f] %6.2f'%(zmin,zmax,zlen))
# how many boxes do we need in each dimension?
# NOTE: math.ceil() returns an int in python3 and a float in
# python2, so we need to convert it to an int for range() to work in
# both versions. See PEP 3141.
self.no_box_x = max(1, int(math.ceil(xlen/box_size)))
self.no_box_y = max(1, int(math.ceil(ylen/box_size)))
self.no_box_z = max(1, int(math.ceil(zlen/box_size)))
#info('No. box x: %6.2f'%self.no_box_x)
#info('No. box y: %6.2f'%self.no_box_y)
#info('No. box z: %6.2f'%self.no_box_z)
# initialize boxes
self.boxes = {} self.boxes = {}
for x in range(self.no_box_x): for x in range(self.num_box_x):
for y in range(self.no_box_y): for y in range(self.num_box_y):
for z in range(self.no_box_z): for z in range(self.num_box_z):
self.boxes[(x,y,z)] = [] self.boxes[(x, y, z)] = []
# put atoms into boxes
for atom in atoms: for atom in atoms:
x = math.floor((atom.x-xmin)/box_size) x = math.floor((atom.x - xmin)/box_size)
y = math.floor((atom.y-ymin)/box_size) y = math.floor((atom.y - ymin)/box_size)
z = math.floor((atom.z-zmin)/box_size) z = math.floor((atom.z - zmin)/box_size)
self.put_atom_in_box(x,y,z,atom) self.put_atom_in_box(x, y, z, atom)
for value in self.boxes.values():
# assign bonds
for key, value in self.boxes.items():
self.find_bonds_for_atoms(value) self.find_bonds_for_atoms(value)
def put_atom_in_box(self, x, y, z, atom):
"""Put an atom in a box.
return Args:
x: box x-coordinates
def put_atom_in_box(self,x,y,z,atom): y: box y-coordinates
# atom in the x,y,z box and the up to 7 neighboring boxes on z: box z-coordinates
# one side of the x,y,z box in each dimension atom: the atom to place in a box
"""
for bx in [x,x+1]: for bx in [x, x+1]:
for by in [y,y+1]: for by in [y, y+1]:
for bz in [z,z+1]: for bz in [z, z+1]:
key = (bx,by,bz) key = (bx, by, bz)
try: try:
self.boxes[key].append(atom) self.boxes[key].append(atom)
except KeyError: except KeyError:
# No box exists for this coordinate
pass pass
#info(atom,'->',key,':',len(self.boxes[key]))
return
def has_bond(self, atom1, atom2): def has_bond(self, atom1, atom2):
"""Look for bond between two atoms.
Args:
atom1: first atom to check
atom2: second atom to check
Returns:
True if there is a bond between atoms
"""
if atom1 in atom2.bonded_atoms or atom2 in atom1.bonded_atoms: if atom1 in atom2.bonded_atoms or atom2 in atom1.bonded_atoms:
return True return True
return False return False
def make_bond(self, atom1, atom2): def make_bond(self, atom1, atom2):
""" Makes a bond between atom1 and atom2 """ """Makes a bond between atom1 and atom2
Args:
atom1: first atom to bond
atom2: second atom to bond
"""
if atom1 == atom2: if atom1 == atom2:
return return
#info('making bond for',atom1,atom2)
if not atom1 in atom2.bonded_atoms: if not atom1 in atom2.bonded_atoms:
atom2.bonded_atoms.append(atom1) atom2.bonded_atoms.append(atom1)
if not atom2 in atom1.bonded_atoms: if not atom2 in atom1.bonded_atoms:
atom1.bonded_atoms.append(atom2) atom1.bonded_atoms.append(atom2)
return
def generate_protein_bond_dictionary(self, atoms): def generate_protein_bond_dictionary(self, atoms):
"""Generate dictionary of protein bonds.
Args:
atoms: list of atoms for bonding
"""
for atom in atoms: for atom in atoms:
for bonded_atom in atom.bonded_atoms: for bonded_atom in atom.bonded_atoms:
resi_i = atom.resName resi_i = atom.res_name
name_i = atom.name name_i = atom.name
resi_j = bonded_atom.resName resi_j = bonded_atom.res_name
name_j = bonded_atom.name name_j = bonded_atom.name
if not name_i in self.backbone_atoms or\ if not name_i in self.backbone_atoms or\
not name_j in self.backbone_atoms: not name_j in self.backbone_atoms:
if not name_i in self.terminal_oxygen_names and\ if not name_i in self.terminal_oxygen_names and\
not name_j in self.terminal_oxygen_names: not name_j in self.terminal_oxygen_names:
if not resi_i in list(self.protein_bonds.keys()): if not resi_i in list(self.protein_bonds.keys()):
self.protein_bonds[resi_i] = {} self.protein_bonds[resi_i] = {}
if not name_i in self.protein_bonds[resi_i]: if not name_i in self.protein_bonds[resi_i]:
self.protein_bonds[resi_i][name_i] = [] self.protein_bonds[resi_i][name_i] = []
if not name_j in self.protein_bonds[resi_i][name_i]: if not name_j in self.protein_bonds[resi_i][name_i]:
self.protein_bonds[resi_i][name_i].append(name_j) self.protein_bonds[resi_i][name_i].append(name_j)
if not resi_j in list(self.protein_bonds.keys()): if not resi_j in list(self.protein_bonds.keys()):
self.protein_bonds[resi_j] = {} self.protein_bonds[resi_j] = {}
if not name_j in self.protein_bonds[resi_j]: if not name_j in self.protein_bonds[resi_j]:
self.protein_bonds[resi_j][name_j] = [] self.protein_bonds[resi_j][name_j] = []
if not name_i in self.protein_bonds[resi_j][name_j]: if not name_i in self.protein_bonds[resi_j][name_j]:
self.protein_bonds[resi_j][name_j].append(name_i) self.protein_bonds[resi_j][name_j].append(name_i)
return
if __name__ == '__main__':
# If called directly, set up protein bond dictionary
import protein, pdb, sys,os
arguments = sys.argv
if len(arguments) != 2:
info('Usage: bonds.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)
pdblist = pdb.readPDB(filename)
my_protein = protein.Protein(pdblist,'test.pdb')
for chain in my_protein.chains:
for residue in chain.residues:
residue.atoms = [atom for atom in residue.atoms if atom.element != 'H']
b = bondmaker()
#b.protein_bonds = {}
atoms = b.find_bonds_for_protein_by_distance(my_protein)
# b.generate_protein_bond_dictionary(atoms)
#file = open(b.data_file_name,'wb')
#pickle.dump(b.protein_bonds, file)
#file.close()

123
propka/calculations.py
View File

@@ -1,98 +1,53 @@
"""PROPKA calculations."""
from __future__ import division import math
from __future__ import print_function import copy
import sys
import math, propka.protonate, propka.bonds,copy, sys import propka.protonate
import propka.bonds
from propka.lib import info, warning from propka.lib import info, warning
#
# methods for setting up atom naming, bonding, and protonation
#
def setup_bonding_and_protonation(parameters, molecular_container): def setup_bonding_and_protonation(parameters, molecular_container):
"""Set up bonding and protonation for a molecule.
Args:
molecular_container: molecule container.
"""
# make bonds # make bonds
my_bond_maker = setup_bonding(parameters, molecular_container) my_bond_maker = setup_bonding(parameters, molecular_container)
# set up ligand atom names # set up ligand atom names
set_ligand_atom_names(molecular_container) set_ligand_atom_names(molecular_container)
# apply information on pi electrons # apply information on pi electrons
my_bond_maker.add_pi_electron_information(molecular_container) my_bond_maker.add_pi_electron_information(molecular_container)
# Protonate atoms # Protonate atoms
if molecular_container.options.protonate_all: if molecular_container.options.protonate_all:
my_protonator = propka.protonate.Protonate(verbose=False) my_protonator = propka.protonate.Protonate(verbose=False)
my_protonator.protonate(molecular_container) my_protonator.protonate(molecular_container)
return
def setup_bonding(parameters, molecular_container): def setup_bonding(parameters, molecular_container):
# make bonds """Set up bonding for a molecular container.
my_bond_maker = propka.bonds.bondmaker()
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
Args:
molecular_container: the molecular container in question
Returns:
BondMaker object
"""
my_bond_maker = propka.bonds.BondMaker()
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
return my_bond_maker return my_bond_maker
def set_ligand_atom_names(molecular_container): def set_ligand_atom_names(molecular_container):
"""Set names for ligands in molecular container.
Args:
molecular_container: molecular container for ligand names
"""
for name in molecular_container.conformation_names: for name in molecular_container.conformation_names:
molecular_container.conformations[name].set_ligand_atom_names() molecular_container.conformations[name].set_ligand_atom_names()
return
#
# The following methods imitates propka3.0 protonation!
#
def setup_bonding_and_protonation_30_style(parameters, molecular_container):
# Protonate atoms
protonate_30_style(molecular_container)
# make bonds
my_bond_maker = propka.bonds.bondmaker()
my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container)
return
def protonate_30_style(molecular_container):
for name in molecular_container.conformation_names:
info('Now protonating', name)
# split atom into residues
curres = -1000000
residue = []
O=None
C=None
for atom in molecular_container.conformations[name].atoms:
if atom.resNumb != curres:
curres = atom.resNumb
if len(residue)>0:
#backbone
[O, C]= addBackBoneHydrogen(residue,O,C)
#arginine
if residue[0].resName == 'ARG':
addArgHydrogen(residue)
#histidine
if residue[0].resName == 'HIS':
addHisHydrogen(residue)
#tryptophan
if residue[0].resName == 'TRP':
addTrpHydrogen(residue)
#amides
if residue[0].resName in ['GLN','ASN']:
addAmdHydrogen(residue)
residue = []
if atom.type=='atom':
residue.append(atom)
return
def addArgHydrogen(residue): def addArgHydrogen(residue):
""" """
Adds Arg hydrogen atoms to residues according to the 'old way'. Adds Arg hydrogen atoms to residues according to the 'old way'.
@@ -159,7 +114,7 @@ def addTrpHydrogen(residue):
elif atom.name == "CE2": elif atom.name == "CE2":
CE = atom CE = atom
if CD == None or NE == None or CE == None: if CD == None or NE == None or CE == None:
str = "Did not find all atoms in %s%4d - in %s" % (self.resName, self.resNumb, "addTrpHydrogen()") str = "Did not find all atoms in %s%4d - in %s" % (self.res_name, self.res_num, "addTrpHydrogen()")
info(str) info(str)
sys.exit(0) sys.exit(0)
@@ -176,15 +131,15 @@ def addAmdHydrogen(residue):
O = None O = None
N = None N = None
for atom in residue: for atom in residue:
if (atom.resName == "GLN" and atom.name == "CD") or (atom.resName == "ASN" and atom.name == "CG"): if (atom.res_name == "GLN" and atom.name == "CD") or (atom.res_name == "ASN" and atom.name == "CG"):
C = atom C = atom
elif (atom.resName == "GLN" and atom.name == "OE1") or (atom.resName == "ASN" and atom.name == "OD1"): elif (atom.res_name == "GLN" and atom.name == "OE1") or (atom.res_name == "ASN" and atom.name == "OD1"):
O = atom O = atom
elif (atom.resName == "GLN" and atom.name == "NE2") or (atom.resName == "ASN" and atom.name == "ND2"): elif (atom.res_name == "GLN" and atom.name == "NE2") or (atom.res_name == "ASN" and atom.name == "ND2"):
N = atom N = atom
if C == None or O == None or N == None: if C == None or O == None or N == None:
str = "Did not find N, C and/or O in %s%4d - in %s" % (atom.resName, atom.resNumb, "addAmdHydrogen()") str = "Did not find N, C and/or O in %s%4d - in %s" % (atom.res_name, atom.res_num, "addAmdHydrogen()")
info(str) info(str)
sys.exit(0) sys.exit(0)
@@ -222,7 +177,7 @@ def addBackBoneHydrogen(residue, O, C):
return [new_O,new_C] return [new_O,new_C]
if N.resName == "PRO": if N.res_name == "PRO":
""" PRO doesn't have an H-atom; do nothing """ """ PRO doesn't have an H-atom; do nothing """
else: else:
H = protonateDirection([N, O, C]) H = protonateDirection([N, O, C])
@@ -311,11 +266,11 @@ def protonateSP2(list):
def make_new_H(atom, x,y,z): def make_new_H(atom, x,y,z):
new_H = propka.atom.Atom() new_H = propka.atom.Atom()
new_H.setProperty(numb = None, new_H.set_property(numb = None,
name = 'H%s'%atom.name[1:], name = 'H%s'%atom.name[1:],
resName = atom.resName, res_name = atom.res_name,
chainID = atom.chainID, chain_id = atom.chain_id,
resNumb = atom.resNumb, res_num = atom.res_num,
x = x, x = x,
y = y, y = y,
z = z, z = z,
@@ -329,7 +284,7 @@ def make_new_H(atom, x,y,z):
new_H.steric_number = 0 new_H.steric_number = 0
new_H.number_of_lone_pairs = 0 new_H.number_of_lone_pairs = 0
new_H.number_of_protons_to_add = 0 new_H.number_of_protons_to_add = 0
new_H.number_of_pi_electrons_in_double_and_triple_bonds = 0 new_H.num_pi_elec_2_3_bonds = 0
atom.bonded_atoms.append(new_H) atom.bonded_atoms.append(new_H)
atom.conformation_container.add_atom(new_H) atom.conformation_container.add_atom(new_H)
@@ -357,7 +312,7 @@ def radial_volume_desolvation(parameters, group):
for atom in all_atoms: for atom in all_atoms:
# ignore atoms in the same residue # ignore atoms in the same residue
if atom.resNumb == group.atom.resNumb and atom.chainID == group.atom.chainID: if atom.res_num == group.atom.res_num and atom.chain_id == group.atom.chain_id:
continue continue
sq_dist = squared_distance(group, atom) sq_dist = squared_distance(group, atom)
@@ -409,15 +364,15 @@ def contactDesolvation(parameters, group):
'N+': 4.5} 'N+': 4.5}
all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms() all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms()
if residue.resName in version.desolvationRadii: if residue.res_name in version.desolvationRadii:
local_cutoff = version.desolvationRadii[residue.resName] local_cutoff = version.desolvationRadii[residue.res_name]
else: else:
local_cutoff = 0.00 local_cutoff = 0.00
residue.Nmass = 0 residue.Nmass = 0
residue.Nlocl = 0 residue.Nlocl = 0
for atom in all_atoms: for atom in all_atoms:
if atom.resNumb != group.atom.resNumb or atom.chainID != group.atom.chainID: if atom.res_num != group.atom.res_num or atom.chain_id != group.atom.chain_id:
dX = atom.x - residue.x dX = atom.x - residue.x
dY = atom.y - residue.y dY = atom.y - residue.y
dZ = atom.z - residue.z dZ = atom.z - residue.z

14
propka/conformation_container.py
View File

@@ -158,7 +158,7 @@ class Conformation_container:
titrate_only = self.molecular_container.options.titrate_only titrate_only = self.molecular_container.options.titrate_only
if titrate_only is not None: if titrate_only is not None:
at = group.atom at = group.atom
if not (at.chainID, at.resNumb, at.icode) in titrate_only: if not (at.chain_id, at.res_num, at.icode) in titrate_only:
group.titratable = False group.titratable = False
if group.residue_type == 'CYS': if group.residue_type == 'CYS':
group.exclude_cys_from_results = True group.exclude_cys_from_results = True
@@ -376,7 +376,7 @@ class Conformation_container:
return [atom for atom in self.atoms if atom.type=='hetatm' and atom.element != 'H'] return [atom for atom in self.atoms if atom.type=='hetatm' and atom.element != 'H']
def get_chain(self,chain): def get_chain(self,chain):
return [atom for atom in self.atoms if atom.chainID != chain] return [atom for atom in self.atoms if atom.chain_id != chain]
def add_atom(self, atom): def add_atom(self, atom):
@@ -388,13 +388,13 @@ class Conformation_container:
atom.molecular_container = self.molecular_container atom.molecular_container = self.molecular_container
# store chain id for bookkeeping # store chain id for bookkeeping
if not atom.chainID in self.chains: if not atom.chain_id in self.chains:
self.chains.append(atom.chainID) self.chains.append(atom.chain_id)
return return
def copy_atom(self, atom): def copy_atom(self, atom):
new_atom = atom.makeCopy() new_atom = atom.make_copy()
self.atoms.append(new_atom) self.atoms.append(new_atom)
new_atom.conformation_container = self new_atom.conformation_container = self
@@ -443,8 +443,8 @@ class Conformation_container:
return return
def sort_atoms_key(self, atom): def sort_atoms_key(self, atom):
key = ord(atom.chainID)*1e7 key = ord(atom.chain_id)*1e7
key += atom.resNumb*1000 key += atom.res_num*1000
if len(atom.name) > len(atom.element): if len(atom.name) > len(atom.element):
key += ord(atom.name[len(atom.element)]) key += ord(atom.name[len(atom.element)])
#info(atom,ord(atom.name[len(atom.element)]), '|%s||%s|'%(atom.name,atom.element)) #info(atom,ord(atom.name[len(atom.element)]), '|%s||%s|'%(atom.name,atom.element))

2
propka/determinants.py
View File

@@ -54,7 +54,7 @@ def addDeterminants(group1, group2, distance, version):
def addSidechainDeterminants(group1, group2, version=None): def addSidechainDeterminants(group1, group2, version=None):
""" """
adding side-chain determinants/perturbations adding side-chain determinants/perturbations
Note, resNumb1 > resNumb2 Note, res_num1 > res_num2
""" """
hbond_interaction = version.hydrogen_bond_interaction(group1, group2) hbond_interaction = version.hydrogen_bond_interaction(group1, group2)

26
propka/group.py
View File

@@ -111,22 +111,22 @@ class Group:
self.common_charge_centre = False self.common_charge_centre = False
self.residue_type = self.atom.resName self.residue_type = self.atom.res_name
if self.atom.terminal: if self.atom.terminal:
self.residue_type = self.atom.terminal self.residue_type = self.atom.terminal
if self.atom.type=='atom': if self.atom.type=='atom':
self.label = '%-3s%4d%2s'%(self.residue_type, atom.resNumb, atom.chainID) self.label = '%-3s%4d%2s'%(self.residue_type, atom.res_num, atom.chain_id)
elif self.atom.resName in ['DA ','DC ','DG ','DT ']: elif self.atom.res_name in ['DA ','DC ','DG ','DT ']:
self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1],
atom.element, atom.element,
atom.name.replace('\'','')[-1], atom.name.replace('\'','')[-1],
atom.resNumb, atom.res_num,
atom.chainID) atom.chain_id)
# self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.resNumb, atom.chainID) # self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.res_num, atom.chain_id)
else: else:
self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chainID) self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chain_id)
# container for squared distances # container for squared distances
@@ -253,7 +253,7 @@ class Group:
return self.label==other.label return self.label==other.label
else: else:
# For heterogene atoms we also need to check the residue number # For heterogene atoms we also need to check the residue number
return self.label==other.label and self.atom.resNumb == other.atom.resNumb return self.label==other.label and self.atom.res_num == other.atom.res_num
def __hash__(self): def __hash__(self):
""" Needed together with __eq__ - otherwise we can't make sets of groups """ """ Needed together with __eq__ - otherwise we can't make sets of groups """
@@ -361,7 +361,7 @@ class Group:
if not self.model_pka_set: if not self.model_pka_set:
self.model_pka = self.parameters.model_pkas[self.residue_type] self.model_pka = self.parameters.model_pkas[self.residue_type]
# check if we should apply a custom model pka # check if we should apply a custom model pka
key = '%s-%s'%(self.atom.resName.strip(), self.atom.name.strip()) key = '%s-%s'%(self.atom.res_name.strip(), self.atom.name.strip())
if key in self.parameters.custom_model_pkas.keys(): if key in self.parameters.custom_model_pkas.keys():
self.model_pka = self.parameters.custom_model_pkas[key] self.model_pka = self.parameters.custom_model_pkas[key]
@@ -1122,7 +1122,7 @@ class Ion_group(Group):
def __init__(self, atom): def __init__(self, atom):
Group.__init__(self,atom) Group.__init__(self,atom)
self.type = 'ION' self.type = 'ION'
self.residue_type = atom.resName.strip() self.residue_type = atom.res_name.strip()
info('Found ion group:', atom) info('Found ion group:', atom)
return return
@@ -1201,7 +1201,7 @@ def is_protein_group(parameters,atom):
### Backbone ### Backbone
if atom.type == 'atom' and atom.name == 'N': if atom.type == 'atom' and atom.name == 'N':
# ignore proline backbone nitrogens # ignore proline backbone nitrogens
if atom.resName != 'PRO': if atom.res_name != 'PRO':
return BBN_group(atom) return BBN_group(atom)
if atom.type == 'atom' and atom.name == 'C': if atom.type == 'atom' and atom.name == 'C':
# ignore C- carboxyl # ignore C- carboxyl
@@ -1209,7 +1209,7 @@ def is_protein_group(parameters,atom):
return BBC_group(atom) return BBC_group(atom)
### Filters for side chains based on PDB protein atom names ### Filters for side chains based on PDB protein atom names
key = '%s-%s'%(atom.resName, atom.name) key = '%s-%s'%(atom.res_name, atom.name)
if key in parameters.protein_group_mapping.keys(): if key in parameters.protein_group_mapping.keys():
return eval('%s_group(atom)'%parameters.protein_group_mapping[key]) return eval('%s_group(atom)'%parameters.protein_group_mapping[key])
@@ -1342,7 +1342,7 @@ def is_ligand_group_by_marvin_pkas(parameters, atom):
def is_ion_group(parameters, atom): def is_ion_group(parameters, atom):
if atom.resName.strip() in parameters.ions.keys(): if atom.res_name.strip() in parameters.ions.keys():
return Ion_group(atom) return Ion_group(atom)
return None return None

4
propka/hybrid36.py
View File

@@ -1,3 +1,7 @@
"""Provides an alternative PDB format that can transparently encode larger atom numbers.
http://cci.lbl.gov/hybrid_36/
"""
import string import string
_hybrid36_upper_chars = set(string.ascii_uppercase) _hybrid36_upper_chars = set(string.ascii_uppercase)

12
propka/iterative.py
View File

@@ -124,9 +124,9 @@ def addIterativeIonPair(object1, object2, interaction, version):
Q2 = object2.Q Q2 = object2.Q
comp1 = object1.pKa_old + annihilation[0] + Q1*coulomb_value comp1 = object1.pKa_old + annihilation[0] + Q1*coulomb_value
comp2 = object2.pKa_old + annihilation[1] + Q2*coulomb_value comp2 = object2.pKa_old + annihilation[1] + Q2*coulomb_value
if object1.resName not in version.parameters.exclude_sidechain_interactions: if object1.res_name not in version.parameters.exclude_sidechain_interactions:
comp1 += Q1*hbond_value comp1 += Q1*hbond_value
if object2.resName not in version.parameters.exclude_sidechain_interactions: if object2.res_name not in version.parameters.exclude_sidechain_interactions:
comp2 += Q2*hbond_value comp2 += Q2*hbond_value
if Q1 == -1.0 and comp1 < comp2: if Q1 == -1.0 and comp1 < comp2:
@@ -155,12 +155,12 @@ def addIterativeIonPair(object1, object2, interaction, version):
# Side-chain # Side-chain
if hbond_value > 0.005: if hbond_value > 0.005:
# residue1 # residue1
if object1.resName not in version.parameters.exclude_sidechain_interactions: if object1.res_name not in version.parameters.exclude_sidechain_interactions:
interaction = [object2, Q1*hbond_value] interaction = [object2, Q1*hbond_value]
annihilation[0] += -Q1*hbond_value annihilation[0] += -Q1*hbond_value
object1.determinants['sidechain'].append(interaction) object1.determinants['sidechain'].append(interaction)
# residue2 # residue2
if object2.resName not in version.parameters.exclude_sidechain_interactions: if object2.res_name not in version.parameters.exclude_sidechain_interactions:
interaction = [object1, Q2*hbond_value] interaction = [object1, Q2*hbond_value]
annihilation[1] += -Q2*hbond_value annihilation[1] += -Q2*hbond_value
object2.determinants['sidechain'].append(interaction) object2.determinants['sidechain'].append(interaction)
@@ -309,7 +309,7 @@ class Iterative:
self.label = group.label self.label = group.label
self.atom = group.atom self.atom = group.atom
self.resName = group.residue_type self.res_name = group.residue_type
self.Q = group.charge self.Q = group.charge
self.pKa_old = None self.pKa_old = None
self.pKa_new = None self.pKa_new = None
@@ -357,7 +357,7 @@ class Iterative:
return self.label==other.label return self.label==other.label
else: else:
# For heterogene atoms we also need to check the residue number # For heterogene atoms we also need to check the residue number
return self.label==other.label and self.atom.resNumb == other.atom.resNumb return self.label==other.label and self.atom.res_num == other.atom.res_num
def __hash__(self): def __hash__(self):
""" Needed together with __eq__ - otherwise we can't make sets of groups """ """ Needed together with __eq__ - otherwise we can't make sets of groups """

10
propka/output.py
View File

@@ -38,13 +38,13 @@ def writePDB(protein, file=None, filename=None, include_hydrogens=False, options
numb = 0 numb = 0
for chain in protein.chains: for chain in protein.chains:
for residue in chain.residues: for residue in chain.residues:
if residue.resName not in ["N+ ", "C- "]: if residue.res_name not in ["N+ ", "C- "]:
for atom in residue.atoms: for atom in residue.atoms:
if include_hydrogens == False and atom.name[0] == "H": if include_hydrogens == False and atom.name[0] == "H":
""" don't print """ """ don't print """
else: else:
numb += 1 numb += 1
line = atom.makePDBLine(numb=numb) line = atom.make_pdb_line2(numb=numb)
line += "\n" line += "\n"
file.write(line) file.write(line)
@@ -130,7 +130,7 @@ def getDeterminantSection(protein, conformation, parameters):
# printing determinants # printing determinants
for chain in protein.conformations[conformation].chains: for chain in protein.conformations[conformation].chains:
for residue_type in parameters.write_out_order: for residue_type in parameters.write_out_order:
groups = [g for g in protein.conformations[conformation].groups if g.atom.chainID == chain] groups = [g for g in protein.conformations[conformation].groups if g.atom.chain_id == chain]
for group in groups: for group in groups:
if group.residue_type == residue_type: if group.residue_type == residue_type:
str += "%s" % ( group.getDeterminantString(parameters.remove_penalised_group) ) str += "%s" % ( group.getDeterminantString(parameters.remove_penalised_group) )
@@ -226,8 +226,8 @@ def writeJackalScapFile(mutationData=None, filename="1xxx_scap.list", options=No
""" """
file = open(filename, 'w') file = open(filename, 'w')
for chainID, code1, resNumb, code2 in mutationData: for chain_id, code1, res_num, code2 in mutationData:
str = "%s, %d, %s\n" % (chainID, resNumb, code2) str = "%s, %d, %s\n" % (chain_id, res_num, code2)
file.write(str) file.write(str)
file.close() file.close()

18
propka/pdb.py
View File

@@ -63,27 +63,27 @@ def protein_precheck(conformations, names):
atoms_by_residue[res_id] = [a] atoms_by_residue[res_id] = [a]
for res_id, res_atoms in atoms_by_residue.items(): for res_id, res_atoms in atoms_by_residue.items():
resname = res_atoms[0].resName res_name = res_atoms[0].res_name
residue_label = '%3s%5s'%(resname, res_id) residue_label = '%3s%5s'%(res_name, res_id)
# ignore ligand residues # ignore ligand residues
if resname not in expected_atom_numbers: if res_name not in expected_atom_numbers:
continue continue
# check for c-terminal # check for c-terminal
if 'C-' in [a.terminal for a in res_atoms]: if 'C-' in [a.terminal for a in res_atoms]:
if len(res_atoms) != expected_atom_numbers[resname]+1: if len(res_atoms) != expected_atom_numbers[res_name]+1:
warning('Unexpected number (%d) of atoms in residue %s in conformation %s' % (len(res_atoms), residue_label, name)) warning('Unexpected number (%d) of atoms in residue %s in conformation %s' % (len(res_atoms), residue_label, name))
continue continue
# check number of atoms in residue # check number of atoms in residue
if len(res_atoms) != expected_atom_numbers[resname]: if len(res_atoms) != expected_atom_numbers[res_name]:
warning('Unexpected number (%d) of atoms in residue %s in conformation %s' % (len(res_atoms), residue_label, name)) warning('Unexpected number (%d) of atoms in residue %s in conformation %s' % (len(res_atoms), residue_label, name))
return return
def resid_from_atom(a): def resid_from_atom(a):
return '%4d %s %s'%(a.resNumb,a.chainID,a.icode) return '%4d %s %s'%(a.res_num,a.chain_id,a.icode)
def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, tags = ['ATOM ', 'HETATM'], chains=None): def get_atom_lines_from_pdb(pdb_file, ignore_residues = [], keep_protons=False, tags = ['ATOM ', 'HETATM'], chains=None):
@@ -197,7 +197,7 @@ def write_mol2_for_atoms(atoms, filename):
substructure_section = '@<TRIPOS>SUBSTRUCTURE\n\n' substructure_section = '@<TRIPOS>SUBSTRUCTURE\n\n'
if len(atoms)>0: if len(atoms)>0:
substructure_section = '@<TRIPOS>SUBSTRUCTURE\n%-7d %10s %7d\n'%(atoms[0].resNumb,atoms[0].resName,atoms[0].numb) substructure_section = '@<TRIPOS>SUBSTRUCTURE\n%-7d %10s %7d\n'%(atoms[0].res_num,atoms[0].res_name,atoms[0].numb)
out = propka.lib.open_file_for_writing(filename) out = propka.lib.open_file_for_writing(filename)
out.write(header%(len(atoms),id-1)) out.write(header%(len(atoms),id-1))
@@ -210,8 +210,8 @@ def write_mol2_for_atoms(atoms, filename):
def get_bond_order(atom1, atom2): def get_bond_order(atom1, atom2):
type = '1' type = '1'
pi_electrons1 = atom1.number_of_pi_electrons_in_double_and_triple_bonds pi_electrons1 = atom1.num_pi_elec_2_3_bonds
pi_electrons2 = atom2.number_of_pi_electrons_in_double_and_triple_bonds pi_electrons2 = atom2.num_pi_elec_2_3_bonds
if '.ar' in atom1.sybyl_type: if '.ar' in atom1.sybyl_type:
pi_electrons1 -=1 pi_electrons1 -=1

32
propka/protonate.py
View File

@@ -121,7 +121,7 @@ class Protonate:
def set_charge(self, atom): def set_charge(self, atom):
# atom is a protein atom # atom is a protein atom
if atom.type=='atom': if atom.type=='atom':
key = '%3s-%s'%(atom.resName, atom.name) key = '%3s-%s'%(atom.res_name, atom.name)
if atom.terminal: if atom.terminal:
debug(atom.terminal) debug(atom.terminal)
key=atom.terminal key=atom.terminal
@@ -172,8 +172,8 @@ class Protonate:
debug('Valence eletrons: %4d'%-self.valence_electrons[atom.element]) debug('Valence eletrons: %4d'%-self.valence_electrons[atom.element])
atom.number_of_protons_to_add -= len(atom.bonded_atoms) atom.number_of_protons_to_add -= len(atom.bonded_atoms)
debug('Number of bonds: %4d'%- len(atom.bonded_atoms)) debug('Number of bonds: %4d'%- len(atom.bonded_atoms))
atom.number_of_protons_to_add -= atom.number_of_pi_electrons_in_double_and_triple_bonds atom.number_of_protons_to_add -= atom.num_pi_elec_2_3_bonds
debug('Pi electrons: %4d'%-atom.number_of_pi_electrons_in_double_and_triple_bonds) debug('Pi electrons: %4d'%-atom.num_pi_elec_2_3_bonds)
atom.number_of_protons_to_add += int(atom.charge) atom.number_of_protons_to_add += int(atom.charge)
debug('Charge: %4.1f'%atom.charge) debug('Charge: %4.1f'%atom.charge)
@@ -185,7 +185,7 @@ class Protonate:
def set_steric_number_and_lone_pairs(self, atom): def set_steric_number_and_lone_pairs(self, atom):
# If we already did this, there is no reason to do it again # If we already did this, there is no reason to do it again
if atom.steric_number_and_lone_pairs_set: if atom.steric_num_lone_pairs_set:
return return
debug('='*10) debug('='*10)
@@ -211,11 +211,11 @@ class Protonate:
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 atom.steric_number += atom.number_of_protons_to_add
debug('%65s: %4d'%('Number of pi-electrons in double and triple bonds(-)',atom.number_of_pi_electrons_in_double_and_triple_bonds)) debug('%65s: %4d'%('Number of pi-electrons in double and triple bonds(-)',atom.num_pi_elec_2_3_bonds))
atom.steric_number -= atom.number_of_pi_electrons_in_double_and_triple_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.number_of_pi_electrons_in_conjugate_double_and_triple_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.number_of_pi_electrons_in_conjugate_double_and_triple_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 atom.steric_number += 0
@@ -231,7 +231,7 @@ class Protonate:
debug('%65s: %4d'%('Steric number',atom.steric_number)) debug('%65s: %4d'%('Steric number',atom.steric_number))
debug('%65s: %4d'%('Number of lone pairs',atom.number_of_lone_pairs)) debug('%65s: %4d'%('Number of lone pairs',atom.number_of_lone_pairs))
atom.steric_number_and_lone_pairs_set = True atom.steric_num_lone_pairs_set = True
return return
@@ -364,11 +364,11 @@ class Protonate:
def add_proton(self, atom, position): def add_proton(self, atom, position):
# Create the new proton # Create the new proton
new_H = propka.atom.Atom() new_H = propka.atom.Atom()
new_H.setProperty(numb = None, new_H.set_property(numb = None,
name = 'H%s'%atom.name[1:], name = 'H%s'%atom.name[1:],
resName = atom.resName, res_name = atom.res_name,
chainID = atom.chainID, chain_id = atom.chain_id,
resNumb = atom.resNumb, res_num = atom.res_num,
x = round(position.x,3), # round of to three digimal points x = round(position.x,3), # round of to three digimal points
y = round(position.y,3), # to avoid round-off differences y = round(position.y,3), # to avoid round-off differences
z = round(position.z,3), # when input file z = round(position.z,3), # when input file
@@ -382,7 +382,7 @@ class Protonate:
new_H.steric_number = 0 new_H.steric_number = 0
new_H.number_of_lone_pairs = 0 new_H.number_of_lone_pairs = 0
new_H.number_of_protons_to_add = 0 new_H.number_of_protons_to_add = 0
new_H.number_of_pi_electrons_in_double_and_triple_bonds = 0 new_H.num_pi_elec_2_3_bonds = 0
new_H.is_protonates = True new_H.is_protonates = True
atom.bonded_atoms.append(new_H) atom.bonded_atoms.append(new_H)
@@ -390,13 +390,13 @@ class Protonate:
atom.conformation_container.add_atom(new_H) atom.conformation_container.add_atom(new_H)
# update names of all protons on this atom # update names of all protons on this atom
new_H.residue_label = "%-3s%4d%2s" % (new_H.name,new_H.resNumb, new_H.chainID) 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') no_protons = atom.count_bonded_elements('H')
if no_protons > 1: if no_protons > 1:
i = 1 i = 1
for proton in atom.get_bonded_elements('H'): for proton in atom.get_bonded_elements('H'):
proton.name = 'H%s%d'%(atom.name[1:],i) proton.name = 'H%s%d'%(atom.name[1:],i)
proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.resNumb, proton.chainID) proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.res_num, proton.chain_id)
i+=1 i+=1