"""
This module defines the InteractionCalculator class, which can be used for
finding interactions between two groups of atoms.
"""
import operator
from collections import Counter
from collections import OrderedDict
from collections import defaultdict
from collections import namedtuple
from past.utils import old_div
import numpy
from schrodinger.application.bioluminate import surfcomp
from schrodinger.application.bioluminate.protein import ALL_RESIDUES
from schrodinger.application.bioluminate.protein import NONPOLAR_RESIDUES
from schrodinger.infra import mm
from schrodinger.infra import mmlist
from schrodinger.infra import structure as infrastructure
from schrodinger.structutils import analyze
from schrodinger.structutils import interactions
from schrodinger.structutils import measure
ORIG_ATOM_NUM_PROP = 'i_temp_orig_atom_num'
# Interaction parameter defaults
HBOND_MAX_DIST = 2.5
HBOND_MIN_A_ANGLE = 90
HBOND_MIN_D_ANGLE = 120
ALLOWABLE_VDW_OVERLAP = 0.4
SALT_BRIDGE_MAX_DIST = 4.0
NEIGHBOR_MAX_DIST = 4.0
MAX_STACK_DIST = 4.0
RESIDUE_HYDROPHOBIC_ATOMS = {
# hydrophobic
'ALA': [' CB '],
'VAL': [' CB ', ' CG1', ' CG2'],
'LEU': [' CB ', ' CG ', ' CD1', ' CD2'],
'ILE': [' CB ', ' CG1', ' CG2', ' CD1'],
'PRO': [' CB ', ' CG ', ' CD '],
'PHE': [' CB ', ' CG ', ' CD1', ' CD2', ' CE1', ' CE2', ' CZ '],
'MET': [' CB ', ' CG ', ' SD ', ' CE '],
'TRP': [
' CB ', ' CG ', ' CD1', ' CD2', ' NE1', ' CE2', ' CE3', ' CZ2', ' CZ3',
' CH2'
],
# polar
'TYR': [' CB ', ' CG ', ' CD1', ' CE1', ' CE2', ' CZ ']
}
HYDROPHOBIC_DIST_CUTOFF = 8.0
EP = 0.00001
# When calculating SASA, floats within this value are considered equal
# (i.e. epsilon)
[docs]def approx_eq(val1, val2):
return (val1 - EP <= val2 <= val1 + EP)
[docs]def lipophilic_ChemScore_value(atom1, atom2, distance=None, cutoff=None):
"""
Calculates a lipophilic score between an atom pair utilizing their van der
Waals radii. An optional argument is available to avoid computation of the
iter-atom distance within the function. The empirical function form is
adopted from ChemScore:
Empirical scoring functions: I. The development of a fast empirical scoring
function to estimate the binding affinity of ligands in receptor complexes.
Eldridge, Murray, Auton, Paolini, and Mee. JCAMD, 1997 (11).
:param atom1: first atom (usually a lipophilic ligand atom)
:type atom1: `structure.StructureAtom`
:param atom2: second atom (usually a lipophilic receptor atom)
:type atom2: `structure.StructureAtom`
:param distance: inter-atom distance
:type distance: float
:param cutoff: distance beyond which to always return 0.0
:type cutoff: float
:return: emperical function value
:rtype: float
"""
if distance is None:
distance = measure.measure_distance(atom1, atom2)
if cutoff and distance >= cutoff:
return 0.0
# Parameters as defined in "Empirical scoring functions" p. 431
R1 = atom1.radius + atom2.radius + 0.5
R2 = R1 + 3.0
# Function form as per Fig. 1
if distance <= R1:
return 1.0
elif distance >= R2:
return 0.0
return old_div((R2 - distance), 3.0)
[docs]class InteractionParams(object):
"""
A class to store settings for interaction parameters
"""
[docs] def __init__(self):
"""
Initialize the class using the default param
"""
self.hbond_max_dist = HBOND_MAX_DIST
self.hbond_min_a_angle = HBOND_MIN_A_ANGLE
self.hbond_min_d_angle = HBOND_MIN_D_ANGLE
self.allowable_vdw_overlap = ALLOWABLE_VDW_OVERLAP
self.salt_bridge_max_dist = SALT_BRIDGE_MAX_DIST
self.neighbor_max_dist = NEIGHBOR_MAX_DIST
self.max_stack_dist = MAX_STACK_DIST
[docs] def paramDict(self):
"""
Return a dictionary of all interaction parameters. Note that changes to
this dictionary will change the instance variables as well.
:return: A dictionary of all interaction parameters
:rtype: dict
"""
return vars(self)
_ResTuple = namedtuple("_ResTuple", ("chain", "resnum", "inscode", "pdbres"))
[docs]class ResTuple(_ResTuple):
"""
A class to store a residue. Unlike schrodinger.structure._Residue, two
ResTuple objects that describe the same residue are equal (and their hashes
are equal). This class will also strip spaces from inscode and pdbres.
"""
def __new__(cls, res):
"""
Create a new instance of the tuple
:param res: The object to take residue data from
:type res: `schrodinger.structure._StructureAtom` or
`schrodinger.structure._Residue` or `ResTuple`
:return: The newly created instance
:rtype: `ResTuple`
"""
chain = res.chain
resnum = res.resnum
inscode = res.inscode
pdbres = res.pdbres.strip().capitalize()
return tuple.__new__(cls, (chain, resnum, inscode, pdbres))
[docs] @classmethod
def init(cls, *args, **kwargs):
"""
Initialize a class instance directly instead of from a _Residue object.
This function allows for the output of __repr__() to be evaluated and
is intended for use in testing.
:return: The initialized instance
:rtype: `ResTuple`
"""
return _ResTuple.__new__(cls, *args, **kwargs)
def __str__(self):
"""
Represent this residue as ex. A:22:Ala
"""
(chain, resnum, inscode, pdbres) = self
if chain == " ":
chain = "_"
if inscode == " ":
inscode = ""
return "%s:%i%s:%s" % (chain, resnum, inscode, pdbres)
def __repr__(self):
"""
Generate a string that can be used to reproduce this object
"""
return self.__class__.__name__ + ".init" + str(tuple(self))
[docs]class Interactions(object):
"""
Store data about all the interactions made by a single residue
"""
# These are the names that will be used in the interaction summary (i.e. the
# Specific Interactions column)
H_BOND = "hb"
PI_STACK = "pi stack"
SALT_BRIDGE = "salt bridge"
CLASH = "clash"
DISULFIDE = "s-s"
[docs] def __init__(self):
"""
Initialize an instance of this class with 0 interactions
"""
self.hbond = Counter()
self.pi_stack = Counter()
self.salt_bridge = Counter()
self.clash = Counter()
self.disulfide = Counter()
self.nearby_res = defaultdict(lambda: float("inf"))
self.interacting_atoms = defaultdict(set)
self.buried_sasa = 0.0
self.vdw_comp = 0.0
self._createInteractionsDict()
def _createInteractionsDict(self):
"""
Create self._interactions that contains all of the specific interaction
counters
"""
interactions = [(self.H_BOND, self.hbond),
(self.PI_STACK, self.pi_stack),
(self.SALT_BRIDGE, self.salt_bridge),
(self.CLASH, self.clash),
(self.DISULFIDE, self.disulfide)]
self.interactions = OrderedDict(interactions)
[docs] @classmethod
def init(cls, **kwargs):
"""
Initialize a class instance from existing interaction dictionaries.
This function allows for the output of __repr__() to be evaluated and
is intended for use in testing.
:return: The initialized instance
:rtype: `Interactions`
"""
interactions = cls()
for (name, attr) in kwargs.items():
setattr(interactions, name, attr)
interactions._createInteractionsDict()
return interactions
def __repr__(self):
"""
Generate a string that can be used to reproduce this object
"""
args = []
for (attr, val) in sorted(vars(self).items()):
if attr == "interactions":
continue
elif isinstance(val, defaultdict):
# defaultdicts don't stringify as executable code due to the
# "<type 'XXX'>" type listing
val = dict(val)
cur_arg = "%s=%s" % (attr, str(val))
args.append(cur_arg)
args = ", ".join(args)
return self.__class__.__name__ + ".init(" + args + ")"
def __eq__(self, other):
"""
Test two Interactions objects for equality. Note that nearby residue
distances are compared for approximate equality since they're floats.
"""
for (attr, val) in sorted(vars(self).items()):
if attr == "interactions":
continue
other_val = getattr(other, attr)
if isinstance(val, float):
if not approx_eq(val, other_val):
return False
elif attr == "nearby_res":
# The distances stored in nearby_res.values() are floats, so we
# need to allow for rounding error in the comparisons
if set(list(val)) != set(list(other_val)):
return False
for (key, self_dist) in val.items():
other_dist = other_val[key]
if not approx_eq(self_dist, other_dist):
return False
elif val != other_val:
return False
return True
def __ne__(self, other):
"""
Test two Interactions objects for inequality.
"""
return not self == other
[docs] def numHbonds(self):
"""
Return the number of hydrogen bonds
:return: The number of hydrogen bonds
:rtype: int
"""
return sum(self.hbond.values())
[docs] def numPiStacks(self):
"""
Return the number of pi stacks
:return: The number of pi stacks
:rtype: int
"""
return sum(self.pi_stack.values())
[docs] def numSaltBridges(self):
"""
Return the number of salt bridges
:return: The number of salt bridges
:rtype: int
"""
return sum(self.salt_bridge.values())
[docs] def numDisulfides(self):
"""
Return the number of disulfide bonds
:return: The number of disulfide bonds
:rtype: int
"""
return sum(self.disulfide.values())
[docs] def numClashs(self):
"""
Return the number of steric clashes
:return: The number of steric clashes
:rtype: int
"""
return sum(self.clash.values())
[docs] def numSpecificInteractions(self):
"""
Return the total number of specific interactions (i.e. ignoring non-
specific interactions such as buried SASA)
:return: The total number of specific interactions
:rtype: int
"""
total = 0
for cur_interaction in self.interactions.values():
total += sum(cur_interaction.values())
return total
[docs] def allInteractingResidues(self):
"""
Return a list of all residues that this one interacts with (not counting
nearby residues)
:return: A set of ResTuple objects
:rtype: set
"""
# Build a set of all residues that are interacted with
interacting_res = set()
for cur_interaction in self.interactions.values():
interacting_res.update(list(cur_interaction))
return interacting_res
[docs] def interactionSummary(self):
"""
Create the interaction summary to display in the specific interactions
column.
:return: The interaction summary
:rtype: str
"""
interacting_res = self.allInteractingResidues()
all_res_inters = []
for cur_res in sorted(interacting_res):
cur_res_inters = []
for (inter_name, cur_inter) in self.interactions.items():
num_inters = cur_inter[cur_res]
if not num_inters:
continue
inter_desc = "%ix %s" % (num_inters, inter_name)
cur_res_inters.append(inter_desc)
if cur_res_inters:
res_desc = ", ".join(cur_res_inters)
res_desc += " to %s" % str(cur_res)
all_res_inters.append(res_desc)
return "\n".join(all_res_inters)
[docs] def nearbyRes(self):
"""
Return the nearby residues
:return: A tuple of
- The nearby residues (as strings), sorted by distance
- The sorted residue distances (floats)
:rtype: tuple
"""
sorted_nearby = sorted(self.nearby_res.items(),
key=operator.itemgetter(1))
nearby_res = [str(res) for (res, dist) in sorted_nearby]
nearby_dist = [dist for (res, dist) in sorted_nearby]
return nearby_res, nearby_dist
[docs]class InteractingResidue():
"""
Store information about a residue and the interactions it makes
"""
[docs] def __init__(self, res, interactions):
"""
Initialize an instance from existing ResTuple and Interactions objects
:param res: The residue
:type res: `ResTuple`
:param interactions: The interactions made by `res`
:type interactions: `Interactions`
"""
# Extract information from the ResTuple object
self.chain = res.chain
self.resnum = res.resnum
self.inscode = res.inscode
self.pdbres = res.pdbres
self.res_str = str(res)
# Extract information from the Interactions object
self.num_hbonds = interactions.numHbonds()
self.num_pi_stacks = interactions.numPiStacks()
self.num_salt_bridges = interactions.numSaltBridges()
self.num_clashs = interactions.numClashs()
self.num_disulfides = interactions.numDisulfides()
self.num_specific_interactions = interactions.numSpecificInteractions()
self.buried_sasa = interactions.buried_sasa
self.vdw_comp = interactions.vdw_comp
self.interaction_summary = interactions.interactionSummary()
self.nearby_res, self.nearby_dist = interactions.nearbyRes()
[docs]class InteractionCalculator(object):
"""
Calculate all interactions between two groups of atoms
"""
[docs] def __init__(
self,
interaction_params=InteractionParams(), # noqa: M511
ignore_backbone=False):
"""
Initialize an instance of the class using the specified parameters
:param interaction_params: The interaction parameters
:type interaction_params: `interaction_calculator.InteractionParams`
:param ignore_backbone: Should the calculations ignore backbone-backbone
interactions?
:type ignore_backbone: bool
"""
self.interaction_params = interaction_params
self._ignore_backbone = ignore_backbone
self.results = None # store interactions for each residue
self.results_overall = None # store overall interface properties
[docs] def calculate(self, struc, asl_expressions):
"""
Calculate all interactions
:param struc: The structure to analyze
:type struc: `schrodinger.structure.Structure`
:param asl_expressions: A list of [asl_expresion for group 1,
asl expression for group2]
:type asl_expressions: list of string
"""
retval = self._prepCalculations(struc, asl_expressions)
(group_atoms, group_strucs, combined_struc) = retval
self._calculateHbonds(struc, group_atoms)
self._calculateSaltBridges(struc, group_atoms)
self._calculatePiStacks(struc, group_strucs)
self._calculateDisulfides(struc, group_atoms)
self._calculateClashes(struc, group_atoms)
# Note that calculateClashes must be called after calculating all
# other specific interactions so that we don't count those
# interactions as clashes
self._calculateHydrophobicScore(struc, group_atoms)
self._calculateNearbyRes(struc, group_atoms)
self._calculateBuriedSASA(group_strucs, combined_struc)
self._calculateVdwComp(struc, group_atoms)
self._removeTempProperty(struc)
def _prepCalculations(self, struc, asl_expressions):
"""
Prepare the specified structure for calculations
:param struc: The structure to analyze
:type struc: `schrodinger.structure.Structure`
:param asl_expressions: A list of [asl_expresion for group 1,
asl expression for group2]
:type asl_expressions: list of string
:return: A tuple of (1) A list of [all protein atom indices in group 1,
all protein atom indices in group 2] (2) A list of [a Structure
object containing only the protein atoms in group 1, a Structure
object containing only the protein atoms in group 2] (3) A single
Structure object containing the protein atoms in groups 1 and 2
"""
self.results = defaultdict(Interactions)
self.results_overall = defaultdict(float)
# Build a list of all atom numbers for all protein atoms in the
# specified areas
group_atoms = [
analyze.evaluate_asl(struc, "( %s ) and ( protein )" % asl)
for asl in asl_expressions
]
self._checkAtomLists(group_atoms)
(group_strucs, combined_struc) = self._prepStrucs(struc, group_atoms)
return group_atoms, group_strucs, combined_struc
def _prepStrucs(self, struc, group_atoms):
"""
Divide up a structure into two new structures, each of which contains
only the specified atoms. In each new structure, the atom indices from
the original structure are recorded in a temporary atom property.
:param struc: The structure to be divided up
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:return: A tuple of two new structures, one for each group.
"""
for atom in struc.atom:
atom.property[ORIG_ATOM_NUM_PROP] = atom.index
all_group_atoms = sorted(set(group_atoms[0] + group_atoms[1]))
combined_struc = struc.extract(all_group_atoms, copy_props=True)
group_strucs = []
for cur_group_atom_nums in group_atoms:
new_struc = struc.extract(cur_group_atom_nums, copy_props=True)
group_strucs.append(new_struc)
return group_strucs, combined_struc
def _removeTempProperty(self, struc):
"""
Remove the ORIG_ATOM_NUM_PROP atom property from all atoms
:param struc: The structure to remove the property from
:type struc: `schrodinger.structure.Structure`
"""
for atom in struc.atom:
del (atom.property[ORIG_ATOM_NUM_PROP])
[docs] def compileOneToManyResults(self):
"""
Compile all of the calculated interactions into a list of
InteractingResidue objects. Each `InteractionResidue` object is a
combined representation of all nearby interacting residues.
:return: A list of InteractingResidue objects, sorted by residue
:rtype: list
"""
compiled_results = []
for (res, interactions) in self.results.items():
interacting_res = InteractingResidue(res, interactions)
compiled_results.append(interacting_res)
return compiled_results
[docs] def compileOneToOneResults(self):
"""
Compile all of the calculated interactions into a list of
InteractingResidue objects. Each `InteractionResidue` object should
only store interaction data of a single nearby residue.
:return: List of individual `InteractionResidue` objects.
:rtype: list[InteractionResidue]
"""
_all_interacting_residues = []
for curr_res, curr_interactions in self.results.items():
for nearby_res in curr_interactions.nearby_res:
interaction = Interactions()
# Only store individual nearby residue interaction information.
interaction.hbond[nearby_res] = \
curr_interactions.hbond[nearby_res]
interaction.pi_stack[nearby_res] = \
curr_interactions.pi_stack[nearby_res]
interaction.salt_bridge[nearby_res] = \
curr_interactions.salt_bridge[nearby_res]
interaction.clash[nearby_res] = \
curr_interactions.clash[nearby_res]
interaction.disulfide[nearby_res] = \
curr_interactions.disulfide[nearby_res]
interaction.nearby_res[nearby_res] = \
curr_interactions.nearby_res[nearby_res]
interaction.buried_sasa = curr_interactions.buried_sasa
interaction.vdw_comp = curr_interactions.vdw_comp
inter_res = InteractingResidue(curr_res, interaction)
_all_interacting_residues.append(inter_res)
# Some residues don't have any neighboring residue interactions
# but may have buried_sasa or vdw_clash values.
if not curr_interactions.nearby_res and (
curr_interactions.vdw_comp or
curr_interactions.buried_sasa):
_all_interacting_residues.append(
InteractingResidue(curr_res, curr_interactions))
return _all_interacting_residues
[docs] @classmethod
def run(
cls,
struc,
asl_expressions,
interaction_params=InteractionParams(), # noqa: M511
ignore_backbone=False,
one_to_one_interactions=False):
"""
A convenience function to initialize this class, calculate all
interactions, and return the compiled results.
:param struc: The structure to analyze
:type struc: `schrodinger.structure.Structure`
:param asl_expressions: A list of [asl_expresion for group 1,
asl expression for group2]
:type asl_expressions: list of string
:param interaction_params: The interaction parameters
:type interaction_params: `interaction_calculator.InteractionParams`
:param ignore_backbone: Should the calculations ignore backbone-backbone
interactions?
:type ignore_backbone: bool
:param one_to_one_interactions: Whether to return a compiled one to
one mapping or one to many mapping of interacting residues.
:type one_to_one_interactions: bool
:return: A list of InteractingResidue objects describing all calculated
interactions, sorted by residue
:rtype: list
"""
calculator = cls(interaction_params, ignore_backbone)
calculator.calculate(struc, asl_expressions)
if one_to_one_interactions:
return calculator.compileOneToOneResults()
else:
return calculator.compileOneToManyResults()
def _checkAtomLists(self, group_atoms):
"""
Make sure that both groups contain protein atoms
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:raise InteractionCalculatorError: If either group does not contain any
protein atoms
"""
if not group_atoms[0]:
bad_group = 1
elif not group_atoms[1]:
bad_group = 2
else:
return
msg = "Group %s does not contain any protein atoms." % (bad_group)
raise InteractionCalculatorError(msg)
def _calculateInteraction(self, iter_func, inter_name, struc, group_atoms,
check_interaction_func):
"""
Record all interactions of the specified type to `self.results`
:param iter_func: A function that produces an iterator that iterates
through all interactions to record
:type iter_func: function
:param inter_name: The interaction name. Must be one of the class
constants from the `Interactions` class.
:type inter_name: str
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:param check_interaction_func: A function that either records
interacting atoms so they're not later counted as clashes
(self._recordInteraction) or checks clashing atoms to see if they've
previously been recorded as interactions (self._checkInteraction). This
function must return False for all interactions to be ignored and True
otherwise.
:type check_interaction_func: function
"""
interaction_iter = iter_func(struc, group_atoms)
count = 0
for (atom1, atom2) in interaction_iter:
ha1_num = self._getHeavyAtomNum(struc, atom1)
ha2_num = self._getHeavyAtomNum(struc, atom2)
if not check_interaction_func(struc, ha1_num, ha2_num):
continue
elif (self._ignore_backbone and
self._backboneInteraction(struc, ha1_num, ha2_num)):
continue
res1 = ResTuple(struc.atom[atom1])
res2 = ResTuple(struc.atom[atom2])
self.results[res1].interactions[inter_name][res2] += 1
self.results[res2].interactions[inter_name][res1] += 1
count += 1
self.results_overall['OVERALL_' + inter_name.upper()] = count
def _calculateHbonds(self, struc, group_atoms):
"""
Record all hydrogen bonds in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
self._calculateInteraction(self._hbondIterator, Interactions.H_BOND,
struc, group_atoms, self._recordInteraction)
def _hbondIterator(self, struc, group_atoms):
"""
Create an iterator that iterates through all hydrogen bonds between two
groups of atoms
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:return: An iterator that produces tuples of two atoms numbers,
representing (the group 1 atom involved in the hydrogen bond, the
group 2 atom involved in the hydrogen bond)
:rtype: iter
:note: This function differs from
`schrodinger.structutils.analyze.hbond_iterator` in that this
function: uses mm.mmct_hbond_get_user_spec_inter() in place of
mm.mmct_hbond_get_inter() so that non-default values for distance
and angle cutoffs can be used. This function Takes two atom lists,
rather than assuming that the second list of atoms is the converse
of the first.
"""
max_dist = self.interaction_params.hbond_max_dist
min_a_angle = self.interaction_params.hbond_min_a_angle
min_d_angle = self.interaction_params.hbond_min_d_angle
atoms1 = mmlist.mmlist(group_atoms[0])
atoms2 = mmlist.mmlist(group_atoms[1])
return_mmlist = mm.mmct_hbond_get_user_spec_inter(
struc, atoms1, struc, atoms2, max_dist, min_a_angle, min_d_angle)
return_pylist = mmlist._mmlist_to_pylist(return_mmlist)
mm.mmlist_delete(return_mmlist)
for i in range(old_div(len(return_pylist), 4)):
atom1 = return_pylist[i * 4 + 1]
atom2 = return_pylist[i * 4 + 3]
yield (atom1, atom2)
def _calculateClashes(self, struc, group_atoms):
"""
Record all steric clashes in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
self._calculateInteraction(self._clashIterator, Interactions.CLASH,
struc, group_atoms, self._checkInteraction)
def _clashIterator(self, struc, group_atoms):
"""
Create an iterator that iterates through all steric clashes between two
groups of atoms
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:return: An iterator that produces tuples of two atoms numbers,
representing (the group 1 atom involved in the steric clash, the group
2 atom involved in the steric clash)
:rtype: iter
"""
allowable_overlap = self.interaction_params.allowable_vdw_overlap
atoms1 = list(group_atoms[0])
atoms2 = list(group_atoms[1])
clash_iter = interactions.clash_iterator(
struc,
atoms1=atoms1,
atoms2=atoms2,
allowable_overlap=allowable_overlap)
return ((atom1.index, atom2.index) for atom1, atom2, dist in clash_iter)
def _calculateSaltBridges(self, struc, group_atoms):
"""
Record all salt bridges in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
self._calculateInteraction(self._saltBridgeIterator,
Interactions.SALT_BRIDGE, struc, group_atoms,
self._recordInteraction)
def _saltBridgeIterator(self, struc, group_atoms):
"""
Create an iterator that iterates through all salt bridges between two
groups of atoms
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:return: An iterator that produces tuples of two atoms numbers,
representing (the group 1 atom involved in the salt bridge, the group
2 atom involved in the salt bridge)
:rtype: iter
"""
cutoff = self.interaction_params.salt_bridge_max_dist
(group1, group2) = group_atoms
salt_bridges = interactions.get_salt_bridges(
struc,
group1,
group2=group2,
cutoff=cutoff,
order_by=interactions.OrderBy.InputOrder)
return [(a1.index, a2.index) for a1, a2 in salt_bridges]
def _calculateHydrophobicScore(self, struc, group_atoms):
"""
Calculate the hydrophobic interactions between the two groups of atoms.
This hydrophobic score is defined in ChemScore and VSGB2.0.
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
#TODO: save the score for each residue in addition to the total score
nearby_iter = self._nearbyAtomIterator(struc, group_atoms,
HYDROPHOBIC_DIST_CUTOFF)
total_score = 0.0
for (atom1, atom2, dist) in nearby_iter:
if self._isHydrophobicAtom(
struc, atom1) and self._isHydrophobicAtom(struc, atom2):
score = lipophilic_ChemScore_value(struc.atom[atom1],
struc.atom[atom2], dist)
# Use -0.5 kcal/mol energy coefficient
total_score += -0.5 * score
self.results_overall[
'OVERALL_INTERFACE_HYDROPHOBIC_SCORE'] = total_score
def _isHydrophobicAtom(self, struc, iatom):
"""
Check if an atom is a hydrophobic atom
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param iatom: the index of the atom being analyzed
:type iatom: integer
:return: whether or not the atom is defined as hydrophobic atom
:rtype: bool
"""
atom = struc.atom[iatom]
resname = atom.pdbres.strip()
flag = False
if resname in ALL_RESIDUES:
if resname in RESIDUE_HYDROPHOBIC_ATOMS:
if atom.pdbname in RESIDUE_HYDROPHOBIC_ATOMS[resname]:
flag = True
else:
if atom.element in ['C', 'S'] and abs(atom.partial_charge) < 0.25:
flag = True
return flag
def _calculateNearbyRes(self, struc, group_atoms):
"""
Record all neighboring residues in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
inter_residues = set()
inter_residues_group1 = set()
inter_atoms = set()
inter_atoms_group1 = set()
nearby_iter = self._nearbyAtomIterator(
struc, group_atoms, self.interaction_params.neighbor_max_dist)
for (atom1, atom2, dist) in nearby_iter:
res1 = ResTuple(struc.atom[atom1])
res2 = ResTuple(struc.atom[atom2])
if dist < self.results[res1].nearby_res[res2]:
self.results[res1].nearby_res[res2] = dist
self.results[res2].nearby_res[res1] = dist
inter_atoms.update([atom1, atom2])
inter_atoms_group1.update([atom1])
inter_residues.update([res1, res2])
inter_residues_group1.update([res1])
self.results_overall['OVERALL_INTERFACE_ATOMS'] = len(inter_atoms)
self.results_overall['OVERALL_INTERFACE_RESIDUES'] = len(inter_residues)
self.results_overall['OVERALL_INTERFACE_ATOMS_GROUP1'] = len(
inter_atoms_group1)
self.results_overall['OVERALL_INTERFACE_ATOMS_GROUP2'] = len(
inter_atoms) - len(inter_atoms_group1)
self.results_overall['OVERALL_INTERFACE_RESIDUES_GROUP1'] = len(
inter_residues_group1)
self.results_overall['OVERALL_INTERFACE_RESIDUES_GROUP2'] = len(
inter_residues) - len(inter_residues_group1)
np = 0
nnp = 0
nlig = 0
for res in inter_residues:
resname = res.pdbres.upper()
if resname in ALL_RESIDUES:
if resname in NONPOLAR_RESIDUES:
nnp += 1
else:
np += 1
else:
nlig += 1
self.results_overall['OVERALL_INTERFACE_RESIDUES_POLAR'] = np
self.results_overall['OVERALL_INTERFACE_RESIDUES_NONPOLAR'] = nnp
self.results_overall['OVERALL_INTERFACE_RESIDUES_OTHER'] = nlig
def _nearbyAtomIterator(self, struc, group_atoms, dist):
"""
Create an iterator that iterates through all neighboring heavy atoms.
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:param dist: the distance cutoff for neighbors
:type dist: float
:return: An iterator that produces tuples of two atoms numbers,
representing (the group 1 atom, the group 2 atom)
:rtype: iter
"""
dist_cell_iter = measure.dist_cell_iterator(struc, dist, group_atoms[0])
for (atom1_num, neighbors) in dist_cell_iter:
atom1 = struc.atom[atom1_num]
if atom1.element == "H":
continue
for neighbor_num in neighbors:
if neighbor_num in group_atoms[1]:
atom2 = struc.atom[neighbor_num]
if atom2.element == "H":
continue
dist = measure.measure_distance(atom1, atom2)
yield (atom1_num, neighbor_num, dist)
def _calculatePiStacks(self, struc, group_strucs):
"""
Record all pi stacking in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_strucs: A list of [Structure object for group 1, Structure
object for group 2]
:type group_strucs: list
"""
max_stack_dist = self.interaction_params.max_stack_dist
stack_iter = PiStackFinder.createIter(struc, group_strucs,
max_stack_dist)
for (ring1_atoms, ring2_atoms) in stack_iter:
# Record interacting atoms combinatorially
ring1_has = [
self._getHeavyAtomNum(struc, atom) for atom in ring1_atoms
]
ring2_has = [
self._getHeavyAtomNum(struc, atom) for atom in ring2_atoms
]
for ha1 in ring1_has:
for ha2 in ring2_has:
self._recordInteraction(struc, ha1, ha2)
# Record the interaction itself
atom_from_ring1 = ring1_atoms[0]
atom_from_ring2 = ring2_atoms[0]
res1 = ResTuple(struc.atom[atom_from_ring1])
res2 = ResTuple(struc.atom[atom_from_ring2])
self.results[res1].pi_stack[res2] += 1
self.results[res2].pi_stack[res1] += 1
def _calculateDisulfides(self, struc, group_atoms):
"""
Record all disulfide bonds in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
self._calculateInteraction(self._disulfideIterator,
Interactions.DISULFIDE, struc, group_atoms,
self._recordInteraction)
def _disulfideIterator(self, struc, group_atoms):
"""
Create an iterator that iterates through all disulfide bonds
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
:return: An iterator that produces tuples of two atoms numbers,
representing (the group 1 atom, the group 2 atom)
:rtype: iter
"""
for atom1_num in group_atoms[0]:
atom1 = struc.atom[atom1_num]
if not self._potentialDisulfideAtom(atom1):
continue
for bond in atom1.bond:
atom2 = bond.atom2
atom2_num = atom2.index
if (self._potentialDisulfideAtom(atom2) and
atom2_num in group_atoms[1]):
yield (atom1_num, atom2_num)
def _potentialDisulfideAtom(self, atom):
"""
Determine if the given atom is a potential disulfide bond participant
:param atom: The atom object to check
:type atom: `schrodinger.structure._StructureAtom`
:return: True if the given atom is a potential disulfide bond
participant, False otherwise
:rtype: bool
"""
return (atom.element == "S" and atom.pdbres.strip() != "MET")
def _calculateBuriedSASA(self, group_strucs, combined_struc):
"""
Record all buried solvent accessible surface area in the structure
:param group_strucs: A list of [Structure object for group 1, Structure
object for group 2]
:type group_strucs: list
:param combined_struc: A single Structure object containing the protein
atoms in groups 1 and 2
:type combined_struc: `schrodinger.structure.Structure`
"""
all_sasas = []
for cur_struc in [combined_struc] + group_strucs:
sasa_by_atom = analyze.calculate_sasa_by_atom(cur_struc)
sasa_by_res = self._sumSASAByRes(cur_struc, sasa_by_atom)
all_sasas.append(sasa_by_res)
group1_residues = [str(res) for res in list(all_sasas[1])]
group2_residues = [str(res) for res in list(all_sasas[2])]
bound_sasa = all_sasas[0]
unbound_sasa = all_sasas[1]
unbound_sasa.update(all_sasas[2])
for (res, bound) in bound_sasa.items():
unbound = unbound_sasa[res]
if unbound == 0 or approx_eq(bound, unbound):
continue
buried_sasa = 1 - (old_div(bound, unbound))
if round(buried_sasa, 3) != 0:
# Don't include values that are to small to be printed in the
# table, since we don't want blank rows to be created. The
# table displays values with precision 3, hence rounding to 3
# decimal places here.
self.results[res].buried_sasa = buried_sasa
delta = unbound - bound
self.results_overall['OVERALL_BURIED_SASA'] += delta
if str(res) in group1_residues:
self.results_overall['OVERALL_BURIED_SASA_GROUP1'] += delta
else:
self.results_overall['OVERALL_BURIED_SASA_GROUP2'] += delta
resname = res.pdbres.upper()
if resname in ALL_RESIDUES:
if resname in NONPOLAR_RESIDUES:
self.results_overall[
'OVERALL_BURIED_SASA_NONPOLAR'] += delta
else:
self.results_overall['OVERALL_BURIED_SASA_POLAR'] += delta
else:
# ligands and non-standard residues
self.results_overall['OVERALL_BURIED_SASA_OTHER'] += delta
def _sumSASAByRes(self, struc, sasa_by_atom):
"""
Sum solvent accessible surface areas by residue
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param sasa_by_atom: A list of SASAs by atom
:type sasa_by_atom: list
:return: A dictionary of {ResTuple: sasa}
:rtype: dict
"""
sasa_by_res = defaultdict(float)
for atom, sasa in zip(struc.atom, sasa_by_atom):
sasa_by_res[ResTuple(atom)] += sasa
return sasa_by_res
def _recordInteraction(self, struc, ha1_num, ha2_num):
"""
Record an interaction between two heavy atoms so that the two atoms
aren't later counted as clashing
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param ha1_num: The atom number of the heavy atom from group 1
:type ha1_num: int
:param ha2_num: The atom number of the heavy atom from group 2
:type ha2_num: int
:return: Always returns True (for compatibility with _checkInteraction
when used in _calculateInteraction)
:rtype: bool
"""
res1 = ResTuple(struc.atom[ha1_num])
res2 = ResTuple(struc.atom[ha2_num])
self.results[res1].interacting_atoms[ha1_num].add(ha2_num)
self.results[res2].interacting_atoms[ha2_num].add(ha1_num)
return True
def _checkInteraction(self, struc, ha1_num, ha2_num):
"""
Find out if two heavy atoms interact so we know whether to ignore their
clashing. (The interactions must have been previously recorded using
_recordInteraction.)
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param ha1_num: The atom number of the heavy atom from group 1
:type ha1_num: int
:param ha2_num: The atom number of the heavy atom from group 2
:type ha2_num: int
:return: False if the heavy atoms interact, True otherwise
:rtype: bool
"""
res = ResTuple(struc.atom[ha1_num])
return ha2_num not in self.results[res].interacting_atoms[ha1_num]
def _getHeavyAtomNum(self, struc, atom_num):
"""
If the specified atom is a hydrogen, get the bound heavy atom
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param atom_num: The atom number of the atom to examine
:type atom_num: int
:return: If the specified atom is a hydrogen, return the atom number of
the bound heavy atom. Otherwise, returns `atom_num`
:rtype: int
"""
atom = struc.atom[atom_num]
if atom.element != "H":
return atom_num
else:
try:
return atom.bond[1].atom2.index
except IndexError:
# If the hydrogen isn't bound to anything, then just return the
# hydrogen itself instead of its bound heavy atom. (If this
# happens, something funny is going on with the structure.)
print("Warning: Atom %i is a hydrogen but isn't bound to "
"anything" % atom_num)
return atom_num
def _backboneInteraction(self, struc, ha1_num, ha2_num):
"""
Determine if the interaction between the specified atoms is a backbone-
backbone interaction. (i.e. Are both of the specified atoms from the
backbone?)
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param ha1_num: The atom number of the heavy atom from group 1
:type ha1_num: int
:param ha2_num: The atom number of the heavy atom from group 2
:type ha2_num: int
:return: True if both heavy atoms are backbone atome, False otherwise
:rtype: bool
"""
BACKBONE_ATOMS = {'C', 'CA', 'N', 'O', 'OXT'}
ha1_pdbname = struc.atom[ha1_num].pdbname.strip()
ha2_pdbname = struc.atom[ha2_num].pdbname.strip()
return (ha1_pdbname in BACKBONE_ATOMS and ha2_pdbname in BACKBONE_ATOMS)
def _calculateVdwComp(self, struc, group_atoms):
"""
Record all van der Waals surface complementarity in the structure
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_atoms: A list of [atom numbers for group 1, atom numbers
for group 2]
:type group_atoms: list
"""
group1_atoms, group2_atoms = group_atoms
comp_by_atom = surfcomp.calc_complementarity_by_atom(
struc, group1_atoms, group2_atoms)
# Sum up the surface complementarity scores by residue
for cur_res in struc.residue:
# Get all the surface values corresponding to this residue
surface_comps = []
for atom_num in cur_res.getAtomIndices():
surface_comps.extend(comp_by_atom.get(atom_num, []))
if not surface_comps:
# This residue isn't on the interaction surface
continue
# Calculate the median to determine the surface complementarity
# value for the residue (see PYTHON-2433)
res_comp = numpy.median(surface_comps)
if round(res_comp, 2) != 0:
# Don't include values that are to small to be printed in the
# table, since we don't want blank rows to be created. The
# table displays values with precision 2, hence rounding to 2
# decimal places here.
self.results[ResTuple(cur_res)].vdw_comp = res_comp
# Record overall values
res_vdw_comps = [
interactions.vdw_comp
for (res, interactions) in self.results.items()
]
self.results_overall[
'OVERALL_SURF_COMPLENTARITY_RESIDUE_AVERAGE'] = numpy.mean(
res_vdw_comps)
self.results_overall['OVERALL_SURF_COMPLENTARITY'] = surfcomp.calc_total_complementarity(\
struc, group1_atoms, group2_atoms)
[docs]class PiStackFinder(object):
"""
Find pi-pi interactions in proteins
:cvar NON_AROMATIC_RES: A set of residue types that don't contain aromatic
side chains
:vartype NON_AROMATIC_RES: set
"""
# A list of residue types with non-aromatic side chains
NON_AROMATIC_RES = {
'ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'ILE', 'LEU',
'LYS', 'MET', 'PRO', 'SER', 'THR', 'VAL'
}
[docs] def __init__(self, max_stack_dist=MAX_STACK_DIST):
"""
Initialize a new object using the specified interaction cutoffs
:param max_stack_dist: The maximum distance between two ring centroids
allowed for face-face interactions.
:type max_stack_dist: float
"""
self.max_stack_dist = max_stack_dist
[docs] @classmethod
def createIter(cls, struc, group_strucs, max_stack_dist=MAX_STACK_DIST):
"""
A convenience function to initalize the class and return an iterator
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_strucs: A list of [Structure object for group 1, Structure
object for group 2]
:type group_strucs: list
:param max_stack_dist: The maximum distance between two ring centroids
allowed for face-face interactions.
:type max_stack_dist: float
:return: An iterator that produces tuples of two atom number lists,
representing (the ring atoms from group 1 involved in the stacking, the
ring atoms from group 2 involved in the stacking
:rtype: iter
"""
pi_stack_finder = cls(max_stack_dist)
return pi_stack_finder.piStacksIterator(struc, group_strucs)
[docs] def piStacksIterator(self, struc, group_strucs):
"""
Create an iterator that iterates through all pi stacking between two
groups of atoms
:param struc: The structure being analyzed
:type struc: `schrodinger.structure.Structure`
:param group_strucs: A list of [Structure object for group 1, Structure
object for group 2]
:type group_strucs: list
:return: An iterator that produces tuples of two atom number lists,
representing (the ring atoms from group 1 involved in the stacking, the
ring atoms from group 2 involved in the stacking
:rtype: iter
"""
max_stack_dist = self.max_stack_dist
atoms1 = [
atom.property[ORIG_ATOM_NUM_PROP] for atom in group_strucs[0].atom
]
atoms2 = [
atom.property[ORIG_ATOM_NUM_PROP] for atom in group_strucs[1].atom
]
params = infrastructure.PiPiParams()
params.setFaceToFaceMaximumDistance(max_stack_dist)
pipis = interactions.find_pi_pi_interactions(struc, struc, atoms1,
atoms2, params)
for cur_pipi in pipis:
yield (cur_pipi.ring1.atoms, cur_pipi.ring2.atoms)
[docs]class InteractionCalculatorError(Exception):
"""
An error that happens during InteractionCalculator calculations
"""
# This class intentionally left blank