From 6bf2442210c90f2da64e7d957fc0cdf9da4bc594 Mon Sep 17 00:00:00 2001 From: Katja Linnemann Date: Thu, 7 Mar 2024 18:25:49 +0100 Subject: [PATCH 1/3] Fixing omega angle calculation in water bridge detection This closes #122 --- plip/structure/detection.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/plip/structure/detection.py b/plip/structure/detection.py index 71fec63..7af526b 100644 --- a/plip/structure/detection.py +++ b/plip/structure/detection.py @@ -291,7 +291,7 @@ def water_bridges(bs_hba, lig_hba, bs_hbd, lig_hbd, water): if not wl.oxy == wd.oxy: continue # Same water molecule and angle within omega - w_angle = vecangle(vector(acc.a.coords, wl.oxy.coords), vector(wl.oxy.coords, don.h.coords)) + w_angle = vecangle(vector(wl.oxy.coords, acc.a.coords), vector(wl.oxy.coords, don.h.coords)) if not config.WATER_BRIDGE_OMEGA_MIN < w_angle < config.WATER_BRIDGE_OMEGA_MAX: continue resnr, reschain, restype = whichresnumber(don.d), whichchain(don.d), whichrestype(don.d) @@ -309,7 +309,7 @@ def water_bridges(bs_hba, lig_hba, bs_hbd, lig_hbd, water): if not wl.oxy == wd.oxy: continue # Same water molecule and angle within omega - w_angle = vecangle(vector(acc.a.coords, wl.oxy.coords), vector(wl.oxy.coords, don.h.coords)) + w_angle = vecangle(vector(wl.oxy.coords, acc.a.coords), vector(wl.oxy.coords, don.h.coords)) if not config.WATER_BRIDGE_OMEGA_MIN < w_angle < config.WATER_BRIDGE_OMEGA_MAX: continue resnr, reschain, restype = whichresnumber(acc.a), whichchain(acc.a), whichrestype(acc.a) From 4bd8648ab988eae8813c0ae9b2c3163badf964b4 Mon Sep 17 00:00:00 2001 From: Katja Linnemann Date: Thu, 7 Mar 2024 18:31:35 +0100 Subject: [PATCH 2/3] Fixing error in water bridge filtering Instead of the bond with angle omega closest to the ideal value the one further away was kept. --- plip/structure/preparation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/plip/structure/preparation.py b/plip/structure/preparation.py index 32326cf..0e77ca2 100644 --- a/plip/structure/preparation.py +++ b/plip/structure/preparation.py @@ -878,7 +878,7 @@ def refine_water_bridges(wbridges, hbonds_ldon, hbonds_pdon): if (wbridge.water.idx, wbridge.a.idx) not in wb_dict: wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge else: - if abs(omega - wb_dict[(wbridge.water.idx, wbridge.a.idx)].w_angle) < abs(omega - wbridge.w_angle): + if abs(omega - wb_dict[(wbridge.water.idx, wbridge.a.idx)].w_angle) > abs(omega - wbridge.w_angle): wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge for wb_tuple in wb_dict: water, acceptor = wb_tuple @@ -888,7 +888,7 @@ def refine_water_bridges(wbridges, hbonds_ldon, hbonds_pdon): wb_dict2[water].append((abs(omega - wb_dict[wb_tuple].w_angle), wb_dict[wb_tuple])) wb_dict2[water] = sorted(wb_dict2[water], key=lambda x: x[0]) else: - if wb_dict2[water][1][0] < abs(omega - wb_dict[wb_tuple].w_angle): + if wb_dict2[water][1][0] > abs(omega - wb_dict[wb_tuple].w_angle): wb_dict2[water] = [wb_dict2[water][0], (wb_dict[wb_tuple].w_angle, wb_dict[wb_tuple])] filtered_wb = [] From 74d7adc98d0aac6cd3383dc236e3fae91a50a613 Mon Sep 17 00:00:00 2001 From: Katja Linnemann Date: Tue, 26 Mar 2024 14:32:13 +0100 Subject: [PATCH 3/3] Fixing another error in water bridge filtering When checking, whether the donor atom already forms a hydrogen bond, the openBabel molecules atom index instead of the original index in the pdb file was used, resulting in random hydrogen bonds in other parts of the structure preventing water bridges from showing up and potentially not filtering out the ones that actually needed filtering. --- plip/structure/preparation.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/plip/structure/preparation.py b/plip/structure/preparation.py index 0e77ca2..2da6fed 100644 --- a/plip/structure/preparation.py +++ b/plip/structure/preparation.py @@ -868,18 +868,24 @@ def refine_pi_cation_laro(all_picat, stacks): def refine_water_bridges(wbridges, hbonds_ldon, hbonds_pdon): """A donor atom already forming a hydrogen bond is not allowed to form a water bridge. Each water molecule can only be donor for two water bridges, selecting the constellation with the omega angle closest to 110 deg.""" - donor_atoms_hbonds = [hb.d.idx for hb in hbonds_ldon + hbonds_pdon] + donor_atoms_hbonds = [hb.d_orig_idx for hb in hbonds_ldon + hbonds_pdon] wb_dict = {} wb_dict2 = {} omega = 110.0 # Just one hydrogen bond per donor atom - for wbridge in [wb for wb in wbridges if wb.d.idx not in donor_atoms_hbonds]: - if (wbridge.water.idx, wbridge.a.idx) not in wb_dict: - wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge + # (TODO: looks wrong, what about donor atoms with multiple H? Does PLIP cover that somehow? + # and only one donor per water-acceptor combination + # (each water-acceptor combo only kept once after this loop, + # the one with angle at water closest to omega) + for wbridge in [wb for wb in wbridges if wb.d_orig_idx not in donor_atoms_hbonds]: + if (wbridge.water_orig_idx, wbridge.a_orig_idx) not in wb_dict: + wb_dict[(wbridge.water_orig_idx, wbridge.a_orig_idx)] = wbridge else: - if abs(omega - wb_dict[(wbridge.water.idx, wbridge.a.idx)].w_angle) > abs(omega - wbridge.w_angle): - wb_dict[(wbridge.water.idx, wbridge.a.idx)] = wbridge + if abs(omega - wb_dict[(wbridge.water_orig_idx, wbridge.a_orig_idx)].w_angle) > abs(omega - wbridge.w_angle): + wb_dict[(wbridge.water_orig_idx, wbridge.a_orig_idx)] = wbridge + # Just two hydrogen bonds per water molecule + # only the two water-acceptor combos with angle closest to omega are kept for wb_tuple in wb_dict: water, acceptor = wb_tuple if water not in wb_dict2: