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

De-lint iterative.py.

Browse Source
This commit is contained in:
Nathan Baker
2020-05-24 12:05:48 -07:00
parent 95e132e520
commit 04b52d458c
2 changed files with 220 additions and 223 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
propka/determinants.py
View File

@@ -42,15 +42,14 @@ def set_determinants(propka_groups, version=None, options=None):
interaction_type = version.parameters.interaction_matrix.get_value(group1.type, interaction_type = version.parameters.interaction_matrix.get_value(group1.type,
group2.type) group2.type)
if interaction_type == 'I': if interaction_type == 'I':
propka.iterative.addtoDeterminantList(group1, group2, propka.iterative.add_to_determinant_list(group1, group2,
distance, distance,
iterative_interactions, iterative_interactions,
version=version) version=version)
elif interaction_type == 'N': elif interaction_type == 'N':
add_determinants(group1, group2, distance, version) add_determinants(group1, group2, distance, version)
# --- Iterative section ---# # --- Iterative section ---#
propka.iterative.add_determinants(iterative_interactions, version, propka.iterative.add_determinants(iterative_interactions, version)
options=options)
def add_determinants(group1, group2, distance, version): def add_determinants(group1, group2, distance, version):

326
propka/iterative.py
View File

@@ -1,57 +1,63 @@
"""Iterative functions for pKa calculations.
from __future__ import division These appear to mostly involve determinants.
from __future__ import print_function """
import math, time
import propka.lib as lib
from propka.determinant import Determinant from propka.determinant import Determinant
import propka.calculations from propka.lib import info, debug
from propka.lib import info, warning, debug
# Some library functions for the interative pKa determinants
def addtoDeterminantList(group1, group2, distance, iterative_interactions, version): # TODO - these are undocumented constants
UNK_MIN_VALUE = 0.005
def add_to_determinant_list(group1, group2, distance, iterative_interactions, version):
"""Add iterative determinantes to the list.
[[R1, R2], [side-chain, coulomb], [A1, A2]], ...
NOTE - the sign is determined when the interaction is added to the iterative object!
NOTE - distance < coulomb_cutoff here
Args:
group1: first group in pair
group2: second group in pair
distance: distance between groups
iterative_interactions: interaction list to modify
version: version object
""" """
Adds 'iterative determinants' to list ..., [[R1, R2], [side-chain, coulomb], [A1, A2]], ...
Note, the sign is determined when the interaction is added to the iterative object!
Note, distance < coulomb_cutoff here
"""
hbond_value = version.hydrogen_bond_interaction(group1, group2) hbond_value = version.hydrogen_bond_interaction(group1, group2)
coulomb_value = version.electrostatic_interaction(group1, group2, distance) coulomb_value = version.electrostatic_interaction(group1, group2, distance)
# adding the interaction to 'iterative_interactions' # adding the interaction to 'iterative_interactions'
if hbond_value or coulomb_value: if hbond_value or coulomb_value:
pair = [group1, group2] pair = [group1, group2]
values = [hbond_value, coulomb_value] values = [hbond_value, coulomb_value]
while None in values: while None in values:
values[values.index(None)] = 0.0 values[values.index(None)] = 0.0
annihilation = [0., 0.] annihilation = [0., 0.]
interaction = [pair, values, annihilation] interaction = [pair, values, annihilation]
iterative_interactions.append(interaction) iterative_interactions.append(interaction)
return
def add_iterative_acid_pair(object1, object2, interaction):
"""Add the Coulomb 'iterative' interaction (an acid pair).
def addIterativeAcidPair(object1, object2, interaction): The higher pKa is raised with QQ+HB
""" The lower pKa is lowered with HB
Adding the Coulomb 'iterative' interaction (an acid pair):
the higher pKa is raised with QQ+HB Args:
the lower pKa is lowered with HB object1: first object in pair
object2: second object in pair
interaction: list with [values, annihilation]
""" """
values = interaction[1] values = interaction[1]
annihilation = interaction[2] annihilation = interaction[2]
hbond_value = values[0] hbond_value = values[0]
coulomb_value = values[1] coulomb_value = values[1]
diff = coulomb_value + 2*hbond_value diff = coulomb_value + 2*hbond_value
comp1 = object1.pKa_old + annihilation[0] + diff comp1 = object1.pka_old + annihilation[0] + diff
comp2 = object2.pKa_old + annihilation[1] + diff comp2 = object2.pka_old + annihilation[1] + diff
annihilation[0] = 0. annihilation[0] = 0.0
annihilation[1] = 0. annihilation[1] = 0.0
if comp1 > comp2: if comp1 > comp2:
# side-chain # side-chain
determinant = [object2, hbond_value] determinant = [object2, hbond_value]
@@ -74,10 +80,15 @@ def addIterativeAcidPair(object1, object2, interaction):
annihilation[1] = -diff annihilation[1] = -diff
def addIterativeBasePair(object1, object2, interaction): def add_iterative_base_pair(object1, object2, interaction):
""" """Add the Coulomb 'iterative' interaction (a base pair).
Adding the Coulomb 'iterative' interaction (a base pair):
the lower pKa is lowered The lower pKa is lowered
Args:
object1: first object in pair
object2: second object in pair
interaction: list with [values, annihilation]
""" """
values = interaction[1] values = interaction[1]
annihilation = interaction[2] annihilation = interaction[2]
@@ -85,10 +96,10 @@ def addIterativeBasePair(object1, object2, interaction):
coulomb_value = values[1] coulomb_value = values[1]
diff = coulomb_value + 2*hbond_value diff = coulomb_value + 2*hbond_value
diff = -diff diff = -diff
comp1 = object1.pKa_old + annihilation[0] + diff comp1 = object1.pka_old + annihilation[0] + diff
comp2 = object2.pKa_old + annihilation[1] + diff comp2 = object2.pka_old + annihilation[1] + diff
annihilation[0] = 0. annihilation[0] = 0.0
annihilation[1] = 0. annihilation[1] = 0.0
if comp1 < comp2: if comp1 < comp2:
# side-chain # side-chain
determinant = [object2, -hbond_value] determinant = [object2, -hbond_value]
@@ -111,254 +122,241 @@ def addIterativeBasePair(object1, object2, interaction):
annihilation[1] = -diff annihilation[1] = -diff
def addIterativeIonPair(object1, object2, interaction, version): def add_iterative_ion_pair(object1, object2, interaction, version):
""" """Add the Coulomb 'iterative' interaction (an acid-base pair)
Adding the Coulomb 'iterative' interaction (an acid-base pair):
the pKa of the acid is lowered & the pKa of the base is raised the pKa of the acid is lowered & the pKa of the base is raised
Args:
object1: first object in pair
object2: second object in pair
interaction: list with [values, annihilation]
version: version object
""" """
values = interaction[1] values = interaction[1]
annihilation = interaction[2] annihilation = interaction[2]
hbond_value = values[0] hbond_value = values[0]
coulomb_value = values[1] coulomb_value = values[1]
Q1 = object1.Q q1 = object1.q
Q2 = object2.Q q2 = object2.q
comp1 = object1.pKa_old + annihilation[0] + Q1*coulomb_value comp1 = object1.pka_old + annihilation[0] + q1*coulomb_value
comp2 = object2.pKa_old + annihilation[1] + Q2*coulomb_value comp2 = object2.pka_old + annihilation[1] + q2*coulomb_value
if object1.res_name not in version.parameters.exclude_sidechain_interactions: if object1.res_name not in version.parameters.exclude_sidechain_interactions:
comp1 += Q1*hbond_value comp1 += q1*hbond_value
if object2.res_name not in version.parameters.exclude_sidechain_interactions: if object2.res_name not in version.parameters.exclude_sidechain_interactions:
comp2 += Q2*hbond_value comp2 += q2*hbond_value
if q1 == -1.0 and comp1 < comp2:
if Q1 == -1.0 and comp1 < comp2:
add_term = True # pKa(acid) < pKa(base) add_term = True # pKa(acid) < pKa(base)
elif Q1 == 1.0 and comp1 > comp2: elif q1 == 1.0 and comp1 > comp2:
add_term = True # pKa(base) > pKa(acid) add_term = True # pKa(base) > pKa(acid)
else: else:
add_term = False add_term = False
annihilation[0] = 0.00 annihilation[0] = 0.00
annihilation[1] = 0.00 annihilation[1] = 0.00
if add_term:
if add_term == True:
# Coulomb # Coulomb
if coulomb_value > 0.005: if coulomb_value > UNK_MIN_VALUE:
# residue1 # residue1
interaction = [object2, Q1*coulomb_value] interaction = [object2, q1*coulomb_value]
annihilation[0] += -Q1*coulomb_value annihilation[0] += -q1*coulomb_value
object1.determinants['coulomb'].append(interaction) object1.determinants['coulomb'].append(interaction)
# residue2 # residue2
interaction = [object1, Q2*coulomb_value] interaction = [object1, q2*coulomb_value]
annihilation[1] += -Q2*coulomb_value annihilation[1] += -q2*coulomb_value
object2.determinants['coulomb'].append(interaction) object2.determinants['coulomb'].append(interaction)
# Side-chain # Side-chain
if hbond_value > 0.005: if hbond_value > UNK_MIN_VALUE:
# residue1 # residue1
if object1.res_name not in version.parameters.exclude_sidechain_interactions: if object1.res_name not in version.parameters.exclude_sidechain_interactions:
interaction = [object2, Q1*hbond_value] interaction = [object2, q1*hbond_value]
annihilation[0] += -Q1*hbond_value annihilation[0] += -q1*hbond_value
object1.determinants['sidechain'].append(interaction) object1.determinants['sidechain'].append(interaction)
# residue2 # residue2
if object2.res_name not in version.parameters.exclude_sidechain_interactions: if object2.res_name not in version.parameters.exclude_sidechain_interactions:
interaction = [object1, Q2*hbond_value] interaction = [object1, q2*hbond_value]
annihilation[1] += -Q2*hbond_value annihilation[1] += -q2*hbond_value
object2.determinants['sidechain'].append(interaction) object2.determinants['sidechain'].append(interaction)
def add_determinants(iterative_interactions, version, options=None): def add_determinants(iterative_interactions, version, _=None):
""" """Add determinants iteratively.
The iterative pKa scheme. Later it is all added in 'calculateTotalPKA' The iterative pKa scheme. Later it is all added in 'calculateTotalPKA'
Args:
iterative_interactions: list of iterative interactions
version: version object
_: options object
""" """
# --- setup --- # --- setup ---
iteratives = [] iteratives = []
done_group = [] done_group = []
# creating iterative objects with references to their real group counterparts # creating iterative objects with references to their real group counterparts
for interaction in iterative_interactions: for interaction in iterative_interactions:
pair = interaction[0] pair = interaction[0]
for group in pair: for group in pair:
if group in done_group: if group in done_group:
#print "done already" # do nothing - already have an iterative object for this group
""" do nothing - already have an iterative object for this group """ pass
else: else:
newIterative = Iterative(group) new_iterative = Iterative(group)
iteratives.append(newIterative) iteratives.append(new_iterative)
done_group.append(group) done_group.append(group)
# Initialize iterative scheme # Initialize iterative scheme
debug("\n --- pKa iterations (%d groups, %d interactions) ---" % debug("\n --- pKa iterations (%d groups, %d interactions) ---" %
(len(iteratives), len(iterative_interactions))) (len(iteratives), len(iterative_interactions)))
converged = False converged = False
iteration = 0 iteration = 0
# set non-iterative pka values as first step # set non-iterative pka values as first step
for itres in iteratives: for iter_ in iteratives:
itres.pKa_iter.append(itres.pKa_NonIterative) iter_.pka_iter.append(iter_.pka_noniterative)
# --- starting pKa iterations --- # --- starting pKa iterations ---
while converged == False: while not converged:
# initialize pka_new
# initialize pKa_new
iteration += 1 iteration += 1
for itres in iteratives: for itres in iteratives:
itres.determinants = {'sidechain':[],'backbone':[],'coulomb':[]} itres.determinants = {'sidechain': [], 'backbone': [],
itres.pKa_new = itres.pKa_NonIterative 'coulomb': []}
itres.pka_new = itres.pka_noniterative
# Adding interactions to temporary determinant container # Adding interactions to temporary determinant container
for interaction in iterative_interactions: for interaction in iterative_interactions:
pair = interaction[0] pair = interaction[0]
values = interaction[1] object1, object2 = find_iterative(pair, iteratives)
annihilation = interaction[2] q1 = object1.q
#print "len(interaction) = %d" % (len(interaction)) q2 = object2.q
object1, object2 = findIterative(pair, iteratives) if q1 < 0.0 and q2 < 0.0:
Q1 = object1.Q # both are acids
Q2 = object2.Q add_iterative_acid_pair(object1, object2, interaction)
if Q1 < 0.0 and Q2 < 0.0: elif q1 > 0.0 and q2 > 0.0:
""" both are acids """ # both are bases
addIterativeAcidPair(object1, object2, interaction) add_iterative_base_pair(object1, object2, interaction)
elif Q1 > 0.0 and Q2 > 0.0:
""" both are bases """
addIterativeBasePair(object1, object2, interaction)
else: else:
""" one of each """ # one of each
addIterativeIonPair(object1, object2, interaction, version) add_iterative_ion_pair(object1, object2, interaction, version)
# Calculating pka_new values
# Calculating pKa_new values
for itres in iteratives: for itres in iteratives:
for type in ['sidechain','backbone','coulomb']: for type_ in ['sidechain', 'backbone', 'coulomb']:
for determinant in itres.determinants[type]: for determinant in itres.determinants[type_]:
itres.pKa_new += determinant[1] itres.pka_new += determinant[1]
# Check convergence # Check convergence
converged = True converged = True
for itres in iteratives: for itres in iteratives:
if itres.pKa_new == itres.pKa_old: if itres.pka_new == itres.pka_old:
itres.converged = True itres.converged = True
else: else:
itres.converged = False itres.converged = False
converged = False converged = False
# reset pKa_old & storing pKa_new in pKa_iter # reset pka_old & storing pka_new in pka_iter
for itres in iteratives: for itres in iteratives:
itres.pKa_old = itres.pKa_new itres.pka_old = itres.pka_new
itres.pKa_iter.append(itres.pKa_new) itres.pka_iter.append(itres.pka_new)
if iteration == 10: if iteration == 10:
info("did not converge in %d iterations" % (iteration)) info("did not converge in %d iterations" % (iteration))
break break
# --- Iterations finished ---
# printing pKa iterations # printing pKa iterations
# formerly was conditioned on if options.verbosity >= 2 - now unnecessary # formerly was conditioned on if options.verbosity >= 2 - now unnecessary
str = "%12s" % (" ") str_ = "%12s" % (" ")
for index in range(0, iteration+1 ): for index in range(iteration+1):
str += "%8d" % (index) str_ += "%8d" % (index)
debug(str) debug(str_)
for itres in iteratives: for itres in iteratives:
str = "%s " % (itres.label) str_ = "%s " % (itres.label)
for pKa in itres.pKa_iter: for pka in itres.pka_iter:
str += "%8.2lf" % (pKa) str_ += "%8.2lf" % (pka)
if itres.converged == False: if not itres.converged:
str += " *" str_ += " *"
debug(str) debug(str_)
# creating real determinants and adding them to group object # creating real determinants and adding them to group object
for itres in iteratives: for itres in iteratives:
for type in ['sidechain','backbone','coulomb']: for type_ in ['sidechain', 'backbone', 'coulomb']:
for interaction in itres.determinants[type]: for interaction in itres.determinants[type_]:
#info('done',itres.group.label,interaction[0],interaction[1]) #info('done',itres.group.label,interaction[0],interaction[1])
value = interaction[1] value = interaction[1]
if value > 0.005 or value < -0.005: if value > UNK_MIN_VALUE or value < -UNK_MIN_VALUE:
g = interaction[0] group = interaction[0]
newDeterminant = Determinant(g, value) new_det = Determinant(group, value)
itres.group.determinants[type].append(newDeterminant) itres.group.determinants[type_].append(new_det)
def find_iterative(pair, iteratives):
"""Find the 'iteratives' that correspond to the groups in 'pair'.
def findIterative(pair, iteratives): Args:
""" pair: groups to match
Function to find the two 'iteratives' that corresponds to the groups in 'pair' iteratives: list of iteratives to search
Returns:
1. first matched iterative
2. second matched iterative
""" """
for iterative in iteratives: for iterative in iteratives:
if iterative.group == pair[0]: if iterative.group == pair[0]:
iterative0 = iterative iterative0 = iterative
elif iterative.group == pair[1]: elif iterative.group == pair[1]:
iterative1 = iterative iterative1 = iterative
return iterative0, iterative1 return iterative0, iterative1
class Iterative: class Iterative:
""" """Iterative class - pKa values and references of iterative groups.
Iterative class - pKa values and references of iterative groups
Note, this class has a fake determinant list, true determinants are NOTE - this class has a fake determinant list, true determinants are made
made after the iterations are finished. after the iterations are finished.
""" """
def __init__(self, group): def __init__(self, group):
""" """Initialize object with group.
Contructer of the iterative object
"""
#print "creating 'iterative object' for %s" % (group.label)
Args:
group: group to use for initialization.
"""
self.label = group.label self.label = group.label
self.atom = group.atom self.atom = group.atom
self.res_name = group.residue_type self.res_name = group.residue_type
self.Q = group.charge self.q = group.charge
self.pKa_old = None self.pka_old = None
self.pKa_new = None self.pka_new = None
self.pKa_iter = [] self.pka_iter = []
self.pKa_NonIterative = 0.00 self.pka_noniterative = 0.00
self.determinants = {'sidechain': [], 'backbone': [], 'coulomb': []} self.determinants = {'sidechain': [], 'backbone': [], 'coulomb': []}
self.group = group self.group = group
self.converged = True self.converged = True
# Calculate the Non-Iterative part of pKa from the group object # Calculate the Non-Iterative part of pKa from the group object
# Side chain # Side chain
side_chain = 0.00 side_chain = 0.00
for determinant in group.determinants['sidechain']: for determinant in group.determinants['sidechain']:
value = determinant.value value = determinant.value
side_chain += value side_chain += value
# Back bone # Back bone
back_bone = 0.00 back_bone = 0.00
for determinant in group.determinants['backbone']: for determinant in group.determinants['backbone']:
value = determinant.value value = determinant.value
back_bone += value back_bone += value
# Coulomb # Coulomb
coulomb = 0.00 coulomb = 0.00
for determinant in group.determinants['coulomb']: for determinant in group.determinants['coulomb']:
value = determinant.value value = determinant.value
coulomb += value coulomb += value
self.pka_noniterative = group.model_pka
self.pKa_NonIterative = group.model_pka self.pka_noniterative += group.energy_volume
self.pKa_NonIterative += group.energy_volume self.pka_noniterative += group.energy_local
self.pKa_NonIterative += group.energy_local self.pka_noniterative += side_chain
self.pKa_NonIterative += side_chain self.pka_noniterative += back_bone
self.pKa_NonIterative += back_bone self.pka_noniterative += coulomb
self.pKa_NonIterative += coulomb self.pka_old = self.pka_noniterative
self.pKa_old = self.pKa_NonIterative
def __eq__(self, other): def __eq__(self, other):
""" """Needed to use objects in sets."""
Check if two groups should be considered identical
"""
if self.atom.type == 'atom': if self.atom.type == 'atom':
# In case of protein atoms we trust the labels # In case of protein atoms we trust the labels
return self.label == other.label return self.label == other.label
else: else:
# For heterogene atoms we also need to check the residue number # For heterogene atoms we also need to check the residue number
return self.label==other.label and self.atom.res_num == other.atom.res_num return self.label == other.label \
and self.atom.res_num == other.atom.res_num
def __hash__(self): def __hash__(self):
""" Needed together with __eq__ - otherwise we can't make sets of groups """ """Needed to use objects in sets."""
return id(self) return id(self)