constava_examples.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Bio2Byte/public_notebooks/blob/main/constava_examples.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CxFL0lgejUTC"
      },
      "source": [
        "# Constava - Examples for common use cases\n",
        "\n",
        "***\n",
        "***\n",
        "\n",
        "In this Jupyter notebook common use cases of the Constava software are described. Code examples are given for the Python interface, as well as their corresponding command line commands.\n",
        "\n",
        "\n",
        "## Table of contents\n",
        "\n",
        "***\n",
        "\n",
        "1. [Preparations](#preparations)\n",
        "2. [Basic usage](#basic-usage)\n",
        "    1. [Extracting backbone dihedrals from a conformationa ensemble](#extracting-backbone-dihedrals)\n",
        "    2. [Analyzing a conformational ensemble](#analyzing-a-conformational-ensemble)\n",
        "    3. [Generating and loading conformational state models](#generating-and-loading-csm)\n",
        "3. [Visualizing the results](#visualization)\n",
        "    1. [Plotting the results](#visualization-plotting)\n",
        "    2. [Mapping a property on the protein structure](##visualization-structure)\n",
        "4. [Specialized use cases](#specialized-use-cases)\n",
        "    1. [Using gmx chi output](#using-gmx-chi-output)\n",
        "    1. [Custom confromational state models](#custom-csm)\n",
        "    1. [Speeding up the analysis](#speeding-up-the-analysis)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GMzWakLMj7cU"
      },
      "source": [
        "<a name=\"preparations\"></a>\n",
        "\n",
        "## 1 Preparations\n",
        "\n",
        "***\n",
        "\n",
        "Execute this cell to get the example data and install necessary packages."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Df-Q94t3HVSx"
      },
      "outputs": [],
      "source": [
        "%%capture\n",
        "# Install constava and necessary dependencies\n",
        "!pip install constava matplotlib py3dmol secstructartist\n",
        "# Download and unzip test data\n",
        "!wget -O constava_examples_data.zip https://zenodo.org/record/8410767/files/constava_examples_data.zip?download=1\n",
        "!unzip -n constava_examples_data.zip"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZEyS5VsJk65f"
      },
      "source": [
        "<a name=\"basic-usage\"></a>\n",
        "\n",
        "## 2 Basic usage\n",
        "\n",
        "***\n",
        "\n",
        "The software provides two modes of interaction. Shell user may use the software from the command line, while users skilled in Python can import it as a module.\n",
        "In the following section the three basic functionalities of the software are described."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tCPWlOMFnYjv"
      },
      "source": [
        "<a name=\"extracting-backbone-dihedrals\"></a>\n",
        "\n",
        "### 2.1 Extracting backbone dihedrals from a conformationa ensemble\n",
        "\n",
        "To analyze a conformational ensemble, the user needs to provide the backbone dihedral angles $(\\phi, \\psi)$ extracted from the ensemble.\n",
        "\n",
        "Below an exemplary code for the extraction of $\\phi$- and $\\psi$-angles from an MD trajectory is given. The `calculate_dihedrals` function returns a `pandas.DataFrame` object with the backbone dihedral information along the trajectory. This information is then stored in a csv file.\n",
        "\n",
        "A warning will inform you, that no dihedral angles are extracted for the first and last residue in a protein. This is expected as for terminal residues only $\\phi$ or $\\psi$ is defined.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xn858E3xIw18"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "from constava.utils.dihedrals import calculate_dihedrals\n",
        "\n",
        "# Calculate dihedrals as a DataFrame\n",
        "dihedrals = calculate_dihedrals(structure=\"./1v66.pdb\", trajectory=\"./1v66.xtc\", degree=True)\n",
        "\n",
        "# Write dihedrals out as a csv\n",
        "dihedrals.to_csv(\"extracted_dihedrals.csv\", index=False, float_format=\"%.1f\")\n",
        "dihedrals.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TomcA2e4obBx"
      },
      "source": [
        "The same behavior can be achieved through the command line:\n",
        "\n",
        "```sh\n",
        "constava dihedrals -O \\\n",
        "    -s \"1v66.pdb\" \\\n",
        "    -f \"1v66.xtc\" \\\n",
        "    -o \"extracted_dihedrals.csv\" \\\n",
        "    --degrees \\\n",
        "    --precision 1\n",
        "```\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T26d_LbxouFp"
      },
      "source": [
        "<a name=\"analyzing-a-conformational-ensemble\"></a>\n",
        "\n",
        "### 2.2 Analyzing a conformational ensemble\n",
        "\n",
        "After the backbone dihedral angles for a conformational ensemble have been\n",
        "calculated, conformational state propensities for each residue can be inferred.\n",
        "\n",
        "The `Constava`-class provides an easy to use interface. Parameters can be\n",
        "altered any time using the `set_param()` method. Conformational state models\n",
        "that are initialized are stored inside the `Constava`-instance. That way, by\n",
        "changing the `input_files` and `output_file` attributes, the  calculation can\n",
        "be run on multiple ensembles without refitting the conformational state model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "BrXZLj6povTG"
      },
      "outputs": [],
      "source": [
        "from constava import Constava\n",
        "\n",
        "# Define input and output files\n",
        "PDBID = \"1v66\"\n",
        "input_file = f\"./{PDBID}_dihedrals.csv\"\n",
        "output_file = \"./constava_results.csv\"\n",
        "\n",
        "# Initialize Constava Python interface with parameters\n",
        "c = Constava(\n",
        "    input_files = input_file,\n",
        "    output_file = output_file,\n",
        "    input_degrees=True,\n",
        "    window = None,\n",
        "    bootstrap = [3,5,10,25],\n",
        "    verbose = 1)\n",
        "\n",
        "# Alter parameters after initialization\n",
        "c.set_param(\"window\", [1,3,5])\n",
        "print(c.show_params())\n",
        "\n",
        "# Run the calculation and write results\n",
        "c.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "twwyJ8Ssovro"
      },
      "source": [
        "An equivalent command from the command line would be:\n",
        "\n",
        "```sh\n",
        "# Run constava with debug-level output\n",
        "constava analyze \\\n",
        "    -i \"1v66_dihedrals.csv\" \\\n",
        "    -o \"constava_results.csv\" \\\n",
        "    --window 1 3 5 \\\n",
        "    --bootstrap 3 5 10 25 \\\n",
        "    --degrees \\\n",
        "    -v\n",
        "```\n",
        "\n",
        "Note, that in the command line each call fits their own conformational state model. To use the same model in multiple runs, first use the `constava fit-model` subcommand to save a model, and then load the model using the `--load-model` flag (see [below](#generating-and-loading-csm))."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HcWZZDfbqpal"
      },
      "source": [
        "<a name=\"generating-and-loading-csm\"></a>\n",
        "\n",
        "### 2.3 Generating and loading conformational state models\n",
        "\n",
        "Conformational state models are usually fitted at runtime. This generates a small but unnecessary computational overhead when running multiple analyses using the same model.\n",
        "\n",
        "To aleviate this issue, a model can be pre-trained and stored to be used in subsequent analyses directly.\n",
        "\n",
        "**Note:** Storing conformational state models is intended for quickly running multiple analyses using the same model. They are **not intended for long-term storage or sharing** your conformational state models.\n",
        "To store or share custom conformational state models, provide the training data as well as the model-fitting parameters."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true
        },
        "id": "wpxYCYPGsNsW"
      },
      "outputs": [],
      "source": [
        "import glob\n",
        "from constava import Constava\n",
        "\n",
        "# Fit a KDE-model on default dataset with bandwidth==0.13\n",
        "c = Constava(verbose = 2)\n",
        "csmodel = c.fit_csmodel(kde_bandwidth = .13)\n",
        "\n",
        "# Write the fitted model out as a pickle\n",
        "csmodel.dump_pickle(\"stored_model.pkl\")\n",
        "\n",
        "# Use the new model to analyze a confromational ensemble\n",
        "PDBID = \"1v66\"\n",
        "input_file = f\"./{PDBID}_dihedrals.csv\"\n",
        "output_file = \"./constava_results.csv\"\n",
        "\n",
        "c = Constava(\n",
        "    input_files = input_file,\n",
        "    output_file = output_file,\n",
        "    model_load = \"stored_model.pkl\",\n",
        "    input_degrees=True,\n",
        "    window = [1, 5, 10, 25],\n",
        "    verbose = 2)\n",
        "c.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Zxfud6VtsN9x"
      },
      "source": [
        "<a name=\"generating-and-loading-csm-shell\"></a>\n",
        "\n",
        "To fit, save and load a model from the command line use the following commands:\n",
        "\n",
        "```sh\n",
        "# Fit the model on default dataset\n",
        "constava fit-model -v \\\n",
        "    -o \"stored_model.pkl\" \\\n",
        "    --kde-bandwidth 0.13\n",
        "\n",
        "# Run a prediction using the saved model\n",
        "constava analyze \\\n",
        "    -i \"1v66_dihedrals.csv\" \\\n",
        "    -o \"1v66_constava.csv\" \\\n",
        "    --load-model \"stored_model.pkl\" \\\n",
        "    --window 1 5 10 25 \\\n",
        "    --degrees \\\n",
        "    -vv\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TSMwVCvr0J0a"
      },
      "source": [
        "<a name=\"visualization\"></a>\n",
        "\n",
        "## 3 Visualizing the results\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eVb6g79_yh6R"
      },
      "source": [
        "<a name=\"visualization-plotting\"></a>\n",
        "\n",
        "### 3.1 Plotting the results\n",
        "\n",
        "Below constava results are plotted for a simple protein structure with four helices."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true
        },
        "id": "k2kAZkNLzMK4"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import secstructartist as ssa\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "# Load the results csv\n",
        "states = [\"coreSheet\", \"surrSheet\", \"Other\", \"Turn\", \"surrHelix\", \"coreHelix\"]\n",
        "results = pd.read_csv(\"1v66_results.csv\")\n",
        "\n",
        "# Select all values, where method was 'window/5/'\n",
        "win5 = results.loc[results[\"#Method\"] == \"window/5/\"]\n",
        "resids = win5[\"ResIndex\"].to_numpy()\n",
        "vari = win5[\"Variability\"].to_numpy()\n",
        "prop = win5[states].to_numpy()\n",
        "secstruct = \"LLHHHHHHHHHHLLHHHHHHHHHHHLLLLLLLHHHHHHHHHHHHHHLLLHHHHHHHHHHHHHL\"\n",
        "\n",
        "# Generate a plot\n",
        "fig, (ax0, ax1, ax2) = plt.subplots(3,1, figsize=(4, 2.5), dpi=120, sharex=True,\n",
        "                                    gridspec_kw={\"hspace\": 0., \"height_ratios\": [.7,2,3]})\n",
        "\n",
        "# Plot secondary structure\n",
        "ssa.draw(secstruct, resids, ax=ax0)\n",
        "ax0.set_ylim([-.7, .55])\n",
        "ax0.set_axis_off()\n",
        "\n",
        "# Plot conformational state variability\n",
        "ax1.grid(True, \"major\", \"y\")\n",
        "ax1.plot(resids, vari)\n",
        "ax1.set_ylabel(\"Variablility\")\n",
        "ax1.set_yticks(np.arange(0, .8, .2))\n",
        "\n",
        "# Plot conformational state propensities as heatmap\n",
        "ax2.imshow(prop.T, vmin=0, vmax=1, cmap=\"Blues\", aspect=\"auto\",\n",
        "           extent=(np.min(resids)-.5, np.max(resids)+.5, len(states) -.5, -.5))\n",
        "ax2.set_yticks(np.arange(len(states)))\n",
        "ax2.set_yticklabels(states)\n",
        "ax2.set_xlabel(\"Residue index\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bf9Lx1_-yijP"
      },
      "source": [
        "<a name=\"visualization-structure\"></a>\n",
        "\n",
        "### 3.2 Mapping a property on the protein structure\n",
        "\n",
        "Analogous to plotting the data, the constava results can also be mapped on a protein structure, using the B-factor column of the PDB file format. Since most values calculated by constava are in the range `[0, 1], it is useful to scale these values up before writing them to a PDB file to use the whole range of the B-factor column."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true
        },
        "id": "rHQFjewM_WoC"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from matplotlib import cm\n",
        "import py3Dmol\n",
        "\n",
        "def constava2pdb(input_pdb, output_pdb, df, state_label, factor = 1):\n",
        "    \"\"\"Write values scaled by `factor` to temp factor column in PDB\"\"\"\n",
        "    with open(input_pdb, \"r\") as ipdb, open(output_pdb, \"w\") as opdb:\n",
        "        for ln in ipdb:\n",
        "            try:\n",
        "                assert ln.startswith(\"ATOM  \")\n",
        "                resn, resi = ln[17:20], int(ln[22:26])\n",
        "                bfac = df.loc[(df[\"ResName\"] == resn) & (df[\"ResIndex\"] == resi), [state_label]].iloc[0,0]\n",
        "            except (AssertionError, IndexError):\n",
        "                opdb.write(ln)\n",
        "            else:\n",
        "                opdb.write(\"{0}{2:6.2f}{1}\".format(ln[:60], ln[66:], factor*bfac))\n",
        "\n",
        "# Try out other properties: coreHelix, surrHelix, coreSheet, surrSheet, Turn, Other, Variability\n",
        "PROPERTY = \"Turn\"\n",
        "\n",
        "# Write \"Conformational State Variability\" values to the pdb file\n",
        "results = pd.read_csv(\"1v66_results.csv\")\n",
        "bstrp5 = results.loc[results[\"#Method\"] == \"bootstrap/5/500//\"]\n",
        "constava2pdb(\"1v66.pdb\", \"1v66_bfactor.pdb\", bstrp5, \"coreHelix\", factor=100)\n",
        "\n",
        "# Make a custom linear gradient\n",
        "mygradient = []\n",
        "for x in np.linspace(.1, .9, 50):\n",
        "    r,g,b,_ = cm.viridis(x)\n",
        "    mygradient.append(f\"#{int(255*r):x}{int(255*g):x}{int(255*b):x}\")\n",
        "\n",
        "# Display protein colored by factor\n",
        "with open(\"1v66_bfactor.pdb\", \"r\") as fhandle:\n",
        "    struct = fhandle.read()\n",
        "\n",
        "view = py3Dmol.view(width=500, height=500)\n",
        "view.setBackgroundColor('black')\n",
        "view.setViewStyle({'style':'outline','color':'#888','width':0.1})\n",
        "view.addModel(struct, \"pdb\")\n",
        "view.setStyle({\"model\": 0}, {'cartoon': {'colorscheme': {\n",
        "    'prop':'b', 'gradient': \"linear\", \"colors\": mygradient, 'min':0,'max':70}}})\n",
        "view.addSurface(py3Dmol.VDW, {\"opacity\": .5, 'colorscheme':{\n",
        "    'prop':'b', 'gradient':\"linear\", \"colors\": mygradient, 'min':0, 'max':70}}, {\"model\": 0})\n",
        "view.zoomTo()\n",
        "view.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lkD3ducwqtc-"
      },
      "source": [
        "<a name=\"specialized-use-cases\"></a>\n",
        "\n",
        "## 4 Specialized use cases\n",
        "\n",
        "***\n",
        "\n",
        "In the following section we will discuss special use cases, that likely only apply to a small number of users."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "imVHn4a9seMs"
      },
      "source": [
        "<a name=\"using-gmx-chi-output\"></a>\n",
        "\n",
        "### 4.1 Using `gmx chi` output\n",
        "\n",
        "Users of the [GROMACS](https://manual.gromacs.org/current/index.html) simulation package may use GROMACS' `gmx chi` module to extract the backbone dihedral angles:\n",
        "\n",
        "```sh\n",
        "gmx chi -s <struct.gro> -f <trajectory.xtc> -rama\n",
        "```\n",
        "\n",
        "This creates a couple of files named `ramaPhiPsi<RESN><RESID>.xvg` with `<RESN>`being the residue's code and `<RESID>` being the residues index. Constava can read those files directly.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true,
          "base_uri": "https://localhost:8080/"
        },
        "id": "Gx5Bkg5AtHIu",
        "outputId": "e34b0fdf-0f1a-4958-8ca6-8d1de6666c02"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:19] Constava: Initializing python interface...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Constava: Initializing python interface...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:19] Setting `window = [1, 3, 5]`\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Setting `window = [1, 3, 5]`\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "ConstavaParameters(input_files=['./gmxchi/ramaPhiPsiLEU23.xvg', './gmxchi/ramaPhiPsiSER50.xvg', './gmxchi/ramaPhiPsiLYS30.xvg', './gmxchi/ramaPhiPsiGLN54.xvg', './gmxchi/ramaPhiPsiARG33.xvg', './gmxchi/ramaPhiPsiTYR25.xvg', './gmxchi/ramaPhiPsiHIS43.xvg', './gmxchi/ramaPhiPsiLEU45.xvg', './gmxchi/ramaPhiPsiGLU6.xvg', './gmxchi/ramaPhiPsiLYS34.xvg', './gmxchi/ramaPhiPsiPRO51.xvg', './gmxchi/ramaPhiPsiALA41.xvg', './gmxchi/ramaPhiPsiARG15.xvg', './gmxchi/ramaPhiPsiVAL11.xvg', './gmxchi/ramaPhiPsiLYS56.xvg', './gmxchi/ramaPhiPsiSER13.xvg', './gmxchi/ramaPhiPsiGLY48.xvg', './gmxchi/ramaPhiPsiALA26.xvg', './gmxchi/ramaPhiPsiLYS40.xvg', './gmxchi/ramaPhiPsiMET10.xvg', './gmxchi/ramaPhiPsiARG63.xvg', './gmxchi/ramaPhiPsiGLU59.xvg', './gmxchi/ramaPhiPsiGLN20.xvg', './gmxchi/ramaPhiPsiMET1.xvg', './gmxchi/ramaPhiPsiGLY24.xvg', './gmxchi/ramaPhiPsiLEU14.xvg', './gmxchi/ramaPhiPsiILE57.xvg', './gmxchi/ramaPhiPsiTHR39.xvg', './gmxchi/ramaPhiPsiHIS35.xvg', './gmxchi/ramaPhiPsiVAL21.xvg', './gmxchi/ramaPhiPsiLYS46.xvg', './gmxchi/ramaPhiPsiGLN9.xvg', './gmxchi/ramaPhiPsiASN29.xvg', './gmxchi/ramaPhiPsiSER4.xvg', './gmxchi/ramaPhiPsiLEU42.xvg', './gmxchi/ramaPhiPsiLYS58.xvg', './gmxchi/ramaPhiPsiALA2.xvg', './gmxchi/ramaPhiPsiLEU22.xvg', './gmxchi/ramaPhiPsiLEU38.xvg', './gmxchi/ramaPhiPsiVAL53.xvg', './gmxchi/ramaPhiPsiLEU44.xvg', './gmxchi/ramaPhiPsiARG64.xvg', './gmxchi/ramaPhiPsiLEU19.xvg', './gmxchi/ramaPhiPsiLEU7.xvg', './gmxchi/ramaPhiPsiGLU36.xvg', './gmxchi/ramaPhiPsiALA52.xvg', './gmxchi/ramaPhiPsiGLY32.xvg', './gmxchi/ramaPhiPsiALA5.xvg', './gmxchi/ramaPhiPsiVAL16.xvg', './gmxchi/ramaPhiPsiLYS8.xvg', './gmxchi/ramaPhiPsiCYS49.xvg', './gmxchi/ramaPhiPsiGLU18.xvg', './gmxchi/ramaPhiPsiALA47.xvg', './gmxchi/ramaPhiPsiARG62.xvg', './gmxchi/ramaPhiPsiMET12.xvg', './gmxchi/ramaPhiPsiHIS31.xvg', './gmxchi/ramaPhiPsiASP3.xvg', './gmxchi/ramaPhiPsiLEU37.xvg', './gmxchi/ramaPhiPsiPHE65.xvg', './gmxchi/ramaPhiPsiGLY27.xvg', './gmxchi/ramaPhiPsiSER17.xvg', './gmxchi/ramaPhiPsiTYR61.xvg', './gmxchi/ramaPhiPsiMET55.xvg', './gmxchi/ramaPhiPsiLEU60.xvg', './gmxchi/ramaPhiPsiARG28.xvg'], input_format='auto', output_file='./constava_results.csv', output_format='auto', model_type='kde', model_load=None, model_data=None, model_dump=None, window=[1, 3, 5], bootstrap=[3, 5, 10, 25], window_series=[], bootstrap_series=[], bootstrap_samples=500, input_degrees=True, model_data_degrees=False, precision=4, kde_bandwidth=0.13, grid_points=10000, seed=None, verbose=1)\n",
            "[2023-10-05 14:52:19] Initializing reader for input file(s)...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Initializing reader for input file(s)...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:19] Initializing writer for results...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Initializing writer for results...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:19] Fitting model to data in: /usr/local/lib/python3.10/dist-packages/constava/data/constava_csdata.json\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Fitting model to data in: /usr/local/lib/python3.10/dist-packages/constava/data/constava_csdata.json\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... model fitted: <kdeModel states=('coreHelix', 'surrHelix', 'coreSheet', 'surrSheet', 'Turn', 'Other')>\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... model fitted: <kdeModel states=('coreHelix', 'surrHelix', 'coreSheet', 'surrSheet', 'Turn', 'Other')>\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] Initializing calculator with <kdeModel states=('coreHelix', 'surrHelix', 'coreSheet', 'surrSheet', 'Turn', 'Other')>...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Initializing calculator with <kdeModel states=('coreHelix', 'surrHelix', 'coreSheet', 'surrSheet', 'Turn', 'Other')>...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: window/1/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/1/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: window/3/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/3/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: window/5/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/5/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: bootstrap/3/500//\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: bootstrap/3/500//\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: bootstrap/5/500//\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: bootstrap/5/500//\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: bootstrap/10/500//\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: bootstrap/10/500//\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] ... adding subsampling method: bootstrap/25/500//\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: bootstrap/25/500//\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] Reading dihedrals from 65 files...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Reading dihedrals from 65 files...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:52:20] Starting inference...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Starting inference...\n",
            "100%|██████████| 65/65 [01:52<00:00,  1.73s/residues]"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Writing results to file: ./constava_results.csv\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "\n",
            "INFO:Constava:Writing results to file: ./constava_results.csv\n"
          ]
        }
      ],
      "source": [
        "import glob\n",
        "from constava import Constava\n",
        "\n",
        "# Define input and output files\n",
        "PDBID = \"1v66\"\n",
        "input_files = glob.glob(f\"./gmxchi/ramaPhiPsi*.xvg\")\n",
        "output_file = \"./constava_results.csv\"\n",
        "\n",
        "# Initialize Constava Python interface with parameters\n",
        "c = Constava(\n",
        "    input_files = input_files,\n",
        "    output_file = output_file,\n",
        "    input_degrees=True,\n",
        "    window = None,\n",
        "    bootstrap = [3,5,10,25],\n",
        "    verbose = 1)\n",
        "\n",
        "# Alter parameters after initialization\n",
        "c.set_param(\"window\", [1,3,5])\n",
        "print(c.show_params())\n",
        "\n",
        "# Run the calculation and write results\n",
        "c.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7twIOhsX4CTr"
      },
      "source": [
        "An equivalent command from the command line would be:\n",
        "\n",
        "```sh\n",
        "# Run constava\n",
        "PDBID=\"1v66\"\n",
        "constava analyze \\\n",
        "    -i ${PDBID}/ramaPhiPsi*.xvg \\\n",
        "    -o constava_results.csv \\\n",
        "    --window 1 3 5 \\\n",
        "    --bootstrap 3 5 10 25 \\\n",
        "    --degrees \\\n",
        "    -v\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G8V-YulcseCN"
      },
      "source": [
        "<a name=\"custom-csm\"></a>\n",
        "\n",
        "### 4.2 Custom confromational state models\n",
        "\n",
        "Besides the default conformational state model from the associated publication,\n",
        "users can also define their own conformational states. This is done by providing\n",
        "a JSON file with $(\\phi, \\psi)$ dihedral pairs that define each conformational state. An example would look like this:\n",
        "\n",
        "> ```python\n",
        "> {\n",
        ">     \"MyHelixState\": [\n",
        ">         [phi1, psi1],\n",
        ">         [phi2, psi2],\n",
        ">         ... ],\n",
        ">     \"MyOtherState\": [\n",
        ">         [phi1, psi1],\n",
        ">         [phi2, psi2],\n",
        ">         ... ],\n",
        ">     ...\n",
        "> }\n",
        "> ```\n",
        "\n",
        "Here, `phi1`, `psi1`, etc. are floats representing the dihedral in radians or degrees (use of --degrees flag required, see below).\n",
        "\n",
        "How to fit and use the probabilistic conformational state model from such a\n",
        "JSON file is described in the code block below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "background_save": true
        },
        "id": "a5bWtei7EMox",
        "outputId": "15848ec7-2f2f-4242-b1e7-9201ef8cdb1a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Constava: Initializing python interface...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Constava: Initializing python interface...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Fitting model to data in: custom_confstates.json\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Fitting model to data in: custom_confstates.json\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] ... model fitted: <kdeModel states=('helix', 'sheet', 'coil', 'N/A')>\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... model fitted: <kdeModel states=('helix', 'sheet', 'coil', 'N/A')>\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Constava: Initializing python interface...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Constava: Initializing python interface...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Initializing reader for input file(s)...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Initializing reader for input file(s)...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Initializing writer for results...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Initializing writer for results...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Loading conformational state models from file: custom_model.pkl\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Loading conformational state models from file: custom_model.pkl\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] ... model loaded: <kdeModel states=('helix', 'sheet', 'coil', 'N/A')>\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... model loaded: <kdeModel states=('helix', 'sheet', 'coil', 'N/A')>\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Initializing calculator with <kdeModel states=('helix', 'sheet', 'coil', 'N/A')>...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Initializing calculator with <kdeModel states=('helix', 'sheet', 'coil', 'N/A')>...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] ... adding subsampling method: window/1/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/1/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] ... adding subsampling method: window/5/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/5/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] ... adding subsampling method: window/10/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/10/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] ... adding subsampling method: window/25/\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:... adding subsampling method: window/25/\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Reading dihedrals from 1 files...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Reading dihedrals from 1 files...\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:13] Starting inference...\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:Constava:Starting inference...\n",
            "100%|██████████| 63/63 [00:17<00:00,  3.61residues/s]"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[2023-10-05 14:54:30] Writing results to file: constava_results.csv\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "\n",
            "INFO:Constava:Writing results to file: constava_results.csv\n"
          ]
        }
      ],
      "source": [
        "from constava import Constava\n",
        "\n",
        "# Fit the grid-interpolation model\n",
        "c = Constava(verbose = 1)\n",
        "csmodel = c.fit_csmodel(model_type = \"kde\",\n",
        "                        model_data = \"custom_confstates.json\",\n",
        "                        model_data_degrees = True,\n",
        "                        kde_bandwidth = .13)\n",
        "\n",
        "# Write the fitted model out as a pickle\n",
        "csmodel.dump_pickle(\"custom_model.pkl\")\n",
        "\n",
        "# Use the new model to analyze a confromational ensemble\n",
        "PDBID = \"1v66\"\n",
        "input_files = glob.glob(f\"./{PDBID}_dihedrals.csv\")\n",
        "output_file = \"constava_results.csv\"\n",
        "c = Constava(\n",
        "    input_files = input_files,\n",
        "    output_file = output_file,\n",
        "    model_load = \"custom_model.pkl\",\n",
        "    input_degrees=True,\n",
        "    window = [1, 5, 10, 25],\n",
        "    verbose = 1)\n",
        "c.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tV0u4DuSENOK"
      },
      "source": [
        "The same actions can be taken from the command line running:\n",
        "\n",
        "```sh\n",
        "# Fit the model on custom dataset\n",
        "constava fit-model -v \\\n",
        "    -i \"custom_confstates.json\" \\\n",
        "    -o \"custom_model.pkl\" \\\n",
        "    --kde-bandwidth 0.13 \\\n",
        "    --degrees\n",
        "\n",
        "# Run a prediction using the saved model\n",
        "constava analyze -v \\\n",
        "    -i \"1v66_dihedrals.csv\" \\\n",
        "    -o \"1v66_constava.csv\" \\\n",
        "    --load-model \"custom_model.pkl\" \\\n",
        "    --window 1 5 10 25 \\\n",
        "    --degrees\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kBQGAby4tHxW"
      },
      "source": [
        "<a name=\"speeding-up-the-analysis\"></a>\n",
        "\n",
        "### 4.3 Speeding up the analysis\n",
        "\n",
        "Analyzing large conformational ensembles with constava can be time consuming. The bottleneck here is the inference of propensities for each conformational state from their respective probability density functions represented by Gaussian kernel density estimators (KDEs).\n",
        "\n",
        "To speed up calculations, users can use the `grid` model type. This replaces the KDEs by a grid representing fixed points on the probability density function. The propensities are consequently inferred using linear interpolation between the closest grid points.\n",
        "\n",
        "While this is less accurate than direct inference from the KDEs, using the `grid` model type has the potential to speed up calculations significantly. Note, that fitting a grid model may take some time since the grid points need first to be inferred from a KDE model.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "w14XdEkctQJH"
      },
      "outputs": [],
      "source": [
        "from constava import Constava\n",
        "\n",
        "# Fit the grid-interpolation model\n",
        "c = Constava(verbose = 1)\n",
        "csmodel = c.fit_csmodel(model_type = \"grid\",\n",
        "                        kde_bandwidth = .13,\n",
        "                        grid_points = 6400)\n",
        "\n",
        "# Write the fitted model out as a pickle\n",
        "csmodel.dump_pickle(\"grid_model.pkl\")\n",
        "\n",
        "# Use the new model to analyze a confromational ensemble\n",
        "PDBID = \"1v66\"\n",
        "input_files = glob.glob(f\"./{PDBID}_dihedrals.csv\")\n",
        "output_file = \"constava_results.csv\"\n",
        "c = Constava(\n",
        "    input_files = input_files,\n",
        "    output_file = output_file,\n",
        "    model_load = \"grid_model.pkl\",\n",
        "    input_degrees=True,\n",
        "    window = [1, 5, 10, 25],\n",
        "    verbose = 1)\n",
        "c.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6qhQ1IzXtZIB"
      },
      "source": [
        "To fit, save and load a grid-interpolation model from the command line use the\n",
        "following commands:\n",
        "\n",
        "```sh\n",
        "# Fit the grid-interpolation model from default dataset\n",
        "constava fit-model -v \\\n",
        "    -o grid_model.pkl \\\n",
        "    --model-type grid \\\n",
        "    --kde-bandwidth 0.13 \\\n",
        "    --grid-points 6400\n",
        "\n",
        "# Run a prediction using the saved model\n",
        "constava analyze -v \\\n",
        "    -i \"1v66_dihedrals.csv\" \\\n",
        "    -o \"constava_results.csv\" \\\n",
        "    --load-model grid_model.pkl \\\n",
        "    --window 1 5 10 25 \\\n",
        "    --degrees\n",
        "```"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}