Untitled8.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyOxuHSS04zcwZDWI+BUJVKa",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "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/Imanjdd/protein_aa_analyze/blob/main/Untitled8.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Bio2Byte - Multiple Sequence Alignment Analysis - Overall behavior.\n",
        "\n",
        "#### **Goal: Compute and visualize the biophysical and sequence conservation of a Multiple Sequence Alignment (MSA) using the b2btools.**\n",
        "\n",
        "*In this notebook we compute the biophysical and sequence conservation of the OVERALL BEHAVIOR of the MSA. To compute them in function of a reference protein, please refer to the following [notebook](https://colab.research.google.com/drive/1BBAGn4PCL9cqYaED5SRe4V8wXoeBibVq#scrollTo=guB0TE1iXZtD).*\n",
        "\n",
        "The following code will:\n",
        "0. Install the b2btools and its dependencies.\n",
        "1. Upload the Multiple Sequence Alignment (MSA) to study.\n",
        "2. Predict the biophysical properties of the input MSA:\n",
        "      *   Backbone and sidechain dynamics (Dynamine)\n",
        "      *   Conformational propensities (sheet, helix, coil, polyproline II) (Dynamine)\n",
        "      *   Early folding propensities (EFoldMine)\n",
        "      *   Disorder propensities (Disomine)\n",
        "3. Compute the sequence conservation (Shannon's entropy) of the input MSA at every postion as well as the MSA occupency.\n",
        "4. Compute 2D plots of the biophysical behaviour, the conservation of the biophysical properties and the sequence conservation.\n",
        "5. Compute the Gaussian Mixture Model scores of the proteins to identify residues with a different behavior compared to the other residues at the same position in the MSA.\n",
        "5. Download your results.\n",
        "\n",
        "**Check out our webserver: [online b2BTools](https://bio2byte.be/b2btools/)**"
      ],
      "metadata": {
        "id": "-i-Lq6IV3JuM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# 1.a Install the Pomegranate Package"
      ],
      "metadata": {
        "id": "1Bp1DcakgHf7"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "cellView": "form",
        "id": "NawWwzLIxdpf"
      },
      "outputs": [],
      "source": [
        "#@title First step : install the cython and scipy dependencies\n",
        "%%capture\n",
        "!pip install \"cython~=0.29.14\" \"scipy>=0.17.0\""
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Second step : Search the pomegranate package on Github\n",
        "!git clone -b v0.14.9 --single-branch https://github.com/jmschrei/pomegranate.git"
      ],
      "metadata": {
        "cellView": "form",
        "id": "JjnksCUjxgta"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Last step : Install the pomegranate package\n",
        "%%capture\n",
        "#@markdown Please be patient, this can take several minutes...\n",
        "%cd pomegranate\n",
        "!python setup.py install"
      ],
      "metadata": {
        "cellView": "form",
        "id": "idJE0vX0xjMS"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 1.b Install the B2BTools Package\n",
        "%%capture\n",
        "!pip install b2bTools==3.0.7 ipympl hmmer\n",
        "#@markdown Please be patient, this can take several minutes...\n",
        "\n",
        "#@markdown ⚠️ Once installed, please restart the session and run step 1.c directly, without rerunning the first step."
      ],
      "metadata": {
        "cellView": "form",
        "id": "TRzH-otixmc2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 1.c Install the dependencies\n",
        "%%capture\n",
        "import math\n",
        "import Bio\n",
        "from Bio import AlignIO\n",
        "import os\n",
        "import json\n",
        "import shutil\n",
        "from google.colab import files\n",
        "from b2bTools.multipleSeq.Predictor import MineSuiteMSA\n",
        "from google.colab import output\n",
        "import numpy as np\n",
        "from sklearn import mixture\n",
        "import sys\n",
        "from os import path"
      ],
      "metadata": {
        "cellView": "form",
        "id": "onUnDIryxp6R"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 1.d Upload your MSA file\n",
        "\n",
        "#@markdown Execute this cell to upload your MSA file from your computer.\n",
        "\n",
        "#@markdown **🚨⚠️ IMPORTANT input file size limits:**\n",
        "#@markdown * Dynamine (dynamics and secondary structure): min 5 residues per sequence\n",
        "#@markdown * Disomine (disorder): min 5 residues per sequence\n",
        "#@markdown * EfoldMine (early folding): min 5 residues per sequence\n",
        "\n",
        "#@markdown **🚨⚠️ IMPORTANT input file format:**\n",
        "#@markdown The file can be in CLUSTAL, FASTA, BaliBase, PSI, A3M, Blast, PHYLIP and STOCKHOLM format\n",
        "\n",
        "msa_filename = None\n",
        "\n",
        "uploaded = files.upload()\n",
        "\n",
        "for fn in uploaded.keys():\n",
        "  print('MSA file \"{name}\" with length {length} bytes uploaded with success'.format(\n",
        "      name=fn, length=len(uploaded[fn])))\n",
        "  msa_filename = fn\n",
        "\n",
        "\n",
        "#@markdown Select the format of your MSA file:\n",
        "fasta = True #@param {type:\"boolean\"}\n",
        "clustal = False #@param {type:\"boolean\"}\n",
        "phylip = False #@param {type:\"boolean\"}\n",
        "stockholm = False #@param {type:\"boolean\"}\n",
        "\n",
        "# BaliBase = False #@param {type:\"boolean\"}\n",
        "# PSI = False #@param {type:\"boolean\"}\n",
        "# A3M = False #@param {type:\"boolean\"}\n",
        "# Blast = False #@param {type:\"boolean\"}\n",
        "\n",
        "extensions=[fasta,clustal,phylip,stockholm]\n",
        "extensions_names=[\"fasta\",\"clustal\",\"phylip\",\"stockholm\"]\n",
        "for ext in extensions:\n",
        "  if ext:\n",
        "    extension = extensions_names[extensions.index(ext)]"
      ],
      "metadata": {
        "cellView": "form",
        "id": "Q2MD5Ei0xuMs"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 2. Select the predictors you want to use and run them\n",
        "#@markdown Select the bio2Byte tools you want to include into the predictions:\n",
        "%%capture\n",
        "\n",
        "#@markdown ### DynaMine predictor tool\n",
        "#@markdown >Fast predictor of protein backbone dynamics using only sequence information as input.\n",
        "#@markdown >The version here also predicts side-chain dynamics and secondary structure predictors\n",
        "#@markdown >using the same principle.\n",
        "\n",
        "#@markdown **Prediction values included**: `backbone`, `sidechain`, `helix`, `ppII`, `coil`, and `sheet`\n",
        "\n",
        "DynaMine = True #@param {type:\"boolean\"}\n",
        "\n",
        "#@markdown ### DisoMine predictor tool\n",
        "#@markdown >Predicts protein disorder with recurrent neural networks not directly\n",
        "#@markdown >from the amino acid sequence, but instead from more generic predictions of key\n",
        "#@markdown >biophysical properties, here protein dynamics, secondary structure\n",
        "#@markdown >and early folding.\n",
        "\n",
        "#@markdown **Prediction values included**: `disoMine`\n",
        "\n",
        "DisoMine = True #@param {type:\"boolean\"}\n",
        "\n",
        "#@markdown ### EFoldMine predictor tool\n",
        "#@markdown >Predicts from the primary amino acid sequence of a protein,\n",
        "#@markdown >which amino acids are likely involved in early folding events.\n",
        "\n",
        "#@markdown **Prediction values included**: `earlyFolding`\n",
        "\n",
        "EFoldMine = True #@param {type:\"boolean\"}\n",
        "\n",
        "#@markdown **Don't forget to run the cell after ticking the boxes.**\n",
        "\n",
        "msaSuite = MineSuiteMSA()\n",
        "\n",
        "if DynaMine and not DisoMine and not EFoldMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('dynamine'))\n",
        "elif DisoMine and not DynaMine and not EFoldMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('disoMine'))\n",
        "elif EFoldMine and not DynaMine and not DisoMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('eFoldMine'))\n",
        "\n",
        "elif DynaMine and DisoMine and not EFoldMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('disoMine', 'dynamine'))\n",
        "elif DynaMine and EFoldMine and not DisoMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('eFoldMine','dynamine'))\n",
        "elif DisoMine and EFoldMine and not DynaMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('eFoldMine', 'disoMine'))\n",
        "\n",
        "elif DisoMine and EFoldMine and DynaMine:\n",
        "  msaSuite.predictAndMapSeqsFromMSA(f\"/content/{msa_filename}\", predTypes = ('eFoldMine', 'disoMine', 'dynamine'))\n",
        "\n",
        "msaSuite.getDistributions()\n",
        "jsondata_list = [msaSuite.alignedPredictionDistribs]"
      ],
      "metadata": {
        "cellView": "form",
        "id": "MlNNYX-r2uae"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 2.b Select a protein from the MSA to study its behavior on top of the behavior of the MSA\n",
        "\n",
        "#@markdown Select the protein you want to include in the study in the dropdown list that appears once you execute this cell.\n",
        "\n",
        "#@markdown In case you don't want to study a protein in particular, select None.\n",
        "\n",
        "\n",
        "import Bio\n",
        "from Bio import AlignIO\n",
        "import ipywidgets as widgets\n",
        "\n",
        "alignment_file=AlignIO.read( msa_filename, extension)\n",
        "\n",
        "ids =[]\n",
        "for prot in alignment_file:\n",
        "  ids.append(prot.id)\n",
        "\n",
        "out = widgets.Dropdown(\n",
        "    options=[\"None\"]+ids,\n",
        "    description='Select protein:',\n",
        "    disabled=False,\n",
        ")\n",
        "display(out)"
      ],
      "metadata": {
        "cellView": "form",
        "id": "Gf8NiZz1zNLn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 2.c. Set the variables in function of the value selected above\n",
        "\n",
        "#selected protein in drop down list\n",
        "selected_prot = []\n",
        "\n",
        "try:\n",
        "  if out.value != \"None\":\n",
        "\n",
        "    selected_prot.append(out.value)\n",
        "\n",
        "    jsondata_list_selected = []\n",
        "    print(f\"The selected protein is {out.value}\")\n",
        "\n",
        "    #prediction selected protein\n",
        "    for prot in selected_prot:\n",
        "      predictions_single_seq = msaSuite.allAlignedPredictions\n",
        "      jsondata_list_selected.append([predictions_single_seq[prot]])\n",
        "except:\n",
        "  pass"
      ],
      "metadata": {
        "cellView": "form",
        "id": "qmJGrA5hzt80"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 4.a. Sequence conservation (Shannon's entropy)\n",
        "#@markdown Once this cell has been executed, the conservation $C$ of the residues\n",
        "#@markdown at a particular column $x$ within the MSA is computed by using\n",
        "#@markdown the Shannon's entropy $E$ and by taking into account the number of gaps $G$ in that column:\n",
        "#@markdown $$ C(x) = (1- E(x))(1-G(x)) $$\n",
        "\n",
        "#@markdown Shannon's entropy at column $x$ equals:\n",
        "\n",
        "#@markdown $$ E(x) = - \\lambda \\sum_{a}^{K} p_{a} log_{2}(p_{a}) $$\n",
        "\n",
        "#@markdown where $K$, the alphabet size which equals 21 as it includes the 20 amino acid types and 1 symbol for the gaps,\n",
        "\n",
        "#@markdown and $ p_{a}$ is the probability of observing the *ath* symbol type at position $x$.\n",
        "\n",
        "#@markdown $\\lambda$ scales the entropy to range [0,1]:\n",
        "\n",
        "#@markdown if the number of residue types in column $x$, $M$ ($1\\leq M \\leq K$) equals 1:\n",
        "#@markdown $$ \\lambda = 1 $$\n",
        "#@markdown else:\n",
        "#@markdown $$ \\lambda = [log_{2}(M)]^{-1} $$\n",
        "\n",
        "#@markdown Moreover, the number of gaps in column $x$:\n",
        "#@markdown $$ G(x) = \\dfrac{n}{N} $$\n",
        "#@markdown where $n$ is the number of gaps in column $x$.\n",
        "\n",
        "#Shannon_entropy calculates the conservation of the amino acids at each position. It takes into account\n",
        "#the gaps in the columns. Check https://doi.org/10.1002/prot.10146 for info formula\n",
        "def shannon_entropy(list_input):\n",
        "    tot = len(list_input) #total number of AA at particular position in MSA\n",
        "    gaps = list_input.count(\"-\") / tot #count frequency of gaps in that column of MSA\n",
        "    unique_base = set(list_input) #remove duplicates, \"-\" is seen as an amino acid type\n",
        "    unique_base_len = len(unique_base) #total number of AA at particular position in MSA with no duplicates\n",
        "\n",
        "\n",
        "    entropy_list = [] # entropy of AA at particular position\n",
        "\n",
        "    for base in unique_base:\n",
        "        n_i = list_input.count(base)\n",
        "        P_i = n_i / tot\n",
        "\n",
        "        entropy_i = P_i * (math.log(P_i, 2))\n",
        "        entropy_list.append(entropy_i)\n",
        "\n",
        "    #sum entropy of every residue at a position and normalize it so it is between 0 and 1\n",
        "    entropy_sum = math.fsum(entropy_list)\n",
        "\n",
        "    if unique_base_len == 1: # log(1, 2) = 0; n/0 throws ZeroDivisionError\n",
        "        shannon_entropy = entropy_sum\n",
        "    else:\n",
        "        unique_base_len_log = math.log(unique_base_len, 2)\n",
        "        shannon_entropy = (-1 / unique_base_len_log) * entropy_sum\n",
        "\n",
        "    #Return entropy AA and entropy gaps at 1 position\n",
        "    #If entropy is high than there are many possible arrangements (high variability)\n",
        "    return shannon_entropy, gaps\n",
        "\n",
        "def conservation(alignment_file):\n",
        "    conservation_AA_list = []\n",
        "\n",
        "    for col_no in range(len(list(alignment_file[0]))):\n",
        "        list_input = list(alignment_file[:, col_no])\n",
        "\n",
        "        sh_entropy_AA, sh_entropy_gaps = shannon_entropy(list_input)\n",
        "\n",
        "        #Translate entropy into conservation. 0 means no conservation, 1 highly conserved\n",
        "        #As e took into account the gaps: if only gaps: entropy=0 freq_gaps=1 conservation = 0\n",
        "        conservation_AA = (1 - sh_entropy_AA) * (1 - sh_entropy_gaps)\n",
        "        conservation_AA_list.append(conservation_AA)\n",
        "\n",
        "    return conservation_AA_list\n",
        "\n",
        "#Read MSA\n",
        "alignment_file = AlignIO.read(msa_filename, extension)\n",
        "conservation_AA_list = conservation(alignment_file)"
      ],
      "metadata": {
        "cellView": "form",
        "id": "vmm_hJp_zuv2"
      },
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 4.b. Sequence occupency in MSA\n",
        "#@markdown Once this cell has been executed, the sequence occupency in the MSA is computed (i.e., the frequency of gaps at every residue position in the MSA)\n",
        "\n",
        "counter = [0] * alignment_file.get_alignment_length()\n",
        "total_seq = len(alignment_file)\n",
        "for prot in alignment_file:\n",
        "  selected_prot_seq=list(prot.seq)\n",
        "  for position,residue in enumerate(selected_prot_seq):\n",
        "      if residue != \"-\":\n",
        "          counter[position] +=1\n",
        "freq_gap = [i/total_seq for i in counter]"
      ],
      "metadata": {
        "cellView": "form",
        "id": "liZRpP6Bz5QP"
      },
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 5. Prepare to plot results\n",
        "#@markdown Run this cell in order to prepare the notebook context to render different plots. Please be patient this can take a few minutes.\n",
        "%%capture\n",
        "%matplotlib widget\n",
        "\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import os\n",
        "import math\n",
        "\n",
        "if not os.path.exists(\"/content/results\"):\n",
        "    os.mkdir(\"/content/results\", )\n",
        "\n",
        "if DynaMine and not DisoMine and not EFoldMine:\n",
        "    NB_SUBPLOTS = 6\n",
        "    PREDICTION_TITLES = {\n",
        "      'backbone': \"DynaMine backbone dynamics\",\n",
        "      'sidechain': \"DynaMine sidechain dynamics\",\n",
        "      'ppII': \"DynaMine conformational propensities: ppII (polyproline II)\",\n",
        "      'coil': \"DynaMine conformational propensities: Coil\",\n",
        "      'sheet': \"DynaMine conformational propensities: Sheet\",\n",
        "      'helix': \"DynaMine conformational propensities: Helix\",\n",
        "  }\n",
        "\n",
        "    PREDICTION_POSITION = {\n",
        "        'backbone':     0,\n",
        "        'sidechain':    1,\n",
        "        'ppII':         2,\n",
        "        'coil':         3,\n",
        "        'sheet':        4,\n",
        "        'helix':        5,\n",
        "    }\n",
        "\n",
        "elif DisoMine and not DynaMine and not EFoldMine:\n",
        "  NB_SUBPLOTS = 1\n",
        "  PREDICTION_TITLES = {\n",
        "       'disoMine': \"Disorder (disoMine)\"\n",
        "  }\n",
        "\n",
        "  PREDICTION_POSITION = {\n",
        "      'disoMine': 0\n",
        "  }\n",
        "\n",
        "elif EFoldMine and not DisoMine and not DynaMine:\n",
        "  NB_SUBPLOTS = 1\n",
        "\n",
        "  PREDICTION_TITLES = {\n",
        "      'earlyFolding': \"Early folding (EFoldMine)\"\n",
        "  }\n",
        "\n",
        "  PREDICTION_POSITION = {\n",
        "      'earlyFolding': 0\n",
        "  }\n",
        "\n",
        "elif DynaMine and DisoMine and not EFoldMine:\n",
        "  NB_SUBPLOTS = 7\n",
        "  PREDICTION_TITLES = {\n",
        "      'backbone': \"DynaMine backbone dynamics\",\n",
        "      'sidechain': \"DynaMine sidechain dynamics\",\n",
        "      'ppII': \"DynaMine conformational propensities: ppII (polyproline II)\",\n",
        "      'coil': \"DynaMine conformational propensities: Coil\",\n",
        "      'sheet': \"DynaMine conformational propensities: Sheet\",\n",
        "      'helix': \"DynaMine conformational propensities: Helix\",\n",
        "      'disoMine': \"Disorder (disoMine)\"\n",
        "  }\n",
        "\n",
        "  PREDICTION_POSITION = {\n",
        "      'backbone':     0,\n",
        "      'sidechain':    1,\n",
        "      'ppII':         2,\n",
        "      'coil':         3,\n",
        "      'sheet':        4,\n",
        "      'helix':        5,\n",
        "      'disoMine':     6\n",
        "  }\n",
        "\n",
        "elif DynaMine and EFoldMine and not DisoMine:\n",
        "  NB_SUBPLOTS = 7\n",
        "  PREDICTION_TITLES = {\n",
        "      'backbone': \"DynaMine backbone dynamics\",\n",
        "      'sidechain': \"DynaMine sidechain dynamics\",\n",
        "      'ppII': \"DynaMine conformational propensities: ppII (polyproline II)\",\n",
        "      'coil': \"DynaMine conformational propensities: Coil\",\n",
        "      'sheet': \"DynaMine conformational propensities: Sheet\",\n",
        "      'helix': \"DynaMine conformational propensities: Helix\",\n",
        "      'earlyFolding': \"Early folding (EFoldMine)\"\n",
        "  }\n",
        "\n",
        "  PREDICTION_POSITION = {\n",
        "      'backbone':     0,\n",
        "      'sidechain':    1,\n",
        "      'ppII':         2,\n",
        "      'coil':         3,\n",
        "      'sheet':        4,\n",
        "      'helix':        5,\n",
        "      'earlyFolding': 6\n",
        "  }\n",
        "elif EFoldMine and DisoMine and not DynaMine:\n",
        "  NB_SUBPLOTS = 2\n",
        "  PREDICTION_TITLES = {\n",
        "      'earlyFolding': \"Early folding (EFoldMine)\",\n",
        "      'disoMine': \"Disorder (disoMine)\"\n",
        "  }\n",
        "\n",
        "  PREDICTION_POSITION = {\n",
        "      'earlyFolding': 0,\n",
        "      'disoMine':     1\n",
        "  }\n",
        "elif DynaMine and DisoMine and EFoldMine:\n",
        "  NB_SUBPLOTS = 8\n",
        "  PREDICTION_TITLES = {\n",
        "      'backbone': \"DynaMine backbone dynamics\",\n",
        "      'sidechain': \"DynaMine sidechain dynamics\",\n",
        "      'ppII': \"DynaMine conformational propensities: ppII (polyproline II)\",\n",
        "      'coil': \"DynaMine conformational propensities: Coil\",\n",
        "      'sheet': \"DynaMine conformational propensities: Sheet\",\n",
        "      'helix': \"DynaMine conformational propensities: Helix\",\n",
        "      'earlyFolding': \"Early folding (EFoldMine)\",\n",
        "      'disoMine': \"Disorder (disoMine)\"\n",
        "  }\n",
        "\n",
        "  PREDICTION_POSITION = {\n",
        "      'backbone':     0,\n",
        "      'sidechain':    1,\n",
        "      'ppII':         2,\n",
        "      'coil':         3,\n",
        "      'sheet':        4,\n",
        "      'helix':        5,\n",
        "      'earlyFolding': 6,\n",
        "      'disoMine':     7\n",
        "  }\n",
        "\n",
        "AXIS_TITLES = {\n",
        "    \"x\": \"Residue position in the MSA\",\n",
        "    \"y\": \"Prediction values\"\n",
        "}\n",
        "\n",
        "\n",
        "import matplotlib as mpl\n",
        "mpl.use(\"Agg\")\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib.ticker import AutoMinorLocator\n",
        "import matplotlib.ticker as ticker\n"
      ],
      "metadata": {
        "cellView": "form",
        "id": "H9wUq47i0E7G"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 5.a. Plot the predicted biophysical conservation of the MSA\n",
        "#@markdown Once this cell has been executed, a file with plots will be saved in the folder \"Results\"\n",
        "#@markdown In case you want to zoom in on the plots, run the cell 5.b below .\n",
        "%%capture\n",
        "\n",
        "#@markdown The plots represent the biophysical behaviour of the proteins in the MSA.\n",
        "#@markdown There is one graph for every biophyisical property that is studied.\n",
        "#@markdown As the median, 1st and 3rd quartile and outliers are shown on the graphs,\n",
        "#@markdown the conservation of the biophysical properties can be observed.\n",
        "\n",
        "#@markdown The dots on the plots show the sequence occupency of the MSA, there are 4 different sizes representing 25%,50%,75% and 100% of the sequence occupency.\n",
        "\n",
        "# from matplotlib.collections import LineCollection\n",
        "\n",
        "def plot_biophysical_msa(jsondata_list_interest,jsondata_list_selected,selected_prots,sequences,freq_gap):\n",
        "\n",
        "    colors = ['blue', 'orange']\n",
        "    residues_count = len(jsondata_list_interest[0]['backbone']['median'])\n",
        "    sequences_count = len(sequences)\n",
        "    #Plot representation\n",
        "    fig, axs = plt.subplots(NB_SUBPLOTS)\n",
        "    fig.set_figwidth(20)\n",
        "    fig.set_figheight(50)\n",
        "\n",
        "    plt.suptitle(f\"Predicted biophysical properties of the MSA: {residues_count} aligned residues from {sequences_count} sequences\", fontsize=14)\n",
        "\n",
        "    # These for loops got too complicated, I have to think\n",
        "    # something simpler to handle the None values in the data\n",
        "    predictions = jsondata_list_interest[0].keys()\n",
        "    for prediction_index, biophys_data in enumerate(predictions):\n",
        "        if biophys_data == 'agmata':\n",
        "            continue\n",
        "\n",
        "        subplot_index_row = PREDICTION_POSITION[biophys_data]\n",
        "\n",
        "        ax = axs[subplot_index_row]\n",
        "        for data, col in zip(jsondata_list_interest, colors):\n",
        "            none_idx = []\n",
        "\n",
        "            for n in range(residues_count):\n",
        "                if data[biophys_data]['median'][n] == None \\\n",
        "                        or data[biophys_data][\n",
        "                    'firstQuartile'][n] == None \\\n",
        "                        or data[biophys_data][\n",
        "                    'thirdQuartile'][n] == None:\n",
        "                    none_idx.append(n)\n",
        "\n",
        "            range_list = []\n",
        "            for n in range(len(none_idx)):\n",
        "                try:\n",
        "                    if none_idx[n] + 1 != none_idx[n + 1]:\n",
        "                        range_list.append(\n",
        "                            (none_idx[n] + 1, none_idx[n + 1]))\n",
        "                    else:\n",
        "                        continue\n",
        "                except:\n",
        "                    if len(none_idx) == 1:\n",
        "                        range_list.append((0, none_idx[0]))\n",
        "                        range_list.append((none_idx[0] + 1, len(\n",
        "                            data[biophys_data][\n",
        "                                'median'])))\n",
        "\n",
        "                    else:\n",
        "                        range_list.append((0, none_idx[0]))\n",
        "                        range_list.append((none_idx[-1] + 1, len(\n",
        "                            data[biophys_data][\n",
        "                                'median'])))\n",
        "\n",
        "            # When there are None values in the data\n",
        "            if range_list:\n",
        "                for tuple in range_list:\n",
        "                    x = np.arange(tuple[0], tuple[1], 1)\n",
        "                    firstq = \\\n",
        "                        data[biophys_data][\n",
        "                            'firstQuartile'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    thirdq = \\\n",
        "                        data[biophys_data][\n",
        "                            'thirdQuartile'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    bottom = \\\n",
        "                        data[biophys_data][\n",
        "                            'bottomOutlier'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    top = \\\n",
        "                        data[biophys_data]['topOutlier'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    ax.fill_between(\n",
        "                        x, firstq, thirdq, alpha=0.3, color=col, label=\"1st-3rd Quartiles\")\n",
        "                    ax.fill_between(\n",
        "                        x, bottom, top, alpha=0.1, color=col, label=\"Outliers\")\n",
        "\n",
        "            # When there aren't None values in the data\n",
        "            else:\n",
        "                x = np.arange(0, len(\n",
        "                    data[biophys_data]['median']), 1)\n",
        "                firstq = data[biophys_data][\n",
        "                    'firstQuartile']\n",
        "                thirdq = data[biophys_data][\n",
        "                    'thirdQuartile']\n",
        "                bottom = data[biophys_data][\n",
        "                    'bottomOutlier']\n",
        "                top = data[biophys_data]['topOutlier']\n",
        "                ax.fill_between(\n",
        "                    x, firstq, thirdq, alpha=0.3, color=col, label=\"1st-3rd Quartiles\")\n",
        "                ax.fill_between(\n",
        "                    x, bottom, top, alpha=0.1, color=col, label=\"Outliers\")\n",
        "\n",
        "            ax.plot(data[biophys_data]['median'], linewidth=1.25, color=col, label=\"Median\")\n",
        "\n",
        "            #Add the selected protein if there are some\n",
        "            colors_2 = [\"magenta\", \"blue\", \"cyan\"] #adapt in function of number of proteins of interest, now max 3 can be studied simultanously\n",
        "            if len(selected_prot)>0:\n",
        "              for count,prot in enumerate(selected_prot):\n",
        "                  ax.plot(jsondata_list_selected[count][0][biophys_data], '-x', linewidth=1.5, color=colors_2[count], label=f\"Prediction {prot}\")\n",
        "\n",
        "            #To represent the amount of sequences at every column of the MSA\n",
        "            limits = [0.25,0.5,0.75,1]\n",
        "            limit_previous = 0\n",
        "            size = 20\n",
        "            for limit in limits:\n",
        "                plot_x = []\n",
        "                plot_y = []\n",
        "                for idx, freq in enumerate(freq_gap):\n",
        "                    if limit_previous < freq <= limit:\n",
        "                        plot_x.append(idx)\n",
        "                        plot_y.append(data[biophys_data]['median'][idx])\n",
        "                ax.scatter(plot_x,plot_y,marker='o',s=size, color=col)\n",
        "                limit_previous = limit\n",
        "                size += 15\n",
        "\n",
        "        ax.set_title(PREDICTION_TITLES[biophys_data])\n",
        "        ax.axis([0, residues_count-1, min(bottom)-0.05, max(top)+0.05])\n",
        "        # #labels\n",
        "        # labels = []\n",
        "        # count = 1\n",
        "        # while count < residues_count:\n",
        "        #     labels.append(str(count))\n",
        "        #     count +=4\n",
        "        # ax.xaxis.set_ticks(ticks=np.arange(0,residues_count+1,4))\n",
        "        # ax.xaxis.set_ticklabels(labels)\n",
        "        # ax.xaxis.set_minor_locator(AutoMinorLocator())\n",
        "        # ax.yaxis.set_minor_locator(AutoMinorLocator())\n",
        "        ax.set_ylabel(AXIS_TITLES['y'])\n",
        "        ax.set_xlabel(AXIS_TITLES['x'])\n",
        "        if biophys_data == 'backbone':\n",
        "            ax.axhline(y=1.0, color='green', linewidth= 1.5, linestyle='-.', label=\"Above: Membrane spaning\") #Membrane spaning\n",
        "            ax.axhline(y=0.8, color='orange', linewidth= 1.5, linestyle='-.', label=\"Above: Rigid\") #Membrane spaning\n",
        "            if min(bottom)-0.05 < 0.69:\n",
        "                ax.axhline(y=0.69, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Context dependent \\nBelow: Flexible\") #context dependent (either rigide or flexible)\n",
        "        if biophys_data == 'earlyFolding':\n",
        "            ax.axhline(y=0.169, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Likely to start folding\") #above: likely start protein folding process\n",
        "        if biophys_data == 'disoMine':\n",
        "            ax.axhline(y=0.5, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Likely to be disordered\") #above: likely disordered\n",
        "        ax.legend(ncol=1, bbox_to_anchor =(1.01,0.5), loc='center left')\n",
        "\n",
        "    plt.tight_layout()\n",
        "    fig.subplots_adjust(top=0.96, hspace = 0.2)\n",
        "    plt.savefig('/content/results/msa_biophysical_conservation.pdf')\n",
        "    plt.show()\n",
        "\n",
        "    return fig, axs\n",
        "\n",
        "sequences =  msaSuite.seqs\n",
        "fig, axs = plot_biophysical_msa(jsondata_list,jsondata_list_selected,selected_prot,sequences,freq_gap)"
      ],
      "metadata": {
        "cellView": "form",
        "id": "BuGj4zMI0TL6"
      },
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 5.b. Zoom in: Plot the predicted biophysical conservation of the MSA\n",
        "#@markdown Same as cell 5.a. but with zooms on specific regions of the MSA.\n",
        "%%capture\n",
        "\n",
        "#@markdown Enter in how many sections you want to cut your plots:\n",
        "\n",
        "n = 5 #@param {type:\"integer\"}\n",
        "#@markdown\n",
        "\n",
        "#@markdown Once this cell has been executed, a file with plots will be saved in the folder \"Results\".\n",
        "\n",
        "#@markdown The plots represent the biophysical behaviour of the proteins in the MSA.\n",
        "#@markdown There is one graph for every biophyisical property that is studied.\n",
        "#@markdown As the median, 1st and 3rd quartile and outliers are shown on the graphs,\n",
        "#@markdown the conservation of the biophysical properties can be observed.\n",
        "\n",
        "#@markdown The dots on the plots show the sequence occupency of the MSA, there are 4 different sizes representing 25%,50%,75% and 100% of the sequence occupency.\n",
        "\n",
        "#create sections\n",
        "residues_count = len(jsondata_list[0]['backbone']['median'])\n",
        "all_regions = []\n",
        "for elem in range(1,n+1):\n",
        "    label = \"Section%d_of_%d\" %(elem, n)\n",
        "    low = int(residues_count * ((1/n)*(elem-1)))\n",
        "    high = int(residues_count * ((1/n)*elem))\n",
        "    if all_regions == []:\n",
        "      low = 1\n",
        "    all_regions.append((label, (low, high)))\n",
        "\n",
        "def plot_biophysical_msa_zoom(jsondata_list_interest,jsondata_list_selected,selected_prot, sequences,freq_gap,all_regions):\n",
        "    colors = ['blue', 'orange']\n",
        "    residues_count = len(jsondata_list_interest[0]['backbone']['median'])\n",
        "    sequences_count = len(sequences)\n",
        "\n",
        "    for region in all_regions:\n",
        "        name_region = region[0]\n",
        "        lower_lim = region[1][0]\n",
        "        upper_lim =region[1][1]\n",
        "\n",
        "        #Plot representation\n",
        "        fig, axs = plt.subplots(NB_SUBPLOTS)\n",
        "        fig.set_figwidth(20)\n",
        "        fig.set_figheight(50)\n",
        "\n",
        "        fig.suptitle(f'Predicted biophysical properties of the MSA: \"{name_region}\" (residue {lower_lim} to residue {upper_lim-1}) from {sequences_count} sequences', fontsize=14)\n",
        "\n",
        "        # These for loops got too complicated, I have to think\n",
        "        # something simpler to handle the None values in the data\n",
        "        predictions = jsondata_list_interest[0].keys()\n",
        "        for prediction_index, biophys_data in enumerate(predictions):\n",
        "            if biophys_data == 'agmata':\n",
        "                continue\n",
        "\n",
        "            subplot_index_row = PREDICTION_POSITION[biophys_data]\n",
        "\n",
        "            ax = axs[subplot_index_row]\n",
        "            for data, col in zip(jsondata_list_interest, colors):\n",
        "                none_idx = []\n",
        "\n",
        "                for n in range(residues_count):\n",
        "                    if data[biophys_data]['median'][n] == None \\\n",
        "                            or data[biophys_data][\n",
        "                        'firstQuartile'][n] == None \\\n",
        "                            or data[biophys_data][\n",
        "                        'thirdQuartile'][n] == None:\n",
        "                        none_idx.append(n)\n",
        "\n",
        "                range_list = []\n",
        "                for n in range(len(none_idx)):\n",
        "                    try:\n",
        "                        if none_idx[n] + 1 != none_idx[n + 1]:\n",
        "                            range_list.append(\n",
        "                                (none_idx[n] + 1, none_idx[n + 1]))\n",
        "                        else:\n",
        "                            continue\n",
        "                    except:\n",
        "                        if len(none_idx) == 1:\n",
        "                            range_list.append((0, none_idx[0]))\n",
        "                            range_list.append((none_idx[0] + 1, len(\n",
        "                                data[biophys_data][\n",
        "                                    'median'])))\n",
        "\n",
        "                        else:\n",
        "                            range_list.append((0, none_idx[0]))\n",
        "                            range_list.append((none_idx[-1] + 1, len(\n",
        "                                data[biophys_data][\n",
        "                                    'median'])))\n",
        "\n",
        "                # When there are None values in the data\n",
        "                if range_list:\n",
        "                    for tuple in range_list:\n",
        "                        x = np.arange(tuple[0], tuple[1], 1)\n",
        "                        firstq = \\\n",
        "                            data[biophys_data][\n",
        "                                'firstQuartile'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        thirdq = \\\n",
        "                            data[biophys_data][\n",
        "                                'thirdQuartile'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        bottom = \\\n",
        "                            data[biophys_data][\n",
        "                                'bottomOutlier'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        top = \\\n",
        "                            data[biophys_data]['topOutlier'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        ax.fill_between(\n",
        "                            x[lower_lim:upper_lim], firstq[lower_lim:upper_lim], thirdq[lower_lim:upper_lim], alpha=0.3, color=col, label=\"1st & 3rd quartile\")\n",
        "                        ax.fill_between(\n",
        "                            x[lower_lim:upper_lim], bottom[lower_lim:upper_lim], top[lower_lim:upper_lim], alpha=0.1, color=col, label=\"Outliers\")\n",
        "\n",
        "                # When there aren't None values in the data\n",
        "                else:\n",
        "                    x = np.arange(lower_lim,upper_lim, 1)\n",
        "                    firstq = data[biophys_data][\n",
        "                        'firstQuartile'][lower_lim:upper_lim]\n",
        "                    thirdq = data[biophys_data][\n",
        "                        'thirdQuartile'][lower_lim:upper_lim]\n",
        "                    bottom = data[biophys_data][\n",
        "                        'bottomOutlier'][lower_lim:upper_lim]\n",
        "                    top = data[biophys_data]['topOutlier'][lower_lim:upper_lim]\n",
        "                    ax.fill_between(\n",
        "                        x, firstq, thirdq, alpha=0.3, color=col, label=\"1st & 3rd quartile\")\n",
        "                    ax.fill_between(\n",
        "                        x, bottom, top, alpha=0.1, color=col, label=\"Outliers\")\n",
        "\n",
        "                    ax.plot(data[biophys_data]['median'],linewidth=1.25, color=\"black\", label= \"Median\")\n",
        "\n",
        "                    #Add selected proteins if there are some\n",
        "                    colors_2 = [\"magenta\", \"blue\", \"cyan\"] #adapt in function of number of proteins of interest, now max 3 can be studied simultanously\n",
        "                    if len(selected_prot)>0:\n",
        "                      for count,prot in enumerate(selected_prot):\n",
        "                          ax.plot(jsondata_list_selected[count][0][biophys_data], '-x', linewidth=1.5, color=colors_2[count], label=f\"Prediction {prot}\")\n",
        "\n",
        "                #To represent the amount of sequences occupency at every column of the MSA\n",
        "                limits = [0.25,0.5,0.75,1] #we divide the occupency in 4 groups: i<25%, 25%<i<50%, 50%<i<75%, i>75%\n",
        "                #depending on the group, the size of the dot will change\n",
        "                limit_previous = 0\n",
        "                size = 20\n",
        "                for limit in limits:\n",
        "                    plot_x = []\n",
        "                    plot_y = []\n",
        "                    for idx, freq in enumerate(freq_gap):\n",
        "                        if limit_previous < freq <= limit:\n",
        "                            plot_x.append(idx)\n",
        "                            plot_y.append(data[biophys_data]['median'][idx])\n",
        "                    if limit == 1:\n",
        "                      ax.scatter(plot_x,plot_y,marker='o',s=size, color=\"grey\", label = \"Size dot ~ freq. seq. occupency\")\n",
        "                    else:\n",
        "                      ax.scatter(plot_x,plot_y,marker='o',s=size, color=\"grey\")\n",
        "                    limit_previous = limit\n",
        "                    size += 15\n",
        "\n",
        "            ax.set_title(PREDICTION_TITLES[biophys_data])\n",
        "            ax.axis([lower_lim,upper_lim, min(bottom)-0.05, max(top)+0.05])\n",
        "            ax.set_ylabel(AXIS_TITLES['y'])\n",
        "            ax.set_xlabel(AXIS_TITLES['x'])\n",
        "            # ax.set_xticks(list(np.arange(lower_lim,upper_lim,10)))\n",
        "            if biophys_data == 'backbone':\n",
        "                ax.axhline(y=1.0, color='green', linewidth= 1.25, linestyle='-.', label=\"Above: Membrane spaning\") #Membrane spaning\n",
        "                ax.axhline(y=0.8, color='orange', linewidth= 1.25, linestyle='-.', label=\"Above: Rigid\") #Membrane spaning\n",
        "                if min(bottom)-0.05 < 0.69:\n",
        "                    ax.axhline(y=0.69, color='red', linewidth= 1.25, linestyle='-.', label=\"Above: Context dependent \\nBelow: Flexible\") #context dependent (either rigide or flexible)\n",
        "            if biophys_data == 'earlyFolding':\n",
        "                ax.axhline(y=0.169, color='red', linewidth= 1.25, linestyle='-.', label=\"Above: Likely to start folding\") #above: likely start protein folding process\n",
        "            if biophys_data == 'disoMine':\n",
        "                ax.axhline(y=0.5, color='red', linewidth= 1.25, linestyle='-.', label=\"Above: Likely to be disordered\") #above: likely disordered\n",
        "            ax.legend(ncol=1, bbox_to_anchor =(1.01,0.5), loc='center left')\n",
        "\n",
        "\n",
        "        plt.tight_layout()\n",
        "        fig.subplots_adjust(top=0.96, hspace = 0.2)\n",
        "        plt.savefig('/content/results/'+ name_region + '_msa_biophysical_conservation.pdf')\n",
        "        plt.show()\n",
        "\n",
        "    return fig, axs\n",
        "\n",
        "sequences =  msaSuite.seqs\n",
        "fig, axs = plot_biophysical_msa_zoom(jsondata_list,jsondata_list_selected,selected_prot, sequences,freq_gap,all_regions)"
      ],
      "metadata": {
        "cellView": "form",
        "id": "dsALZYId0dgC"
      },
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 5.c. Plot the biophysical and sequence conservation of the MSA\n",
        "#@markdown Once this cell has been executed, a file with plots will be saved in the folder \"Results\"\n",
        "#@markdown In case you want to zoom in on the plots, run the cell 5.d below.\n",
        "%%capture\n",
        "\n",
        "#@markdown The plots represent the biophysical behaviour of the proteins in the MSA as well as their sequence conservation.\n",
        "#@markdown There is one graph for every biophyisical property that is studied.\n",
        "#@markdown As the median, 1st and 3rd quartile and outliers are shown on the graphs,\n",
        "#@markdown the conservation of the biophysical properties can be observed.\n",
        "\n",
        "#@markdown The dots on the plots show the sequence occupency of the MSA, there are 4 different sizes representing 25%,50%,75% and 100% of the sequence occupency.\n",
        "\n",
        "def plot_entropy_biphysical_msa(jsondata_list, jsondata_list_selected, selected_prot, sequences, conservation_AA_list, freq_gap):\n",
        "    sequences_count = len(sequences)\n",
        "    residues_count = len(jsondata_list[0]['backbone']['median'])\n",
        "\n",
        "    #Color map for conservation sequence\n",
        "    cmap = mpl.cm.Blues(np.linspace(0,1,20))\n",
        "    cmap = mpl.colors.ListedColormap(cmap[:5])\n",
        "    colors = ['blue', 'orange']\n",
        "\n",
        "    #Plot representation\n",
        "    fig, axs = plt.subplots(NB_SUBPLOTS)\n",
        "    fig.set_figwidth(20)\n",
        "    fig.set_figheight(50)\n",
        "\n",
        "    fig.suptitle(f\"Predicted biophysical properties of the MSA: {residues_count} aligned residues from {len(sequences)} sequences\", fontsize=16)\n",
        "\n",
        "    predictions = jsondata_list[0].keys()\n",
        "    for prediction_index, biophys_data in enumerate(predictions):\n",
        "        if biophys_data == 'agmata':\n",
        "            continue\n",
        "\n",
        "        subplot_index_row = PREDICTION_POSITION[biophys_data]\n",
        "\n",
        "        ax = axs[subplot_index_row]\n",
        "        for data, col in zip(jsondata_list, colors):\n",
        "            none_idx = []\n",
        "\n",
        "            for n in range(residues_count):\n",
        "                if data[biophys_data]['median'][n] == None \\\n",
        "                        or data[biophys_data][\n",
        "                    'firstQuartile'][n] == None \\\n",
        "                        or data[biophys_data][\n",
        "                    'thirdQuartile'][n] == None:\n",
        "                    none_idx.append(n)\n",
        "\n",
        "            range_list = []\n",
        "            for n in range(len(none_idx)):\n",
        "                try:\n",
        "                    if none_idx[n] + 1 != none_idx[n + 1]:\n",
        "                        range_list.append(\n",
        "                            (none_idx[n] + 1, none_idx[n + 1]))\n",
        "                    else:\n",
        "                        continue\n",
        "                except:\n",
        "                    if len(none_idx) == 1:\n",
        "                        range_list.append((0, none_idx[0]))\n",
        "                        range_list.append((none_idx[0] + 1, len(\n",
        "                            data[biophys_data][\n",
        "                                'median'])))\n",
        "\n",
        "                    else:\n",
        "                        range_list.append((0, none_idx[0]))\n",
        "                        range_list.append((none_idx[-1] + 1, len(\n",
        "                            data[biophys_data][\n",
        "                                'median'])))\n",
        "\n",
        "            # When there are None values in the data\n",
        "            if range_list:\n",
        "                for tuple in range_list:\n",
        "                    x = np.arange(tuple[0], tuple[1], 1)\n",
        "                    firstq = \\\n",
        "                        data[biophys_data][\n",
        "                            'firstQuartile'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    thirdq = \\\n",
        "                        data[biophys_data][\n",
        "                            'thirdQuartile'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    bottom = \\\n",
        "                        data[biophys_data][\n",
        "                            'bottomOutlier'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    top = \\\n",
        "                        data[biophys_data]['topOutlier'][\n",
        "                        tuple[0]:tuple[1]]\n",
        "                    ax.plot(x,firstq, linewidth=0.5, color=\"black\", label=\"1st & 3rd quartile\")\n",
        "                    ax.plot(x,thirdq, linewidth=0.5, color=\"black\")\n",
        "                    ax.plot(x,bottom, alpha=0.25, linestyle =\"--\", color=\"black\", label=\"Bottom & Top outliers\")\n",
        "                    ax.plot(x,top, alpha=0.25, linestyle =\"--\", color=\"black\")\n",
        "\n",
        "            # When there aren't None values in the data\n",
        "            else:\n",
        "                x = np.arange(0, len(\n",
        "                    data[biophys_data]['median']), 1)\n",
        "                firstq = data[biophys_data][\n",
        "                    'firstQuartile']\n",
        "                thirdq = data[biophys_data][\n",
        "                    'thirdQuartile']\n",
        "                bottom = data[biophys_data][\n",
        "                    'bottomOutlier']\n",
        "                top = data[biophys_data]['topOutlier']\n",
        "\n",
        "                ax.plot(x,firstq, linewidth=0.5, color=\"black\", label=\"1st & 3rd quartile\")\n",
        "                ax.plot(x,thirdq, linewidth=0.5, color=\"black\")\n",
        "                ax.plot(x,bottom, alpha=0.25, linestyle =\"--\", color=\"black\", label=\"Bottom & Top outliers\")\n",
        "                ax.plot(x,top, alpha=0.25, linestyle =\"--\", color=\"black\")\n",
        "\n",
        "                ax.plot(data[biophys_data]['median'], linewidth=1.25, color=\"black\", label= \"Median\")\n",
        "\n",
        "                #Add selected proteins if there are some\n",
        "                colors_2 = [\"magenta\", \"blue\", \"cyan\"] #adapt in function of number of proteins of interest, now max 3 can be studied simultanously\n",
        "                if len(selected_prot)>0:\n",
        "                  for count,prot in enumerate(selected_prot):\n",
        "                      ax.plot(jsondata_list_selected[count][0][biophys_data], '-x', linewidth=1.5, color=colors_2[count], label=f\"Prediction {prot}\")\n",
        "\n",
        "                #plot sequence conservation\n",
        "                entropy_values = conservation_AA_list\n",
        "\n",
        "                extent = [0, residues_count, min(bottom)-0.05, max(top)+0.05]\n",
        "                ax.imshow(np.array(entropy_values)[np.newaxis,:], cmap=cmap, extent=extent,aspect = \"auto\", vmin=0, vmax=1)\n",
        "                cbar = ax.figure.colorbar(\n",
        "                    mpl.cm.ScalarMappable(norm=None, cmap=cmap), shrink=1, pad = 0.03, ax=ax, ticks=[0, 1])\n",
        "                cbar.ax.set_yticklabels(['0%', '100%'])\n",
        "                cbar.ax.set_ylabel('Sequence conservation', rotation=270)\n",
        "\n",
        "            #To represent the amount of sequences occupency at every column of the MSA\n",
        "            limits = [0.25,0.5,0.75,1] #we divide the occupency in 4 groups: i<25%, 25%<i<50%, 50%<i<75%, i>75%\n",
        "            #depending on the group, the size of the dot will change\n",
        "            limit_previous = 0\n",
        "            size = 20\n",
        "            for limit in limits:\n",
        "                plot_x = []\n",
        "                plot_y = []\n",
        "                for idx, freq in enumerate(freq_gap):\n",
        "                    if limit_previous < freq <= limit:\n",
        "                        plot_x.append(idx)\n",
        "                        plot_y.append(data[biophys_data]['median'][idx])\n",
        "                if limit == 1:\n",
        "                  ax.scatter(plot_x,plot_y,marker='o',s=size, color=\"grey\", label = \"Size dot ~ freq. seq. occupency\")\n",
        "                else:\n",
        "                  ax.scatter(plot_x,plot_y,marker='o',s=size, color=\"grey\")\n",
        "                limit_previous = limit\n",
        "                size += 15\n",
        "\n",
        "        ax.set_title(PREDICTION_TITLES[biophys_data] , fontsize=14)\n",
        "        ax.axis([0, residues_count, min(bottom)-0.05, max(top)+0.05])\n",
        "        # ax.xaxis.set_major_locator(ticker.FixedLocator(np.arange(0,residues_count,4)))\n",
        "\n",
        "        # #labels\n",
        "        # labels = []\n",
        "        # count = 1\n",
        "        # while count < residues_count+1:\n",
        "        #     labels.append(str(count))\n",
        "        #     count +=4\n",
        "        # ax.xaxis.set_major_formatter(ticker.FixedFormatter(labels))\n",
        "        # ax.xaxis.set_minor_locator(ticker.FixedLocator(np.arange(0.99,residues_count+0.99,1)))\n",
        "        # ax.yaxis.set_minor_locator(AutoMinorLocator())\n",
        "        ax.set_ylabel(AXIS_TITLES['y'], fontsize=10)\n",
        "        ax.set_xlabel(AXIS_TITLES['x'], fontsize=10)\n",
        "\n",
        "        if biophys_data == 'backbone':\n",
        "            ax.axhline(y=1.0, color='green', linewidth= 1.5, linestyle='-.', label=\"Above: Membrane spaning\") #Membrane spaning\n",
        "            ax.axhline(y=0.8, color='orange', linewidth= 1.5, linestyle='-.', label=\"Above: Rigid\") #Membrane spaning\n",
        "            if min(bottom)-0.05 < 0.69:\n",
        "                ax.axhline(y=0.69, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Context dependent \\nBelow: Flexible\") #context dependent (either rigide or flexible)\n",
        "        if biophys_data == 'earlyFolding':\n",
        "            ax.axhline(y=0.169, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Likely to start folding\") #above: likely start protein folding process\n",
        "        if biophys_data == 'disoMine':\n",
        "            ax.axhline(y=0.5, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Likely to be disordered\") #above: likely disordered\n",
        "        ax.legend(ncol=1, bbox_to_anchor =(1.1,0.5), loc='center left')\n",
        "\n",
        "\n",
        "    plt.tight_layout()\n",
        "    fig.subplots_adjust(top=0.96, hspace = 0.6)\n",
        "\n",
        "    plt.savefig('/content/results/msa_sequence_biophysical_conservation.pdf')\n",
        "\n",
        "    return fig, axs\n",
        "\n",
        "sequences =  msaSuite.seqs\n",
        "fig, axs = plot_entropy_biphysical_msa(jsondata_list,jsondata_list_selected,selected_prot,sequences,conservation_AA_list,freq_gap)"
      ],
      "metadata": {
        "cellView": "form",
        "id": "BfW0yTyl0pUy"
      },
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 5.d. Zoom in: Plot the sequence and biophysical conservation of the MSA\n",
        "#@markdown Same as cell 5.c. but with zooms on specific regions of the MSA.\n",
        "%%capture\n",
        "\n",
        "#@markdown Enter in how many sections you want to cut your plots:\n",
        "n = 5 #@param {type:\"integer\"}\n",
        "\n",
        "#@markdown\n",
        "\n",
        "#@markdown The plots represent the biophysical behaviour of the proteins in the MSA as well as their sequence conservation.\n",
        "#@markdown There is one graph for every biophyisical property that is studied.\n",
        "#@markdown As the median, 1st and 3rd quartile and outliers are shown on the graphs,\n",
        "#@markdown the conservation of the biophysical properties can be observed.\n",
        "\n",
        "#@markdown The dots on the plots show the sequence occupency of the MSA, there are 4 different sizes representing 25%,50%,75% and 100% of the sequence occupency.\n",
        "\n",
        "#create sections\n",
        "residues_count = len(jsondata_list[0]['backbone']['median'])\n",
        "all_regions = []\n",
        "for elem in range(1,n+1):\n",
        "    label = \"Section%d_of_%d\" %(elem, n)\n",
        "    low = int(residues_count * ((1/n)*(elem-1)))\n",
        "    high = int(residues_count * ((1/n)*elem))\n",
        "    if all_regions == []:\n",
        "      low = 1\n",
        "    all_regions.append((label, (low, high)))\n",
        "print(all_regions)\n",
        "#@markdown The dots on the plots show the sequence occupency of the MSA, there are 4 different sizes representing 25%,50%,75% and 100% of the sequence occupency.\n",
        "\n",
        "def plot_entropy_biphysical_msa_zoom(jsondata_list,jsondata_list_selected,selected_prot, sequences, all_regions, conservation_AA_list, freq_gap):\n",
        "    sequences_count = len(sequences)\n",
        "    residues_count = len(jsondata_list[0]['backbone']['median'])\n",
        "\n",
        "    #Color map for conservation sequence\n",
        "    cmap = mpl.cm.Blues(np.linspace(0,1,20))\n",
        "    cmap = mpl.colors.ListedColormap(cmap[:5])\n",
        "\n",
        "    colors = ['blue', 'orange']\n",
        "\n",
        "    for region in all_regions:\n",
        "        name_region = region[0]\n",
        "        lower_lim = region[1][0]\n",
        "        upper_lim = region[1][1]\n",
        "\n",
        "        #Plot representation\n",
        "        fig, axs = plt.subplots(8)\n",
        "        fig.set_figwidth(20)\n",
        "        fig.set_figheight(50)\n",
        "\n",
        "        fig.suptitle(f'Predicted biophysical properties of the MSA: \"{name_region}\" (residue {lower_lim} to residue {upper_lim-1}) from {sequences_count} sequences', fontsize=14)\n",
        "\n",
        "        # These for loops got too complicated, I have to think\n",
        "        # something simpler to handle the None values in the data\n",
        "        predictions = jsondata_list[0].keys()\n",
        "        for prediction_index, biophys_data in enumerate(predictions):\n",
        "            if biophys_data == 'agmata':\n",
        "                continue\n",
        "\n",
        "            subplot_index_row = PREDICTION_POSITION[biophys_data]\n",
        "\n",
        "            ax = axs[subplot_index_row]\n",
        "            for data, col in zip(jsondata_list, colors):\n",
        "                none_idx = []\n",
        "\n",
        "                for n in range(residues_count):\n",
        "                    if data[biophys_data]['median'][n] == None \\\n",
        "                            or data[biophys_data][\n",
        "                        'firstQuartile'][n] == None \\\n",
        "                            or data[biophys_data][\n",
        "                        'thirdQuartile'][n] == None:\n",
        "                        none_idx.append(n)\n",
        "\n",
        "                range_list = []\n",
        "                for n in range(len(none_idx)):\n",
        "                    try:\n",
        "                        if none_idx[n] + 1 != none_idx[n + 1]:\n",
        "                            range_list.append(\n",
        "                                (none_idx[n] + 1, none_idx[n + 1]))\n",
        "                        else:\n",
        "                            continue\n",
        "                    except:\n",
        "                        if len(none_idx) == 1:\n",
        "                            range_list.append((0, none_idx[0]))\n",
        "                            range_list.append((none_idx[0] + 1, len(\n",
        "                                data[biophys_data][\n",
        "                                    'median'])))\n",
        "\n",
        "                        else:\n",
        "                            range_list.append((0, none_idx[0]))\n",
        "                            range_list.append((none_idx[-1] + 1, len(\n",
        "                                data[biophys_data][\n",
        "                                    'median'])))\n",
        "\n",
        "                # When there are None values in the data\n",
        "                if range_list:\n",
        "                    for tuple in range_list:\n",
        "                        x = np.arange(tuple[0], tuple[1], 1)\n",
        "                        firstq = \\\n",
        "                            data[biophys_data][\n",
        "                                'firstQuartile'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        thirdq = \\\n",
        "                            data[biophys_data][\n",
        "                                'thirdQuartile'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        bottom = \\\n",
        "                            data[biophys_data][\n",
        "                                'bottomOutlier'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        top = \\\n",
        "                            data[biophys_data]['topOutlier'][\n",
        "                            tuple[0]:tuple[1]]\n",
        "                        ax.plot(x,firstq, linewidth=0.25, color=\"black\", label=\"1st & 3rd quartile\")\n",
        "                        ax.plot(x,thirdq, linewidth=0.25, color=\"black\")\n",
        "                        ax.plot(x,bottom, alpha=0.25, linestyle =\"--\", color=\"black\", label=\"Bottom & Top outliers\")\n",
        "                        ax.plot(x,top, alpha=0.25, linestyle =\"--\", color=\"black\")\n",
        "\n",
        "                # When there aren't None values in the data\n",
        "                else:\n",
        "                    x = np.arange(0, len(\n",
        "                        data[biophys_data]['median']), 1)\n",
        "                    firstq = data[biophys_data][\n",
        "                        'firstQuartile']\n",
        "                    thirdq = data[biophys_data][\n",
        "                        'thirdQuartile']\n",
        "                    bottom = data[biophys_data][\n",
        "                        'bottomOutlier']\n",
        "                    top = data[biophys_data]['topOutlier']\n",
        "\n",
        "                    ax.plot(x,firstq, linewidth=0.5, color=\"black\", label=\"1st & 3rd quartile\")\n",
        "                    ax.plot(x,thirdq, linewidth=0.5, color=\"black\")\n",
        "                    ax.plot(x,bottom, alpha=0.25, linestyle =\"--\", color=\"black\", label=\"Bottom & Top outliers\")\n",
        "                    ax.plot(x,top, alpha=0.25, linestyle =\"--\", color=\"black\")\n",
        "\n",
        "                    ax.plot(data[biophys_data]['median'], linewidth=1.25, color=\"black\", label= \"Median\")\n",
        "\n",
        "                    #Add selected proteins if there are some\n",
        "                    colors_2 = [\"magenta\", \"blue\", \"cyan\"] #adapt in function of number of proteins of interest, now max 3 can be studied simultanously\n",
        "                    if len(selected_prot)>0:\n",
        "                      for count,prot in enumerate(selected_prot):\n",
        "                          ax.plot(jsondata_list_selected[count][0][biophys_data], '-x', linewidth=1.5, color=colors_2[count], label=f\"Prediction {prot}\")\n",
        "\n",
        "                    #Sequence conservation\n",
        "                    entropy_values = conservation_AA_list\n",
        "                    extent = [0, residues_count, min(bottom)-0.05, max(top)+0.05]\n",
        "                    ax.imshow(np.array(entropy_values)[np.newaxis,:], cmap=cmap, extent=extent,aspect = \"auto\", vmin=0, vmax=1)\n",
        "                    cbar = ax.figure.colorbar(\n",
        "                        mpl.cm.ScalarMappable(norm=None, cmap=cmap), shrink=1, pad = 0.03, ax=ax, ticks=[0, 1])\n",
        "                    cbar.ax.set_yticklabels(['0%', '100%'])\n",
        "                    cbar.ax.set_ylabel('Sequence conservation', rotation=270)\n",
        "\n",
        "                #To represent the amount of sequences at every column of the MSA\n",
        "                limits = [0.25,0.5,0.75,1]\n",
        "                limit_previous = 0\n",
        "                size = 20\n",
        "                for limit in limits:\n",
        "                    plot_x = []\n",
        "                    plot_y = []\n",
        "                    for idx, freq in enumerate(freq_gap):\n",
        "                        if limit_previous < freq <= limit:\n",
        "                            plot_x.append(idx)\n",
        "                            plot_y.append(data[biophys_data]['median'][idx])\n",
        "                    if limit == 1:\n",
        "                      ax.scatter(plot_x,plot_y,marker='o',s=size, color=\"grey\", label = \"Size dot ~ freq. seq. occupency\")\n",
        "                    else:\n",
        "                      ax.scatter(plot_x,plot_y,marker='o',s=size, color=\"grey\")\n",
        "                    limit_previous = limit\n",
        "                    size += 15\n",
        "\n",
        "            ax.set_title(PREDICTION_TITLES[biophys_data] , fontsize=14)\n",
        "            ax.axis([lower_lim, upper_lim, min(bottom)-0.05, max(top)+0.05])\n",
        "            # ax.xaxis.set_major_locator(ticker.FixedLocator(np.arange(lower_lim,upper_lim-1,4)))\n",
        "            # #labels\n",
        "            # labels = []\n",
        "            # count = lower_lim + 1\n",
        "            # while count < upper_lim+1:\n",
        "            #     labels.append(str(count))\n",
        "            #     count +=4\n",
        "            # ax.xaxis.set_major_formatter(ticker.FixedFormatter(labels))\n",
        "            # ax.xaxis.set_minor_locator(ticker.FixedLocator(np.arange(0.99,residues_count+0.99,1)))\n",
        "            # ax.yaxis.set_minor_locator(AutoMinorLocator())\n",
        "            ax.set_ylabel(AXIS_TITLES['y'], fontsize=10)\n",
        "            ax.set_xlabel(AXIS_TITLES['x'], fontsize=10)\n",
        "\n",
        "\n",
        "            if biophys_data == 'backbone':\n",
        "                ax.axhline(y=1.0, color='green', linewidth= 1.5, linestyle='-.', label=\"Above: Membrane spaning\") #Membrane spaning\n",
        "                ax.axhline(y=0.8, color='orange', linewidth= 1.5, linestyle='-.', label=\"Above: Rigid\") #Membrane spaning\n",
        "                if min(bottom)-0.05 < 0.69:\n",
        "                    ax.axhline(y=0.69, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Context dependent \\nBelow: Flexible\") #context dependent (either rigide or flexible)\n",
        "            if biophys_data == 'earlyFolding':\n",
        "                ax.axhline(y=0.169, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Likely to start folding\") #above: likely start protein folding process\n",
        "            if biophys_data == 'disoMine':\n",
        "                ax.axhline(y=0.5, color='red', linewidth= 1.5, linestyle='-.', label=\"Above: Likely to be disordered\") #above: likely disordered\n",
        "            ax.legend(ncol=1, bbox_to_anchor =(1.1,0.5), loc='center left')\n",
        "\n",
        "\n",
        "        plt.tight_layout()\n",
        "        fig.subplots_adjust(top=0.96, hspace = 0.6)\n",
        "\n",
        "        plt.savefig('/content/results/'+ name_region + '_msa_sequence_biophysical_conservation.pdf')\n",
        "\n",
        "        plt.show()\n",
        "\n",
        "    return fig, axs\n",
        "\n",
        "sequences =  msaSuite.seqs\n",
        "fig, axs = plot_entropy_biphysical_msa_zoom(jsondata_list,jsondata_list_selected,selected_prot, sequences, all_regions, conservation_AA_list, freq_gap)"
      ],
      "metadata": {
        "cellView": "form",
        "id": "EEkb2bPs1ANW"
      },
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 6. Get the Gaussian Mixture Model (GMM) scores\n",
        "%%capture\n",
        "#@markdown Once this cell has been executed, the GMM scores will be computed and the most different residues will be shown.\n",
        "#@markdown The results are also saved in the Results folder.\n",
        "\n",
        "#@markdown Note: only the following biophysical properties are considered to calculate\n",
        "#@markdown the GMM scores: backbone & sidechain dynamics, coil, sheet & helix conformation propensity,\n",
        "#@markdown early folding propsensity and disorder propensity.\n",
        "\n",
        "def GMM_scores (jsondata_list, single_preds):\n",
        "    #Train predictor\n",
        "    data_full = []\n",
        "    for aln_pos in range(len(jsondata_list[0]['backbone']['median'])):\n",
        "        pred_vector = []\n",
        "        for biophys in ['backbone', 'sidechain', 'coil', 'sheet', 'helix',\n",
        "                        'earlyFolding', 'disoMine']:\n",
        "            pred_vector.append(jsondata_list[0][biophys]['median'][aln_pos])\n",
        "        data_full.append(pred_vector)\n",
        "    X_train = np.vstack(data_full)\n",
        "    clf = mixture.GaussianMixture(n_components = 1, covariance_type = 'full',\n",
        "                                      verbose = 2, verbose_interval = 1)\n",
        "    clf.fit(X_train)\n",
        "\n",
        "    #Predict GMM scores\n",
        "    gmm_dict = {}\n",
        "    scores = []\n",
        "    gmm_info = {}\n",
        "    sequences = {}\n",
        "\n",
        "    for prot in single_preds.keys():\n",
        "      if prot != \"sequence\":\n",
        "        full_pred = []\n",
        "        for res in range(len(single_preds[prot]['backbone'])):\n",
        "            pred_vector = []\n",
        "            for biophys in ['backbone', 'sidechain', 'coil', 'sheet', 'helix',\n",
        "                            'earlyFolding', \"disoMine\"]:\n",
        "                if single_preds[prot][biophys][res] != None:\n",
        "                    pred_vector.append(single_preds[prot][biophys][res])\n",
        "\n",
        "            full_pred.append(pred_vector)\n",
        "\n",
        "        full_pred = [list for list in full_pred if list]\n",
        "\n",
        "        preds = np.vstack(full_pred)\n",
        "        gmm_info[prot] = clf.score_samples(preds).tolist()\n",
        "        scores.extend(clf.score_samples(preds).tolist())\n",
        "        gmm_info[prot] = clf.score_samples(preds).tolist()\n",
        "        scores.extend(clf.score_samples(preds).tolist())\n",
        "\n",
        "\n",
        "    for prot in predictions_single_seq[\"sequence\"].keys():\n",
        "      sequences[prot]=predictions_single_seq[\"sequence\"][prot]\n",
        "\n",
        "    #most different biophysical behaviour\n",
        "    perc_5 = np.percentile(scores, 5)\n",
        "    perc_1 = np.percentile(scores, 1)\n",
        "\n",
        "    for key in gmm_info.keys():\n",
        "\n",
        "        worst_5perc = []\n",
        "        worst_1perc = []\n",
        "\n",
        "        #find the residues that are the most different compared to the other residues at the same position in the MSA\n",
        "        #we analyze the difference in biophysical behavior between the residues\n",
        "        for i in gmm_info[key]:\n",
        "            if i <= perc_5:\n",
        "                idx = gmm_info[key].index(i)\n",
        "                AA = sequences[key][idx]\n",
        "                worst_5perc.append(AA+str(idx+1))\n",
        "            if i <= perc_1:\n",
        "                idx = gmm_info[key].index(i)\n",
        "                AA = sequences[key][idx]\n",
        "                worst_1perc.append(AA+str(idx+1))\n",
        "\n",
        "        #To avoid to have the same AA in both lists\n",
        "        for AA in worst_5perc:\n",
        "            if AA in worst_1perc:\n",
        "                worst_5perc.remove(AA)\n",
        "\n",
        "        worst = [worst_5perc,worst_1perc]\n",
        "        new_list = []\n",
        "        for worst_ in worst:\n",
        "            consecutive = []\n",
        "            for position in range(len(worst_)-1):\n",
        "                idx = worst_[position][1:]\n",
        "                next = worst_[position+1][1:]\n",
        "                if int(idx) +1 == int(next):\n",
        "                    consecutive.append(True) #True means there is a consecutive one coming\n",
        "                else:\n",
        "                    consecutive.append(False)\n",
        "            consecutive.append(False)\n",
        "\n",
        "            found_first = False\n",
        "            new = []\n",
        "            for i in range(len(worst_)):\n",
        "                if consecutive[i] == True and found_first == False:\n",
        "                    found_first = True\n",
        "                    AA_first = worst_[i]\n",
        "                if consecutive[i] == False and found_first == True:\n",
        "                    AA_last = worst_[i]\n",
        "                    concatenated = AA_first+\"-\"+AA_last\n",
        "                    new.append(concatenated)\n",
        "                    found_first = False\n",
        "                elif consecutive[i] == False and found_first == False:\n",
        "                    new.append(worst_[i])\n",
        "            new_list.append(new)\n",
        "        worst_5perc_new = new_list[0]\n",
        "        worst_1perc_new = new_list[1]\n",
        "\n",
        "\n",
        "        #Worst 5% means the residues are lower than the 5th percentille\n",
        "        #Worst 1% means the residues are lower than the 1th percentill\n",
        "        gmm_dict[key] = {\"GMMscore\": gmm_info[key], \"Worst 5%\":worst_5perc_new, \"Worst 1%\": worst_1perc_new}\n",
        "\n",
        "    with open('/content/results/GMM_quantification.json', 'w') as fp:\n",
        "        json.dump(gmm_dict, fp)\n",
        "\n",
        "    return gmm_dict\n",
        "\n",
        "jsondata_list = [msaSuite.alignedPredictionDistribs]\n",
        "single_preds = msaSuite.allAlignedPredictions\n",
        "gmm_scores = GMM_scores(jsondata_list, single_preds)\n",
        "\n",
        "\n",
        "print(\"\\n\")\n",
        "print(\"GMM score analysis to find the residues that are the most different:\")\n",
        "for prot in gmm_scores.keys():\n",
        "    print(prot)\n",
        "    print(\"Worst 5%:\")\n",
        "    print(str(gmm_scores[prot][\"Worst 5%\"]))\n",
        "    print(\"Worst 1%:\")\n",
        "    print(str(gmm_scores[prot][\"Worst 1%\"]))\n",
        "    print(\"\\n\")"
      ],
      "metadata": {
        "cellView": "form",
        "id": "Su7Su3SC1T6z"
      },
      "execution_count": 13,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title 7. Download the predictions\n",
        "\n",
        "#@markdown Once this cell has been executed, a zip-archive with\n",
        "#@markdown the predictions, the corresponding plots and the GMM scores will created and is ready to be downloaded.\n",
        "\n",
        "final_results = {}\n",
        "\n",
        "predictions = msaSuite.allAlignedPredictions\n",
        "\n",
        "for prot in predictions.keys():\n",
        "  if prot != \"sequence\":\n",
        "      predictions[prot][\"sequence\"] = predictions_single_seq[\"sequence\"][prot]\n",
        "      final_results[prot] = {**predictions[prot]}\n",
        "final_results[\"statistics\"] = jsondata_list[0]\n",
        "\n",
        "json.dump(final_results, open('/content/results/predicted_biophysical_features.json', 'w'), indent=4)\n",
        "\n",
        "!zip -r /content/b2b_tools_results.zip /content/results\n",
        "\n",
        "files.download(\"/content/b2b_tools_results.zip\")"
      ],
      "metadata": {
        "cellView": "form",
        "id": "hexhk1b61Zxj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Questions? Feedback?\n",
        "\n",
        "Contact us through our [Feedback](https://www.bio2byte.be/b2btools/feedback) page on the bio2Byte website or via our email to bio2byte@vub.be.\n",
        "\n",
        "\n",
        "## PyPi repository\n",
        "\n",
        "Want to try out the b2btools? Download the Bio2Byte's tools package from our PyPi repository: https://pypi.org/project/b2bTools/.\n",
        "\n",
        "## Citations\n",
        "\n",
        "- Implementation of the b2btools to study the protein biophysical features and their conservation:\n",
        "> Kagami, L. P., Orlando, G., Raimondi, D., Ancien, F., Dixit, B., Gavaldá-García, J., Ramasamy, P., Roca-Martínez, J., Tzavella, K., & Vranken, W. (2021). b2bTools: Online predictions for protein biophysical features and their conservation. Nucleic Acids Research, 49(W1), W52–W59. https://doi.org/10.1093/nar/gkab425\n",
        "\n",
        "- DynaMine:\n",
        "> Cilia, E., Pancsa, R., Tompa, P., Lenaerts, T., & Vranken, W. F. (2013). From protein sequence to dynamics and disorder with DynaMine. Nature Communications, 4(1), 2741–2741. https://doi.org/10.1038/ncomms3741\n",
        "\n",
        "- EFoldMine:\n",
        "> Raimondi, D., Orlando, G., Pancsa, R., Khan, T., & Vranken, W. F. (2017). Exploring the Sequence-based Prediction of Folding Initiation Sites in Proteins. Scientific Reports, 7(1), 8826–8826. https://doi.org/10.1038/s41598-017-08366-3\n",
        "\n",
        "- Disomine:\n",
        "> Orlando, G., Raimondi, D., Codicè, F., Tabaro, F., & Vranken, W. (2022). Prediction of Disordered Regions in Proteins with Recurrent Neural Networks and Protein Dynamics. Journal of Molecular Biology, 434(12), 167579. https://doi.org/10.1016/j.jmb.2022.167579\n",
        "\n",
        "- AgMata:\n",
        "> Orlando, G., Silva, A., Macedo-Ribeiro, S., Raimondi, D., & Vranken, W. (2020). Accurate prediction of protein beta-aggregation with generalized statistical potentials. Bioinformatics, 36(7), 2076–2081. https://doi.org/10.1093/bioinformatics/btz912\n",
        "\n",
        "- PSPer:\n",
        "> Orlando, G., Raimondi, D., Tabaro, F., Codicè, F., Moreau, Y., & Vranken, W. F. (2019). Computational identification of prion-like RNA-binding proteins that form liquid phase-separated condensates. Bioinformatics, 35(22), 4617–4623. https://doi.org/10.1093/bioinformatics/btz274"
      ],
      "metadata": {
        "id": "y1bfRVpC3fwy"
      }
    }
  ]
}