"""
PDB molecular container
=======================
Molecular container for storing all contents of PDB files.
"""
import logging
import os
from typing import Dict, List, Optional, Tuple
from propka.parameters import Parameters
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, Options
_LOGGER = logging.getLogger(__name__)
[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.
"""
conformation_names: List[str]
conformations: Dict[str, ConformationContainer]
name: Optional[str]
version: propka.version.Version
def __init__(self, parameters: Parameters, options: 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) -> None:
"""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) -> None:
"""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) -> None:
"""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) -> 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: Tuple[float, float, float] = (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: List[Tuple[float, float]] = []
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: Tuple[Optional[float], float] = (None, 1e6)
for point in profile:
opt = min(opt, point, key=lambda v: v[1])
# find values within 80 % of optimum
range_80pct: Tuple[Optional[float], Optional[float]] = (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: Tuple[Optional[float], Optional[float]] = (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: str = '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: List[List[float]] = []
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: str = 'AVR', grid=(0., 14.), *,
precision: float = 1e-4) -> Tuple[float, float]:
"""Get the isoelectric points for folded and unfolded states.
Args:
conformation: conformation to test
grid: pH window [min, max]
precision: Compute pI up to this precision
Returns:
1. Folded state PI
2. Unfolded state PI
"""
conf = self.conformations[conformation]
WHICH_UNFOLDED = 0
WHICH_FOLDED = 1
def pi(which, pH, min_, max_):
charge = conf.calculate_charge(
self.version.parameters, ph=pH)[which]
if max_ - min_ > precision:
if charge > 0.0:
min_ = pH
else:
max_ = pH
next_pH = (min_ + max_) / 2
return pi(which, next_pH, min_, max_)
return pH
start = (grid[0] + grid[1]) / 2, grid[0], grid[1]
return (
pi(WHICH_FOLDED, *start),
pi(WHICH_UNFOLDED, *start),
)