diff --git a/propka/atom.py b/propka/atom.py index 4b71b2a..b93f051 100644 --- a/propka/atom.py +++ b/propka/atom.py @@ -1,66 +1,74 @@ - -from __future__ import division -from __future__ import print_function - -import string, propka.lib, propka.group - +"""Atom class - contains all atom information found in the PDB file""" +import string +import propka.lib +import propka.group from . import hybrid36 -class Atom: - """ - Atom class - contains all atom information found in the pdbfile - """ + +class Atom(object): + """Atom class - contains all atom information found in the PDB file""" def __init__(self, line=None, verbose=False): + """Initialize Atom object. - self.set_properties(line) - - self.residue_label = "%-3s%4d%2s" % (self.name,self.resNumb, self.chainID) - - self.groups_extracted = 0 + Args: + line: Line from a PDB file to set properties of atom. + verbose: TODO - this does not appear to be used. Can we remove it? + """ + 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_type = None self.number_of_bonded_elements = {} - self.cysteine_bridge = False - self.bonded_atoms = [] self.residue = None self.conformation_container = None self.molecular_container = None self.is_protonated = False - self.steric_number_and_lone_pairs_set = False - + self.steric_num_lone_pairs_set = False self.terminal = None self.charge = 0 self.charge_set = False self.steric_number = 0 self.number_of_lone_pairs = 0 self.number_of_protons_to_add = 0 - self.number_of_pi_electrons_in_double_and_triple_bonds = 0 - self.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = 0 + self.num_pi_elec_2_3_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 self.sybyl_type = '' self.sybyl_assigned = False self.marvin_pka = False - - return - - - - def set_properties(self, line): + """Line from PDB file to set properties of atom. + Args: + line: PDB file line + """ self.name = '' self.numb = 0 self.x = 0.0 self.y = 0.0 self.z = 0.0 - self.resNumb = 0 - self.resName = '' - self.chainID = 'A' + self.res_num = 0 + self.res_name = '' + self.chain_id = 'A' self.type = '' self.occ = '1.0' self.beta = '0.0' @@ -69,19 +77,19 @@ class Atom: if line: self.name = line[12:16].strip() - self.numb = int( hybrid36.decode(line[ 6:11]) ) - self.x = float( line[30:38].strip() ) - self.y = float( line[38:46].strip() ) - self.z = float( line[46:54].strip() ) - self.resNumb = int( line[22:26].strip() ) - self.resName = "%-3s" % (line[17:20].strip()) - self.chainID = line[21] + self.numb = int(hybrid36.decode(line[6:11])) + self.x = float(line[30:38].strip()) + self.y = float(line[38:46].strip()) + self.z = float(line[46:54].strip()) + self.res_num = int(line[22:26].strip()) + self.res_name = "%-3s" % (line[17:20].strip()) + self.chain_id = line[21] # Set chain id to "_" if it is just white space. - if not self.chainID.strip(): - self.chainID = '_' + if not self.chain_id.strip(): + self.chain_id = '_' 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.occ = line[55:60].strip() @@ -92,104 +100,136 @@ class Atom: self.element = line[12:14].strip().strip(string.digits) if len(self.name) == 4: self.element = self.element[0] - if len(self.element)==2: - self.element = '%1s%1s'%(self.element[0], self.element[1].lower()) + if len(self.element) == 2: + 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): - self.group_type = type - return + Args: + type_: group type of atom + """ + self.group_type = type_ def count_bonded_elements(self, element): - res = 0 - for ba in self.bonded_atoms: - if element == ba.element: - res +=1 - return res + """Count number of bonded atoms with same element. + Args: + element: element type for test. + Returns: + number of bonded atoms. + """ + return len(self.get_bonded_elements(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 = [] for ba in self.bonded_atoms: - if ba.element == element: - res.append(ba) + if ba.element == element: + res.append(ba) return res 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): - """ check if is found within bonds of self """ + """Check if is found within 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: if ba == other_atom: return True if max_bonds > cur_bond: if ba.is_atom_within_bond_distance(other_atom, max_bonds, cur_bond+1): return True - 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. - - def setProperty(self, - numb = None, - name = None, - resName = None, - chainID = None, - resNumb = None, - x = None, - y = None, - z = None, - occ = None, - beta = None): + Args: + numb: Atom number + name: Atom name + res_name: residue name + chain_id: chain ID + res_num: residue number + x: atom x-coordinate + y: atom y-coordinate + z: atom z-coordinate + occ: atom occupancy + beta: atom temperature factor """ - 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. """ - - if numb != None: self.numb = numb - if name != None: self.name = name - if resName != None: self.resName = resName - if chainID != None: self.chainID = chainID - if resNumb != None: self.resNumb = resNumb - if x != None: self.x = x - if y != None: self.y = y - if z != None: self.z = z - if occ != None: self.occ = occ - if beta != None: self.beta = beta - - - - def makeCopy(self): - """ - making a copy of this 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 - + new_atom = Atom() + new_atom.type = self.type + new_atom.numb = self.numb + new_atom.name = self.name + new_atom.element = self.element + new_atom.res_name = self.res_name + new_atom.res_num = self.res_num + new_atom.chain_id = self.chain_id + new_atom.x = self.x + new_atom.y = self.y + new_atom.z = self.z + new_atom.occ = self.occ + new_atom.beta = self.beta + new_atom.terminal = self.terminal + new_atom.residue_label = self.residue_label + new_atom.icode = self.icode + return new_atom 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 = '-' model_pka = '-' if self.group: @@ -199,197 +239,168 @@ class Atom: if self.group.titratable: model_pka = '%6.2f'%self.group.model_pka - - - str = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s \n" % (self.type.upper(), - self.numb, - propka.lib.makeTidyAtomLabel(self.name, - self.element), - self.resName, - self.chainID, - self.resNumb, - self.x, - self.y, - self.z, - group, - model_pka) - - - - - return str + str_ = "%-6s%5d %s " % (self.type.upper(), self.numb, + propka.lib.makeTidyAtomLabel(self.name, self.element)) + str_ += "%s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s \n" % (self.res_name, self.chain_id, + self.res_num, self.x, self.y, + self.z, group, model_pka) + return str_ 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 = [] for atom in self.bonded_atoms: bonded.append(atom.numb) bonded.sort() - # write 'em out for b in bonded: res += '%5d'%b res += '\n' return res - def get_input_parameters(self): - """ Method for getting the input parameters stored in the - occupancy and b-factor fields in input files""" - + """Extract the input parameters stored in the occupancy and b-factor + fields in input files""" # Set the group type if self.occ != '-': # make sure to set the terminal - if self.occ in ['N+','C-']: + if self.occ in ['N+', 'C-']: self.terminal = self.occ - # save the ligand group charge - if self.occ =='BLG': - self.charge=+1 - elif self.occ =='ALG': - self.charge=-1 - + if self.occ == 'BLG': + self.charge = +1 + elif self.occ == 'ALG': + self.charge = -1 # generic ions - if self.occ in ['1P','2P','1N','2N']: - self.resName=self.occ - self.occ='Ion' - + if self.occ in ['1P', '2P', '1N', '2N']: + self.res_name = self.occ + self.occ = 'Ion' # correct the group type - self.occ = self.occ.replace('N+','Nterm') - self.occ = self.occ.replace('C-','Cterm') - self.occ = self.occ.replace('ION','Ion') - self.occ = self.occ.replace('ALG','titratable_ligand') - self.occ = self.occ.replace('BLG','titratable_ligand') - self.occ = self.occ.replace('LG','non_titratable_ligand') - + self.occ = self.occ.replace('N+', 'Nterm') + self.occ = self.occ.replace('C-', 'Cterm') + self.occ = self.occ.replace('ION', 'Ion') + self.occ = self.occ.replace('ALG', 'titratable_ligand') + self.occ = self.occ.replace('BLG', 'titratable_ligand') + self.occ = self.occ.replace('LG', 'non_titratable_ligand') # try to initialise the group 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: - 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 if self.beta != '-': self.group.model_pka = float(self.beta) self.group.model_pka_set = True - # set occ and beta to standard values self.occ = '1.00' self.beta = '0.00' - - - return - - - def make_pdb_line(self): + """Create PDB line. - # making string - str = "%-6s%5d %s %s%2s%4d%12.3lf%8.3lf%8.3lf%6s%6s\n" % (self.type.upper(), - 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) + TODO - this could/should be a @property method/attribute + TODO - figure out difference between make_pdb_line, make_input_line, and make_pdb_line2 - - return str - - - - 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: + String with PDB line. """ - returns a pdb ATOM-line for various purposes; - specifying arguments over-writes. + str_ = "%-6s%5d " % (self.type.upper(), self.numb) + 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 - if name == None: name = self.name - if resName == None: resName = self.resName - if chainID == None: chainID = self.chainID - if resNumb == None: resNumb = self.resNumb - if x == None: x = self.x - if y == None: y = self.y - if z == None: z = self.z - if occ == None: occ = self.occ - if beta == None: beta = self.beta + str_ = "%-4d %-4s " % (id_, propka.lib.makeTidyAtomLabel(self.name, + self.element)) + str_ += "%10.4f %10.4f %10.4f " % (self.x, self.y, self.z) + str_ += "%6s %6d %10s %10.4f\n" % (self.sybyl_type.replace('-', ''), + self.res_num, self.res_name, 0.0) + return str_ - # making string - str = "ATOM " - str += "%6d" % (numb) - 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 + def make_pdb_line2(self, numb=None, name=None, res_name=None, chain_id=None, + res_num=None, x=None, y=None, z=None, occ=None, beta=None): + """Create a PDB line. + 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 - """ - return propka.lib.makeTidyAtomLabel(self.name,self.element) + if numb is None: + numb = self.numb + 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): - 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) - - -# def get_element(self): -# """ try to extract element if not already done""" -# if self.element == '': -# self.element = self.name.strip(string.digits)[0] -# return 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, + self.chain_id, self.x, self.y, + self.z, self.element) def set_residue(self, residue): - """ Makes a references to the parent residue""" - if self.residue == None: + """ Makes a reference to the parent residue + + Args: + residue: the parent residue + """ + if self.residue is None: self.residue = residue - - - - diff --git a/propka/bonds.py b/propka/bonds.py index 4d58dec..8e7ef66 100644 --- a/propka/bonds.py +++ b/propka/bonds.py @@ -1,336 +1,334 @@ - -from __future__ import division -from __future__ import print_function - -import pickle,sys,os,math,propka.calculations +"""PROPKA representation of bonds.""" +# TODO - is pickle still used? +import pickle +# TODO - eliminate use of sys +import sys +# TODO - eliminate use of os +import os +import math import json - import pkg_resources +import propka.calculations +# TODO - replace the info/warning imports with logging functionality 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): - # predefined bonding distances - self.distances = { - 'S-S':2.5, - 'F-F':1.7} - + self.distances = {'S-S' : DISULFIDE_DISTANCE, 'F-F' : FLUORIDE_DISTANCE} self.distances_squared = {} for k in self.distances.keys(): - self.distances_squared[k]=self.distances[k]*self.distances[k] - - self.H_dist = 1.5; - self.default_dist = 2.0; - - self.H_dist_squared = self.H_dist * self.H_dist + self.distances_squared[k] = self.distances[k] * self.distances[k] + self.H_dist = HYDROGEN_DISTANCE + self.default_dist = DEFAULT_DISTANCE + self.H_dist_squared = self.H_dist * self.H_dist self.default_dist_squared = self.default_dist * self.default_dist - - self.max_sq_distance = max(list(self.distances_squared.values())+[self.default_dist_squared]) - + distances = list(self.distances_squared.values()) + [self.default_dist_squared] + self.max_sq_distance = max(distances) # protein bonding data 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.intra_residue_backbone_bonds = {'N': ['CA'], - 'CA':['N','C'], - 'C': ['CA','O'], - 'O': ['C']} - - self.number_of_pi_electrons_in_bonds_in_backbone = {'C':1, - 'O':1} - - self.number_of_pi_electrons_in_conjugate_bonds_in_backbone = {'N':1} - - self.number_of_pi_electrons_in_bonds_in_sidechains = {'ARG-CZ' :1, - 'ARG-NH1':1, - 'ASN-OD1':1, - 'ASN-CG' :1, - 'ASP-OD1':1, - 'ASP-CG' :1, - 'GLU-OE1':1, - 'GLU-CD' :1, - 'GLN-OE1':1, - 'GLN-CD' :1, - 'HIS-CG' :1, - 'HIS-CD2':1, - 'HIS-ND1':1, - 'HIS-CE1':1, - 'PHE-CG' :1, - 'PHE-CD1':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.intra_residue_backbone_bonds = {'N': ['CA'], 'CA': ['N', 'C'], + 'C': ['CA', 'O'], 'O': ['C']} + self.num_pi_elec_bonds_backbone = {'C': 1, 'O': 1} + self.num_pi_elec_conj_bonds_backbone = {'N': 1} + self.num_pi_elec_bonds_sidechains = {'ARG-CZ' : 1, 'ARG-NH1': 1, + 'ASN-OD1': 1, 'ASN-CG' : 1, + 'ASP-OD1': 1, 'ASP-CG' : 1, + 'GLU-OE1': 1, 'GLU-CD' : 1, + 'GLN-OE1': 1, 'GLN-CD' : 1, + 'HIS-CG' : 1, 'HIS-CD2': 1, + 'HIS-ND1': 1, 'HIS-CE1': 1, + 'PHE-CG' : 1, 'PHE-CD1': 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.num_pi_elec_conj_bonds_sidechains = {'ARG-NE': 1, 'ARG-NH2': 1, + 'ASN-ND2': 1, 'GLN-NE2': 1, + 'HIS-NE2': 1, 'TRP-NE1': 1} + self.num_pi_elec_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.num_pi_elec_conj_bonds_ligands = {'N.am': 1, 'N.pl3': 1} self.backbone_atoms = list(self.intra_residue_backbone_bonds.keys()) - - self.terminal_oxygen_names = ['OXT','O\'\''] - - - return - + self.terminal_oxygen_names = ['OXT', 'O\'\''] def find_bonds_for_protein(self, protein): - """ Finds bonds proteins based on the way atoms - normally bond in proteins""" + """Finds bonds proteins based on the way atoms normally bond in proteins. + Args: + protein: the protein to search for bonds + """ info('++++ Side chains ++++') # side chains for chain in protein.chains: 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) - info('++++ Backbones ++++') # backbone last_residues = [] for chain in protein.chains: - for i in range(1,len(chain.residues)): - if chain.residues[i-1].resName.replace(' ','') not in ['N+','C-']: - if chain.residues[i].resName.replace(' ','') not in ['N+','C-']: + for i in range(1, len(chain.residues)): + if chain.residues[i-1].res_name.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]) last_residues.append(chain.residues[i]) - info('++++ terminal oxygen ++++') # terminal OXT for last_residue in last_residues: self.find_bonds_for_terminal_oxygen(last_residue) - info('++++ cysteines ++++') # Cysteines for chain in protein.chains: - for i in range(0,len(chain.residues)): - if chain.residues[i].resName == 'CYS': - for j in range(0,len(chain.residues)): - if chain.residues[j].resName == 'CYS' and j != i: + for i in range(0, len(chain.residues)): + if chain.residues[i].res_name == 'CYS': + for j in range(0, len(chain.residues)): + if chain.residues[j].res_name == 'CYS' and j != i: self.check_for_cysteine_bonds(chain.residues[i], chain.residues[j]) - return 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: if atom1.name == 'SG': for atom2 in cys2.atoms: 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) - - return - def find_bonds_for_terminal_oxygen(self, residue): + """Look for bonds for terminal oxygen. + + Args: + residue - test residue + """ for atom1 in residue.atoms: if atom1.name in self.terminal_oxygen_names: for atom2 in residue.atoms: if atom2.name == 'C': self.make_bond(atom1, atom2) - return - - + # TODO - stopped here. def connect_backbone(self, residue1, residue2): - """ Sets up bonds in the backbone """ - # residue 1 + """Sets up bonds in the backbone + + Args: + residue1: first residue to connect + residue2: second residue to connect + """ self.find_bonds_for_residue_backbone(residue1) - - # residue 2 self.find_bonds_for_residue_backbone(residue2) - - # inter-residue bond for atom1 in residue1.atoms: if atom1.name == 'C': for atom2 in residue2.atoms: 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) - return - 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: - if atom1.name in list(self.number_of_pi_electrons_in_bonds_in_backbone.keys()): - atom1.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_in_backbone[atom1.name] - 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 - 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_bonds_backbone.keys()): + atom1.num_pi_elec_2_3_bonds \ + = self.num_pi_elec_bonds_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: for atom2 in residue.atoms: if atom2.name in self.intra_residue_backbone_bonds[atom1.name]: self.make_bond(atom1, atom2) - return - - def find_bonds_for_side_chain(self, atoms): - """ Finds bonds for a side chain """ - for atom1 in atoms: + """Finds bonds for a side chain. - key = '%s-%s'%(atom1.resName,atom1.name) - if key in list(self.number_of_pi_electrons_in_bonds_in_sidechains.keys()): - 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()): - atom1.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_sidechains[key] + Args: + atoms: list of atoms to check for bonds + """ + for atom1 in atoms: + 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.terminal_oxygen_names: for atom2 in atoms: - if atom2.name in self.protein_bonds[atom1.resName][atom1.name]: - self.make_bond(atom1,atom2) - - return - + if atom2.name in self.protein_bonds[atom1.res_name][atom1.name]: + self.make_bond(atom1, atom2) 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 self.find_bonds_for_atoms(ligand.atoms) self.add_pi_electron_table_info(ligand.atoms) - return - 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 for atom in atoms: # for ligands if atom.type == 'hetatm': - if atom.sybyl_type in self.number_of_pi_electrons_in_bonds_ligands.keys(): - atom.number_of_pi_electrons_in_double_and_triple_bonds = self.number_of_pi_electrons_in_bonds_ligands[atom.sybyl_type] - if atom.sybyl_type in self.number_of_pi_electrons_in_conjugate_bonds_in_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] - + if atom.sybyl_type in self.num_pi_elec_bonds_ligands.keys(): + atom.num_pi_elec_2_3_bonds = self.num_pi_elec_bonds_ligands[atom.sybyl_type] + if atom.sybyl_type in self.num_pi_elec_conj_bonds_ligands.keys(): + atom.num_pi_elec_conj_2_3_bonds \ + = self.num_pi_elec_conj_bonds_ligands[atom.sybyl_type] # for protein if atom.type == 'atom': - key = '%s-%s'%(atom.resName,atom.name) - if key in list(self.number_of_pi_electrons_in_bonds_in_sidechains.keys()): - atom.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()): - atom.number_of_pi_electrons_in_conjugate_double_and_triple_bonds = self.number_of_pi_electrons_in_conjugate_bonds_in_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 + key = '%s-%s' % (atom.res_name, atom.name) + if key in list(self.num_pi_elec_bonds_sidechains.keys()): + atom.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()): + atom.num_pi_elec_conj_2_3_bonds = self.num_pi_elec_conj_bonds_sidechains[key] + 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): - """ 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) atoms = [] for chain in molecule.chains: 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: atoms.append(atom) - - self.find_bonds_for_atoms_using_boxes(atoms) ##### - #self.find_bonds_for_atoms(atoms) ##### + self.find_bonds_for_atoms_using_boxes(atoms) return 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) - 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: continue - 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 if atoms[i].element == 'S' and atoms[j].element == 'S': atoms[i].cysteine_bridge = True atoms[j].cysteine_bridge = True - - return - - 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: return False - - key = '%s-%s'%(atom1.element,atom2.element) + key = '%s-%s' % (atom1.element, atom2.element) 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 - if sq_dist < self.default_dist_squared and h_count==0: + if sq_dist < self.default_dist_squared and h_count == 0: return True if key in self.distances_squared.keys(): if sq_dist < self.distances_squared[key]: return True - - - return False 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: 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): + """Add pi electron information to a molecule. + + Args: + molecules: list of molecules for adding pi electron information. + """ for name in molecules.conformation_names: self.add_pi_electron_table_info(molecules.conformations[name].atoms) - return - 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 - - # find min and max coordinates - xmin = 1e6; xmax = -1e6 - ymin = 1e6; ymax = -1e6 - zmin = 1e6; zmax = -1e6 + Args: + atoms: list of atoms for finding bonds + """ + box_size = BOX_SIZE + xmin = POS_MAX + 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: if atom.x > xmax: xmax = atom.x @@ -344,147 +342,95 @@ class bondmaker: ymin = atom.y if atom.z < zmin: zmin = atom.z - - xlen = xmax-xmin - ylen = ymax-ymin - zlen = zmax-zmin - - #info('x range: [%6.2f;%6.2f] %6.2f'%(xmin,xmax,xlen)) - #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 + xlen = xmax - xmin + ylen = ymax - ymin + 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))) + self.num_box_z = max(1, int(math.ceil(zlen/box_size))) self.boxes = {} - for x in range(self.no_box_x): - for y in range(self.no_box_y): - for z in range(self.no_box_z): - self.boxes[(x,y,z)] = [] - - # put atoms into boxes + for x in range(self.num_box_x): + for y in range(self.num_box_y): + for z in range(self.num_box_z): + self.boxes[(x, y, z)] = [] for atom in atoms: - x = math.floor((atom.x-xmin)/box_size) - y = math.floor((atom.y-ymin)/box_size) - z = math.floor((atom.z-zmin)/box_size) - self.put_atom_in_box(x,y,z,atom) - - # assign bonds - for key, value in self.boxes.items(): + x = math.floor((atom.x - xmin)/box_size) + y = math.floor((atom.y - ymin)/box_size) + z = math.floor((atom.z - zmin)/box_size) + self.put_atom_in_box(x, y, z, atom) + for value in self.boxes.values(): self.find_bonds_for_atoms(value) + def put_atom_in_box(self, x, y, z, atom): + """Put an atom in a box. - return - - def put_atom_in_box(self,x,y,z,atom): - # atom in the x,y,z box and the up to 7 neighboring boxes on - # one side of the x,y,z box in each dimension - - for bx in [x,x+1]: - for by in [y,y+1]: - for bz in [z,z+1]: - key = (bx,by,bz) + Args: + x: box x-coordinates + y: box y-coordinates + z: box z-coordinates + atom: the atom to place in a box + """ + for bx in [x, x+1]: + for by in [y, y+1]: + for bz in [z, z+1]: + key = (bx, by, bz) try: self.boxes[key].append(atom) except KeyError: - # No box exists for this coordinate pass - #info(atom,'->',key,':',len(self.boxes[key])) - - return - 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: return True return False - 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: return - - #info('making bond for',atom1,atom2) if not atom1 in atom2.bonded_atoms: atom2.bonded_atoms.append(atom1) - if not atom2 in atom1.bonded_atoms: atom1.bonded_atoms.append(atom2) - return - - def generate_protein_bond_dictionary(self, atoms): + """Generate dictionary of protein bonds. + Args: + atoms: list of atoms for bonding + """ for atom in atoms: for bonded_atom in atom.bonded_atoms: - resi_i = atom.resName + resi_i = atom.res_name name_i = atom.name - resi_j = bonded_atom.resName + resi_j = bonded_atom.res_name name_j = bonded_atom.name - if not name_i in self.backbone_atoms or\ not name_j in self.backbone_atoms: if not name_i in self.terminal_oxygen_names and\ not name_j in self.terminal_oxygen_names: - if not resi_i in list(self.protein_bonds.keys()): self.protein_bonds[resi_i] = {} if not name_i in self.protein_bonds[resi_i]: 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) - - if not resi_j in list(self.protein_bonds.keys()): self.protein_bonds[resi_j] = {} if not name_j in self.protein_bonds[resi_j]: 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) - - - 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 ') - 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() diff --git a/propka/calculations.py b/propka/calculations.py index 4cc075d..ca0f5a6 100644 --- a/propka/calculations.py +++ b/propka/calculations.py @@ -1,98 +1,53 @@ - -from __future__ import division -from __future__ import print_function - -import math, propka.protonate, propka.bonds,copy, sys +"""PROPKA calculations.""" +import math +import copy +import sys +import propka.protonate +import propka.bonds from propka.lib import info, warning -# -# methods for setting up atom naming, bonding, and protonation -# - def setup_bonding_and_protonation(parameters, molecular_container): + """Set up bonding and protonation for a molecule. + + Args: + molecular_container: molecule container. + """ # make bonds my_bond_maker = setup_bonding(parameters, molecular_container) - # set up ligand atom names set_ligand_atom_names(molecular_container) - # apply information on pi electrons my_bond_maker.add_pi_electron_information(molecular_container) - # Protonate atoms if molecular_container.options.protonate_all: my_protonator = propka.protonate.Protonate(verbose=False) my_protonator.protonate(molecular_container) - return - - - def setup_bonding(parameters, molecular_container): - # make bonds - my_bond_maker = propka.bonds.bondmaker() - my_bond_maker.find_bonds_for_molecules_using_boxes(molecular_container) + """Set up bonding for a 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 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: 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): """ Adds Arg hydrogen atoms to residues according to the 'old way'. @@ -159,7 +114,7 @@ def addTrpHydrogen(residue): elif atom.name == "CE2": CE = atom 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) sys.exit(0) @@ -176,15 +131,15 @@ def addAmdHydrogen(residue): O = None N = None 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 - 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 - 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 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) sys.exit(0) @@ -222,7 +177,7 @@ def addBackBoneHydrogen(residue, O, C): return [new_O,new_C] - if N.resName == "PRO": + if N.res_name == "PRO": """ PRO doesn't have an H-atom; do nothing """ else: H = protonateDirection([N, O, C]) @@ -311,11 +266,11 @@ def protonateSP2(list): def make_new_H(atom, x,y,z): new_H = propka.atom.Atom() - new_H.setProperty(numb = None, + new_H.set_property(numb = None, name = 'H%s'%atom.name[1:], - resName = atom.resName, - chainID = atom.chainID, - resNumb = atom.resNumb, + res_name = atom.res_name, + chain_id = atom.chain_id, + res_num = atom.res_num, x = x, y = y, z = z, @@ -329,7 +284,7 @@ def make_new_H(atom, x,y,z): new_H.steric_number = 0 new_H.number_of_lone_pairs = 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.conformation_container.add_atom(new_H) @@ -357,7 +312,7 @@ def radial_volume_desolvation(parameters, group): for atom in all_atoms: # 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 sq_dist = squared_distance(group, atom) @@ -409,15 +364,15 @@ def contactDesolvation(parameters, group): 'N+': 4.5} all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms() - if residue.resName in version.desolvationRadii: - local_cutoff = version.desolvationRadii[residue.resName] + if residue.res_name in version.desolvationRadii: + local_cutoff = version.desolvationRadii[residue.res_name] else: local_cutoff = 0.00 residue.Nmass = 0 residue.Nlocl = 0 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 dY = atom.y - residue.y dZ = atom.z - residue.z diff --git a/propka/conformation_container.py b/propka/conformation_container.py index ce78e9a..ce1a9c7 100644 --- a/propka/conformation_container.py +++ b/propka/conformation_container.py @@ -158,7 +158,7 @@ class Conformation_container: titrate_only = self.molecular_container.options.titrate_only if titrate_only is not None: 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 if group.residue_type == 'CYS': 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'] 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): @@ -388,13 +388,13 @@ class Conformation_container: atom.molecular_container = self.molecular_container # store chain id for bookkeeping - if not atom.chainID in self.chains: - self.chains.append(atom.chainID) + if not atom.chain_id in self.chains: + self.chains.append(atom.chain_id) return def copy_atom(self, atom): - new_atom = atom.makeCopy() + new_atom = atom.make_copy() self.atoms.append(new_atom) new_atom.conformation_container = self @@ -443,8 +443,8 @@ class Conformation_container: return def sort_atoms_key(self, atom): - key = ord(atom.chainID)*1e7 - key += atom.resNumb*1000 + key = ord(atom.chain_id)*1e7 + key += atom.res_num*1000 if len(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)) diff --git a/propka/determinants.py b/propka/determinants.py index 4299432..3154875 100644 --- a/propka/determinants.py +++ b/propka/determinants.py @@ -54,7 +54,7 @@ def addDeterminants(group1, group2, distance, version): def addSidechainDeterminants(group1, group2, version=None): """ adding side-chain determinants/perturbations - Note, resNumb1 > resNumb2 + Note, res_num1 > res_num2 """ hbond_interaction = version.hydrogen_bond_interaction(group1, group2) diff --git a/propka/group.py b/propka/group.py index 475e892..14cc5ec 100644 --- a/propka/group.py +++ b/propka/group.py @@ -111,22 +111,22 @@ class Group: self.common_charge_centre = False - self.residue_type = self.atom.resName + self.residue_type = self.atom.res_name if self.atom.terminal: self.residue_type = self.atom.terminal if self.atom.type=='atom': - self.label = '%-3s%4d%2s'%(self.residue_type, atom.resNumb, atom.chainID) - elif self.atom.resName in ['DA ','DC ','DG ','DT ']: + self.label = '%-3s%4d%2s'%(self.residue_type, atom.res_num, atom.chain_id) + elif self.atom.res_name in ['DA ','DC ','DG ','DT ']: self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element, atom.name.replace('\'','')[-1], - atom.resNumb, - atom.chainID) + atom.res_num, + 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: - 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 @@ -253,7 +253,7 @@ class Group: return self.label==other.label else: # 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): """ Needed together with __eq__ - otherwise we can't make sets of groups """ @@ -361,7 +361,7 @@ class Group: if not self.model_pka_set: self.model_pka = self.parameters.model_pkas[self.residue_type] # 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(): self.model_pka = self.parameters.custom_model_pkas[key] @@ -1122,7 +1122,7 @@ class Ion_group(Group): def __init__(self, atom): Group.__init__(self,atom) self.type = 'ION' - self.residue_type = atom.resName.strip() + self.residue_type = atom.res_name.strip() info('Found ion group:', atom) return @@ -1201,7 +1201,7 @@ def is_protein_group(parameters,atom): ### Backbone if atom.type == 'atom' and atom.name == 'N': # ignore proline backbone nitrogens - if atom.resName != 'PRO': + if atom.res_name != 'PRO': return BBN_group(atom) if atom.type == 'atom' and atom.name == 'C': # ignore C- carboxyl @@ -1209,7 +1209,7 @@ def is_protein_group(parameters,atom): return BBC_group(atom) ### 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(): 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): - if atom.resName.strip() in parameters.ions.keys(): + if atom.res_name.strip() in parameters.ions.keys(): return Ion_group(atom) return None diff --git a/propka/hybrid36.py b/propka/hybrid36.py index f5b96af..e0af2dc 100644 --- a/propka/hybrid36.py +++ b/propka/hybrid36.py @@ -1,3 +1,7 @@ +"""Provides an alternative PDB format that can transparently encode larger atom numbers. + +http://cci.lbl.gov/hybrid_36/ +""" import string _hybrid36_upper_chars = set(string.ascii_uppercase) diff --git a/propka/iterative.py b/propka/iterative.py index 8aee97a..3614d93 100644 --- a/propka/iterative.py +++ b/propka/iterative.py @@ -124,9 +124,9 @@ def addIterativeIonPair(object1, object2, interaction, version): Q2 = object2.Q comp1 = object1.pKa_old + annihilation[0] + Q1*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 - 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 if Q1 == -1.0 and comp1 < comp2: @@ -155,12 +155,12 @@ def addIterativeIonPair(object1, object2, interaction, version): # Side-chain if hbond_value > 0.005: # 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] annihilation[0] += -Q1*hbond_value object1.determinants['sidechain'].append(interaction) # 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] annihilation[1] += -Q2*hbond_value object2.determinants['sidechain'].append(interaction) @@ -309,7 +309,7 @@ class Iterative: self.label = group.label self.atom = group.atom - self.resName = group.residue_type + self.res_name = group.residue_type self.Q = group.charge self.pKa_old = None self.pKa_new = None @@ -357,7 +357,7 @@ class Iterative: return self.label==other.label else: # 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): """ Needed together with __eq__ - otherwise we can't make sets of groups """ diff --git a/propka/output.py b/propka/output.py index 7087233..531af02 100644 --- a/propka/output.py +++ b/propka/output.py @@ -38,13 +38,13 @@ def writePDB(protein, file=None, filename=None, include_hydrogens=False, options numb = 0 for chain in protein.chains: 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: if include_hydrogens == False and atom.name[0] == "H": """ don't print """ else: numb += 1 - line = atom.makePDBLine(numb=numb) + line = atom.make_pdb_line2(numb=numb) line += "\n" file.write(line) @@ -130,7 +130,7 @@ def getDeterminantSection(protein, conformation, parameters): # printing determinants for chain in protein.conformations[conformation].chains: 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: if group.residue_type == residue_type: 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') - for chainID, code1, resNumb, code2 in mutationData: - str = "%s, %d, %s\n" % (chainID, resNumb, code2) + for chain_id, code1, res_num, code2 in mutationData: + str = "%s, %d, %s\n" % (chain_id, res_num, code2) file.write(str) file.close() diff --git a/propka/pdb.py b/propka/pdb.py index a2abfba..0b01588 100644 --- a/propka/pdb.py +++ b/propka/pdb.py @@ -63,27 +63,27 @@ def protein_precheck(conformations, names): atoms_by_residue[res_id] = [a] for res_id, res_atoms in atoms_by_residue.items(): - resname = res_atoms[0].resName - residue_label = '%3s%5s'%(resname, res_id) + res_name = res_atoms[0].res_name + residue_label = '%3s%5s'%(res_name, res_id) # ignore ligand residues - if resname not in expected_atom_numbers: + if res_name not in expected_atom_numbers: continue # check for c-terminal 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)) continue # 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)) return 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): @@ -197,7 +197,7 @@ def write_mol2_for_atoms(atoms, filename): substructure_section = '@SUBSTRUCTURE\n\n' if len(atoms)>0: - substructure_section = '@SUBSTRUCTURE\n%-7d %10s %7d\n'%(atoms[0].resNumb,atoms[0].resName,atoms[0].numb) + substructure_section = '@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.write(header%(len(atoms),id-1)) @@ -210,8 +210,8 @@ def write_mol2_for_atoms(atoms, filename): def get_bond_order(atom1, atom2): type = '1' - pi_electrons1 = atom1.number_of_pi_electrons_in_double_and_triple_bonds - pi_electrons2 = atom2.number_of_pi_electrons_in_double_and_triple_bonds + pi_electrons1 = atom1.num_pi_elec_2_3_bonds + pi_electrons2 = atom2.num_pi_elec_2_3_bonds if '.ar' in atom1.sybyl_type: pi_electrons1 -=1 diff --git a/propka/protonate.py b/propka/protonate.py index 6a4927e..cef3216 100644 --- a/propka/protonate.py +++ b/propka/protonate.py @@ -121,7 +121,7 @@ class Protonate: def set_charge(self, atom): # atom is a protein atom if atom.type=='atom': - key = '%3s-%s'%(atom.resName, atom.name) + key = '%3s-%s'%(atom.res_name, atom.name) if atom.terminal: debug(atom.terminal) key=atom.terminal @@ -172,8 +172,8 @@ class Protonate: 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.number_of_pi_electrons_in_double_and_triple_bonds - debug('Pi electrons: %4d'%-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.num_pi_elec_2_3_bonds) atom.number_of_protons_to_add += int(atom.charge) debug('Charge: %4.1f'%atom.charge) @@ -185,7 +185,7 @@ class Protonate: def set_steric_number_and_lone_pairs(self, atom): # 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 debug('='*10) @@ -211,11 +211,11 @@ class Protonate: 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.number_of_pi_electrons_in_double_and_triple_bonds)) - atom.steric_number -= 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.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)) - atom.steric_number -= 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.num_pi_elec_conj_2_3_bonds debug('%65s: %4d'%('Number of donated co-ordinated bonds',0)) atom.steric_number += 0 @@ -231,7 +231,7 @@ class Protonate: debug('%65s: %4d'%('Steric number',atom.steric_number)) 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 @@ -364,11 +364,11 @@ class Protonate: def add_proton(self, atom, position): # Create the new proton new_H = propka.atom.Atom() - new_H.setProperty(numb = None, + new_H.set_property(numb = None, name = 'H%s'%atom.name[1:], - resName = atom.resName, - chainID = atom.chainID, - resNumb = atom.resNumb, + 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 @@ -382,7 +382,7 @@ class Protonate: new_H.steric_number = 0 new_H.number_of_lone_pairs = 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 atom.bonded_atoms.append(new_H) @@ -390,13 +390,13 @@ class Protonate: 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.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') 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.resNumb, proton.chainID) + proton.residue_label = "%-3s%4d%2s" % (proton.name,proton.res_num, proton.chain_id) i+=1