{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from CalcLigRMSD import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# CalcLigRMSD" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "CalcLigRMSD calculates the Root-mean-square deviation (RMSD) between two prealigned molecules. \n", "\n", "This jupyter notebook shows examples of how to use the CalcLigRMSD function. \n", "This function is particularly useful in 2 cases:\n", "\n", "1. when the atom names in the compared structures do not match, i.e. when the same atoms present different names in the coordinate files. \n", "2. when one or both structures have missing atoms. This occurs, for example, in crystallographic structures when some atoms are not well defined in the electron density map.\n", "\n", "Below, there are also examples of how to align protein structures and extract ligand coordinates using a python/pymol function. The code below may be used both when the protein structures have multiple chains/domains and when protein structures present multiple copies of the ligand of interest, either in the binding pocket or on the surface." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### How to use the CalcLigRMSD function:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function CalcLigRMSD in module CalcLigRMSD:\n", "\n", "CalcLigRMSD(lig1, lig2, rename_lig2=True, output_filename='tmp.pdb')\n", " Calculate the Root-mean-square deviation (RMSD) between two prealigned ligands, \n", " even when atom names between the two ligands are not matching.\n", " The symmetry of the molecules is taken into consideration (e.g. tri-methyl groups). \n", " Moreover, if one ligand structure has missing atoms (e.g. undefined electron density in the crystal structure), \n", " the RMSD is calculated for the maximum common substructure (MCS).\n", " \n", " Parameters\n", " ----------\n", " lig1 : RDKit molecule\n", " lig2 : RDKit molecule\n", " rename_lig2 : bool, optional\n", " True to rename the atoms of lig2 according to the atom names of lig1\n", " output_filename : str, optional\n", " If rename_lig2 is set to True, a PDB file with the renamed lig2 atoms is written as output.\n", " This may be useful to check that the RMSD has been \"properly\" calculated, \n", " i.e. that the atoms have been properly matched for the calculation of the RMSD.\n", " \n", " Returns\n", " -------\n", " rmsd : float\n", " Root-mean-square deviation between the two input molecules\n", "\n" ] } ], "source": [ "help(CalcLigRMSD)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# RMSD between Pre-aligned Ligands with not Matching Atom Names " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic Example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Complex between the kinase protein AURKA and the JVE ligand\n", "\n", "Comparison between the crystal structure of JVE (pink, PDB code: 4UZH) and the docked pose of JVE (yellow).\n", "JVE was docked into the the AURKA structure 2C6E using AutodockVina. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 1.66\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(\"data/docked_2c6e_JVE_pH74_netcharge1.pdb\")\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/4uzh_JVE.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/JVE_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The atom names in the docked structure and in the crystal structure are different. However, the function is able to match the atoms and calculate the RMSD.\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Examples where symmetry needs to be considered" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Complex between the kinase protein AURKA and the N15 ligand" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 1.49\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(\"data/docked_2c6e_N15_pH74_netcharge1.pdb\")\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/3w2c_N15.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/N15_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The N15 ligand contains two symmetric groups (highlighted in the Figure below). The algorithm takes symmetry into consideration: it calculates the RMSD between all possible \"symmetrical\" combinations and returns the minimum RMSD. The structures with matched atom names corresponding to the lowest RMSD are reported in the figure below.\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Complex between the kinase protein AURKA and the AKI ligand" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 4.55\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(\"data/docked_2c6e_AKI_pH74_netcharge1.pdb\")\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/3m11_AKI.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/AKI_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The AKI ligand contains four symmetric groups (highlighted in the Figure below). Again, the code takes symmetry into consideration!\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example for a very symmetrical molecule\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "BEG is a very symmetrical molecule (see [this ref](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-019-0362-7)). Not only it does it contain two benzene groups but the entire molecule is symmetrical. We have downloaded the crystal structure 1D4I and extracted the coordinates of the BEG molecule (1d4i_BEG.pdb). In a copy (BEG_test.pdb), we have then renamed the atoms randomly. If the alghoritm is able to correctly match the atoms, then it will return an RMSD of zero. And indeed:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 0.00\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(\"data/BEG_test.pdb\")\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/1d4i_BEG.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = False, output_filename = 'data/AKI_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example where there are missing atoms in the crystal structure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Two crystal structures are available for the EGFR/1C9 complex. However, the electron density of the 1C9 is not well defined in any of the structures (see figure below). In such cases, one can calculate the RMSD between the docked and crystal structures by excluding the atoms which are not defined in the crystal structure.\n", "\n", "This is achived by calculating the RMSD for the maximum common substructure." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 6.47\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(\"data/docked_3w2p_1C9_pH74_netcharge1.pdb\")\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/4i23_1C9.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/1C9_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 7.71\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(\"data/docked_3w2p_1C9_pH74_netcharge1.pdb\")\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/4i24_1C9.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/1C9_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this figure, the missing groups are highlighted in green in the docked pose while the truncation point in the crystal structure is highlighted in red." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Align docking and crystal structures prior to the RMSD calculation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The function above takes as input pre-aligned ligand structures. In this section, we show how one could align crystallographic protein-ligand structure to the docked structure and extract the ligand coordinates from the crystal structure using Pymol. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "from pymol import cmd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Load the crystallographic and docked protein-ligand structures in pymol\n", "\n", "Docking programs often return as output only the coordinates of the docked ligand.\n", "Therefore, in order to be able to perform the protein alignment,\n", "we first need to load both the protein structure used for docking and the docked ligand structure in pymol\n", "and merge the coordinates in a temporary PDB file, complex.pdb\n", "\n", "\n", "Input files:\n", "- `docked_lig_file` : name of file containing the coordinates of the docked ligand\n", "- `protein_file` : name of file containing the the protein structure used for docking\n", "- `crystal_file`: name of file containing the crystal structure of the complex" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Function to merge the ligand and protein PDB files in a unique PDB file \n", "def merge_prot_lig_pdbs(protein_file, docked_lig_file, output_file = 'complex.pdb'):\n", " cmd.reinitialize()\n", " cmd.load(protein_file, 'protein')\n", " cmd.load(docked_lig_file, 'docked_lig')\n", " cmd.save(output_file, \"protein or docked_lig\")\n", " \n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " PyMOL not running, entering library mode (experimental)\n" ] } ], "source": [ "docked_lig_file = \"data/docked_2c6e_SKE_pH74_netcharge1.pdb\"\n", "protein_file = \"data/aurka_protein_2c6e.pdb\"\n", "crystal_file = \"data/5dt0.pdb\"\n", "\n", "merge_prot_lig_pdbs(protein_file, docked_lig_file)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can reinitialize Pymol and load both the docked and crystal complexes in Pymol. The step above can be skipped if the docked complex is directly available." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# load both the docked and crystallographic protein-ligand structures \n", "cmd.reinitialize()\n", "cmd.load(\"complex.pdb\", \"docked\")\n", "cmd.load(crystal_file, \"crystal\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Proteins Alignment" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simple Example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Align the proteins on the backbone alpha carbons and extract the coordinates of the ligand from the crystal structure:\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# align the protein structures on the CA \n", "cmd.align(\"crystal and name CA\", \"docked and name CA\")\n", "\n", "# save the coordinates of the ligand in a PDB file\n", "cmd.save(\"data/lig_crystal_aligned.pdb\", \"crystal and resn SKE\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, you can calculate the RMSD:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 8.21\n" ] } ], "source": [ "docked_ligand = Chem.MolFromPDBFile(docked_lig_file)\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/lig_crystal_aligned.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/1C9_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Alignment for Structures with Multiple Chains/Domains and Different/Multiple Ligands" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The structure alignment and ligand extraction becomes more complicated for crystal structures that contain dimers, or higher oligomers, multi-protein complexes, multiple and different ligands!\n", "In these case, the code below can be used. It reads as input the crystallographic and docked structures.\n", "The function does the following:\n", "\n", " - check for the presence of organic molecules in the crystal structure.\n", " - check for the presence of organic molecules in each of the chains\n", " - if a chain contains an organic molecule, then that chain is aligned to the docked structure based on the C-alpha backbone atoms.\n", " - calculate the distance between the organic molecule in the aligned crystal structure and the organic molecule in the aligned docked structure.\n", " - if the distance is lower than a given threshold (max_dist), the coordinates of the organic molecule from the aligned crystal structure are stored in a PDB file.\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "from pymol import stored\n", "import numpy as np \n", "\n", "def find_organic_ligs(obj_name):\n", " cmd.select(f\"ligs\", f\"{obj_name} and organic\")\n", " stored.residues = set()\n", " cmd.iterate(f\"ligs\", 'stored.residues.add(resn)')\n", " ligs_resname = list(stored.residues)\n", " return list(set(ligs_resname) - set(['ACE', 'NME']))\n", "\n", "\n", "def extract_lig_from_aligned_crystal(crystal_file, docked_file, max_dist = 7):\n", " count = 0\n", " # load structures\n", " cmd.reinitialize()\n", " cmd.load(docked_file, \"docked\")\n", " cmd.load(crystal_file, \"crystal\")\n", " # identify organic molecules in both the docked and crystal structures\n", " lig_docked_name = find_organic_ligs('docked')[0]\n", " lig_crystal_names = find_organic_ligs('crystal')\n", " # check the number of chains in the crystal structure\n", " chains = cmd.get_chains(\"crystal\")\n", " print(f\"The crystal structure {crystal_file} contains {len(chains)} chains and {len(lig_crystal_names)} ligands {lig_crystal_names}\")\n", " # loop over the organic molecules\n", " for lig_crystal_name in lig_crystal_names:\n", " print(f\"\\n{lig_crystal_name}\")\n", " # loop over the chains\n", " for chain in chains:\n", " sel = cmd.select(f\"ligs{chain}\", f\"crystal and chain {chain} and resn {lig_crystal_name}\")\n", " if sel == 0:\n", " continue\n", " else:\n", " cmd.align(f\"crystal and chain {chain} and name CA\", \"docked and name CA\")\n", " cmd.save(\"data/crystal_aligned.pdb\", \"crystal\")\n", " cmd.select(f\"ligs{chain}\", f\"crystal and chain {chain} and resn {lig_crystal_name}\")\n", " stored.residues = set()\n", " cmd.iterate(f\"ligs{chain}\", 'stored.residues.add(resv)')\n", " lig_res = list(stored.residues)\n", " # measure the distance between the center of mass of the ligands in docked and crystal structures.\n", " center2 = cmd.centerofmass(f\"docked and resn {lig_docked_name}\")\n", " # there may be multiple ligand copies in the crystal structure. Loop over them.\n", " for lig1 in lig_res:\n", " center1 = cmd.centerofmass(f\"crystal and chain {chain} and resi {lig1}\")\n", " # calculate distance\n", " dist = np.linalg.norm(np.array(center1)-np.array(center2)) \n", " # if the distance is lower than the threshold, the ligand is in the binding pocket.\n", " # Therefore, the coordinates of this crystal ligand are written out in a PDB file\n", " if dist < max_dist:\n", " count += 1\n", " cmd.save(f\"data/{lig_crystal_name}_crystal_aligned_{count}.pdb\", f\"crystal and chain {chain} and resi {lig1}\")\n", " print(f\"The coordinates of {lig_crystal_name} extracted from the crystal structure and chain {chain} are stored in data/{lig_crystal_name}_crystal_aligned_{count}.pdb\")\n", " else:\n", " print(f\"PDB {crystal_file} cmpd {lig_crystal_name} res {lig1} is out of the binding pocket\")\n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example where the crystal structure contains multiple chains\n", "\n", "The PDB structure 6C83 contains four protein chains, two of the target of interest AURKA and two of another protein.\n", "The two monomers of AURKA, both contain the compound of interest in the binding pocket. The other two proteins do not contain any organic molecule.\n", "\n", "In this case, we can use the function described above to perform the protein alignment and extract the coordinates of the ligand. As both monomers contain the ligand, two PDB files will be outputted.\n", "\n", "" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# merge protein and docked ligand\n", "docked_lig_file = \"data/docked_2c6e_ACP_pH74_netcharge1.pdb\"\n", "protein_file = \"data/aurka_protein_2c6e.pdb\"\n", "crystal_file = \"data/6c83.pdb\"\n", "\n", "merge_prot_lig_pdbs(protein_file, docked_lig_file)\n" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The crystal structure data/6c83.pdb contains 4 chains and 1 ligands ['ACP']\n", "\n", "ACP\n", "The coordinates of ACP extracted from the crystal structure and chain A are stored in data/ACP_crystal_aligned_1.pdb\n", "The coordinates of ACP extracted from the crystal structure and chain B are stored in data/ACP_crystal_aligned_2.pdb\n" ] } ], "source": [ "# extract the coordinated of the ligand from the crystal structure\n", "extract_lig_from_aligned_crystal(crystal_file, \"complex.pdb\", max_dist = 7)\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 3.82\n" ] } ], "source": [ "# Calculate RMSD for ACP extracted from chain A:\n", "docked_ligand = Chem.MolFromPDBFile(docked_lig_file)\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/ACP_crystal_aligned_1.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/ACP_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 3.02\n" ] } ], "source": [ "# Calculate RMSD for ACP extracted from chain B:\n", "docked_ligand = Chem.MolFromPDBFile(docked_lig_file)\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/ACP_crystal_aligned_2.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/ACP_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Example where the crystal structure contains different ligands\n", "\n", "The PDB structure 5DPV of AURKA contains two different ligands, 5DN and SKE. Tthe figure below shows that SKE (blue) is located in the binding pocket while 5DN (orange) lies on the surface of the protein. \n", "\n", "It occurs very often that a crystal structure contains different ligands or multiple copies of the same ligand. It is important to extract only the coordinates of the ligand located in the binding pocket, corresponding to the ligand of interest.\n", "\n", "The function described above checks the location of the ligand before extracting the coordinates. And this is achived by measuring the distance between the center of mass of the ligand in the docked structure (assuming that it is located in the binding pocket) and that of the ligand in the crystal structure. \n", "\n", "" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "docked_lig_file = \"data/docked_2c6e_SKE_pH74_netcharge1.pdb\"\n", "protein_file = \"data/aurka_protein_2c6e.pdb\"\n", "crystal_file = \"data/5dpv.pdb\"\n", "\n", "merge_prot_lig_pdbs(protein_file, docked_lig_file)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The crystal structure data/5dpv.pdb contains 1 chains and 2 ligands ['SKE', '5DN']\n", "\n", "SKE\n", "The coordinates of SKE extracted from the crystal structure and chain A are stored in data/SKE_crystal_aligned_1.pdb\n", "\n", "5DN\n", "PDB data/5dpv.pdb cmpd 5DN res 401 is out of the binding pocket\n" ] } ], "source": [ "extract_lig_from_aligned_crystal(crystal_file, \"complex.pdb\", max_dist = 7)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMSD: 8.17\n" ] } ], "source": [ "# Calculare RMSD:\n", "docked_ligand = Chem.MolFromPDBFile(docked_lig_file)\n", "crystal_ligand = Chem.MolFromPDBFile(\"data/SKE_crystal_aligned_1.pdb\")\n", "rmsd2 = CalcLigRMSD(docked_ligand, crystal_ligand, rename_lig2 = True, output_filename = 'data/SKE_renamed.pdb')\n", "print(f\"RMSD: {rmsd2:.2f}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.4" } }, "nbformat": 4, "nbformat_minor": 4 }