"""
Module to generate cross link restraints by analyzing a trajectory.
Copyright Schrodinger, LLC. All rights reserved.
"""
import copy
import itertools
import math
from typing import List
from typing import Optional
import numpy as np
from scipy.stats import circmean
from scipy.stats import circvar
from scipy.stats import ttest_ind
from schrodinger.application.desmond import cms
from schrodinger.application.desmond.packages import msys
from schrodinger.application.desmond.packages import topo
from schrodinger.application.desmond.packages import traj
from schrodinger.structutils import analyze
from . import utils
__all__ = [
'CrossLinkRestraint', 'CrossLinkRestraintCentroid',
'CrossLinkRestraintInteraction', 'FragmentLinkingRestraint',
'gen_cross_link_restraint', 'gen_cross_link_restraint_centroid',
'gen_cross_link_restraint_interaction'
]
[docs]class CrossLinkRestraint:
[docs] def __init__(self, A=None, B=None, C=None, a=None, b=None, c=None):
"""
`A`, `B`, and `C` are indices of the ligand's atoms, and `a`, `b`, and
`c` are of the receptor's.
"""
self.A = A
self.B = B
self.C = C
self.a = a
self.b = b
self.c = c
# yapf: disable
# The values are (mean, var).
self.Aa = (None, None) # NOQA: E221
self.BAa = (None, None) # NOQA: E221
self.Aab = (None, None) # NOQA: E221
self.BAab = (None, None) # NOQA: E221
self.CBAa = (None, None) # NOQA: E221
self.Aabc = (None, None) # NOQA: E221
# The extra angles are for excluding colinear structures, not part of
# the cross link restraint.
self.CBA = (None, None)
self.abc = (None, None)
# Time series of the corresponding measurements
self.ts_Aa = []
self.ts_BAa = []
self.ts_Aab = []
self.ts_BAab = []
self.ts_CBAa = []
self.ts_Aabc = []
self.ts_CBA = []
self.ts_abc = []
# yapf: enable
def __str__(self):
s = f"CBAabc = {self.C} {self.B} {self.A} {self.a} {self.b} {self.c}\n"
if (None, None) in [
self.Aa, self.BAa, self.Aab, self.BAab, self.CBAa, self.Aabc,
self.CBA, self.abc
]:
s += "incomplete"
else:
s += "Aa = %.5g +- %.5g\n" % (
self.Aa[0],
math.sqrt(self.Aa[1]),
)
s += "BAa = %.5g +- %.5g\n" % (
self.BAa[0],
math.sqrt(self.BAa[1]),
)
s += "Aab = %.5g +- %.5g\n" % (
self.Aab[0],
math.sqrt(self.Aab[1]),
)
s += "BAab = %.5g +- %.5g\n" % (
self.BAab[0],
math.sqrt(self.BAab[1]),
)
s += "Aabc = %.5g +- %.5g\n" % (
self.Aabc[0],
math.sqrt(self.Aabc[1]),
)
s += "CBAa = %.5g +- %.5g\n" % (
self.CBAa[0],
math.sqrt(self.CBAa[1]),
)
s += "CBA = %.5g +- %.5g\n" % (
self.CBA[0],
math.sqrt(self.CBA[1]),
)
s += "abc = %.5g +- %.5g\n" % (
self.abc[0],
math.sqrt(self.abc[1]),
)
return s
[docs] def asMsjSetting(self,
fc,
sigma,
schedule_name,
lambda_state,
schedule_dihed_name=None):
"""
:type fc: list or tuple with 3 elements
:param fc: Force constants of the stretch, the angle, and the dihedral
restraints.
:type sigma: list or tuple with 3 elements
:param sigma: Sigmas of the stretch, the angle, and the dihedral
restraints.
:type schedule_name: str
:param schedule_name: Lambda schedule name of the restraints
:type lambda_state: int
:param lambda_state: Must be either 0 or 1. `0` represent the reference
state, whereas `1` the mutant state.
:param schedule_dihed_name: Name of the dihed schedule or
use `schedule_name` if not set.
:rtype: str
:return: Msj settings of the restraints.
"""
schedule_dihed_name = schedule_dihed_name or schedule_name
# N.B.: This works only for alchemical FEP.
A_state = 1 - lambda_state
B_state = lambda_state
fcA = [e * A_state for e in fc]
fcB = [e * B_state for e in fc]
return f"""
{{name = alchemical_stretch_fbhw
atoms = ["a.n {self.A}" "a.n {self.a}"]
r0A = {self.Aa[0]}
r0B = {self.Aa[0]}
sigmaA = {sigma[0]}
sigmaB = {sigma[0]}
force_constants = [{fcA[0]} {fcB[0]}]
schedule = {schedule_name}
}}
{{name = alchemical_angle_fbhw
atoms = ["a.n {self.B}" "a.n {self.A}" "a.n {self.a}"]
theta0A = {self.BAa[0]}
theta0B = {self.BAa[0]}
sigmaA = {sigma[1]}
sigmaB = {sigma[1]}
force_constants = [{fcA[1]} {fcB[1]}]
schedule = {schedule_name}
}}
{{name = alchemical_angle_fbhw
atoms = ["a.n {self.A}" "a.n {self.a}" "a.n {self.b}"]
theta0A = {self.Aab[0]}
theta0B = {self.Aab[0]}
sigmaA = {sigma[1]}
sigmaB = {sigma[1]}
force_constants = [{fcA[1]} {fcB[1]}]
schedule = {schedule_name}
}}
{{name = alchemical_improper_fbhw
atoms = ["a.n {self.B}" "a.n {self.A}" "a.n {self.a}" "a.n {self.b}"]
phi0A = {self.BAab[0]}
phi0B = {self.BAab[0]}
sigmaA = {sigma[2]}
sigmaB = {sigma[2]}
force_constants = [{fcA[2]} {fcB[2]}]
schedule = {schedule_dihed_name}
}}
{{name = alchemical_improper_fbhw
atoms = ["a.n {self.C}" "a.n {self.B}" "a.n {self.A}" "a.n {self.a}"]
phi0A = {self.CBAa[0]}
phi0B = {self.CBAa[0]}
sigmaA = {sigma[2]}
sigmaB = {sigma[2]}
force_constants = [{fcA[2]} {fcB[2]}]
schedule = {schedule_dihed_name}
}}
{{name = alchemical_improper_fbhw
atoms = ["a.n {self.A}" "a.n {self.a}" "a.n {self.b}" "a.n {self.c}"]
phi0A = {self.Aabc[0]}
phi0B = {self.Aabc[0]}
sigmaA = {sigma[2]}
sigmaB = {sigma[2]}
force_constants = [{fcA[2]} {fcB[2]}]
schedule = {schedule_dihed_name}
}}
"""
@property
def atoms(self):
"""
Returns indices of atoms currently included in this restraint.
"""
return list(
filter(None, [self.A, self.B, self.C, self.a, self.b, self.c]))
@property
def variance(self):
"""
Returns sum of all variances.
"""
# FIXME: Shall we weight them? It seems working fine without any
# weightings. But note these values are NOT of the same units.
tmp = [
self.Aa, self.BAa, self.Aab, self.BAab, self.CBAa, self.Aabc,
self.CBA, self.abc
]
return sum(filter(None, list(zip(*tmp))[1]))
def _calcRefvalues(self):
"""
Calculates the means and variances for the six geometries from the
corresponding time series.
"""
for r in ("Aa", "BAa", "Aab", "CBA", "abc"):
ts = np.asarray(getattr(self, "ts_" + r))
mean, var = (ts.mean(), ts.var()) if len(ts) else (None, None)
setattr(self, r, (mean, var))
for r in ("BAab", "Aabc", "CBAa"):
ts = np.asarray(getattr(self, "ts_" + r))
mean, var = (circmean(ts, -180, 180), circvar(ts, -180, 180)) \
if len(ts) else (None, None)
setattr(self, r, (mean, var))
def _newRestraints(self, atoms: List[int], which_atom: str):
"""
Generates a number of new `CrossLinkRestraint` objects by replacing the
current value of the atom field, as specified by `which_atom`, with each
element of the given `atoms`.
"""
for a in atoms:
new_restraint = copy.deepcopy(self)
setattr(new_restraint, which_atom, a)
yield new_restraint
@staticmethod
def _findAtoms(ct, cross_link_restraint, r_lazy, asl):
"""
Finds atoms that can be considered for the given `cross_link_restraint`,
which is partially completed.
"""
# All selected atoms must be within `r_lazy` distance to each other.
for i in cross_link_restraint.atoms:
asl = "(%s) and (within %f (a.n %d))" % (asl, r_lazy, i)
# Excludes the atoms already selected for the cross-link restraint.
restraint_atom_str = " ".join(map(str, cross_link_restraint.atoms))
asl = "(%s) and (not (a.n %s))" % (asl, restraint_atom_str)
return analyze.evaluate_asl(ct, asl)
@staticmethod
def _violateMidpoint(ct, atoms, r_clone2):
"""
"Midpoint" is the method used by Desmond to calculate the neighbor
lists. Without going too much technical detail, here we ensure the
following condition:
1. Calculate the centroid of the given atoms (`atoms`) in the structure
(`ct`).
2. All given atoms must be within a sphere that is centered at the
centroid and has a radius of `sqrt(r_clone2)`.
:type atoms: list[Union(int, None)]
:type r_clone2: float
:rtype: bool
"""
atoms_xyz = ct.getXYZ(copy=False)[[i - 1 for i in atoms]]
centroid = atoms_xyz.mean(axis=0)
for xyz in atoms_xyz:
dx, dy, dz = xyz - centroid
d2 = dx * dx + dy * dy + dz * dz
if (d2 > r_clone2):
return True
return False
@staticmethod
def _distilCandidates(candidates,
model: cms.Cms,
tr: List[traj.Frame],
r_clone2,
min_angle=10.0) -> List['CrossLinkRestraint']:
"""
Removes candidates that violate the following rules:
1. Satisfy the midpoint rule (see `_violateMidpoint` method above for
detail), in any of the given structures (`cts`).
2. Calculate the average reference values of the six geometries. The
reference values of BAa and Aab angles should NOT be close to either
0 or 180 degrees.
Returns the surviving candidates.
:type candidates: `list[CrossLinkRestraint]`
:type model: cms.Cms
:param model: Simulation system
:type tr: list[traj.Frame]
:param tr: The trajectory to examine
:type r_clone2: `float`
:param r_clone2: = r_clone * r_clone. (See `findBest` docstring for
the definition of "r_clone")
:type min_angle: `float`
:param min_angle: Mininum angle (in degrees) away from 0 or 180 degrees.
"""
# For each restraint candidate, calculates the time series for the six
# geometries. Eliminates the candidates that violate the midpoint rule.
ct = model.fsys_ct.copy()
# Read in the frames only when needed to reduce memory usage
for fr in tr:
topo.update_fsys_ct_from_frame(ct, model, fr)
for c in list(candidates):
if CrossLinkRestraint._violateMidpoint(ct, c.atoms, r_clone2):
candidates.remove(c)
else:
if c.Aa == (None, None):
c.ts_Aa.append(ct.measure(c.A, c.a))
if c.B is not None:
if c.BAa == (None, None):
c.ts_BAa.append(ct.measure(c.B, c.A, c.a))
if c.C is not None:
if c.CBAa == (None, None):
c.ts_CBAa.append(ct.measure(c.C, c.B, c.A, c.a))
if c.CBA == (None, None):
c.ts_CBA.append(ct.measure(c.C, c.B, c.A))
if c.b is not None:
if c.Aab == (None, None):
c.ts_Aab.append(ct.measure(c.A, c.a, c.b))
if c.c is not None:
if c.Aabc == (None, None):
c.ts_Aabc.append(ct.measure(c.A, c.a, c.b, c.c))
if c.abc == (None, None):
c.ts_abc.append(ct.measure(c.a, c.b, c.c))
if c.B is not None and c.BAab == (None, None):
c.ts_BAab.append(ct.measure(c.B, c.A, c.a, c.b))
# For the surviving candidates, calculates the reference values of the
# six geometries. Eliminates the candidates that have angles close to 0
# or 180 degrees.
for c in list(candidates):
c._calcRefvalues()
for angle in c._angles():
if (angle is not None and
(angle < min_angle or angle > 180.0 - min_angle)):
candidates.remove(c)
break
return candidates
def _angles(self):
return (self.BAa[0], self.Aab[0], self.CBA[0], self.abc[0])
[docs] @staticmethod
def findBest(model,
tr,
*,
ligand_asl,
receptor_asl,
r_clone=4,
min_angle=10.0,
verbose=False,
use_bonded_atoms=False) -> Optional['CrossLinkRestraint']:
"""
Examine the trajectory, and find and return the "best" cross-link
restraint between the ligand molecule and the receptor molecule. Three
atoms from the ligand molecule and three atoms from the receptor
molecule will be selected. We denote the three ligand atoms as A, B,
and C, the three receptor atoms as a, b, and c. Relative restraints
will be applied to the following geometries:
- 1 stretch restraint: A-a distance
- 2 angle restraints: BAa and Aab angles
- 3 dihedral restraints: CBAa, BAab, Aabc dihedral angles
Only the given ligand atoms (`lig_atoms`) will be considered for the
ligand atoms of the restraint, whereas all "protein" atoms will be
considered for the receptor atoms of the restraint. We use the ASL
expression "protein" to find protein atoms. It's OK if the `lig_atoms`
is part of the result of that expression.
Caveats:
- This won't work if the receptor is NOT a protein.
- We assumed that there are at least 3 atoms in the ligand molecule
and 3 atoms in the receptor molecule(s).
The "best" criteria are the following:
1. To satisfy the requirement of Desmond's midpoint algorithm, all
atoms are close to each other (i.e., within a sphere of radius of
`r_clone`). This is mainly for the efficiency of the Desmond
backend.
2. No ill geometries such as colinear structures
3. The variance of the restraint is the least.
:type model: cms.Cms
:param model: Simulation system
:type tr: list[traj.Frame]
:param tr: The trajectory to examine
:type ligand_asl: str
:parma ligand_asl: ASL expression to specify candidate ligand atoms for
the restraint
:type receptor_asl: str
:parma receptor_asl: ASL expression to specify candidate receptor atoms
for the restraint
:type r_clone: float
:param r_clone: Radius of particle / home box visibility (Desmond's
definition). Its value is half of the real space cutoff distance.
We use 4 (Angstroms) as a safe default, assuming the cutoff
distance is 8 which is a bit less than the typical value. The
exact value doesn't matter, but it's better to be less as opposed
(to greater) than the actual r_clone used in the simulation.
:type min_angle: `float`
:param min_angle: Mininum angle (in degrees) away from 0 or 180 degrees.
:param use_bonded_atoms: Set to True to choose 'ABC' and 'abc' atoms
that are bonded together and form an angle. The order of the atoms
does not matter, just that they are bonded together.
Default of False will not place restrictions on how these
atoms are connected.
:rtype: `CrossLinkRestraint`
"""
r_clone2 = r_clone * r_clone
r_lazy = r_clone * 2
lig_atoms = model.select_atom(ligand_asl)
if use_bonded_atoms:
candidates = []
# Iterate three ligand atoms are bonded together in any order
for A, B, C in analyze.angle_iterator(model.fsys_ct, lig_atoms):
receptor_atoms = model.select_atom(
f'{receptor_asl} and (within {r_lazy} (a.n {A} {B} {C}))')
if len(receptor_atoms) == 0:
continue
# Iterate over three receptor atoms are bonded together in any order
for a, b, c in analyze.angle_iterator(model.fsys_ct,
receptor_atoms):
# Consider the equivalent angles starting from either terminal
# atom. The ordering can affect the variance.
for A_, B_, C_ in ((A, B, C), (C, B, A)):
for a_, b_, c_ in ((a, b, c), (c, b, a)):
candidates.append(
CrossLinkRestraint(A=A_,
B=B_,
C=C_,
a=a_,
b=b_,
c=c_))
# Distil and sort restraint candiates
candidates = CrossLinkRestraint._distilCandidates(
candidates, model, tr, r_clone2, min_angle)
candidates.sort(key=lambda r: r.variance)
else:
candidates = [CrossLinkRestraint(A=atom) for atom in lig_atoms]
# We find "A" and "a", then "B" and "b", and then "C" and "c".
verbose and print("%d candidates for 'A'..." % len(candidates))
for which_atom in ["a", "B", "b", "C", "c"]:
new_candidates = []
for c in candidates:
atoms = CrossLinkRestraint._findAtoms(
model, c, r_lazy, receptor_asl if
(which_atom > "Z") else ligand_asl)
new_candidates += c._newRestraints(atoms, which_atom)
candidates = CrossLinkRestraint._distilCandidates(
new_candidates, model, tr, r_clone2, min_angle)
verbose and print("%d candidates for '%s'..." %
(len(candidates), which_atom))
# We only look further at the currently best 128 candidates, which
# should be more than enough. Larger numbers will expand the search
# but at the same time slow it down.
candidates.sort(key=lambda r: r.variance)
candidates = candidates[:128]
if not candidates:
return None
return candidates[0]
[docs]class CrossLinkRestraintCentroid(CrossLinkRestraint):
"""
Generate restraints using the centroid method.
See `findBest` for more details.
"""
@property
def variance(self):
"""
Returns sum of all variances.
"""
tmp = [
self.Aa,
self.BAa,
self.Aab,
self.BAab,
self.CBAa,
self.Aabc,
]
return sum(filter(None, list(zip(*tmp))[1]))
def _angles(self):
return (self.Aab[0], self.BAa[0])
[docs] @staticmethod
def findBest(cms_model: cms.Cms,
tr: List["traj.TrajFrame"],
*,
ligand_asl: str,
receptor_asl: str,
r_clone=9999,
min_angle=45.0) -> Optional['CrossLinkRestraint']:
"""
Generate restraints using the centroid of the binding pocket
and the ligand as the basis for the restraint candidates.
These candidates are then filtered to pick the one
that results in the minimal variance for the restraint
terms over the input trajectory.
See `CrossLinkRestraint.findBest` for more details.
"""
# ASLs
ligand_heavy_asl = utils._get_heavy_asl(ligand_asl)
receptor_ca_asl = f'({receptor_asl}) and a.ptype CA and within 20 ({ligand_asl})' # region to define two axes
binding_pocket_asl = f'({receptor_asl}) and a.ptype CA and within 10 ({ligand_asl})'
ct = cms_model.fsys_ct
# Check for case where ligand is far from protein (should be rare)
if not analyze.evaluate_asl(ct, receptor_ca_asl):
return None
if not analyze.evaluate_asl(ct, binding_pocket_asl):
return None
# Find centroid of binding pocket
bp_centroid_aid = utils._find_centroid_aid(ct, binding_pocket_asl)
candidate = None
min_var = 9999
# Find the ligand centroid
lig_centroid_aids = utils._find_aids_within_cutoff_of_centroid(
ct, ligand_heavy_asl)
for lig_centroid_aid in lig_centroid_aids:
# Find the ligand atom at the farthest end to the centroid
lig_end_aid = utils._find_max_dist_aid(ct, ligand_heavy_asl,
lig_centroid_aid)
# Find the centroid of the first group of protein atoms for axis 1
group1_aids = utils._find_axis_aids(ct, receptor_ca_asl,
lig_end_aid, lig_centroid_aid,
bp_centroid_aid)
if not group1_aids:
continue
group1_asl = f'''a.n {' '.join(map(str, group1_aids))}'''
group1_centroid_aid = utils._find_centroid_aid(ct, group1_asl)
# Find the centroid of the second group of protein atoms for axis 2
group2_aids = utils._find_axis_aids(ct, receptor_ca_asl,
group1_centroid_aid,
lig_centroid_aid,
bp_centroid_aid)
if not group2_aids:
continue
group2_asl = f'''a.n {' '.join(map(str, group2_aids))}'''
group2_centroid_aid = utils._find_centroid_aid(ct, group2_asl)
# A B C ligand atoms
# a b c protein atoms
A = lig_centroid_aid
a = bp_centroid_aid
b = group1_centroid_aid
c = group2_centroid_aid
# Calculate distance, angle and torsion from trajectory
for lig_aid2, lig_aid3 in utils._get_two_atoms(cms_model, A):
for C, B in itertools.permutations([lig_aid2, lig_aid3], 2):
cur_candidate = CrossLinkRestraintCentroid(A=A,
B=B,
C=C,
a=a,
b=b,
c=c)
cur_candidate = CrossLinkRestraintCentroid._distilCandidates(
[cur_candidate], cms_model, tr, r_clone, min_angle)
if cur_candidate and cur_candidate[0].variance < min_var:
min_var = cur_candidate[0].variance
candidate = cur_candidate[0]
return candidate
[docs]class CrossLinkRestraintInteraction(CrossLinkRestraint):
"""
Generate restraints using the interactions between
the receptor and ligand. See `findBest` for more details.
"""
@property
def variance(self):
"""
Returns sum of all variances.
"""
tmp = [
self.Aa,
self.BAa,
self.Aab,
self.BAab,
self.CBAa,
self.Aabc,
]
return sum(filter(None, list(zip(*tmp))[1]))
def _angles(self):
return (self.Aab[0], self.BAa[0])
[docs] @staticmethod
def findBest(
msys_model: "msys.System", # noqa: F821
cms_model: cms.Cms,
tr: List["traj.TrajFrame"],
*,
ligand_asl: str,
receptor_asl: str,
min_angle=45.0,
freq_cutoff=0.5,
max_variance=9999,
max_dist_variance=4.0**2,
num_traj_segments=3) -> Optional['CrossLinkRestraint']:
"""
Uses the protein-ligand interactions to determine the cross link restraints.
These interactions include hydrogen bonds and salt bridge interactions.
The candidates are considered by picking the ones with the
most consistent interactions throughout the input trajectory.
Then the candidates are filtered based on the sum of the
variance for all terms, and the A-a distance variance.
Then the candidates are filtered to pick the one
that results in the minimal overall variance.
See `CrossLinkRestraint.findBest` for more details.
:param freq_cutoff: Interactions must be present at least
this fraction of the trajectory to be considered for
the restraint. The default is 0.5.
:param max_variance: The maximum variance to allow the interaction
to be used as the restraint. The default is 9999 (mixed units).
:param max_dist_variance: The maximum variance to allow for the
A-a distance. Default is 4.0**2 (in Angstrom**2).
:param num_traj_segments: Split the trajectory into `num_traj_segments`
segments for the analysis.
"""
# Set to a large value to disable the r_clone check
r_clone = 9999.0
# Interaction frequencies for each segment of the trajectory
freq = utils._get_protein_ligand_interaction_freq_dict(
msys_model, cms_model, tr, ligand_asl, receptor_asl,
num_traj_segments)
if not freq:
print("No interactions found.")
return None
# Get centroid of the ligand heavy atoms
ct = cms_model.fsys_ct
lig_centroid_data = utils._get_centroid_data(
ct, utils._get_heavy_asl(ligand_asl))
# Loop through all interactions, try the one with the
# smallest distance to the ligand heavy atom centroid.
# If the variance for the selected interaction is too large,
# remove it from consideration and try with the next interaction.
found_candidate = False
while not found_candidate:
selected_freq = freq_cutoff
selected_key = None
selected_data = None
for k in freq:
data = freq[k]
avg_freq = np.average(data)
if avg_freq >= selected_freq:
selected_freq = avg_freq
selected_key = k
selected_data = data
if selected_key is None:
print("No interactions found.")
return None
selected_interactions = [selected_key]
# compare interaction frequency of other sites
# with the site having largest interaction frequency
# based on the t-test. if the difference is not
# significant compared with largest interaction
# frequency site, the site will also be selected.
for k in freq:
if k != selected_key:
_, pval = ttest_ind(selected_data, freq[k])
# pval is nan if both have the same exact values
if math.isnan(pval) or pval > 0.05:
selected_interactions.append(k)
dist = None
selected_freq = None
selected_key = None
for k in selected_interactions:
r = np.linalg.norm(ct.atom[k.ligand_aid].xyz -
lig_centroid_data.centroid)
if dist is None or r < dist:
dist = r
selected_freq = np.average(freq[k])
selected_key = k
if dist is None:
print("No strong interactions found.")
return None
candidate = None
cur_min_variance = max_variance
lig_aid1 = selected_key.ligand_aid
# Loop through the ligand atoms and the protein atoms corresponding
# to the chosen interaction. Check that the variance is below
# the cut-offs and select as a candidate. Repeat for all atom
# permutations, picking the one that results in the lowest overall variance.
for lig_aid2, lig_aid3 in utils._get_two_atoms(cms_model, lig_aid1):
# The SO2 group can easily flip leading to problems, so this is explicitly
# checked for and excluded here.
if not utils._is_so2(cms_model, lig_aid1) and not utils._is_so2(
cms_model, lig_aid2) and not utils._is_so2(
cms_model, lig_aid3):
for C, B, A in itertools.permutations(
[lig_aid1, lig_aid2, lig_aid3], 3):
for a, b, c in itertools.permutations(
selected_key.receptor_aids, 3):
cur_candidate = CrossLinkRestraintInteraction(A=A,
B=B,
C=C,
a=a,
b=b,
c=c)
cur_candidate = CrossLinkRestraintInteraction._distilCandidates(
[cur_candidate], cms_model, tr, r_clone,
min_angle)
if cur_candidate and cur_candidate[0].variance < cur_min_variance and \
cur_candidate[0].Aa[1] < max_dist_variance:
candidate = cur_candidate[0]
cur_min_variance = candidate.variance
found_candidate = True
if found_candidate:
print(f'Selected interaction frequency: {selected_freq}')
print(f'Selected interaction: {selected_key}')
print(f'Selected candidate: {candidate}')
print(f'Selected distance: {dist}')
print(f'Selected variance sum: {candidate.variance}')
print(f'Selected r0 deviation: {candidate.Aa[1]}')
return candidate
# Remove this interaction from consideration
# and try again.
del freq[selected_key]
if not freq:
# Nothing left
print('No good interactions after filtering.')
return None
[docs]class FragmentLinkingRestraint:
[docs] def __init__(self, A: int, B: int, a: int, b: int, Aa: float, BAa: float,
Aab: float):
"""
Container for the stretch restraint Aa and the angle
restraints BAa and Aab. The parameters `A` `B` `a` and `b`
are the atom indicies for the restraints.
:param Aa: Value for the stretch restraint in Angstrom.
:param BAa: Value for the first angle restraint in degrees.
:param Aab: Value for the second angle restraint in degrees.
"""
self.A = A
self.B = B
self.a = a
self.b = b
# Restraint parameters.
self.Aa = Aa
self.BAa = BAa
self.Aab = Aab
def __str__(self):
return f'BAab = {self.B} {self.A} {self.a} {self.b}'
[docs] def asMsjSetting(self, fc, sigma, schedule_name, lambda_state,
atom_asl_dict):
"""
:type fc: list or tuple with 2 elements
:param fc: Force constants of the stretch and the angle restraints.
:type sigma: list or tuple with 2 elements
:param sigma: Sigmas of the stretch and the angle restraints.
:type schedule_name: str
:param schedule_name: Lambda schedule name of the restraints
:type lambda_state: int
:param lambda_state: Must be 0, 1 or None. `0` means to apply the
restraints to the reference state, whereas `1` means to apply
the restraints to the mutant state. `None` means to use regular
nonalchemical restraints, which are fixed regardless of lambda.
:type atom_asl_dict: dict
:param atom_asl_dict: If not None, specify the atom names
'A', 'B', 'a', 'b', as the keys and the corresponding restraint asl
as the values. Default of None means to use the atom numbers as
the asl.
:rtype: str
:return: Msj settings of the restraints.
"""
if atom_asl_dict is None:
atom_asl_dict = {
'A': f"a.n {self.A}",
'B': f"a.n {self.B}",
'a': f"a.n {self.a}",
'b': f"a.n {self.b}",
}
if lambda_state is None:
return f"""{{name = stretch_fbhw
atoms = ["{atom_asl_dict['A']}" "{atom_asl_dict['a']}"]
r0 = {self.Aa}
sigma = {sigma[0]}
force_constants = [{fc[0]}]
}}
{{name = angle_fbhw
atoms = ["{atom_asl_dict['B']}" "{atom_asl_dict['A']}" "{atom_asl_dict['a']}"]
theta0 = {self.BAa}
sigma = {sigma[1]}
force_constants = [{fc[1]}]
}}
{{name = angle_fbhw
atoms = ["{atom_asl_dict['A']}" "{atom_asl_dict['a']}" "{atom_asl_dict['b']}"]
theta0 = {self.Aab}
sigma = {sigma[1]}
force_constants = [{fc[1]}]
}}
"""
else:
A_state = 1 - lambda_state
B_state = lambda_state
fcA = [e * A_state for e in fc]
fcB = [e * B_state for e in fc]
return f"""{{name = alchemical_stretch_fbhw
atoms = ["{atom_asl_dict['A']}" "{atom_asl_dict['a']}"]
r0A = {self.Aa}
r0B = {self.Aa}
sigmaA = {sigma[0]}
sigmaB = {sigma[0]}
force_constants = [{fcA[0]} {fcB[0]}]
schedule = {schedule_name}
}}
{{name = alchemical_angle_fbhw
atoms = ["{atom_asl_dict['B']}" "{atom_asl_dict['A']}" "{atom_asl_dict['a']}"]
theta0A = {self.BAa}
theta0B = {self.BAa}
sigmaA = {sigma[1]}
sigmaB = {sigma[1]}
force_constants = [{fcA[1]} {fcB[1]}]
schedule = {schedule_name}
}}
{{name = alchemical_angle_fbhw
atoms = ["{atom_asl_dict['A']}" "{atom_asl_dict['a']}" "{atom_asl_dict['b']}"]
theta0A = {self.Aab}
theta0B = {self.Aab}
sigmaA = {sigma[1]}
sigmaB = {sigma[1]}
force_constants = [{fcA[1]} {fcB[1]}]
schedule = {schedule_name}
}}
"""
[docs] @classmethod
def findRestraint(cls, ct, fragment0_asl,
fragment1_asl) -> "FragmentLinkingRestraint":
"""
Given a structure and the fragment0_asl/fragment1_asl asls,
find the restraint that
1) minimizes the distance for the stretch term between fragments
2) maximizes the distance between the stretch term
atom, and the atom used for the angle term within each fragment.
"""
frag0_atom_idxs = analyze.evaluate_asl(ct, fragment0_asl)
frag1_atom_idxs = analyze.evaluate_asl(ct, fragment1_asl)
atom_coords = ct.getXYZ()
def distance_squared(atom0, atom1):
# atom_coords are 0 indexed, atoms are 1 indexed
delta = atom_coords[atom0 - 1, :] - atom_coords[atom1 - 1, :]
return np.inner(delta, delta)
# For the stretch term, find the minimal
# distance between the matched atoms
A, a = min(itertools.product(frag0_atom_idxs, frag1_atom_idxs),
key=lambda atom_idxs: distance_squared(*atom_idxs))
Aa = math.sqrt(distance_squared(A, a))
# For the angle terms, find the atom of the same fragment
# that is the farthest away.
def find_maximum_distance_index(ref_atom_idx, frag_atom_idxs) -> int:
return max(
frag_atom_idxs,
key=lambda atom_idx: distance_squared(ref_atom_idx, atom_idx))
B = find_maximum_distance_index(A, frag0_atom_idxs)
b = find_maximum_distance_index(a, frag1_atom_idxs)
BAa = ct.measure(B, A, a)
Aab = ct.measure(A, a, b)
return cls(A=A, B=B, a=a, b=b, Aa=Aa, BAa=BAa, Aab=Aab)
# FIXME: The `stage.SystemBuilder` class is dying. All still useful logic for
# restraint settings should be scavanged into this module.
# -YW, 12/17/2018
[docs]def gen_cross_link_restraint(model,
tr,
*,
ligand_asl,
receptor_asl,
r_clone=4,
min_angle=10.0,
use_bonded_atoms=False) -> CrossLinkRestraint:
"""
A cross-link restraint comprises one stretch restraint, two angle
restraints, and three dihedral restraints between two molecules, which
typically are the ligand and the receptor. A cross link restraint will
completely restrain the relative distance and orientation of the two
molecules.
N.B.: A cross link restraint requires to identify three atoms from the
receptor (denoted as a, b, and c), and three from the ligand (denoted as A,
B, and C). If any of the two molecules are too small to not even have three
non-terminal atoms, this function will NOT work, and a `RuntimeError`
exception will be raised. Also the ligand atoms and the receptor atoms
should NOT have overlaps, otherwise a `RuntimeError` exception will be
raised.
:type model: `cms.Cms`
:type tr: `list` of `traj.Frame`
:param tr: A MD simulation trajectory that we will analyze to find out the
optimal 6 atoms and the equilibrium values for the cross-link restraint.
:type ligand_asl: `str`
:parma ligand_asl: ASL expression to specify candidate ligand atoms for the
restraint
:type receptor_asl: `str`
:parma receptor_asl: ASL expression to specify candidate receptor atoms for
the restraint
:type r_clone: `int`
:type min_angle: `float`
:param use_bonded_atoms: See `CrossLinkRestraint.findBest` for information on
this parameter.
:rtype: `CrossLinkRestraint`
"""
utils._check_asls(model, ligand_asl, receptor_asl)
best_restraint = CrossLinkRestraint.findBest(
model,
tr,
ligand_asl=ligand_asl,
receptor_asl=receptor_asl,
r_clone=r_clone,
min_angle=min_angle,
use_bonded_atoms=use_bonded_atoms)
utils._check_restraint(best_restraint)
return best_restraint
[docs]def gen_cross_link_restraint_centroid(
model: cms.Cms,
tr: List["traj.TrajFrame"],
*,
ligand_asl: str,
receptor_asl: str,
min_angle=45.0) -> CrossLinkRestraintCentroid:
"""
Generate cross link restraints using the centroid method.
See `gen_cross_link_restraint` for more detail.
"""
utils._check_asls(model, ligand_asl, receptor_asl)
best_restraint = CrossLinkRestraintCentroid.findBest(
model,
tr,
ligand_asl=ligand_asl,
receptor_asl=receptor_asl,
min_angle=min_angle)
utils._check_restraint(best_restraint)
return best_restraint
[docs]def gen_cross_link_restraint_interaction(
msys_model: "msys.System", # noqa: F821
model: cms.Cms,
tr: List["traj.TrajFrame"],
*,
ligand_asl: str,
receptor_asl: str,
min_angle=45.0) -> CrossLinkRestraintInteraction:
"""
Generate cross link restraints looking at hydrogen bonds
and salt bridges to determine the atoms to use for the restraint.
See `gen_cross_link_restraint` for more detail.
:param msys_model: Msys model used to analyze the trajectory.
:raise RuntimeError: If the input asls or structures are invalid.
:raise CrossLinkGenerationError: If a suitable restraint could
not be generated.
"""
utils._check_asls(model, ligand_asl, receptor_asl)
best_restraint = CrossLinkRestraintInteraction.findBest(
msys_model,
model,
tr,
ligand_asl=ligand_asl,
receptor_asl=receptor_asl,
min_angle=min_angle)
utils._check_restraint(best_restraint)
return best_restraint