"""
PDB molecular container
=======================
Molecular container for storing all contents of PDB files.
"""
import logging
import os
import propka.version
from propka.output import write_pka, print_header, print_result
from propka.conformation_container import ConformationContainer
from propka.lib import make_grid
_LOGGER = logging.getLogger(__name__)
# TODO - these are constants whose origins are a little murky
UNK_PI_CUTOFF = 0.01
# Maximum number of iterations for finding PI
MAX_ITERATION = 4
[docs]class MolecularContainer:
"""Container for storing molecular contents of PDB files.
TODO - this class name does not conform to PEP8 but has external use.
We should deprecate and change eventually.
.. versionchanged:: 3.4.0
Removed :meth:`write_propka` and
:meth:`additional_setup_when_reading_input_file` as reading and writing
PROPKA input files is no longer supported.
"""
def __init__(self, parameters, options=None):
"""Initialize molecular container.
Args:
parameters: Parameters() object
options: options object
"""
# printing out header before parsing input
print_header()
self.conformation_names = []
self.conformations = {}
self.options = options
self.name = None
try:
version_class = getattr(propka.version, parameters.version)
self.version = version_class(parameters)
except AttributeError as err:
print(err)
errstr = 'Error: Version {0:s} does not exist'.format(
parameters.version)
raise Exception(errstr)
[docs] def find_covalently_coupled_groups(self):
"""Find covalently coupled groups."""
for name in self.conformation_names:
self.conformations[name].find_covalently_coupled_groups()
[docs] def find_non_covalently_coupled_groups(self):
"""Find non-covalently coupled groups."""
verbose = self.options.display_coupled_residues
for name in self.conformation_names:
self.conformations[name].find_non_covalently_coupled_groups(
verbose=verbose)
[docs] def calculate_pka(self):
"""Calculate pKa values."""
# calculate for each conformation
for name in self.conformation_names:
self.conformations[name].calculate_pka(
self.version, self.options)
# find non-covalently coupled groups
self.find_non_covalently_coupled_groups()
# find the average of the conformations
self.average_of_conformations()
# print out the conformation-average results
print_result(self, 'AVR', self.version.parameters)
[docs] def write_pka(self, filename=None, reference="neutral",
direction="folding", options=None):
"""Write pKa information to a file.
Args:
filename: file to write to
reference: reference state
direction: folding vs. unfolding
options: options object
"""
if filename is None:
filename = os.path.join('{0:s}.pka'.format(self.name))
# if the display_coupled_residues option is true, write the results out
# to an alternative pka file
if self.options.display_coupled_residues:
filename = os.path.join('{0:s}_alt_state.pka'.format(self.name))
if (hasattr(self.version.parameters, 'output_file_tag')
and len(self.version.parameters.output_file_tag) > 0):
filename = os.path.join(
'{0:s}_{1:s}.pka'.format(
self.name, self.version.parameters.output_file_tag))
write_pka(
self, self.version.parameters, filename=filename,
conformation='AVR', reference=reference)
[docs] def get_folding_profile(self, conformation='AVR', reference="neutral",
grid=[0., 14., 0.1]):
"""Get a folding profile.
Args:
conformation: conformation to select
reference: reference state
direction: folding direction (folding)
grid: the grid of pH values [min, max, step_size]
options: options object
Returns:
TODO - figure out what these are
1. profile
2. opt
3. range_80pct
4. stability_range
"""
# calculate stability profile
profile = []
for ph in make_grid(*grid):
conf = self.conformations[conformation]
ddg = conf.calculate_folding_energy(ph=ph, reference=reference)
profile.append([ph, ddg])
# find optimum
opt = [None, 1e6]
for point in profile:
opt = min(opt, point, key=lambda v: v[1])
# find values within 80 % of optimum
range_80pct = [None, None]
values_within_80pct = [p[0] for p in profile if p[1] < 0.8*opt[1]]
if len(values_within_80pct) > 0:
range_80pct = [min(values_within_80pct), max(values_within_80pct)]
# find stability range
stability_range = [None, None]
stable_values = [p[0] for p in profile if p[1] < 0.0]
if len(stable_values) > 0:
stability_range = [min(stable_values), max(stable_values)]
return profile, opt, range_80pct, stability_range
[docs] def get_charge_profile(self, conformation='AVR', grid=[0., 14., .1]):
"""Get charge profile for conformation as function of pH.
Args:
conformation: conformation to test
grid: grid of pH values [min, max, step]
Returns:
list of charge state values
"""
charge_profile = []
for ph in make_grid(*grid):
conf = self.conformations[conformation]
q_unfolded, q_folded = conf.calculate_charge(
self.version.parameters, ph=ph)
charge_profile.append([ph, q_unfolded, q_folded])
return charge_profile
[docs] def get_pi(self, conformation='AVR', grid=[0., 14., 1], iteration=0):
"""Get the isoelectric points for folded and unfolded states.
Args:
conformation: conformation to test
grid: grid of pH values [min, max, step]
iteration: iteration number of process
Returns:
1. Folded state PI
2. Unfolded state PI
"""
charge_profile = self.get_charge_profile(
conformation=conformation, grid=grid)
pi_folded = pi_unfolded = [None, 1e6, 1e6]
for point in charge_profile:
pi_folded = min(pi_folded, point, key=lambda v: abs(v[2]))
pi_unfolded = min(pi_unfolded, point, key=lambda v: abs(v[1]))
# If results are not good enough, do it again with a higher sampling
# resolution
pi_folded_value = pi_folded[0]
pi_unfolded_value = pi_unfolded[0]
step = grid[2]
# TODO - need to warn if maximum number of iterations is exceeded
if ((pi_folded[2] > UNK_PI_CUTOFF
or pi_unfolded[1] > UNK_PI_CUTOFF) and iteration < MAX_ITERATION):
pi_folded_value, _ = self.get_pi(
conformation=conformation,
grid=[pi_folded[0]-step, pi_folded[0]+step, step/10.0],
iteration=iteration+1)
_, pi_unfolded_value = self.get_pi(
conformation=conformation,
grid=[pi_unfolded[0]-step, pi_unfolded[0]+step, step/10.0],
iteration=iteration+1)
return pi_folded_value, pi_unfolded_value