diff options
Diffstat (limited to 'plip/structure/detection.py')
| -rw-r--r-- | plip/structure/detection.py | 479 |
1 files changed, 479 insertions, 0 deletions
diff --git a/plip/structure/detection.py b/plip/structure/detection.py new file mode 100644 index 0000000..71fec63 --- /dev/null +++ b/plip/structure/detection.py @@ -0,0 +1,479 @@ +import itertools +from collections import defaultdict +from collections import namedtuple + +import numpy as np +from openbabel.openbabel import OBAtomAtomIter + +from plip.basic import config, logger +from plip.basic.supplemental import vecangle, vector, euclidean3d, projection +from plip.basic.supplemental import whichresnumber, whichrestype, whichchain + +logger = logger.get_logger() + +def filter_contacts(pairings): + """Filter interactions by two criteria: + 1. No interactions between the same residue (important for intra mode). + 2. No duplicate interactions (A with B and B with A, also important for intra mode).""" + if not config.INTRA: + return pairings + filtered1_pairings = [p for p in pairings if (p.resnr, p.reschain) != (p.resnr_l, p.reschain_l)] + already_considered = [] + filtered2_pairings = [] + for contact in filtered1_pairings: + try: + dist = 'D{}'.format(round(contact.distance, 2)) + except AttributeError: + try: + dist = 'D{}'.format(round(contact.distance_ah, 2)) + except AttributeError: + dist = 'D{}'.format(round(contact.distance_aw, 2)) + res1, res2 = ''.join([str(contact.resnr), contact.reschain]), ''.join( + [str(contact.resnr_l), contact.reschain_l]) + data = {res1, res2, dist} + if data not in already_considered: + filtered2_pairings.append(contact) + already_considered.append(data) + return filtered2_pairings + + +################################################## +# FUNCTIONS FOR DETECTION OF SPECIFIC INTERACTIONS +################################################## + +def hydrophobic_interactions(atom_set_a, atom_set_b): + """Detection of hydrophobic pliprofiler between atom_set_a (binding site) and atom_set_b (ligand). + Definition: All pairs of qualified carbon atoms within a distance of HYDROPH_DIST_MAX + """ + data = namedtuple('hydroph_interaction', 'bsatom bsatom_orig_idx ligatom ligatom_orig_idx ' + 'distance restype resnr reschain restype_l, resnr_l, reschain_l') + pairings = [] + for a, b in itertools.product(atom_set_a, atom_set_b): + if a.orig_idx == b.orig_idx: + continue + e = euclidean3d(a.atom.coords, b.atom.coords) + if not config.MIN_DIST < e < config.HYDROPH_DIST_MAX: + continue + restype, resnr, reschain = whichrestype(a.atom), whichresnumber(a.atom), whichchain(a.atom) + restype_l, resnr_l, reschain_l = whichrestype(b.orig_atom), whichresnumber(b.orig_atom), whichchain(b.orig_atom) + contact = data(bsatom=a.atom, bsatom_orig_idx=a.orig_idx, ligatom=b.atom, ligatom_orig_idx=b.orig_idx, + distance=e, restype=restype, resnr=resnr, + reschain=reschain, restype_l=restype_l, + resnr_l=resnr_l, reschain_l=reschain_l) + pairings.append(contact) + return filter_contacts(pairings) + + +def hbonds(acceptors, donor_pairs, protisdon, typ): + """Detection of hydrogen bonds between sets of acceptors and donor pairs. + Definition: All pairs of hydrogen bond acceptor and donors with + donor hydrogens and acceptor showing a distance within HBOND DIST MIN and HBOND DIST MAX + and donor angles above HBOND_DON_ANGLE_MIN + """ + data = namedtuple('hbond', 'a a_orig_idx d d_orig_idx h distance_ah distance_ad angle type protisdon resnr ' + 'restype reschain resnr_l restype_l reschain_l sidechain atype dtype') + pairings = [] + for acc, don in itertools.product(acceptors, donor_pairs): + if not typ == 'strong': + continue + # Regular (strong) hydrogen bonds + dist_ah = euclidean3d(acc.a.coords, don.h.coords) + dist_ad = euclidean3d(acc.a.coords, don.d.coords) + if not config.MIN_DIST < dist_ad < config.HBOND_DIST_MAX: + continue + vec1, vec2 = vector(don.h.coords, don.d.coords), vector(don.h.coords, acc.a.coords) + v = vecangle(vec1, vec2) + if not v > config.HBOND_DON_ANGLE_MIN: + continue + protatom = don.d.OBAtom if protisdon else acc.a.OBAtom + ligatom = don.d.OBAtom if not protisdon else acc.a.OBAtom + is_sidechain_hbond = protatom.GetResidue().GetAtomProperty(protatom, 8) # Check if sidechain atom + resnr = whichresnumber(don.d) if protisdon else whichresnumber(acc.a) + resnr_l = whichresnumber(acc.a_orig_atom) if protisdon else whichresnumber(don.d_orig_atom) + restype = whichrestype(don.d) if protisdon else whichrestype(acc.a) + restype_l = whichrestype(acc.a_orig_atom) if protisdon else whichrestype(don.d_orig_atom) + reschain = whichchain(don.d) if protisdon else whichchain(acc.a) + rechain_l = whichchain(acc.a_orig_atom) if protisdon else whichchain(don.d_orig_atom) + # Next line prevents H-Bonds within amino acids in intermolecular interactions + if config.INTRA is not None and whichresnumber(don.d) == whichresnumber(acc.a): + continue + # Next line prevents backbone-backbone H-Bonds + if config.INTRA is not None and protatom.GetResidue().GetAtomProperty(protatom, + 8) and ligatom.GetResidue().GetAtomProperty( + ligatom, 8): + continue + contact = data(a=acc.a, a_orig_idx=acc.a_orig_idx, d=don.d, d_orig_idx=don.d_orig_idx, h=don.h, + distance_ah=dist_ah, distance_ad=dist_ad, angle=v, type=typ, protisdon=protisdon, + resnr=resnr, restype=restype, reschain=reschain, resnr_l=resnr_l, + restype_l=restype_l, reschain_l=rechain_l, sidechain=is_sidechain_hbond, + atype=acc.a.type, dtype=don.d.type) + pairings.append(contact) + return filter_contacts(pairings) + + +def pistacking(rings_bs, rings_lig): + """Return all pi-stackings between the given aromatic ring systems in receptor and ligand.""" + data = namedtuple( + 'pistack', + 'proteinring ligandring distance angle offset type restype resnr reschain restype_l resnr_l reschain_l') + pairings = [] + for r, l in itertools.product(rings_bs, rings_lig): + # DISTANCE AND RING ANGLE CALCULATION + d = euclidean3d(r.center, l.center) + b = vecangle(r.normal, l.normal) + a = min(b, 180 - b if not 180 - b < 0 else b) # Smallest of two angles, depending on direction of normal + + # RING CENTER OFFSET CALCULATION (project each ring center into the other ring) + proj1 = projection(l.normal, l.center, r.center) + proj2 = projection(r.normal, r.center, l.center) + offset = min(euclidean3d(proj1, l.center), euclidean3d(proj2, r.center)) + + # RECEPTOR DATA + resnr, restype, reschain = whichresnumber(r.atoms[0]), whichrestype(r.atoms[0]), whichchain(r.atoms[0]) + resnr_l, restype_l, reschain_l = whichresnumber(l.orig_atoms[0]), whichrestype( + l.orig_atoms[0]), whichchain(l.orig_atoms[0]) + + # SELECTION BY DISTANCE, ANGLE AND OFFSET + passed = False + if not config.MIN_DIST < d < config.PISTACK_DIST_MAX: + continue + if 0 < a < config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX: + ptype = 'P' + passed = True + if 90 - config.PISTACK_ANG_DEV < a < 90 + config.PISTACK_ANG_DEV and offset < config.PISTACK_OFFSET_MAX: + ptype = 'T' + passed = True + if passed: + contact = data(proteinring=r, ligandring=l, distance=d, angle=a, offset=offset, + type=ptype, resnr=resnr, restype=restype, reschain=reschain, + resnr_l=resnr_l, restype_l=restype_l, reschain_l=reschain_l) + pairings.append(contact) + return filter_contacts(pairings) + + +def pication(rings, pos_charged, protcharged): + """Return all pi-Cation interaction between aromatic rings and positively charged groups. + For tertiary and quaternary amines, check also the angle between the ring and the nitrogen. + """ + data = namedtuple( + 'pication', 'ring charge distance offset type restype resnr reschain restype_l resnr_l reschain_l protcharged') + pairings = [] + if len(rings) == 0 or len(pos_charged) == 0: + return pairings + for ring in rings: + c = ring.center + for p in pos_charged: + d = euclidean3d(c, p.center) + # Project the center of charge into the ring and measure distance to ring center + proj = projection(ring.normal, ring.center, p.center) + offset = euclidean3d(proj, ring.center) + if not config.MIN_DIST < d < config.PICATION_DIST_MAX or not offset < config.PISTACK_OFFSET_MAX: + continue + if type(p).__name__ == 'lcharge' and p.fgroup == 'tertamine': + # Special case here if the ligand has a tertiary amine, check an additional angle + # Otherwise, we might have have a pi-cation interaction 'through' the ligand + n_atoms = [a_neighbor for a_neighbor in OBAtomAtomIter(p.atoms[0].OBAtom)] + n_atoms_coords = [(a.x(), a.y(), a.z()) for a in n_atoms] + amine_normal = np.cross(vector(n_atoms_coords[0], n_atoms_coords[1]), + vector(n_atoms_coords[2], n_atoms_coords[0])) + b = vecangle(ring.normal, amine_normal) + # Smallest of two angles, depending on direction of normal + a = min(b, 180 - b if not 180 - b < 0 else b) + if not a > 30.0: + resnr, restype = whichresnumber(ring.atoms[0]), whichrestype(ring.atoms[0]) + reschain = whichchain(ring.atoms[0]) + resnr_l, restype_l = whichresnumber(p.orig_atoms[0]), whichrestype(p.orig_atoms[0]) + reschain_l = whichchain(p.orig_atoms[0]) + contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular', + restype=restype, resnr=resnr, reschain=reschain, + restype_l=restype_l, resnr_l=resnr_l, reschain_l=reschain_l, + protcharged=protcharged) + pairings.append(contact) + break + resnr = whichresnumber(p.atoms[0]) if protcharged else whichresnumber(ring.atoms[0]) + resnr_l = whichresnumber(ring.orig_atoms[0]) if protcharged else whichresnumber(p.orig_atoms[0]) + restype = whichrestype(p.atoms[0]) if protcharged else whichrestype(ring.atoms[0]) + restype_l = whichrestype(ring.orig_atoms[0]) if protcharged else whichrestype(p.orig_atoms[0]) + reschain = whichchain(p.atoms[0]) if protcharged else whichchain(ring.atoms[0]) + reschain_l = whichchain(ring.orig_atoms[0]) if protcharged else whichchain(p.orig_atoms[0]) + contact = data(ring=ring, charge=p, distance=d, offset=offset, type='regular', restype=restype, + resnr=resnr, reschain=reschain, restype_l=restype_l, resnr_l=resnr_l, + reschain_l=reschain_l, protcharged=protcharged) + pairings.append(contact) + return filter_contacts(pairings) + + +def saltbridge(poscenter, negcenter, protispos): + """Detect all salt bridges (pliprofiler between centers of positive and negative charge)""" + data = namedtuple( + 'saltbridge', 'positive negative distance protispos resnr restype reschain resnr_l restype_l reschain_l') + pairings = [] + for pc, nc in itertools.product(poscenter, negcenter): + if not config.MIN_DIST < euclidean3d(pc.center, nc.center) < config.SALTBRIDGE_DIST_MAX: + continue + resnr = pc.resnr if protispos else nc.resnr + resnr_l = whichresnumber(nc.orig_atoms[0]) if protispos else whichresnumber(pc.orig_atoms[0]) + restype = pc.restype if protispos else nc.restype + restype_l = whichrestype(nc.orig_atoms[0]) if protispos else whichrestype(pc.orig_atoms[0]) + reschain = pc.reschain if protispos else nc.reschain + reschain_l = whichchain(nc.orig_atoms[0]) if protispos else whichchain(pc.orig_atoms[0]) + contact = data(positive=pc, negative=nc, distance=euclidean3d(pc.center, nc.center), protispos=protispos, + resnr=resnr, restype=restype, reschain=reschain, resnr_l=resnr_l, restype_l=restype_l, + reschain_l=reschain_l) + pairings.append(contact) + return filter_contacts(pairings) + + +def halogen(acceptor, donor): + """Detect all halogen bonds of the type Y-O...X-C""" + data = namedtuple('halogenbond', 'acc acc_orig_idx don don_orig_idx distance don_angle acc_angle restype ' + 'resnr reschain restype_l resnr_l reschain_l donortype acctype sidechain') + pairings = [] + for acc, don in itertools.product(acceptor, donor): + dist = euclidean3d(acc.o.coords, don.x.coords) + if not config.MIN_DIST < dist < config.HALOGEN_DIST_MAX: + continue + vec1, vec2 = vector(acc.o.coords, acc.y.coords), vector(acc.o.coords, don.x.coords) + vec3, vec4 = vector(don.x.coords, acc.o.coords), vector(don.x.coords, don.c.coords) + acc_angle, don_angle = vecangle(vec1, vec2), vecangle(vec3, vec4) + is_sidechain_hal = acc.o.OBAtom.GetResidue().GetAtomProperty(acc.o.OBAtom, 8) # Check if sidechain atom + if not config.HALOGEN_ACC_ANGLE - config.HALOGEN_ANGLE_DEV < acc_angle \ + < config.HALOGEN_ACC_ANGLE + config.HALOGEN_ANGLE_DEV: + continue + if not config.HALOGEN_DON_ANGLE - config.HALOGEN_ANGLE_DEV < don_angle \ + < config.HALOGEN_DON_ANGLE + config.HALOGEN_ANGLE_DEV: + continue + restype, reschain, resnr = whichrestype(acc.o), whichchain(acc.o), whichresnumber(acc.o) + restype_l, reschain_l, resnr_l = whichrestype(don.orig_x), whichchain(don.orig_x), whichresnumber(don.orig_x) + contact = data(acc=acc, acc_orig_idx=acc.o_orig_idx, don=don, don_orig_idx=don.x_orig_idx, + distance=dist, don_angle=don_angle, acc_angle=acc_angle, + restype=restype, resnr=resnr, + reschain=reschain, restype_l=restype_l, + reschain_l=reschain_l, resnr_l=resnr_l, donortype=don.x.OBAtom.GetType(), acctype=acc.o.type, + sidechain=is_sidechain_hal) + pairings.append(contact) + return filter_contacts(pairings) + + +def water_bridges(bs_hba, lig_hba, bs_hbd, lig_hbd, water): + """Find water-bridged hydrogen bonds between ligand and protein. For now only considers bridged of first degree.""" + data = namedtuple('waterbridge', 'a a_orig_idx atype d d_orig_idx dtype h water water_orig_idx distance_aw ' + 'distance_dw d_angle w_angle type resnr restype reschain resnr_l restype_l reschain_l protisdon') + pairings = [] + # First find all acceptor-water pairs with distance within d + # and all donor-water pairs with distance within d and angle greater theta + lig_aw, prot_aw, lig_dw, prot_hw = [], [], [], [] + for w in water: + for acc1 in lig_hba: + dist = euclidean3d(acc1.a.coords, w.oxy.coords) + if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST: + lig_aw.append((acc1, w, dist)) + for acc2 in bs_hba: + dist = euclidean3d(acc2.a.coords, w.oxy.coords) + if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST: + prot_aw.append((acc2, w, dist)) + for don1 in lig_hbd: + dist = euclidean3d(don1.d.coords, w.oxy.coords) + d_angle = vecangle(vector(don1.h.coords, don1.d.coords), vector(don1.h.coords, w.oxy.coords)) + if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST \ + and d_angle > config.WATER_BRIDGE_THETA_MIN: + lig_dw.append((don1, w, dist, d_angle)) + for don2 in bs_hbd: + dist = euclidean3d(don2.d.coords, w.oxy.coords) + d_angle = vecangle(vector(don2.h.coords, don2.d.coords), vector(don2.h.coords, w.oxy.coords)) + if config.WATER_BRIDGE_MINDIST <= dist <= config.WATER_BRIDGE_MAXDIST \ + and d_angle > config.WATER_BRIDGE_THETA_MIN: + prot_hw.append((don2, w, dist, d_angle)) + + for l, p in itertools.product(lig_aw, prot_hw): + acc, wl, distance_aw = l + don, wd, distance_dw, d_angle = p + 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)) + 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) + resnr_l, reschain_l, restype_l = whichresnumber(acc.a_orig_atom), whichchain( + acc.a_orig_atom), whichrestype(acc.a_orig_atom) + contact = data(a=acc.a, a_orig_idx=acc.a_orig_idx, atype=acc.a.type, d=don.d, d_orig_idx=don.d_orig_idx, + dtype=don.d.type, h=don.h, water=wl.oxy, water_orig_idx=wl.oxy_orig_idx, + distance_aw=distance_aw, distance_dw=distance_dw, d_angle=d_angle, w_angle=w_angle, + type='first_deg', resnr=resnr, restype=restype, + reschain=reschain, restype_l=restype_l, resnr_l=resnr_l, reschain_l=reschain_l, protisdon=True) + pairings.append(contact) + for p, l in itertools.product(prot_aw, lig_dw): + acc, wl, distance_aw = p + don, wd, distance_dw, d_angle = l + 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)) + 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) + resnr_l, reschain_l, restype_l = whichresnumber(don.d_orig_atom), whichchain( + don.d_orig_atom), whichrestype(don.d_orig_atom) + contact = data(a=acc.a, a_orig_idx=acc.a_orig_idx, atype=acc.a.type, d=don.d, d_orig_idx=don.d_orig_idx, + dtype=don.d.type, h=don.h, water=wl.oxy, water_orig_idx=wl.oxy_orig_idx, + distance_aw=distance_aw, distance_dw=distance_dw, + d_angle=d_angle, w_angle=w_angle, type='first_deg', resnr=resnr, + restype=restype, reschain=reschain, + restype_l=restype_l, reschain_l=reschain_l, resnr_l=resnr_l, protisdon=False) + pairings.append(contact) + return filter_contacts(pairings) + + +def metal_complexation(metals, metal_binding_lig, metal_binding_bs): + """Find all metal complexes between metals and appropriate groups in both protein and ligand, as well as water""" + data = namedtuple('metal_complex', 'metal metal_orig_idx metal_type target target_orig_idx target_type ' + 'coordination_num distance resnr restype ' + 'reschain restype_l reschain_l resnr_l location rms, geometry num_partners complexnum') + pairings_dict = {} + pairings = [] + # #@todo Refactor + metal_to_id = {} + metal_to_orig_atom = {} + for metal, target in itertools.product(metals, metal_binding_lig + metal_binding_bs): + distance = euclidean3d(metal.m.coords, target.atom.coords) + if not distance < config.METAL_DIST_MAX: + continue + if metal.m not in pairings_dict: + pairings_dict[metal.m] = [(target, distance), ] + metal_to_id[metal.m] = metal.m_orig_idx + metal_to_orig_atom[metal.m] = metal.orig_m + else: + pairings_dict[metal.m].append((target, distance)) + for cnum, metal in enumerate(pairings_dict): + rms = 0.0 + excluded = [] + # cnum +1 being the complex number + contact_pairs = pairings_dict[metal] + num_targets = len(contact_pairs) + vectors_dict = defaultdict(list) + for contact_pair in contact_pairs: + target, distance = contact_pair + vectors_dict[target.atom.idx].append(vector(metal.coords, target.atom.coords)) + + # Listing of coordination numbers and their geometries + configs = {2: ['linear', ], + 3: ['trigonal.planar', 'trigonal.pyramidal'], + 4: ['tetrahedral', 'square.planar'], + 5: ['trigonal.bipyramidal', 'square.pyramidal'], + 6: ['octahedral', ]} + + # Angle signatures for each geometry (as seen from each target atom) + ideal_angles = {'linear': [[180.0]] * 2, + 'trigonal.planar': [[120.0, 120.0]] * 3, + 'trigonal.pyramidal': [[109.5, 109.5]] * 3, + 'tetrahedral': [[109.5, 109.5, 109.5, 109.5]] * 4, + 'square.planar': [[90.0, 90.0, 90.0, 90.0]] * 4, + 'trigonal.bipyramidal': [[120.0, 120.0, 90.0, 90.0]] * 3 + [[90.0, 90.0, 90.0, 180.0]] * 2, + 'square.pyramidal': [[90.0, 90.0, 90.0, 180.0]] * 4 + [[90.0, 90.0, 90.0, 90.0]], + 'octahedral': [[90.0, 90.0, 90.0, 90.0, 180.0]] * 6} + angles_dict = {} + + for target in vectors_dict: + cur_vector = vectors_dict[target] + other_vectors = [] + for t in vectors_dict: + if not t == target: + [other_vectors.append(x) for x in vectors_dict[t]] + angles = [vecangle(pair[0], pair[1]) for pair in itertools.product(cur_vector, other_vectors)] + angles_dict[target] = angles + + all_total = [] # Record fit information for each geometry tested + gdata = namedtuple('gdata', 'geometry rms coordination excluded diff_targets') # Geometry Data + # Can't specify geometry with only one target + if num_targets == 1: + final_geom = 'NA' + final_coo = 1 + excluded = [] + rms = 0.0 + else: + for coo in sorted(configs, reverse=True): # Start with highest coordination number + geometries = configs[coo] + for geometry in geometries: + signature = ideal_angles[geometry] # Set of ideal angles for geometry, from each perspective + geometry_total = 0 + geometry_scores = [] # All scores for one geometry (from all subsignatures) + used_up_targets = [] # Use each target just once for a subsignature + not_used = [] + coo_diff = num_targets - coo # How many more observed targets are there? + + # Find best match for each subsignature + for subsignature in signature: # Ideal angles from one perspective + best_target = None # There's one best-matching target for each subsignature + best_target_score = 999 + + for k, target in enumerate(angles_dict): + if target not in used_up_targets: + observed_angles = angles_dict[target] # Observed angles from perspective of one target + single_target_scores = [] + used_up_observed_angles = [] + for i, ideal_angle in enumerate(subsignature): + # For each angle in the signature, find the best-matching observed angle + best_match = None + best_match_diff = 999 + for j, observed_angle in enumerate(observed_angles): + if j not in used_up_observed_angles: + diff = abs(ideal_angle - observed_angle) + if diff < best_match_diff: + best_match_diff = diff + best_match = j + if best_match is not None: + used_up_observed_angles.append(best_match) + single_target_scores.append(best_match_diff) + # Calculate RMS for target angles + target_total = sum([x ** 2 for x in single_target_scores]) ** 0.5 # Tot. score targ/sig + if target_total < best_target_score: + best_target_score = target_total + best_target = target + + used_up_targets.append(best_target) + geometry_scores.append(best_target_score) + # Total score is mean of RMS values + geometry_total = np.mean(geometry_scores) + # Record the targets not used for excluding them when deciding for a final geometry + [not_used.append(target) for target in angles_dict if target not in used_up_targets] + all_total.append(gdata(geometry=geometry, rms=geometry_total, coordination=coo, + excluded=not_used, diff_targets=coo_diff)) + + # Make a decision here. Starting with the geometry with lowest difference in ideal and observed partners ... + # Check if the difference between the RMS to the next best solution is not larger than 0.5 + if not num_targets == 1: # Can't decide for any geoemtry in that case + all_total = sorted(all_total, key=lambda x: abs(x.diff_targets)) + for i, total in enumerate(all_total): + next_total = all_total[i + 1] + this_rms, next_rms = total.rms, next_total.rms + diff_to_next = next_rms - this_rms + if diff_to_next > 0.5: + final_geom, final_coo, rms, excluded = total.geometry, total.coordination, total.rms, total.excluded + break + elif next_total.rms < 3.5: + final_geom, final_coo, = next_total.geometry, next_total.coordination + rms, excluded = next_total.rms, next_total.excluded + break + elif i == len(all_total) - 2: + final_geom, final_coo, rms, excluded = "NA", "NA", float('nan'), [] + break + + # Record all contact pairing, excluding those with targets superfluous for chosen geometry + only_water = set([x[0].location for x in contact_pairs]) == {'water'} + if not only_water: # No complex if just with water as targets + logger.info(f'metal ion {metal.type} complexed with {final_geom} geometry (coo. number {final_coo}/ {num_targets} observed)') + for contact_pair in contact_pairs: + target, distance = contact_pair + if target.atom.idx not in excluded: + metal_orig_atom = metal_to_orig_atom[metal] + restype_l, reschain_l, resnr_l = whichrestype(metal_orig_atom), whichchain( + metal_orig_atom), whichresnumber(metal_orig_atom) + contact = data(metal=metal, metal_orig_idx=metal_to_id[metal], metal_type=metal.type, + target=target, target_orig_idx=target.atom_orig_idx, target_type=target.type, + coordination_num=final_coo, distance=distance, resnr=target.resnr, + restype=target.restype, reschain=target.reschain, location=target.location, + rms=rms, geometry=final_geom, num_partners=num_targets, complexnum=cnum + 1, + resnr_l=resnr_l, restype_l=restype_l, reschain_l=reschain_l) + pairings.append(contact) + return filter_contacts(pairings) |