De-lint coupled_groups.py

Some methods/attributes were changed but were checked in Google to make
sure other codes were not affected.
This commit is contained in:
Nathan Baker
2020-05-23 17:57:12 -07:00
parent 83c54ec153
commit c985e02713
2 changed files with 216 additions and 173 deletions

View File

@@ -3,7 +3,7 @@ import functools
import propka.ligand import propka.ligand
from propka.output import make_interaction_map from propka.output import make_interaction_map
from propka.determinant import Determinant from propka.determinant import Determinant
from propka.coupled_groups import nccg from propka.coupled_groups import NCCG
from propka.determinants import setBackBoneDeterminants, setIonDeterminants from propka.determinants import setBackBoneDeterminants, setIonDeterminants
from propka.determinants import setDeterminants from propka.determinants import setDeterminants
from propka.group import Group, is_group from propka.group import Group, is_group
@@ -112,7 +112,7 @@ class ConformationContainer:
if len(list(filter(lambda g: len(g.non_covalently_coupled_groups) > 0, if len(list(filter(lambda g: len(g.non_covalently_coupled_groups) > 0,
self.get_titratable_groups()))) > 0: self.get_titratable_groups()))) > 0:
self.non_covalently_coupled_groups = True self.non_covalently_coupled_groups = True
nccg.identify_non_covalently_coupled_groups(self, verbose=verbose) NCCG.identify_non_covalently_coupled_groups(self, verbose=verbose)
# re-do the check # re-do the check
if len(list(filter(lambda g: len(g.non_covalently_coupled_groups) > 0, if len(list(filter(lambda g: len(g.non_covalently_coupled_groups) > 0,
self.get_titratable_groups()))) > 0: self.get_titratable_groups()))) > 0:

View File

