diff --git a/propka.cfg b/propka.cfg index a12e38a..58b1aec 100644 --- a/propka.cfg +++ b/propka.cfg @@ -342,7 +342,7 @@ common_charge_centre 0 remove_penalised_group 1 # non-covalent coupling -max_intrinsic_pKa_diff 2.0 +max_intrinsic_pka_diff 2.0 min_interaction_energy 0.5 max_free_energy_diff 1.0 min_swap_pka_shift 1.0 diff --git a/propka/calculations.py b/propka/calculations.py index 5bf2449..c2f15b5 100644 --- a/propka/calculations.py +++ b/propka/calculations.py @@ -89,8 +89,8 @@ def setup_bonding_and_protonation(parameters, molecular_container): 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) + PROTONATOR = propka.protonate.Protonate(verbose=False) + PROTONATOR.protonate(molecular_container) def setup_bonding(molecular_container): @@ -386,10 +386,10 @@ def radial_volume_desolvation(parameters, group): """ all_atoms = group.atom.conformation_container.get_non_hydrogen_atoms() volume = 0.0 - # TODO - Nathan really wants to rename the Nmass variable. + # TODO - Nathan really wants to rename the num_volume variable. # He had to re-read the original paper to figure out what it was. # A better name would be num_volume. - group.Nmass = 0 + group.num_volume = 0 min_dist_4th = MIN_DISTANCE_4TH for atom in all_atoms: # ignore atoms in the same residue @@ -410,11 +410,11 @@ def radial_volume_desolvation(parameters, group): volume += dv_inc # buried if sq_dist < parameters.buried_cutoff_squared: - group.Nmass += 1 - group.buried = calculate_weight(parameters, group.Nmass) + group.num_volume += 1 + group.buried = calculate_weight(parameters, group.num_volume) scale_factor = calculate_scale_factor(parameters, group.buried) volume_after_allowance = max(0.00, volume-parameters.desolvationAllowance) - group.Emass = group.charge * parameters.desolvationPrefactor \ + group.energy_volume = group.charge * parameters.desolvationPrefactor \ * volume_after_allowance * scale_factor @@ -589,7 +589,7 @@ def hydrogen_bond_interaction(group1, group2, version): # is closer to the titratable atom than the hydrogen. In either # case, we set the angle factor to 0 f_angle = 0.0 - weight = version.calculate_pair_weight(group1.Nmass, group2.Nmass) + weight = version.calculate_pair_weight(group1.num_volume, group2.num_volume) exception, value = version.check_exceptions(group1, group2) if exception: # Do nothing, value should have been assigned. @@ -614,7 +614,7 @@ def electrostatic_interaction(group1, group2, dist, version): # check if we should do coulomb interaction at all if not version.check_coulomb_pair(group1, group2, dist): return None - weight = version.calculate_pair_weight(group1.Nmass, group2.Nmass) + weight = version.calculate_pair_weight(group1.num_volume, group2.num_volume) value = version.calculate_coulomb_energy(dist, weight) return value @@ -633,7 +633,7 @@ def check_coulomb_pair(parameters, group1, group2, dist): Returns: Boolean """ - num_volume = group1.Nmass + group2.Nmass + num_volume = group1.num_volume + group2.num_volume do_coulomb = True # check if both groups are titratable (ions are taken care of in # determinants::set_ion_determinants) @@ -695,7 +695,7 @@ def backbone_reorganization(parameters, conformation): value = 1.0 - (dist-UNK_BACKBONE_DISTANCE2) \ / (UNK_BACKBONE_DISTANCE1-UNK_BACKBONE_DISTANCE2) dpka += UNK_PKA_SCALING2*min(1.0, value) - titratable_group.Elocl = dpka*weight + titratable_group.energy_local = dpka*weight def check_exceptions(version, group1, group2): @@ -799,7 +799,7 @@ def check_coo_coo_exception(group1, group2, version): closest_atom2) f_angle = 1.00 value = hydrogen_bond_energy(dist, dpka_max, cutoff, f_angle) - weight = calculate_pair_weight(version.parameters, group1.Nmass, group2.Nmass) + weight = calculate_pair_weight(version.parameters, group1.num_volume, group2.num_volume) value = value * (1.0 + weight) return exception, value @@ -816,7 +816,7 @@ def check_coo_his_exception(group1, group2, version): 2. value associated with atypical interaction (None if Boolean is False) """ exception = False - if check_buried(group1.Nmass, group2.Nmass): + if check_buried(group1.num_volume, group2.num_volume): exception = True return exception, version.parameters.COO_HIS_exception @@ -833,7 +833,7 @@ def check_oco_his_exception(group1, group2, version): 2. value associated with atypical interaction (None if Boolean is False) """ exception = False - if check_buried(group1.Nmass, group2.Nmass): + if check_buried(group1.num_volume, group2.num_volume): exception = True return exception, version.parameters.OCO_HIS_exception @@ -850,7 +850,7 @@ def check_cys_his_exception(group1, group2, version): 2. value associated with atypical interaction (None if Boolean is False) """ exception = False - if check_buried(group1.Nmass, group2.Nmass): + if check_buried(group1.num_volume, group2.num_volume): exception = True return exception, version.parameters.CYS_HIS_exception @@ -867,7 +867,7 @@ def check_cys_cys_exception(group1, group2, version): 2. value associated with atypical interaction (None if Boolean is False) """ exception = False - if check_buried(group1.Nmass, group2.Nmass): + if check_buried(group1.num_volume, group2.num_volume): exception = True return exception, version.parameters.CYS_CYS_exception diff --git a/propka/coupled_groups.py b/propka/coupled_groups.py index 49777f8..54197ab 100644 --- a/propka/coupled_groups.py +++ b/propka/coupled_groups.py @@ -33,7 +33,7 @@ class NonCovalentlyCoupledGroups: return {'coupling_factor': -1.0} # calculate intrinsic pKa's, if not already done for group in [group1, group2]: - if not hasattr(group, 'intrinsic_pKa'): + if group.intrinsic_pka is None: group.calculate_intrinsic_pka() use_ph = self.parameters.pH if self.parameters.pH == 'variable': @@ -70,12 +70,12 @@ class NonCovalentlyCoupledGroups: and return_on_fail: return {'coupling_factor': -1.0} # check intrinsic pka diff - if abs(group1.intrinsic_pKa - group2.intrinsic_pKa) \ - > self.parameters.max_intrinsic_pKa_diff and return_on_fail: + if abs(group1.intrinsic_pka - group2.intrinsic_pka) \ + > self.parameters.max_intrinsic_pka_diff and return_on_fail: return {'coupling_factor': -1.0} # if everything is OK, calculate the coupling factor and return all info factor = self.get_free_energy_diff_factor(default_energy, swapped_energy) \ - * self.get_pka_diff_factor(group1.intrinsic_pKa, group2.intrinsic_pKa) \ + * self.get_pka_diff_factor(group1.intrinsic_pka, group2.intrinsic_pka) \ * self.get_interaction_factor(interaction_energy) return {'coupling_factor': factor, 'default_energy': default_energy, 'swapped_energy': swapped_energy, @@ -95,8 +95,8 @@ class NonCovalentlyCoupledGroups: """ intrinsic_pka_diff = abs(pka1-pka2) res = 0.0 - if intrinsic_pka_diff <= self.parameters.max_intrinsic_pKa_diff: - res = 1-(intrinsic_pka_diff/self.parameters.max_intrinsic_pKa_diff)**2 + if intrinsic_pka_diff <= self.parameters.max_intrinsic_pka_diff: + res = 1-(intrinsic_pka_diff/self.parameters.max_intrinsic_pka_diff)**2 return res def get_free_energy_diff_factor(self, energy1, energy2): @@ -147,7 +147,7 @@ class NonCovalentlyCoupledGroups: info(' Detecting non-covalently coupled residues') info('-' * 103) info(' Maximum pKa difference: %4.2f pKa units' \ - % self.parameters.max_intrinsic_pKa_diff) + % self.parameters.max_intrinsic_pka_diff) info(' Minimum interaction energy: %4.2f pKa units' \ % self.parameters.min_interaction_energy) info(' Maximum free energy diff.: %4.2f pKa units' \ @@ -263,7 +263,7 @@ class NonCovalentlyCoupledGroups: str_ = ' ' + '-' * 113 + '\n' for group in system: str_ += self.tagged_format(' %-8s|' % tag, - group.getDeterminantString(), + group.get_determinant_string(), all_labels) return str_ + '\n' @@ -355,8 +355,8 @@ class NonCovalentlyCoupledGroups: (group1.label, group2.label, data['coupling_factor'], data['default_energy'], data['swapped_energy'], data['default_energy'] - data['swapped_energy'], data['pH'], - data['interaction_energy'], group1.intrinsic_pKa, group2.intrinsic_pKa, - group1.intrinsic_pKa-group2.intrinsic_pKa, data['swapped_pka1'], + data['interaction_energy'], group1.intrinsic_pka, group2.intrinsic_pka, + group1.intrinsic_pka-group2.intrinsic_pka, data['swapped_pka1'], data['swapped_pka2'], data['pka_shift1'], data['pka_shift2']) return str_ diff --git a/propka/determinants.py b/propka/determinants.py index 7b75d3f..101547a 100644 --- a/propka/determinants.py +++ b/propka/determinants.py @@ -188,8 +188,8 @@ def set_ion_determinants(conformation_container, version): for ion_group in conformation_container.get_ions(): dist_sq = squared_distance(titratable_group, ion_group) if dist_sq < version.parameters.coulomb_cutoff2_squared: - weight = version.calculate_pair_weight(titratable_group.Nmass, - ion_group.Nmass) + weight = version.calculate_pair_weight(titratable_group.num_volume, + ion_group.num_volume) # the pKa of both acids and bases are shifted up by negative # ions (and vice versa) value = (-ion_group.charge) \ diff --git a/propka/group.py b/propka/group.py index c524dac..49e2552 100644 --- a/propka/group.py +++ b/propka/group.py @@ -1,336 +1,332 @@ -# -# Class for storing groups important for propka calculations -# - -from __future__ import division -from __future__ import print_function - -import propka.ligand, propka.determinant, propka.ligand_pka_values, math, propka.protonate +"""Routines and classes for storing groups important to PROPKA calculations.""" +import math +import propka.ligand +import propka.protonate +from propka.ligand_pka_values import ligand_pka_values +from propka.determinant import Determinant from propka.lib import info, warning -my_protonator = propka.protonate.Protonate(verbose=False) - -expected_atoms_acid_interactions = { - 'COO':{'O':2}, - 'HIS':{'H':2, 'N':2}, - 'CYS':{'S':1}, - 'TYR':{'O':1}, - 'LYS':{'N':1}, - 'ARG':{'H':5, 'N':3}, - 'ROH':{'O':1}, - 'AMD':{'H':2, 'N':1}, - 'TRP':{'H':1, 'N':1}, - 'N+': {'N':1}, - 'C-': {'O':2}, - 'BBN':{'H':1, 'N':1,}, - 'BBC':{'O':1}, - 'NAR':{'H':1, 'N':1}, - 'NAM':{'H':1, 'N':1}, - 'F': {'F':1}, - 'Cl': {'Cl':1}, - 'OH': {'H':1, 'O':1}, - 'OP': {'O':1}, - 'O3': {'O':1}, - 'O2': {'O':1}, - 'SH': {'S':1}, - 'CG': {'H':5, 'N':3}, - 'C2N':{'H':4, 'N':2}, - 'OCO':{'O':2}, - 'N30':{'H':4, 'N':1}, - 'N31':{'H':3, 'N':1}, - 'N32':{'H':2, 'N':1}, - 'N33':{'H':1, 'N':1}, - 'NP1':{'H':2, 'N':1}, - 'N1' :{'N':1} -} - -expected_atoms_base_interactions = { - 'COO':{'O':2}, - 'HIS':{'N':2}, - 'CYS':{'S':1}, - 'TYR':{'O':1}, - 'LYS':{'N':1}, - 'ARG':{'N':3}, - 'ROH':{'O':1}, - 'AMD':{'O':1}, - 'TRP':{'N':1}, - 'N+': {'N':1}, - 'C-': {'O':2}, - 'BBN':{'H':1, 'N':1,}, - 'BBC':{'O':1}, - 'NAR':{'H':1, 'N':1}, - 'NAM':{'H':1, 'N':1}, - 'F': {'F':1}, - 'Cl': {'Cl':1}, - 'OH': {'H':1, 'O':1}, - 'OP': {'O':1}, - 'O3': {'O':1}, - 'O2': {'O':1}, - 'SH': {'S':1}, - 'CG': {'N':3}, - 'C2N':{'N':2}, - 'OCO':{'O':2}, - 'N30':{'N':1}, - 'N31':{'N':1}, - 'N32':{'N':1}, - 'N33':{'N':1}, - 'NP1':{'N':1}, - 'N1' :{'N':1} -} +# Constants that start with "UNK_" are a mystery to me +UNK_PKA_SCALING = -1.36 +PROTONATOR = propka.protonate.Protonate(verbose=False) +EXPECTED_ATOMS_ACID_INTERACTIONS = {'COO': {'O': 2}, 'HIS': {'H': 2, 'N': 2}, + 'CYS': {'S': 1}, 'TYR': {'O': 1}, + 'LYS': {'N': 1}, 'ARG': {'H': 5, 'N': 3}, + 'ROH': {'O': 1}, 'AMD': {'H': 2, 'N': 1}, + 'TRP': {'H': 1, 'N': 1}, 'N+': {'N': 1}, + 'C-': {'O': 2}, 'BBN': {'H': 1, 'N': 1,}, + 'BBC': {'O': 1}, 'NAR': {'H': 1, 'N': 1}, + 'NAM': {'H': 1, 'N': 1}, 'F': {'F': 1}, + 'Cl': {'Cl': 1}, 'OH': {'H': 1, 'O': 1}, + 'OP': {'O': 1}, 'O3': {'O': 1}, + 'O2': {'O': 1}, 'SH': {'S': 1}, + 'CG': {'H': 5, 'N': 3}, + 'C2N': {'H': 4, 'N': 2}, 'OCO': {'O': 2}, + 'N30': {'H': 4, 'N': 1}, + 'N31': {'H': 3, 'N': 1}, + 'N32': {'H': 2, 'N': 1}, + 'N33': {'H': 1, 'N': 1}, + 'NP1': {'H': 2, 'N': 1}, 'N1': {'N': 1}} +EXPECTED_ATOMS_BASE_INTERACTIONS = {'COO': {'O': 2}, 'HIS': {'N': 2}, + 'CYS': {'S': 1}, 'TYR': {'O': 1}, + 'LYS': {'N': 1}, 'ARG': {'N': 3}, + 'ROH': {'O': 1}, 'AMD': {'O': 1}, + 'TRP': {'N': 1}, 'N+': {'N': 1}, + 'C-': {'O': 2}, 'BBN': {'H': 1, 'N': 1}, + 'BBC': {'O': 1}, 'NAR': {'H': 1, 'N': 1}, + 'NAM': {'H': 1, 'N': 1}, 'F': {'F': 1}, + 'Cl': {'Cl': 1}, 'OH': {'H': 1, 'O': 1}, + 'OP': {'O': 1}, 'O3': {'O': 1}, + 'O2': {'O': 1}, 'SH': {'S': 1}, + 'CG': {'N': 3}, 'C2N': {'N': 2}, + 'OCO': {'O': 2}, 'N30': {'N': 1}, + 'N31': {'N': 1}, 'N32': {'N': 1}, + 'N33': {'N': 1}, 'NP1': {'N': 1}, + 'N1': {'N': 1}} class Group: + """Class for storing groups important to pKa calculations.""" + def __init__(self, atom): - #info('Made new %s group from %s'%(type,atom)) + """Initialize with an atom. + + Args: + atom: atom object + """ self.atom = atom self.type = '' atom.group = self - # set up data structures - self.determinants = {'sidechain':[],'backbone':[],'coulomb':[]} + self.determinants = {'sidechain': [], 'backbone': [], 'coulomb': []} self.pka_value = 0.0 self.model_pka = 0.0 - - self.Emass = 0.0 - self.Nmass = 0.0 - self.Elocl = 0.0 - self.Nlocl = 0.0 + # Energy associated with volume interactions + self.energy_volume = 0.0 + # Number of atoms associated with volume interactions + self.num_volume = 0.0 + # Energy associated with local interactions + self.energy_local = 0.0 + # Number of atoms associated with local interactions + self.num_local = 0.0 self.buried = 0.0 self.x = 0.0 self.y = 0.0 self.z = 0.0 self.charge = 0 + self.parameters = None + self.exclude_cys_from_results = None self.interaction_atoms_for_acids = [] self.interaction_atoms_for_bases = [] self.model_pka_set = False - + self.intrinsic_pka = None + self.titratable = None # information on covalent and non-covalent coupling self.non_covalently_coupled_groups = [] self.covalently_coupled_groups = [] self.coupled_titrating_group = None self.common_charge_centre = False - - 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.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.res_num, - atom.chain_id) - -# self.label = '%1s%1s%1s%4d%2s'%(self.residue_type[1], atom.element,atom.name[-1], atom.res_num, atom.chain_id) + if self.atom.type == 'atom': + 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.res_num, + atom.chain_id) else: - self.label = '%-3s%4s%2s'%(self.residue_type, atom.name, atom.chain_id) - - + self.label = '%-3s%4s%2s' % (self.residue_type, atom.name, + atom.chain_id) # container for squared distances self.squared_distances = {} - return - - - # - # Coupling-related methods - # def couple_covalently(self, other): - """ Couple this group with another group """ + """Couple this group with another group. + + Args: + other: other group for coupling + """ # do the coupling if not other in self.covalently_coupled_groups: self.covalently_coupled_groups.append(other) - if not self in other.covalently_coupled_groups: other.covalently_coupled_groups.append(self) - return - def couple_non_covalently(self, other): - """ Couple this group with another group """ + """Non-covalenthly couple this group with another group. + + Args: + other: other group for coupling + """ # do the coupling if not other in self.non_covalently_coupled_groups: self.non_covalently_coupled_groups.append(other) - if not self in other.non_covalently_coupled_groups: other.non_covalently_coupled_groups.append(self) - return + def get_covalently_coupled_groups(self): + """Get covalently coupled groups. - def get_covalently_coupled_groups(self): return self.covalently_coupled_groups - def get_non_covalently_coupled_groups(self): return self.non_covalently_coupled_groups + Returns: + list of covalently coupled groups. + """ + return self.covalently_coupled_groups + def get_non_covalently_coupled_groups(self): + """Get non-covalently coupled groups. + + Returns: + list of covalently coupled groups. + """ + return self.non_covalently_coupled_groups def share_determinants(self, others): + """Share determinants between this group and others. + Args: + others: list of other groups + """ # for each determinant type for other in others: if other == self: + the_other = other continue - - for type in ['sidechain','backbone','coulomb']: - for g in other.determinants[type]: self.share_determinant(g,type) - + for type_ in ['sidechain', 'backbone', 'coulomb']: + for det in other.determinants[type_]: + self.share_determinant(det, type_) # recalculate pka values self.calculate_total_pka() - other.calculate_total_pka() + the_other.calculate_total_pka() - return + def share_determinant(self, new_determinant, type_): + """Add determinant to this group's list of determinants. - - def share_determinant(self, new_determinant, type): + Args: + new_determinant: determinant to add + type_: type of determinant + """ added = False # first check if we already have a determinant with this label - for own_determinant in self.determinants[type]: + for own_determinant in self.determinants[type_]: if own_determinant.group == new_determinant.group: # if so, find the average value - avr = (own_determinant.value + new_determinant.value)/2.0 + avr = 0.5*(own_determinant.value + new_determinant.value) own_determinant.value = avr new_determinant.value = avr added = True - # otherwise we just add the determinant to our list if not added: - self.determinants[type].append(propka.determinant.Determinant(new_determinant.group, - new_determinant.value)) - - return + self.determinants[type_].append(Determinant(new_determinant.group, + new_determinant.value)) def make_covalently_coupled_line(self): + """Create line for covalent coupling. + Returns: + string + """ # first check if there are any coupled groups at all if len(self.covalently_coupled_groups) == 0: return '' - - line = 'CCOUPL%5d'%self.atom.numb - + line = 'CCOUPL%5d' % self.atom.numb # extract and sort numbers of coupled groups coupled = [] - for g in self.covalently_coupled_groups: - coupled.append(g.atom.numb) + for group in self.covalently_coupled_groups: + coupled.append(group.atom.numb) coupled.sort() - # write 'em out - for b in coupled: - line += '%5d'%b + for num in coupled: + line += '%5d' % num line += '\n' - return line def make_non_covalently_coupled_line(self): + """Create line for non-covalent coupling. + + Returns: + string + """ # first check if there are any coupled groups at all if len(self.non_covalently_coupled_groups) == 0: return '' - - line = 'NCOUPL%5d'%self.atom.numb - + line = 'NCOUPL%5d' % self.atom.numb # extract and sort numbers of coupled groups coupled = [] - for g in self.non_covalently_coupled_groups: - coupled.append(g.atom.numb) + for group in self.non_covalently_coupled_groups: + coupled.append(group.atom.numb) coupled.sort() - # write 'em out - for b in coupled: - line += '%5d'%b + for num in coupled: + line += '%5d' % num line += '\n' - - return line - # - # Bookkeeping methods - # def __eq__(self, other): - """ - Check if two groups should be considered identical - """ + """Needed for creating sets of groups.""" if self.atom.type == 'atom': # In case of protein atoms we trust the labels - return self.label==other.label + 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.res_num == other.atom.res_num + 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 """ + """Needed for creating sets of groups.""" return id(self) def __iadd__(self, other): if self.type != other.type: - raise Exception('Cannot add groups of different types (%s and %s)'%(self.type,other.type)) - + errstr = 'Cannot add groups of different types (%s and %s)' \ + % (self.type, other.type) + raise Exception(errstr) # add all values self.pka_value += other.pka_value - self.Nmass += other.Nmass - self.Emass += other.Emass - self.Nlocl += other.Nlocl - self.Elocl += other.Elocl + self.num_volume += other.num_volume + self.energy_volume += other.energy_volume + self.num_local += other.num_local + self.energy_local += other.energy_local self.buried += other.buried # and add all determinants - for type in ['sidechain','backbone','coulomb']: - for determinant in other.determinants[type]: - self.add_determinant(determinant, type) + # TODO - list ['sidechain', 'backbone', 'coulomb'] should be constant + # This list appears all over the code and should be moved to a constant + # higher in the package + for type_ in ['sidechain', 'backbone', 'coulomb']: + for determinant in other.determinants[type_]: + self.add_determinant(determinant, type_) return self + def add_determinant(self, new_determinant, type_): + """Add to current and creates non-present determinants. - def add_determinant(self, new_determinant, type): - """ Adds to current and creates non-present determinants""" + Args: + new_determinant: new determinant to add + type_: determinant type + """ # first check if we already have a determinant with this label - for own_determinant in self.determinants[type]: + for own_determinant in self.determinants[type_]: if own_determinant.group == new_determinant.group: # if so, add the value own_determinant.value += new_determinant.value return - # otherwise we just add the determinant to our list - self.determinants[type].append(propka.determinant.Determinant(new_determinant.group, - new_determinant.value)) + self.determinants[type_].append(Determinant(new_determinant.group, + new_determinant.value)) - return + def set_determinant(self, new_determinant, type_): + """Overwrite current and create non-present determinants. - def set_determinant(self, new_determinant, type): - """ Overwrites current and creates non-present determinants""" + Args: + new_determinant: new determinant to add + type_: determinant type + """ # first check if we already have a determinant with this label - for own_determinant in self.determinants[type]: + for own_determinant in self.determinants[type_]: if own_determinant.group == new_determinant.group: # if so, overwrite the value own_determinant.value = new_determinant.value return - # otherwise we just add the determinant to our list - self.determinants[type].append(propka.determinant.Determinant(new_determinant.group, - new_determinant.value)) - - return + self.determinants[type_].append(Determinant(new_determinant.group, + new_determinant.value)) def remove_determinants(self, labels): - """ removes all determinants with label in labels """ - for type in ['sidechain','backbone','coulomb']: - matches = list(filter(lambda d: d.label in labels, [d for d in self.determinants[type]])) - for m in matches: self.determinants[type].remove(m) + """Remove all determinants with specified labels. - return + Args: + labels: list of labels to remove + """ + for type_ in ['sidechain', 'backbone', 'coulomb']: + matches = list(filter(lambda d: d.label in labels, \ + [d for d in self.determinants[type_]])) + for match in matches: + self.determinants[type_].remove(match) def __truediv__(self, value): value = float(value) # divide all values self.pka_value /= value - self.Nmass /= value - self.Emass /= value - self.Nlocl /= value - self.Elocl /= value + self.num_volume /= value + self.energy_volume /= value + self.num_local /= value + self.energy_local /= value self.buried /= value # and all determinants - for type in ['sidechain','backbone','coulomb']: - for determinant in self.determinants[type]: + for type_ in ['sidechain', 'backbone', 'coulomb']: + for determinant in self.determinants[type_]: determinant.value /= value return self def clone(self): + """Create a copy of this group. + + Returns: + Copy of this group. + """ res = Group(self.atom) res.type = self.type res.residue_type = self.residue_type @@ -344,17 +340,14 @@ class Group: return res def setup(self): + """Set up a group.""" # set the charges if self.type in self.parameters.charge.keys(): self.charge = self.parameters.charge[self.type] if self.residue_type in self.parameters.ions.keys(): self.charge = self.parameters.ions[self.residue_type] - - #info('ION setup',self,self.residue_type, self.charge) - # find the center and the interaction atoms self.setup_atoms() - # set the model pka value self.titratable = False if self.residue_type in self.parameters.model_pkas.keys(): @@ -364,783 +357,841 @@ class Group: 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] - self.model_pka_set = True - if self.model_pka_set and not self.atom.cysteine_bridge: self.titratable = True self.exclude_cys_from_results = False - return - def setup_atoms(self): - # This method is overwritten for some types of groups + """Set up atoms in group. + + This method is overwritten for some types of groups + """ # set the center at the position of the main atom self.set_center([self.atom]) # set the main atom as interaction atom self.set_interaction_atoms([self.atom], [self.atom]) - return - def set_interaction_atoms(self, interaction_atoms_for_acids, interaction_atoms_for_bases): - [a.set_group_type(self.type) for a in interaction_atoms_for_acids+interaction_atoms_for_bases] + def set_interaction_atoms(self, interaction_atoms_for_acids, + interaction_atoms_for_bases): + """Set interacting atoms and group types. + Args: + interaction_atoms_for_acids: list of atoms for acid interactions + interaction_atoms_for_base: list of atoms for base interactions + """ + for atom in interaction_atoms_for_acids + interaction_atoms_for_bases: + atom.set_group_type(self.type) self.interaction_atoms_for_acids = interaction_atoms_for_acids self.interaction_atoms_for_bases = interaction_atoms_for_bases - # check if all atoms have been identified ok = True - for [expect, found, t] in [[expected_atoms_acid_interactions, self.interaction_atoms_for_acids, 'acid'], - [expected_atoms_base_interactions, self.interaction_atoms_for_bases, 'base']]: + for [expect, found, _] in [[EXPECTED_ATOMS_ACID_INTERACTIONS, + self.interaction_atoms_for_acids, 'acid'], + [EXPECTED_ATOMS_BASE_INTERACTIONS, + self.interaction_atoms_for_bases, 'base']]: if self.type in expect.keys(): - for e in expect[self.type].keys(): - if len([a for a in found if a.element==e]) != expect[self.type][e]: + for elem in expect[self.type].keys(): + if len([a for a in found if a.element == elem]) \ + != expect[self.type][elem]: ok = False - if not ok: - warning('Missing atoms or failed protonation for %s (%s) -- please check the structure' % (self.label, self.type)) + str_ = 'Missing atoms or failed protonation for ' + str_ += '%s (%s) -- please check the structure' % (self.label, + self.type) + warning(str_) warning('%s' % self) - Na = sum([expected_atoms_acid_interactions[self.type][e] for e in expected_atoms_acid_interactions[self.type].keys()]) - Nb = sum([expected_atoms_base_interactions[self.type][e] for e in expected_atoms_base_interactions[self.type].keys()]) - - warning('Expected %d interaction atoms for acids, found:' % Na) + num_acid = sum([EXPECTED_ATOMS_ACID_INTERACTIONS[self.type][e] \ + for e in EXPECTED_ATOMS_ACID_INTERACTIONS[self.type].keys()]) + num_base = sum([EXPECTED_ATOMS_BASE_INTERACTIONS[self.type][e] \ + for e in EXPECTED_ATOMS_BASE_INTERACTIONS[self.type].keys()]) + warning('Expected %d interaction atoms for acids, found:' % num_acid) for i in range(len(self.interaction_atoms_for_acids)): - warning(' %s' % self.interaction_atoms_for_acids[i]) - - warning('Expected %d interaction atoms for bases, found:' % Nb) + warning(' %s' % self.interaction_atoms_for_acids[i]) + warning('Expected %d interaction atoms for bases, found:' % num_base) for i in range(len(self.interaction_atoms_for_bases)): - warning(' %s' % self.interaction_atoms_for_bases[i]) - - - #return - - return + warning(' %s' % self.interaction_atoms_for_bases[i]) def get_interaction_atoms(self, interacting_group): + """Get atoms involved in interaction with other group. + + Args: + interacting_group: other group + Returns: + list of atoms + """ if interacting_group.residue_type in self.parameters.base_list: return self.interaction_atoms_for_bases else: - return self.interaction_atoms_for_acids #default is acid interaction atoms - cf. 3.0 + # default is acid interaction atoms - cf. 3.0 + return self.interaction_atoms_for_acids def set_center(self, atoms): + """Set center of group based on atoms. + + Args: + atoms: list of atoms + """ if not atoms: raise ValueError("At least one atom must be specified") - # reset center - self.x = 0.0; self.y = 0.0; self.z = 0.0 - - # find the average positon of atoms + self.x = 0.0 + self.y = 0.0 + self.z = 0.0 + # find the average position of atoms for atom in atoms: - self.x += atom.x; self.y += atom.y; self.z += atom.z - + self.x += atom.x + self.y += atom.y + self.z += atom.z self.x /= float(len(atoms)) self.y /= float(len(atoms)) self.z /= float(len(atoms)) - return + def get_determinant_string(self, remove_penalised_group=False): + """Create a string to identify this determinant. - def getDeterminantString(self, remove_penalised_group=False): + Args: + remove_penalised_group: Boolean flag to remove penalized groups + Returns: + string + """ if self.coupled_titrating_group and remove_penalised_group: return '' - number_of_sidechain = len(self.determinants['sidechain']) - number_of_backbone = len(self.determinants['backbone']) - number_of_coulomb = len(self.determinants['coulomb']) - - number_of_lines = max(1, number_of_sidechain, number_of_backbone, number_of_coulomb) - str = "" + number_of_backbone = len(self.determinants['backbone']) + number_of_coulomb = len(self.determinants['coulomb']) + number_of_lines = max(1, number_of_sidechain, number_of_backbone, + number_of_coulomb) + str_ = "" for line_number in range(number_of_lines): - str += "%s" % (self.label) + str_ += "%s" % (self.label) if line_number == 0: - str += " %6.2lf" %(self.pka_value) - if len(self.non_covalently_coupled_groups)>0: - str+='*' + str_ += " %6.2lf" %(self.pka_value) + if len(self.non_covalently_coupled_groups) > 0: + str_ += '*' else: - str+=' ' - - # if self.atom.cysteine_bridge: - # str += " BONDED " - # else: - str += " %4d%2s " % (int(100.0*self.buried), "%") - - str += " %6.2lf %4d" % (self.Emass, self.Nmass) - str += " %6.2lf %4d" % (self.Elocl, self.Nlocl) + str_ += ' ' + str_ += " %4d%2s " % (int(100.0*self.buried), "%") + str_ += " %6.2lf %4d" % (self.energy_volume, self.num_volume) + str_ += " %6.2lf %4d" % (self.energy_local, self.num_local) else: - str += "%40s" % (" ") - + str_ += "%40s" % (" ") # add the determinants - for type in ['sidechain','backbone','coulomb']: - str += self.get_determinant_for_string(type,line_number) - + for type_ in ['sidechain', 'backbone', 'coulomb']: + str_ += self.get_determinant_for_string(type_, line_number) # adding end-of-line - str += "\n" + str_ += "\n" + str_ += "\n" + return str_ - str += "\n" + def get_determinant_for_string(self, type_, number): + """Return a string describing determinant. - return str - - def get_determinant_for_string(self, type, number): - if number >= len(self.determinants[type]): + Args: + type_: determinant type + number: determinant index number + Returns: + string + """ + if number >= len(self.determinants[type_]): empty_determinant = "%s%4d%2s" % ("XXX", 0, "X") return "%8.2lf %s" % (0.0, empty_determinant) else: - determinant = self.determinants[type][number] + determinant = self.determinants[type_][number] return "%8.2lf %s" % (determinant.value, determinant.label) - def calculate_total_pka(self): + """Calculate total pKa based on determinants associated with this + group.""" # if this is a cysteine involved in a di-sulphide bond if self.atom.cysteine_bridge: self.pka_value = 99.99 return - - - self.pka_value = self.model_pka + self.Emass + self.Elocl - + self.pka_value = self.model_pka + self.energy_volume + self.energy_local for determinant_type in ['sidechain', 'backbone', 'coulomb']: for determinant in self.determinants[determinant_type]: self.pka_value += determinant.value - return - - - - def calculate_intrinsic_pka(self): + """Calculate the intrinsic pKa values from the desolvation + determinants, back-bone hydrogen bonds, and side-chain hydrogen bonds + to non-titratable residues. """ - Calculates the intrinsic pKa values from the desolvation determinants, back-bone hydrogen bonds, - and side-chain hydrogen bond to non-titratable residues - """ - back_bone = 0.0 + back_bone = 0.0 for determinant in self.determinants['backbone']: value = determinant.value back_bone += value - side_chain = 0.0 for determinant in self.determinants['sidechain']: - if determinant.label[0:3] not in ['ASP','GLU','LYS','ARG','HIS','CYS','TYR','C- ','N+ ']: + if determinant.label[0:3] not in ['ASP', 'GLU', 'LYS', 'ARG', + 'HIS', 'CYS', 'TYR', 'C- ', + 'N+ ']: value = determinant.value side_chain += value + self.intrinsic_pka = self.model_pka + self.energy_volume \ + + self.energy_local + back_bone + side_chain - self.intrinsic_pKa = self.model_pka + self.Emass + self.Elocl + back_bone + side_chain + def get_summary_string(self, remove_penalised_group=False): + """Create summary string for this group. - return - - - - - - def getSummaryString(self, remove_penalised_group=False): + Args: + remove_penalised_group: Boolean to ignore penalized groups + Returns: + string + """ if self.coupled_titrating_group and remove_penalised_group: return '' - ligand_type = '' if self.atom.type == 'hetatm': ligand_type = self.type - penalty = '' if self.coupled_titrating_group: - penalty = ' NB: Discarded due to coupling with %s'%self.coupled_titrating_group.label - - str = " %9s %8.2lf %10.2lf %18s %s\n" % (self.label, - self.pka_value, - self.model_pka,ligand_type, - penalty) - - return str + penalty = ' NB: Discarded due to coupling with %s' \ + % self.coupled_titrating_group.label + str_ = " %9s %8.2lf %10.2lf %18s %s\n" % (self.label, + self.pka_value, + self.model_pka, + ligand_type, + penalty) + return str_ def __str__(self): - return 'Group (%s) for %s'%(self.type,self.atom) - - - - # - # Energy-related methods - # + return 'Group (%s) for %s' % (self.type, self.atom) def calculate_folding_energy(self, parameters, ph=None, reference=None): + """Return the electrostatic energy of this residue at specified pH. + + Args: + parameters: parameters for energy calculation + ph: pH value for calculation + reference: reference state for calculation + Returns: + float describing energy """ - returning the electrostatic energy of this residue at pH 'pH' - """ - if ph == None: - pH = parameters.pH - if reference == None: + if ph is None: + ph = parameters.pH + if reference is None: reference = parameters.reference - # If not titratable, the contribution is zero - if not self.titratable: return 0.00 - - # calculating the ddG(neutral --> low-pH) contribution - ddG_neutral = 0.00 + # calculating the ddg(neutral --> low-pH) contribution + ddg_neutral = 0.00 if reference == 'neutral' and self.charge > 0.00: pka_prime = self.pka_value for determinant in self.determinants['coulomb']: if determinant.value > 0.00: pka_prime -= determinant.value - ddG_neutral = -1.36*(pka_prime - self.model_pka) - - # calculating the ddG(low-pH --> pH) contribution + ddg_neutral = UNK_PKA_SCALING*(pka_prime - self.model_pka) + # calculating the ddg(low-pH --> pH) contribution # folded - x = ph - self.pka_value - y = 10**x - Q_pro = math.log10(1+y) - + dpka = ph - self.pka_value + conc_ratio = 10**dpka + q_pro = math.log10(1+conc_ratio) # unfolded - x = ph - self.model_pka - y = 10**x - Q_mod = math.log10(1+y) + dpka = ph - self.model_pka + conc_ratio = 10**dpka + q_mod = math.log10(1+conc_ratio) + ddg_low = UNK_PKA_SCALING*(q_pro - q_mod) + ddg = ddg_neutral + ddg_low + return ddg - ddG_low = -1.36*(Q_pro - Q_mod) - ddG = ddG_neutral + ddG_low - - return ddG - - def calculate_charge(self, parmaeters, ph=7.0, state='folded'): + def calculate_charge(self, _, ph=7.0, state='folded'): + """Calculate the charge of the specified state at the specified pH. + Args: + _: parameters for calculation + ph: pH value + state: "folded" or "unfolded" + Returns: + float with charge + """ if state == "unfolded": - x = self.charge * (self.model_pka - ph) + q_dpka = self.charge * (self.model_pka - ph) else: - x = self.charge * (self.pka_value - ph) - - y = 10**x - charge = self.charge*(y/(1.0+y)) - + q_dpka = self.charge * (self.pka_value - ph) + conc_ratio = 10**q_dpka + charge = self.charge*(conc_ratio/(1.0+conc_ratio)) return charge def use_in_calculations(self): - """ - Whether this group should be included in the results report. If - --titrate_only option is specified, only residues that are titratable - and are in that list are included; otherwise all titratable residues - and CYS residues are included. + """Indicate whether group should be included in results report. + + If --titrate_only option is specified, only residues that are + titratable and are in that list are included; otherwise all titratable + residues and CYS residues are included. """ return self.titratable or (self.residue_type == 'CYS' and \ not self.exclude_cys_from_results) -class COO_group(Group): +class COOGroup(Group): + """Carboxyl group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'COO' def setup_atoms(self): + """Set up group.""" # Identify the two caroxyl oxygen atoms the_oxygens = self.atom.get_bonded_elements('O') - # set the center using the two oxygen carboxyl atoms (if present) if the_oxygens: self.set_center(the_oxygens) else: self.set_center([self.atom]) - # FIXME perhaps it would be better to ignore this group completely + # TODO - perhaps it would be better to ignore this group completely # if the oxygen is missing from this residue? - self.set_interaction_atoms(the_oxygens, the_oxygens) - return -class HIS_group(Group): +class HISGroup(Group): + """Histidine group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'HIS' def setup_atoms(self): + """Set up atoms in group.""" # Find the atoms in the histidine ring ring_atoms = propka.ligand.is_ring_member(self.atom) if len(ring_atoms) != 5: warning('His group does not seem to contain a ring', self) - # protonate ring - for r in ring_atoms: - my_protonator.protonate_atom(r) - + for ring_atom in ring_atoms: + PROTONATOR.protonate_atom(ring_atom) # set the center using the ring atoms if ring_atoms: self.set_center(ring_atoms) else: # Missing side-chain atoms self.set_center([self.atom]) - # FIXME perhaps it would be better to ignore this group completely? - + # TODO - perhaps it would be better to ignore this group completely? # find the hydrogens on the ring-nitrogens hydrogens = [] nitrogens = [ra for ra in ring_atoms if ra.element == 'N'] - for nitrogen in nitrogens: hydrogens.extend(nitrogen.get_bonded_elements('H')) - self.set_interaction_atoms(hydrogens+nitrogens, nitrogens) - return +class CYSGroup(Group): + """Cysteine group.""" -class CYS_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'CYS' -class TYR_group(Group): +class TYRGroup(Group): + """Tyrosine group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'TYR' -class LYS_group(Group): +class LYSGroup(Group): + """Lysine group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'LYS' -class ARG_group(Group): +class ARGGroup(Group): + """Arginine group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'ARG' def setup_atoms(self): + """Set up group.""" # set the center at the position of the main atom self.set_center([self.atom]) - # find all the hydrogens on the nitrogen atoms nitrogens = self.atom.get_bonded_elements('N') - for n in nitrogens: - my_protonator.protonate_atom(n) - + for nitrogen in nitrogens: + PROTONATOR.protonate_atom(nitrogen) hydrogens = [] for nitrogen in nitrogens: hydrogens.extend(nitrogen.get_bonded_elements('H')) self.set_interaction_atoms(nitrogens+hydrogens, nitrogens) - return +class ROHGroup(Group): + """Alcohol group.""" -class ROH_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'ROH' -class SER_group(Group): + +class SERGroup(Group): + """Serine group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'SER' -class AMD_group(Group): + +class AMDGroup(Group): + """Amide group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'AMD' def setup_atoms(self): + """Setup group.""" # Identify the oxygen and nitrogen amide atoms the_oxygen = self.atom.get_bonded_elements('O') the_nitrogen = self.atom.get_bonded_elements('N') - # add protons to the nitrogen - my_protonator.protonate_atom(the_nitrogen[0]) + PROTONATOR.protonate_atom(the_nitrogen[0]) the_hydrogens = the_nitrogen[0].get_bonded_elements('H') - # set the center using the oxygen and nitrogen amide atoms self.set_center(the_oxygen+the_nitrogen) - - self.set_interaction_atoms(the_nitrogen+the_hydrogens,the_oxygen) - - return + self.set_interaction_atoms(the_nitrogen + the_hydrogens, the_oxygen) -class TRP_group(Group): +class TRPGroup(Group): + """Tryptophan group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'TRP' def setup_atoms(self): + """Set up atoms in group.""" # set the center at the position of the main atom self.set_center([self.atom]) - # find the hydrogen on the nitrogen atom - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogen = self.atom.get_bonded_elements('H') self.set_interaction_atoms(the_hydrogen+[self.atom], [self.atom]) - return -class Nterm_group(Group): +class NtermGroup(Group): + """N-terminus group.""" def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'N+' -class Cterm_group(Group): +class CtermGroup(Group): + """C-terminus group.""" def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'COO' # this is to deal with the COO-C- parameter unification. def setup_atoms(self): + """Set up atoms in group.""" # Identify the carbon and other oxygen carboxyl atoms the_carbons = self.atom.get_bonded_elements('C') if not the_carbons: self.set_center([self.atom]) - # FIXME perhaps it would be better to ignore this group completely + # TODO - perhaps it would be better to ignore this group completely # if the carbon is missing from this residue? else: the_other_oxygen = the_carbons[0].get_bonded_elements('O') the_other_oxygen.remove(self.atom) - # set the center and interaction atoms the_oxygens = [self.atom]+ the_other_oxygen self.set_center(the_oxygens) self.set_interaction_atoms(the_oxygens, the_oxygens) - return +class BBNGroup(Group): + """Backbone nitrogen group.""" - -class BBN_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'BBN' self.residue_type = 'BBN' - def setup_atoms(self): + """Set up atoms in group.""" # Identify the hydrogen - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogen = self.atom.get_bonded_elements('H') - # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogen+[self.atom], the_hydrogen+[self.atom]) - return -class BBC_group(Group): +class BBCGroup(Group): + """Backbone carbon group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'BBC' self.residue_type = 'BBC' def setup_atoms(self): + """Set up atoms in group.""" # Identify the oxygen the_oxygen = self.atom.get_bonded_elements('O') - # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_oxygen, the_oxygen) - return -class NAR_group(Group): +class NARGroup(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'NAR' self.residue_type = 'NAR' info('Found NAR group:', atom) - return - def setup_atoms(self): + """Set up atoms in group.""" # Identify the hydrogen - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogen = self.atom.get_bonded_elements('H') - # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogen+[self.atom], the_hydrogen+[self.atom]) - return +class NAMGroup(Group): + """Unknown group. + TODO - identify this group. + """ - -class NAM_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'NAM' self.residue_type = 'NAM' info('Found NAM group:', atom) - return - def setup_atoms(self): + """Set up atoms in this group.""" # Identify the hydrogen - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogen = self.atom.get_bonded_elements('H') - # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogen+[self.atom], the_hydrogen+[self.atom]) - return +class FGroup(Group): + """Fluoride group.""" -class F_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'F' self.residue_type = 'F' info('Found F group:', atom) - return -class Cl_group(Group): + +class ClGroup(Group): + """Chloride group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'Cl' self.residue_type = 'Cl' info('Found Cl group:', atom) - return -class OH_group(Group): + +class OHGroup(Group): + """Hydroxide group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'OH' self.residue_type = 'OH' info('Found OH group:', atom) - return - def setup_atoms(self): + """Set up atoms in this group.""" # Identify the hydrogen - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogen = self.atom.get_bonded_elements('H') - # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogen+[self.atom], the_hydrogen+[self.atom]) - return -class OP_group(Group): + +class OPGroup(Group): + """Phosphate group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'OP' self.residue_type = 'OP' info('Found OP group:', atom) - return - def setup_atoms(self): + """Set up atoms in this group.""" # Identify the hydrogen - my_protonator.protonate_atom(self.atom) - #the_hydrogen = self.atom.get_bonded_elements('H') - + PROTONATOR.protonate_atom(self.atom) # set the center using the oxygen self.set_center([self.atom]) #self.set_interaction_atoms(the_hydrogen+[self.atom], the_hydrogen+[self.atom]) self.set_interaction_atoms([self.atom], [self.atom]) - return -class O3_group(Group): +class O3Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'O3' self.residue_type = 'O3' info('Found O3 group:', atom) - return -class O2_group(Group): +class O2Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'O2' self.residue_type = 'O2' info('Found O2 group:', atom) - return -class SH_group(Group): + +class SHGroup(Group): + """Sulfhydryl group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'SH' self.residue_type = 'SH' info('Found SH group:', atom) - return -class CG_group(Group): - """Guadinium group""" +class CGGroup(Group): + """Guadinium group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'CG' self.residue_type = 'CG' info('Found CG group:', atom) - return def setup_atoms(self): + """Set up atoms in this group.""" # Identify the nitrogens the_nitrogens = self.atom.get_bonded_elements('N') - # set the center using the nitrogen self.set_center([self.atom]) - the_hydrogens = [] - for n in the_nitrogens: - my_protonator.protonate_atom(n) - the_hydrogens += n.get_bonded_elements('H') + for nitrogen in the_nitrogens: + PROTONATOR.protonate_atom(nitrogen) + the_hydrogens += nitrogen.get_bonded_elements('H') self.set_interaction_atoms(the_hydrogens+the_nitrogens, the_nitrogens) - return -class C2N_group(Group): - """Amidinium group""" +class C2NGroup(Group): + """Amidinium group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'C2N' self.residue_type = 'C2N' info('Found C2N group:', atom) - return def setup_atoms(self): + """Set up atoms in this group.""" # Identify the nitrogens the_nitrogens = self.atom.get_bonded_elements('N') - the_nitrogens = [n for n in the_nitrogens if len(n.get_bonded_heavy_atoms())==1] - + the_nitrogens = [n for n in the_nitrogens \ + if len(n.get_bonded_heavy_atoms()) == 1] # set the center using the nitrogen self.set_center([self.atom]) - the_hydrogens = [] - for n in the_nitrogens: - my_protonator.protonate_atom(n) - the_hydrogens += n.get_bonded_elements('H') - + for nitrogen in the_nitrogens: + PROTONATOR.protonate_atom(nitrogen) + the_hydrogens += nitrogen.get_bonded_elements('H') self.set_interaction_atoms(the_hydrogens+the_nitrogens, the_nitrogens) - return -class OCO_group(Group): + +class OCOGroup(Group): + """Carboxyl group.""" + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'OCO' self.residue_type = 'OCO' info('Found OCO group:', atom) - return def setup_atoms(self): - # Identify the two caroxyl oxygen atoms + """Set up atoms in group.""" + # Identify the two carboxyl oxygen atoms the_oxygens = self.atom.get_bonded_elements('O') - # set the center using the two oxygen carboxyl atoms self.set_center(the_oxygens) self.set_interaction_atoms(the_oxygens, the_oxygens) - return +class N30Group(Group): + """Unknown group. + + TODO - identify this group. + """ -class N30_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'N30' self.residue_type = 'N30' info('Found N30 group:', atom) - return def setup_atoms(self): + """Set up atoms in this group.""" # Identify the nitrogens - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogens = self.atom.get_bonded_elements('H') # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogens+[self.atom], [self.atom]) - return -class N31_group(Group): + +class N31Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'N31' self.residue_type = 'N31' info('Found N31 group:', atom) - return def setup_atoms(self): + """Set up atoms in this group.""" # Identify the nitrogens - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogens = self.atom.get_bonded_elements('H') # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogens+[self.atom], [self.atom]) - return -class N32_group(Group): + +class N32Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'N32' self.residue_type = 'N32' info('Found N32 group:', atom) - return def setup_atoms(self): + """Set up atoms in this group.""" # Identify the nitrogens - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogens = self.atom.get_bonded_elements('H') # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogens+[self.atom], [self.atom]) - return -class N33_group(Group): + +class N33Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'N33' self.residue_type = 'N33' info('Found N33 group:', atom) - return def setup_atoms(self): + """Set up atoms in this group.""" # Identify the nitrogens - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogens = self.atom.get_bonded_elements('H') # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogens+[self.atom], [self.atom]) - return -class NP1_group(Group): + +class NP1Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'NP1' self.residue_type = 'NP1' info('Found NP1 group:', atom) - return - def setup_atoms(self): + """Set up atoms in group.""" # Identify the nitrogens - my_protonator.protonate_atom(self.atom) + PROTONATOR.protonate_atom(self.atom) the_hydrogens = self.atom.get_bonded_elements('H') # set the center using the nitrogen self.set_center([self.atom]) self.set_interaction_atoms(the_hydrogens+[self.atom], [self.atom]) - return -class N1_group(Group): + +class N1Group(Group): + """Unknown group. + + TODO - identify this group. + """ + def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'N1' self.residue_type = 'N1' info('Found N1 group:', atom) - return +class IonGroup(Group): + """Ion group.""" -class Ion_group(Group): def __init__(self, atom): - Group.__init__(self,atom) + Group.__init__(self, atom) self.type = 'ION' self.residue_type = atom.res_name.strip() info('Found ion group:', atom) - return -class non_titratable_ligand_group(Group): +class NonTitratableLigandGroup(Group): + """Non-titratable ligand group.""" + def __init__(self, atom): Group.__init__(self, atom) self.type = 'LG' self.residue_type = 'LG' -# info('Non-titratable ligand group',atom) - return -class titratable_ligand_group(Group): + +class TitratableLigandGroup(Group): + """Titratable ligand group.""" + def __init__(self, atom): Group.__init__(self, atom) # set the charge and determine type (acid or base) self.charge = atom.charge - if self.charge <0: + if self.charge < 0: self.type = 'ALG' self.residue_type = 'ALG' elif self.charge > 0: @@ -1148,8 +1199,6 @@ class titratable_ligand_group(Group): self.residue_type = 'BLG' else: raise Exception('Unable to determine type of ligand group - charge not set?') - - # check if marvin model pka has been calculated # this is not true if we are reading an input file if atom.marvin_pka: @@ -1157,192 +1206,196 @@ class titratable_ligand_group(Group): info('Titratable ligand group ', atom, self.model_pka, self.charge) self.model_pka_set = True - return - def is_group(parameters, atom): - atom.groups_extracted = 1 + """Identify whether the atom belongs to a group. + Args: + parameters: parameters for check + atom: atom to check + Returns: + group for atom or None + """ + atom.groups_extracted = 1 # check if this atom belongs to a protein group protein_group = is_protein_group(parameters, atom) - if protein_group: return protein_group - + if protein_group: + return protein_group # check if this atom belongs to a ion group ion_group = is_ion_group(parameters, atom) - if ion_group: return ion_group - + if ion_group: + return ion_group # check if this atom belongs to a ligand group if parameters.ligand_typing == 'marvin': ligand_group = is_ligand_group_by_marvin_pkas(parameters, atom) elif parameters.ligand_typing == 'sybyl': - ligand_group = is_ligand_group_by_sybyl_types(parameters, atom) + ligand_group = None elif parameters.ligand_typing == 'groups': ligand_group = is_ligand_group_by_groups(parameters, atom) else: - raise Exception('Unknown ligand typing method \'%s\''%parameters.ligand_typing) - - if ligand_group: return ligand_group - - - + errstr = 'Unknown ligand typing method \'%s\'' % parameters.ligand_typing + raise Exception(errstr) + if ligand_group: + return ligand_group return None -def is_protein_group(parameters,atom): +def is_protein_group(parameters, atom): + """Identify whether the atom belongs to a protein group. + + Args: + parameters: parameters for check + atom: atom to check + Returns: + group for atom or None + """ if atom.type != 'atom': return None - ### Check for termial groups if atom.terminal == 'N+': - return Nterm_group(atom) + return NtermGroup(atom) elif atom.terminal == 'C-': - return Cterm_group(atom) - + return CtermGroup(atom) ### Backbone if atom.type == 'atom' and atom.name == 'N': # ignore proline backbone nitrogens if atom.res_name != 'PRO': - return BBN_group(atom) + return BBNGroup(atom) if atom.type == 'atom' and atom.name == 'C': # ignore C- carboxyl if atom.count_bonded_elements('O') == 1: - return BBC_group(atom) - + return BBCGroup(atom) ### Filters for side chains based on PDB protein atom names - key = '%s-%s'%(atom.res_name, 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]) - - return None - -def is_ligand_group_by_sybyl_types(parameters, atom): - - + class_str = "%sGroup" % parameters.protein_group_mapping[key] + group_class = globals()[class_str] + return group_class(atom) return None -def is_ligand_group_by_groups(parameters, atom): +def is_ligand_group_by_groups(_, atom): + """Identify whether the atom belongs to a ligand group by checking groups. + + Args: + _: parameters for check + atom: atom to check + Returns: + group for atom or None + """ ### Ligand group filters if atom.type != 'hetatm': return None - - my_protonator.protonate_atom(atom) - + PROTONATOR.protonate_atom(atom) if atom.sybyl_type == 'N.ar': - if len(atom.get_bonded_heavy_atoms())==2: - return NAR_group(atom) - + if len(atom.get_bonded_heavy_atoms()) == 2: + return NARGroup(atom) if atom.sybyl_type == 'N.am': - return NAM_group(atom) - + return NAMGroup(atom) if atom.sybyl_type in ['N.3', 'N.4']: heavy_bonded = atom.get_bonded_heavy_atoms() if len(heavy_bonded) == 0: - return N30_group(atom) + return N30Group(atom) elif len(heavy_bonded) == 1: - return N31_group(atom) + return N31Group(atom) elif len(heavy_bonded) == 2: - return N32_group(atom) + return N32Group(atom) elif len(heavy_bonded) == 3: - return N33_group(atom) - + return N33Group(atom) if atom.sybyl_type == 'N.1': - return N1_group(atom) - + return N1Group(atom) if atom.sybyl_type == 'N.pl3': # make sure that this atom is not part of a guadinium or amidinium group bonded_carbons = atom.get_bonded_elements('C') if len(bonded_carbons) == 1: bonded_nitrogens = bonded_carbons[0].get_bonded_elements('N') if len(bonded_nitrogens) == 1: - return NP1_group(atom) - - + return NP1Group(atom) if atom.sybyl_type == 'C.2': # Guadinium and amidinium groups bonded_nitrogens = atom.get_bonded_elements('N') - npls = [n for n in bonded_nitrogens if (n.sybyl_type == 'N.pl3' and len(n.get_bonded_heavy_atoms())==1)] + npls = [n for n in bonded_nitrogens if (n.sybyl_type == 'N.pl3' \ + and len(n.get_bonded_heavy_atoms()) == 1)] if len(npls) == 2: - n_with_max_two_heavy_atom_bonds = [n for n in bonded_nitrogens if len(n.get_bonded_heavy_atoms())<3] + n_with_max_two_heavy_atom_bonds = [n for n in bonded_nitrogens \ + if len(n.get_bonded_heavy_atoms()) < 3] if len(n_with_max_two_heavy_atom_bonds) == 2: - return C2N_group(atom) + return C2NGroup(atom) if len(bonded_nitrogens) == 3: - return CG_group(atom) + return CGGroup(atom) # carboxyl group bonded_oxygens = atom.get_bonded_elements('O') bonded_oxygens = [b for b in bonded_oxygens if 'O.co2' in b.sybyl_type] if len(bonded_oxygens) == 2: - return OCO_group(atom) - - + return OCOGroup(atom) if atom.sybyl_type == 'F': - return F_group(atom) - + return FGroup(atom) if atom.sybyl_type == 'Cl': - return Cl_group(atom) - + return ClGroup(atom) if atom.sybyl_type == 'O.3': if len(atom.get_bonded_heavy_atoms()) == 1: # phosphate group if atom.count_bonded_elements('P') == 1: - return OP_group(atom) + return OPGroup(atom) # hydroxyl group else: - return OH_group(atom) + return OHGroup(atom) # other SP3 oxygen else: - return O3_group(atom) - + return O3Group(atom) if atom.sybyl_type == 'O.2': - return O2_group(atom) - - + return O2Group(atom) if atom.sybyl_type == 'S.3': # sulphide group if len(atom.get_bonded_heavy_atoms()) == 1: - return SH_group(atom) - # other SP3 oxygen - #else: - # return S3_group(atom) - - + return SHGroup(atom) return None def is_ligand_group_by_marvin_pkas(parameters, atom): + """Identify whether the atom belongs to a ligand group by calculating + 'Marvin pKas'. + + Args: + parameters: parameters for check + atom: atom to check + Returns: + group for atom or None + """ if atom.type != 'hetatm': return None - # calculate Marvin ligand pkas for this conformation container # if not already done + # TODO - double-check testing coverage of these functions. if not atom.conformation_container.marvin_pkas_calculated: - lpka = propka.ligand_pka_values.ligand_pka_values(parameters) + lpka = ligand_pka_values(parameters) lpka.get_marvin_pkas_for_molecular_container(atom.molecular_container, - min_ph=parameters.min_ligand_model_pka, - max_ph=parameters.max_ligand_model_pka) - - + min_pH=parameters.min_ligand_model_pka, + max_pH=parameters.max_ligand_model_pka) if atom.marvin_pka: - return titratable_ligand_group(atom) - + return TitratableLigandGroup(atom) # Special case of oxygen in carboxyl group not assigned a pka value by marvin if atom.sybyl_type == 'O.co2': atom.charge = -1.0 - other_oxygen = [o for o in atom.get_bonded_elements('C')[0].get_bonded_elements('O') if not o==atom][0] + other_oxygen = [o for o \ + in atom.get_bonded_elements('C')[0].get_bonded_elements('O') \ + if not o == atom][0] atom.marvin_pka = other_oxygen.marvin_pka - return titratable_ligand_group(atom) - - + return TitratableLigandGroup(atom) if atom.element in parameters.hydrogen_bonds.elements: - return non_titratable_ligand_group(atom) - + return NonTitratableLigandGroup(atom) return None def is_ion_group(parameters, atom): + """Identify whether the atom belongs to an ion group. + Args: + parameters: parameters for check + atom: atom to check + Returns: + group for atom or None + """ if atom.res_name.strip() in parameters.ions.keys(): - return Ion_group(atom) - + return IonGroup(atom) return None diff --git a/propka/iterative.py b/propka/iterative.py index 132d42b..e3bebf9 100644 --- a/propka/iterative.py +++ b/propka/iterative.py @@ -339,8 +339,8 @@ class Iterative: coulomb += value self.pKa_NonIterative = group.model_pka - self.pKa_NonIterative += group.Emass - self.pKa_NonIterative += group.Elocl + self.pKa_NonIterative += group.energy_volume + self.pKa_NonIterative += group.energy_local self.pKa_NonIterative += side_chain self.pKa_NonIterative += back_bone self.pKa_NonIterative += coulomb diff --git a/propka/output.py b/propka/output.py index 531af02..6e15d32 100644 --- a/propka/output.py +++ b/propka/output.py @@ -133,7 +133,7 @@ def getDeterminantSection(protein, conformation, parameters): 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) ) + str += "%s" % ( group.get_determinant_string(parameters.remove_penalised_group) ) # Add a warning in case of coupled residues if protein.conformations[conformation].non_covalently_coupled_groups and not protein.options.display_coupled_residues: @@ -151,7 +151,7 @@ def getSummarySection(protein, conformation, parameters): for residue_type in parameters.write_out_order: for group in protein.conformations[conformation].groups: if group.residue_type == residue_type: - str += "%s" % ( group.getSummaryString(parameters.remove_penalised_group) ) + str += "%s" % ( group.get_summary_string(parameters.remove_penalised_group) ) return str diff --git a/propka/parameters.py b/propka/parameters.py index 143a59b..90e47b5 100644 --- a/propka/parameters.py +++ b/propka/parameters.py @@ -33,7 +33,7 @@ parameters = ['Nmin','Nmax','desolvationSurfaceScalingFactor','desolvat 'include_H_in_interactions','coupling_max_number_of_bonds', 'min_bond_distance_for_hydrogen_bonds','coupling_penalty', 'shared_determinants','common_charge_centre','hide_penalised_group', 'remove_penalised_group', - 'max_intrinsic_pKa_diff','min_interaction_energy','max_free_energy_diff','min_swap_pka_shift', + 'max_intrinsic_pka_diff','min_interaction_energy','max_free_energy_diff','min_swap_pka_shift', 'min_pka','max_pka','sidechain_interaction'] strings = ['version','output_file_tag','ligand_typing','pH','reference'] diff --git a/propka/propka.cfg b/propka/propka.cfg index a12e38a..58b1aec 100644 --- a/propka/propka.cfg +++ b/propka/propka.cfg @@ -342,7 +342,7 @@ common_charge_centre 0 remove_penalised_group 1 # non-covalent coupling -max_intrinsic_pKa_diff 2.0 +max_intrinsic_pka_diff 2.0 min_interaction_energy 0.5 max_free_energy_diff 1.0 min_swap_pka_shift 1.0 diff --git a/propka/version.py b/propka/version.py index 861fc5b..86c4268 100644 --- a/propka/version.py +++ b/propka/version.py @@ -18,8 +18,8 @@ class version: def calculate_desolvation(self, group): return self.desolvation_model(self.parameters, group) - def calculate_pair_weight(self, Nmass1, Nmass2): - return self.weight_pair_method(self.parameters, Nmass1, Nmass2) + def calculate_pair_weight(self, num_volume1, num_volume2): + return self.weight_pair_method(self.parameters, num_volume1, num_volume2) # side chains def hydrogen_bond_interaction(self, group1, group2):