CCTBX_Demo.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "CCTBX Demo.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/cschlick/CryoEM-ReRefine-Ligands/blob/master/CCTBX_Demo.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jDU-VY8eo2cG"
      },
      "source": [
        "# Introduction to CCTBX\n",
        "\n",
        "The CCTBX project (computational crystallography toolbox) initially started as a library containing algorithms for the handling of unit cells, space groups, and atomic scatterers. Since then it has grown into a large project to support a variety of structural biology workflows, including x-ray and neutron crystallography, cryo-electron microscopy, and x-ray free electron lasers. \n",
        "\n",
        "The CCTBX library is written in Python, with many of the lower level algorithms accelerated using C++ via the Boost Python bindings.\n",
        "\n",
        "The principal high-level cctbx objects are:\n",
        "\n",
        "  - Data Manager: Read and write files and keep track of them\n",
        "  - Model Manager: A molecular model or models and restraints that go with them\n",
        "  - Map Manager: A 3D map or maps and information about gridding and map symmetry\n",
        "  - Map Model Manager: One or more related maps and models as a group\n",
        "  - F-Model Manager: Combines a model and a related Miller array (map coefficeints, or other data for each index of a map)\n",
        "\n",
        "\\\n",
        "These high-level cctbx objects are composed of lower-level objects that implement various concepts in macromolecular crystallography, which are in turn composed of CCTBX flex arrays for efficient computation of numerical data. Finally, user defined parameters are managed and shared betwen objects using the CCTBX Phil module. \n",
        "\n",
        "\\\n",
        "The examples presented here aim to illustrate working with the core CCTBX data types using example data. Much of the content is adapted from the official CCTBX documentation: http://cci.lbl.gov/docs/cctbx/\n",
        "\n",
        "\\\n",
        "Many of the examples are drawn directly from those created by Dorothee Liebschner and Tom Terwilliger. Additional credit goes to Billy Poon for packaging CCTBX as a conda package, and to Billy and Claudia for helping to get CCTBX working in a Colab notebook."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2zb2Szc1NfWN"
      },
      "source": [
        "#Installation\n",
        "### Installation is in three steps:\n",
        "####    1.  the Anaconda Python package manager, CCTBX, and py3Dmol\n",
        "####    2. Retrieve the chemical data repositories and some example data\n",
        "####    3. Restart the Runtime and load CCTBX environment variables"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7ZyxiGwSLcrs",
        "cellView": "form"
      },
      "source": [
        "#@title Install software\n",
        "#@markdown Please execute this cell by pressing the _Play_ button \n",
        "#@markdown on the left. Installation may take a few minutes.\n",
        "#@markdown Double click on this text to show/hide the installation script\n",
        "\n",
        "from IPython.utils import io\n",
        "import tqdm.notebook\n",
        "\n",
        "total = 50\n",
        "\n",
        "MINICONDA_PREFIX=\"/usr/local\"\n",
        "MINICONDA_INSTALLER_SCRIPT=\"Miniconda3-py37_4.9.2-Linux-x86_64.sh\"\n",
        "\n",
        "with tqdm.notebook.tqdm(total=total) as pbar:\n",
        "    with io.capture_output() as captured:\n",
        "      # install anaconda\n",
        "      %shell wget https://repo.continuum.io/miniconda/$MINICONDA_INSTALLER_SCRIPT\n",
        "      pbar.update(10)\n",
        "      %shell chmod +x {MINICONDA_INSTALLER_SCRIPT}\n",
        "      %shell ./{MINICONDA_INSTALLER_SCRIPT} -b -f -p {MINICONDA_PREFIX}\n",
        "      %shell conda install -c conda-forge cctbx -y\n",
        "      pbar.update(30)\n",
        "      \n",
        "      # provide location of libtbx environment\n",
        "      %shell ln -s /usr/local/share/cctbx /usr/share/cctbx\n",
        "\n",
        "      # replace system libstdc++ with conda version\n",
        "      # restart runtime after running\n",
        "      # to keep changes persistent, modify the filesystem, not environment variables\n",
        "      %shell sudo cp /usr/local/lib/libstdc++.so.6.0.29 /usr/lib/x86_64-linux-gnu/libstdc++.so.6 \n",
        "\n",
        "\n",
        "      # non cctbx installs\n",
        "      %shell pip install py3dmol\n",
        "      %shell apt install -y subversion\n",
        "\n",
        "      # python scripts for py3dmol + cctbx\n",
        "      %shell mkdir cctbx_py3dmol\n",
        "      %shell cd cctbx_py3dmol\n",
        "      %shell wget https://raw.githubusercontent.com/cschlick/cctbx-notebooks/master/cctbx_py3dmol_wrapper.py\n",
        "      %shell wget https://raw.githubusercontent.com/cschlick/cctbx-notebooks/master/cubetools.py\n",
        "      %shell cd ../\n",
        "      pbar.update(10)\n",
        "\n",
        "      "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SKOxhUpP0N5k",
        "cellView": "form"
      },
      "source": [
        "#@title Download data\n",
        "#@markdown Please execute this cell by pressing the _Play_ button \n",
        "#@markdown on the left. This should take a minute.\n",
        "%%capture\n",
        "%%bash\n",
        "\n",
        "# get chemical data\n",
        "cd  /usr/local/lib/python3.7/site-packages\n",
        "mkdir chem_data\n",
        "cd chem_data\n",
        "\n",
        "svn --quiet --non-interactive --trust-server-cert co https://github.com/phenix-project/geostd.git/trunk geostd\n",
        "svn --quiet --non-interactive --trust-server-cert co https://github.com/rlabduke/mon_lib.git/trunk mon_lib\n",
        "svn --quiet --non-interactive --trust-server-cert export https://github.com/rlabduke/reference_data.git/trunk/Top8000/Top8000_rotamer_pct_contour_grids rotarama_data\n",
        "svn --quiet --non-interactive --trust-server-cert --force export https://github.com/rlabduke/reference_data.git/trunk/Top8000/Top8000_ramachandran_pct_contour_grids rotarama_data\n",
        "svn --quiet --non-interactive --trust-server-cert co https://github.com/rlabduke/reference_data.git/trunk/Top8000/Top8000_cablam_pct_contour_grids cablam_data\n",
        "svn --quiet --non-interactive --trust-server-cert co https://github.com/rlabduke/reference_data.git/trunk/Top8000/rama_z rama_z\n",
        "\n",
        "# update rotamer and cablam cache\n",
        "/usr/local/bin/mmtbx.rebuild_rotarama_cache\n",
        "/usr/local/bin/mmtbx.rebuild_cablam_cache\n",
        "\n",
        "\n",
        "# get probe\n",
        "git clone https://github.com/cschlick/cctbx-notebooks.git\n",
        "mkdir -p /usr/local/share/cctbx/probe/exe\n",
        "cp cctbx-notebooks/probe /usr/local/share/cctbx/probe/exe\n",
        "chmod +x /usr/local/share/cctbx/probe/exe/probe\n",
        "mkdir -p /usr/lib/python3.7/site-packages/probe\n",
        "mkdir -p /usr/lib/python3.7/site-packages/reduce\n",
        "\n",
        "# tutorial data\n",
        "cd /content/ \n",
        "wget https://gitcdn.link/repo/phenix-lbl/cctbx_tutorial_files/master/2019_melk/1aba_pieces.pdb\n",
        "wget https://gitcdn.link/repo/phenix-lbl/cctbx_tutorial_files/master/2019_melk/1aba_model.pdb\n",
        "wget https://gitcdn.link/repo/phenix-lbl/cctbx_tutorial_files/master/2019_melk/4zyp.mtz\n",
        "wget https://gitcdn.link/repo/phenix-lbl/cctbx_tutorial_files/master/2019_melk/1aba_reflections.mtz\n",
        "wget https://gitcdn.link/repo/phenix-lbl/cctbx_tutorial_files/master/2019_melk/resname_mix.pdb\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2-vBTJbT-dt7"
      },
      "source": [
        "**WAIT!**...don't click the next button yet...Go to \"Runtime\" and select \"Restart Runtime\" first.  Then click the button below as noted."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GojrMInlLosF",
        "cellView": "form"
      },
      "source": [
        "#@title Set environment variables and do imports\n",
        "#@markdown Execute this cell by clicking the run button the left\n",
        "%%capture\n",
        "#  Manually update sys.path\n",
        "%env PYTHONPATH=\"$/env/python:/usr/local/lib:/usr/local/lib/python3.7/lib-dynload:/usr/local/lib/python3.7/site-packages\"\n",
        "\n",
        "# add conda paths to sys.path\n",
        "import sys, os\n",
        "sys.path.extend(['/usr/local/lib', '/usr/local/lib/python3.7/lib-dynload', '/usr/local/lib/python3.7/site-packages',\"/content/cctbx_py3dmol\"])\n",
        "os.environ[\"MMTBX_CCP4_MONOMER_LIB\"]=\"/usr/local/lib/python3.7/site-packages/chem_data/geostd\"\n",
        "os.environ[\"LIBTBX_BUILD\"]= \"/usr/local/share/cctbx\"\n",
        "\n",
        "import cctbx\n",
        "import libtbx\n",
        "libtbx.env.add_repository(\"/usr/local/lib/python3.7/site-packages/\")\n",
        "libtbx.env.module_dist_paths['probe']=\"\"\n",
        "\n",
        "from cubetools import *\n",
        "from cctbx_py3dmol_wrapper import CCTBX3dmolWrapper\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WTtUdcU_NOHi"
      },
      "source": [
        "#### Test\n",
        "We should (hopefully) not see any errors. Use the Run button or Shift+Enter to execute the cells"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "to8tosuwM7t4"
      },
      "source": [
        "import cctbx\n",
        "from scitbx.array_family import flex\n",
        "a = flex.random_double(100)\n",
        "print(flex.min(a))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8QOMLb5LiH1u"
      },
      "source": [
        "We can also run linux shell commands from Colab cells, either using the %shell syntax (using during installation), or the ! prefix for a single line. Lets check that the downloads were successful and that we have some files to work with. For example, we hope to see 1aba_model.pdb in the current working directory."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tgUQut9fiEX6"
      },
      "source": [
        "! pwd"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Vk_xLQGU5kU5"
      },
      "source": [
        "! ls"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lx9V36w3NnGF"
      },
      "source": [
        "# File IO via the DataManager\n",
        "\n",
        "When you want to perform some actions on data or a model, you have to start with reading the file. Once the action is done, you probably want to write the files. \n",
        "\n",
        "To do this conveniently, you can use the cctbx DataManager. The DataManager lets you read and write files describing atomic models, restraints, reflection data files, symmetry files, sequence files, and MRC/CCP4 format map files."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6NVbHFPRwoky"
      },
      "source": [
        "from iotbx.data_manager import DataManager\n",
        "\n",
        "dm = DataManager()\n",
        "filename_pdb = \"1aba_model.pdb\"\n",
        "dm.process_model_file(filename_pdb)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZMsYsx85oWzK"
      },
      "source": [
        "The data manager has now read our file. If desired, a single data manager instance can read many files and store them simultaneously. To access a specific file, we can request it using the filename as the key..."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8w6APmZpolDO"
      },
      "source": [
        "model = dm.get_model(filename=filename_pdb)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AjeWyuCZoml9"
      },
      "source": [
        "The model object is an instance of the CCTBX Model Manager. The Model Manager stores all information related to a molecular model. This includes the molecular composition, space group and symmetry information, information about geometry restraints, and any metadata that was available in the input file.\n",
        "\n",
        "\\\n",
        "We will use the py3dmol package: https://3dmol.org/ to display structures in the notebook. We can customize the depictions by writing a custom show function (already downloaded during the data download step). Run the next cell to import the show function and show your model. You can manipulate the model with your mouse..."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8wk23bK7PMN-"
      },
      "source": [
        "from cctbx_py3dmol_wrapper import CCTBX3dmolWrapper\n",
        "show = CCTBX3dmolWrapper()\n",
        "show(model)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xrpcArdSwQp8"
      },
      "source": [
        "To save a model to the disk, we also use the CCTBX data manager. One of the many advantages of the data manager is the ability to handle multiple file formats transparently."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EHnHRVwdgKm9"
      },
      "source": [
        "filename_cif = \"1aba_model.cif\"\n",
        "dm.write_model_file(model,filename=filename_cif,format=\"cif\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FmudB0Yskdyp"
      },
      "source": [
        "Lets check to make sure it worked"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qDF-8kuATyy6"
      },
      "source": [
        "! head 1aba_model.cif -n 100"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nSbof8NlT28v"
      },
      "source": [
        "# Models"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fW4SpgeExQU5"
      },
      "source": [
        "## Model Input\n",
        "\n",
        "The model manager object provides a a container for all information related to a model. The categories of data inclue the model metadata, the model composition, and geometry geometry restraints for the model. By maintaining a clear separation between these categories of data, we can avoid passing unneccessary information to low level functions.\n",
        "\n",
        "\\\n",
        "First we will consider the model metadata which is stored as a model input object. In contrast to the model composition which is highly structured data, the model input can vary substantially depending on the source of the file and what the authors chose to include. Lets take a look at the model input for this file. Execute the cell below to get the model_input object, then the button below that to print out some characteristics.|"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RlkLxSJ8sdGa"
      },
      "source": [
        "model_input= model.get_model_input()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O7DBf-GEw4Ib"
      },
      "source": [
        "Some fields are structured, and are directly accessible by member functions"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KgDim35XsdJC"
      },
      "source": [
        "print(model_input.file_type())\n",
        "print(model_input.get_experiment_type())\n",
        "print(model_input.get_program_name())\n",
        "print(model_input.deposition_date())\n",
        "print(model_input.resolution())\n",
        "print(model_input.sequence_from_SEQRES())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iLZy5mQEw_Oc"
      },
      "source": [
        "The entire contents of the file are retained, and can be parsed to extract miscellaneous fields. Note that many fields relating to the data are not populated. However, we can determine them later when we look at the data for this structure."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gE5wmcwxudGw"
      },
      "source": [
        "# Get the whole title section of the pdb\n",
        "print(\"\\n\".join(model_input.title_section()))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JpgyqW_wKpat"
      },
      "source": [
        "# print the remark section\n",
        "print(\"\\n\".join(model_input.remark_section())[:2000])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n76Hi28SuKvE"
      },
      "source": [
        "## Crystal symmetry\n",
        "\n",
        "The crystal symmetry object stores structured data about the unit cell and space group"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dbWcoVaCf7L7"
      },
      "source": [
        "cs = model.crystal_symmetry()\n",
        "print(cs)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PQLrv9p3xU6k"
      },
      "source": [
        "## Model composition\n",
        "The core componnet of the Model Manager object is the information about molecular composition, stored as a hierarchy object. The hierarchy is written as a C++ object for speed, and structured to reflect the hierarchical nature of macromolecules. It also aims to be independent of the original file format. That is, the hierarchy object doesn't know or care if the original file was a PDB file, mmCIF file, or any other molecular file format.\n",
        "\n",
        "\\\n",
        "Lets take a closer look at the composition of this protein"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SFwvomYt0fBh"
      },
      "source": [
        "hierarchy = model.get_hierarchy()\n",
        "print(hierarchy.composition())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GLUFPuOu2I3q"
      },
      "source": [
        "#### Group args\n",
        "The composition method returns a group_args object. This is a CCTBX object that is commonly returned as the result from CCTBX programs and methods. It essentially functions the same as a dictionary, holding data as a set of key,value pairs. The advantage of the group_args object is that the keys can be accessed using dot syntax, ie: \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DmfVYbfFx-4d"
      },
      "source": [
        "result = hierarchy.composition()\n",
        "n_atoms = result.n_atoms\n",
        "print(n_atoms)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QbmfDcJeyIjz"
      },
      "source": [
        "This is a convienent way to work with program results, especially if the results contain multiple nested group_args. Alternatively, the fields can be accessed by iterating through the keys"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kjiEC6XDyFG8"
      },
      "source": [
        "for key in result.keys():\n",
        "  print(key,\":\",result.get(key))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZjFwXZkL5Z4V"
      },
      "source": [
        "### Architecture of the hierarchy object\n",
        "\n",
        "```\n",
        "model(s)\n",
        "  id\n",
        "  chain(s)\n",
        "    id\n",
        "    residue_group(s)\n",
        "      resseq\n",
        "      icode\n",
        "      atom_group(s)\n",
        "        altloc,\n",
        "        resname\n",
        "          atom(s)\n",
        "            name\n",
        "            element\n",
        "            charge\n",
        "            xyz\n",
        "            occ\n",
        "            b\n",
        "```\n",
        " - The model, chain, and atom levels of the hierarchy object are probably immediately obvious to someone familiar with the content of model files (such as a PDB file). Note that that is no 'residue type' in the data structure. Instead, there are the two types residue_group and atom_group. They are related to alternative conformations. \n",
        "\n",
        " - If there are no alternative conformations in the model, all residue groups contain exactly one atom group, which contains all the atoms of a residue. A file with alternative conformations will lead to residue groups with multiple atom groups, one for each conformer. (Note: about a quarter of the files in the PDB database contain alternative conformations).\n",
        "\n",
        "\n",
        "#### Looping over hierarchy"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KM6CDHow8HnY"
      },
      "source": [
        "for m in hierarchy.models(): \n",
        "  for chain in m.chains():\n",
        "    for rg in chain.residue_groups():\n",
        "      for ag in rg.atom_groups():\n",
        "        for atom in ag.atoms():\n",
        "          # your code here\n",
        "          pass\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CnPQpHD8Koq_"
      },
      "source": [
        "Many hierarchy objects have methods to directly access objects lower down in the hierarchy as an iterable. Lets see what is in the first four residues (atom groups)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2ofZ8oTzKgfp"
      },
      "source": [
        "atom_groups = list(hierarchy.atom_groups())\n",
        "\n",
        "for ag in atom_groups[:4]:\n",
        "  rg = ag.parent() # access the residue group\n",
        "  chain = rg.parent().parent() # access the chain\n",
        "  print(\"Chain ID:\",chain.id)\n",
        "  print(\" Residue Num:\",rg.resseq)\n",
        "  print(\" Residue Name:\",ag.resname)\n",
        "  for atom in ag.atoms():\n",
        "    print(\"     Atom:\",atom.name)\n",
        "  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sVCslaZtMpRQ"
      },
      "source": [
        "###Example: Truncate to Poly-alanine\n",
        "We have residue names and atom names now, but to perform the truncation, we need to know which residues are amino acids (because we only want to truncate amino acids), and which atom names we want to keep.\n",
        "The iotbx.pdb module contains the sub-module amino_acid_codes. This sub-module contains two Python dictionaries, one of which is (shortened):\n",
        "```\n",
        "one_letter_given_three_letter = {\n",
        "\"ALA\": \"A\",\n",
        "\"ARG\": \"R\",\n",
        "...\n",
        "\"TYR\": \"Y\",\n",
        "\"VAL\": \"V\"}\n",
        "```\n",
        "\n",
        "We don't need the one-letter codes, but we can re-use the keys of this dictionary to efficently decide if a residue name corresponds to an amino acid."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kc9d4CcnKgsa"
      },
      "source": [
        "import iotbx.pdb.amino_acid_codes\n",
        "aa_resnames = iotbx.pdb.amino_acid_codes.one_letter_given_three_letter\n",
        "ala_atom_names = set([\" N  \", \" CA \", \" C  \", \" O  \", \" CB \"])\n",
        "\n",
        "\n",
        "def truncate_polyala(model):\n",
        "  model_polyala = model.deep_copy()\n",
        "\n",
        "  for m in model_polyala.get_hierarchy().models():\n",
        "    for chain in m.chains():\n",
        "      for rg in chain.residue_groups():\n",
        "        for ag in rg.atom_groups():\n",
        "          if ag.resname in aa_resnames:\n",
        "            ag.resname = \"ALA\"\n",
        "            for atom in ag.atoms():\n",
        "              if atom.name not in ala_atom_names:\n",
        "                ag.remove_atom(atom=atom)\n",
        "      \n",
        "    return model_polyala\n",
        "\n",
        "model_polyala = truncate_polyala(model)\n",
        "show(model_polyala)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5Sb0jXht6wZD"
      },
      "source": [
        "#### Improved truncation for rare cases\n",
        "For most practical purposes, the previous the previous method is sufficient. However, there are some files in the PDB for which this is not true. One example is the structure with the PDB ID 1ysl. Lets see what happens when we run the previous method on a fragment from 1ysl:\n",
        "\n",
        "**Note:** We are using the same instance of Data Manager here that we used before. The datamanager can store multiple files, accessed by their filename"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N9E4fvC7O4QJ"
      },
      "source": [
        "filename_1ysl = \"resname_mix.pdb\"\n",
        "dm.process_model_file(filename=filename_1ysl)\n",
        "model_1ysl = dm.get_model(filename=filename_1ysl)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UXk5HzYn8W2U"
      },
      "source": [
        "show(model_1ysl)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "traI-wAZ-RFI"
      },
      "source": [
        "Lets try to run our poly-alanine truncation function on this fragment:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BMHyvEpu8wSQ"
      },
      "source": [
        "model_1ysl_trunc= truncate_polyala(model_1ysl)\n",
        "show(model_1ysl_trunc)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xL9xdn258-zs"
      },
      "source": [
        "It is not all alanine. The reason is because this fragment contains a modified residue CSD (oxidized cysteine). In fact, this residue also has an alternate conformation. Conformation \"A\" is the modified amino acid CSD, conformation \"B\" is a regular CYS."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "D_LrhCG0-L0w"
      },
      "source": [
        "print(model_1ysl_trunc.model_as_pdb())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E-B9rPWR-TCS"
      },
      "source": [
        "Rare cases like this are the reason why the residue_group and atom_group levels are needed in the hierarchy. Here are two different residue names for the same member of a chain. Even though this sitution is rare, it is entirely plausible and valid, and a comprehensive PDB processing library has to be able to handle it.\n",
        "Our script will only truncate the CYS residue, but it would be better if it also truncated the non-standard CSD residue in the A alternative conformation. Let's find out what it takes to achieve this.\n",
        "\n",
        "One way would be to add the 3 letter code for all modified amino acids to the dictionary. This however relies on the completeness of the dictionary; what happens if the 3 letter code changes or if other entries are added?\n",
        "Another possibility is to check if there is at least one amino acid in a residue group, and if so, apply the truncation to all residues in the group, even if they don't have a standard residue name. This means, for each residue group we have to loop over the atom groups twice, first to scan for at least one standard amino-acid residue name, and if there is one, a second time to do the truncation. This will double effort (and computing time), but it is required to accomodate all possible cases."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UDMOwpfkCMo_"
      },
      "source": [
        "import iotbx.pdb.amino_acid_codes\n",
        "aa_resnames = iotbx.pdb.amino_acid_codes.one_letter_given_three_letter\n",
        "ala_atom_names = set([\" N  \", \" CA \", \" C  \", \" O  \", \" CB \"])\n",
        "\n",
        "\n",
        "def truncate_polyala_improved(model):\n",
        "  model_polyala = model.deep_copy()\n",
        "\n",
        "  def has_amino_acid(residue_group):\n",
        "    for ag in residue_group.atom_groups():\n",
        "      if ag.resname in aa_resnames:\n",
        "        return True\n",
        "    return False\n",
        "\n",
        "\n",
        "  for m in model_polyala.get_hierarchy().models():\n",
        "    for chain in m.chains():\n",
        "      for rg in chain.residue_groups():\n",
        "        if has_amino_acid(rg):\n",
        "          for ag in rg.atom_groups():\n",
        "              ag.resname = \"ALA\"\n",
        "              for atom in ag.atoms():\n",
        "                if atom.name not in ala_atom_names:\n",
        "                  ag.remove_atom(atom=atom)\n",
        "      \n",
        "    return model_polyala\n",
        "\n",
        "model_polyala = truncate_polyala(model)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OD4jcuNeDZie"
      },
      "source": [
        "model_1ysl_trunc = truncate_polyala_improved(model_1ysl)\n",
        "show(model_1ysl_trunc)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nSvNwE8WCERH"
      },
      "source": [
        "Success! The fragment was truncated to poly-alanine. Hopefully this section emphasizes the distinction between residue groups and atom groups in the CCTBX hierarchy."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hxdooC4AEjob"
      },
      "source": [
        "### Selections\n",
        "In many cases it is desireable to process a subset of a model. Fortunately, the model object supports a variety of string selections. In many cases this can enable you to skip looping over the hierarchy. \n",
        "\n",
        "\\\n",
        "For example, to select the first 20 residues:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1_Vo7i_eLDMv"
      },
      "source": [
        "model_subset = model.select(model.selection(\"resseq 1:20\"))\n",
        "show(model_subset)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UuBsic9eLaRZ"
      },
      "source": [
        "What is happening here? The inner statement parses the string selection, and turns it into a boolean flex array. Flex arrays are CCTBX numerical arrays implemented in C++. Each value in the array relates to an atom in the model, and the value True/False corresponds to whether or not the atom was selected by the selection string."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3FFvTCq0Llqz"
      },
      "source": [
        "sel = model.selection(\"resseq 1:20\")\n",
        "print(sel)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rEXdaypcNeeA"
      },
      "source": [
        "We can see how many values are set to True using the count() method on the flex array. We can see the total length of the array using the normal Python len() function."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jmua7ExILp-g"
      },
      "source": [
        "print(sel.count(True))\n",
        "print(len(sel))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eIExWt6yNjl0"
      },
      "source": [
        "This corresponds to 160 atoms selected, out of 880 atoms total. You may recall from the hierarchy.composition() function, there are 880 atoms in this structure.\n",
        "\n",
        "\\\n",
        "The outer component of our selection statement above, the select() function, takes the boolean flex array and converts it into a new model that only includes the selected atoms. This means that you can construct your own selection array if desired instead of using the string selection.\n",
        "\n",
        "\\\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VhtvDl00QYTx"
      },
      "source": [
        "Other options for string selections are documented on this Phenix site:\\\n",
        " https://phenix-online.org/documentation/reference/atom_selections.html\n",
        "\n",
        " \\\n",
        " For example, say we want to select all aromatic amino acids. This will demonstrate using the logical \"or\" expression in string selections:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2Mfodf5QOo7E"
      },
      "source": [
        "model_aromatic = model.select(model.selection(\"resname TYR or resname PHE or resname TRP\"))\n",
        "show(model_aromatic)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6-mWbzpIRUqi"
      },
      "source": [
        "There were not any Tryptophans, but we efficiently selected all the Tyrosines and Phenylalanines. We could have done this by looping over the hierarchy, but doing that manually is significantly slower"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fURf9mj7TGhu"
      },
      "source": [
        "%%time\n",
        "sel = [False]*model.get_number_of_atoms()\n",
        "\n",
        "for m in model.get_hierarchy().models():\n",
        "  for chain in m.chains():\n",
        "    for rg in chain.residue_groups():\n",
        "      for ag in rg.atom_groups():\n",
        "        if ag.resname in [\"PHE\",\"TYR\",\"TRP\"]:\n",
        "          for atom in rg.atoms():\n",
        "            sel[atom.i_seq]=True\n",
        "\n",
        "sel = flex.bool(sel)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mtmdr0sFS7cG"
      },
      "source": [
        "%%time \n",
        "sel = model.selection(\"resname TYR or resname PHE or resname TRP\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V6QVTcs9Tdeb"
      },
      "source": [
        "Note that the first example, looping over the hierarchy, uses the atom.i_seq property. This is the unique integer sequence value for each atom, and this is the position in the selection array that the atom corresponds to."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v8EmKUfaTXBT"
      },
      "source": [
        "#### Flex and Numpy: An example using B-factors and coordinates\n",
        "\n",
        "Flex arrays are numerical arrays implemented in C++. They have a lot in common with the widely used Numpy arrays. In fact, we can directly access a Numpy array from a Flex array. This is useful for passing data on to other libraries based on Numpy."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MznfYcQBT4xx"
      },
      "source": [
        "b_iso = model.get_b_iso()\n",
        "print(b_iso)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7nkbkSP7VGAF"
      },
      "source": [
        "b_iso_np = b_iso.as_numpy_array()\n",
        "print(type(b_iso_np))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5mV9CoLWE_UI"
      },
      "source": [
        "b_iso_np"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0NsXm0RNVJxc"
      },
      "source": [
        "list(b_iso_np)==list(b_iso)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dg0Wl7QiWBrP"
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "np.all(np.array(b_iso) == b_iso_np)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "keQq8EmzWIGa"
      },
      "source": [
        "all(flex.double(b_iso_np)==b_iso)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DEHmzmskVR44"
      },
      "source": [
        "Both can easily be converted back and forth directly, or through an intermediate Python list. Because of this interoperabiliy, we can easily use any of the Python libraries based on Numpy, for example Matplotlib."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7zbLuoWOT9Ng"
      },
      "source": [
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OHLAYXHpU3_u"
      },
      "source": [
        "_ = plt.hist(b_iso)\n",
        "plt.title(\"B-factor histogram\")\n",
        "plt.xlabel(\"B-factor\")\n",
        "plt.ylabel(\"Count\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IdnzsD_EYKF3"
      },
      "source": [
        "##### Cartesian coordinates\n",
        "Lets extract the cartesian coordinates for the model. It will be a flex array of N 3D vectors, where N is the number of atoms."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QaNKpRoojwf1"
      },
      "source": [
        "model_copy = model.deep_copy()\n",
        "xyz = model_copy.get_sites_cart()\n",
        "xyz"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iEXBlyKOlD45"
      },
      "source": [
        "We will move a fragment of 100 atoms by a translation vector to illustrate setting new coordinates for a model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8XublZH_j0W-"
      },
      "source": [
        "translation = (20,20,20)\n",
        "\n",
        "for i in range(100):\n",
        "  x,y,z = xyz[i]\n",
        "  a,b,c = translation\n",
        "  xyz[i] = x+a, y+b, z+c"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hS0t1EankjsQ"
      },
      "source": [
        "model_copy.set_sites_cart(xyz)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ahWee-QJknku"
      },
      "source": [
        "show(model_copy)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ae61pOYWMvbz"
      },
      "source": [
        "## Model geometry\n",
        "The final major category of information stored by the model manager is information related to geometry restraints. Geometry restraints provide a framework for evaluating how far away a model is from an idealized model.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u7Xi-_JsUP6R"
      },
      "source": [
        "model.process_input_model(make_restraints=True)\n",
        "model.geometry_statistics().show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "njoeCt55UuX3"
      },
      "source": [
        "geometry_restraints = model.get_restraints_manager()\n",
        "geometry_energies = geometry_restraints.geometry.energies_sites(sites_cart = model.get_sites_cart())\n",
        "angle_deltas = geometry_energies.angle_proxies.deltas(sites_cart=model.get_sites_cart())\n",
        "bond_deltas = geometry_energies.bond_proxies.deltas(sites_cart=model.get_sites_cart())\n",
        "\n",
        "fig,axes = plt.subplots(1,2,figsize=(12,5))\n",
        "ax1,ax2 = axes\n",
        "ax1.hist(bond_deltas,bins=50)\n",
        "ax1.set_xlabel(\"Bond delta (angstroms)\")\n",
        "\n",
        "ax2.hist(angle_deltas,bins=50)\n",
        "ax2.set_xlabel(\"Angle delta (degrees)\")\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EbA64ye8nVTD"
      },
      "source": [
        "Knowing how far a feature deviates from the ideal provides a way to incorporate molecular geometry into a target function. A large number of restraint types can be managed using the geometry restraints manager:\n",
        "  - Bonds\n",
        "  - angles\n",
        "  - chirality\n",
        "  - dihedrals\n",
        "  - planarity\n",
        "  - secondary structure\n",
        "  - non-bonded overlaps\n",
        "  - ramachandran\n",
        "  - rotamers\n",
        "  - etc"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fERaKyGztXxW"
      },
      "source": [
        "# Maps"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YKblpLp-C8Sd"
      },
      "source": [
        "##The Map Model Manager\n",
        "After models, probably the second most integral component of the CCTBX data structures are density maps. The two data structures we will introduce here are the Map manager and the Map Model Manager. The Map Model Manager in particular is a very high level object with many convienence functions. For example, if we want to generate some test data without reading in files, we can do that."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T6xqgx1BC9GJ"
      },
      "source": [
        "from iotbx.map_model_manager import map_model_manager      \n",
        "mmm=map_model_manager()\n",
        "mmm.generate_map()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sCU0VwQEC-Gs"
      },
      "source": [
        "show(mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lpKPZagKhpQ7"
      },
      "source": [
        "This small fragment was taken from a library, and the map was calculated using the model. The generate_map() method has many options that let you choose parameters such as:\n",
        "  - The number of residues (n_residues)\n",
        "  - The B-factor (b_iso)\n",
        "  - The high resolution limit (d_min)\n",
        "  - The scattering table (scattering_table)\n",
        "  \n",
        "It is also possible to supply a model or model file directly. Note that the default choice for the scattering table is “electron”, so it will create a cryo-EM map. If you want an electron density or nuclear density map, the value will have to be modified accordingly. \n",
        "\n",
        "\\\n",
        "Let's experiment with the effect of changing the B-factor and resolution for the simulated map."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x8-RClxGhnTJ"
      },
      "source": [
        "d_min = 8\n",
        "mmm=map_model_manager()\n",
        "mmm.set_resolution(d_min)\n",
        "mmm.generate_map(d_min=d_min)\n",
        "show(mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T2njRWffnYUO"
      },
      "source": [
        "This will be very familiar to all cryo-electron microscopists. Lets try to get a very nice looking map."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "usXuGMIIlZ8r"
      },
      "source": [
        "d_min = 1.0\n",
        "mmm=map_model_manager()\n",
        "mmm.set_resolution(d_min)\n",
        "mmm.generate_map(d_min=d_min,b_iso=10)\n",
        "show(mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JXsK7iagm0mg"
      },
      "source": [
        "We can play with different values for the isotropic B-factor..."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QNAK8PLVl9mB"
      },
      "source": [
        "d_min = 1.0\n",
        "mmm=map_model_manager()\n",
        "mmm.set_resolution(d_min)\n",
        "mmm.generate_map(d_min=d_min,b_iso=80)\n",
        "show(mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kXDUe6al28S_"
      },
      "source": [
        "When calculating the map using a higher isotropic B-factor, the density is visually degraded, which is what we expect."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LYlrRfB8pJsl"
      },
      "source": [
        "Lets generate a map for the 1aba model we have been working on. We will aim to generate a map at the resolution typical for cryo-EM data. Because we are trying to simulate a cryo-EM map, we need to change the space group to P1\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DTUYkEuE89_o"
      },
      "source": [
        "# 1aba crystal symmetry (Not P1)\n",
        "print(\"\\n1aba:\\n\",model.crystal_symmetry())\n",
        "\n",
        "from cctbx.maptbx.box import shift_and_box_model\n",
        "model_p1 = shift_and_box_model(model)\n",
        "\n",
        "# should be P1\n",
        "print(\"\\n1aba reboxed as P1:\\n\",model_p1.crystal_symmetry())\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QFxfzPTJJcBV"
      },
      "source": [
        "You can see that the shift_and_box_model() function made a P1 Unit cell that is big enough to accommodate the model. This is commonly used as the crystal symmetry for cryo-EM maps in CCTBX. Now we will simulate the density."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4fDf-tigqRCL"
      },
      "source": [
        "mmm=map_model_manager()\n",
        "mmm.generate_map(model=model_p1,d_min=3)\n",
        "show(mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0ivVNkOuepeZ"
      },
      "source": [
        "##The Map Manager"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JzlRrWIiKCc3"
      },
      "source": [
        "To access the raw numerical density values, we can access the map_data flex on the Map Manager object"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "MWBu5wqZMN-m"
      },
      "source": [
        "mm = mmm.map_manager()\n",
        "map_data = mm.map_data()\n",
        "map_data"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0Sj5CqMoMUNt"
      },
      "source": [
        "The shape of the map_data array can be accessed using the all() method"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uxY7S9xZMw6c"
      },
      "source": [
        "map_data.all()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kg87PjSwMzK9"
      },
      "source": [
        "These are the dimensions of the volume data. They were determined automatically when we generated the map, based on the resolution supplied and the unit cell dimensions. We can retrieve the resulting pixel sizes from the Map Manager as well:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qYfoWWxJMyep"
      },
      "source": [
        "mm.pixel_sizes()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PQ4U_xhDO_0v"
      },
      "source": [
        "### Extracting a sub map\n",
        "\n",
        "Often times it is desireable to only work with the density in one local area, for example, around a ligand. One of the convienent features of the Map Model Manager is that it can do this for you. Lets try to look at the density for the MES ligand."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sj8qEXsTQf3_"
      },
      "source": [
        "ligand_mmm = mmm.extract_all_maps_around_model(selection_string=\"resname MES\")\n",
        "show(ligand_mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4v-ersL1KldQ"
      },
      "source": [
        "The box_cushion argument controls how large the extracted density is. We can try with a smaller value if we don't need so much density in the neighborhood of our ligand"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rWfniX_hTxXx"
      },
      "source": [
        "ligand_mmm = mmm.extract_all_maps_around_model(selection_string=\"resname MES\",\n",
        "                                               box_cushion=2.0)\n",
        "show(ligand_mmm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jwmSORCiViJq"
      },
      "source": [
        "Now we can more easily evaluate the local density around the ligand. To adjust the contour level, we can pass different values to the show() function. First lets histogram the density values:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ByIOa8-eTV9F"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "_ = plt.hist(ligand_mmm.map_data(),bins=30)\n",
        "plt.xlabel(\"Map value\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D9T78bhUWheK"
      },
      "source": [
        "Lets set the value to 0.4"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "okngXvDTTb9P"
      },
      "source": [
        "show(ligand_mmm,map_threshold=0.4)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0yyF4PECMF5B"
      },
      "source": [
        "To finish, we will now save the map to illustrate how to save real space maps using the Data Manager"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AgxA1COPKBkr"
      },
      "source": [
        "dm.write_real_map_file(mm,filename=\"1aba_calculated_at_3A.mrc\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h0GitqEdtEUD"
      },
      "source": [
        "# Reflections"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uIUnEcvPC-Xj"
      },
      "source": [
        "\n",
        "## Miller array\n",
        "Miller arrays are containers for experimental (or calculated data). A miller array contains the crystal symmetry, an array of Miller indices (h,k,l), a boolean flag indicating anomalous pairs and a flex array containing data (X-ray amplitudes or intensities, experimental sigmas, etc.).\n",
        "\n",
        "\\\n",
        "Lets open the reflections for a new structure, 4zyp\n",
        " "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xyzZfy2zDKzL"
      },
      "source": [
        "filename = \"/content/4zyp.mtz\"\n",
        "dm.process_miller_array_file(filename)\n",
        "reflections_reader = dm.get_miller_array(filename=filename)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D30HUjGhXU3f"
      },
      "source": [
        "The reflections reader is a container which can store multiple miller arrays. Lets build a dictionary of all the miller arrays and the data labels there were present in the file"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EoqGS5kDD2Ps"
      },
      "source": [
        "miller_array_dict = {ma.info().label_string():ma for ma in reflections_reader.as_miller_arrays()}\n",
        "\n",
        "for key,value in miller_array_dict.items():\n",
        "  print(key,\":\",type(value))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cigzUnzrXyon"
      },
      "source": [
        "We will choose the FOBS,SIGFOBS array and look at some of the data stored in the miller array"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3jnR529jEYiA"
      },
      "source": [
        "f_obs = miller_array_dict[\"FOBS,SIGFOBS\"]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-iJIKrx2D0g2"
      },
      "source": [
        "print(\"Resolution limits:\", f_obs.d_max_min())\n",
        "print(\"\\nSpace group info:\",f_obs.space_group_info())\n",
        "print(\"\\nCompleteness:\",f_obs.completeness())\n",
        "print(\"\\nSize:\",f_obs.size())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3U3OG7YeYJul"
      },
      "source": [
        "We can also select the R-free in the MTZ file. Let's have a closer look at that array."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aeQPQIhYRQUC"
      },
      "source": [
        "r_free = miller_array_dict[\"R-free-flags\"]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OMRO5Z72YWWd"
      },
      "source": [
        "print(\"Resolution limits:\", r_free.d_max_min())\n",
        "print(\"\\nSpace group info:\",r_free.space_group_info())\n",
        "print(\"\\nCompleteness:\",r_free.completeness())\n",
        "print(\"\\nSize:\",r_free.size())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_q9_WmMsYg7o"
      },
      "source": [
        "The size of the miller arrays FOBS,SIGFOBS and R-free-flags is not the same. In other words, the number of hkl indices with corresponding data or R-free-flag differs in the two arrays. This situation occurs rather frequently, even with deposited data in the PDB. For example, the mismatch can occur when R-free-flags are copied from another data set.\n",
        "\n",
        "\\\n",
        "Programs dealing with data have to address this mismatch. Let's say we want to calculate R-factors. If a reflection with an associated intensity does not have an R-free flag, it can't be assigned neither to the working set nor the the test set of reflections. It is therefore useful to obtain a common set of reflections. Let's try to obtain such a common set for out test data set.\n",
        "\n",
        "\\\n",
        "It is instructive to determine how many reflections are mismatched. The miller class has the method match_indices, which will do all the work."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IuzbHBCjSLH3"
      },
      "source": [
        "match_result = f_obs.match_indices(r_free)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_VABzUIhSEyT"
      },
      "source": [
        "n_common = match_result.pairs().size()\n",
        "n_fobs_only = f_obs.size()-n_common\n",
        "n_rfree_only = r_free.size()-n_common\n",
        "print(\"Common:\",n_common)\n",
        "print(\"F-obs only:\",n_fobs_only)\n",
        "print(\"R-free only:\",n_rfree_only)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zidGAeSzY3N-"
      },
      "source": [
        "We just found out that the R-free-flags array has about 500 additional miller indices. Note that match_indices is a low level method and is typicially not used. It is nevertheless instructive for this example. Practically, to obtain a common set of reflections, the common_sets() method can be used:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zM0EphJKRheM"
      },
      "source": [
        "f_obs, r_free = f_obs.common_sets(r_free)\n",
        "print(f_obs.size())\n",
        "print(r_free.size())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6VX0WqUqVxIo"
      },
      "source": [
        "##F-model: Synthesis of model and data\n",
        "Let's go back to our structure, 1aba. We will include code for imports and to open the model and data mtz file together, to enable starting the examples here."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "d5ZzeNttWE_G"
      },
      "source": [
        "from iotbx.data_manager import DataManager\n",
        "\n",
        "dm = DataManager()\n",
        "filename_pdb = \"1aba_model.pdb\"\n",
        "dm.process_model_file(filename_pdb)\n",
        "model = dm.get_model(filename=filename_pdb)\n",
        "\n",
        "filename_mtz = \"/content/1aba_reflections.mtz\"\n",
        "dm.process_miller_array_file(filename=filename_mtz)\n",
        "reflections_reader = dm.get_miller_array(filename=filename_mtz)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OFeA776NR6_5"
      },
      "source": [
        "And select F-obs and R-free miller arrays"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "G6dPrV60WQYY"
      },
      "source": [
        "miller_array_dict = {ma.info().label_string():ma for ma in reflections_reader.as_miller_arrays()}\n",
        "print(miller_array_dict.keys())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6jmpSVIZWnnG"
      },
      "source": [
        "# select fobs and rfree\n",
        "f_obs = miller_array_dict['FOBS,SIGFOBS']\n",
        "r_free = miller_array_dict['FreeR_flag']\n",
        "f_obs, r_free = f_obs.common_sets(r_free)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iu9BJJjqSKdA"
      },
      "source": [
        "We will instantiate an F Model Manager object. This combines the data and the structure together. Note that the model we will used is an xray structure object. The xray structure object is a reduced form of the model, containing only the information about the scatterers which is needed for calculations."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tAvikit-P-gb"
      },
      "source": [
        "import mmtbx\n",
        "\n",
        "xrs = model.get_xray_structure()\n",
        "r_free_flags = r_free.array(data = r_free.data()==0)\n",
        "\n",
        "f_model = mmtbx.f_model.manager(\n",
        "  f_obs          = f_obs,\n",
        "  r_free_flags   = r_free_flags,\n",
        "  xray_structure = xrs)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5RIxriv0aeK2"
      },
      "source": [
        "The update_all_scales() method is necessary to calculate scale factors and add a bulk solvent model. This procedure puts the calculated and observed structure factor amplitudes to a common scale and adjusts calculated low resolution structure factors."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yv2fAp38R3iQ"
      },
      "source": [
        "_ = f_model.update_all_scales()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E-qpKMeHakhW"
      },
      "source": [
        "Now we will print out the R-work and R-free values"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RgkX5BqqXNao"
      },
      "source": [
        "print(\"R-work:\",f_model.r_work())\n",
        "print(\"R-free:\",f_model.r_free())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EkJKX-64bFtr"
      },
      "source": [
        "Now perform a fast-Fourier-transform (FFT) to convert the structure factors into a real map"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qW9P33t4ZVKH"
      },
      "source": [
        "f_calc = f_model.f_calc()\n",
        "fft_map = f_calc.fft_map()\n",
        "fft_map.apply_sigma_scaling()\n",
        "map_data = fft_map.real_map_unpadded()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bYkK-TlPb6Q7"
      },
      "source": [
        "Now we can instantiate a map manager manually from the three basic components:\n",
        "  - Map data\n",
        "  - Unit cell grid\n",
        "  - Unit cell crystal symmetry"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2ghE_VIMaI08"
      },
      "source": [
        "from iotbx.map_manager import map_manager\n",
        "\n",
        "mm=map_manager(map_data=map_data,\n",
        "                unit_cell_grid=map_data.all(), # the shape of the map array\n",
        "                unit_cell_crystal_symmetry=xrs.crystal_symmetry(),\n",
        "                wrapping=True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IPp-XpCMZTTr"
      },
      "source": [
        "from cctbx_py3dmol_wrapper import CCTBX3dmolWrapper\n",
        "show = CCTBX3dmolWrapper()\n",
        "show(mm)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bcg4sHUz5-bl"
      },
      "source": [
        "Of course, this is an F-calc map, so it is similar to what we got before using the generate_map() method. Let's try to calculate the 2mFo-DFc map."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zZHAW4h56bHo"
      },
      "source": [
        "from mmtbx import map_tools\n",
        "\n",
        "map_coeff = map_tools.electron_density_map(fmodel = f_model).map_coefficients(map_type=\"2mFo-DFc\")\n",
        "fft_map = map_coeff.fft_map()\n",
        "fft_map.apply_sigma_scaling()\n",
        "map_data = fft_map.real_map_unpadded()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OIe35R5Hc5yx"
      },
      "source": [
        "We can instantiate a Map Model Manager with the maps and the model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AYBA2K-Yahn0"
      },
      "source": [
        "from iotbx.map_manager import map_manager\n",
        "from iotbx.map_model_manager import map_model_manager\n",
        "\n",
        "mm_2fofc=map_manager(map_data=map_data,\n",
        "                unit_cell_grid=map_data.all(),\n",
        "                unit_cell_crystal_symmetry=xrs.crystal_symmetry(),\n",
        "                wrapping=True)\n",
        "\n",
        "mmm_2fofc = map_model_manager(model=model,map_manager=mm_2fofc)\n",
        "ligand_mmm_2fofc = mmm_2fofc.extract_all_maps_around_model(selection_string=\"resname MES\")\n",
        "show(ligand_mmm_2fofc,map_threshold=0.4)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "e0UnzkCVF_cR"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}