@@ -1,113 +1,113 @@
"""Describe coupling between groups."""
from __future__ import division import itertools
from __future__ import print_function import propka.lib
from propka.group import Group
import math, propka.output, propka.group, propka.lib, itertools from propka.output import make_interaction_map
from propka.lib import info, warning from propka.lib import info
class non_covalently_couple_groups: class NonCovalentlyCoupledGroups:
"""Groups that are coupled without covalent bonding."""
def __init__(self): def __init__(self):
self.parameters = None self.parameters = None
# self.do_intrinsic = False
# self.do_pair_wise = False
self.do_prot_stat = True self.do_prot_stat = True
return def is_coupled_protonation_state_probability(self, group1, group2,
energy_method,
return_on_fail=True):
# """Check whether two groups are energetically coupled.
# Methods for finding coupled groups
#
def is_coupled_protonation_state_probability(self, group1, group2, energy_method, return_on_fail=True):
Args:
group1: first group for interaction
group2: second group for interaction
energy_method: function for calculating energy
return_on_fail: return if part of the calculation fails
Returns:
dictionary describing coupling
"""
# check if the interaction energy is high enough # check if the interaction energy is high enough
interaction_energy = max(self.get_interaction(group1,group2), self.get_interaction(group2,group1)) interaction_energy = max(self.get_interaction(group1, group2),
self.get_interaction(group2, group1))
if interaction_energy<=self.parameters.min_interaction_energy and return_on_fail: if interaction_energy <= self.parameters.min_interaction_energy \
return {'coupling_factor':-1.0} and return_on_fail:
return {'coupling_factor': -1.0}
# calculate intrinsic pKa's, if not already done # calculate intrinsic pKa's, if not already done
for group in [group1, group2]: for group in [group1, group2]:
if not hasattr(group, 'intrinsic_pKa'): if not hasattr(group, 'intrinsic_pKa'):
group.calculate_intrinsic_pka() group.calculate_intrinsic_pka()
use_ph = self.parameters.pH
use_pH = self.parameters.pH
if self.parameters.pH == 'variable': if self.parameters.pH == 'variable':
use_pH = min(group1.pka_value, group2.pka_value) use_ph = min(group1.pka_value, group2.pka_value)
default_energy = energy_method(ph=use_ph,
default_energy = energy_method(ph=use_pH, reference=self.parameters.reference) reference=self.parameters.reference)
default_pka1 = group1.pka_value default_pka1 = group1.pka_value
default_pka2 = group2.pka_value default_pka2 = group2.pka_value
# check that pka values are within relevant limits # check that pka values are within relevant limits
if max(default_pka1, default_pka2) < self.parameters.min_pka or \ if max(default_pka1, default_pka2) < self.parameters.min_pka or \
min(default_pka1, default_pka2) > self.parameters.max_pka: min(default_pka1, default_pka2) > self.parameters.max_pka:
if return_on_fail: if return_on_fail:
return {'coupling_factor':-1.0} return {'coupling_factor': -1.0}
# Swap interactions and re-calculate pKa values # Swap interactions and re-calculate pKa values
self.swap_interactions([group1], [group2]) self.swap_interactions([group1], [group2])
group1.calculate_total_pka() group1.calculate_total_pka()
group2.calculate_total_pka() group2.calculate_total_pka()
# store swapped energy and pka's # store swapped energy and pka's
swapped_energy = energy_method(ph=use_pH, reference=self.parameters.reference) swapped_energy = energy_method(ph=use_ph, reference=self.parameters.reference)
swapped_pka1 = group1.pka_value swapped_pka1 = group1.pka_value
swapped_pka2 = group2.pka_value swapped_pka2 = group2.pka_value
pka_shift1 = swapped_pka1 - default_pka1 pka_shift1 = swapped_pka1 - default_pka1
pka_shift2 = swapped_pka2 - default_pka2 pka_shift2 = swapped_pka2 - default_pka2
# Swap back to original protonation state # Swap back to original protonation state
self.swap_interactions([group1], [group2]) self.swap_interactions([group1], [group2])
group1.calculate_total_pka() group1.calculate_total_pka()
group2.calculate_total_pka() group2.calculate_total_pka()
# check difference in free energy # check difference in free energy
if abs(default_energy - swapped_energy) > self.parameters.max_free_energy_diff and return_on_fail: if abs(default_energy - swapped_energy) > self.parameters.max_free_energy_diff \
return {'coupling_factor':-1.0} and return_on_fail:
return {'coupling_factor': -1.0}
# check pka shift # check pka shift
if max(abs(pka_shift1), abs(pka_shift2)) < self.parameters.min_swap_pka_shift and return_on_fail: if max(abs(pka_shift1), abs(pka_shift2)) < self.parameters.min_swap_pka_shift \
return {'coupling_factor':-1.0} and return_on_fail:
return {'coupling_factor': -1.0}
# check intrinsic pka diff # 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) \
return {'coupling_factor':-1.0} > 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 # if everything is OK, calculate the coupling factor and return all info
factor = self.get_free_energy_diff_factor(default_energy, swapped_energy)*\ 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) * self.get_interaction_factor(interaction_energy)
return {'coupling_factor': factor, 'default_energy': default_energy,
return {'coupling_factor':factor, 'swapped_energy': swapped_energy,
'default_energy':default_energy, 'interaction_energy': interaction_energy,
'swapped_energy':swapped_energy, 'swapped_pka1': swapped_pka1, 'swapped_pka2': swapped_pka2,
'interaction_energy':interaction_energy, 'pka_shift1': pka_shift1, 'pka_shift2': pka_shift2,
'swapped_pka1':swapped_pka1, 'pH': use_ph}
'swapped_pka2':swapped_pka2,
'pka_shift1':pka_shift1,
'pka_shift2':pka_shift2,
'pH':use_pH}
#
# Methods for calculating the coupling factor
#
def get_pka_diff_factor(self, pka1, pka2): def get_pka_diff_factor(self, pka1, pka2):
"""Get scaling factor for difference between intrinsic pKa values.
Args:
pka1: first pKa to compare
pka2: second pKa to compare
Returns:
float value of scaling factor
"""
intrinsic_pka_diff = abs(pka1-pka2) intrinsic_pka_diff = abs(pka1-pka2)
res = 0.0 res = 0.0
if intrinsic_pka_diff <= self.parameters.max_intrinsic_pKa_diff: if intrinsic_pka_diff <= self.parameters.max_intrinsic_pKa_diff:
res = 1-(intrinsic_pka_diff/self.parameters.max_intrinsic_pKa_diff)**2 res = 1-(intrinsic_pka_diff/self.parameters.max_intrinsic_pKa_diff)**2
return res return res
def get_free_energy_diff_factor(self, energy1, energy2): def get_free_energy_diff_factor(self, energy1, energy2):
"""Get scaling factor for difference between free energies.
Args:
energy1: first energy to compare
energy2: second energy to compare
Returns:
float value of scaling factor
"""
free_energy_diff = abs(energy1-energy2) free_energy_diff = abs(energy1-energy2)
res = 0.0 res = 0.0
if free_energy_diff <= self.parameters.max_free_energy_diff: if free_energy_diff <= self.parameters.max_free_energy_diff:
@@ -115,19 +115,28 @@ class non_covalently_couple_groups:
return res return res
def get_interaction_factor(self, interaction_energy): def get_interaction_factor(self, interaction_energy):
"""Get scaling factor related to interaction energy.
Args:
interaction_energy: interaction energy
Returns:
float value of scaling factor
"""
res = 0.0 res = 0.0
interaction_energy = abs(interaction_energy) interaction_energy = abs(interaction_energy)
if interaction_energy >= self.parameters.min_interaction_energy: if interaction_energy >= self.parameters.min_interaction_energy:
res = (interaction_energy-self.parameters.min_interaction_energy)/(1.0+interaction_energy-self.parameters.min_interaction_energy) res = (interaction_energy-self.parameters.min_interaction_energy) \
/ (1.0+interaction_energy-self.parameters.min_interaction_energy)
return res return res
def identify_non_covalently_coupled_groups(self, conformation,
verbose=True):
"""Find coupled residues in protein.
Args:
conformation: protein conformation to test
def identify_non_covalently_coupled_groups(self, conformation, verbose=True): verbose: verbose output (boolean)
""" Finds coupled residues in protein """ """
self.parameters = conformation.parameters self.parameters = conformation.parameters
if verbose: if verbose:
info('') info('')
@@ -137,186 +146,220 @@ class non_covalently_couple_groups:
info('-' * 103) info('-' * 103)
info(' Detecting non-covalently coupled residues') info(' Detecting non-covalently coupled residues')
info('-' * 103) info('-' * 103)
info(' Maximum pKa difference: %4.2f pKa units' % self.parameters.max_intrinsic_pKa_diff) info(' Maximum pKa difference: %4.2f pKa units' \
info(' Minimum interaction energy: %4.2f pKa units' % self.parameters.min_interaction_energy) % self.parameters.max_intrinsic_pKa_diff)
info(' Maximum free energy diff.: %4.2f pKa units' % self.parameters.max_free_energy_diff) info(' Minimum interaction energy: %4.2f pKa units' \
info(' Minimum swap pKa shift: %4.2f pKa units' % self.parameters.min_swap_pka_shift) % self.parameters.min_interaction_energy)
info(' Maximum free energy diff.: %4.2f pKa units' \
% self.parameters.max_free_energy_diff)
info(' Minimum swap pKa shift: %4.2f pKa units' \
% self.parameters.min_swap_pka_shift)
info(' pH: %6s ' % str(self.parameters.pH)) info(' pH: %6s ' % str(self.parameters.pH))
info(' Reference: %s' % self.parameters.reference) info(' Reference: %s' % self.parameters.reference)
info(' Min pKa: %4.2f' % self.parameters.min_pka) info(' Min pKa: %4.2f' % self.parameters.min_pka)
info(' Max pKa: %4.2f' % self.parameters.max_pka) info(' Max pKa: %4.2f' % self.parameters.max_pka)
info('') info('')
# find coupled residues # find coupled residues
titratable_groups = conformation.get_titratable_groups() titratable_groups = conformation.get_titratable_groups()
if not conformation.non_covalently_coupled_groups: if not conformation.non_covalently_coupled_groups:
for g1 in titratable_groups: for group1 in titratable_groups:
for g2 in titratable_groups: for group2 in titratable_groups:
if g1==g2: if group1 == group2:
break break
if not group1 in group2.non_covalently_coupled_groups \
if not g1 in g2.non_covalently_coupled_groups and self.do_prot_stat: and self.do_prot_stat:
data = self.is_coupled_protonation_state_probability(g1, g2,conformation.calculate_folding_energy) data = self.\
if data['coupling_factor'] >0.0: is_coupled_protonation_state_probability(group1,
g1.couple_non_covalently(g2) group2,
conformation.\
calculate_folding_energy)
if data['coupling_factor'] > 0.0:
group1.couple_non_covalently(group2)
if verbose: if verbose:
self.print_out_swaps(conformation) self.print_out_swaps(conformation)
return def print_out_swaps(self, conformation):
"""Print out something having to do with coupling interactions.
def print_out_swaps(self, conformation, verbose=True): Args:
map = propka.output.make_interaction_map('Non-covalent coupling map for %s'%conformation, conformation: conformation to print
"""
map_ = make_interaction_map('Non-covalent coupling map for %s' % conformation,
conformation.get_non_covalently_coupled_groups(), conformation.get_non_covalently_coupled_groups(),
lambda g1,g2: g1 in g2.non_covalently_coupled_groups) lambda g1, g2: g1 in g2.non_covalently_coupled_groups)
info(map) info(map_)
for system in conformation.get_coupled_systems(conformation.\
for system in conformation.get_coupled_systems(conformation.get_non_covalently_coupled_groups(),propka.group.Group.get_non_covalently_coupled_groups): get_non_covalently_coupled_groups(), \
Group.get_non_covalently_coupled_groups):
self.print_system(conformation, list(system)) self.print_system(conformation, list(system))
return
def print_system(self, conformation, system): def print_system(self, conformation, system):
"""Print out something about the system.
Args:
conformation: conformation to print
system: system to print
"""
info('System containing %d groups:' % len(system)) info('System containing %d groups:' % len(system))
# make list of interactions within this system
# make list of interactions withi this system interactions = list(itertools.combinations(system, 2))
interactions = list(itertools.combinations(system,2))
# print out coupling info for each interaction # print out coupling info for each interaction
coup_info = '' coup_info = ''
for interaction in interactions: for interaction in interactions:
data = self.is_coupled_protonation_state_probability(interaction[0], interaction[1],conformation.calculate_folding_energy, return_on_fail=False) data = self.is_coupled_protonation_state_probability(interaction[0], \
coup_info += self.make_data_to_string(data,interaction[0],interaction[1])+'\n\n' interaction[1], conformation.calculate_folding_energy, \
return_on_fail=False)
coup_info += self.make_data_to_string(data, interaction[0], \
interaction[1]) + '\n\n'
info(coup_info) info(coup_info)
# make list of possible combinations of swap to try out # make list of possible combinations of swap to try out
combinations = propka.lib.generate_combinations(interactions) combinations = propka.lib.generate_combinations(interactions)
# Make possible swap combinations # Make possible swap combinations
swap_info = '' swap_info = ''
swap_info += self.print_determinants_section(system,'Original') swap_info += self.print_determinants_section(system, 'Original')
for combination in combinations: for combination in combinations:
# Tell the user what is swap in this combination # Tell the user what is swap in this combination
swap_info += 'Swapping the following interactions:\n' swap_info += 'Swapping the following interactions:\n'
for interaction in combination: for interaction in combination:
swap_info += ' %s %s\n'%(interaction[0].label, interaction[1].label) swap_info += ' %s %s\n' % (interaction[0].label,
interaction[1].label)
# swap... # swap...
for interaction in combination: for interaction in combination:
self.swap_interactions([interaction[0]],[interaction[1]]) self.swap_interactions([interaction[0]], [interaction[1]])
swap_info += self.print_determinants_section(system, 'Swapped')
swap_info += self.print_determinants_section(system,'Swapped')
# ...and swap back
#for interaction in combination:
# self.swap_interactions([interaction[0]], [interaction[1]])
info(swap_info) info(swap_info)
return @staticmethod
# def get_interaction(group1, group2, include_side_chain_hbs=True):
# Interaction and swapping methods """Get interaction energy between two groups.
#
def get_interaction(self, group1, group2, include_side_chain_hbs = True): Args:
group1: first group for interaction
group2: second group for interaction
include_side_chain_hbs: include side-chain hydrogen bonds in energy
Returns:
interaction energy (float)
"""
determinants = group1.determinants['coulomb'] determinants = group1.determinants['coulomb']
if include_side_chain_hbs: if include_side_chain_hbs:
determinants = group1.determinants['sidechain'] + group1.determinants['coulomb'] determinants = group1.determinants['sidechain'] \
+ group1.determinants['coulomb']
interaction_energy = 0.0 interaction_energy = 0.0
for det in determinants: for det in determinants:
if group2 == det.group: if group2 == det.group:
interaction_energy += det.value interaction_energy += det.value
return interaction_energy return interaction_energy
def print_determinants_section(self, system, tag): def print_determinants_section(self, system, tag):
"""Print determinants of system.
Args:
system: set of groups
tag: something to add to output
Returns:
string with summary
"""
all_labels = [g.label for g in system] all_labels = [g.label for g in system]
s = ' '+'-'*113+'\n' str_ = ' ' + '-' * 113 + '\n'
for group in system: for group in system:
s += self.tagged_format(' %-8s|' % tag, group.getDeterminantString(), all_labels) str_ += self.tagged_format(' %-8s|' % tag,
group.getDeterminantString(),
all_labels)
return str_ + '\n'
return s+'\n' def swap_interactions(self, groups1, groups2, include_side_chain_hbs=True):
"""Swap interactions between two groups.
def swap_interactions(self, groups1, groups2, include_side_chain_hbs = True): Args:
group1: first group to swap
for i in range(len(groups1)): group2: second group to swap
group1 = groups1[i] """
for i, group1 in enumerate(groups1):
group2 = groups2[i] group2 = groups2[i]
# swap the interactions! # swap the interactions!
self.transfer_determinant(group1.determinants['coulomb'], group2.determinants['coulomb'], group1.label, group2.label) self.transfer_determinant(group1.determinants['coulomb'],
group2.determinants['coulomb'],
group1.label, group2.label)
if include_side_chain_hbs: if include_side_chain_hbs:
self.transfer_determinant(group1.determinants['sidechain'], group2.determinants['sidechain'], group1.label, group2.label) self.transfer_determinant(group1.determinants['sidechain'],
group2.determinants['sidechain'],
group1.label, group2.label)
# re-calculate pKa values # re-calculate pKa values
group1.calculate_total_pka() group1.calculate_total_pka()
group2.calculate_total_pka() group2.calculate_total_pka()
return @staticmethod
def transfer_determinant(determinants1, determinants2,
label1, label2):
"""Transfer information between two sets of determinants.
Args:
def transfer_determinant(self, determinants1, determinants2, label1, label2): determinants1: determinant list
determinants2: determinant list
label1: label for list 1
label2: label for list 2
"""
# find out what to transfer... # find out what to transfer...
from1to2 = [] from1to2 = []
from2to1 = [] from2to1 = []
for det in determinants1: for det in determinants1:
if det.label == label2: if det.label == label2:
from1to2.append(det) from1to2.append(det)
for det in determinants2: for det in determinants2:
if det.label == label1: if det.label == label1:
from2to1.append(det) from2to1.append(det)
# ...and transfer it! # ...and transfer it!
for det in from1to2: for det in from1to2:
det.label = label1 det.label = label1
determinants2.append(det) determinants2.append(det)
determinants1.remove(det) determinants1.remove(det)
for det in from2to1: for det in from2to1:
det.label = label2 det.label = label2
determinants1.append(det) determinants1.append(det)
determinants2.remove(det) determinants2.remove(det)
return @staticmethod
def tagged_format(tag, str_, labels):
"""Tag a string.
# Args:
# Output methods tag: tag to add
# str_: string to tag
labels: labels to replace
def tagged_format(self, tag, s, labels): Returns:
s = "%s %s"%(tag,s) tagged string
s = s.replace('\n','\n%s '%tag) """
str_ = "%s %s" % (tag, str_)
str_ = str_.replace('\n', '\n%s ' % tag)
for label in labels: for label in labels:
s = s.replace(label, '\033[31m%s\033[30m'%label) str_ = str_.replace(label, '\033[31m%s\033[30m' % label)
return s+'\n' return str_ + '\n'
@staticmethod
def make_data_to_string(data, group1, group2):
"""Describe interaction between groups.
def make_data_to_string(self, data, group1, group2): Args:
s = """ %s and %s coupled (prot.state): %5.2f data: data about interactions
group1: first group
group2: second group
Returns:
formatted string with information.
"""
str_ = \
""" %s and %s coupled (prot.state): %5.2f
Energy levels: %6.2f, %6.2f (difference: %6.2f) at pH %6.2f Energy levels: %6.2f, %6.2f (difference: %6.2f) at pH %6.2f
Interaction energy: %6.2f Interaction energy: %6.2f
Intrinsic pka's: %6.2f, %6.2f (difference: %6.2f) Intrinsic pka's: %6.2f, %6.2f (difference: %6.2f)
Swapped pKa's: %6.2f, %6.2f (difference: %6.2f, %6.2f)"""%(group1.label, Swapped pKa's: %6.2f, %6.2f (difference: %6.2f, %6.2f)""" % \
group2.label, (group1.label, group2.label, data['coupling_factor'],
data['coupling_factor'],
data['default_energy'], data['swapped_energy'], data['default_energy'], data['swapped_energy'],
data['default_energy'] - data['swapped_energy'], data['default_energy'] - data['swapped_energy'], data['pH'],
data['pH'], data['interaction_energy'], group1.intrinsic_pKa, group2.intrinsic_pKa,
data['interaction_energy'], group1.intrinsic_pKa-group2.intrinsic_pKa, data['swapped_pka1'],
group1.intrinsic_pKa, data['swapped_pka2'], data['pka_shift1'], data['pka_shift2'])
group2.intrinsic_pKa,
group1.intrinsic_pKa-group2.intrinsic_pKa,
data['swapped_pka1'],
data['swapped_pka2'],
data['pka_shift1'],
data['pka_shift2'])
return s return str_
nccg = non_covalently_couple_groups() NCCG = NonCovalentlyCoupledGroups()