Commit 063873b2 authored by jdurrant's avatar jdurrant
Browse files

Before removing cruft from old hydrogen/halogen-bond code.

parent 6c2f57df
...@@ -137,8 +137,7 @@ def dihedral(point1, point2, point3, point4): # never tested ...@@ -137,8 +137,7 @@ def dihedral(point1, point2, point3, point4): # never tested
b1Xb2 = cross_product(b1, b2) b1Xb2 = cross_product(b1, b2)
b1XMagb2 = vector_scalar_multiply(b1, b2.magnitude()) b1XMagb2 = vector_scalar_multiply(b1, b2.magnitude())
radians = math.atan2(dot_product(b1XMagb2, b2Xb3), dot_product(b1Xb2, b2Xb3)) return math.atan2(dot_product(b1XMagb2, b2Xb3), dot_product(b1Xb2, b2Xb3))
return radians
def angle_between_three_points(point1, point2, point3): # As in three connected atoms def angle_between_three_points(point1, point2, point3): # As in three connected atoms
...@@ -151,10 +150,8 @@ def angle_between_points(point1, point2): ...@@ -151,10 +150,8 @@ def angle_between_points(point1, point2):
new_point1 = return_normalized_vector(point1) new_point1 = return_normalized_vector(point1)
new_point2 = return_normalized_vector(point2) new_point2 = return_normalized_vector(point2)
dot_prod = dot_product(new_point1, new_point2) dot_prod = dot_product(new_point1, new_point2)
if dot_prod > 1.0: dot_prod = min(dot_prod, 1.0)
dot_prod = 1.0 # to prevent errors that can rarely occur dot_prod = max(dot_prod, -1.0)
if dot_prod < -1.0:
dot_prod = -1.0
return math.acos(dot_prod) return math.acos(dot_prod)
......
...@@ -31,17 +31,181 @@ from binana._utils.shim import fabs ...@@ -31,17 +31,181 @@ from binana._utils.shim import fabs
# Be sure to update the corresponding function in # Be sure to update the corresponding function in
# binana.interactions.__init__.py as well! # binana.interactions.__init__.py as well!
# O-H distance is 0.96 A, N-H is 1.01 A. See # # O-H distance is 0.96 A, N-H is 1.01 A. See
# http://www.science.uwaterloo.ca/~cchieh/cact/c120/bondel.html # # http://www.science.uwaterloo.ca/~cchieh/cact/c120/bondel.html
max_donor_X_dist = { # _max_donor_X_dist = {
"H": 1.3, # "H": 1.3,
"I": 2.04 * 1.4, # O-I: 2.04 per avogadro # "I": 2.04 * 1.4, # O-I: 2.04 per avogadro
"BR": 1.86 * 1.4, # O-Br: 1.86 # "BR": 1.86 * 1.4, # O-Br: 1.86
"Br": 1.86 * 1.4, # "Br": 1.86 * 1.4,
"CL": 1.71 * 1.4, # O-Cl: 1.71 # "CL": 1.71 * 1.4, # O-Cl: 1.71
"Cl": 1.71 * 1.4, # "Cl": 1.71 * 1.4,
"F": 1.33 * 1.4, # O-F: 1.33 # "F": 1.33 * 1.4, # O-F: 1.33
} # }
# # def _mimic_set_minus(set1, set2):
# # for s1 in set1:
# # while s1 in set2:
# # set2.remove(s1)
# # return set2
# def _detect_hydrogen_bonds_without_hydrogens(close_donors_acceptors, ligand, receptor):
# # In case the user has not added hydrogen bonds to their structures. Not as
# # accurate.
# hbonds = {}
# pdb_hbonds = Mol()
# hbonds_labels = []
# for lig_atm, recep_atm in close_donors_acceptors:
# lig_categories = ligand._categorize_donor_acceptor_without_hydrogens(lig_atm)
# recep_categories = receptor._categorize_donor_acceptor_without_hydrogens(
# recep_atm
# )
# if len(lig_categories) < len(recep_categories):
# recep_categories = _mimic_set_minus(lig_categories, recep_categories)
# # Note that transcrypt doesn't support set subtraction.
# # recep_categories = recep_categories - lig_categories
# if len(recep_categories) < len(lig_categories):
# lig_categories = _mimic_set_minus(recep_categories, lig_categories)
# # lig_categories = lig_categories - recep_categories
# # If below occurs, there can't be donor-acceptor pair
# if len(lig_categories) == 0:
# continue
# if len(recep_categories) == 0:
# continue
# if len(lig_categories) > 0 and len(recep_categories) > 0:
# # They could be complementary (donor-acceptor), but they could also
# # both be donor or both be acceptor.
# # Pick first lig one (in case there are many)
# lig_category = list(lig_categories)[0]
# # Remove that one from recep options.
# recep_categories = _mimic_set_minus([lig_category], recep_categories)
# if len(recep_categories) == 0:
# # There are no complementary donor-acceptor pairs.
# continue
# # Pick first of remaining receptor ones
# recep_category = list(recep_categories)[0]
# # *****
# comment = "RECEPTOR" if lig_category == "ACCEPTOR" else "LIGAND"
# hbonds_key = (
# "HDONOR_"
# + comment
# + "_"
# + recep_atm.side_chain_or_backbone()
# + "_"
# + recep_atm.structure
# )
# pdb_hbonds.add_new_atom(lig_atm.copy_of())
# pdb_hbonds.add_new_atom(recep_atm.copy_of())
# hashtable_entry_add_one(hbonds, hbonds_key)
# hbonds_labels.append(
# (
# lig_atm.string_id(),
# # There is no middle hydrogen in this case, so just repeat
# # one of the heteroatoms.
# lig_atm.string_id()
# if comment == "LIGAND"
# else recep_atm.string_id(),
# recep_atm.string_id(),
# comment,
# )
# )
# return hbonds, pdb_hbonds, hbonds_labels
# def _detect_hbonds_by_angle(
# close_donors_acceptors, ligand, receptor, angle_cutoff, hydrogen_bond
# ):
# # Use this function if detecting halogen bonds or hydrogen bonds when the
# # models both include hydrogen atoms.
# hbonds = {}
# pdb_hbonds = Mol()
# hbonds_labels = []
# central_atom_names = ["H"] if hydrogen_bond else ["I", "BR", "Br", "CL", "Cl", "F"]
# # Now see if there's some sort of hydrogen (or halogen) bond between these
# # two atoms. default distance cutoff = 4, angle cutoff = 40. Note that these
# # cutoffs are generous.
# for lig_atm, recep_atm in close_donors_acceptors:
# # now build a list of all the hydrogens (or halogens) close to these
# # atoms
# h_or_hal = []
# for mol, hetero_atom, lbl in [
# [ligand, lig_atm, "LIGAND"],
# [receptor, recep_atm, "RECEPTOR"],
# ]:
# for atm_index in mol.all_atoms.keys():
# central_atom = mol.all_atoms[atm_index]
# element = central_atom.element
# if element in central_atom_names:
# # so it's a hydrogen (or halogen)
# dist = central_atom.coordinates.dist_to(hetero_atom.coordinates)
# if dist < _max_donor_X_dist[element]: # 1.3 for H
# central_atom.comment = lbl
# h_or_hal.append(central_atom)
# # now we need to check the angles
# for center_atm in h_or_hal:
# angle = angle_between_three_points(
# lig_atm.coordinates,
# center_atm.coordinates,
# recep_atm.coordinates,
# )
# if fabs(180 - angle * 180.0 / math.pi) > angle_cutoff:
# # Angle is too big.
# continue
# # If you get here, angle is ok.
# hbonds_key = (
# "HDONOR_"
# + center_atm.comment
# + "_"
# + recep_atm.side_chain_or_backbone()
# + "_"
# + recep_atm.structure
# )
# pdb_hbonds.add_new_atom(lig_atm.copy_of())
# pdb_hbonds.add_new_atom(center_atm.copy_of())
# pdb_hbonds.add_new_atom(recep_atm.copy_of())
# hashtable_entry_add_one(hbonds, hbonds_key)
# hbonds_labels.append(
# (
# lig_atm.string_id(),
# center_atm.string_id(),
# recep_atm.string_id(),
# center_atm.comment,
# )
# )
# return hbonds, pdb_hbonds, hbonds_labels
def _mimic_set_minus(set1, set2):
for s2 in set2:
set1 = [s for s in set1 if s[0] != s2[0]]
return set1
def get_hydrogen_or_halogen_bonds( def get_hydrogen_or_halogen_bonds(
...@@ -79,89 +243,206 @@ def get_hydrogen_or_halogen_bonds( ...@@ -79,89 +243,206 @@ def get_hydrogen_or_halogen_bonds(
the log file ("labels"). the log file ("labels").
""" """
central_atoms = ["H"] if hydrogen_bond else ["I", "BR", "Br", "CL", "Cl", "F"] hbonds = {}
pdb_hbonds = Mol()
hbonds_labels = []
dist_cutoff = _set_default(dist_cutoff, HYDROGEN_HALOGEN_BOND_DIST_CUTOFF) dist_cutoff = _set_default(dist_cutoff, HYDROGEN_HALOGEN_BOND_DIST_CUTOFF)
angle_cutoff = _set_default(angle_cutoff, HYDROGEN_HALOGEN_BOND_ANGLE_CUTOFF) angle_cutoff = _set_default(angle_cutoff, HYDROGEN_HALOGEN_BOND_ANGLE_CUTOFF)
hbonds = {} # Check if hydrogen atoms added.
pdb_hbonds = Mol() lig_and_recep_have_hydrogens = ligand.has_hydrogens and receptor.has_hydrogens
hbonds_labels = []
# Calculate the distances. # Calculate the distances.
ligand_receptor_dists = _get_ligand_receptor_dists(ligand, receptor) ligand_receptor_dists = _get_ligand_receptor_dists(ligand, receptor)
# Now see if there's some sort of hydrogen (or halogen) bond between these # Get all donor-acceptor pairs that are near each other.
# two atoms. distance cutoff = 4, angle cutoff = 40. Note that these cutoffs close_donors_acceptors = [
# are generous. [ligand_atom, receptor_atom]
for ligand_atom, receptor_atom, dist in ligand_receptor_dists: for ligand_atom, receptor_atom, dist in ligand_receptor_dists
if ( if (
dist < dist_cutoff dist < dist_cutoff
and ligand_atom.element in ["O", "N"] and ligand_atom.element in ["O", "N"]
and receptor_atom.element in ["O", "N"] and receptor_atom.element in ["O", "N"]
): )
# now build a list of all the hydrogens (or halogens) close to these ]
# atoms
hydrogens = [] # Go through those pairs and find ones with complementary receptor/ligand
# labels.
for atm_index in ligand.all_atoms.keys(): for ligand_atom, receptor_atom in close_donors_acceptors:
element = ligand.all_atoms[atm_index].element hbond_detected = False
if element in central_atoms: lig_atm_hbond_infs = ligand.is_hbond_donor_acceptor(ligand_atom, hydrogen_bond)
dist = ligand.all_atoms[atm_index].coordinates.dist_to( recep_atm_hbond_infs = receptor.is_hbond_donor_acceptor(
ligand_atom.coordinates receptor_atom, hydrogen_bond
) )
if dist < max_donor_X_dist[element]: # 1.3 for H
# so it's a hydrogen (or halogen) # print(
ligand.all_atoms[atm_index].comment = "LIGAND" # "Note that above codewasworking before refactoring. Can always try restoring if you can't fix below."
hydrogens.append(ligand.all_atoms[atm_index]) # )
for atm_index in receptor.all_atoms.keys(): # print(
element = receptor.all_atoms[atm_index].element # "This whole process of winnowing down acceptable donor-acceptor pairs is falwed. You need to combinatorial check them all, especially since angles can be used to eliminate them in some cases."
if element in central_atoms: # )
dist = receptor.all_atoms[atm_index].coordinates.dist_to(
receptor_atom.coordinates for lig_atm_hbond_inf in lig_atm_hbond_infs:
) if hbond_detected:
if dist < max_donor_X_dist[element]: # 1.3 for H break
# so it's a hydrogen (or halogen)
receptor.all_atoms[atm_index].comment = "RECEPTOR" lig_donor_or_accept, lig_center_atom = lig_atm_hbond_inf
hydrogens.append(receptor.all_atoms[atm_index]) for recep_atm_hbond_inf in recep_atm_hbond_infs:
accept_donor_or_accept, accept_center_atom = recep_atm_hbond_inf
# now we need to check the angles
for hydrogen in hydrogens: if lig_donor_or_accept == accept_donor_or_accept:
if ( # Both acceptors or both donors. Doesn't work.
fabs( continue
180
- angle_between_three_points( center_atom = (
ligand_atom.coordinates, lig_center_atom
hydrogen.coordinates, if lig_donor_or_accept == "DONOR"
receptor_atom.coordinates, else accept_center_atom
) )
* 180.0
/ math.pi # Now that you've got the atoms, check the angles if appropriate.
) if lig_and_recep_have_hydrogens or not hydrogen_bond:
<= angle_cutoff # Hydrogens present and you're detecting hydrogen bonds, or
): # you're detecting halogen bonds (so hydrogens don't
hbonds_key = ( # matter).
"HDONOR_"
+ hydrogen.comment angle = angle_between_three_points(
+ "_" ligand_atom.coordinates,
+ receptor_atom.side_chain_or_backbone() center_atom.coordinates,
+ "_" receptor_atom.coordinates,
+ receptor_atom.structure
) )
pdb_hbonds.add_new_atom(ligand_atom.copy_of())
pdb_hbonds.add_new_atom(hydrogen.copy_of()) if fabs(180 - angle * 180.0 / math.pi) > angle_cutoff:
pdb_hbonds.add_new_atom(receptor_atom.copy_of()) # Angle is too big.
hashtable_entry_add_one(hbonds, hbonds_key) continue
hbonds_labels.append( # Now collect that data.
( comment = "RECEPTOR" if lig_donor_or_accept == "ACCEPTOR" else "LIGAND"
ligand_atom.string_id(),
hydrogen.string_id(), hbonds_key = (
receptor_atom.string_id(), "HDONOR_"
hydrogen.comment, + comment
) + "_"
+ receptor_atom.side_chain_or_backbone()
+ "_"
+ receptor_atom.structure
)
pdb_hbonds.add_new_atom(ligand_atom.copy_of())
pdb_hbonds.add_new_atom(center_atom.copy_of())
pdb_hbonds.add_new_atom(receptor_atom.copy_of())
hashtable_entry_add_one(hbonds, hbonds_key)
hbonds_labels.append(
(
ligand_atom.string_id(),
center_atom.string_id(),
receptor_atom.string_id(),
comment,
) )
)
# If you get here, it's identified a hydrogen bond. No need to
# keep checking for a hydrogen bond between these two atoms.
hbond_detected = True
break
# if len(lig_atm_hbond_infs) < len(recep_atm_hbond_infs):
# recep_atm_hbond_infs = _mimic_set_minus(
# recep_atm_hbond_infs, lig_atm_hbond_infs
# )
# # Note that transcrypt doesn't support set subtraction.
# # recep_atm_hbond_inf = recep_atm_hbond_inf - lig_atm_hbond_inf
# if len(recep_atm_hbond_infs) < len(lig_atm_hbond_infs):
# lig_atm_hbond_infs = _mimic_set_minus(
# lig_atm_hbond_infs, recep_atm_hbond_infs
# )
# # lig_atm_hbond_inf = lig_atm_hbond_inf - recep_atm_hbond_inf
# # If below occurs, there can't be donor-acceptor pair
# if len(lig_atm_hbond_infs) == 0:
# continue
# if len(recep_atm_hbond_infs) == 0:
# continue
# if len(lig_atm_hbond_infs) > 0 and len(recep_atm_hbond_infs) > 0:
# # They could be complementary (donor-acceptor), but they could also
# # both be donor or both be acceptor.
# # Pick first lig one (in case there are many)
# lig_category = lig_atm_hbond_infs[0]
# # Remove that one from recep options.
# recep_atm_hbond_infs = _mimic_set_minus(
# recep_atm_hbond_infs, [lig_category]
# )
# if len(recep_atm_hbond_infs) == 0:
# # There are no complementary donor-acceptor pairs.
# continue
# # Pick first of remaining receptor ones
# recep_category = recep_atm_hbond_infs[0]
# # print(
# # "Should therebe mutliple center atoms, and you need to pick between them?"
# # )
# center_atom = (
# recep_category[1] if lig_category[0] == "ACCEPTOR" else lig_category[1]
# )
# # Now that you've got the atoms, check the angles if appropriate.
# # TODO: Conditional
# angle = angle_between_three_points(
# ligand_atom.coordinates,
# center_atom.coordinates,
# receptor_atom.coordinates,
# )
# if fabs(180 - angle * 180.0 / math.pi) > angle_cutoff:
# # Angle is too big.
# continue
# # Now collect that data.
# comment = "RECEPTOR" if lig_category[0] == "ACCEPTOR" else "LIGAND"
# hbonds_key = (
# "HDONOR_"
# + comment
# + "_"
# + receptor_atom.side_chain_or_backbone()
# + "_"
# + receptor_atom.structure
# )
# pdb_hbonds.add_new_atom(ligand_atom.copy_of())
# pdb_hbonds.add_new_atom(center_atom.copy_of())
# pdb_hbonds.add_new_atom(receptor_atom.copy_of())
# hashtable_entry_add_one(hbonds, hbonds_key)
# # import pdb; pdb.set_trace()
# hbonds_labels.append(
# (
# ligand_atom.string_id(),
# center_atom.string_id(),
# receptor_atom.string_id(),
# comment,
# )
# )
# if not hydrogen_bond or (ligand.has_hydrogens and receptor.has_hydrogens):
# hbonds, pdb_hbonds, hbonds_labels = _detect_hbonds_by_angle(
# close_donors_acceptors, ligand, receptor, angle_cutoff, hydrogen_bond
# )
# else:
# hbonds, pdb_hbonds, hbonds_labels = _detect_hydrogen_bonds_without_hydrogens(
# close_donors_acceptors, ligand, receptor
# )
return { return {
"counts": hbonds, "counts": hbonds,
......
...@@ -64,33 +64,43 @@ def get_salt_bridges(ligand, receptor, cutoff=None): ...@@ -64,33 +64,43 @@ def get_salt_bridges(ligand, receptor, cutoff=None):
key = "SALT-BRIDGE_" + structure key = "SALT-BRIDGE_" + structure
for index in receptor_charge.indices: for index in receptor_charge.indices:
pdb_salt_bridges.add_new_atom(receptor.all_atoms[index].copy_of()) idx = int(index)
atom = receptor.all_atoms[idx]
pdb_salt_bridges.add_new_atom(atom.copy_of())
for index in ligand_charge.indices: for index in ligand_charge.indices:
pdb_salt_bridges.add_new_atom(ligand.all_atoms[index].copy_of()) idx = int(index)
atom = ligand.all_atoms[idx]
pdb_salt_bridges.add_new_atom(atom.copy_of())
hashtable_entry_add_one(salt_bridges, key) hashtable_entry_add_one(salt_bridges, key)
salt_bridge_labels.append( salt_bridge_labels.append(
( (
"[" (
+ " / ".join( (
[ "["
ligand.all_atoms[index].string_id() + " / ".join(
for index in ligand_charge.indices ligand.all_atoms[int(index)].string_id()
] for index in ligand_charge.indices
) )
+ "]", )
"[" + "]"
+ " / ".join( ),
[ (
receptor.all_atoms[index].string_id() (
for index in receptor_charge.indices "["
] + " / ".join(
) receptor.all_atoms[int(index)].string_id()
+ "]", for index in receptor_charge.indices
)
)
+ "]"
),
) )
) )
return { return {
"counts": salt_bridges, "counts": salt_bridges,
"mol": pdb_salt_bridges, "mol": pdb_salt_bridges,
......
...@@ -220,7 +220,7 @@ def _get_active_site_flexibility(active_site_flexibility, output): ...@@ -220,7 +220,7 @@ def _get_active_site_flexibility(active_site_flexibility, output):
return output return output
def _get_hbonds(hbonds, output, hydrogen_bond = True): def _get_hbonds(hbonds, output, hydrogen_bond=True):
name = "Hydrogen" if hydrogen_bond else "Halogen" name = "Hydrogen" if hydrogen_bond else "Halogen"
...@@ -525,7 +525,7 @@ def collect( ...@@ -525,7 +525,7 @@ def collect(
output = _get_active_site_flexibility(active_site_flexibility, output) output = _get_active_site_flexibility(active_site_flexibility, output)
output = _get_hbonds(hydrogen_bonds, output, True) # hydrogen bonds output = _get_hbonds(hydrogen_bonds, output, True) # hydrogen bonds
output = _get_hbonds(hydrogen_bonds, output, False) # halogen bonds output = _get_hbonds(halogen_bonds, output, False) # halogen bonds
output = _get_hydrophobics(hydrophobics, output) output = _get_hydrophobics(hydrophobics, output)
......