"""
Data structures for groups
==========================
Routines and classes for storing groups important to PROPKA calculations.
.. versionchanged:: 3.4.0
Removed :func:`initialize_atom_group` as reading PROPKA inputs is no longer
supported.
"""
import logging
import math
import propka.ligand
import propka.protonate
from propka.ligand_pka_values import LigandPkaValues
from propka.determinant import Determinant
_LOGGER = logging.getLogger(__name__)
# Constants that start with "UNK_" are a mystery to me
UNK_PKA_SCALING = -1.36
PROTONATOR = propka.protonate.Protonate(verbose=False)
#: acids
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}}
#: bases
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}}
[docs]class Group:
"""Class for storing groups important to pKa calculations.
.. versionchanged:: 3.4.0
Removed :meth:`make_covalently_coupled_line` and
:meth:`make_non_covalently_coupled_line` as writing PROPKA inputs is no
longer supported.
"""
def __init__(self, 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.pka_value = 0.0
self.model_pka = 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':
fmt = "{g.residue_type:<3s}{a.res_num:>4d}{a.chain_id:>2s}"
self.label = fmt.format(g=self, a=atom)
elif self.atom.res_name in ['DA ', 'DC ', 'DG ', 'DT ']:
fmt = "{type:1s}{elem:1s}{name:1s}{res_num:>4d}{chain:>2s}"
self.label = fmt.format(
type=self.residue_type[1], elem=atom.element,
name=atom.name.replace('\'', '')[-1], res_num=atom.res_num,
chain=atom.chain_id)
else:
fmt = "{type:<3s}{name:>4s}{chain:>2s}"
self.label = fmt.format(
type=self.residue_type, name=atom.name, chain=atom.chain_id)
# container for squared distances
self.squared_distances = {}
[docs] def couple_covalently(self, other):
"""Couple this group with another group.
Args:
other: other group for coupling
"""
# do the coupling
if other not in self.covalently_coupled_groups:
self.covalently_coupled_groups.append(other)
if self not in other.covalently_coupled_groups:
other.covalently_coupled_groups.append(self)
[docs] def couple_non_covalently(self, other):
"""Non-covalenthly couple this group with another group.
Args:
other: other group for coupling
"""
# do the coupling
if other not in self.non_covalently_coupled_groups:
self.non_covalently_coupled_groups.append(other)
if self not in other.non_covalently_coupled_groups:
other.non_covalently_coupled_groups.append(self)
[docs] def get_covalently_coupled_groups(self):
"""Get covalently coupled groups.
Returns:
list of covalently coupled groups.
"""
return self.covalently_coupled_groups
[docs] def get_non_covalently_coupled_groups(self):
"""Get non-covalently coupled groups.
Returns:
list of covalently coupled groups.
"""
return self.non_covalently_coupled_groups
[docs] 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 det in other.determinants[type_]:
self.share_determinant(det, type_)
# recalculate pka values
self.calculate_total_pka()
the_other.calculate_total_pka()
[docs] def share_determinant(self, new_determinant, type_):
"""Add determinant to this group's list of determinants.
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_]:
if own_determinant.group == new_determinant.group:
# if so, find the average value
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(
Determinant(new_determinant.group, new_determinant.value))
def __eq__(self, other):
"""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
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))
def __hash__(self):
"""Needed for creating sets of groups."""
return id(self)
def __iadd__(self, other):
if self.type != other.type:
str_ = (
'Cannot add groups of different types '
'({0:s} and {1:s})'.format(self.type, other.type))
raise Exception(str_)
# add all values
self.pka_value += other.pka_value
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
# 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
[docs] def add_determinant(self, new_determinant, type_):
"""Add 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_]:
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(Determinant(new_determinant.group,
new_determinant.value))
[docs] def set_determinant(self, new_determinant, type_):
"""Overwrite current and create 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_]:
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(Determinant(new_determinant.group,
new_determinant.value))
[docs] def remove_determinants(self, labels):
"""Remove all determinants with specified labels.
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.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_]:
determinant.value /= value
return self
[docs] 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
res.model_pka = self.model_pka
res.coupled_titrating_group = self.coupled_titrating_group
res.covalently_coupled_groups = self.covalently_coupled_groups
res.non_covalently_coupled_groups = self.non_covalently_coupled_groups
res.titratable = self.titratable
res.exclude_cys_from_results = self.exclude_cys_from_results
res.charge = self.charge
return res
[docs] 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]
# 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():
if not self.model_pka_set:
self.model_pka = self.parameters.model_pkas[self.residue_type]
# check if we should apply a custom model pka
key = '{0:s}-{1:s}'.format(
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
[docs] def setup_atoms(self):
"""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])
[docs] 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, _] 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 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:
str_ = 'Missing atoms or failed protonation for '
str_ += '{0:s} ({1:s}) -- please check the structure'.format(
self.label, self.type)
_LOGGER.warning(str_)
_LOGGER.warning('{0:s}'.format(str(self)))
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()])
_LOGGER.warning(
'Expected {0:d} interaction atoms for acids, found:'.format(
num_acid))
for i in range(len(self.interaction_atoms_for_acids)):
_LOGGER.warning(
' {0:s}'.format(
str(self.interaction_atoms_for_acids[i])))
_LOGGER.warning(
'Expected {0:d} interaction atoms for bases, found:'.format(
num_base))
for i in range(len(self.interaction_atoms_for_bases)):
_LOGGER.warning(
' {0:s}'.format(
str(self.interaction_atoms_for_bases[i])))
[docs] 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:
# default is acid interaction atoms - cf. 3.0
return self.interaction_atoms_for_acids
[docs] 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 position of atoms
for atom in atoms:
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))
[docs] def get_determinant_string(self, remove_penalised_group=False):
"""Create a string to identify this determinant.
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_ = ""
for line_number in range(number_of_lines):
str_ += "{0:s}".format(self.label)
if line_number == 0:
str_ += " {0:6.2f}".format(self.pka_value)
if len(self.non_covalently_coupled_groups) > 0:
str_ += '*'
else:
str_ += ' '
str_ += " {0:4d}{1:>2s} ".format(int(100.0*self.buried), "%")
str_ += " {0:6.2f} {1:4d}".format(
self.energy_volume, int(self.num_volume))
str_ += " {0:6.2f} {1:4d}".format(
self.energy_local, int(self.num_local))
else:
str_ += "{0:>40s}".format(" ")
# add the determinants
for type_ in ['sidechain', 'backbone', 'coulomb']:
str_ += self.get_determinant_for_string(type_, line_number)
# adding end-of-line
str_ += "\n"
str_ += "\n"
return str_
[docs] def get_determinant_for_string(self, type_, number):
"""Return a string describing determinant.
Args:
type_: determinant type
number: determinant index number
Returns:
string
"""
if number >= len(self.determinants[type_]):
return " 0.00 XXX 0 X"
else:
determinant = self.determinants[type_][number]
return "{0:8.2f} {1:s}".format(
determinant.value, determinant.label)
[docs] 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.energy_volume + self.energy_local)
for determinant_type in ['sidechain', 'backbone', 'coulomb']:
for determinant in self.determinants[determinant_type]:
self.pka_value += determinant.value
[docs] 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.
"""
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+ ']:
value = determinant.value
side_chain += value
self.intrinsic_pka = (
self.model_pka + self.energy_volume + self.energy_local
+ back_bone + side_chain)
[docs] def get_summary_string(self, remove_penalised_group=False):
"""Create summary string for this group.
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 {0:s}'.format(
self.coupled_titrating_group.label))
fmt = (
" {g.label:>9s} {g.pka_value:8.2f} {g.model_pka:10.2f} "
"{type:>18s} {penalty:s}\n")
return fmt.format(g=self, type=ligand_type, penalty=penalty)
def __str__(self):
str_ = 'Group ({0:s}) for {1:s}'.format(self.type, str(self.atom))
return str_
[docs] 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
"""
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
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 = UNK_PKA_SCALING*(pka_prime - self.model_pka)
# calculating the ddg(low-pH --> pH) contribution
# folded
dpka = ph - self.pka_value
conc_ratio = 10**dpka
q_pro = math.log10(1+conc_ratio)
# unfolded
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
[docs] 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":
q_dpka = self.charge * (self.model_pka - ph)
else:
q_dpka = self.charge * (self.pka_value - ph)
conc_ratio = 10**q_dpka
charge = self.charge*(conc_ratio/(1.0+conc_ratio))
return charge
[docs] def use_in_calculations(self):
"""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)
[docs]class COOGroup(Group):
"""Carboxyl group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'COO'
[docs] 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])
# 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)
[docs]class HISGroup(Group):
"""Histidine group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'HIS'
[docs] 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:
_LOGGER.warning(
'His group does not seem to contain a ring %s', self
)
# protonate ring
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])
# 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)
[docs]class CYSGroup(Group):
"""Cysteine group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'CYS'
[docs]class TYRGroup(Group):
"""Tyrosine group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'TYR'
[docs]class LYSGroup(Group):
"""Lysine group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'LYS'
[docs]class ARGGroup(Group):
"""Arginine group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'ARG'
[docs] 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 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)
[docs]class ROHGroup(Group):
"""Alcohol group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'ROH'
[docs]class SERGroup(Group):
"""Serine group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'SER'
[docs]class AMDGroup(Group):
"""Amide group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'AMD'
[docs] 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')
if not (the_oxygen and the_nitrogen):
_LOGGER.warning(f"Missing N or O atom: {self}")
self.set_center([self.atom])
return
# add protons to the nitrogen
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)
[docs]class TRPGroup(Group):
"""Tryptophan group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'TRP'
[docs] 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
PROTONATOR.protonate_atom(self.atom)
the_hydrogen = self.atom.get_bonded_elements('H')
self.set_interaction_atoms(the_hydrogen+[self.atom], [self.atom])
[docs]class NtermGroup(Group):
"""N-terminus group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'N+'
[docs]class CtermGroup(Group):
"""C-terminus group."""
def __init__(self, atom):
Group.__init__(self, atom)
# this is to deal with the COO-C- parameter unification.
self.type = 'COO'
[docs] 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])
# 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)
[docs]class BBNGroup(Group):
"""Backbone nitrogen group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'BBN'
self.residue_type = 'BBN'
[docs] def setup_atoms(self):
"""Set up atoms in group."""
# Identify the hydrogen
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])
[docs]class BBCGroup(Group):
"""Backbone carbon group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'BBC'
self.residue_type = 'BBC'
[docs] 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)
[docs]class NARGroup(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'NAR'
self.residue_type = 'NAR'
_LOGGER.info('Found NAR group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in group."""
# Identify the hydrogen
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])
[docs]class NAMGroup(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'NAM'
self.residue_type = 'NAM'
_LOGGER.info('Found NAM group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the hydrogen
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])
[docs]class FGroup(Group):
"""Fluoride group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'F'
self.residue_type = 'F'
_LOGGER.info('Found F group: %s', atom)
[docs]class ClGroup(Group):
"""Chloride group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'Cl'
self.residue_type = 'Cl'
_LOGGER.info('Found Cl group: %s', atom)
[docs]class OHGroup(Group):
"""Hydroxide group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'OH'
self.residue_type = 'OH'
_LOGGER.info('Found OH group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the hydrogen
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])
[docs]class OPGroup(Group):
"""Phosphate group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'OP'
self.residue_type = 'OP'
_LOGGER.info('Found OP group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the hydrogen
PROTONATOR.protonate_atom(self.atom)
# set the center using the oxygen
self.set_center([self.atom])
self.set_interaction_atoms([self.atom], [self.atom])
[docs]class O3Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'O3'
self.residue_type = 'O3'
_LOGGER.info('Found O3 group: %s', atom)
[docs]class O2Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'O2'
self.residue_type = 'O2'
_LOGGER.info('Found O2 group: %s', atom)
[docs]class SHGroup(Group):
"""Sulfhydryl group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'SH'
self.residue_type = 'SH'
_LOGGER.info('Found SH group: %s', atom)
[docs]class CGGroup(Group):
"""Guadinium group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'CG'
self.residue_type = 'CG'
_LOGGER.info('Found CG group: %s', atom)
[docs] 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 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)
[docs]class C2NGroup(Group):
"""Amidinium group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'C2N'
self.residue_type = 'C2N'
_LOGGER.info('Found C2N group: %s', atom)
[docs] 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]
# set the center using the nitrogen
self.set_center([self.atom])
the_hydrogens = []
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)
[docs]class OCOGroup(Group):
"""Carboxyl group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'OCO'
self.residue_type = 'OCO'
_LOGGER.info('Found OCO group: %s', atom)
[docs] def setup_atoms(self):
"""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)
[docs]class N30Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'N30'
self.residue_type = 'N30'
_LOGGER.info('Found N30 group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the nitrogens
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])
[docs]class N31Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'N31'
self.residue_type = 'N31'
_LOGGER.info('Found N31 group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the nitrogens
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])
[docs]class N32Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'N32'
self.residue_type = 'N32'
_LOGGER.info('Found N32 group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the nitrogens
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])
[docs]class N33Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'N33'
self.residue_type = 'N33'
_LOGGER.info('Found N33 group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in this group."""
# Identify the nitrogens
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])
[docs]class NP1Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'NP1'
self.residue_type = 'NP1'
_LOGGER.info('Found NP1 group: %s', atom)
[docs] def setup_atoms(self):
"""Set up atoms in group."""
# Identify the nitrogens
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])
[docs]class N1Group(Group):
"""Unknown group.
TODO - identify this group.
"""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'N1'
self.residue_type = 'N1'
_LOGGER.info('Found N1 group: %s', atom)
[docs]class IonGroup(Group):
"""Ion group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'ION'
self.residue_type = atom.res_name.strip()
_LOGGER.info('Found ion group: %s', atom)
[docs]class NonTitratableLigandGroup(Group):
"""Non-titratable ligand group."""
def __init__(self, atom):
Group.__init__(self, atom)
self.type = 'LG'
self.residue_type = 'LG'
[docs]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:
self.type = 'ALG'
self.residue_type = 'ALG'
elif self.charge > 0:
self.type = 'BLG'
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:
self.model_pka = atom.marvin_pka
_LOGGER.info(
'Titratable ligand group %s %s %s',
atom, self.model_pka, self.charge
)
self.model_pka_set = True
[docs]def is_group(parameters, atom):
"""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
# check if this atom belongs to a ion group
ion_group = is_ion_group(parameters, atom)
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 = None
elif parameters.ligand_typing == 'groups':
ligand_group = is_ligand_group_by_groups(parameters, atom)
else:
raise Exception(
'Unknown ligand typing method \'{0.s}\''.format(
parameters.ligand_typing))
if ligand_group:
return ligand_group
return None
[docs]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 NtermGroup(atom)
elif atom.terminal == 'C-':
return CtermGroup(atom)
# Backbone
if atom.type == 'atom' and atom.name == 'N':
# ignore proline backbone nitrogens
if atom.res_name != 'PRO':
return BBNGroup(atom)
if atom.type == 'atom' and atom.name == 'C':
# ignore C- carboxyl
if atom.count_bonded_elements('O') == 1:
return BBCGroup(atom)
# Filters for side chains based on PDB protein atom names
key = '{0:s}-{1:s}'.format(atom.res_name, atom.name)
if key in parameters.protein_group_mapping.keys():
class_str = "{0:s}Group".format(parameters.protein_group_mapping[key])
group_class = globals()[class_str]
return group_class(atom)
return None
[docs]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
PROTONATOR.protonate_atom(atom)
if atom.sybyl_type == 'N.ar':
if len(atom.get_bonded_heavy_atoms()) == 2:
return NARGroup(atom)
if atom.sybyl_type == 'N.am':
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 N30Group(atom)
elif len(heavy_bonded) == 1:
return N31Group(atom)
elif len(heavy_bonded) == 2:
return N32Group(atom)
elif len(heavy_bonded) == 3:
return N33Group(atom)
if atom.sybyl_type == 'N.1':
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 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)]
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]
if len(n_with_max_two_heavy_atom_bonds) == 2:
return C2NGroup(atom)
if len(bonded_nitrogens) == 3:
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 OCOGroup(atom)
if atom.sybyl_type == 'F':
return FGroup(atom)
if atom.sybyl_type == 'Cl':
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 OPGroup(atom)
# hydroxyl group
else:
return OHGroup(atom)
# other SP3 oxygen
else:
return O3Group(atom)
if atom.sybyl_type == 'O.2':
return O2Group(atom)
if atom.sybyl_type == 'S.3':
# sulphide group
if len(atom.get_bonded_heavy_atoms()) == 1:
return SHGroup(atom)
return None
[docs]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 = LigandPkaValues(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)
if atom.marvin_pka:
return TitratableLigandGroup(atom)
# Special case of oxygen in carboxyl group not assigned 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]
atom.marvin_pka = other_oxygen.marvin_pka
return TitratableLigandGroup(atom)
if atom.element in parameters.hydrogen_bonds.elements:
return NonTitratableLigandGroup(atom)
return None
[docs]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 IonGroup(atom)
return None