From ebd0451db8dedc8003d94127d8753bd314413daa Mon Sep 17 00:00:00 2001 From: clami66 Date: Thu, 16 May 2024 11:23:51 +0200 Subject: [PATCH] variable bond length based on element types --- src/DockQ/DockQ.py | 23 ++++++-- src/DockQ/constants.py | 124 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 141 insertions(+), 6 deletions(-) create mode 100644 src/DockQ/constants.py diff --git a/src/DockQ/DockQ.py b/src/DockQ/DockQ.py index 5bfa62b..c3cd0ac 100755 --- a/src/DockQ/DockQ.py +++ b/src/DockQ/DockQ.py @@ -18,6 +18,7 @@ try: from .operations import residue_distances, get_fnat_stats from .parsers import PDBParser, MMCIFParser + from .constants import * except ImportError: warnings.warn( """WARNING: It looks like cython is not working, @@ -25,6 +26,7 @@ ) from operations_nocy import residue_distances, get_fnat_stats from parsers import PDBParser, MMCIFParser + from constants import * def parse_args(): @@ -197,8 +199,11 @@ def calc_sym_corrected_lrmsd( get_atoms_per_residue((aligned_ref_receptor, aligned_sample_receptor), what="receptor") ) - sample_ligand_atoms_ids = [atom.id for atom in sample_ligand.get_atoms() if atom.id != "H"] - ref_ligand_atoms_ids = [atom.id for atom in ref_ligand.get_atoms() if atom.id != "H"] + sample_ligand_atoms_ids = [atom.id for atom in sample_ligand.get_atoms()] + sample_ligand_atoms_ele = [atom.element for atom in sample_ligand.get_atoms()] + + ref_ligand_atoms_ids = [atom.id for atom in ref_ligand.get_atoms()] + ref_ligand_atoms_ele = [atom.element for atom in ref_ligand.get_atoms()] sample_ligand_atoms = np.array([atom.coord for atom in sample_ligand.get_atoms() if atom.id in ref_ligand_atoms_ids]) ref_ligand_atoms = np.array([atom.coord for atom in ref_ligand.get_atoms() if atom.id in sample_ligand_atoms_ids]) @@ -211,13 +216,13 @@ def calc_sym_corrected_lrmsd( sample_rotated_ligand_atoms = np.dot(sample_ligand_atoms, rot) + tran - sample_graph = create_graph(sample_ligand_atoms) - ref_graph = create_graph(ref_ligand_atoms) + sample_graph = create_graph(sample_ligand_atoms, sample_ligand_atoms_ele) + ref_graph = create_graph(ref_ligand_atoms, ref_ligand_atoms_ele) min_Lrms = float("inf") best_mapping = None + print(sample_graph, ref_graph) for isomorphism in nx.vf2pp_all_isomorphisms(sample_graph, ref_graph): - model_i = list(isomorphism.keys()) native_i = list(isomorphism.values()) @@ -593,14 +598,20 @@ def run_on_chains( return info -def create_graph(atom_list, threshold=2.0): +def create_graph(atom_list, atom_ids): import networkx as nx G = nx.Graph() + for i, atom_i in enumerate(atom_list): + cr_i = covalent_radius[atom_ids[i]] for j, atom_j in enumerate(atom_list): + cr_j = covalent_radius[atom_ids[j]] distance = np.linalg.norm(atom_i - atom_j) + threshold = (cr_i + cr_j + bond_tolerance) if i != j else 0.0 + #print(atom_ids[i], atom_ids[j], threshold) if distance < threshold: # Adjust threshold as needed G.add_edge(i, j) + return G diff --git a/src/DockQ/constants.py b/src/DockQ/constants.py new file mode 100644 index 0000000..6fcd5f4 --- /dev/null +++ b/src/DockQ/constants.py @@ -0,0 +1,124 @@ +from typing import Dict + +bond_tolerance: float = 0.4 + +covalent_radius: Dict = { + "H": 0.31, + "HE": 0.28, + "LI": 1.28, + "BE": 0.96, + "B": 0.84, + "C": 0.76, + "N": 0.71, + "O": 0.66, + "F": 0.57, + "NE": 0.58, + "NA": 1.66, + "MG": 1.41, + "AL": 1.21, + "SI": 1.11, + "P": 1.07, + "S": 1.05, + "CL": 1.02, + "AR": 1.06, + "K": 2.03, + "CA": 1.76, + "SC": 1.7, + "TI": 1.6, + "V": 1.53, + "CR": 1.39, + "MN": 1.5, + "FE": 1.42, + "CO": 1.38, + "NI": 1.24, + "CU": 1.32, + "ZN": 1.22, + "GA": 1.22, + "GE": 1.2, + "AS": 1.19, + "SE": 1.2, + "BR": 1.2, + "KR": 1.16, + "RB": 2.2, + "SR": 1.95, + "Y": 1.9, + "ZR": 1.75, + "NB": 1.64, + "MO": 1.54, + "TC": 1.47, + "RU": 1.46, + "RH": 1.42, + "PD": 1.39, + "AG": 1.45, + "CD": 1.44, + "IN": 1.42, + "SN": 1.39, + "SB": 1.39, + "TE": 1.38, + "I": 1.39, + "XE": 1.4, + "CS": 2.44, + "BA": 2.15, + "LA": 2.07, + "CE": 2.04, + "PR": 2.03, + "ND": 2.01, + "PM": 1.99, + "SM": 1.98, + "EU": 1.98, + "GD": 1.96, + "TB": 1.94, + "DY": 1.92, + "HO": 1.92, + "ER": 1.89, + "TM": 1.9, + "YB": 1.87, + "LU": 1.87, + "HF": 1.75, + "TA": 1.7, + "W": 1.62, + "RE": 1.51, + "OS": 1.44, + "IR": 1.41, + "PT": 1.36, + "AU": 1.36, + "HG": 1.32, + "TL": 1.45, + "PB": 1.46, + "BI": 1.48, + "PO": 1.4, + "AT": 1.5, + "RN": 1.5, + "FR": 2.6, + "RA": 2.21, + "AC": 2.15, + "TH": 2.06, + "PA": 2.0, + "U": 1.96, + "NP": 1.9, + "PU": 1.87, + "AM": 1.8, + "CM": 1.69, + "BK": 2.0, + "CF": 2.0, + "ES": 2.0, + "FM": 2.0, + "MD": 2.0, + "NO": 2.0, + "LR": 2.0, + "RF": 2.0, + "DB": 2.0, + "SG": 2.0, + "BH": 2.0, + "HS": 2.0, + "MT": 2.0, + "DS": 2.0, + "RG": 2.0, + "CN": 2.0, + "UUT": 2.0, + "UUQ": 2.0, + "UUP": 2.0, + "UUH": 2.0, + "UUS": 2.0, + "UUO": 2.0, +}