Notebooks/json_GPCRchimeraDB_natural.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "AQCGa7gnNVZ3"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import requests, sys\n",
    "import json\n",
    "import os\n",
    "import urllib.request\n",
    "import numpy as np\n",
    "from collections import defaultdict\n",
    "from Bio import AlignIO, SeqIO\n",
    "import re\n",
    "import html\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LinearSegmentedColormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "VZnsxbIFNDzL"
   },
   "outputs": [],
   "source": [
    "# Download the needed files\n",
    "\n",
    "#Download uniprot ID of all proteins on GPCRdb and their classification\n",
    "\n",
    "#Download file with classification info\n",
    "filename_listGPCRdb = \"../data/240923_Classification_GPCRdb.xlsx\"\n",
    "listGPCRdb_df = pd.read_excel(filename_listGPCRdb)\n",
    "\n",
    "#Load the excel file of your protein entry from https://gpcrdb.org/mutational_landscape/ #check for updates?!\n",
    "\n",
    "filename_mutagenesis = \"../data/GPCRdb_variants.xlsx\"\n",
    "mutagenesis_GPCRdb_raw_data = pd.read_excel(filename_mutagenesis)\n",
    "\n",
    "#Load the excel file with Gprot and Barr coupling data from https://gproteindb.org/signprot/statistics_venn (filtered with own code)\n",
    "#this info is also on GtoP but not as complete\n",
    "filename_Gprot = \"../data/GproteinDB_table.xlsx\"\n",
    "Gprot_GPCRdb = pd.read_excel(filename_Gprot)\n",
    "filename_Barr = \"../data/BarrDB_table.xlsx\"\n",
    "Barr_GPCRdb = pd.read_excel(filename_Barr)\n",
    "\n",
    "#chimeric design info\n",
    "filename_chimeric_designs = \"../data/previous_designs.xlsx\"\n",
    "chimeric_design_df = pd.read_excel(filename_chimeric_designs)\n",
    "\n",
    "# #Download all structures on GPCRdb to access the GPCR state\n",
    "# requestURL = \"https://gpcrdb.org/services/structure/\"\n",
    "# r = requests.get(requestURL)\n",
    "\n",
    "# if not r.ok:\n",
    "#   r.raise_for_status()\n",
    "#   sys.exit()\n",
    "\n",
    "# structures_chain = json.loads(r.text)\n",
    "# json.dump(structures_chain, open(\"../data/structures_data.json\", \"w\"), indent=2)\n",
    "structures_chain= json.load(open(\"../data/structures_data.json\"))\n",
    "\n",
    "#Download all endogenous on GPCRdb to access the GPCR state\n",
    "# requestURL = \"https://gpcrdb.org/services/ligands/endogenousligands/\"\n",
    "# r = requests.get(requestURL)\n",
    "\n",
    "# if not r.ok:\n",
    "#   r.raise_for_status()\n",
    "#   sys.exit()\n",
    "\n",
    "# endogenous_ligands = json.loads(r.text)\n",
    "# json.dump(endogenous_ligands, open(\"../data/endogenous_ligands.json\", \"w\"), indent=2)\n",
    "endogenous_ligands = json.load(open(\"../data/endogenous_ligands.json\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n",
      "{'Q8IYL9': 'MNSTCIEEQHDLDHYLFPIVYIFVIIVSIPANIGSLCVSFLQAKKESELGIYLFSLSLSDLLYALTLPLWIDYTWNKDNWTFSPALCKGSAFLMYMNFYSSTAFLTCIAVDRYLAVVYPLKFFFLRTRRFALMVSLSIWILETIFNAVMLWEDETVVEYCDAEKSNFTLCYDKYPLEKWQINLNLFRTCTGYAIPLVTILICNRKVYQAVRHNKATENKEKKRIIKLLVSITVTFVLCFTPFHVMLLIRCILEHAVNFEDHSNSGKRTYTMYRITVALTSLNCVADPILYCFVTETGRYDMWNILKFCTGRCNTSQRQRKRILSVSTKDTMELEVLE', 'Q8TDS4': 'MNRHHLQDHFLEIDKKNCCVFRDDFIVKVLPPVLGLEFIFGLLGNGLALWIFCFHLKSWKSSRIFLFNLAVADFLLIICLPFLMDNYVRRWDWKFGDIPCRLMLFMLAMNRQGSIIFLTVVAVDRYFRVVHPHHALNKISNRTAAIISCLLWGITIGLTVHLLKKKMPIQNGGANLCSSFSICHTFQWHEAMFLLEFFLPLGIILFCSARIIWSLRQRQMDRHAKIKRAITFIMVVAIVFVICFLPSVVVRIRIFWLLHTSGTQNCEVYRSVDLAFFITLSFTYMNSMLDPVVYYFSSPSFPNFFSTLINRCLQRKMTGEPDNNRSTSVELTGDPNKTRGAPEALMANSGEPWSPSYLGPTSP', 'P02699': 'MNGTEGPNFYVPFSNKTGVVRSPFEAPQYYLAEPWQFSMLAAYMFLLIMLGFPINFLTLYVTVQHKKLRTPLNYILLNLAVADLFMVFGGFTTTLYTSLHGYFVFGPTGCNLEGFFATLGGEIALWSLVVLAIERYVVVCKPMSNFRFGENHAIMGVAFTWVMALACAAPPLVGWSRYIPEGMQCSCGIDYYTPHEETNNESFVIYMFVVHFIIPLIVIFFCYGQLVFTVKEAAAQQQESATTQKAEKEVTRMVIIMVIAFLICWLPYAGVAFYIFTHQGSDFGPIFMTIPAFFAKTSAVYNPVIYIMMNKQFRNCMVTTLCCGKNPLGDDEASTTVSKTETSQVAPA', 'P11617': 'MSTMGSWVYITVELAIAVLAILGNVLVCWAVWLNSNLQNVTNYFVVSLAAADIAVGVLAIPFAITISTGFCAACHNCLFFACFVLVLTQSSIFSLLAIAIDRYIAIRIPLRYNGLVTGTRAKGIIAVCWVLSFAIGLTPMLGWNNCSQPKEGRNYSQGCGEGQVACLFEDVVPMNYMVYYNFFAFVLVPLLLMLGVYLRIFLAARRQLKQMESQPLPGERARSTLQKEVHAAKSLAIIVGLFALCWLPLHIINCFTFFCPECSHAPLWLMYLTIVLSHTNSVVNPFIYAYRIREFRQTFRKIIRSHVLRRREPFKAGGTSARALAAHGSDGEQISLRLNGHPPGVWANGSAPHPERRPNGYTLGLVSGGIAPESHGDMGLPDVELLSHELKGACPESPGLEGPLAQDGAGVS', 'P30542': 'MPPSISAFQAAYIGIEVLIALVSVPGNVLVIWAVKVNQALRDATFCFIVSLAVADVAVGALVIPLAILINIGPQTYFHTCLMVACPVLILTQSSILALLAIAVDRYLRVKIPLRYKMVVTPRRAAVAIAGCWILSFVVGLTPMFGWNNLSAVERAWAANGSMGEPVIKCEFEKVISMEYMVYFNFFVWVLPPLLLMVLIYLEVFYLIRKQLNKKVSASSGDPQKYYGKELKIAKSLALILFLFALSWLPLHILNCITLFCPSCHKPSILTYIAIFLTHGNSAMNPIVYAFRIQKFRVTFLKIWNDHFRCQPAPPIDEDLPEERPDD', 'P04274': 'MGPPGNDSDFLLTTNGSHVPDHDVTEERDEAWVVGMAILMSVIVLAIVFGNVLVITAIAKFERLQTVTNYFITSLACADLVMGLAVVPFGASHILMKMWNFGNFWCEFWTSIDVLCVTASIETLCVIAVDRYIAITSPFKYQSLLTKNKARMVILMVWIVSGLTSFLPIQMHWYRATHQKAIDCYHKETCCDFFTNQAYAIASSIVSFYVPLVVMVFVYSRVFQVAKRQLQKIDKSEGRFHSPNLGQVEQDGRSGHGLRRSSKFCLKEHKALKTLGIIMGTFTLCWLPFFIVNIVHVIQDNLIPKEVYILLNWLGYVNSAFNPLIYCRSPDFRIAFQELLCLRRSSSKAYGNGYSSNSNGKTDYMGEASGCQLGQEKESERLCEDPPGTESFVNCQGTVPSLSLDSQGRNCSTNDSPL', 'P07550': 'MGQPGNGSAFLLAPNGSHAPDHDVTQERDEVWVVGMGIVMSLIVLAIVFGNVLVITAIAKFERLQTVTNYFITSLACADLVMGLAVVPFGAAHILMKMWTFGNFWCEFWTSIDVLCVTASIETLCVIAVDRYFAITSPFKYQSLLTKNKARVIILMVWIVSGLTSFLPIQMHWYRATHQEAINCYANETCCDFFTNQAYAIASSIVSFYVPLVIMVFVYSRVFQEAKRQLQKIDKSEGRFHVQNLSQVEQDGRTGHGLRRSSKFCLKEHKALKTLGIIMGTFTLCWLPFFIVNIVHVIQDNLIRKEVYILLNWIGYVNSGFNPLIYCRSPDFRIAFQELLCLRRSSLKAYGNGYSSNGNTGEQSGYHVEQEKENKLLCEDLPGTEDFVGHQGTVPSDNIDSQGRNCSTNDSLL', 'P20309': 'MTLHNNSTTSPLFPNISSSWIHSPSDAGLPPGTVTHFGSYNVSRAAGNFSSPDGTTDDPLGGHTVWQVVFIAFLTGILALVTIIGNILVIVSFKVNKQLKTVNNYFLLSLACADLIIGVISMNLFTTYIIMNRWALGNLACDLWLAIDYVASNASVMNLLVISFDRYFSITRPLTYRAKRTTKRAGVMIGLAWVISFVLWAPAILFWQYFVGKRTVPPGECFIQFLSEPTITFGTAIAAFYMPVTIMTILYWRIYKETEKRTKELAGLQASGTEAETENFVHPTGSSRSCSSYELQQQSMKRSNRRKYGRCHFWFTTKSWKPSSEQMDQDHSSSDSWNNNDAAASLENSASSDEEDIGSETRAIYSIVLKLPGHSTILNSTKLPSSDNLQVPEEELGMVDLERKADKLQAQKSVDDGGSFPKSFSKLPIQLESAVDTAKTSDVNSSVGKSTATLPLSFKEATLAKRFALKTRSQITKRKRMSLVKEKKAAQTLSAILLAFIITWTPYNIMVLVNTFCDSCIPKTFWNLGYWLCYINSTVNPVCYALCNKTFRTTFKMLLLCQCDKKKRRKQQYQQRQSVIFHKRAPEQAL'}\n"
     ]
    }
   ],
   "source": [
    "#for now select only 6 (demonstration purposes)\n",
    "of_interest=[\"\"]\n",
    "\n",
    "naturals_entry_data = \"../data/240923_cleaned_all_mammals_classA.fasta\"\n",
    "entry_uniprotID_seq = {}\n",
    "for record in SeqIO.parse(naturals_entry_data,\"fasta\"):\n",
    "    # if record.id in of_interest:\n",
    "        entry_uniprotID_seq[record.id]=str(record.seq)\n",
    "print(len(entry_uniprotID_seq))\n",
    "print(entry_uniprotID_seq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def coupling_Gprot_Barr(uniprot_id):\n",
    "    Gprot_coupling_data = []\n",
    "    Gprot_coupling_data_prot = {}\n",
    "    gprot = \"\"\n",
    "    for i in Gprot_GPCRdb[Gprot_GPCRdb['Uniprot ID'] == uniprot_id].iloc[0][1:]:\n",
    "        if not i is np.NaN:\n",
    "            gprot +=i\n",
    "            gprot += \", \"\n",
    "    Gprot_coupling_data_prot[\"value\"]=gprot[:-2]\n",
    "    Gprot_coupling_data_prot[\"reference\"]=\"https://gproteindb.org/signprot/statistics_venn\"\n",
    "    Gprot_coupling_data.append(Gprot_coupling_data_prot)\n",
    "\n",
    "    Barr_coupling_data = []\n",
    "    Barr_coupling_data_prot = {}\n",
    "    barr = \"\"\n",
    "    for i in Barr_GPCRdb[Barr_GPCRdb['Uniprot ID'] == uniprot_id].iloc[0][1:]:\n",
    "        if not i is np.NaN:\n",
    "            barr +=i\n",
    "            barr += \", \"\n",
    "    Barr_coupling_data_prot[\"value\"]=barr[:-2]\n",
    "    Barr_coupling_data_prot[\"reference\"]=\"https://arrestindb.org/signprot/arrestin_venn\"\n",
    "    Barr_coupling_data.append(Barr_coupling_data_prot)\n",
    "\n",
    "    return Gprot_coupling_data,Barr_coupling_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_mutagenesis_info_Uniprot(variations_uniprot_json,uniprot_id):\n",
    "    mutations_Uniprot = []\n",
    "    for variant in range(len(variations_uniprot_json[\"features\"])):\n",
    "        mutation = {}\n",
    "        mutation[\"start\"] = int(variations_uniprot_json[\"features\"][variant][\"begin\"])\n",
    "        mutation[\"end\"] = int(variations_uniprot_json[\"features\"][variant][\"end\"])\n",
    "        mutation_type = variations_uniprot_json[\"features\"][variant][\"consequenceType\"]\n",
    "        mutation[\"type\"] = mutation_type\n",
    "\n",
    "        typical_AA = \"AVLIPMCFYWSTQNHKRDEG\"\n",
    "        mutation[\"original residue\"] = variations_uniprot_json[\"features\"][variant][\"wildType\"]\n",
    "        if not mutation[\"original residue\"] in typical_AA:\n",
    "            continue\n",
    "        predictions = []\n",
    "        prediction = {}\n",
    "        consensus = []\n",
    "        if mutation_type.lower() == \"missense\":\n",
    "            mutation[\"alternative residue\"] = variations_uniprot_json[\"features\"][variant][\"mutatedType\"]\n",
    "            try:\n",
    "                nb_predictors = len(variations_uniprot_json[\"features\"][variant][\"predictions\"])\n",
    "                for predictor in range(nb_predictors):\n",
    "                    prediction = {}\n",
    "                    algorithm = variations_uniprot_json[\"features\"][variant][\"predictions\"][predictor][\"predAlgorithmNameType\"].lower()\n",
    "                    score = variations_uniprot_json[\"features\"][variant][\"predictions\"][predictor][\"score\"]\n",
    "                    prediction_value = variations_uniprot_json[\"features\"][variant][\"predictions\"][predictor][\"predictionValType\"]\n",
    "                    if algorithm.lower() == \"polyphen\":\n",
    "                        prediction[\"predictor\"] = \"polyphen\"\n",
    "                        prediction[\"value\"] = str(score)\n",
    "                        if score >= 0.2:\n",
    "                            prediction[\"prediction\"] = \"probably damaging\"\n",
    "                            consensus.append(\"-\")\n",
    "                        elif score >= 0.1:\n",
    "                            prediction[\"prediction\"] = \"possibly damaging\"\n",
    "                            consensus.append(\"-\")\n",
    "                        else:\n",
    "                            prediction[\"prediction\"] = \"benign\"\n",
    "                            consensus.append(\"+\")\n",
    "                        predictions.append(prediction)\n",
    "                    elif algorithm.lower() == \"sift\":\n",
    "                        prediction[\"predictor\"] = \"SIFT\"\n",
    "                        prediction[\"value\"] = str(score)\n",
    "                        if score <= 0.05:\n",
    "                            prediction[\"prediction\"] = \"deleterious\"\n",
    "                            consensus.append(\"-\")\n",
    "                        else:\n",
    "                            prediction[\"prediction\"] = \"tolerated\"\n",
    "                            consensus.append(\"+\")\n",
    "                        predictions.append(prediction)\n",
    "            except:\n",
    "                prediction[\"predictor\"] = \"\"\n",
    "                prediction[\"value\"] = \"\"\n",
    "                try:\n",
    "                    effect = variations_uniprot_json[\"features\"][variant][\"clinicalSignificances\"][0][\"type\"]\n",
    "                    if \"benign\" in effect.lower() or \"tolerated\" in effect.lower():\n",
    "                        prediction[\"prediction\"] = \"tolerated\"\n",
    "                        consensus.append(\"+\")\n",
    "                    elif \"deleterious\" in effect.lower() or \"damaging\" in effect.lower():\n",
    "                        prediction[\"prediction\"] = \"deleterious\"\n",
    "                        consensus.append(\"-\")\n",
    "                    else:\n",
    "                        prediction[\"prediction\"] = \"not indicated\"\n",
    "                except:\n",
    "                    prediction[\"prediction\"] = \"not indicated\"\n",
    "                predictions.append(prediction)\n",
    "\n",
    "        elif mutation_type.lower() == \"stop gained\":\n",
    "            mutation[\"alternative residue\"] = \"termination\"\n",
    "            prediction[\"predictor\"] = \"\"\n",
    "            prediction[\"value\"] = \"\"\n",
    "            try:\n",
    "                effect = variations_uniprot_json[\"features\"][variant][\"clinicalSignificances\"][0][\"type\"]\n",
    "                if \"benign\" in effect.lower() or \"tolerated\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"tolerated\"\n",
    "                    consensus.append(\"+\")\n",
    "                elif \"deleterious\" in effect.lower() or \"damaging\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"deleterious\"\n",
    "                    consensus.append(\"-\")\n",
    "                else:\n",
    "                    prediction[\"prediction\"] = \"not indicated\"\n",
    "            except:\n",
    "                prediction[\"prediction\"] = \"not indicated\"\n",
    "            predictions.append(prediction)\n",
    "\n",
    "        elif mutation_type.lower() == \"inframe deletion\":\n",
    "            mutation[\"alternative residue\"] = \"missing\"\n",
    "            prediction[\"predictor\"] = \"\"\n",
    "            prediction[\"value\"] = \"\"\n",
    "            try:\n",
    "                effect = variations_uniprot_json[\"features\"][variant][\"clinicalSignificances\"][0][\"type\"]\n",
    "                if \"benign\" in effect.lower() or \"tolerated\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"tolerated\"\n",
    "                    consensus.append(\"+\")\n",
    "                elif \"deleterious\" in effect.lower() or \"damaging\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"deleterious\"\n",
    "                    consensus.append(\"-\")\n",
    "                else:\n",
    "                    prediction[\"prediction\"] = \"not indicated\"\n",
    "            except:\n",
    "                prediction[\"prediction\"] = \"not indicated\"\n",
    "            predictions.append(prediction)\n",
    "\n",
    "        elif mutation_type.lower() == \"frameshift\":\n",
    "            mutation[\"alternative residue\"] = \"\"\n",
    "            prediction[\"predictor\"] = \"\"\n",
    "            prediction[\"value\"] = \"\"\n",
    "            try:\n",
    "                effect = variations_uniprot_json[\"features\"][variant][\"clinicalSignificances\"][0][\"type\"]\n",
    "                if \"benign\" in effect.lower() or \"tolerated\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"tolerated\"\n",
    "                    consensus.append(\"+\")\n",
    "                elif \"deleterious\" in effect.lower() or \"damaging\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"deleterious\"\n",
    "                    consensus.append(\"-\")\n",
    "                else:\n",
    "                        prediction[\"prediction\"] = \"not indicated\"\n",
    "            except:\n",
    "                prediction[\"prediction\"] = \"not indicated\"\n",
    "            predictions.append(prediction)\n",
    "\n",
    "        elif mutation_type.lower() == \"stop lost\":\n",
    "            mutation[\"alternative residue\"] = variations_uniprot_json[\"features\"][variant][\"mutatedType\"]\n",
    "            prediction[\"predictor\"] = \"\"\n",
    "            prediction[\"value\"] = \"\"\n",
    "            try:\n",
    "                effect = variations_uniprot_json[\"features\"][variant][\"clinicalSignificances\"][0][\"type\"]\n",
    "                if \"benign\" in effect.lower() or \"tolerated\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"tolerated\"\n",
    "                    consensus.append(\"+\")\n",
    "                elif \"deleterious\" in effect.lower() or \"damaging\" in effect.lower():\n",
    "                    prediction[\"prediction\"] = \"deleterious\"\n",
    "                    consensus.append(\"-\")\n",
    "                else:\n",
    "                    prediction[\"prediction\"] = \"not indicated\"\n",
    "            except:\n",
    "                prediction[\"prediction\"] = \"not indicated\"\n",
    "            predictions.append(prediction)\n",
    "\n",
    "        else: #skip if it not a missense, frameshift, stop gained, stop lost, inframe deletion\n",
    "            continue\n",
    "\n",
    "        mutation[\"effect(s)\"] = predictions\n",
    "\n",
    "        if len(consensus) == 0:\n",
    "                consensus_ = \"no consensus\"\n",
    "        else:\n",
    "            for element in consensus:\n",
    "                if len(list(set(consensus))) == 1:\n",
    "                    if list(set(consensus))[0] == \"+\":\n",
    "                        consensus_ = \"tolerated\"\n",
    "                    elif list(set(consensus))[0] == \"-\":\n",
    "                        consensus_ = \"deleterious\"\n",
    "                else:\n",
    "                    consensus_ = \"no consensus\"\n",
    "\n",
    "        mutation[\"consensus effect\"] = consensus_\n",
    "\n",
    "        mutation[\"reference\"] = f\"https://www.uniprot.org/uniprotkb/{uniprot_id}/variant-viewer\"\n",
    "        mutations_Uniprot.append(mutation)\n",
    "\n",
    "    return mutations_Uniprot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_pharmaco_info_GPCRdb(pharmaco_data_json):\n",
    "    pharmaco = []\n",
    "    for dic in pharmaco_data_json:\n",
    "        pharmaco_dic = {}\n",
    "        pharmaco_dic[\"start\"]= dic[\"mutation_pos\"]\n",
    "        pharmaco_dic[\"end\"]= dic[\"mutation_pos\"]\n",
    "        pharmaco_dic[\"original residue\"]= dic[\"mutation_from\"]\n",
    "        pharmaco_dic[\"alternative residue\"]= dic[\"mutation_to\"]\n",
    "        pharmaco_dic[\"studied parameter\"]= dic[\"exp_type\"]\n",
    "        DB = dic[\"ligand_id\"]\n",
    "        pharmaco_dic[\"ligand\"]= dic[\"ligand_name\"]\n",
    "        if \"CHEMBL\" in DB: #chembl\n",
    "            pharmaco_dic[\"link ligand\"] = f\"https://www.ebi.ac.uk/chembl/compound_report_card/{DB}/\"\n",
    "        elif DB.isnumeric(): #pubchem\n",
    "            pharmaco_dic[\"link ligand\"] = f\"https://pubchem.ncbi.nlm.nih.gov/compound/{DB}\"\n",
    "        else:\n",
    "            pharmaco_dic[\"link ligand\"] = \"\"\n",
    "        pharmaco_dic[\"ligand type\"]= \"\"\n",
    "        effect_value = round(dic[\"exp_fold_change\"],1)\n",
    "        if effect_value < 0:\n",
    "            impact = \"increase\"\n",
    "            effect_value = effect_value*-1\n",
    "        else:\n",
    "            impact = \"decrease\"\n",
    "\n",
    "        if effect_value == 0.0:\n",
    "            pharmaco_dic[\"effect\"]=\"None\"\n",
    "        else:\n",
    "            pharmaco_dic[\"effect\"]= str(effect_value) + \" fold \" + impact\n",
    "\n",
    "        pharmaco_dic[\"reference\"]= dic[\"reference\"]\n",
    "\n",
    "        pharmaco.append(pharmaco_dic)\n",
    "    return pharmaco"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_sq_atom_distance(i, j):\n",
    "    \"\"\"Squared euclidean distance between two 3d points\"\"\"\n",
    "    return (i[0] - j[0]) * (i[0] - j[0]) + \\\n",
    "            (i[1] - j[1]) * (i[1] - j[1]) + \\\n",
    "            (i[2] - j[2]) * (i[2] - j[2])\n",
    "\n",
    "def identify_gaps(pdb_file, chain_pdb, offset, end): #code modified from pdb_gap.py file from pdbtools Copyright 2018 João Pedro Rodrigues\n",
    "    fhandle = open(pdb_file, 'r')\n",
    "    centroid = ' CA '  # respect spacing. 'CA  ' != ' CA '\n",
    "    distance_threshold = 4.0 * 4.0\n",
    "    prev_at = (None, None, None, None, (None, None, None))\n",
    "    model = 0\n",
    "    n_gaps = 0\n",
    "    gap = []\n",
    "    for line in fhandle:\n",
    "\n",
    "        if line.startswith('MODEL'):\n",
    "            model = int(line[10:14])\n",
    "\n",
    "        elif line.startswith('ATOM'):\n",
    "            atom_name = line[12:16]\n",
    "            if atom_name != centroid:\n",
    "                continue\n",
    "\n",
    "            resn = line[17:20]\n",
    "            resi = int(line[22:26])\n",
    "            chain = line[21]\n",
    "            x = float(line[30:38])\n",
    "            y = float(line[38:46])\n",
    "            z = float(line[46:54])\n",
    "\n",
    "            at_uid = (model, chain, resi, resn, atom_name, (x, y, z))\n",
    "            if prev_at[0] == at_uid[0] and prev_at[1] == at_uid[1]:\n",
    "                d = calculate_sq_atom_distance(at_uid[5], prev_at[5])\n",
    "                if d > distance_threshold:\n",
    "                    gap.append([prev_at[1],prev_at[2],at_uid[1],at_uid[2]])\n",
    "                    n_gaps += 1\n",
    "                elif prev_at[2] + 1 != at_uid[2]:\n",
    "                    gap.append([prev_at[1],prev_at[2],at_uid[1],at_uid[2]])\n",
    "                    n_gaps += 1\n",
    "\n",
    "            prev_at = at_uid\n",
    "\n",
    "    gaps_cleaned = []\n",
    "    start = offset\n",
    "    for section in gap:\n",
    "        if section[0] == chain_pdb and section[2] == chain_pdb:\n",
    "            stop = section[1]\n",
    "            if start < 1000 and stop < 1000:\n",
    "                gaps_cleaned.append([start,stop])\n",
    "            start = section[3]\n",
    "    gaps_cleaned.append([start,end])\n",
    "    return gaps_cleaned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def merge_duplicates(dicts):\n",
    "    # Step 1: Group dictionaries by the 'start' key\n",
    "    grouped = defaultdict(list)\n",
    "    for d in dicts:\n",
    "        grouped[d['start']].append(d)\n",
    "    \n",
    "    result = []\n",
    "    conflicts = []\n",
    "\n",
    "    # Step 2: Process each group\n",
    "    for start, items in grouped.items():\n",
    "        if len(items) > 1:\n",
    "            # Check if all 'type' values are the same\n",
    "            types = set(d['type'] for d in items)\n",
    "            if len(types) == 1:\n",
    "                # Merge 'reference' values\n",
    "                merged_references = \"PDBePISA \"\n",
    "                for d in items:\n",
    "                    merged_references += d['description'].split('.')[0] + \" \"\n",
    "                    merged_references += d['description'].split(' ')[-1] + \" \"\n",
    "                # Create a new dictionary with merged references\n",
    "                new_dict = items[0].copy()\n",
    "                new_dict['description'] = merged_references[:-1]\n",
    "                new_dict['reference'] = \"https://www.ebi.ac.uk/pdbe/pisa/\"\n",
    "                result.append(new_dict)\n",
    "            else:\n",
    "                # Print dictionaries with different 'type' values\n",
    "                for d in items:\n",
    "                    conflicts.append(d)\n",
    "        else:\n",
    "            result.append(items[0])\n",
    "    \n",
    "    return result, conflicts\n",
    "\n",
    "def retrieve_interacting_residues_PDB(pdb_id, chain_pdb,uniprot_pdb_start,pdb_start):\n",
    "    # retrieve the interacting residues in the PDBs from PISA, need to make sure it doesn't take into accound the interactions between 2 sym GPCRs\n",
    "    # interacting residues is defined by a bsa > 0\n",
    "    #https://github.com/PDBe-KB/pdbe-pisa-json/blob/main/PISA-APIs.ipynb\n",
    "\n",
    "    # interacting_residues_list = []\n",
    "    # binders_chain= []\n",
    "    # for pdb_id, chain_pdb,uniprot_pdb_start,pdb_start in zip(pdb_ids,chain_pdbs,uniprot_pdb_starts,pdb_starts):\n",
    "        try: #when its just 1 chain or 1 chain and a ligand PISA doesn't work\n",
    "            difference_pdb_uniprot_start = pdb_start-uniprot_pdb_start\n",
    "            interacting_residues = []\n",
    "            response = requests.get(f\"https://www.ebi.ac.uk/pdbe/api/pisa/assembly/{pdb_id.lower()}/1\")\n",
    "            interface_count = response.json()[pdb_id.lower()][\"assembly\"][\"interface_count\"]\n",
    "            for i in range(1,interface_count+1):\n",
    "                response_single_interface = requests.get(f\"https://www.ebi.ac.uk/pdbe/api/pisa/interface/{pdb_id.lower()}/1/{i}/\")\n",
    "                data = response_single_interface.json()\n",
    "                if \"/\" in chain_pdb:\n",
    "                    chain_pdb = chain_pdb.split(\"/\")\n",
    "                for j in range(len(data[\"molecules\"])):\n",
    "                    if isinstance(chain_pdb,str):\n",
    "                        if data[\"molecules\"][j][\"chain_id\"]==chain_pdb:\n",
    "                            for bsa,position in zip(data[\"molecules\"][j][\"buried_surface_areas\"],data[\"molecules\"][j]['residue_seq_ids']):\n",
    "                                if bsa >0.0:\n",
    "                                    interacting_residues.append(int(position)-difference_pdb_uniprot_start)\n",
    "                            if j == 0: #there is supposed to be only 2 molecules, the GPCR and the interacting molecule\n",
    "                                chain_interacting_molecule = data[\"molecules\"][1][\"chain_id\"]\n",
    "                            else: \n",
    "                                chain_interacting_molecule = data[\"molecules\"][0][\"chain_id\"]\n",
    "                    elif isinstance(chain_pdb,list):\n",
    "                        if chain_pdb[0] in data[\"molecules\"][j][\"chain_id\"] and chain_pdb[1] in data[\"molecules\"][j+1][\"chain_id\"]:\n",
    "                            break\n",
    "                        else:\n",
    "                            if chain_pdb[0] in data[\"molecules\"][j][\"chain_id\"] or chain_pdb[1] in data[\"molecules\"][j][\"chain_id\"]:\n",
    "                                for bsa,position in zip(data[\"molecules\"][j][\"buried_surface_areas\"],data[\"molecules\"][j]['residue_seq_ids']):\n",
    "                                    if bsa >0.0:\n",
    "                                        interacting_residues.append(int(position)-difference_pdb_uniprot_start)\n",
    "                            if j == 0: #there is supposed to be only 2 molecules, the GPCR and the interacting molecule\n",
    "                                chain_interacting_molecule = data[\"molecules\"][1][\"chain_id\"]\n",
    "                            else: \n",
    "                                chain_interacting_molecule = data[\"molecules\"][0][\"chain_id\"]\n",
    "        except:\n",
    "            pass\n",
    "\n",
    "        return interacting_residues,chain_interacting_molecule\n",
    "\n",
    "def extract_name_binders(chain_of_interest,pdb_file_path):\n",
    "\n",
    "    molecule_name = None\n",
    "    current_molecule = \"\"\n",
    "    reading_molecule = False\n",
    "    found_chain = False\n",
    "    \n",
    "    # Open and read the PDB file\n",
    "    with open(pdb_file_path, 'r') as file:\n",
    "        lines = file.readlines()\n",
    "    \n",
    "    for line in lines:\n",
    "        if line.startswith(\"COMPND\"):\n",
    "            # Start reading the molecule name if \"MOLECULE\" is in the line\n",
    "            if \"MOLECULE\" in line:\n",
    "                reading_molecule = True\n",
    "                current_molecule = line.split(\":\")[1].strip().rstrip(\";\")  # Extract initial part of the molecule name\n",
    "            # If the molecule name is being read and it continues on the next line\n",
    "            elif reading_molecule and \"CHAIN\" not in line:\n",
    "                pattern = r\"\\d+\\s+(.+)\"\n",
    "                match = re.search(pattern, line)\n",
    "                current_molecule += \" \"+match.group(1).strip().rstrip(\";\")\n",
    "            # Once we reach the chain of interest\n",
    "            if f\"CHAIN: {chain_of_interest}\" in line:\n",
    "                found_chain = True\n",
    "            # If molecule and chain have been found, stop reading\n",
    "            if found_chain and current_molecule and \";\" in line:\n",
    "                molecule_name = current_molecule\n",
    "                break\n",
    "    return molecule_name\n",
    " \n",
    "    \n",
    "def translate_interacting_residues_IC_EC(interacting_residues,binders,pdb_list,all_regions):\n",
    "    interactions_pdb_list_EC = []\n",
    "    interactions_pdb_list_IC = []\n",
    "    for residues,binder,pdb_id in zip(interacting_residues,binders,pdb_list):\n",
    "        for residue in residues:\n",
    "            for region in all_regions:\n",
    "                if residue >= region[\"start\"] and residue <= region[\"end\"]:\n",
    "                    if region[\"name\"] in [\"Nterm\",\"ECL1\",\"ECL2\",\"ECL3\"]:\n",
    "                        region_residue = \"EC\"\n",
    "                    elif region[\"name\"] in [\"Cterm\",\"ICL1\",\"ICL2\",\"ICL3\"]:\n",
    "                        region_residue = \"IC\"\n",
    "                    elif region[\"name\"] in [\"TM1\",\"TM3\",\"TM5\",\"TM7\"]:\n",
    "                        if residue <= region[\"start\"] + round(((region[\"end\"]-region[\"start\"])/2)):\n",
    "                            region_residue = \"EC\"\n",
    "                        else:\n",
    "                            region_residue = \"IC\"\n",
    "                    elif region[\"name\"] in [\"TM2\",\"TM4\",\"TM6\"]:\n",
    "                        if residue <= region[\"start\"] + round(((region[\"end\"]-region[\"start\"])/2)):\n",
    "                            region_residue = \"IC\"\n",
    "                        else:\n",
    "                            region_residue = \"EC\"\n",
    "\n",
    "                    interactions_pdb = {}\n",
    "                    interactions_pdb[\"start\"]=int(residue)\n",
    "                    interactions_pdb[\"end\"]=int(residue)\n",
    "                    if region_residue == \"IC\":\n",
    "                        interactions_pdb[\"type\"]=\"Intracellular binding pocket residue\"\n",
    "                    else:\n",
    "                        interactions_pdb[\"type\"]=\"Extracellular binding pocket residue\"\n",
    "                    \n",
    "                    interactions_pdb[\"description\"]=f\"{binder}. Inferred from {pdb_id}\"\n",
    "                    interactions_pdb[\"reference\"]=\"https://www.ebi.ac.uk/pdbe/pisa/\"\n",
    "\n",
    "                    if region_residue == \"IC\":\n",
    "                        interactions_pdb_list_IC.append(interactions_pdb)\n",
    "                    else:\n",
    "                        interactions_pdb_list_EC.append(interactions_pdb)\n",
    "\n",
    "    #merge duplicates\n",
    "    interactions_pdb_list_IC_no_duplicates,_ = merge_duplicates(interactions_pdb_list_IC)\n",
    "    interactions_pdb_list_EC_no_duplicates,_ = merge_duplicates(interactions_pdb_list_EC)\n",
    "\n",
    "    return interactions_pdb_list_IC_no_duplicates, interactions_pdb_list_EC_no_duplicates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pad_line(line):\n",
    "    \"\"\"Helper function to pad line to 80 characters in case it is shorter\"\"\"\n",
    "    size_of_line = len(line)\n",
    "    if size_of_line < 80:\n",
    "        padding = 80 - size_of_line + 1\n",
    "        line = line.strip('\\n') + ' ' * padding + '\\n'\n",
    "    return line[:81]  # 80 + newline character\n",
    "\n",
    "def check_first_residue_pdb_uniprot(pdb_file,chain_interest,sequence):\n",
    "    #SIFT mapping identifies 1st residue from PDB in uniprot. But that 1st residue is not necessarly the first residue that is resolved in the PDB\n",
    "    #Here check if the first residue that is resolved corresponds to a residue at the same position in sequence\n",
    "    #if yes, no need to renumber\n",
    "    Three_to_One_AA = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',\n",
    "     'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N', \n",
    "     'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W', \n",
    "     'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}\n",
    "\n",
    "    fhandle = open(pdb_file, 'r')\n",
    "    _pad_line = pad_line\n",
    "    for line in fhandle:\n",
    "        line = _pad_line(line)\n",
    "        if line.startswith(\"ATOM\") and line[21]==chain_interest:\n",
    "            residue = Three_to_One_AA[line[17:20].upper()]\n",
    "            pos = int(line[24:27])\n",
    "            break\n",
    "    if sequence[pos-1] == residue:\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "        \n",
    "def SIFT_mapping(pdb_id,uniprot_id):\n",
    "    #find equivalence start uniprot and pdb \n",
    "    # https://www.ebi.ac.uk/pdbe/api/doc/sifts.html\n",
    "    response = requests.get(f\"https://www.ebi.ac.uk/pdbe/api/mappings/uniprot/{pdb_id}\")\n",
    "    data = response.json()\n",
    "    uniprot_pdb_start = data[pdb_id.lower()]['UniProt'][uniprot_id]['mappings'][0]['unp_start']\n",
    "    pdb_start = data[pdb_id.lower()]['UniProt'][uniprot_id]['mappings'][0]['start']['author_residue_number']\n",
    "    if not pdb_start:\n",
    "        pdb_start = data[pdb_id.lower()]['UniProt'][uniprot_id]['mappings'][0]['start']['residue_number']\n",
    "    return uniprot_pdb_start, pdb_start\n",
    "\n",
    "def renumber_pdb_uniprot_start(pdb_file,pdb_id,chain,uniprot_id,sequence):\n",
    "    uniprot_pdb_start, pdb_start=SIFT_mapping(pdb_id,uniprot_id)\n",
    "    new_pdb_file_path = f\"../examples/3Dstructures/pdb_renumbered/{pdb_id}.pdb\"\n",
    "    modify = False\n",
    "    if not check_first_residue_pdb_uniprot(pdb_file,chain,sequence) and uniprot_pdb_start != pdb_start :\n",
    "    #if both SIFT and the manually matching with the sequence doesn't work, renumber\n",
    "        modify = True\n",
    "    if modify:\n",
    "        !python pdb_reres.py -$uniprot_pdb_start -$chain $pdb_file > $new_pdb_file_path\n",
    "    else:\n",
    "        uniprot_pdb_start, pdb_start = 0, 0\n",
    "    return new_pdb_file_path, uniprot_pdb_start, pdb_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def retrieve_pdb_dsbonds_interactions(uniprot_json,uniprot_id,sequence):\n",
    "    structures = []\n",
    "    ds_bonds = []\n",
    "    pdbs = []\n",
    "    chains = []\n",
    "    uniprot_pdb_starts = []\n",
    "    pdb_starts = []\n",
    "    interacting_residues_list = []\n",
    "    binders_list= []\n",
    "    for i in range(len(uniprot_json['uniProtKBCrossReferences'])):\n",
    "        if uniprot_json['uniProtKBCrossReferences'][i]['database'] == 'PDB':\n",
    "            pdb_id = uniprot_json['uniProtKBCrossReferences'][i]['id']\n",
    "            length_chain = 0\n",
    "            if ',' in uniprot_json['uniProtKBCrossReferences'][i][\"properties\"][2][\"value\"]: #when there is a \",\" it means that there are multiple fragments, let's assume it's longer than 200 residues then\n",
    "                sections = uniprot_json['uniProtKBCrossReferences'][i][\"properties\"][2][\"value\"].split(\",\")\n",
    "            else:\n",
    "                sections = [uniprot_json['uniProtKBCrossReferences'][i][\"properties\"][2][\"value\"]]\n",
    "            for j in range(len(sections)):\n",
    "                range_chain=sections[j].split(\"=\")[1].split(\"-\")\n",
    "                length_chain += int(range_chain[1])-int(range_chain[0])\n",
    "            if length_chain > 200:\n",
    "                #download pdb file\n",
    "                try:\n",
    "                    pdb_file_path = f'../data/tmp/{pdb_id}.pdb'\n",
    "                    urllib.request.urlretrieve(f'https://files.rcsb.org/download/{pdb_id}.pdb', pdb_file_path)\n",
    "\n",
    "                    #get conformational state structure => provided by GPCRdb\n",
    "                    state = \"Undetermined\" #default\n",
    "                    for structure in structures_chain:\n",
    "                        if structure[\"pdb_code\"] == pdb_id:\n",
    "                            state = structure[\"state\"]\n",
    "                            break\n",
    "                    chain_pdb =  uniprot_json['uniProtKBCrossReferences'][i][\"properties\"][2][\"value\"][0]\n",
    "                    full_chain_pdb = uniprot_json['uniProtKBCrossReferences'][i][\"properties\"][2][\"value\"].split(\"=\")[0]\n",
    "\n",
    "                    #get SIFTS mapping PDB seq > UniProt mapping and renumber pdb. New pdb being written\n",
    "                    new_pdb_file_path, uniprot_pdb_start, pdb_start=renumber_pdb_uniprot_start(pdb_file_path,pdb_id,chain_pdb,uniprot_id,sequence)\n",
    "                    if os.path.exists(new_pdb_file_path):\n",
    "                        path_DB_pdb = f\"file:///examples/3Dstructures/pdb_renumbered/{pdb_id}.pdb\"\n",
    "                    else:\n",
    "                        path_DB_pdb = f\"https://files.rcsb.org/download/{pdb_id}.pdb\"\n",
    "                        new_pdb_file_path = pdb_file_path\n",
    "                    #this is needed to find the interactions within the pdb file\n",
    "                    pdbs.append(pdb_id)\n",
    "                    chains.append(full_chain_pdb)\n",
    "                    uniprot_pdb_starts.append(uniprot_pdb_start)\n",
    "                    pdb_starts.append(pdb_start)\n",
    "\n",
    "                    #get offset and gaps in structure\n",
    "                    offset = 0\n",
    "                    start_found = False\n",
    "                    gaps = []\n",
    "                    gap = False\n",
    "                    with open(new_pdb_file_path, 'r') as f:\n",
    "                        for line in f:\n",
    "                            line_list = line.split()\n",
    "                            if not start_found:\n",
    "                                if line_list[0]==\"ATOM\" and line_list[4]==chain_pdb:\n",
    "                                    if int(line_list[5]) >= 1:\n",
    "                                        offset = int(line_list[5]) - 1 #line_list[5] gives the 1st position so the offset is line_list[5] -1\n",
    "                                        start_found = True\n",
    "                            else:\n",
    "                                if line_list[0]==\"ATOM\" and line_list[4]==chain_pdb:\n",
    "                                    end = int(line_list[5])\n",
    "                            \n",
    "                    gaps = identify_gaps (new_pdb_file_path, chain_pdb, offset+1, end)\n",
    "                    structures.append({\"value\":pdb_id,\"chain\": chain_pdb, \"state\":state, \"offset\":  offset, \"gaps\": gaps, \"url\":path_DB_pdb, \"reference\":f\"https://www.rcsb.org/structure/{pdb_id.upper()}\"})\n",
    "\n",
    "                    #Look for disulfide bridges in PDB files\n",
    "                    with open(new_pdb_file_path, 'r') as f:\n",
    "                        for line in f:\n",
    "                            line_list = line.split()\n",
    "                            if line_list[0]==\"SSBOND\":\n",
    "                                if line_list[3] == chain_pdb and line_list[6] == chain_pdb:\n",
    "                                    ssbond = {}\n",
    "                                    ssbond[\"start\"] = int(line_list[4])\n",
    "                                    ssbond[\"end\"]= int(line_list[7])\n",
    "                                    ssbond[\"description\"] = \"Disulfide bond\"\n",
    "                                    ssbond[\"reference\"] = \"Extracted from PDB files\"\n",
    "                                    ds_bonds.append(ssbond)\n",
    "\n",
    "\n",
    "                    #find interacting residues at ligand binding site and G protein binding site\n",
    "                    interacting_residues,binder_chain = retrieve_interacting_residues_PDB(pdb_id,chain_pdb,uniprot_pdb_start,pdb_start)\n",
    "                    interacting_residues_list.append(list(set(interacting_residues)))\n",
    "                    binder = extract_name_binders(binder_chain,pdb_file_path)\n",
    "                    if binder != None:\n",
    "                        binders_list.append(binder)\n",
    "\n",
    "                    #remove pdb file\n",
    "                    os.remove(pdb_file_path)\n",
    "                \n",
    "                except:\n",
    "                    continue\n",
    "            else:\n",
    "                with open(\"to_look_at_pdbs.txt\",\"a\") as f:\n",
    "                    f.write(pdb_id+\"\\n\")\n",
    "                    f.close()\n",
    "    \n",
    "\n",
    "    #Keep only the DS bonds that are in every PDB file (sometimes add DS bond to stabilize structure to be crystalized)\n",
    "    ds_bonds_noduplicates=[dict(t) for t in {tuple(d.items()) for d in ds_bonds}]\n",
    "    for single in ds_bonds_noduplicates:\n",
    "        count = 0\n",
    "        for multiple in ds_bonds:\n",
    "            if single == multiple:\n",
    "                count +=1\n",
    "        if count != len(structures):\n",
    "            ds_bonds_noduplicates.remove(single)   \n",
    "            \n",
    "    return structures, ds_bonds_noduplicates, interacting_residues_list, binders_list, pdbs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#equivalence between positions in sequence and in MSA\n",
    "#dictionary with list of list. In every sublist, 2 elements, 1st is the position in sequence, the 2nd the position in MSA\n",
    "def map_seq_MSA(sequence_aligned):\n",
    "    previous = 0\n",
    "    translate = {}\n",
    "    sequence_nogaps = sequence_aligned.replace(\"-\",\"\")\n",
    "    for res in range(len(sequence_nogaps)):\n",
    "        idx_msa = previous + sequence_aligned[previous:].index(sequence_nogaps[res])\n",
    "        translate[res+1]=idx_msa+1\n",
    "        previous = idx_msa + 1\n",
    "    return translate\n",
    "\n",
    "#Microswitches/motifs - identify them based on their defined columns in mammalian MSA\n",
    "def motifs_microswitches_literature(MSA,uniprot_id):\n",
    "    #GPCRdb finds sodium pockets\n",
    "    #As microswitches are well defined in literature we can check ourselves if these well knwon microswitches are present in our gpcrs\n",
    "    #All known microswitches in literature for class A\n",
    "    E_DRY_W = {\"positions\":[\"3.49\", \"3.50\", \"3.51\"],\"residues\":[\"ED\", \"R\", \"WY\"], \"name\": \"E/DRY/W motif (ionic lock switch)\"}\n",
    "    CWxP = {\"positions\":[\"6.47\", \"6.48\", \"6.50\"], \"residues\":[\"C\", \"W\", \"P\"], \"name\": \"CWxP motif (transmission toggle switch)\"}\n",
    "    NPxxY = {\"positions\":[\"7.49\", \"7.50\", \"7.53\"], \"residues\":[\"N\",\"P\",\"Y\"], \"name\": \"NPxxY motif (tyr toggle switch)\"}\n",
    "    PIF = {\"positions\": [\"5.50\", \"3.40\", \"6.44\"], \"residues\":[\"P\",\"I\",\"F\"], \"name\": \"PIF motif\"}\n",
    "    hydrophobic_lock = {\"positions\":[\"3.43\",\"6.40\"], \"residues\":[\"LVIM\", \"LVIM\"], \"name\": \"hydrophobic lock\"}\n",
    "    ionic_lock = {\"positions\":[\"6.30\"], \"residues\":[\"DE\"], \"name\": \"ionic lock\"}\n",
    "    #disulfide bond between TM3 and ECL2 is already identified by Uniprot in the \"Disulfide bonds\" section\n",
    "    #Sodium binding pocket (allosteric action): middle of the 7TMs. Identified by GPCRdb but are the identified ones all of them???\n",
    "\n",
    "    #the positions are in human readable format (not pyton - starts at 0)\n",
    "    MSA_E_DRY_W = {\"positions\":[694,695,696],\"residues\":[\"ED\", \"R\", \"WY\"], \"name\": \"E/DRY/W motif (ionic lock switch)\"}\n",
    "    MSA_CWxP = {\"positions\":[1211,1212,1214], \"residues\":[\"C\", \"W\", \"P\"], \"name\": \"CWxP motif (transmission toggle switch)\"}\n",
    "    MSA_NPxxY = {\"positions\":[1290,1291,1294], \"residues\":[\"N\",\"P\",\"Y\"], \"name\": \"NPxxY motif (tyr toggle switch)\"}\n",
    "    MSA_PIF = {\"positions\": [947,685,1204], \"residues\":[\"P\",\"I\",\"F\"], \"name\": \"PIF motif\"}\n",
    "    MSA_hydrophobic_lock = {\"positions\":[688,1200], \"residues\":[\"LVIM\", \"LVIM\"], \"name\": \"Hydrophobic lock\"}\n",
    "    MSA_ionic_lock = {\"positions\":[1190], \"residues\":[\"DE\"], \"name\": \"Ionic lock\"}\n",
    "    MSA_sodium_pocket = {\"positions\":[619,684], \"residues\":[\"D\",\"S\"], \"name\": \"Sodium binding pocket\"}\n",
    "\n",
    "    TM1x50={\"positions\":[579],\"residues\":[\"N\"], \"name\": \"1.50 (BW numbering)\"}\n",
    "    TM2x50={\"positions\":[619],\"residues\":[\"D\"], \"name\": \"2.50 (BW numbering)\"}\n",
    "    TM3x50={\"positions\":[695],\"residues\":[\"R\"], \"name\": \"3.50 (BW numbering)\"}\n",
    "    TM4x50={\"positions\":[751],\"residues\":[\"P\"], \"name\": \"4.50 (BW numbering)\"}\n",
    "    TM5x50={\"positions\":[947],\"residues\":[\"P\"], \"name\": \"5.50 (BW numbering)\"}\n",
    "    TM6x50={\"positions\":[1214],\"residues\":[\"P\"], \"name\": \"6.50 (BW numbering)\"}\n",
    "    TM7x50={\"positions\":[1291],\"residues\":[\"P\"], \"name\": \"7.50 (BW numbering)\"}\n",
    "\n",
    "    alignment = AlignIO.read(open(MSA), \"fasta\")\n",
    "    len_MSA=alignment.get_alignment_length()\n",
    "    record_dict = SeqIO.index(MSA, \"fasta\")\n",
    "    aligned_seq_interest = str(record_dict[uniprot_id].seq)\n",
    "    translate_seq_MSA = map_seq_MSA(aligned_seq_interest) #gives position of a unaligned res in msa\n",
    "    translate_MSA_seq = {v: k for k, v in translate_seq_MSA.items()} #gives position of a aligned res in unaligned seq\n",
    "    microswitch_types = [MSA_E_DRY_W, MSA_CWxP, MSA_NPxxY, MSA_PIF, MSA_hydrophobic_lock, MSA_ionic_lock,\n",
    "                         TM1x50,TM2x50,TM3x50,TM4x50,TM5x50,TM6x50,TM7x50]\n",
    "    microswitches = []\n",
    "    microswitches_residues = []\n",
    "\n",
    "    for microswitch_type in microswitch_types:\n",
    "        are_there = []\n",
    "        for position, residue in zip(microswitch_type[\"positions\"], microswitch_type[\"residues\"]):\n",
    "            if aligned_seq_interest[position-1] in residue:\n",
    "                are_there.append(True)\n",
    "            else:\n",
    "                are_there.append(False)\n",
    "        for i, (position, residue) in enumerate(zip(microswitch_type[\"positions\"], microswitch_type[\"residues\"])):\n",
    "            microswitch_residue = {}\n",
    "            \n",
    "            #take into account the possibility that there is a gap at that position in the MSA\n",
    "            if position in translate_MSA_seq:\n",
    "                microswitch_residue[\"start\"] = translate_MSA_seq[position]\n",
    "                microswitch_residue[\"end\"] = translate_MSA_seq[position]\n",
    "                residue_motif = aligned_seq_interest[position-1]\n",
    "            else:\n",
    "                for next in range(position+1,len_MSA):\n",
    "                    if next in translate_MSA_seq:\n",
    "                        microswitch_residue[\"start\"] = translate_MSA_seq[next]\n",
    "                        microswitch_residue[\"end\"] = translate_MSA_seq[next]\n",
    "                        residue_motif = aligned_seq_interest[next-1]\n",
    "                        break\n",
    "\n",
    "            if not all(are_there) and not are_there[i]:\n",
    "                if residue_motif == \"F\" and microswitch_type[\"name\"]==\"PIF motif\":\n",
    "                    microswitch_residue[\"description\"] = residue_motif + \" instead of \"+ residue+ \" from \" + \" (part of \" + microswitch_type[\"name\"]+ \")\"\n",
    "                elif residue_motif == \"R\" and microswitch_type[\"name\"]==\"E/DRY/W motif (ionic lock switch)\":\n",
    "                    microswitch_residue[\"description\"] = residue_motif + \" instead of \"+ residue + \" (part of \" + microswitch_type[\"name\"]+ \")\"\n",
    "                elif \"(BW numbering)\" in microswitch_type[\"name\"]:\n",
    "                    microswitch_residue[\"description\"] = residue_motif+ \" instead of \"+ residue + \" \" + microswitch_type[\"name\"]\n",
    "                else:\n",
    "                    microswitch_residue[\"description\"] = residue_motif + \" instead of \"+ residue + \" (part of \" + microswitch_type[\"name\"]+ \")\"\n",
    "            else:\n",
    "                if residue_motif == \"F\" and microswitch_type[\"name\"]==\"PIF motif\":\n",
    "                    microswitch_residue[\"description\"] = residue_motif + \" part of \" + microswitch_type[\"name\"] +\" and hydrophobic lock\"\n",
    "                elif residue_motif == \"R\" and microswitch_type[\"name\"]==\"E/DRY/W motif (ionic lock switch)\":\n",
    "                    microswitch_residue[\"description\"] = residue_motif+ \" part of \" + microswitch_type[\"name\"] +\" and ionic lock\"\n",
    "                elif \"(BW numbering)\" in microswitch_type[\"name\"]:\n",
    "                    microswitch_residue[\"description\"] = residue_motif+ \" \" + microswitch_type[\"name\"]\n",
    "                else:\n",
    "                    microswitch_residue[\"description\"] = residue_motif+ \" part of \" + microswitch_type[\"name\"]\n",
    "            if are_there[i]:\n",
    "                microswitch_residue[\"conserved\"] = \"yes\"\n",
    "            else:\n",
    "                microswitch_residue[\"conserved\"] = \"no\"\n",
    "            microswitch_residue[\"reference\"] = \"Based on alignment\"\n",
    "            microswitches_residues.append(microswitch_residue)\n",
    "\n",
    "    return microswitches_residues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#features from Uniprot binding site, PTM, natural variants\n",
    "def features_uniprot(uniprot_json,uniprot_id):\n",
    "    binding_sites = []\n",
    "    PTMs = []\n",
    "    disulfide_bonds = []\n",
    "    mutagenesiss = []\n",
    "\n",
    "    for i in range(len(uniprot_json['features'])):\n",
    "\n",
    "        #Motifs/Microswitches Uniprot are already identified by the self written code above where we check if all well defined microswitches knwon in literature are present/absent in the GPCR of interest\n",
    "\n",
    "        #Binding site (orthosteric & allosteric)\n",
    "        if uniprot_json['features'][i]['type'] == 'Binding site':\n",
    "            binding_site = {}\n",
    "            binding_site['start'] = uniprot_json['features'][i]['location']['start']['value']\n",
    "            binding_site['end'] = uniprot_json['features'][i]['location']['end']['value']\n",
    "            binding_site['type'] = \"Ligand binding residue\"\n",
    "            binding_site['description'] = uniprot_json['features'][i]['ligand']['name']\n",
    "            binding_site[\"reference\"] = f\"https://www.uniprot.org/uniprotkb/{uniprot_id}/entry\"\n",
    "            binding_sites.append(binding_site)\n",
    "        ##PTMs\n",
    "        #Glycosylation\n",
    "        elif uniprot_json['features'][i]['type'] == 'Glycosylation' or uniprot_json['features'][i]['type'] == 'Lipidation' or uniprot_json['features'][i]['type'] == 'Modified residue':\n",
    "            PTM = {}\n",
    "            if uniprot_json['features'][i]['type'] == 'Glycosylation':\n",
    "                PTM['start'] = uniprot_json['features'][i]['location']['start']['value']\n",
    "                PTM['end'] = uniprot_json['features'][i]['location']['end']['value']\n",
    "                PTM['description'] = 'Glycosylation'\n",
    "            elif uniprot_json['features'][i]['type'] == 'Lipidation':\n",
    "                PTM['start'] = uniprot_json['features'][i]['location']['start']['value']\n",
    "                PTM['end'] = uniprot_json['features'][i]['location']['end']['value']\n",
    "                PTM['description'] = 'Lipidation'\n",
    "            else:\n",
    "                ptm_types = [\"phospho\",\"methyl\",\"acetyl\", \"amid\", \"pyrro\", \"hydroxy\", \"l-\", \"d-\", \"sulf\",\"nitro\"]\n",
    "                full_ptm_types = [\"Phosphorylation\", \"Methylation\", \"Acetylation\", \"Amidation\", \"Pyrrolidone carboxylic acid\", \"Hydroxylation\", \"Isomerization\", \"Isomerization\", \"Sulfation\", \"Nitrosylation\"]\n",
    "                description = uniprot_json['features'][i]['description']\n",
    "                for idx, ptm_type in enumerate(ptm_types):\n",
    "                    if ptm_type in uniprot_json['features'][i]['description'].lower():\n",
    "                        description = full_ptm_types[idx]\n",
    "                        break\n",
    "                PTM['start'] = uniprot_json['features'][i]['location']['start']['value']\n",
    "                PTM['end'] = uniprot_json['features'][i]['location']['end']['value']\n",
    "                PTM['description'] = description\n",
    "            PTM['reference'] = f\"https://www.uniprot.org/uniprotkb/{uniprot_id}/entry\"\n",
    "            PTMs.append(PTM)\n",
    "\n",
    "    return binding_sites,PTMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#retrieve PTMs on Scop3P and compare with PTMs we have already found\n",
    "def retrieve_PTM_Scop3P(uniprot_id):\n",
    "    requestURL = f\"https://iomics.ugent.be/scop3p/api/modifications?accession={uniprot_id}\"\n",
    "    try: \n",
    "        r = requests.get(requestURL)\n",
    "        if not r.ok:\n",
    "            r.raise_for_status()\n",
    "            sys.exit()\n",
    "        scop3P_PTM = json.loads(r.text)\n",
    "    except:\n",
    "        scop3P_PTM = []\n",
    "\n",
    "    if len(scop3P_PTM)>0:\n",
    "\n",
    "        #add if not alrady there\n",
    "        ptms = []\n",
    "        for ptm in scop3P_PTM[\"modifications\"]:\n",
    "            position = ptm['position']\n",
    "            # if not position in positions_previous:\n",
    "            ptm_dict = {'start': position , 'end': position, 'description': 'Phosphorylation', 'reference': \"https://iomics.ugent.be/scop3p/index?protein={uniprot_id}\"}\n",
    "            ptms.append(ptm_dict)\n",
    "    else:\n",
    "        ptms=[]\n",
    "    return ptms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find TM regions based on TMbed predictions\n",
    "#find the TM regions\n",
    "def TMbed_retrieve_allregions(TM_regions_file,uniprotID):\n",
    "    \n",
    "    regions = []\n",
    "    region_names = ['Nterm','TM1', 'ICL1', 'TM2', 'ECL1', 'TM3', 'ICL2', 'TM4', 'ECL2', 'TM5', 'ICL3', 'TM6', 'ECL3', 'TM7','Cterm']\n",
    "\n",
    "    with open(TM_regions_file, 'r') as f:\n",
    "        count = 0\n",
    "        TM_counter = 1\n",
    "        name =\"\"\n",
    "        ok=False\n",
    "        error = False\n",
    "        for line in f:\n",
    "            line=line.replace('\\n',\"\").replace('\\t',\" \")\n",
    "            line=list(line.split())\n",
    "            count +=1\n",
    "            if \">\" in line[0]:\n",
    "                try:\n",
    "                    name = line[0].split('|')[1]\n",
    "                except:\n",
    "                    name = line[0][1:]\n",
    "                count = 1\n",
    "            if name == uniprotID and (count%3==0):\n",
    "                counter_region = 0\n",
    "                start = False\n",
    "                \n",
    "                previous = 1\n",
    "                for idx,element in enumerate(line[0]):\n",
    "                    if start == False:\n",
    "                        if element == \"H\" or element == \"h\":\n",
    "                            start = True\n",
    "                            lower_lim_TM = idx+1 #human readble limits\n",
    "                            region = {}\n",
    "                            region[\"name\"] = region_names[counter_region]\n",
    "                            region[\"start\"] = previous\n",
    "                            region[\"end\"] = lower_lim_TM-1\n",
    "                            region[\"reference\"] = \"TMbed\"\n",
    "                            regions.append(region)\n",
    "                            counter_region +=1\n",
    "                    else:\n",
    "                        if element != \"H\" and element != \"h\":\n",
    "                            upper_lim_TM = idx  #human readble limits\n",
    "                            start = False\n",
    "                            region = {}\n",
    "                            region[\"name\"] = region_names[counter_region]\n",
    "                            region[\"start\"] = lower_lim_TM\n",
    "                            region[\"end\"] = upper_lim_TM\n",
    "                            region[\"reference\"] = \"https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-022-04873-x\"\n",
    "                            regions.append(region)\n",
    "                            counter_region +=1\n",
    "                            previous=upper_lim_TM+1\n",
    "                region = {}\n",
    "                region[\"name\"] = region_names[counter_region]\n",
    "                region[\"start\"] = previous\n",
    "                region[\"end\"] = len(line[0])\n",
    "                region[\"reference\"] = \"https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-022-04873-x\"\n",
    "                regions.append(region)\n",
    "                break\n",
    "\n",
    "    return regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutting_pts_2_ss_region(dict_regions,all_regions=None):\n",
    "    translated_regions = []\n",
    "    if isinstance(dict_regions,dict):\n",
    "        for region in dict_regions.keys():\n",
    "            lower_lim = dict_regions[region][0]\n",
    "            upper_lim = dict_regions[region][1]\n",
    "            for ss_region in all_regions:\n",
    "                if lower_lim >= ss_region[\"start\"] and lower_lim <= ss_region[\"end\"]:\n",
    "                    ss_lower_lim = ss_region[\"name\"]\n",
    "                if upper_lim >= ss_region[\"start\"] and upper_lim <= ss_region[\"end\"]:\n",
    "                    ss_upper_lim = ss_region[\"name\"]\n",
    "            translated_regions.append(str(lower_lim)+\"-\"+str(upper_lim)+f\" ({ss_lower_lim}-{ss_upper_lim})\")\n",
    "    elif isinstance(dict_regions,list):\n",
    "        if all_regions:\n",
    "            for section in dict_regions:\n",
    "                info = []\n",
    "                for pos in section:\n",
    "                    for ss_region in all_regions:\n",
    "                        if int(pos) >= ss_region[\"start\"] and int(pos) <= ss_region[\"end\"]:\n",
    "                            ss_pos = ss_region[\"name\"]\n",
    "                            break\n",
    "                    info.append(pos)\n",
    "                    info.append(ss_pos)\n",
    "                translated_regions.append(str(info[0])+\"-\"+str(info[2])+f\" ({info[1]}-{info[3]})\")\n",
    "        else:\n",
    "            for region in dict_regions:\n",
    "                lower_lim = region[0]\n",
    "                upper_lim = region[1]\n",
    "                translated_regions.append(str(lower_lim)+\"-\"+str(upper_lim))\n",
    "    return translated_regions\n",
    "\n",
    "def often_used_cutting_pts(cutting_points_parent,chimera,sequence,all_regions_cutting_pts,related_chimeras,dict_regions):\n",
    "    cutting_pts = chimera[cutting_points_parent]\n",
    "    for region in cutting_pts:\n",
    "        positions=region.split(\" \")[0]\n",
    "        name=region.split(\" \")[1][1:-1]\n",
    "        if positions.split(\"-\")[0] != \"1\":\n",
    "            position = int(positions.split(\"-\")[0])\n",
    "            name_region= name.split(\"-\")[0]\n",
    "            lower_lim = [d for d in dict_regions if d.get(\"name\") == name_region][0][\"start\"]\n",
    "            upper_lim = [d for d in dict_regions if d.get(\"name\") == name_region][0][\"end\"]\n",
    "            a_third = round((upper_lim-lower_lim)/3)\n",
    "            if position < (lower_lim + a_third):\n",
    "                idx = 0\n",
    "            elif position < (lower_lim + 2*a_third):\n",
    "                idx = 1\n",
    "            else:\n",
    "                idx = 2\n",
    "            try:\n",
    "                all_regions_cutting_pts[name_region].append([sequence[position-1]+str(position),idx])\n",
    "            except:\n",
    "                all_regions_cutting_pts[name_region]=[[sequence[position-1]+str(position),idx]]\n",
    "            try:\n",
    "                related_chimeras[sequence[position-1]+str(position)].append(chimera[\"name\"])\n",
    "            except:\n",
    "                related_chimeras[sequence[position-1]+str(position)]=[chimera[\"name\"]]\n",
    "        if positions.split(\"-\")[1] != str(len(sequence)):\n",
    "            position = int(positions.split(\"-\")[1])\n",
    "            name_region= name.split(\"-\")[1]\n",
    "            lower_lim = [d for d in dict_regions if d.get(\"name\") == name_region][0][\"start\"]\n",
    "            upper_lim = [d for d in dict_regions if d.get(\"name\") == name_region][0][\"end\"]\n",
    "            a_third = lower_lim+round((upper_lim-lower_lim)/3)\n",
    "            if position < (lower_lim + a_third):\n",
    "                idx = 0\n",
    "            elif position < (lower_lim + 2*a_third):\n",
    "                idx = 1\n",
    "            else:\n",
    "                idx = 2\n",
    "            try:\n",
    "                all_regions_cutting_pts[name_region].append([sequence[position-1]+str(position),idx])\n",
    "            except:\n",
    "                all_regions_cutting_pts[name_region]=[[sequence[position-1]+str(position),idx]]\n",
    "            try:\n",
    "                related_chimeras[sequence[position-1]+str(position)].append(chimera[\"name\"])\n",
    "            except:\n",
    "                related_chimeras[sequence[position-1]+str(position)]=[chimera[\"name\"]]\n",
    "    return all_regions_cutting_pts,related_chimeras\n",
    "\n",
    "\n",
    "def retrieve_involvement_natural_chimeric_design(uniprot_id,abb_name,sequence,chimeric_design_df):\n",
    "\n",
    "    involvement = []\n",
    "\n",
    "    all_regions_cutting_pts_all = {}\n",
    "    related_chimeras_all = {}\n",
    "    for parent_column_id in ['Reference_id','Target_id']:\n",
    "\n",
    "        #find rows that have uniprot as ref id or target id\n",
    "        designs_parent = chimeric_design_df[chimeric_design_df[parent_column_id] == uniprot_id]\n",
    "\n",
    "        #Info from rows\n",
    "        names_chimeras = designs_parent['Chimera_name'].tolist()\n",
    "        ids_chimeras = designs_parent['Chimera_name_ids'].tolist()\n",
    "        regions_chimera = designs_parent['Chimera_parts'].tolist()\n",
    "        name_target_chimeras = designs_parent['Target_name'].tolist()\n",
    "        id_target_chimeras = designs_parent['Target_id'].tolist()\n",
    "\n",
    "        name_ref_chimeras = designs_parent['Reference_name'].tolist()\n",
    "        id_ref_chimeras = designs_parent['Reference_id'].tolist()\n",
    "\n",
    "        regions_ref_chimeras = designs_parent['Reference_cutting_points'].tolist()\n",
    "        regions_target_chimeras = designs_parent['Target_cutting_points'].tolist()\n",
    "\n",
    "        expression = designs_parent['Expression binary'].tolist()\n",
    "        fct = designs_parent['Function binary'].tolist()\n",
    "\n",
    "        application = designs_parent['Application'].tolist()\n",
    "        type_chimera = designs_parent['Chimera Type (1/2/3)'].tolist()\n",
    "        Gprot = designs_parent['G-protein'].tolist()\n",
    "        Ligand = designs_parent['Ligand'].tolist()\n",
    "        structures = designs_parent['3D structure PDB'].tolist()\n",
    "        biblio = designs_parent['DOI'].tolist()\n",
    "\n",
    "        for i,name in enumerate(names_chimeras):\n",
    "\n",
    "            all_regions = TMbed_retrieve_allregions(\"../data/TM_regions_chimeras_TMbed.txt\",name)\n",
    "            cutting_pt_chimera = cutting_pts_2_ss_region(eval(regions_chimera[i]),all_regions) \n",
    "\n",
    "            all_regions_ref = TMbed_retrieve_allregions(\"../data/Uniprot_TMbed.txt\",id_ref_chimeras[i])\n",
    "            cutting_pt_ref = cutting_pts_2_ss_region(eval(regions_ref_chimeras[i]),all_regions_ref)\n",
    "\n",
    "            all_regions_target = TMbed_retrieve_allregions(\"../data/Uniprot_TMbed.txt\",id_target_chimeras[i])\n",
    "            cutting_pt_target = cutting_pts_2_ss_region(eval(regions_target_chimeras[i]),all_regions_target)\n",
    "\n",
    "            pharma_name_ref = html.unescape(get_pharma_name(id_ref_chimeras[i],name_ref_chimeras[i]))\n",
    "            pharma_name_target = html.unescape(get_pharma_name(id_target_chimeras[i], name_target_chimeras[i]))\n",
    "\n",
    "            pharma_name_ref_ = pharma_name_ref\n",
    "            pharma_name_target_ = pharma_name_target\n",
    "            if \"receptor\" in pharma_name_ref.lower():\n",
    "                pharma_name_ref_ = pharma_name_ref.replace(\" receptor\",\"\")\n",
    "            if \"receptor\" in pharma_name_target.lower():\n",
    "                pharma_name_target_ = pharma_name_target.replace(\" receptor\",\"\")\n",
    "            pharma_name = pharma_name_ref_ + \" \" + pharma_name_target_ + \" receptor\"\n",
    "            if \"adrenoceptor\" in pharma_name:\n",
    "                pharma_name = pharma_name.replace(\" receptor\",\"\")\n",
    "\n",
    "            if isinstance(structures[i],str):\n",
    "                pdb = structures[i]\n",
    "            else:\n",
    "                pdb = \"\"\n",
    "\n",
    "            if isinstance(Gprot[i],str):\n",
    "                gprot=Gprot[i]\n",
    "            else:\n",
    "                gprot=\"\"\n",
    "\n",
    "            if isinstance(Ligand[i],str):\n",
    "                ligand=Ligand[i]\n",
    "            else:\n",
    "                ligand=\"\"\n",
    "\n",
    "            chimera={\n",
    "            \"name\":name,\n",
    "            \"name_pharma\":pharma_name,\n",
    "            \"id\":ids_chimeras[i],\n",
    "            \"ref\": name_ref_chimeras[i],\n",
    "            \"ref_pharma_name\": pharma_name_ref,\n",
    "            \"target\": name_target_chimeras[i],\n",
    "            \"target_pharma_name\":pharma_name_target,\n",
    "            \"cutting_point_chimera\": cutting_pt_chimera,\n",
    "            \"cutting_point_ref\": cutting_pt_ref,\n",
    "            \"cutting_point_target\": cutting_pt_target,\n",
    "            \"expression_function\": fct[i],\n",
    "            \"type\":type_chimera[i],\n",
    "            \"GprotLigand\": gprot+\" \"+ligand,\n",
    "            \"application\": application[i]+\" \"+pdb,\n",
    "            \"reference\": biblio[i]\n",
    "            }\n",
    "\n",
    "            involvement.append(chimera)\n",
    "\n",
    "            if \"_\".join(name.split(\"_\")[:2]) == abb_name:\n",
    "                all_regions_cutting_pts_all,related_chimeras_all=often_used_cutting_pts(\"cutting_point_ref\",chimera,sequence,all_regions_cutting_pts_all,related_chimeras_all,all_regions_ref)\n",
    "            else: \n",
    "                all_regions_cutting_pts_all,related_chimeras_all=often_used_cutting_pts(\"cutting_point_target\",chimera,sequence,all_regions_cutting_pts_all,related_chimeras_all,all_regions_target)                \n",
    "    \n",
    "    for key,value in all_regions_cutting_pts_all.items():\n",
    "        unique_items = list(map(list, set(map(tuple, value))))\n",
    "        all_regions_cutting_pts_all[key] = unique_items\n",
    "\n",
    "    for key,value in related_chimeras_all.items():\n",
    "        unique_items = list(set(value))\n",
    "        related_chimeras_all[key] = unique_items\n",
    "\n",
    "    return involvement,all_regions_cutting_pts_all,related_chimeras_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pharma_name(uniprotID,abb_name):\n",
    "    #GtoP or gpcrdb_name or pharmacological name\n",
    "    try:\n",
    "        requestURL = f\"https://gpcrdb.org/services/protein/accession/{uniprotID}\"\n",
    "\n",
    "        r = requests.get(requestURL, headers={ \"Accept\" : \"application/json\"})\n",
    "        if not r.ok:\n",
    "            info_entry = None\n",
    "        else:\n",
    "            info_entry = json.loads(r.text)\n",
    "    except:\n",
    "        info_entry = None\n",
    "\n",
    "    if not info_entry is None:\n",
    "        clean_html_tags = re.compile('<.*?>')\n",
    "        pharma_name = re.sub(clean_html_tags, '', info_entry[\"name\"])\n",
    "    else:\n",
    "        pharma_name = abb_name\n",
    "\n",
    "    return pharma_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_duplicates(list_dictionaries,uniprot_id, descriminator1, descriminator2=None):\n",
    "    # Initialize a dictionary to count occurrences of (start, description) pairs\n",
    "    count = {}\n",
    "\n",
    "    # First pass: Count occurrences of each (start, description) pair\n",
    "    for d in list_dictionaries:\n",
    "        if descriminator2 != None:\n",
    "            identifier = (d[descriminator1], d[descriminator2])\n",
    "        else:\n",
    "            identifier = (d[descriminator1])\n",
    "        if identifier in count:\n",
    "            count[identifier] += 1\n",
    "        else:\n",
    "            count[identifier] = 1\n",
    "\n",
    "    # Initialize a set to track seen (start, description) pairs\n",
    "    seen = set()\n",
    "    # Initialize a list to store the filtered dictionaries\n",
    "    unique_dict_list = []\n",
    "\n",
    "    # Second pass: Filter dictionaries and update \"other_key\" for duplicates\n",
    "    for d in list_dictionaries:\n",
    "        if descriminator2 != None:\n",
    "            identifier = (d[descriminator1], d[descriminator2])\n",
    "        else:\n",
    "            identifier = (d[descriminator1])\n",
    "        if identifier not in seen:\n",
    "            # If this (start, description) pair is a duplicate (appears more than once)\n",
    "            if count[identifier] > 1:\n",
    "                # don't choose randomly the reference, set it to the preferred reference type\n",
    "                if descriminator2 != None:\n",
    "                    d[\"reference\"] = f\"https://iomics.ugent.be/scop3p/index?protein={uniprot_id}\"\n",
    "                else:\n",
    "                    d[\"reference\"] = \"https://www.ebi.ac.uk/pdbe/pisa/\"\n",
    "            # Add it to the seen set\n",
    "            seen.add(identifier)\n",
    "            # Add the dictionary to the unique list\n",
    "            unique_dict_list.append(d)\n",
    "    return unique_dict_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heatmap_cutting_pts(pharma_name,abbreviated_name,chimeras,allregions):\n",
    "\n",
    "    # Number of heatmaps\n",
    "    num_heatmaps = len(allregions)\n",
    "\n",
    "    # Create a custom colormap starting with white and transitioning to red\n",
    "    colors = [\"white\", \"#ffcccc\", \"#ff9999\", \"#ff6666\", \"#ff3333\", \"#ff0000\", \"#cc0000\", \"#990000\", \"#660000\", \"#330000\"]\n",
    "    cmap_reds_white = LinearSegmentedColormap.from_list(\"RedsWhite\", colors)\n",
    "\n",
    "    for functionality in [True,False]:\n",
    "        for parent_nameDB,parent_name,parent_pts in zip([\"EC\",\"IC\"],[\"ref\",\"target\"],[\"cutting_point_ref\",\"cutting_point_target\"]):\n",
    "            heatmaps_data = []\n",
    "            for i,region in enumerate(allregions):\n",
    "                heatmap = [0]* (region[\"end\"]-region[\"start\"]+1)\n",
    "                nb_designs = 0\n",
    "                for design in chimeras:\n",
    "                    if functionality:\n",
    "                        if design[\"expression_function\"].lower() == \"yes\" or design[\"expression_function\"].lower() == \"not assessed\":\n",
    "                            if design[parent_name] == abbreviated_name:\n",
    "                                for cutting_pt in design[parent_pts]:\n",
    "                                    if cutting_pt.split(\" \")[1].split(\"-\")[0][1:] == region[\"name\"]:\n",
    "                                        position = int(cutting_pt.split(\" \")[0].split(\"-\")[0]) - region[\"start\"]+ 1 \n",
    "                                        heatmap[position-1] +=1\n",
    "                                    if cutting_pt.split(\" \")[1].split(\"-\")[1][:-1] == region[\"name\"]:\n",
    "                                        position = int(cutting_pt.split(\" \")[0].split(\"-\")[1]) - region[\"start\"] + 1 \n",
    "                                        heatmap[position-1] +=1\n",
    "                                nb_designs +=1\n",
    "                    else:\n",
    "                        if design[\"expression_function\"].lower() != \"yes\" and design[\"expression_function\"].lower() != \"not assessed\":\n",
    "                            if design[parent_name] == abbreviated_name:\n",
    "                                for cutting_pt in design[parent_pts]:\n",
    "                                    if cutting_pt.split(\" \")[1].split(\"-\")[0][1:] == region[\"name\"]:\n",
    "                                        position = int(cutting_pt.split(\" \")[0].split(\"-\")[0]) - region[\"start\"]+ 1 \n",
    "                                        heatmap[position-1] +=1\n",
    "                                    if cutting_pt.split(\" \")[1].split(\"-\")[1][:-1] == region[\"name\"]:\n",
    "                                        position = int(cutting_pt.split(\" \")[0].split(\"-\")[1]) - region[\"start\"] + 1 \n",
    "                                        heatmap[position-1] +=1\n",
    "                                nb_designs +=1\n",
    "                heatmaps_data.append((np.array(heatmap)/nb_designs).reshape(-1, 1))\n",
    "\n",
    "            # Create a figure with a grid of subplots arranged horizontally\n",
    "            fig, axes = plt.subplots(1, num_heatmaps, figsize=(10,5))\n",
    "            if functionality:\n",
    "                fig.suptitle(f\"Functional Designs ({nb_designs}) with {html.unescape(pharma_name)} ({abbreviated_name}) as {parent_nameDB} side parent\",fontsize=12)\n",
    "            else:\n",
    "                fig.suptitle(f\"Non-functional Designs ({nb_designs}) with {html.unescape(pharma_name)} ({abbreviated_name}) as {parent_nameDB} side parent\",fontsize=12)\n",
    "            \n",
    "            # Plot each heatmap\n",
    "            for i, ax in enumerate(axes):\n",
    "                cax = ax.matshow(heatmaps_data[i], aspect='auto', cmap=cmap_reds_white, vmin=0, vmax=1)  # Ensure vmin and vmax are set to 0 and 1\n",
    "                ax.yaxis.set_ticks_position('left')\n",
    "                ax.yaxis.set_tick_params(labelleft=True)\n",
    "\n",
    "                # Major ticks every 4 elements\n",
    "                major_tick_positions = range(0, heatmaps_data[i].shape[0], 5)\n",
    "                major_tick_labels = np.arange(allregions[i][\"start\"], allregions[i][\"end\"] + 1)[::5]\n",
    "                ax.set_yticks(major_tick_positions)\n",
    "                ax.set_yticklabels(major_tick_labels, fontsize=6)\n",
    "                \n",
    "                # Minor ticks at every element\n",
    "                minor_tick_positions = range(heatmaps_data[i].shape[0])\n",
    "                ax.set_yticks(minor_tick_positions, minor=True)\n",
    "                ax.set_title(allregions[i][\"name\"],fontsize=10)\n",
    "                ax.set_xticks([])  # Hide x-axis ticks\n",
    "                # for y in range(heatmaps_data[i].shape[0]):\n",
    "                #     ax.axhline(y - 0.5, color='black', linewidth=0.5)  # Line between each row\n",
    "\n",
    "            # Adjust layout to make room for colorbar\n",
    "            plt.subplots_adjust(right=0.85, wspace=1.2)\n",
    "\n",
    "            # Add color bar completely to the right\n",
    "            cbar_ax = fig.add_axes([0.87, 0.15, 0.02, 0.7])\n",
    "            cbar = fig.colorbar(cax, cax=cbar_ax, orientation='vertical')\n",
    "\n",
    "            # Add label to the color bar\n",
    "            cbar.set_label('Frequency of designs with these cutting points', fontsize=10)\n",
    "            # Add ticks and labels to the color bar\n",
    "            # cbar_ticks = [0, 0.25, 0.5, 0.75, 1.0]\n",
    "            # cbar.set_ticks(cbar_ticks)\n",
    "            # cbar.set_ticklabels(['0%', '25%', '50%', '75%', '100%'])\n",
    "\n",
    "            if functionality:\n",
    "                plt.savefig(f\"../examples/heatmap_cutting_pts/{abbreviated_name}_{parent_nameDB}_functional.png\",dpi=200)\n",
    "            else:\n",
    "                plt.savefig(f\"../examples/heatmap_cutting_pts/{abbreviated_name}_{parent_nameDB}_non_functional.png\",dpi=200)\n",
    "\n",
    "            plt.close()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert scientific name UniProt to common name\n",
    "def parse_species_file(file_path):\n",
    "    species_dict = {}\n",
    "    \n",
    "    with open(file_path, 'r') as file:\n",
    "        lines = file.readlines()\n",
    "    \n",
    "    scientific_name = None\n",
    "    \n",
    "    for line in lines:\n",
    "        if \"N=\" in line:\n",
    "            scientific_name = line.split(\"N=\")[1].strip()\n",
    "        elif \"C=\" in line and scientific_name:\n",
    "            common_name = line.split(\"C=\")[1].strip()\n",
    "            species_dict[scientific_name] = common_name\n",
    "            scientific_name = None  # Reset for the next entry\n",
    "    \n",
    "    return species_dict\n",
    "\n",
    "file_path = '../data/UniProt_names_scientific_common.txt' \n",
    "species_dict = parse_species_file(file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q8IYL9\n",
      "Q8TDS4\n",
      "P02699\n",
      "P11617\n",
      "P30542\n",
      "P04274\n",
      "P07550\n",
      "P20309\n"
     ]
    }
   ],
   "source": [
    "for uniprotID, prot_seq in entry_uniprotID_seq.items():\n",
    "    requestURL = f\"https://rest.uniprot.org/uniprotkb/{uniprotID}.json\"\n",
    "    r = requests.get(requestURL, headers={ \"Accept\" : \"application/json\"})\n",
    "\n",
    "    if not r.ok:\n",
    "        r.raise_for_status()\n",
    "        sys.exit()\n",
    "\n",
    "    uniprot_json = json.loads(r.text)\n",
    "\n",
    "    class_ = 'A' #always class A for now\n",
    "    abbreviated_name = uniprot_json[\"uniProtkbId\"]\n",
    "    names = []\n",
    "    names.append({\"value\":uniprot_json['proteinDescription']['recommendedName']['fullName']['value'], \"reference\":\"UniProt\"})\n",
    "    try:\n",
    "        if 'alternativeNames' in uniprot_json['proteinDescription']:\n",
    "            for i in range(len(uniprot_json['proteinDescription']['alternativeNames'])):\n",
    "                names.append({\"value\":uniprot_json['proteinDescription']['alternativeNames'][i]['fullName']['value'],\"reference\":\"UniProt\"})\n",
    "    except:\n",
    "        pass\n",
    "\n",
    "    # if \"a\"==\"a\":\n",
    "    if not os.path.exists(f'../examples/json_entries/all_mammals_latest/{abbreviated_name.upper()}.json'):\n",
    "        print(uniprotID)\n",
    "        species_scientific = uniprot_json['organism']['scientificName']\n",
    "        if species_scientific in species_dict:\n",
    "            species = species_dict[species_scientific]\n",
    "        \n",
    "        #find classfication based on human classification in GPCRdb\n",
    "        #find human ortholog\n",
    "        family = \"\"\n",
    "        subclass_ligand = \"\"\n",
    "        subclass_phylo = \"\"\n",
    "        try:\n",
    "            if listGPCRdb_df[listGPCRdb_df['Uniprot ID'] == uniprotID][\"Phylogenetically-based\"].values[0] == \"A-other\":\n",
    "                if not abbreviated_name.endswith(\"HUMAN\"):\n",
    "                    abbreviated_name_human = (abbreviated_name.split('_')[0]+\"_\"+\"HUMAN\").lower()\n",
    "                    uniprot_id_human = listGPCRdb_df[listGPCRdb_df['Name'] == abbreviated_name_human]['Uniprot ID'].values[0] \n",
    "            else:\n",
    "                uniprot_id_human = uniprotID\n",
    "            family = listGPCRdb_df[listGPCRdb_df['Uniprot ID'] == uniprot_id_human][\"Subclass\"].values[0].rstrip() #need to change this to interpro API for those not on GPCRdb\n",
    "            if \"Class A\" in family:\n",
    "                family = family.replace(\"Class A \",\"\")\n",
    "                family = family[0].upper()+family[1:]\n",
    "            if \"receptors\" in family:\n",
    "                family = family.replace(\"receptors\",\"\").rstrip()\n",
    "            elif \"receptor\" in family:\n",
    "                family = family.replace(\"receptor\",\"\").rstrip()\n",
    "            subclass_ligand = listGPCRdb_df[listGPCRdb_df['Uniprot ID'] == uniprot_id_human][\"Ligand-based\"].values[0].rstrip() #need to change so it works for all mammals (put same as what we have for humans? What with those not on GPCRdb?)\n",
    "            if \"receptors\" in subclass_ligand:\n",
    "                subclass_ligand = subclass_ligand.replace(\"receptors\",\"\").rstrip()\n",
    "            elif \"receptor\" in subclass_ligand:\n",
    "                subclass_ligand = subclass_ligand.replace(\"receptor\",\"\").rstrip()\n",
    "            subclass_phylo = listGPCRdb_df[listGPCRdb_df['Uniprot ID'] == uniprot_id_human][\"Phylogenetically-based\"].values[0]\n",
    "            if \"A-\" in subclass_phylo:\n",
    "                subclass_phylo = subclass_phylo.split('-')[1].rstrip()\n",
    "                subclass_phylo = subclass_phylo[0].upper()+subclass_phylo[1:]\n",
    "        except:\n",
    "            if any(\"olfactory\" in entry[\"value\"].lower() for entry in names): #many of the olfactory humans are not in GPCRdb so no classification. We took the olf human ones from GROSS\n",
    "                family = \"Olfactory\"\n",
    "                subclass_ligand = \"Olfactory\"\n",
    "                subclass_phylo = \"Olfactory\"\n",
    "            else:\n",
    "                print(\"Problem with classification (family,subfamily...)\")\n",
    "\n",
    "        #GtoP or gpcrdb_name or pharmacological name\n",
    "        try:\n",
    "            requestURL = f\"https://gpcrdb.org/services/protein/accession/{uniprotID}\"\n",
    "\n",
    "            r = requests.get(requestURL, headers={ \"Accept\" : \"application/json\"})\n",
    "            if not r.ok:\n",
    "                info_entry = None\n",
    "            else:\n",
    "                info_entry = json.loads(r.text)\n",
    "        except:\n",
    "            info_entry = None\n",
    "\n",
    "        if not info_entry is None:\n",
    "            clean_html_tags = re.compile('<.*?>')\n",
    "            pharma_name = re.sub(clean_html_tags, '', info_entry[\"name\"])\n",
    "        else:\n",
    "            pharma_name = abbreviated_name\n",
    "        \n",
    "        #Gprot and Barr coupling data from GPCRdb\n",
    "        try:\n",
    "            Gprot_coupling_data,Barr_coupling_data=coupling_Gprot_Barr(uniprotID)\n",
    "        except:\n",
    "            Gprot_coupling_data,Barr_coupling_data=[],[]\n",
    "        \n",
    "        #Variations Uniprot\n",
    "        requestURL = f\"https://www.ebi.ac.uk/proteins/api/variation/{uniprotID}\"\n",
    "\n",
    "        r = requests.get(requestURL, headers={ \"Accept\" : \"application/json\"})\n",
    "        if not r.ok:\n",
    "            variations_uniprot_json = None\n",
    "        else:\n",
    "            variations_uniprot_json = json.loads(r.text)\n",
    "\n",
    "        #Pharmacological info from GPCRdb mutants (API service)\n",
    "        requestURL = f\"https://gpcrdb.org/services/mutants/{abbreviated_name.lower()}\"\n",
    "        try:\n",
    "            r = requests.get(requestURL, headers={ \"Accept\" : \"application/json\"})\n",
    "            pharmaco_data_json = json.loads(r.text)\n",
    "        except:\n",
    "            pharmaco_data_json = None\n",
    "\n",
    "        #name endogenous ligands\n",
    "        ligands = []\n",
    "        for value in endogenous_ligands:\n",
    "            if value[\"receptor\"] == abbreviated_name.lower():\n",
    "                # ligands.append({\"value\":  value[\"ligand_name\"].replace(\"&\", \"\").replace(\";\", \"\").replace(\"<sub>\",\"_\"), \"reference\": \"GPCRdb\"})\n",
    "                ligands.append({\"value\":  value[\"ligand_name\"].replace(\"&\", \"\").replace(\";\", \"\").replace(\"<sup>\",\"\").replace(\"</sup>\",\"\").replace(\"<sub>\",\"\").replace(\"</sub>\",\"\"), \"reference\": \"GPCRdb\"})\n",
    "\n",
    "        seen_values = set()\n",
    "        ligands = [item for item in ligands if item[\"value\"] not in seen_values and not seen_values.add(item[\"value\"])]\n",
    "    \n",
    "        MSA = \"../data/MSA_all_mammalian.fasta\"\n",
    "        microswitches_literature = motifs_microswitches_literature(MSA,uniprotID)\n",
    "\n",
    "        #Retrieve the mutagenesis information from Uniprot related to the entry\n",
    "        if variations_uniprot_json:\n",
    "            mutations_Uniprot = retrieve_mutagenesis_info_Uniprot(variations_uniprot_json,uniprotID)\n",
    "        else:\n",
    "            mutations_Uniprot = []\n",
    "\n",
    "        #Pharmacological data - mutants info from GPCRdb API\n",
    "        if pharmaco_data_json:\n",
    "            pharmaco = retrieve_pharmaco_info_GPCRdb(pharmaco_data_json)\n",
    "        else:\n",
    "            pharmaco = []\n",
    "\n",
    "        #If no conformational biosensor just put confo_biosensor = []\n",
    "        confo_biosensor = []\n",
    "\n",
    "        #find chimeras this parent is involved in\n",
    "        #find the cutting points used for this parent in all designs that we know so we can learn form that\n",
    "        chimeras,all_regions_cutting_pts,related_chimeras=retrieve_involvement_natural_chimeric_design(uniprotID,abbreviated_name,prot_seq,chimeric_design_df)\n",
    "        \n",
    "        cutting_point_values = [] #leave empty if it's a natural\n",
    "\n",
    "        #get structures: pdb id, chain, state and offset\n",
    "        #From the PDB (exp structures)\n",
    "        structures, ds_bonds,interacting_residues, binders, PDBs = retrieve_pdb_dsbonds_interactions(uniprot_json,uniprotID,prot_seq)\n",
    "        #AlphaFold2\n",
    "        structures.append({\"value\":f\"AlphaFold2\",\"chain\": \"A\",  \"state\":\"Undetermined\", \"offset\":  0, \"gaps\": [], \"url\": f\"https://alphafold.ebi.ac.uk/files/AF-{uniprotID}-F1-model_v4.pdb\", \"reference\":\"AFDB\"})\n",
    "        \n",
    "        #AlphaFold multistate. Don't have a AF ms for every GPCR (only humans). Need to check if file exist:\n",
    "        af_ms_active = f\"../examples/3Dstructures/AFms/Active/{uniprotID}.pdb\"\n",
    "        af_ms_inactive = f\"../examples/3Dstructures/AFms/Inactive/{uniprotID}.pdb\"\n",
    "        if os.path.exists(af_ms_active):\n",
    "            af_ms_active = f\"file:///examples/3Dstructures/AFms/Active/{uniprotID}.pdb\"\n",
    "            structures.append({\"value\":f\"AlphaFold2 Active\",\"chain\": \"A\",  \"state\":\"Active\", \"offset\":  0, \"gaps\": [], \"url\": af_ms_active, \"reference\":\"AlphaFold multistate\"})\n",
    "        if os.path.exists(af_ms_inactive):\n",
    "            af_ms_inactive = f\"file:///examples/3Dstructures/AFms/Inactive/{uniprotID}.pdb\"\n",
    "            structures.append({\"value\":f\"AlphaFold2 Inactive\",\"chain\": \"A\",  \"state\":\"Inactive\", \"offset\":  0, \"gaps\": [], \"url\": af_ms_inactive, \"reference\":\"AlphaFold multistate\"})\n",
    "\n",
    "        #retrieve all limits of the secondary structure elements\n",
    "        TM_regions_file = \"../data/Uniprot_TMbed.txt\"\n",
    "        allregions = TMbed_retrieve_allregions(TM_regions_file,uniprotID)\n",
    "\n",
    "        #retrieve the residues interacting with ligand/Gprot/Nb/Ab in PDB and link it to region\n",
    "        #Add manually extra IC and EC contacts\n",
    "        #should follow the following structure: list regrouping all dictionaries with 1 dict per contact\n",
    "        #in dictionary: {\"start\":,\"end\",\"type\",\"description\",\"reference\"}\n",
    "        #for the EC contacts the types can be \"orthosteric\",\"allosteric\",\"VHH EC\"\n",
    "        #for the IC contacts the types can be \"G-protein\",\"VHH IC\"\n",
    "        manual_ICs,manual_ECs = translate_interacting_residues_IC_EC(interacting_residues,binders,PDBs,allregions)\n",
    "\n",
    "        #Features uniprot PTM, binding site uniptor\n",
    "        features = {}\n",
    "        ligand_BS_uniprot,PTMs_uniprot = features_uniprot(uniprot_json,uniprotID)\n",
    "\n",
    "        #Scop3P phosphorylations\n",
    "        PTMs_scop3P=retrieve_PTM_Scop3P(uniprotID)\n",
    "\n",
    "        info = {}\n",
    "\n",
    "        #Abbreviated name\n",
    "        info[\"Abbreviated name\"] = [{\"value\": abbreviated_name.upper(), \"reference\": \"UniProt\"}]\n",
    "\n",
    "        #pharma name\n",
    "        info[\"Pharmacological name\"] = [{\"value\": pharma_name, \"reference\": \"GPCRdb\"}]\n",
    "\n",
    "        #Name\n",
    "        info[\"Name(s)\"] = names\n",
    "\n",
    "        #Uniprot ID\n",
    "        info[\"Uniprot ID\"] = [{\"value\": uniprotID, \"reference\": \"UniProt\"}]\n",
    "\n",
    "        #Species\n",
    "        info[\"Organism\"] =  [{\"value\":species, \"reference\": \"UniProt\"}]\n",
    "\n",
    "        #Class\n",
    "        info[\"Class\"] = [{\"value\":class_, \"reference\": \"GPCRdb\"}]\n",
    "\n",
    "        #Family\n",
    "        info[\"Family\"] = [{\"value\": family, \"reference\": \"GPCRdb\"}]\n",
    "\n",
    "        #Subclass\n",
    "        #Phylogenetically based & Ligand based\n",
    "        info[\"Subclass\"] = {\"Phylogenetically based\": [{\"value\": subclass_phylo, \"reference\": \"10.1124/mol.63.6.1256\"}],\n",
    "                            \"Ligand based\": [{\"value\":subclass_ligand, \"reference\": \"GPCRdb\"}]}\n",
    "\n",
    "        #Endogenous ligand \n",
    "        info[\"Endogenous ligand\"]=ligands\n",
    "\n",
    "        #Gport and Barr coupling data\n",
    "        info[\"G-protein coupling\"]=Gprot_coupling_data\n",
    "        info[\"Beta-arrestin coupling\"] = Barr_coupling_data\n",
    "\n",
    "        #Structures\n",
    "        info[\"Structures\"] = structures\n",
    "\n",
    "        #Info related to chimeric design\n",
    "        info[\"Conformational biosensor\"] = confo_biosensor\n",
    "        info[\"Involvement in chimeric design\"] = chimeras\n",
    "        info[\"Cutting point values\"] = cutting_point_values\n",
    "        info[\"Known cutting points and designs\"] = {\"Known cutting points\":all_regions_cutting_pts,\"Known designs\":related_chimeras}\n",
    "\n",
    "        features['Microswitches'] = microswitches_literature\n",
    "\n",
    "        #remove duplicate PTMs, keep Scop3P ref\n",
    "        PTMs_tot = PTMs_uniprot + PTMs_scop3P\n",
    "        PTMs_tot_unique = remove_duplicates(PTMs_tot,uniprotID,\"start\",descriminator2=\"description\")\n",
    "\n",
    "        features['PTMs'] = PTMs_tot_unique\n",
    "        features['Disulfide bonds'] = ds_bonds\n",
    "        features['Mutagenesis'] = mutations_Uniprot\n",
    "        features['Pharmacological mutagenesis'] = pharmaco\n",
    "\n",
    "        #remove duplicate contacts, keep uniprot\n",
    "        Contacts_EC_tot = manual_ECs + ligand_BS_uniprot\n",
    "        Contacts_EC_unique = remove_duplicates(Contacts_EC_tot,uniprotID,\"start\")\n",
    "\n",
    "        features[\"Contacts\"] = Contacts_EC_unique +  manual_ICs\n",
    "\n",
    "        info[\"Features\"] = features\n",
    "\n",
    "        #Sequence\n",
    "        info[\"Sequence\"] = [{\"value\":prot_seq, \"reference\": \"UniProt\"}]\n",
    "\n",
    "        #Secondary structure info\n",
    "        info[\"Limits regions\"] = allregions\n",
    "\n",
    "        #Gather info that could be useful for chimeric design\n",
    "        known_info = []\n",
    "        if len(confo_biosensor) > 0:\n",
    "            known_info.append({\"value\": \"Confo biosensor\"})\n",
    "        if len(chimeras) > 0:\n",
    "            known_info.append({\"value\": \"Parent chimera\"})\n",
    "        info[\"Known info chimeric design\"] = known_info\n",
    "\n",
    "        json.dump(info, open(f'../examples/json_entries/updated_naturals/{abbreviated_name.upper()}.json', 'w'), indent=2)\n",
    "\n",
    "        # heatmap_cutting_pts(pharma_name,abbreviated_name,chimeras,allregions)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "private_outputs": true,
   "provenance": []
  },
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}