diff --git a/kinoml/core/kinase.py b/kinoml/core/kinase.py new file mode 100644 index 00000000..219e8479 --- /dev/null +++ b/kinoml/core/kinase.py @@ -0,0 +1,74 @@ +""" +Defines the Kinase class + +[WIP] +""" + + +class Kinase(object): + def __init__( + self, + pdb, + chain, + kinase_id, + name, + struct_id, + ligand, + pocket_seq, + numbering, + key_res, + dihedrals, + distances, + mean_dist, + ): + """This script defines a Kinase class of which any kinase can be represented as an object with the + following parameters: + + Parameters + ---------- + pdb: str + The PDB code of the structure. + chain: str + The chain index of the structure. + kinase_id: int + The standard ID of a kinase enforced by the KLIFS database. + name: str + The standard name of the kinase used by the KLIFS database. + struct_id: int + The ID associated with a specific chain in the pdb structure of a kinase. + ligand: str + The ligand name as it appears in the pdb file. + pocket_seq: str + The 85 discontinuous residues (from multi-sequence alignment) that define the binding pocket of a kinase. + numbering: list of int + The residue indices of the 85 pocket residues specific to the structure. + key_res: list of int + A list of residue indices that are relevant to the collective variables. + dihedrals: list of floats + A list (one frame) or lists (multiple frames) of dihedrals relevant to kinase conformation. + distances: list of floats + A list (one frame) or lists (multiple frames) of intramolecular distances relevant to kinase conformation. + mean_dist: float + A float (one frame) or a list of floats (multiple frames), which is the mean pairwise distance between + ligand heavy atoms and the CAs of the 85 pocket residues. + + .. todo :: + + This is WAY too many positional arguments. Can we use kwargs instead, or somehow simplify the positional arguments into logical groups? + Many of these will be optional if we want to represent aspects of a structure, so there's no need to make them all requiredself. + Also, we will likely not want to mix features (distances, dihedrals) with structural information directly. + + """ + + self.pdb = pdb + self.chain = chain + self.kinase_id = kinase_id + self.name = name + self.struct_id = struct_id + self.ligand = ligand + self.pocket_seq = pocket_seq + self.numbering = numbering + self.key_res = key_res + self.dihedrals = dihedrals + self.distances = distances + self.mean_dist = mean_dist diff --git a/kinoml/features/dunbrack_cluster.py b/kinoml/features/dunbrack_cluster.py new file mode 100644 index 00000000..68d4a117 --- /dev/null +++ b/kinoml/features/dunbrack_cluster.py @@ -0,0 +1,233 @@ +""" +Tools to assign a structure or a trajectory of structures into +conformational clusters based on Modi and Dunbrack, 2019 (https://pubmed.ncbi.nlm.nih.gov/30867294/) +""" +from pathlib import Path +import tempfile +from typing import Union + +from appdirs import user_cache_dir +import pandas as pd + + +def assign(dihedrals, distances): + from math import cos + import numpy as np + + # define the centroid values for Dunbrack features + centroid = dict() + centroid[(0, "x_phi")] = -129.0 + centroid[(0, "x_psi")] = 179.0 + centroid[(0, "d_phi")] = 61.0 + centroid[(0, "d_psi")] = 81.0 + centroid[(0, "f_phi")] = -97.0 + centroid[(0, "f_psi")] = 20.0 + centroid[(0, "f_chi1")] = -71.0 + + centroid[(1, "x_phi")] = -119.0 + centroid[(1, "x_psi")] = 168.0 + centroid[(1, "d_phi")] = 59.0 + centroid[(1, "d_psi")] = 34.0 + centroid[(1, "f_phi")] = -89.0 + centroid[(1, "f_psi")] = -8.0 + centroid[(1, "f_chi1")] = 56.0 + + centroid[(2, "x_phi")] = -112.0 + centroid[(2, "x_psi")] = -8.0 + centroid[(2, "d_phi")] = -141.0 + centroid[(2, "d_psi")] = 148.0 + centroid[(2, "f_phi")] = -128.0 + centroid[(2, "f_psi")] = 23.0 + centroid[(2, "f_chi1")] = -64.0 + + centroid[(3, "x_phi")] = -135.0 + centroid[(3, "x_psi")] = 175.0 + centroid[(3, "d_phi")] = 60.0 + centroid[(3, "d_psi")] = 65.0 + centroid[(3, "f_phi")] = -79.0 + centroid[(3, "f_psi")] = 145.0 + centroid[(3, "f_chi1")] = -73.0 + + centroid[(4, "x_phi")] = -125.0 + centroid[(4, "x_psi")] = 172.0 + centroid[(4, "d_phi")] = 60.0 + centroid[(4, "d_psi")] = 33.0 + centroid[(4, "f_phi")] = -85.0 + centroid[(4, "f_psi")] = 145.0 + centroid[(4, "f_chi1")] = 49.0 + + centroid[(5, "x_phi")] = -106.0 + centroid[(5, "x_psi")] = 157.0 + centroid[(5, "d_phi")] = 69.0 + centroid[(5, "d_psi")] = 21.0 + centroid[(5, "f_phi")] = -62.0 + centroid[(5, "f_psi")] = 134.0 + centroid[(5, "f_chi1")] = -145.0 + + assignment = list() + for i in range(len(distances)): + ## reproduce the Dunbrack clustering + ## level1: define the DFG positions + if distances[i][0] <= 11.0 and distances[i][1] <= 11.0: # angstroms + ## can only be BABtrans + assignment.append(7) + elif distances[i][0] > 11.0 and distances[i][1] < 14.0: + ## can only be BBAminus + assignment.append(6) + else: + ## belong to DFGin and possibly clusters 0 - 5 + mindist = 10000.0 + cluster_assign = 0 + + for c in range(6): + total_dist = ( + float( + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][0] - centroid[(c, "x_phi")]) * np.pi / 180.0) + ) + ) + + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][1] - centroid[(c, "x_psi")]) * np.pi / 180.0) + ) + ) + + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][2] - centroid[(c, "d_phi")]) * np.pi / 180.0) + ) + ) + + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][3] - centroid[(c, "d_psi")]) * np.pi / 180.0) + ) + ) + + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][4] - centroid[(c, "f_phi")]) * np.pi / 180.0) + ) + ) + + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][5] - centroid[(c, "f_psi")]) * np.pi / 180.0) + ) + ) + + ( + 2.0 + * ( + 1.0 + - cos((dihedrals[i][6] - centroid[(c, "f_chi1")]) * np.pi / 180.0) + ) + ) + ) + / 7 + ) + if total_dist < mindist: + mindist = total_dist + clust_assign = c + assignment.append(clust_assign) + return assignment + + +class PDBDunbrack: + + _PDB_DUNBRACK_LIBRARY = Path(f"{user_cache_dir()}/pdb_dunbrack_library.csv") + + def __init__(self): + self.pdb_dunbrack_library = self.update() + + def __repr__(self): + return f"" + + def update(self, reinitialize: bool = False) -> pd.DataFrame: + """ + Update DataFrame holding information about kinases from the KLIFS database and the corresponding Dunbrack + cluster. + Parameters + ---------- + reinitialize: bool + If the DataFrame should be built from scratch. + Returns + ------- + pdb_dunbrack_library: pd.DataFrame + DataFrame holding information about kinases from KLIFS and the corresponding Dunbrack cluster. + """ + from .klifs import query_klifs_database + import klifs_utils + import MDAnalysis as mda + from .protein_struct_features import key_klifs_residues, compute_simple_protein_features + from tqdm import tqdm + + # get available kinase information from KLIFS + klifs_kinase_ids = klifs_utils.remote.kinases.kinase_names().kinase_ID.to_list() + klifs_kinase_df = klifs_utils.remote.structures.structures_from_kinase_ids( + klifs_kinase_ids + ) + + # initialize library + if not self._PDB_DUNBRACK_LIBRARY.is_file() or reinitialize is True: + columns = list(klifs_kinase_df.columns) + ["dunbrack_cluster"] + pdb_dunbrack_library = pd.DataFrame(columns=columns) + pdb_dunbrack_library.to_csv(self._PDB_DUNBRACK_LIBRARY, index=False) + + pdb_dunbrack_library = pd.read_csv(self._PDB_DUNBRACK_LIBRARY) + + counter = 0 + for index, row in tqdm(klifs_kinase_df.iterrows(), total=klifs_kinase_df.shape[0]): + structure_id = row["structure_ID"] + if structure_id not in list(pdb_dunbrack_library["structure_ID"]): + counter += 1 + try: # assign dunbrack cluster + with tempfile.NamedTemporaryFile(suffix=".pdb", mode="w+t") as temp_file: + pdb_text = klifs_utils.remote.coordinates.complex._complex_pdb_text( + structure_id + ) + temp_file.write(pdb_text) + u = mda.Universe(temp_file.name) + klifs = query_klifs_database(row["pdb"], row["chain"]) + key_res = key_klifs_residues(klifs["numbering"]) + dihedrals, distances = compute_simple_protein_features(u, key_res) + assignment = assign(dihedrals, distances)[0] + except: # catch all errors and assign None + assignment = None + row["dunbrack_cluster"] = assignment + pdb_dunbrack_library = pdb_dunbrack_library.append(row, ignore_index=True) + if counter % 10 == 0: # save every 10th structure, so one can pause in between + pdb_dunbrack_library.to_csv(self._PDB_DUNBRACK_LIBRARY, index=False) + + pdb_dunbrack_library.to_csv(self._PDB_DUNBRACK_LIBRARY, index=False) + return pdb_dunbrack_library + + def structures_by_cluster(self, cluster_id: Union[int, None]) -> pd.DataFrame: + """ + Get KLIFS kinase information of structures matching the given Dunbrack cluster ID. + Parameters + ---------- + cluster_id: int or None + Dunbrack cluser ID of interest. None for structures that were not assigned to any cluster. + Returns + ------- + structures: pd.DataFrame + KLIFS kinase information about matching structures. + """ + if cluster_id is None: + structures = self.pdb_dunbrack_library[ + self.pdb_dunbrack_library["dunbrack_cluster"].isnull() + ] + else: + structures = self.pdb_dunbrack_library[ + self.pdb_dunbrack_library["dunbrack_cluster"] == cluster_id + ] + return structures diff --git a/kinoml/features/kinase.py b/kinoml/features/kinase.py new file mode 100644 index 00000000..32256feb --- /dev/null +++ b/kinoml/features/kinase.py @@ -0,0 +1,74 @@ +""" +kinase_model.py +Defines the Kinase class + +""" + + +class Kinase(object): + def __init__( + self, + pdb, + chain, + kinase_id, + name, + struct_id, + ligand, + pocket_seq, + numbering, + key_res, + dihedrals, + distances, + mean_dist, + ): + """This script defines a Kinase class of which any kinase can be represented as an object with the + following parameters: + + Parameters + ---------- + pdb: str + The PDB code of the structure. + chain: str + The chain index of the structure. + kinase_id: int + The standard ID of a kinase enforced by the KLIFS database. + name: str + The standard name of the kinase used by the KLIFS database. + struct_id: int + The ID associated with a specific chain in the pdb structure of a kinase. + ligand: str + The ligand name as it appears in the pdb file. + pocket_seq: str + The 85 discontinuous residues (from multi-sequence alignment) that define the binding pocket of a kinase. + numbering: list of int + The residue indices of the 85 pocket residues specific to the structure. + key_res: list of int + A list of residue indices that are relevant to the collective variables. + dihedrals: list of floats + A list (one frame) or lists (multiple frames) of dihedrals relevant to kinase conformation. + distances: list of floats + A list (one frame) or lists (multiple frames) of intramolecular distances relevant to kinase conformation. + mean_dist: float + A float (one frame) or a list of floats (multiple frames), which is the mean pairwise distance between + ligand heavy atoms and the CAs of the 85 pocket residues. + + .. todo :: + + This is WAY too many positional arguments. Can we use kwargs instead, or somehow simplify the positional arguments into logical groups? + Many of these will be optional if we want to represent aspects of a structure, so there's no need to make them all requiredself. + Also, we will likely not want to mix features (distances, dihedrals) with structural information directly. + + """ + + self.pdb = pdb + self.chain = chain + self.kinase_id = kinase_id + self.name = name + self.struct_id = struct_id + self.ligand = ligand + self.pocket_seq = pocket_seq + self.numbering = numbering + self.key_res = key_res + self.dihedrals = dihedrals + self.distances = distances + self.mean_dist = mean_dist diff --git a/kinoml/features/klifs.py b/kinoml/features/klifs.py new file mode 100644 index 00000000..1973318e --- /dev/null +++ b/kinoml/features/klifs.py @@ -0,0 +1,141 @@ +""" +Tools for querying the KLIFS database + +http://klifs.vu-compmedchem.nl/ + +""" + +# Setup general logging (guarantee output/error message in case of interruption) +# TODO: Can we log to the terminal instead? +import logging + +logger = logging.getLogger(__name__) +logging.root.setLevel(logging.INFO) +logging.basicConfig(level=logging.INFO, format="%(message)s") +logging.getLogger("urllib3").setLevel(logging.WARNING) + + +def query_klifs_database(pdbid, chainid): + """ + Retrieve KLIFS information from the KLIFTS database. + + Parameters + ---------- + pdbid: str + The PDB code of the inquiry kinase. + chainid: str + The chain index of the inquiry kinase. + + Returns + ------- + klifs : dict + Relevant KLIFS information for the desired pdbid and chain, with the following keys:value pairs + + kinase_id: int + The standard ID of a kinase enforced by the KLIFS database. + name: str + The standard name of the kinase used by the KLIFS database. + pocket_seq: str + The 85 discontinuous residues (from multisequence alignment) that define the binding pocket of a kinase. + struct_id: int + The ID associated with a specific chain in the pdb structure of a kinase. + ligand: str + The ligand name as it appears in the pdb file. + numbering: list of int + The residue indices of the 85 pocket residues specific to the structure. + + .. todo :: Why not have this return a KLIFS object instead? A Python object model for the KLIFS data would likely be very useful. + + """ + import urllib, requests + + # get information of the query kinase from the KLIFS database and gives values + # of kinase_id, name and pocket_seq (numbering) + url = "http://klifs.vu-compmedchem.nl/api/structures_pdb_list?pdb-codes=" + str(pdbid) + + # check to make to sure the search returns valid info + # if return is empty + if len(requests.get(url).text) == 0: + raise ValueError("No data found in KLIFS for pdbid '{}'.".format(pdbid)) + else: + # clean up the info from KLIFS + clean = requests.get(url).text.replace("true", "True").replace("false", "False") + + # each pdb code corresponds to multiple structures + chain_found = False + import ast + + for structure in ast.literal_eval(clean): + numbering = None + # find the specific chain + if isinstance(structure, int): ## if the stucture is not found in the klifs database + kinase_id = None + name = None + pocket_seq = None + struct_id = None + ligand = None + return { + "kinase_id": kinase_id, + "name": name, + "struct_id": struct_id, + "ligand": ligand, + "pocket_seq": pocket_seq, + "numbering": numbering, + } + else: + if structure["chain"] == str(chainid): + kinase_id = int(structure["kinase_ID"]) + name = str(structure["kinase"]) + pocket_seq = str(structure["pocket"]) + struct_id = int(structure["structure_ID"]) + # make sure the specified structure is not an apo structure + ligand = None + if structure["ligand"] != 0: + ligand = str(structure["ligand"]) + chain_found = True + if not chain_found: + raise ValueError( + "No data found for chainid '{}'." + "Please make sure you provide a capital letter (A, B, C, ...) as a chain ID.".format( + chainid + ) + ) + + # Get the numbering of the 85 pocket residues + cmd = "http://klifs.vu-compmedchem.nl/details.php?structure_id=" + str(struct_id) + preload = urllib.request.urlopen(cmd) + info = urllib.request.urlopen(cmd) + for line_number, line in enumerate(info): + line = line.decode() + if "pocketResidues=[" in line: + numbering = ast.literal_eval((line[line.find("=") + 1 : line.find(";")])) + # check if there is gaps/missing residues among the pocket residues. + # If so, enforce their indices as 0 and avoid using them to compute collective variables. + if numbering != None and len(numbering) > 0: + for i in range(len(numbering)): + if numbering[i] == -1: + # logging.info( + # "Warning: There is a gap/missing residue at position: " + + # str(i + 1) + + # ". Its index will be enforced as 0 and it will not be used to compute collective variables." + # ) + numbering[i] = 0 + # print("numbering:") + # print(numbering) + # print out kinase information + # logging.info("Kinase ID: " + str(kinase_id)) + # logging.info("Kinase name: " + str(name)) + # logging.info("Pocket residues: " + str(pocket_seq)) + # logging.info("Structure ID: " + str(struct_id)) + # logging.info("Ligand name: " + str(ligand)) + # logging.info("Numbering of the 85 pocket residues: " + str(numbering)) + + # TODO: Return an object (or potentially a dict) containing this information, rather than just a list of arguments. + return { + "kinase_id": kinase_id, + "name": name, + "struct_id": struct_id, + "ligand": ligand, + "pocket_seq": pocket_seq, + "numbering": numbering, + } diff --git a/kinoml/features/protein_struct_features.py b/kinoml/features/protein_struct_features.py new file mode 100644 index 00000000..d4bfdc2a --- /dev/null +++ b/kinoml/features/protein_struct_features.py @@ -0,0 +1,199 @@ +""" +protein.py +This is a tool to featurize kinase conformational changes through the entire Kinome. + +""" + + +def key_klifs_residues(numbering): + """ + Retrieve a list of PDB residue indices relevant to key kinase conformations mapped via KLIFS. + + Define indices of the residues relevant to a list of 12 collective variables relevant to + kinase conformational changes. These variables include: angle between aC and aE helices, + the key K-E salt bridge, DFG-Phe conformation (two distances), X-DFG-Phi, X-DFG-Psi, + DFG-Asp-Phi, DFG-Asp-Psi, DFG-Phe-Phi, DFG-Phe-Psi, DFG-Phe-Chi1, and the FRET L-S distance. + All features are under the current numbering of the structure provided. + + Parameters + ---------- + numbering : list of int + numbering[klifs_index] is the residue number for the given PDB file corresponding to KLIFS residue index 'klifs_index' + + Returns + ------- + key_res : list of int + Key residue indices + + """ + if numbering == None: + print("The structure was not found in the klifs database.") + key_res = None + return key_res + + key_res = dict() # initialize key_res (which read from the 0-based numbering list) + for i in range(5): + key_res[f"group{i}"] = list() + ## feature group 0: A-loop backbone dihedrals + key_res["group0"].append(numbering[83]) # start of A-loop + + ## feature group 1: P-loop backbone dihedrals + key_res["group1"].append(numbering[3]) # res0 in P-loop + key_res["group1"].append(numbering[4]) # res1 in P-loop + key_res["group1"].append(numbering[5]) # res2 in P-loop + key_res["group1"].append(numbering[6]) # res3 in P-loop + key_res["group1"].append(numbering[7]) # res4 in P-loop + key_res["group1"].append(numbering[8]) # res5 in P-loop + + ## feature group 2: aC-related features + # angle between aC and aE helices and the key salt bridge + key_res["group2"].append(numbering[19]) # res0 in aC + key_res["group2"].append(numbering[29]) # res10 in aC + key_res["group2"].append(numbering[62]) # end of aE + key_res["group2"].append(numbering[16]) # K in beta III + key_res["group2"].append(numbering[23]) # E in aC + + ## feature group 3: DFG-related features + key_res["group3"].append(numbering[79]) # X-DFG + key_res["group3"].append(numbering[80]) # DFG-Asp + key_res["group3"].append(numbering[81]) # DFG-Phe + key_res["group3"].append(numbering[27]) # ExxxX + + ## feature group 4: the FRET distance + # not in the list of 85 (equivalent to Aura"S284"), use the 100% conserved beta III K as a reference + key_res["group4"].append(numbering[16] + 120) + + # not in the list of 85 (equivalent to Aura"L225"), use the 100% conserved beta III K as a reference + key_res["group4"].append(numbering[16] + 61) + + return key_res + + +def compute_simple_protein_features(u, key_res): + """ + This function takes the PDB code, chain id and certain coordinates of a kinase from + a command line and returns its structural features. + + Parameters + ---------- + u : object + A MDAnalysis.core.universe.Universe object of the input structure (a pdb file or a simulation trajectory). + key_res : dict of int + A dictionary (with keys 'group0' ... 'group4') of feature-related residue indices in five feature groups. + Returns + ------- + features: list of floats + A list (single structure) or lists (multiple frames in a trajectory) of 72 features in 5 groups (A-loop, P-loop, aC, DFG, FRET) + + .. todo :: Use kwargs with sensible defaults instead of relying only on positional arguments. + + + """ + from MDAnalysis.core.groups import AtomGroup + from MDAnalysis.analysis.dihedrals import Dihedral + from MDAnalysis.analysis.distances import dist + import numpy as np + import pandas as pd + + # get the array of atom indices for the calculation of: + # * seven dihedrals (a 7*4 array where each row contains indices of the four atoms for each dihedral) + # * two ditances (a 2*2 array where each row contains indices of the two atoms for each dihedral) + dih = np.zeros(shape=(7, 4), dtype=int, order="C") + dis = np.zeros(shape=(2, 2), dtype=int, order="C") + + # name list of the dihedrals and distances + dih_names = ["xDFG_phi", "xDFG_psi", "dFG_phi", "dFG_psi", "DfG_phi", "DfG_psi", "DfG_chi1"] + dis_names = ["DFG_conf1", "DFG_conf2", "DFG_conf3", "DFG_conf4"] + + # parse the topology info (0-based atom indices) + + ### dihedrals (feature group 3) + # dihedral 0 & 1: X-DFG Phi & Psi + dih[0][0] = int(u.select_atoms(f"resid {key_res['group3'][0]-1} and name C")[0].ix) # xxDFG C + dih[0][1] = int(u.select_atoms(f"resid {key_res['group3'][0]} and name N")[0].ix) # xDFG N + dih[0][2] = int(u.select_atoms(f"resid {key_res['group3'][0]} and name CA")[0].ix) # xDFG CA + dih[0][3] = int(u.select_atoms(f"resid {key_res['group3'][0]} and name C")[0].ix) # xDFG C + dih[1][0] = dih[0][1] # xDFG N + dih[1][1] = dih[0][2] # xDFG CA + dih[1][2] = dih[0][3] # xDFG C + dih[1][3] = int(u.select_atoms(f"resid {key_res['group3'][1]} and name N")[0].ix) # DFG-Asp N + + # dihedral 2 & 3: DFG-Asp Phi & Psi + dih[2][0] = dih[0][3] # xDFG C + dih[2][1] = dih[1][3] # DFG-Asp N + dih[2][2] = int( + u.select_atoms(f"resid {key_res['group3'][1]} and name CA")[0].ix + ) # DFG-Asp CA + dih[2][3] = int(u.select_atoms(f"resid {key_res['group3'][1]} and name C")[0].ix) # DFG-Asp C + dih[3][0] = dih[2][1] # DFG-Asp N + dih[3][1] = dih[2][2] # DFG-Asp CA + dih[3][2] = dih[2][3] # DFG-Asp C + dih[3][3] = int(u.select_atoms(f"resid {key_res['group3'][2]} and name N")[0].ix) # DFG-Phe N + + # dihedral 4 & 5: DFG-Phe Phi & Psi + dih[4][0] = dih[2][3] # DFG-Asp C + dih[4][1] = dih[3][3] # DFG-Phe N + dih[4][2] = int( + u.select_atoms(f"resid {key_res['group3'][2]} and name CA")[0].ix + ) # DFG-Phe CA + dih[4][3] = int(u.select_atoms(f"resid {key_res['group3'][2]} and name C")[0].ix) # DFG-Phe C + dih[5][0] = dih[4][1] # DFG-Phe N + dih[5][1] = dih[4][2] # DFG-Phe CA + dih[5][2] = dih[4][3] # DFG-Phe C + dih[5][3] = int( + u.select_atoms(f"resid {key_res['group3'][2]+1} and name N")[0].ix + ) # DFG-Gly N + + # dihedral 6: DFG-Phe Chi1 + dih[6][0] = dih[3][3] # DFG-Phe N + dih[6][1] = dih[4][2] # DFG-Phe CA + dih[6][2] = int( + u.select_atoms(f"resid {key_res['group3'][2]} and name CB")[0].ix + ) # DFG-Phe CB + dih[6][3] = int( + u.select_atoms(f"resid {key_res['group3'][2]} and name CG")[0].ix + ) # DFG-Phe CG + + ### distances + ## Dunbrack distances D1, D2 + dis[0][0] = int(u.select_atoms(f"resid {key_res['group3'][3]} and name CA")[0].ix) # ExxxX CA + dis[0][1] = int( + u.select_atoms(f"resid {key_res['group3'][2]} and name CZ")[0].ix + ) # DFG-Phe CZ + dis[1][0] = int( + u.select_atoms(f"resid {key_res['group2'][3]} and name CA")[0].ix + ) # K in beta III CA + dis[1][1] = dis[0][1] # DFG-Phe CZ + + # check if there is any missing coordinates; if so, skip dihedral/distance calculation for those residues + check_flag = 1 + for i in range(len(dih)): + if 0 in dih[i]: + dih[i] = [0, 0, 0, 0] + check_flag = 0 + + for i in range(len(dis)): + if 0 in dis[i]: + dis[i] = [0, 0] + check_flag = 0 + if check_flag: + print("There is no missing coordinates. All dihedrals and distances will be computed.") + + # compute dihedrals and distances + distances = list() + dih_ags = list() + for i in range(7): # for each of the dihedrals + dih_ags.append(AtomGroup(dih[i], u)) + dihedrals = Dihedral(dih_ags).run().angles + + each_frame = list() + for i in range(2): + ag0 = AtomGroup([dis[i][0]], u) # first atom in each atom pair + ag1 = AtomGroup([dis[i][1]], u) # second atom in each atom pair + each_frame.append(dist(ag0, ag1)[-1][0]) + each_frame = np.array(each_frame) + distances.append(each_frame) + + # clean up + del u, dih, dis + return dihedrals, distances