OpenFold.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/aqlaboratory/openfold/blob/main/notebooks/OpenFold.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pc5-mbsX9PZC"
      },
      "source": [
        "# OpenFold Colab\n",
        "\n",
        "Runs a simplified version of [OpenFold](https://github.com/aqlaboratory/openfold) on a target sequence. Adapted from DeepMind's [official AlphaFold Colab](https://colab.research.google.com/github/deepmind/alphafold/blob/main/notebooks/AlphaFold.ipynb).\n",
        "\n",
        "**Differences to AlphaFold v2.0**\n",
        "\n",
        "OpenFold is a trainable PyTorch reimplementation of AlphaFold 2. For the purposes of inference, it is practically identical to the original (\"practically\" because ensembling is excluded from OpenFold (recycling is enabled, however)).\n",
        "\n",
        "In this notebook, OpenFold is run with your choice of our original OpenFold parameters or DeepMind's publicly released parameters for AlphaFold 2.\n",
        "\n",
        "**Note**\n",
        "\n",
        "Like DeepMind's official Colab, this notebook uses **no templates (homologous structures)** and a selected portion of the full [BFD database](https://bfd.mmseqs.com/).\n",
        "\n",
        "**Citing this work**\n",
        "\n",
        "Any publication that discloses findings arising from using this notebook should [cite](https://github.com/deepmind/alphafold/#citing-this-work) DeepMind's [AlphaFold paper](https://doi.org/10.1038/s41586-021-03819-2).\n",
        "\n",
        "**Licenses**\n",
        "\n",
        "This Colab supports inference with the [AlphaFold model parameters](https://github.com/deepmind/alphafold/#model-parameters-license), made available under the Creative Commons Attribution 4.0 International ([CC BY 4.0](https://creativecommons.org/licenses/by/4.0/legalcode)) license. The Colab itself is provided under the [Apache 2.0 license](https://www.apache.org/licenses/LICENSE-2.0). See the full license statement below.\n",
        "\n",
        "**More information**\n",
        "\n",
        "You can find more information about how AlphaFold/OpenFold works in DeepMind's two Nature papers:\n",
        "\n",
        "*   [AlphaFold methods paper](https://www.nature.com/articles/s41586-021-03819-2)\n",
        "*   [AlphaFold predictions of the human proteome paper](https://www.nature.com/articles/s41586-021-03828-1)\n",
        "\n",
        "FAQ on how to interpret AlphaFold/OpenFold predictions are [here](https://alphafold.ebi.ac.uk/faq)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "rowN0bVYLe9n",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "8a4c9c7d-a555-460f-de39-d6aff2272d36"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Length of input sequence : 716\n"
          ]
        }
      ],
      "source": [
        "#@markdown ### Enter the amino acid sequence to fold ⬇️\n",
        "#@markdown For multiple sequences, separate sequences with a colon `:`\n",
        "input_sequence = 'MVDATRVPMDERFRTLKKKLEEGMVFTEYEQIPKKKANGIFSTAALPENAERSRIREVVPYEENRVELIPTKENNTGYINASHIKVVVGGAEWHYIATQGPLPHTCHDFWQMVWEQGVNVIAMVTAEEEGGRTKSHRYWPKLGSKHSSATYGKFKVTTKFRTDSVCYATTGLKVKHLLSGQERTVWHLQYTDWPDHGCPEDVQGFLSYLEEIQSVRRHTNSMLEGTKNRHPPIVVHCSAGVGRTGVLILSELMIYCLEHNEKVEVPMMLRLLREQRMFMIQTIAQYKFVYQVLIQFLQNSRLI:MVDATRVPMDERFRTLKKKLEEGMVFTEYEQIPKKKANGIFSTAALPENAERSRIREVVPYEENRVELIPTKENNTGYINASHIKVVVGGAEWHYIATQGPLPHTCHDFWQMVWEQGVNVIAMVTAEEEGGRTKSHRYWPKLGSKHSSATYGKFKVTTKFRTDSVCYATTGLKVKHLLSGQERTVWHLQYTDWPDHGCPEDVQGFLSYLEEIQSVRRHTNSMLEGTKNRHPPIVVHCSAGVGRTGVLILSELMIYCLEHNEKVEVPMMLRLLREQRMFMIQTIAQYKFVYQVLIQFLQNSRLI:GHMAEPQRHTMLCMCCKCEARIELVVESSADDLRAFQQLFLNTLSFVCPWCASQQ:GHMAEPQRHTMLCMCCKCEARIELVVESSADDLRAFQQLFLNTLSFVCPWCASQQ'  #@param {type:\"string\"}\n",
        "\n",
        "#@markdown ### Configure the model ⬇️\n",
        "\n",
        "weight_set = 'AlphaFold' #@param [\"OpenFold\", \"AlphaFold\"]\n",
        "model_mode = 'multimer' #@param [\"monomer\", \"multimer\"]\n",
        "relax_prediction = True #@param {type:\"boolean\"}\n",
        "\n",
        "\n",
        "# Remove all whitespaces, tabs and end lines; upper-case\n",
        "input_sequence = input_sequence.translate(str.maketrans('', '', ' \\n\\t')).upper()\n",
        "aatypes = set('ACDEFGHIKLMNPQRSTVWY')  # 20 standard aatypes\n",
        "allowed_chars = aatypes.union({':'})\n",
        "if not set(input_sequence).issubset(allowed_chars):\n",
        "  raise Exception(f'Input sequence contains non-amino acid letters: {set(input_sequence) - allowed_chars}. OpenFold only supports 20 standard amino acids as inputs.')\n",
        "\n",
        "if ':' in input_sequence and weight_set != 'AlphaFold':\n",
        "  raise ValueError('Input sequence is a multimer, must select Alphafold weight set')\n",
        "\n",
        "import enum\n",
        "@enum.unique\n",
        "class ModelType(enum.Enum):\n",
        "  MONOMER = 0\n",
        "  MULTIMER = 1\n",
        "\n",
        "model_type_dict = {\n",
        "    'monomer': ModelType.MONOMER,\n",
        "    'multimer': ModelType.MULTIMER,\n",
        "}\n",
        "\n",
        "model_type = model_type_dict[model_mode]\n",
        "print(f'Length of input sequence : {len(input_sequence.replace(\":\", \"\"))}')\n",
        "#@markdown After making your selections, execute this cell by pressing the\n",
        "#@markdown *Play* button on the left."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "woIxeCPygt7K"
      },
      "outputs": [],
      "source": [
        "#@title Install third-party software\n",
        "#@markdown Please execute this cell by pressing the *Play* button on\n",
        "#@markdown the left.\n",
        "\n",
        "\n",
        "#@markdown **Note**: This installs the software on the Colab\n",
        "#@markdown notebook in the cloud and not on your computer.\n",
        "\n",
        "import os, time\n",
        "from IPython.utils import io\n",
        "from sys import version_info\n",
        "import subprocess\n",
        "\n",
        "python_version = f\"{version_info.major}.{version_info.minor}\"\n",
        "\n",
        "\n",
        "os.system(\"wget -qnc https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-Linux-x86_64.sh\")\n",
        "os.system(\"bash Mambaforge-Linux-x86_64.sh -bfp /usr/local\")\n",
        "os.system(\"mamba config --set auto_update_conda false\")\n",
        "os.system(f\"mamba install -y -c conda-forge -c bioconda kalign2=2.04 hhsuite=3.3.0 openmm=7.7.0 python={python_version} pdbfixer biopython=1.79\")\n",
        "os.system(\"pip install -q torch ml_collections py3Dmol modelcif\")\n",
        "\n",
        "try:\n",
        "  with io.capture_output() as captured:\n",
        "\n",
        "    # Create a ramdisk to store a database chunk to make Jackhmmer run fast.\n",
        "    %shell sudo apt install --quiet --yes hmmer\n",
        "    %shell sudo mkdir -m 777 --parents /tmp/ramdisk\n",
        "    %shell sudo mount -t tmpfs -o size=9G ramdisk /tmp/ramdisk\n",
        "\n",
        "    %shell wget -q -P /content \\\n",
        "      https://git.scicore.unibas.ch/schwede/openstructure/-/raw/7102c63615b64735c4941278d92b554ec94415f8/modules/mol/alg/src/stereo_chemical_props.txt\n",
        "\n",
        "    %shell mkdir -p /content/openfold/openfold/resources\n",
        "\n",
        "    commit = \"e2e19f16676b1a409f9ba3a6f69b11ee7f5887c2\"\n",
        "    os.system(f\"pip install -q git+https://github.com/aqlaboratory/openfold.git@{commit}\")\n",
        "\n",
        "    os.system(f\"cp -f -p /content/stereo_chemical_props.txt /usr/local/lib/python{python_version}/site-packages/openfold/resources/\")\n",
        "\n",
        "except subprocess.CalledProcessError as captured:\n",
        "  print(captured)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "VzJ5iMjTtoZw"
      },
      "outputs": [],
      "source": [
        "#@title Download model weights\n",
        "#@markdown Please execute this cell by pressing the *Play* button on\n",
        "#@markdown the left.\n",
        "\n",
        "# Define constants\n",
        "GIT_REPO='https://github.com/aqlaboratory/openfold'\n",
        "ALPHAFOLD_PARAM_SOURCE_URL = 'https://storage.googleapis.com/alphafold/alphafold_params_2022-12-06.tar'\n",
        "OPENFOLD_PARAMS_DIR = './openfold/openfold/resources/openfold_params'\n",
        "ALPHAFOLD_PARAMS_DIR = './openfold/openfold/resources/params'\n",
        "ALPHAFOLD_PARAMS_PATH = os.path.join(\n",
        "  ALPHAFOLD_PARAMS_DIR, os.path.basename(ALPHAFOLD_PARAM_SOURCE_URL)\n",
        ")\n",
        "\n",
        "try:\n",
        "  with io.capture_output() as captured:\n",
        "    if(weight_set == 'AlphaFold'):\n",
        "      %shell mkdir --parents \"{ALPHAFOLD_PARAMS_DIR}\"\n",
        "      %shell wget -O {ALPHAFOLD_PARAMS_PATH} {ALPHAFOLD_PARAM_SOURCE_URL}\n",
        "      %shell tar --extract --verbose --file=\"{ALPHAFOLD_PARAMS_PATH}\" \\\n",
        "        --directory=\"{ALPHAFOLD_PARAMS_DIR}\" --preserve-permissions\n",
        "      %shell rm \"{ALPHAFOLD_PARAMS_PATH}\"\n",
        "    elif(weight_set == 'OpenFold'):\n",
        "      # Install AWS CLI\n",
        "      %shell curl \"https://awscli.amazonaws.com/awscli-exe-linux-x86_64.zip\" -o \"awscliv2.zip\"\n",
        "      %shell unzip -qq awscliv2.zip\n",
        "      %shell sudo ./aws/install\n",
        "      %shell rm awscliv2.zip\n",
        "      %shell rm -rf ./aws\n",
        "      %shell mkdir --parents \"{OPENFOLD_PARAMS_DIR}\"\n",
        "\n",
        "      %shell aws s3 cp \\\n",
        "        --no-sign-request \\\n",
        "        --region us-east-1 \\\n",
        "        s3://openfold/openfold_params \"{OPENFOLD_PARAMS_DIR}\" \\\n",
        "        --recursive\n",
        "    else:\n",
        "      raise ValueError(\"Invalid weight set\")\n",
        "except subprocess.CalledProcessError as captured:\n",
        "  print(captured)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "_FpxxMo-mvcP"
      },
      "outputs": [],
      "source": [
        "#@title Import Python packages\n",
        "#@markdown Please execute this cell by pressing the *Play* button on\n",
        "#@markdown the left.\n",
        "\n",
        "import unittest.mock\n",
        "import sys\n",
        "from typing import Dict, Sequence\n",
        "\n",
        "sys.path.insert(0, f'/usr/local/lib/python{python_version}/dist-packages/')\n",
        "sys.path.insert(0, f'/usr/local/lib/python{python_version}/site-packages/')\n",
        "\n",
        "# Allows us to skip installing these packages\n",
        "unnecessary_modules = [\n",
        "  \"dllogger\",\n",
        "  \"pytorch_lightning\",\n",
        "  \"pytorch_lightning.utilities\",\n",
        "  \"pytorch_lightning.callbacks.early_stopping\",\n",
        "  \"pytorch_lightning.utilities.seed\",\n",
        "]\n",
        "for unnecessary_module in unnecessary_modules:\n",
        "  sys.modules[unnecessary_module] = unittest.mock.MagicMock()\n",
        "\n",
        "import os\n",
        "\n",
        "from urllib import request\n",
        "from concurrent import futures\n",
        "from google.colab import files\n",
        "import json\n",
        "from matplotlib import gridspec\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import py3Dmol\n",
        "import torch\n",
        "import shutil\n",
        "import tqdm\n",
        "import tqdm.notebook\n",
        "\n",
        "TQDM_BAR_FORMAT = '{l_bar}{bar}| {n_fmt}/{total_fmt} [elapsed: {elapsed} remaining: {remaining}]'\n",
        "\n",
        "# Prevent shell magic being broken by openmm, prevent this cryptic error:\n",
        "# \"NotImplementedError: A UTF-8 locale is required. Got ANSI_X3.4-1968\"\n",
        "import locale\n",
        "def getpreferredencoding(do_setlocale = True):\n",
        "    return \"UTF-8\"\n",
        "locale.getpreferredencoding = getpreferredencoding\n",
        "\n",
        "from openfold import config\n",
        "from openfold.data import feature_pipeline\n",
        "from openfold.data import parsers\n",
        "from openfold.data import data_pipeline\n",
        "from openfold.data import msa_pairing\n",
        "from openfold.data import feature_processing_multimer\n",
        "from openfold.data.tools import jackhmmer\n",
        "from openfold.model import model\n",
        "from openfold.np import protein\n",
        "from openfold.np.relax import relax\n",
        "from openfold.np.relax.utils import overwrite_b_factors\n",
        "from openfold.utils.import_weights import import_jax_weights_\n",
        "from openfold.utils.tensor_utils import tensor_tree_map\n",
        "\n",
        "from IPython import display\n",
        "from ipywidgets import GridspecLayout\n",
        "from ipywidgets import Output"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W4JpOs6oA-QS"
      },
      "source": [
        "## Making a prediction\n",
        "\n",
        "Note that the search against databases and the actual prediction can take some time, from minutes to hours, depending on the length of the protein and what type of GPU you are allocated by Colab (see FAQ below)."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Search against genetic databases\n",
        "\n",
        "#@markdown Once this cell has been executed, you will see\n",
        "#@markdown statistics about the multiple sequence alignment\n",
        "#@markdown (MSA) that will be used by OpenFold. In particular,\n",
        "#@markdown you’ll see how well each residue is covered by similar\n",
        "#@markdown sequences in the MSA.\n",
        "\n",
        "# --- Find the closest source --\n",
        "test_url_pattern = 'https://storage.googleapis.com/alphafold-colab{:s}/latest/uniref90_2021_03.fasta.1'\n",
        "ex = futures.ThreadPoolExecutor(3)\n",
        "def fetch(source):\n",
        "  request.urlretrieve(test_url_pattern.format(source))\n",
        "  return source\n",
        "fs = [ex.submit(fetch, source) for source in ['', '-europe', '-asia']]\n",
        "source = None\n",
        "for f in futures.as_completed(fs):\n",
        "  source = f.result()\n",
        "  ex.shutdown()\n",
        "  break\n",
        "\n",
        "# Run the search against chunks of genetic databases (since the genetic\n",
        "# databases don't fit in Colab ramdisk).\n",
        "\n",
        "jackhmmer_binary_path = '/usr/bin/jackhmmer'\n",
        "\n",
        "# --- Parse multiple sequences, if there are any ---\n",
        "def split_multiple_sequences(sequence):\n",
        "  seqs = sequence.split(':')\n",
        "  sorted_seqs = sorted(seqs, key=lambda s: len(s))\n",
        "\n",
        "  # TODO: Handle the homomer case when writing fasta sequences\n",
        "  fasta_path_tuples = []\n",
        "  for idx, seq in enumerate(set(sorted_seqs)):\n",
        "    fasta_path = f'target_{idx+1}.fasta'\n",
        "    with open(fasta_path, 'wt') as f:\n",
        "      f.write(f'>query\\n{seq}\\n')\n",
        "    fasta_path_tuples.append((seq, fasta_path))\n",
        "  fasta_path_by_seq = dict(fasta_path_tuples)\n",
        "\n",
        "  return sorted_seqs, fasta_path_by_seq\n",
        "\n",
        "sequences, fasta_path_by_sequence = split_multiple_sequences(input_sequence)\n",
        "db_results_by_sequence = {seq: {} for seq in fasta_path_by_sequence.keys()}\n",
        "\n",
        "DB_ROOT_PATH = f'https://storage.googleapis.com/alphafold-colab{source}/latest/'\n",
        "db_configs = {}\n",
        "db_configs['smallbfd'] = {\n",
        "    'database_path': f'{DB_ROOT_PATH}uniref90_2021_03.fasta',\n",
        "    'z_value': 65984053,\n",
        "    'num_jackhmmer_chunks': 17,\n",
        "}\n",
        "db_configs['mgnify'] = {\n",
        "    'database_path': f'{DB_ROOT_PATH}mgy_clusters_2022_05.fasta',\n",
        "    'z_value': 304820129,\n",
        "    'num_jackhmmer_chunks': 120,\n",
        "}\n",
        "db_configs['uniref90'] = {\n",
        "     'database_path': f'{DB_ROOT_PATH}uniref90_2022_01.fasta',\n",
        "     'z_value': 144113457,\n",
        "     'num_jackhmmer_chunks': 62,\n",
        "}\n",
        "\n",
        "# Search UniProt and construct the all_seq features only for heteromers, not homomers.\n",
        "if model_type == ModelType.MULTIMER and len(set(sequences)) > 1:\n",
        "  db_configs['uniprot'] = {\n",
        "       'database_path': f'{DB_ROOT_PATH}uniprot_2021_04.fasta',\n",
        "       'z_value': 225013025 + 565928,\n",
        "       'num_jackhmmer_chunks': 101,\n",
        "       }\n",
        "\n",
        "total_jackhmmer_chunks = sum([d['num_jackhmmer_chunks'] for d in db_configs.values()])\n",
        "with tqdm.notebook.tqdm(total=total_jackhmmer_chunks, bar_format=TQDM_BAR_FORMAT) as pbar:\n",
        "  def jackhmmer_chunk_callback(i):\n",
        "    pbar.update(n=1)\n",
        "\n",
        "  for db_name, db_config in db_configs.items():\n",
        "    pbar.set_description(f'Searching {db_name}')\n",
        "    jackhmmer_runner = jackhmmer.Jackhmmer(\n",
        "        binary_path=jackhmmer_binary_path,\n",
        "        database_path=db_config['database_path'],\n",
        "        get_tblout=True,\n",
        "        num_streamed_chunks=db_config['num_jackhmmer_chunks'],\n",
        "        streaming_callback=jackhmmer_chunk_callback,\n",
        "        z_value=db_config['z_value'])\n",
        "\n",
        "    db_results = jackhmmer_runner.query_multiple(fasta_path_by_sequence.values())\n",
        "    for seq, result in zip(fasta_path_by_sequence.keys(), db_results):\n",
        "      db_results_by_sequence[seq][db_name] = result\n",
        "\n",
        "\n",
        "# --- Extract the MSAs and visualize ---\n",
        "# Extract the MSAs from the Stockholm files.\n",
        "# NB: deduplication happens later in data_pipeline.make_msa_features.\n",
        "\n",
        "MAX_HITS_BY_DB = {\n",
        "    'uniref90': 10000,\n",
        "    'smallbfd': 5000,\n",
        "    'mgnify': 501,\n",
        "    'uniprot': 50000,\n",
        "}\n",
        "\n",
        "msas_by_seq_by_db = {seq: {} for seq in sequences}\n",
        "full_msa_by_seq = {seq: [] for seq in sequences}\n",
        "\n",
        "for seq, sequence_result in db_results_by_sequence.items():\n",
        "  print(f'parsing_results_for_sequence {seq}')\n",
        "  for db_name, db_results in sequence_result.items():\n",
        "    unsorted_results = []\n",
        "    for i, result in enumerate(db_results):\n",
        "      msa_obj = parsers.parse_stockholm(result['sto'])\n",
        "      e_values_dict = parsers.parse_e_values_from_tblout(result['tbl'])\n",
        "      target_names = msa_obj.descriptions\n",
        "      e_values = [e_values_dict[t.split('/')[0]] for t in target_names]\n",
        "      zipped_results = zip(msa_obj.sequences, msa_obj.deletion_matrix, target_names, e_values)\n",
        "      if i != 0:\n",
        "        # Only take query from the first chunk\n",
        "        zipped_results = [x for x in zipped_results if x[2] != 'query']\n",
        "      unsorted_results.extend(zipped_results)\n",
        "    sorted_by_evalue = sorted(unsorted_results, key=lambda x: x[3])\n",
        "    msas, del_matrix, targets, _ = zip(*sorted_by_evalue)\n",
        "    db_msas = parsers.Msa(msas, del_matrix, targets)\n",
        "    if db_msas:\n",
        "      if db_name in MAX_HITS_BY_DB:\n",
        "        db_msas.truncate(MAX_HITS_BY_DB[db_name])\n",
        "      msas_by_seq_by_db[seq][db_name] = db_msas\n",
        "      full_msa_by_seq[seq].extend(db_msas.sequences)\n",
        "      msa_size = len(set(db_msas.sequences))\n",
        "      print(f'{msa_size} Sequences Found in {db_name}')\n",
        "\n",
        "\n",
        "fig = plt.figure(figsize=(12, 3))\n",
        "max_num_alignments = 0\n",
        "\n",
        "for seq_idx, seq in enumerate(set(sequences)):\n",
        "  full_msas = full_msa_by_seq[seq]\n",
        "  deduped_full_msa = list(dict.fromkeys(full_msas))\n",
        "  total_msa_size = len(deduped_full_msa)\n",
        "  print(f'\\n{total_msa_size} Sequences Found in Total\\n')\n",
        "\n",
        "  aa_map = {restype: i for i, restype in enumerate('ABCDEFGHIJKLMNOPQRSTUVWXYZ-')}\n",
        "  msa_arr = np.array([[aa_map[aa] for aa in seq] for seq in deduped_full_msa])\n",
        "  num_alignments, num_res = msa_arr.shape\n",
        "  plt.plot(np.sum(msa_arr != aa_map['-'], axis=0), label=f'Chain {seq_idx}')\n",
        "  max_num_alignments = max(num_alignments, max_num_alignments)\n",
        "\n",
        "\n",
        "plt.title('Per-Residue Count of Non-Gap Amino Acids in the MSA')\n",
        "plt.ylabel('Non-Gap Count')\n",
        "plt.yticks(range(0, max_num_alignments + 1, max(1, int(max_num_alignments / 3))))\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "o7BqQN_gfYtq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XUo6foMQxwS2"
      },
      "outputs": [],
      "source": [
        "#@title Run OpenFold and download prediction\n",
        "\n",
        "#@markdown Once this cell has been executed, a zip-archive with\n",
        "#@markdown the obtained prediction will be automatically downloaded\n",
        "#@markdown to your computer.\n",
        "\n",
        "# Color bands for visualizing plddt\n",
        "PLDDT_BANDS = [\n",
        "  (0, 50, '#FF7D45'),\n",
        "  (50, 70, '#FFDB13'),\n",
        "  (70, 90, '#65CBF3'),\n",
        "  (90, 100, '#0053D6')\n",
        "]\n",
        "\n",
        "# --- Run the model ---\n",
        "if model_type == ModelType.MONOMER:\n",
        "  model_names = [\n",
        "   'finetuning_3.pt',\n",
        "   'finetuning_4.pt',\n",
        "   'finetuning_5.pt',\n",
        "   'finetuning_ptm_2.pt',\n",
        "   'finetuning_no_templ_ptm_1.pt'\n",
        "  ]\n",
        "elif model_type == ModelType.MULTIMER:\n",
        "  model_names = [\n",
        "    'model_1_multimer_v3',\n",
        "    'model_2_multimer_v3',\n",
        "    'model_3_multimer_v3',\n",
        "    'model_4_multimer_v3',\n",
        "    'model_5_multimer_v3',\n",
        "  ]\n",
        "\n",
        "def _placeholder_template_feats(num_templates_, num_res_):\n",
        "  return {\n",
        "      'template_aatype': np.zeros((num_templates_, num_res_, 22), dtype=np.int64),\n",
        "      'template_all_atom_positions': np.zeros((num_templates_, num_res_, 37, 3), dtype=np.float32),\n",
        "      'template_all_atom_mask': np.zeros((num_templates_, num_res_, 37), dtype=np.float32),\n",
        "      'template_domain_names': np.zeros((num_templates_,), dtype=np.float32),\n",
        "      'template_sum_probs': np.zeros((num_templates_, 1), dtype=np.float32),\n",
        "  }\n",
        "\n",
        "\n",
        "def make_features(\n",
        "    sequences: Sequence[str],\n",
        "    msas_by_seq_by_db: Dict[str, Dict[str, parsers.Msa]],\n",
        "    model_type: ModelType):\n",
        "  num_templates = 1 # Placeholder for generating fake template features\n",
        "  feature_dict = {}\n",
        "\n",
        "  for idx, seq in enumerate(sequences, start=1):\n",
        "    _chain_id = f'chain_{idx}'\n",
        "    num_res = len(seq)\n",
        "\n",
        "    feats = data_pipeline.make_sequence_features(seq, _chain_id, num_res)\n",
        "    msas_without_uniprot = [msas_by_seq_by_db[seq][db] for db in db_configs.keys() if db != 'uniprot']\n",
        "    msa_feats = data_pipeline.make_msa_features(msas_without_uniprot)\n",
        "    feats.update(msa_feats)\n",
        "    feats.update(_placeholder_template_feats(num_templates, num_res))\n",
        "\n",
        "    if model_type == ModelType.MONOMER:\n",
        "      feature_dict[seq] = feats\n",
        "    if model_type == ModelType.MULTIMER:\n",
        "      # Perform extra pair processing steps for heteromers\n",
        "      if len(set(sequences)) > 1:\n",
        "        uniprot_msa = msas_by_seq_by_db[seq]['uniprot']\n",
        "        uniprot_msa_features = data_pipeline.make_msa_features([uniprot_msa])\n",
        "        valid_feat_names = msa_pairing.MSA_FEATURES + (\n",
        "                'msa_species_identifiers',\n",
        "            )\n",
        "        pair_feats = {\n",
        "                f'{k}_all_seq': v for k, v in uniprot_msa_features.items()\n",
        "                if k in valid_feat_names\n",
        "            }\n",
        "        feats.update(pair_feats)\n",
        "\n",
        "      feats = data_pipeline.convert_monomer_features(feats, _chain_id)\n",
        "      feature_dict[_chain_id] = feats\n",
        "\n",
        "  if model_type == ModelType.MONOMER:\n",
        "    np_example = feature_dict[sequences[0]]\n",
        "  elif model_type == ModelType.MULTIMER:\n",
        "    all_chain_feats = data_pipeline.add_assembly_features(feature_dict)\n",
        "    features = feature_processing_multimer.pair_and_merge(all_chain_features=all_chain_feats)\n",
        "    np_example = data_pipeline.pad_msa(features, 512)\n",
        "\n",
        "  return np_example\n",
        "\n",
        "\n",
        "output_dir = 'prediction'\n",
        "os.makedirs(output_dir, exist_ok=True)\n",
        "\n",
        "plddts = {}\n",
        "pae_outputs = {}\n",
        "weighted_ptms = {}\n",
        "unrelaxed_proteins = {}\n",
        "\n",
        "with tqdm.notebook.tqdm(total=len(model_names), bar_format=TQDM_BAR_FORMAT) as pbar:\n",
        "  for i, model_name in enumerate(model_names, start = 1):\n",
        "    pbar.set_description(f'Running {model_name}')\n",
        "\n",
        "    feature_dict = make_features(sequences, msas_by_seq_by_db, model_type)\n",
        "\n",
        "    if(weight_set == \"AlphaFold\"):\n",
        "      if model_type == ModelType.MONOMER:\n",
        "        config_preset = f\"model_{i}\"\n",
        "      elif model_type == ModelType.MULTIMER:\n",
        "        config_preset = f'model_{i}_multimer_v3'\n",
        "    else:\n",
        "      if(\"_no_templ_\" in model_name):\n",
        "        config_preset = \"model_3\"\n",
        "      else:\n",
        "        config_preset = \"model_1\"\n",
        "      if(\"_ptm_\" in model_name):\n",
        "        config_preset += \"_ptm\"\n",
        "\n",
        "    cfg = config.model_config(config_preset)\n",
        "\n",
        "    # Force the model to only use 3 recycling updates\n",
        "    cfg.data.common.max_recycling_iters = 3\n",
        "    cfg.model.recycle_early_stop_tolerance = -1\n",
        "\n",
        "    openfold_model = model.AlphaFold(cfg)\n",
        "    openfold_model = openfold_model.eval()\n",
        "    if(weight_set == \"AlphaFold\"):\n",
        "      params_name = os.path.join(\n",
        "        ALPHAFOLD_PARAMS_DIR, f\"params_{config_preset}.npz\"\n",
        "      )\n",
        "      import_jax_weights_(openfold_model, params_name, version=config_preset)\n",
        "    elif(weight_set == \"OpenFold\"):\n",
        "      params_name = os.path.join(\n",
        "        OPENFOLD_PARAMS_DIR,\n",
        "        model_name,\n",
        "      )\n",
        "      d = torch.load(params_name)\n",
        "      openfold_model.load_state_dict(d)\n",
        "    else:\n",
        "      raise ValueError(f\"Invalid weight set: {weight_set}\")\n",
        "\n",
        "    openfold_model = openfold_model.cuda()\n",
        "\n",
        "    pipeline = feature_pipeline.FeaturePipeline(cfg.data)\n",
        "    processed_feature_dict = pipeline.process_features(\n",
        "      feature_dict,\n",
        "      mode='predict',\n",
        "      is_multimer = (model_type == ModelType.MULTIMER),\n",
        "    )\n",
        "\n",
        "    processed_feature_dict = tensor_tree_map(\n",
        "        lambda t: t.cuda(), processed_feature_dict\n",
        "    )\n",
        "\n",
        "    with torch.no_grad():\n",
        "      prediction_result = openfold_model(processed_feature_dict)\n",
        "\n",
        "    # Move the batch and output to np for further processing\n",
        "    processed_feature_dict = tensor_tree_map(\n",
        "      lambda t: np.array(t[..., -1].cpu()), processed_feature_dict\n",
        "    )\n",
        "    prediction_result = tensor_tree_map(\n",
        "      lambda t: np.array(t.cpu()), prediction_result\n",
        "    )\n",
        "\n",
        "    mean_plddt = prediction_result['plddt'].mean()\n",
        "\n",
        "    if model_type == ModelType.MONOMER:\n",
        "      if 'predicted_aligned_error' in prediction_result:\n",
        "        pae_outputs[model_name] = (\n",
        "            prediction_result['predicted_aligned_error'],\n",
        "            prediction_result['max_predicted_aligned_error']\n",
        "        )\n",
        "      else:\n",
        "        # Get the pLDDT confidence metrics. Do not put pTM models here as they\n",
        "        # should never get selected.\n",
        "        plddts[model_name] = prediction_result['plddt']\n",
        "    elif model_type == ModelType.MULTIMER:\n",
        "      # Multimer models are sorted by pTM+ipTM.\n",
        "      plddts[model_name] = prediction_result['plddt']\n",
        "      pae_outputs[model_name] = (prediction_result['predicted_aligned_error'],\n",
        "                                 prediction_result['max_predicted_aligned_error'])\n",
        "\n",
        "      weighted_ptms[model_name] = prediction_result['weighted_ptm_score']\n",
        "\n",
        "    # Set the b-factors to the per-residue plddt.\n",
        "    final_atom_mask = prediction_result['final_atom_mask']\n",
        "    b_factors = prediction_result['plddt'][:, None] * final_atom_mask\n",
        "    unrelaxed_protein = protein.from_prediction(\n",
        "      processed_feature_dict,\n",
        "      prediction_result,\n",
        "      remove_leading_feature_dimension=False,\n",
        "      b_factors=b_factors,\n",
        "    )\n",
        "    unrelaxed_proteins[model_name] = unrelaxed_protein\n",
        "\n",
        "    # Delete unused outputs to save memory.\n",
        "    del openfold_model\n",
        "    del processed_feature_dict\n",
        "    del prediction_result\n",
        "    pbar.update(n=1)\n",
        "\n",
        "  # Find the best model according to the mean pLDDT.\n",
        "  if model_type == ModelType.MONOMER:\n",
        "    best_model_name = max(plddts.keys(), key=lambda x: plddts[x].mean())\n",
        "  elif model_type == ModelType.MULTIMER:\n",
        "    best_model_name = max(weighted_ptms.keys(), key=lambda x: weighted_ptms[x])\n",
        "  best_pdb = protein.to_pdb(unrelaxed_proteins[best_model_name])\n",
        "\n",
        "  # --- AMBER relax the best model ---\n",
        "  if(relax_prediction):\n",
        "    pbar.set_description(f'AMBER relaxation')\n",
        "    amber_relaxer = relax.AmberRelaxation(\n",
        "        max_iterations=0,\n",
        "        tolerance=2.39,\n",
        "        stiffness=10.0,\n",
        "        exclude_residues=[],\n",
        "        max_outer_iterations=20,\n",
        "        use_gpu=True,\n",
        "    )\n",
        "    relaxed_pdb, _, _ = amber_relaxer.process(\n",
        "        prot=unrelaxed_proteins[best_model_name]\n",
        "    )\n",
        "    best_pdb = relaxed_pdb\n",
        "\n",
        "  # Write out the prediction\n",
        "  pred_output_path = os.path.join(output_dir, 'selected_prediction.pdb')\n",
        "  with open(pred_output_path, 'w') as f:\n",
        "    f.write(best_pdb)\n",
        "\n",
        "  pbar.update(n=1)  # Finished AMBER relax.\n",
        "\n",
        "# Construct multiclass b-factors to indicate confidence bands\n",
        "# 0=very low, 1=low, 2=confident, 3=very high\n",
        "banded_b_factors = []\n",
        "for plddt in plddts[best_model_name]:\n",
        "  for idx, (min_val, max_val, _) in enumerate(PLDDT_BANDS):\n",
        "    if plddt >= min_val and plddt <= max_val:\n",
        "      banded_b_factors.append(idx)\n",
        "      break\n",
        "banded_b_factors = np.array(banded_b_factors)[:, None] * final_atom_mask\n",
        "to_visualize_pdb = overwrite_b_factors(best_pdb, banded_b_factors)\n",
        "\n",
        "# --- Visualise the prediction & confidence ---\n",
        "show_sidechains = True\n",
        "def plot_plddt_legend():\n",
        "  \"\"\"Plots the legend for pLDDT.\"\"\"\n",
        "  thresh = [\n",
        "            'Very low (pLDDT < 50)',\n",
        "            'Low (70 > pLDDT > 50)',\n",
        "            'Confident (90 > pLDDT > 70)',\n",
        "            'Very high (pLDDT > 90)']\n",
        "\n",
        "  colors = [x[2] for x in PLDDT_BANDS]\n",
        "\n",
        "  plt.figure(figsize=(2, 2))\n",
        "  for c in colors:\n",
        "    plt.bar(0, 0, color=c)\n",
        "  plt.legend(thresh, frameon=False, loc='center', fontsize=20)\n",
        "  plt.xticks([])\n",
        "  plt.yticks([])\n",
        "  ax = plt.gca()\n",
        "  ax.spines['right'].set_visible(False)\n",
        "  ax.spines['top'].set_visible(False)\n",
        "  ax.spines['left'].set_visible(False)\n",
        "  ax.spines['bottom'].set_visible(False)\n",
        "  plt.title('Model Confidence', fontsize=20, pad=20)\n",
        "  return plt\n",
        "\n",
        "# Show the structure coloured by chain if the multimer model has been used.\n",
        "if model_type == ModelType.MULTIMER:\n",
        "  multichain_view = py3Dmol.view(width=800, height=600)\n",
        "  multichain_view.addModelsAsFrames(to_visualize_pdb)\n",
        "  multichain_style = {'cartoon': {'colorscheme': 'chain'}}\n",
        "  multichain_view.setStyle({'model': -1}, multichain_style)\n",
        "  multichain_view.zoomTo()\n",
        "  multichain_view.show()\n",
        "\n",
        "# Color the structure by per-residue pLDDT\n",
        "color_map = {i: bands[2] for i, bands in enumerate(PLDDT_BANDS)}\n",
        "view = py3Dmol.view(width=800, height=600)\n",
        "view.addModelsAsFrames(to_visualize_pdb)\n",
        "style = {'cartoon': {\n",
        "    'colorscheme': {\n",
        "        'prop': 'b',\n",
        "        'map': color_map}\n",
        "        }}\n",
        "if show_sidechains:\n",
        "  style['stick'] = {}\n",
        "view.setStyle({'model': -1}, style)\n",
        "view.zoomTo()\n",
        "\n",
        "grid = GridspecLayout(1, 2)\n",
        "out = Output()\n",
        "with out:\n",
        "  view.show()\n",
        "grid[0, 0] = out\n",
        "\n",
        "out = Output()\n",
        "with out:\n",
        "  plot_plddt_legend().show()\n",
        "grid[0, 1] = out\n",
        "\n",
        "display.display(grid)\n",
        "\n",
        "# Display pLDDT and predicted aligned error (if output by the model).\n",
        "if pae_outputs:\n",
        "  num_plots = 2\n",
        "else:\n",
        "  num_plots = 1\n",
        "\n",
        "plt.figure(figsize=[8 * num_plots, 6])\n",
        "plt.subplot(1, num_plots, 1)\n",
        "plt.plot(plddts[best_model_name])\n",
        "plt.title('Predicted LDDT')\n",
        "plt.xlabel('Residue')\n",
        "plt.ylabel('pLDDT')\n",
        "\n",
        "if num_plots == 2:\n",
        "  plt.subplot(1, 2, 2)\n",
        "  pae, max_pae = list(pae_outputs.values())[0]\n",
        "  plt.imshow(pae, vmin=0., vmax=max_pae, cmap='Greens_r')\n",
        "  plt.colorbar(fraction=0.046, pad=0.04)\n",
        "\n",
        "  # Display lines at chain boundaries.\n",
        "  best_unrelaxed_prot = unrelaxed_proteins[best_model_name]\n",
        "  total_num_res = best_unrelaxed_prot.residue_index.shape[-1]\n",
        "  chain_ids = best_unrelaxed_prot.chain_index\n",
        "  for chain_boundary in np.nonzero(chain_ids[:-1] - chain_ids[1:]):\n",
        "    if chain_boundary.size:\n",
        "      plt.plot([0, total_num_res], [chain_boundary, chain_boundary], color='red')\n",
        "      plt.plot([chain_boundary, chain_boundary], [0, total_num_res], color='red')\n",
        "  plt.title('Predicted Aligned Error')\n",
        "  plt.xlabel('Scored residue')\n",
        "  plt.ylabel('Aligned residue')\n",
        "\n",
        "# Save pLDDT and predicted aligned error (if it exists)\n",
        "pae_output_path = os.path.join(output_dir, 'predicted_aligned_error.json')\n",
        "if pae_outputs:\n",
        "  # Save predicted aligned error in the same format as the AF EMBL DB\n",
        "  rounded_errors = np.round(pae.astype(np.float64), decimals=1)\n",
        "  indices = np.indices((len(rounded_errors), len(rounded_errors))) + 1\n",
        "  indices_1 = indices[0].flatten().tolist()\n",
        "  indices_2 = indices[1].flatten().tolist()\n",
        "  pae_data = json.dumps([{\n",
        "      'residue1': indices_1,\n",
        "      'residue2': indices_2,\n",
        "      'distance': rounded_errors.flatten().tolist(),\n",
        "      'max_predicted_aligned_error': max_pae.item()\n",
        "  }],\n",
        "                        indent=None,\n",
        "                        separators=(',', ':'))\n",
        "  with open(pae_output_path, 'w') as f:\n",
        "    f.write(pae_data)\n",
        "\n",
        "\n",
        "# --- Download the predictions ---\n",
        "shutil.make_archive(base_name='prediction', format='zip', root_dir=output_dir)\n",
        "files.download(f'{output_dir}.zip')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lUQAn5LYC5n4"
      },
      "source": [
        "### Interpreting the prediction\n",
        "\n",
        "Please see the [AlphaFold methods paper](https://www.nature.com/articles/s41586-021-03819-2) and the [AlphaFold predictions of the human proteome paper](https://www.nature.com/articles/s41586-021-03828-1), as well as [DeepMind's FAQ](https://alphafold.ebi.ac.uk/faq) on how to interpret AlphaFold/OpenFold predictions. More information about the predictions of the AlphaFold Multimer model can be found in the [Alphafold Multimer paper](https://www.biorxiv.org/content/10.1101/2022.03.11.484043v3.full.pdf)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jeb2z8DIA4om"
      },
      "source": [
        "## FAQ & Troubleshooting\n",
        "\n",
        "\n",
        "*   How do I get a predicted protein structure for my protein?\n",
        "    *   Click on the _Connect_ button on the top right to get started.\n",
        "    *   Paste the amino acid sequence of your protein (without any headers) into the “Enter the amino acid sequence to fold”.\n",
        "    *   Run all cells in the Colab, either by running them individually (with the play button on the left side) or via _Runtime_ > _Run all._\n",
        "    *   The predicted protein structure will be downloaded once all cells have been executed. Note: This can take minutes to hours - see below.\n",
        "*   How long will this take?\n",
        "    *   Downloading the OpenFold source code can take up to a few minutes.\n",
        "    *   Downloading and installing the third-party software can take up to a few minutes.\n",
        "    *   The search against genetic databases can take minutes to hours.\n",
        "    *   Running OpenFold and generating the prediction can take minutes to hours, depending on the length of your protein and on which GPU-type Colab has assigned you.\n",
        "*   My Colab no longer seems to be doing anything, what should I do?\n",
        "    *   Some steps may take minutes to hours to complete.\n",
        "    *   Sometimes, running the \"installation\" cells more than once can corrupt the OpenFold installation.\n",
        "    *   If nothing happens or if you receive an error message, try restarting your Colab runtime via _Runtime_ > _Restart runtime_.\n",
        "    *   If this doesn’t help, reset your Colab runtime via _Runtime_ > _Factory reset runtime_.\n",
        "*   How does what's run in this notebook compare to the full versions of Alphafold/Openfold?\n",
        "    *   This Colab version of OpenFold searches a selected portion of the BFD dataset and currently doesn’t use templates, so its accuracy is reduced in comparison to the full version, which is analogous to what's described in the [AlphaFold paper](https://doi.org/10.1038/s41586-021-03819-2) and [Github repo](https://github.com/deepmind/alphafold/). The full version of OpenFold can be run from our own [GitHub repo](https://github.com/aqlaboratory/openfold).\n",
        "*   What is a Colab?\n",
        "    *   See the [Colab FAQ](https://research.google.com/colaboratory/faq.html).\n",
        "*   I received a warning “Notebook requires high RAM”, what do I do?\n",
        "    *   The resources allocated to your Colab vary. See the [Colab FAQ](https://research.google.com/colaboratory/faq.html) for more details.\n",
        "    *   You can execute the Colab nonetheless.\n",
        "*   I received an error “Colab CPU runtime not supported” or “No GPU/TPU found”, what do I do?\n",
        "    *   Colab CPU runtime is not supported. Try changing your runtime via _Runtime_ > _Change runtime type_ > _Hardware accelerator_ > _GPU_.\n",
        "    *   The type of GPU allocated to your Colab varies. See the [Colab FAQ](https://research.google.com/colaboratory/faq.html) for more details.\n",
        "    *   If you receive “Cannot connect to GPU backend”, you can try again later to see if Colab allocates you a GPU.\n",
        "    *   [Colab Pro](https://colab.research.google.com/signup) offers priority access to GPUs.\n",
        "*   Does this tool install anything on my computer?\n",
        "    *   No, everything happens in the cloud on Google Colab.\n",
        "    *   At the end of the Colab execution a zip-archive with the obtained prediction will be automatically downloaded to your computer.\n",
        "*   How should I share feedback and bug reports?\n",
        "    *   Please share any feedback and bug reports as an [issue](https://github.com/aqlaboratory/openfold/issues) on Github.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YfPhvYgKC81B"
      },
      "source": [
        "# License and Disclaimer\n",
        "\n",
        "This Colab notebook and other information provided is for theoretical modelling only, caution should be exercised in its use. It is provided ‘as-is’ without any warranty of any kind, whether expressed or implied. Information is not intended to be a substitute for professional medical advice, diagnosis, or treatment, and does not constitute medical or other professional advice.\n",
        "\n",
        "## AlphaFold/OpenFold Code License\n",
        "\n",
        "Copyright 2021 AlQuraishi Laboratory\n",
        "\n",
        "Copyright 2021 DeepMind Technologies Limited.\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\"); you may not use this file except in compliance with the License. You may obtain a copy of the License at https://www.apache.org/licenses/LICENSE-2.0.\n",
        "\n",
        "Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.\n",
        "\n",
        "## Model Parameters License\n",
        "\n",
        "DeepMind's AlphaFold parameters are made available under the terms of the Creative Commons Attribution 4.0 International (CC BY 4.0) license. You can find details at: https://creativecommons.org/licenses/by/4.0/legalcode\n",
        "\n",
        "\n",
        "## Third-party software\n",
        "\n",
        "Use of the third-party software, libraries or code referred to in this notebook may be governed by separate terms and conditions or license provisions. Your use of the third-party software, libraries or code is subject to any such terms and you should check that you can comply with any applicable restrictions or terms and conditions before use.\n",
        "\n",
        "\n",
        "## Mirrored Databases\n",
        "\n",
        "The following databases have been mirrored by DeepMind, and are available with reference to the following:\n",
        "* UniRef90: v2021\\_03 (unmodified), by The UniProt Consortium, available under a [Creative Commons Attribution-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nd/4.0/).\n",
        "* MGnify: v2019\\_05 (unmodified), by Mitchell AL et al., available free of all copyright restrictions and made fully and freely available for both non-commercial and commercial use under [CC0 1.0 Universal (CC0 1.0) Public Domain Dedication](https://creativecommons.org/publicdomain/zero/1.0/).\n",
        "* BFD: (modified), by Steinegger M. and Söding J., modified by DeepMind, available under a [Creative Commons Attribution-ShareAlike 4.0 International License](https://creativecommons.org/licenses/by/4.0/). See the Methods section of the [AlphaFold proteome paper](https://www.nature.com/articles/s41586-021-03828-1) for details."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "gpuType": "T4"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "nbformat": 4,
  "nbformat_minor": 0
}