From 353b66b478b3a06bcbe29bf0a0e599e53df46f1f Mon Sep 17 00:00:00 2001 From: Andrea Volkamer Date: Wed, 18 Aug 2021 12:40:29 +0200 Subject: [PATCH] re-add HM functionalities --- examples/homology.ipynb | 371 +++++++++++++++++++++++++++++++++++ kinoml/modeling/alignment.py | 208 ++++++++++++++++++++ kinoml/modeling/homology.py | 149 ++++++++++++++ 3 files changed, 728 insertions(+) create mode 100644 examples/homology.ipynb create mode 100644 kinoml/modeling/homology.py diff --git a/examples/homology.ipynb b/examples/homology.ipynb new file mode 100644 index 00000000..7e927bc3 --- /dev/null +++ b/examples/homology.ipynb @@ -0,0 +1,371 @@ +{ + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.8-final" + }, + "orig_nbformat": 2, + "kernelspec": { + "name": "python37864bitkinomltestconda53eeb1b28d7e44c880a287b87d4ac530", + "display_name": "Python 3.7.8 64-bit ('kinoml_test': conda)" + } + }, + "nbformat": 4, + "nbformat_minor": 2, + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Homology Modelling" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Author: William Glass\n", + "\n", + "This notebook shows examples of how to perform homology modelling in KinoML. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# get relevant imports\n", + "from kinoml.modeling.homology import HomologyModel\n", + "from kinoml.modeling.alignment import Alignment\n", + "from kinoml.core.proteins import ProteinStructure" + ] + }, + { + "source": [ + "### Basic usage" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To start, we need a template from which to base our homology model on." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "hm = HomologyModel()\n", + "\n", + "# If we already have our model (e.g. if prepared using Spruce TK etc), so we can load it easily\n", + "structure = ProteinStructure.from_file(\"./4yne_protein.pdb\")\n", + "\n", + "# If we just want to grab a structure from the PDB we could use the `from_name` attribute in `ProteinStructure`, e.g:\n", + "structure_from_pdb = ProteinStructure.from_name('4yne')\n", + "\n", + "# Once we have our structure, we need to extract it's sequence. Note, this will often not be the canonical sequence.\n", + "sequence = structure_from_pdb.sequence" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now have our template structure and its sequence. If, for some reason, we didn't have access to the structure but did have access to the sequence we could run a BLAST search to find a PDB structure for our template:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stderr", + "text": "@> Blast searching NCBI PDB database for \"GSGIR...\"\n@> Blast search completed in 47.5s.\n" + } + ], + "source": [ + "model_templates = hm.get_pdb_template(sequence)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "dict_keys(['4yne', '4aoj', '5h3q', '6d1y', '6iqn', '6nss', '5kml', '6nsp', '6npt', '5kmi', '4gt5', '6pl1', '5jfs', '6d22', '5wr7', '4pmm', '5i8a', '5kvt', '4f0i', '4ymj', '3v5q', '6kzc', '4asz', '1luf', '3zos', '5fdp', '6few', '5bvk', '6y23', '6brj', '6tu9', '6fer', '3zzw', '4gt4', '3eta', '1irk', '4ibm', '5hhw', '5e1s', '1p14', '3ekk', '6pxv', '5xff', '5xfj', '5jkg', '4uxq', '1i44', '4tye', '5nud', '4qqj', '6jpe', '4qqt', '2z8c', '1gag', '1rqq', '4qq5', '6iuo', '4qqc', '4xlv', '4k33', '4fnw', '2yjr', '3lco', '3lw0', '1p4o', '6pnx', '4ux0', '1agw', '3rhx', '4zsa', '3kxx', '5a4c', '4rwi', '4f63', '4wun', '3js2', '5zv2', '3c4f', '3gql', '5aa8', '5aa9', '4fnx', '2xb7', '6mx8', '4fnz', '6e0r', '4tt7', '4dce', '5a9u', '2xp2', '3lcs', '2yfx', '3aox', '2yjs', '4z55', '4fob', '2yhv', '3l9p', '4hw7', '6t2w', '4xcu', '6lvm', '3tt0', '1jqh', '3o23', '3qqu', '2oj9', '3i81', '5fxq', '5fxr', '4d2r', '1m7n', '3lvp', '3d94', '6jk8', '6mzw', '6nvl', '5flf', '5vnd', '4anl', '4ans', '2ogv', '5i9u', '1mqb', '3zbf', '4p2k', '5ek7', '4trl', '5a46', '5fxs', '6pyh', '3gqi', '6vg3', '4j99', '2zm3', '1k3a', '6fek', '6i83', '2ivt', '2ivv', '2ivs', '6i82', '6vhg', '5fm2', '4csv', '2pvy', '5ugl', '2pzp', '4j97', '4j96', '4j98', '2pwl', '2pz5', '2rei', '5amn', '6nja', '2pzr', '6lvk', '6lvl', '1gjo', '3ri1', '5ugx', '2q0b', '5ui0', '5uhn', '2py3', '2psq', '2pvf', '6in0', '3b2t', '5eg3', '3kul', '4xli', '3cd3', '3bkb', '6jmf', '2xyn', '5to8', '3et7', '3cc6', '4h1j', '2qoc', '3dzq', '3zfx', '2qob', '2qoo', '5mja', '3zfm', '1jpa', '2xyu', '2y6m', '2hel', '3sxr', '6i99', '3fxx', '2qoi', '2gsf', '2qof', '2qod', '3miy', '3qgw', '4hct', '5mjb', '2hen', '2r2p', '4ot5', '2qok', '2qol', '2qo7', '2qon', '4nwm', '1sm2', '4kio', '3v5j', '3zfy', '3qrj', '2hyy', '2gqg', '2hzi', '3qri', '5hu9', '3t9t', '2hz0', '3dk6', '2g2f', '2f4j', '4xey', '3dk7', '4twp', '3oy3', '6bl8', '2v7a', '2z60', '2hiw', '4zog', '4wa9', '3dk3', '2g1t', '3oxz', '2e2b', '1fpu', '2qoh', '3pyy', '5mo4', '1opk', '2fo0', '1opl', '2yn8', '6fni'])" + }, + "metadata": {}, + "execution_count": 4 + } + ], + "source": [ + "model_templates.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "model_template = list(model_templates.keys())[0]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this toy example our BLAST search using the `4YNE` sequence has, as expected, returned the `4YNE` structure in the PDB as the \"best\" model for us to use as our template.\n", + "\n", + "Typically, we will want to search the PDB with a query sequence and find the most relevant PDB structure to use as our template. We can use the `get_sequence` attribute to obtain the full canonical sequence from a database based on the a unique ID (either using the default `backend=\"uniprot\"` or `backend=\"ncbi\"`). Defaults: `backend=True` and `kinase=True`, the latter ensures we refine our search to sequences that contain a kinase domain." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "uniprot_seq = hm.get_sequence('P04629', kinase=True) # Get sequence based on the UNIPROT ID = P04629. Default: `kinase=True` is shown here for clarity." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then run a blast search using this canonical sequence to search the PDB for the best template model to use." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stderr", + "text": "@> Blast searching NCBI PDB database for \"IVLKW...\"\n@> Blast search completed in 68.7s.\n" + } + ], + "source": [ + "model_templates = hm.get_pdb_template(uniprot_seq) # Need to fix timeouts" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "dict_keys(['4aoj', '5h3q', '5kml', '6d1y', '6iqn', '5kmi', '4gt5', '4pmm', '5i8a', '4yne', '6nsp', '5kvt', '6nss', '4f0i', '5jfs', '6d22', '5wr7', '6npt', '6pl1', '3v5q', '6kzc', '4ymj', '4asz', '1luf', '3zos', '5fdp', '5bvk', '6few', '6y23', '6brj', '6fer', '6tu9', '3eta', '1irk', '4ibm', '5hhw', '5e1s', '1p14', '3ekk', '1i44', '6pxv', '3zzw', '4gt4', '4fnw', '2yjr', '2z8c', '1gag', '1rqq', '4xlv', '5aa9', '5aa8', '6mx8', '2xb7', '4fnz', '4dce', '2xp2', '6e0r', '4tt7', '5a9u', '3aox', '2yfx', '2yjs', '4fob', '2yhv', '3lcs', '4z55', '3l9p', '3lw0', '1p4o', '4fnx', '4ans', '4anl', '6pyh', '3lco', '2oj9', '5fxq', '1jqh', '3o23', '3i81', '5fxr', '3qqu', '4d2r', '3lvp', '1m7n', '3d94', '6jk8', '5xfj', '5xff', '4uxq', '5jkg', '6t2w', '3lcd', '3zbf', '4tye', '5nud', '4qqj', '6jpe', '4qqt', '3gql', '4hw7', '2ogv', '5fxs', '4qq5', '4qqc', '6iuo', '3tt0', '4k33', '2zm3', '1k3a', '5a46', '4zsa', '3rhx', '5a4c', '1agw', '4ux0', '4f63', '4rwi', '4wun', '3kxx', '5zv2', '3js2', '6pnx', '3c4f', '3gqi', '2qoo', '3kul', '6nvl', '6mzw', '5flf', '5vnd', '4hvs', '2ivs', '6i82', '5fm2', '4csv', '4xcu', '2gsf', '3fxx', '2qoi', '2qof', '2qod', '6vg3', '6lvm', '5i9u', '4p2k', '5ek7', '1mqb', '6gqj', '6gqk', '3bea', '2i0v', '6fek', '6i83', '2ivt', '2ivv', '6vhg', '2qok', '2qol', '2qo7', '2qon', '4trl', '2i1m', '5amn', '6nja', '4j96', '4j99', '2pvy', '2pzp', '4j97', '4j98', '2pwl', '2pz5', '2qoc', '6ne7', '4cki', '6in0', '2rei', '3qgw', '5ui0', '2pzr', '3ri1', '6lvl', '1gjo', '6lvk', '2q0b', '5uhn', '2py3', '5ugl', '2pvf', '2qob', '3zfx', '3b2t', '5ugx', '3cly', '2psq', '5mja', '3zfy', '3dzq', '3miy', '3v5j', '1sm2', '4hct', '6agx', '5eg3', '5mjb', '3zfm', '1jpa', '2y6m', '2xyu', '2hel', '3cc6', '3et7', '4kio', '3t9t', '4xli', '2xyn', '2hen', '6jmf', '3bkb', '6i99', '3sxr', '2r2p', '4twp', '3qrj', '3oy3', '2z60', '3qri', '2hyy', '2hz0', '5hu9', '4zog', '2hzi', '3oxz', '2hiw', '3pyy', '2g1t', '2e2b', '4wa9', '1fpu', '2f4j', '2g2f', '2qoh', '3dk6', '4xey', '6bl8', '3dk3', '3dk7', '2v7a', '2gqg'])" + }, + "metadata": {}, + "execution_count": 9 + } + ], + "source": [ + "model_templates.keys()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, `model_templates.keys()` has returned a number of PDB structures we could use as templates for our homology model. " + ] + }, + { + "source": [ + "### Generating the alignment" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A structure (either from a user made one or one downloaded from the PDB with `ProteinStructure.from_name()`) can be used as template on to which we can build our homology model. \n", + "\n", + "We also need a target sequence, this is downloaded from the UniProt server using `HomologyModel().get_sequence()` as shown above.\n", + "\n", + "This information can be used in `Alignment.get_alignment()` to produce an alignment of the two sequences. \n", + "\n", + "The user can also use `Alignment.make_ali_file()` that generates an alignment file based on MODELLER formatting." + ] + }, + { + "source": [ + "#### Generating an alignment and homology model: no ligand" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "source": [ + "In this example we shall use the `4yne` PDB structure as a template and use the Uniprot sequence for NTRK1 to build a homolgy model." + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "structure_from_pdb = ProteinStructure.from_name('4yne') \n", + "uniprot_seq = hm.get_sequence('P04629', kinase=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "output_type": "execute_result", + "data": { + "text/plain": "dict_keys(['score', 'sequence_identity', 'symbols', 'codes'])" + }, + "metadata": {}, + "execution_count": 16 + } + ], + "source": [ + "# Generate the sequence alignment between the template and target sequence\n", + "seq_alignment = Alignment.get_alignment(structure_from_pdb.sequence.sequence, up_seq)\n", + "# To access the alignment information you can use the metadata keys:\n", + "seq_alignment.metadata.keys()" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Now, we can generate an alignment file (based on Modeller)\n", + "seq_alignment.make_ali_file(\n", + " seq_alignment.metadata['symbols'][0], # the aligned template sequence\n", + " seq_alignment.metadata['symbols'][1], # the aligned target sequence\n", + " structure_from_pdb,\n", + " uniprot_seq\n", + " )" + ] + }, + { + "source": [ + "The alignment file can now be used to generate a homolgy model using Modeller. In a real use case the user will need to generate a large number of models and score them (e.g. using DOPE, QMEAN etc). Here, we shall only generate one. " + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "tags": [ + "outputPrepend" + ] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": ".0000 1.000\n 9 Distance restraints 1 (CA-CA) : 0 0 0 0.000 0.000 0.0000 1.000\n10 Distance restraints 2 (N-O) : 0 0 0 0.000 0.000 0.0000 1.000\n11 Mainchain Phi dihedral restraints : 0 0 0 0.000 0.000 0.0000 1.000\n12 Mainchain Psi dihedral restraints : 0 0 0 0.000 0.000 0.0000 1.000\n13 Mainchain Omega dihedral restraints: 36 4 6 0.261 0.261 95.357 1.000\n14 Sidechain Chi_1 dihedral restraints: 28 0 2 1.616 1.616 20.504 1.000\n15 Sidechain Chi_2 dihedral restraints: 21 0 3 1.358 1.358 19.340 1.000\n16 Sidechain Chi_3 dihedral restraints: 14 0 0 1.143 1.143 9.9072 1.000\n17 Sidechain Chi_4 dihedral restraints: 6 0 0 1.778 1.778 5.2542 1.000\n18 Disulfide distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n19 Disulfide angle restraints : 0 0 0 0.000 0.000 0.0000 1.000\n20 Disulfide dihedral angle restraints: 0 0 0 0.000 0.000 0.0000 1.000\n21 Lower bound distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n22 Upper bound distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n23 Distance restraints 3 (SDCH-MNCH) : 0 0 0 0.000 0.000 0.0000 1.000\n24 Sidechain Chi_5 dihedral restraints: 0 0 0 0.000 0.000 0.0000 1.000\n25 Phi/Psi pair of dihedral restraints: 0 0 0 0.000 0.000 0.0000 1.000\n26 Distance restraints 4 (SDCH-SDCH) : 0 0 0 0.000 0.000 0.0000 1.000\n27 Distance restraints 5 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n28 NMR distance restraints 6 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n29 NMR distance restraints 7 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n30 Minimal distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n31 Non-bonded restraints : 10831 0 63 2.683 2.683 -481.17 0.600\n32 Atomic accessibility restraints : 0 0 0 0.000 0.000 0.0000 1.000\n33 Atomic density restraints : 0 0 0 0.000 0.000 0.0000 1.000\n34 Absolute position restraints : 0 0 0 0.000 0.000 0.0000 1.000\n35 Dihedral angle difference restraint: 0 0 0 0.000 0.000 0.0000 1.000\n36 GBSA implicit solvent potential : 0 0 0 0.000 0.000 -416.02 1.000\n37 EM density fitting potential : 0 0 0 0.000 0.000 0.0000 1.000\n38 SAXS restraints : 0 0 0 0.000 0.000 0.0000 1.000\n39 Symmetry restraints : 0 0 0 0.000 0.000 0.0000 1.000\n\n\n\nList of the violated restraints:\n A restraint is violated when the relative difference\n from the best value (RVIOL) is larger than CUTOFF.\n\n ICSR ... index of a restraint in the current set.\n RESNO ... residue numbers of the first two atoms.\n ATM ... IUPAC atom names of the first two atoms.\n FEAT ... the value of the feature in the model.\n restr ... the mean of the basis restraint with the smallest\n difference from the model (local minimum).\n viol ... difference from the local minimum.\n rviol ... relative difference from the local minimum.\n RESTR ... the best value (global minimum).\n VIOL ... difference from the best value.\n RVIOL ... relative difference from the best value.\n\n\n-------------------------------------------------------------------------------------------------\n\nFeature 1 : Bond length potential \nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 36 28D 28D CB CA 217 216 1.73 1.54 0.20 5.38 1.54 0.20 5.38\n 2 42 28D 28D C CA 221 216 1.93 1.49 0.44 12.74 1.49 0.44 12.74\n 3 43 28D 29K C N 221 223 1.91 1.35 0.57 20.17 1.35 0.57 20.17\n 4 49 29K 29K N CA 223 224 2.02 1.43 0.59 19.33 1.43 0.59 19.33\n 5 51 29K 29K C CA 230 224 2.16 1.49 0.67 19.50 1.49 0.67 19.50\n 6 52 29K 30M C N 230 232 1.77 1.35 0.42 14.89 1.35 0.42 14.89\n 7 140 164R 164R CD CG 1258 1257 1.78 1.53 0.25 6.96 1.53 0.25 6.96\n 8 141 164R 164R NE CD 1259 1258 1.70 1.49 0.21 6.10 1.49 0.21 6.10\n\n-------------------------------------------------------------------------------------------------\n\nFeature 2 : Bond angle potential \nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 280 25P 25P N CA 190 191 135.55 108.20 27.35 6.21 108.20 27.35 6.21\n 2 317 28D 28D CA CB 216 217 132.78 108.00 24.78 5.74 108.00 24.78 5.74\n 3 319 28D 28D CA C 216 221 137.08 116.50 20.58 5.92 116.50 20.58 5.92\n 4 323 28D 28D C CA 221 216 134.17 108.00 26.17 6.06 108.00 26.17 6.06\n 5 324 28D 29K C N 221 223 150.19 120.00 30.19 6.86 120.00 30.19 6.86\n 6 327 29K 29K N CA 223 224 143.42 107.00 36.42 10.47 107.00 36.42 10.47\n 7 330 29K 29K CA C 224 230 151.43 116.50 34.93 10.04 116.50 34.93 10.04\n 8 457 164R 164R CB CG 1256 1257 135.02 113.60 21.42 5.26 113.60 21.42 5.26\n 9 458 164R 164R CG CD 1257 1258 143.15 107.50 35.65 9.42 107.50 35.65 9.42\n 10 463 164R 164R C CA 1263 1255 140.82 108.00 32.82 7.60 108.00 32.82 7.60\n\n-------------------------------------------------------------------------------------------------\n\nFeature 13 : Mainchain Omega dihedral restraints \nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 999 163S 163S CA C 1249 1252 -152.81 180.00 27.19 5.44 180.00 27.19 5.44\n\n\nreport______> Distribution of short non-bonded contacts:\n\n serious non-bonded atom clash: 149 222 2.223\n\nDISTANCE1: 0.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40\nDISTANCE2: 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50\nFREQUENCY: 0 0 1 0 6 15 30 70 145 196 214 247 328 350 347\n\n\n<< end of ENERGY.\niupac_m_485_> OD1/2 will be swapped: 120.2134 28 28\niupac_m_486_> OE1/2 will be swapped: -133.7539 39 39\niupac_m_486_> OE1/2 will be swapped: -178.8265 42 42\niupac_m_487_> NH1/2 swapped: 179.9657 45 45\niupac_m_485_> OD1/2 will be swapped: 122.0238 165 165\niupac_m_483_> CD1/2 CE1/2 swapped: 96.0942 167 167\niupac_m_397W> Atoms were not swapped because of the uncertainty of how to handle the H atom.\niupac_m_483_> CD1/2 CE1/2 swapped: -90.8451 172 172\niupac_m_397W> Atoms were not swapped because of the uncertainty of how to handle the H atom.\niupac_m_484W> Dihedral still outside +-90: 93.9427\nopenf___224_> Open ${LIB}/atmcls-mf.lib\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10650794 10401.166 10.157\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10650858 10401.229 10.157\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10650916 10401.285 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10650980 10401.348 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10651046 10401.412 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10651110 10401.475 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10651266 10401.627 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10651394 10401.752 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10651688 10402.039 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10651880 10402.227 10.158\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10652384 10402.719 10.159\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10652640 10402.969 10.159\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10653636 10403.941 10.160\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10654020 10404.316 10.160\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10654596 10404.879 10.161\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10656630 10406.865 10.163\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10657526 10407.740 10.164\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10661978 10412.088 10.168\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10663322 10413.400 10.169\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10672940 10422.793 10.179\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10674924 10424.730 10.180\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10695570 10444.893 10.200\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10698578 10447.830 10.203\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10744544 10492.719 10.247\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10749024 10497.094 10.251\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10850244 10595.941 10.348\nopenf___224_> Open ${LIB}/dist-mf.lib\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10850640 10596.328 10.348\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10851436 10597.105 10.349\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10852232 10597.883 10.349\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10853420 10599.043 10.351\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10853820 10599.434 10.351\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10854420 10600.020 10.352\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10855320 10600.898 10.352\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10856668 10602.215 10.354\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10856868 10602.410 10.354\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10858892 10604.387 10.356\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10859192 10604.680 10.356\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10862228 10607.645 10.359\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10862676 10608.082 10.359\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10867228 10612.527 10.364\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10867900 10613.184 10.364\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10874728 10619.852 10.371\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10875736 10620.836 10.372\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10885980 10630.840 10.382\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10886580 10631.426 10.382\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10886980 10631.816 10.383\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10888492 10633.293 10.384\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10903856 10648.297 10.399\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10904756 10649.176 10.400\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10905356 10649.762 10.400\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10907624 10651.977 10.402\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10930672 10674.484 10.424\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10932016 10675.797 10.426\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10932916 10676.676 10.426\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10936320 10680.000 10.430\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10970892 10713.762 10.463\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10972908 10715.730 10.465\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10974256 10717.047 10.466\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 10979360 10722.031 10.471\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11031216 10772.672 10.520\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11034240 10775.625 10.523\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11036264 10777.602 10.525\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11043920 10785.078 10.532\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11121704 10861.039 10.606\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11126240 10865.469 10.611\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11129276 10868.434 10.614\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11140760 10879.648 10.625\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11257436 10993.590 10.736\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11264240 11000.234 10.742\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11268792 11004.680 10.747\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11286020 11021.504 10.763\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11461036 11192.418 10.930\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11471248 11202.391 10.940\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11478076 11209.059 10.946\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11503916 11234.293 10.971\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11766440 11490.664 11.221\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11781752 11505.617 11.236\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11791996 11515.621 11.246\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 11830756 11553.473 11.283\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12224540 11938.027 11.658\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12247508 11960.457 11.680\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12262872 11975.461 11.695\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12321012 12032.238 11.750\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12911688 12609.070 12.314\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12946140 12642.715 12.346\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 12969188 12665.223 12.368\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 13056400 12750.391 12.452\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 13942416 13615.641 13.297\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 13994100 13666.113 13.346\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 14028672 13699.875 13.379\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 14159488 13827.625 13.504\n\nDynamically allocated memory at amaxgroup_restr [B,KiB,MiB]: 15488512 15125.500 14.771\nread_pa_232_> parameters BONDS ANGLS DIHEDS IMPROPS MODE\n 0 0 0 0 12561\n>> Model assessment by DOPE potential\niatmcls_286W> MODEL atom not classified: LEU:OXT LEU\n\nDynamically allocated memory at amaxrestraints [B,KiB,MiB]: 15990084 15615.316 15.249\n\nDynamically allocated memory at amaxrestraints [B,KiB,MiB]: 16407684 16023.129 15.648\n\n\n>> ENERGY; Differences between the model's features and restraints:\nNumber of all residues in MODEL : 272\nNumber of all, selected real atoms : 2175 257\nNumber of all, selected pseudo atoms : 0 0\nNumber of all static, selected restraints : 1084 1084\nCOVALENT_CYS : F\nNONBONDED_SEL_ATOMS : 1\nNumber of non-bonded pairs (excluding 1-2,1-3,1-4): 64807\nDynamic pairs routine : 5, NATM_SEL x NATM double loop\nAtomic shift for contacts update (UPDATE_DYNAMIC) : 0.390\nLENNARD_JONES_SWITCH : 6.500 7.500\nCOULOMB_JONES_SWITCH : 6.500 7.500\nRESIDUE_SPAN_RANGE : 1 9999\nNLOGN_USE : 15\nCONTACT_SHELL : 15.000\nDYNAMIC_PAIRS,_SPHERE,_COULOMB,_LENNARD,_MODELLER : T F F F T\nSPHERE_STDV : 0.050\nRADII_FACTOR : 0.820\nCurrent energy : -2364.1306\n\n\n\n\n<< end of ENERGY.\niatmcls_286W> MODEL atom not classified: LEU:OXT LEU\nDOPE score : -2364.130615\n>> Model assessment by GA341 potential\n\nSurface library : /Users/williamglass/opt/anaconda3/envs/kinoml_test/lib/modeller-9.24/modlib/surf5.de\nPair library : /Users/williamglass/opt/anaconda3/envs/kinoml_test/lib/modeller-9.24/modlib/pair9.de\nChain identifier : _\n% sequence identity : 100.000000\nSequence length : 272\nCompactness : 0.345183\nNative energy (pair) : -160.943322\nNative energy (surface) : -10.143724\nNative energy (combined) : -5.629863\nZ score (pair) : -6.208156\nZ score (surface) : -6.076228\nZ score (combined) : -8.785641\nGA341 score : 1.000000\n\nopenf___224_> Open P04629.BL00010001.pdb\nwrpdb___568_> Residues, atoms, selected atoms: 272 2175 2175\n\n>> Summary of successfully produced models:\nFilename molpdf\n----------------------------------------\nP04629.B99990001.pdb 2008.34644\n\n\n>> Summary of successfully produced loop models:\nFilename molpdf DOPE score GA341 score\n----------------------------------------------------------------------\nP04629.BL00010001.pdb 725.45709 -2364.13062 1.00000\n\n" + } + ], + "source": [ + "hm.make_model(structure_from_pdb, uniprot_seq, seq_alignment, num_models=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Generating an alignment and homology model: with ligand" + ] + }, + { + "source": [ + "Often, structures will contain the coordinates of bound ligands. If these are of interest and need to be retained this is easy to do:" + ], + "cell_type": "markdown", + "metadata": {} + }, + { + "source": [ + "# We generate the same alignment file (based on Modeller), but now retain the bound ligand\n", + "seq_alignment.make_ali_file(\n", + " seq_alignment.metadata['symbols'][0], # the aligned template sequence\n", + " seq_alignment.metadata['symbols'][1], # the aligned target sequence\n", + " structure_from_pdb,\n", + " uniprot_seq,\n", + " ligand=True\n", + " )" + ], + "cell_type": "code", + "metadata": {}, + "execution_count": 19, + "outputs": [] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "tags": [ + "outputPrepend" + ] + }, + "outputs": [ + { + "output_type": "stream", + "name": "stdout", + "text": " 24.805 1.000\n14 Sidechain Chi_1 dihedral restraints: 29 0 1 1.336 1.336 10.207 1.000\n15 Sidechain Chi_2 dihedral restraints: 22 0 0 1.424 1.424 14.724 1.000\n16 Sidechain Chi_3 dihedral restraints: 14 0 0 1.901 1.901 11.528 1.000\n17 Sidechain Chi_4 dihedral restraints: 6 0 0 1.474 1.474 2.8384 1.000\n18 Disulfide distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n19 Disulfide angle restraints : 0 0 0 0.000 0.000 0.0000 1.000\n20 Disulfide dihedral angle restraints: 0 0 0 0.000 0.000 0.0000 1.000\n21 Lower bound distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n22 Upper bound distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n23 Distance restraints 3 (SDCH-MNCH) : 0 0 0 0.000 0.000 0.0000 1.000\n24 Sidechain Chi_5 dihedral restraints: 0 0 0 0.000 0.000 0.0000 1.000\n25 Phi/Psi pair of dihedral restraints: 0 0 0 0.000 0.000 0.0000 1.000\n26 Distance restraints 4 (SDCH-SDCH) : 0 0 0 0.000 0.000 0.0000 1.000\n27 Distance restraints 5 (X-Y) : 2 0 0 0.000 0.000 0.0000 1.000\n28 NMR distance restraints 6 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n29 NMR distance restraints 7 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n30 Minimal distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n31 Non-bonded restraints : 12747 0 152 2.725 2.725 -143.55 0.600\n32 Atomic accessibility restraints : 0 0 0 0.000 0.000 0.0000 1.000\n33 Atomic density restraints : 0 0 0 0.000 0.000 0.0000 1.000\n34 Absolute position restraints : 0 0 0 0.000 0.000 0.0000 1.000\n35 Dihedral angle difference restraint: 0 0 0 0.000 0.000 0.0000 1.000\n36 GBSA implicit solvent potential : 0 0 0 0.000 0.000 -423.62 1.000\n37 EM density fitting potential : 0 0 0 0.000 0.000 0.0000 1.000\n38 SAXS restraints : 0 0 0 0.000 0.000 0.0000 1.000\n39 Symmetry restraints : 0 0 0 0.000 0.000 0.0000 1.000\n\n\n\nList of the violated restraints:\n A restraint is violated when the relative difference\n from the best value (RVIOL) is larger than CUTOFF.\n\n ICSR ... index of a restraint in the current set.\n RESNO ... residue numbers of the first two atoms.\n ATM ... IUPAC atom names of the first two atoms.\n FEAT ... the value of the feature in the model.\n restr ... the mean of the basis restraint with the smallest\n difference from the model (local minimum).\n viol ... difference from the local minimum.\n rviol ... relative difference from the local minimum.\n RESTR ... the best value (global minimum).\n VIOL ... difference from the best value.\n RVIOL ... relative difference from the best value.\n\n\n-------------------------------------------------------------------------------------------------\n\nFeature 2 : Bond angle potential \nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 454 163S 163S N CA 1248 1249 127.18 107.00 20.18 5.80 107.00 20.18 5.80\n\n\nreport______> Distribution of short non-bonded contacts:\n\n\nDISTANCE1: 0.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40\nDISTANCE2: 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50\nFREQUENCY: 0 0 0 0 1 17 45 122 178 262 216 278 349 379 424\n\n\n<< end of ENERGY.\nopenf___224_> Open P04629.IL00000001.pdb\nwrpdb___568_> Residues, atoms, selected atoms: 273 2202 2202\nopenf___224_> Open P04629.IL00000001.pdb\n\nDynamically allocated memory at amaxmodel [B,KiB,MiB]: 17402241 16994.375 16.596\nread_mo_297_> Segments, residues, atoms: 2 273 2202\nread_mo_298_> Segment: 1 1 A 272 A 2175\nread_mo_298_> Segment: 2 273 B 273 B 27\niatmcls_286W> MODEL atom not classified: LEU:OXT LEU\nrandomi_498_> Atoms,selected atoms,random_seed,amplitude: 2202 266 1 5.0000\nrandomi_496_> Amplitude is > 0; randomization is done.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\n\nDynamically allocated memory at amaxrestraints [B,KiB,MiB]: 17507913 17097.570 16.697\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\npreppdf_458W> Both Lennard-Jones and statistical potential terms selected.\n\n\n>> ENERGY; Differences between the model's features and restraints:\nNumber of all residues in MODEL : 273\nNumber of all, selected real atoms : 2202 266\nNumber of all, selected pseudo atoms : 0 0\nNumber of all static, selected restraints : 1125 1125\nCOVALENT_CYS : F\nNONBONDED_SEL_ATOMS : 1\nNumber of non-bonded pairs (excluding 1-2,1-3,1-4): 10490\nDynamic pairs routine : 5, NATM_SEL x NATM double loop\nAtomic shift for contacts update (UPDATE_DYNAMIC) : 0.390\nLENNARD_JONES_SWITCH : 6.500 7.500\nCOULOMB_JONES_SWITCH : 6.500 7.500\nRESIDUE_SPAN_RANGE : 0 99999\nNLOGN_USE : 15\nCONTACT_SHELL : 8.000\nDYNAMIC_PAIRS,_SPHERE,_COULOMB,_LENNARD,_MODELLER : T F F T T\nSPHERE_STDV : 0.050\nRADII_FACTOR : 0.820\nCurrent energy : -43.0811\n\n\n\n\n\nSummary of the restraint violations: \n\n NUM ... number of restraints.\n NUMVI ... number of restraints with RVIOL > VIOL_REPORT_CUT[i].\n RVIOL ... relative difference from the best value.\n NUMVP ... number of restraints with -Ln(pdf) > VIOL_REPORT_CUT2[i].\n RMS_1 ... RMS(feature, minimally_violated_basis_restraint, NUMB).\n RMS_2 ... RMS(feature, best_value, NUMB).\n MOL.PDF ... scaled contribution to -Ln(Molecular pdf).\n\n # RESTRAINT_GROUP NUM NUMVI NUMVP RMS_1 RMS_2 MOL.PDF S_i\n------------------------------------------------------------------------------------------------------\n 1 Bond length potential : 273 5 8 0.037 0.037 99.309 1.000\n 2 Bond angle potential : 372 7 16 0.109 0.109 302.53 1.000\n 3 Stereochemical cosine torsion poten: 265 0 18 1.169 1.169 107.41 1.000\n 4 Stereochemical improper torsion pot: 105 1 3 0.062 0.062 28.496 1.000\n 5 Soft-sphere overlap restraints : 0 0 0 0.000 0.000 0.0000 1.000\n 6 Lennard-Jones 6-12 potential : 10490 0 43 2.550 2.550 480.94 1.000\n 7 Coulomb point-point electrostatic p: 0 0 0 0.000 0.000 0.0000 1.000\n 8 H-bonding potential : 0 0 0 0.000 0.000 0.0000 1.000\n 9 Distance restraints 1 (CA-CA) : 0 0 0 0.000 0.000 0.0000 1.000\n10 Distance restraints 2 (N-O) : 0 0 0 0.000 0.000 0.0000 1.000\n11 Mainchain Phi dihedral restraints : 0 0 0 0.000 0.000 0.0000 1.000\n12 Mainchain Psi dihedral restraints : 0 0 0 0.000 0.000 0.0000 1.000\n13 Mainchain Omega dihedral restraints: 37 0 1 0.106 0.106 16.216 1.000\n14 Sidechain Chi_1 dihedral restraints: 29 0 1 1.260 1.260 14.721 1.000\n15 Sidechain Chi_2 dihedral restraints: 22 0 4 0.982 0.982 14.266 1.000\n16 Sidechain Chi_3 dihedral restraints: 14 0 0 1.336 1.336 11.121 1.000\n17 Sidechain Chi_4 dihedral restraints: 6 0 0 1.525 1.525 5.9976 1.000\n18 Disulfide distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n19 Disulfide angle restraints : 0 0 0 0.000 0.000 0.0000 1.000\n20 Disulfide dihedral angle restraints: 0 0 0 0.000 0.000 0.0000 1.000\n21 Lower bound distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n22 Upper bound distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n23 Distance restraints 3 (SDCH-MNCH) : 0 0 0 0.000 0.000 0.0000 1.000\n24 Sidechain Chi_5 dihedral restraints: 0 0 0 0.000 0.000 0.0000 1.000\n25 Phi/Psi pair of dihedral restraints: 0 0 0 0.000 0.000 0.0000 1.000\n26 Distance restraints 4 (SDCH-SDCH) : 0 0 0 0.000 0.000 0.0000 1.000\n27 Distance restraints 5 (X-Y) : 2 0 0 0.204 0.204 0.61168 1.000\n28 NMR distance restraints 6 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n29 NMR distance restraints 7 (X-Y) : 0 0 0 0.000 0.000 0.0000 1.000\n30 Minimal distance restraints : 0 0 0 0.000 0.000 0.0000 1.000\n31 Non-bonded restraints : 10436 0 62 2.658 2.658 -702.71 0.600\n32 Atomic accessibility restraints : 0 0 0 0.000 0.000 0.0000 1.000\n33 Atomic density restraints : 0 0 0 0.000 0.000 0.0000 1.000\n34 Absolute position restraints : 0 0 0 0.000 0.000 0.0000 1.000\n35 Dihedral angle difference restraint: 0 0 0 0.000 0.000 0.0000 1.000\n36 GBSA implicit solvent potential : 0 0 0 0.000 0.000 -422.00 1.000\n37 EM density fitting potential : 0 0 0 0.000 0.000 0.0000 1.000\n38 SAXS restraints : 0 0 0 0.000 0.000 0.0000 1.000\n39 Symmetry restraints : 0 0 0 0.000 0.000 0.0000 1.000\n\n\n\nList of the violated restraints:\n A restraint is violated when the relative difference\n from the best value (RVIOL) is larger than CUTOFF.\n\n ICSR ... index of a restraint in the current set.\n RESNO ... residue numbers of the first two atoms.\n ATM ... IUPAC atom names of the first two atoms.\n FEAT ... the value of the feature in the model.\n restr ... the mean of the basis restraint with the smallest\n difference from the model (local minimum).\n viol ... difference from the local minimum.\n rviol ... relative difference from the local minimum.\n RESTR ... the best value (global minimum).\n VIOL ... difference from the best value.\n RVIOL ... relative difference from the best value.\n\n\n-------------------------------------------------------------------------------------------------\n\nFeature 1 : Bond length potential \nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 72 40A 40A CB CA 309 308 1.72 1.54 0.18 5.02 1.54 0.18 5.02\n 2 73 40A 40A N CA 307 308 1.62 1.43 0.19 6.16 1.43 0.19 6.16\n 3 75 40A 40A C CA 310 308 1.86 1.49 0.37 10.66 1.49 0.37 10.66\n 4 76 40A 41S C N 310 312 1.55 1.35 0.21 7.33 1.35 0.21 7.33\n 5 85 42E 42E CD CG 322 321 1.75 1.52 0.23 5.92 1.52 0.23 5.92\n\n-------------------------------------------------------------------------------------------------\n\nFeature 2 : Bond angle potential \nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 275 23L 24L C N 180 182 141.46 120.00 21.46 4.88 120.00 21.46 4.88\n 2 278 24L 24L N CA 182 183 124.07 107.00 17.07 4.91 107.00 17.07 4.91\n 3 370 39E 40A C N 305 307 143.17 120.00 23.17 5.26 120.00 23.17 5.26\n 4 392 42E 42E CB CG 320 321 141.08 108.00 33.08 7.66 108.00 33.08 7.66\n\n-------------------------------------------------------------------------------------------------\n\nFeature 4 : Stereochemical improper torsion potentia\nList of the RVIOL violations larger than : 4.5000\n\n # ICSR RESNO1/2 ATM1/2 INDATM1/2 FEAT restr viol rviol RESTR VIOL RVIOL\n 1 646 23L 23L C CA 180 175 14.76 0.00 14.76 5.20 0.00 14.76 5.20\n\n\nreport______> Distribution of short non-bonded contacts:\n\n serious non-bonded atom clash: 288 310 2.190\n\nDISTANCE1: 0.00 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40\nDISTANCE2: 2.10 2.20 2.30 2.40 2.50 2.60 2.70 2.80 2.90 3.00 3.10 3.20 3.30 3.40 3.50\nFREQUENCY: 0 1 0 5 1 18 37 87 151 222 221 278 323 365 424\n\n\n<< end of ENERGY.\niupac_m_486_> OE1/2 will be swapped: 176.8874 26 26\niupac_m_485_> OD1/2 will be swapped: -116.1402 28 28\niupac_m_486_> OE1/2 will be swapped: -159.3526 39 39\niupac_m_486_> OE1/2 will be swapped: -177.0422 42 42\niupac_m_487_> NH1/2 swapped: 177.0632 164 164\niupac_m_487_> NH1/2 swapped: 178.8781 173 173\n>> Model assessment by DOPE potential\niatmcls_286W> MODEL atom not classified: LEU:OXT LEU\n\nDynamically allocated memory at amaxrestraints [B,KiB,MiB]: 18012077 17589.918 17.178\n\nDynamically allocated memory at amaxrestraints [B,KiB,MiB]: 18434861 18002.793 17.581\n\n\n>> ENERGY; Differences between the model's features and restraints:\nNumber of all residues in MODEL : 273\nNumber of all, selected real atoms : 2202 266\nNumber of all, selected pseudo atoms : 0 0\nNumber of all static, selected restraints : 1125 1125\nCOVALENT_CYS : F\nNONBONDED_SEL_ATOMS : 1\nNumber of non-bonded pairs (excluding 1-2,1-3,1-4): 61556\nDynamic pairs routine : 5, NATM_SEL x NATM double loop\nAtomic shift for contacts update (UPDATE_DYNAMIC) : 0.390\nLENNARD_JONES_SWITCH : 6.500 7.500\nCOULOMB_JONES_SWITCH : 6.500 7.500\nRESIDUE_SPAN_RANGE : 1 9999\nNLOGN_USE : 15\nCONTACT_SHELL : 15.000\nDYNAMIC_PAIRS,_SPHERE,_COULOMB,_LENNARD,_MODELLER : T F F F T\nSPHERE_STDV : 0.050\nRADII_FACTOR : 0.820\nCurrent energy : -3064.3940\n\n\n\n\n<< end of ENERGY.\niatmcls_286W> MODEL atom not classified: LEU:OXT LEU\nDOPE score : -3064.394043\n>> Model assessment by GA341 potential\n\nSurface library : /Users/williamglass/opt/anaconda3/envs/kinoml_test/lib/modeller-9.24/modlib/surf5.de\nPair library : /Users/williamglass/opt/anaconda3/envs/kinoml_test/lib/modeller-9.24/modlib/pair9.de\nChain identifier : _\n% sequence identity : 95.539001\nSequence length : 272\nCompactness : 0.346651\nNative energy (pair) : -175.982511\nNative energy (surface) : -10.258501\nNative energy (combined) : -6.008461\nZ score (pair) : -6.470797\nZ score (surface) : -6.205995\nZ score (combined) : -8.523598\nGA341 score : 1.000000\n\nopenf___224_> Open P04629.BL00010001.pdb\nwrpdb___568_> Residues, atoms, selected atoms: 273 2202 2202\n\n>> Summary of successfully produced models:\nFilename molpdf\n----------------------------------------\nP04629.B99990001.pdb 2415.46313\n\n\n>> Summary of successfully produced loop models:\nFilename molpdf DOPE score GA341 score\n----------------------------------------------------------------------\nP04629.BL00010001.pdb -43.08109 -3064.39404 1.00000\n\n" + } + ], + "source": [ + "hm.make_model(structure_from_pdb, uniprot_seq, seq_alignment, num_models=1, ligand=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ] +} \ No newline at end of file diff --git a/kinoml/modeling/alignment.py b/kinoml/modeling/alignment.py index ae9e5225..1e465a02 100644 --- a/kinoml/modeling/alignment.py +++ b/kinoml/modeling/alignment.py @@ -1,3 +1,211 @@ +from ..core.proteins import ProteinStructure +from ..core.sequences import KinaseDomainAminoAcidSequence +from typing import Union + + +class Alignment: + """Alignment representation of protein sequences""" + + def __init__(self, metadata=None, alignment=None, alignment_file_path=None, *args, **kwargs): + + from appdirs import user_cache_dir + + if metadata is None: + metadata = {} + self.metadata = metadata + self.alignment = alignment + self.alignment_file_path = f"{user_cache_dir()}/alignment.ali" + + @classmethod + def get_alignment(cls, seq1: str, seq2: str, local: bool = True): + """ + Generate an alignment between two sequences + + Parameters + ---------- + seq1: str + The first sequence to be aligned + seq1: str + The second sequence to be aligned + local: bool + If false, a global alignment is performed + (based on the Needleman-Wunsch algorithm), + otherwise a local alignment is performed + (based on the Smith–Waterman algorithm). + (Default: True) + + Returns + ------- + Alignment + """ + + import biotite.sequence as seq + import biotite.sequence.align as align + import numpy as np + + # create the default matrix + # TODO add more options for the choice of matrix + matrix = align.SubstitutionMatrix.std_protein_matrix() + + alignments = align.align_optimal( + seq.ProteinSequence(seq1), + seq.ProteinSequence(seq2), + matrix, + local=local, + ) + + alignment = alignments[0] + + score = alignment.score + seq_identity = align.get_sequence_identity(alignment) + symbols = align.get_symbols(alignment) + codes = align.get_codes(alignment) + + return cls( + alignment=alignment, + metadata={ + "score": score, + "sequence_identity": seq_identity, + "symbols": symbols, + "codes": codes, + }, + ) + + def _format_alignment_with_ligands(self, template, aligned_template_seq): + + from MDAnalysis.lib.util import convert_aa_code + + sel = template.universe.select_atoms("protein and name CA") + resnames = sel.resnames.tolist() + resnames_one_letter = [convert_aa_code(x) for x in resnames] + + # get the last four AAs in the aligned target sequence + end_chunk = aligned_template_seq[-4:][::-1] + aa_store = [] + # loop over residues in reverse + for i, letter in enumerate(resnames_one_letter[::-1]): + # read a chunk of four AAs + current_chunk = resnames_one_letter[::-1][i : i + 4] + + # if the end of the aligned target sequence is hit, stop + if current_chunk == end_chunk: + break + # otherwise store the one letter AAs that need to be added + else: + aa_store.append(letter) + + # create a one element list of dashes + dashes = "".join(["-" for x in aa_store]) + # create one element list of AAs in the canonical order + amino_acids = "".join(aa_store[::-1]) + + return amino_acids, dashes + + def make_ali_file( + self, + aligned_template_seq: str, + aligned_target_seq: str, + template: ProteinStructure, + target: Union[str, KinaseDomainAminoAcidSequence], + ligand: bool = False, + ): + """ + Generate an alignment file in MODELLER format + + Parameters + ---------- + aligned_template_seq: str + The aligned template sequence + aligned_target_seq: str + The aligned target sequence + template: ProteinStructure + The template structure + target: list of str or KinaseDomainAminoAcidSequence + The original target sequence + ligand: bool + Specify whether to retain a ligand in the alignment. + (Default: False) + """ + + # Convert None entries into dashes + seq1_dashed = [i or "-" for i in aligned_template_seq] + seq2_dashed = [i or "-" for i in aligned_target_seq] + + # Setup formatting for MODELLER alignment file + max_length = 75 + + # TODO handle if using backend='ncbi' + + # handle if target is KinaseDomainAminoAcidSequence vs. str + try: + sequence_id = getattr(target, "metadata")["uniprot_id"] + sequence_begin = getattr(target, "metadata")["begin"] + sequence_end = getattr(target, "metadata")["end"] + except: + sequence_id = "sequence_id" + sequence_begin = "1" + sequence_end = len(target) + + protein_id = template.metadata["id"] + # get index for the first residue in the structure sequence + # ignoring negative and 0 resids + resids = template.universe.residues.resids + protein_start_index = next( + (i for i, x in enumerate(template.universe.residues.resids) if x > 0), + None, + ) + protein_start = resids[protein_start_index] + protein_end = "" + + if ligand: + # extract the AAs and No. dashes that need to be added in the alignment + tail_amino_acids, tail_dashes = self._format_alignment_with_ligands( + template, aligned_template_seq + ) + + # add on bufffer AAs / dashes so that Modeller finds the ligand correctly + seq1_dashed += tail_amino_acids + seq2_dashed += tail_dashes + + # write alignment file in MODELLER format + with open(self.alignment_file_path, "w") as ali_file: + for i in range(len(seq1_dashed)): + if i == 0: + ali_file.write(f">P1;{protein_id}\n") + ali_file.write( + f"structure:{protein_id}:{protein_start}:A:{protein_end}: ::: : \n" + ) + ali_file.write(seq1_dashed[i]) + if i == len(seq1_dashed) - 1: + if ligand: + ali_file.write("/.*") + else: + ali_file.write("*") + if (i + 1) % max_length == 0: + ali_file.write("\n") + + for i in range(len(seq2_dashed)): + # start new line below first sequence + if i == 0: + ali_file.write(f"\n>P1;{sequence_id}\n") + if ligand: + ali_file.write( + f"sequence:{sequence_id}:{sequence_begin}: :{sequence_end + len(tail_dashes)}: ::: : \n" + ) + else: + ali_file.write( + f"sequence:{sequence_id}:{sequence_begin}: :{sequence_end}: ::: : \n" + ) + ali_file.write(seq2_dashed[i]) + if i == len(seq2_dashed) - 1: + if ligand: + ali_file.write("/.*") + else: + ali_file.write("*") + if (i + 1) % max_length == 0: + ali_file.write("\n") + + def sequence_similarity( sequence1: str, sequence2: str, diff --git a/kinoml/modeling/homology.py b/kinoml/modeling/homology.py new file mode 100644 index 00000000..a03d6163 --- /dev/null +++ b/kinoml/modeling/homology.py @@ -0,0 +1,149 @@ +"""WIP""" +from ..core.proteins import ProteinStructure +from ..core.sequences import KinaseDomainAminoAcidSequence +from .alignment import Alignment + + +class HomologyModel: # TODO inherent a Base class? + def __init__(self, metadata=None, *args, **kwargs): + # TODO specify id, template, and sequence here? e.g.: + + if metadata is None: + metadata = {} + self.metadata = metadata + + # self.identifier = identifier + # self.template = template + # self.sequence = sequence + + def get_pdb_template( + self, + sequence, + ): + """ + Retrieve a template structure from PDB from a BLAST search + + Parameters + ---------- + sequence: str + A string of the protein sequence + + Returns + ------- + hits: dict + A dictionary generated from ProDy with PDB hits. + """ + + from prody import blastPDB + import tempfile + import pickle + + blast_record = blastPDB(sequence) + hits = blast_record.getHits() + + # TODO need to address when BLAST search times out + # TODO add option based on sequency similarity cut off + + return hits + + def get_sequence(self, identifier: str, kinase: bool = True, backend: str = "uniprot"): + """ + Retrieve a sequence based on an identifier code + + Parameters + ---------- + identifier: str + A string of the identifier to query (e.g 'P04629' from UniProt) + kinase: bool + Specify whether the sequence search should query kinase domains. + (Default: True) + backend: str + Specify which database to query. Options = ["uniprot", "ncbi"] + + Returns + ------- + sequence: str + A protein sequence. + """ + + import requests + from io import StringIO + + if kinase: + try: + from_method = getattr(KinaseDomainAminoAcidSequence, f"from_{backend}") + + except AttributeError: + # TODO implement ncbi backend + raise NotImplementedError( + f"Backend {backend} not supported. Please choose from ['uniprot', 'ncbi']" + ) + + else: + sequence = from_method(identifier) + + else: + params = {"query": identifier, "format": "fasta"} + response = requests.get("http://www.uniprot.org/uniprot/", params) + sequence = response.text.split("\n", 1)[1].replace("\n", "") + + return sequence + + def make_model( + self, + template_system: ProteinStructure, + target_sequence: KinaseDomainAminoAcidSequence, + alignment: Alignment, + num_models: int = 100, + ligand: bool = False, + ): + """ + Generate homology model(s) based on a template and alignment with MODELLER + + Parameters + ---------- + template_system: ProteinStructure + The template system. + target_sequence: KinaseDomainAminoAcidSequence + The target sequence + alignment: Alignment + An Alignment object containing information on aligned sequences. These + sequences must be the same as present in template_system and target_sequence + num_models: int + The number of homology models to generate. default = 100 + ligand: bool, optional=False + Specify whether to include a bound ligand present in the template structure. + """ + + from modeller import log, environ + from modeller.automodel import dope_loopmodel, assess + + log.verbose() + env = environ() + + pdb_id = template_system.metadata["id"] + + env.io.atom_files_directory = [ + template_system.metadata["path"].split(".")[0].split(pdb_id)[0] + ] + + if ligand: + # Read in HETATM records from template + env.io.hetatm = True + + # TODO handle usage of "ncbi" instead of "uniprot_id" + # specify model paramters and score via DOPE + a = dope_loopmodel( + env, + alnfile=alignment.alignment_file_path, + knowns=template_system.metadata["id"], + sequence=target_sequence.metadata["uniprot_id"], + loop_assess_methods=(assess.DOPE, assess.GA341), + ) + + a.starting_model = 1 + a.ending_model = num_models + + # TODO could add loop refinement option here + # TODO need to specify output directory for models + a.make()