make_alphafold_jobs_af23.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Prepare Sequence Data"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pygtop\n",
    "import requests\n",
    "from Bio import SeqIO\n",
    "from Bio.Seq import Seq\n",
    "from Bio.SeqRecord import SeqRecord\n",
    "import os\n",
    "import json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Receptor Name Receptor Accession    Ligand Name  \\\n",
      "0      5ht1b_human             P28222  5-HT-moduline   \n",
      "1      5ht1b_human             P28222  5-HT-moduline   \n",
      "2      5ht1b_human             P28222  5-HT-moduline   \n",
      "3      5ht1d_human             P28221  5-HT-moduline   \n",
      "4      5ht1d_human             P28221  5-HT-moduline   \n",
      "...            ...                ...            ...   \n",
      "2877      oxyr_rat             P70536       oxytocin   \n",
      "2878      oxyr_rat             P70536       oxytocin   \n",
      "2879      oxyr_rat             P70536       oxytocin   \n",
      "2880      oxyr_rat             P70536       oxytocin   \n",
      "2881      oxyr_rat             P70536       oxytocin   \n",
      "\n",
      "                         InChIKey  Smiles   Sequence  Alternate ID  \\\n",
      "0     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN       LSAL       7048565   \n",
      "1     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN       LSAL           137   \n",
      "2     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN       LSAL       7048566   \n",
      "3     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN       LSAL       7048565   \n",
      "4     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN       LSAL           137   \n",
      "...                           ...     ...        ...           ...   \n",
      "2877  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN  CYIQNCPLG          2042   \n",
      "2878  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN  CYIQNCPLG       DB00107   \n",
      "2879  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN  CYIQNCPLG  CHEMBL395429   \n",
      "2880  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN  CYIQNCPLG          2174   \n",
      "2881  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN  CYIQNCPLG        439302   \n",
      "\n",
      "                     Source  \n",
      "0                   PubChem  \n",
      "1     Guide To Pharmacology  \n",
      "2                   PubChem  \n",
      "3                   PubChem  \n",
      "4     Guide To Pharmacology  \n",
      "...                     ...  \n",
      "2877           Drug Central  \n",
      "2878               DrugBank  \n",
      "2879    ChEMBL_compound_ids  \n",
      "2880  Guide To Pharmacology  \n",
      "2881                PubChem  \n",
      "\n",
      "[2882 rows x 8 columns]\n"
     ]
    }
   ],
   "source": [
    "endo_data = pd.read_csv(\"input_files/peptide_ligands_Arman.csv\")\n",
    "print(endo_data)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             gtp_id                        chainA_name  \\\n",
      "ligand_name                                              \n",
      "INSL3          1995              Insuline-like A chain   \n",
      "LH             1159              Lutropin subunit beta   \n",
      "relaxin-3      1990                            A chain   \n",
      "relaxin-1      1988                            A chain   \n",
      "TSH            3920           Thyrotropin subunit beta   \n",
      "INSL5          2000       Insulin-like peptide B chain   \n",
      "hCG            1160  Choriogonadotropin subunit beta 3   \n",
      "thrombin       4453                        light chain   \n",
      "relaxin        1989                            A chain   \n",
      "\n",
      "                                                    chainA_seq chainA_id  \\\n",
      "ligand_name                                                                \n",
      "INSL3                               AAATNPARYCCLSGCTQQDLLTLCPY    P51460   \n",
      "LH           SREPLRPWCHPINAILAVEKEGCPVCITVNTTICAGYCPTMMRVLQ...    P01229   \n",
      "relaxin-3                             DVLAGLSSSCCKWGCSKSEISSLC    Q8WXF3   \n",
      "relaxin-1                              PYVALFEKCCLIGCTKRSLAKYC    P04808   \n",
      "TSH          FCIPTEYTMHIERRECAYCLTINTTICAGYCMTRDINGKLFLPKYA...    P01222   \n",
      "INSL5                                 KESVRLCGLEYIRTVIYICASSRW    Q9Y5Q6   \n",
      "hCG          SKEPLRPRCRPINATLAVEKEGCPVCITVNTTICAGYCPTMTRVLQ...    P0DN86   \n",
      "thrombin                  TFGSGEADCGLRPLFEKKSLEDKTERELLESYIDGR    P00734   \n",
      "relaxin                               QLYSALANKCCHVGCTKRSLARFC    P04090   \n",
      "\n",
      "                                   chainB_name  \\\n",
      "ligand_name                                      \n",
      "INSL3                    Insuline-like B chain   \n",
      "LH           Glycoprotein hormones alpha chain   \n",
      "relaxin-3                              B chain   \n",
      "relaxin-1                              B chain   \n",
      "TSH          Glycoprotein hormones alpha chain   \n",
      "INSL5                    Insuline-like B chain   \n",
      "hCG          Glycoprotein hormones alpha chain   \n",
      "thrombin                           heavy chain   \n",
      "relaxin                                B chain   \n",
      "\n",
      "                                                    chainB_seq chainB_id  \n",
      "ligand_name                                                               \n",
      "INSL3                      LGPAPTPEMREKLCGHHFVRALVRVCGGPRWSTEA    P51460  \n",
      "LH           XAPDVQDCPECTLQENPFFSQPGAPILQCMGCCFSRAYPTPLRSKK...    P01215  \n",
      "relaxin-3                          RAAPYGVRLCGREFIRAVIFTCGGSRW    Q8WXF3  \n",
      "relaxin-1                      VAAKWKDDVIKLCGRELVRAQIAICGMSTWS    P04808  \n",
      "TSH          XAPDVQDCPECTLQENPFFSQPGAPILQCMGCCFSRAYPTPLRSKK...    P01215  \n",
      "INSL5                                    QDLQTLCCTDGCSMTDLSALC    Q9Y5Q6  \n",
      "hCG          XAPDVQDCPECTLQENPFFSQPGAPILQCMGCCFSRAYPTPLRSKK...    P01215  \n",
      "thrombin     IVEGSDAEIGMSPWQVMLFRKSPQELLCGASLISDRWVLTAAHCLL...    P00734  \n",
      "relaxin                          DSWMEEVIKLCGRELVRAQIAICGMSTWS    P04090  \n"
     ]
    }
   ],
   "source": [
    "spec_lig_data = pd.read_csv(\"input_files/peptide_specligands_Arman.csv\", index_col=0)\n",
    "print(spec_lig_data)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Receptor Name                         gpr39_human\n",
      "Receptor Accession                         O43194\n",
      "Ligand Name                           neurotensin\n",
      "InChIKey              PCJGZPGTCUMMOT-ISULXFBGSA-N\n",
      "Smiles                                        NaN\n",
      "Sequence                            XLYENKPRRPYIL\n",
      "Alternate ID                                 1579\n",
      "Source                      Guide To Pharmacology\n",
      "Name: 644, dtype: object\n",
      "\n",
      "Receptor Name                   gpr15_human\n",
      "Receptor Accession                   P49685\n",
      "Ligand Name                          GPR15L\n",
      "InChIKey                                NaN\n",
      "Smiles                                  NaN\n",
      "Sequence                                NaN\n",
      "Alternate ID                          10567\n",
      "Source                Guide To Pharmacology\n",
      "dtype: object\n",
      "\n",
      "Receptor Name                   gpr15_human\n",
      "Receptor Accession                   P49685\n",
      "Ligand Name                    GPR15L_71-81\n",
      "InChIKey                                NaN\n",
      "Smiles                                  NaN\n",
      "Sequence                                NaN\n",
      "Alternate ID                          10568\n",
      "Source                Guide To Pharmacology\n",
      "Name: 2883, dtype: object\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/ls/ljr1g7ps2xd3yb1pnstmg4440000gn/T/ipykernel_13447/1094284396.py:15: FutureWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "  rec_row = pd.Series(index=endo_data.columns)\n",
      "/var/folders/ls/ljr1g7ps2xd3yb1pnstmg4440000gn/T/ipykernel_13447/1094284396.py:21: FutureWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "  lig_row = pd.Series(index=endo_data.columns)\n",
      "/var/folders/ls/ljr1g7ps2xd3yb1pnstmg4440000gn/T/ipykernel_13447/1094284396.py:21: FutureWarning: The default dtype for empty Series will be 'object' instead of 'float64' in a future version. Specify a dtype explicitly to silence this warning.\n",
      "  lig_row = pd.Series(index=endo_data.columns)\n"
     ]
    }
   ],
   "source": [
    "# manually adding GPR39_human/neurotensin complex and GRP15/GRP15L\n",
    "\n",
    "fields = {'Receptor Name': \"rec\", 'Receptor Accession': \"rec\", 'Ligand Name': \"lig\",\n",
    "          'InChIKey': \"lig\", 'Smiles': \"lig\", 'Sequence': \"lig\", 'Alternate ID': \"lig\",\n",
    "          'Source': \"lig\", #'Receptor Sequence': \"rec\"\n",
    "          }\n",
    "\n",
    "for rec_name, rec_acc, lig_name, lig_gtp in [(\"gpr39_human\", \"O43194\", \"neurotensin\", \"1579\"), # neurotensin\n",
    "                          (\"gpr15_human\", \"P49685\", \"GPR15L\", \"10567\"), # GRP15L full length\n",
    "                          (\"gpr15_human\", \"P49685\", \"GPR15L_71-81\", \"10568\") # GRP15L shorter version\n",
    "                          ]:\n",
    "    try:\n",
    "        rec_row = endo_data[endo_data[\"Receptor Name\"] == rec_name].iloc[0, :]\n",
    "    except IndexError:\n",
    "        rec_row = pd.Series(index=endo_data.columns)\n",
    "        rec_row[\"Receptor Name\"] = rec_name\n",
    "        rec_row[\"Receptor Accession\"] = rec_acc\n",
    "    try:\n",
    "        lig_row = endo_data[endo_data[\"Alternate ID\"] == lig_gtp].iloc[0, :]\n",
    "    except IndexError:\n",
    "        lig_row = pd.Series(index=endo_data.columns)\n",
    "        lig_row[\"Ligand Name\"] = lig_name\n",
    "        lig_row[\"Alternate ID\"] = lig_gtp\n",
    "        lig_row[\"Source\"] = \"Guide To Pharmacology\"\n",
    "    new_row = rec_row.copy()\n",
    "    for key, from_row in fields.items():\n",
    "        if from_row == \"rec\":\n",
    "            new_row[key] = rec_row[key]\n",
    "        else:\n",
    "            new_row[key] = lig_row[key]\n",
    "    print()\n",
    "    print(new_row)\n",
    "    endo_data.loc[len(endo_data.index)] = new_row"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">  Sequences for ligand 'GPR15L_71-81' missing \n",
      "   > Guide to Pharmacology ID is 10568\n",
      "> Ligand 'relaxin' is in special set (has multiple chains). Getting sequences from there\n",
      "> Added partially missing sequences to ligand 'FSH' entries \n",
      "> Ligand 'INSL3' is in special set (has multiple chains). Getting sequences from there\n",
      "> Ligand 'thrombin' is in special set (has multiple chains). Getting sequences from there\n",
      "> Ligand 'relaxin-1' is in special set (has multiple chains). Getting sequences from there\n",
      "> Ligand 'LH' is in special set (has multiple chains). Getting sequences from there\n",
      "> Ligand 'INSL5' is in special set (has multiple chains). Getting sequences from there\n",
      "> Ligand 'TSH' is in special set (has multiple chains). Getting sequences from there\n",
      "> Ligand 'hCG' is in special set (has multiple chains). Getting sequences from there\n",
      ">  Sequences for ligand 'GPR15L' missing \n",
      "   > Guide to Pharmacology ID is 10567\n",
      "   > UniProt accession is Q6UWK7\n",
      "> Ligand 'relaxin-3' is in special set (has multiple chains). Getting sequences from there\n",
      "\n",
      "> 2 ligands have missing sequences\n"
     ]
    }
   ],
   "source": [
    "column_name = \"Sequence\"\n",
    "counter = 0\n",
    "ligand_list = list(set(endo_data[\"Ligand Name\"]))\n",
    "\n",
    "# loop over all ligands to check for sequence\n",
    "for lig in ligand_list:\n",
    "    lig_data = endo_data[endo_data[\"Ligand Name\"] == lig]\n",
    "\n",
    "    # ligand has sequence but not for all entries\n",
    "    if pd.isna(lig_data[\"Sequence\"]).any() and not pd.isna(lig_data[\"Sequence\"]).all():\n",
    "        print(\"> Added partially missing sequences to ligand '%s' entries \" % lig)\n",
    "        # get sequence from first ligand entry with a sequence\n",
    "        sequence = lig_data[~pd.isna(lig_data[\"Sequence\"])][\"Sequence\"].iloc[0]\n",
    "        # set sequence for all ligand entries without sequence\n",
    "        endo_data.loc[endo_data[\"Ligand Name\"] == lig, \"Sequence\"] = sequence\n",
    "    if pd.isna(lig_data[\"Sequence\"]).all():\n",
    "        if lig in spec_lig_data.index:\n",
    "            print(\"> Ligand '%s' is in special set (has multiple chains). Getting sequences from there\" % lig)\n",
    "            lig_entry = str(spec_lig_data.loc[lig, :].to_dict())\n",
    "            endo_data.loc[endo_data[\"Ligand Name\"] == lig, \"Sequence\"] = lig_entry\n",
    "            continue\n",
    "        print(\">  Sequences for ligand '%s' missing \" % lig)\n",
    "        counter += 1\n",
    "        gtp_id = int(lig_data[lig_data[\"Source\"] == \"Guide To Pharmacology\"][\"Alternate ID\"].values[0])\n",
    "\n",
    "        print(\"   > Guide to Pharmacology ID is %s\" % gtp_id)\n",
    "        my_drug = pygtop.get_ligand_by_id(gtp_id)\n",
    "\n",
    "        sequence = my_drug.one_letter_sequence()\n",
    "        endo_data.loc[endo_data[\"Ligand Name\"] == lig, \"Sequence\"] = sequence\n",
    "\n",
    "        db_links = my_drug.database_links()\n",
    "        for db_link in db_links:\n",
    "            db_link = db_link.json_data\n",
    "            if db_link[\"database\"] == \"UniProtKB\":\n",
    "                uniprot_id = db_link[\"accession\"]\n",
    "                print(\"   > UniProt accession is %s\" % uniprot_id)\n",
    "            elif db_link[\"database\"] == \"PubChem\":\n",
    "                pubchem_id = db_link[\"accession\"]\n",
    "                print(\"   > PubChem accession is %s\" % pubchem_id)\n",
    "            else:\n",
    "                continue\n",
    "\n",
    "print(\"\\n> %i ligands have missing sequences\" % counter)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> Sequence of 'ednrb_rat' could not be retrieved.\n"
     ]
    }
   ],
   "source": [
    "receptor_list = list(set(endo_data[\"Receptor Name\"].values))\n",
    "\n",
    "def getSequence(protein_id):\n",
    "    #print(\" > Querying GPCRdb for protein '%s' \" % protein_id)\n",
    "    url = \"https://gpcrdb.org/services/protein/\"+protein_id\n",
    "    #print(\"   \"+url)\n",
    "    response = requests.get(url)\n",
    "    protein_data = response.json()\n",
    "    return protein_data['sequence']\n",
    "\n",
    "# Getting receptor sequences from GPCRdb\n",
    "for rec in receptor_list:\n",
    "    try:\n",
    "        sequence = getSequence(rec)\n",
    "    except KeyError:\n",
    "        print(\"> Sequence of '%s' could not be retrieved.\" % rec)\n",
    "        pass\n",
    "    endo_data.loc[endo_data[\"Receptor Name\"] == rec, \"Receptor Sequence\"] = sequence"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "# manually adding receptor sequence for one missing entry on GPCRdb\n",
    "for entry, fastafile in [(\"ntr2_mouse\", \"input_files/P70310_NTR2_mouse.fasta\")]:\n",
    "    record = SeqIO.parse(open(fastafile, \"r\"), \"fasta\")\n",
    "    sequence = str(list(record)[0].seq)\n",
    "    endo_data.loc[endo_data[\"Receptor Name\"] == entry, \"Receptor Sequence\"] = sequence"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> 0 receptors have missing sequences\n"
     ]
    }
   ],
   "source": [
    "print(\"> %i receptors have missing sequences\" % list(pd.isna(endo_data[\"Receptor Sequence\"])).count(True))"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " > 736 unique complex pairs to model\n"
     ]
    }
   ],
   "source": [
    "endo_data_unique = endo_data[endo_data[\"Source\"] == \"Guide To Pharmacology\"]\n",
    "endo_data_unique = endo_data_unique.drop_duplicates(subset=[\"Receptor Name\", \"Sequence\"])\n",
    "print(\" > %i unique complex pairs to model\" % len(endo_data_unique))\n",
    "endo_data_unique.to_csv(\"input_files/peptide_ligands_Arman_completed.csv\")"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Prepare AlphaFold Jobs"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [],
   "source": [
    "# helper function to (quick and dirty) generate AF scripts for ILF cluster\n",
    "def af_ilf_jobscript(jobname,\n",
    "                     fastafile,\n",
    "                     jobscript_file,\n",
    "                     only_compute_MSAs,\n",
    "                     use_precompute_MSAs,\n",
    "                     working_dir=\"$PWD\",\n",
    "                     singleton_dependency=False,\n",
    "                     do_amber=True,\n",
    "                     ncores=6,\n",
    "                     mem=20):\n",
    "    af_dir = \"/projects/ilfgrid/apps/alphafold-2.3.1\"\n",
    "    af_db_dir = \"/projects/ilfgrid/data/alphafold-genetic-databases/\"\n",
    "    if only_compute_MSAs:\n",
    "        cp_cmd = \"\"\n",
    "        gswitch = \"#SBATCH --exclude=ilfgridgpun02fl\"\n",
    "        af_cmd = f\"cd $AF_DIR && /projects/ilfgrid/apps/alphafold-2.2.0/prep_MSAs.sh -d $AFDB_DIR -n {ncores} -m {mem} -o $WDIR/ -f $FASTAFILE\"\n",
    "\n",
    "    else:\n",
    "        cp_cmd = \"cp $WDIR/$BASENAME/ranked_0.pdb $WDIR/${BASENAME}.pdb\"\n",
    "        gswitch = \"#SBATCH --gres=gpu:1                                                         # no. of gpus (can only use one)\"\n",
    "        gpu_check = \"\"\n",
    "        af_cmd = f\"cd $AF_DIR && $AF_DIR/run_alphafold.sh -d $AFDB_DIR -n {ncores} -o $WDIR/ -f $FASTAFILE -t 2023-01-01 -m multimer -r {str(do_amber).lower()} -l 1 -e false -p {str(use_precompute_MSAs).lower()} -a $SLURM_JOB_GPUS\"\n",
    "\n",
    "\n",
    "    if singleton_dependency:\n",
    "        dependency = \"#SBATCH --dependency=singleton\"\n",
    "    else:\n",
    "        dependency = \"\"\n",
    "\n",
    "    command = f'''#!/bin/bash -l\n",
    "#SBATCH --job-name={jobname}\n",
    "###SBATCH --output={working_dir}/slurm-%A.out\n",
    "###SBATCH --chdir={af_dir}/                    # change location to the install directory\n",
    "#SBATCH --partition=standard                                                 # the queue you submit to\n",
    "#SBATCH --mem={mem}G                                                            # the amount of memory to reserve\n",
    "#SBATCH --ntasks={ncores}                                                          # combined with --cpus-per-task it determines\n",
    "#SBATCH --cpus-per-task=1                                                    # the number of cpus to use (here: 12)\n",
    "#SBATCH --nodes=1                                                            # no. of nodes to use\n",
    "{gswitch}\n",
    "{dependency}\n",
    "echo \">> Job started at: $(date)\"\n",
    "echo \"   On Machine:     $(hostname)\"\n",
    "echo -e \"   In directory:   $(pwd)\\\\n\"\n",
    "echo -e \">> GPU is \\\"${{SLURM_JOB_GPUS}}\\\"\\\\n\"\n",
    "start=`date +%s`\n",
    "\n",
    "WDIR=\"{working_dir}\"\n",
    "FASTAFILE=\"$WDIR/{fastafile}\"\n",
    "\n",
    "AF_DIR=\"{af_dir}\"\n",
    "AFDB_DIR=\"{af_db_dir}/\"\n",
    "module load miniconda\n",
    "conda activate $AF_DIR/AF_2.3.1_conda_env\n",
    "\n",
    "'''\n",
    "\n",
    "    command = command+af_cmd\n",
    "\n",
    "    command = command+f'''\n",
    "\n",
    "BASENAME=`basename $FASTAFILE .fasta`\n",
    "{cp_cmd}\n",
    "\n",
    "end=`date +%s`\n",
    "runtime=$((end-start))\n",
    "echo -e \\\"\\\\n>> Job finished at: $(date)\\\"\n",
    "echo -e \\\"\\\\n>> Runtime: $runtime s\\\"'''\n",
    "\n",
    "    with open(os.path.join(jobscript_file), \"w+\") as f:\n",
    "        f.write(command)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Receptor Name Receptor Accession    Ligand Name  \\\n",
      "1      5ht1b_human             P28222  5-HT-moduline   \n",
      "4      5ht1d_human             P28221  5-HT-moduline   \n",
      "6       c3ar_human             Q16581            C3a   \n",
      "7       c3ar_human             Q16581            C3a   \n",
      "8       c3ar_human             Q16581            C5a   \n",
      "...            ...                ...            ...   \n",
      "2846    oxyr_human             P30559       oxytocin   \n",
      "2864      oxyr_rat             P70536       oxytocin   \n",
      "2882   gpr39_human             O43194    neurotensin   \n",
      "2883   gpr15_human             P49685         GPR15L   \n",
      "2884   gpr15_human             P49685   GPR15L_71-81   \n",
      "\n",
      "                         InChIKey  Smiles  \\\n",
      "1     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN   \n",
      "4     IMIVWAUMTAIVPJ-XUXIUFHCSA-N     NaN   \n",
      "6                             NaN     NaN   \n",
      "7                             NaN     NaN   \n",
      "8                             NaN     NaN   \n",
      "...                           ...     ...   \n",
      "2846  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN   \n",
      "2864  XNOPRXBHLZRZKH-DSZYJQQASA-N     NaN   \n",
      "2882  PCJGZPGTCUMMOT-ISULXFBGSA-N     NaN   \n",
      "2883                          NaN     NaN   \n",
      "2884                          NaN     NaN   \n",
      "\n",
      "                                               Sequence  Alternate ID  \\\n",
      "1                                                  LSAL           137   \n",
      "4                                                  LSAL           137   \n",
      "6     SVQLTEKRMDKVGKYPKELRKCCEDGMRENPMRFSCQRRTRFISLG...          3640   \n",
      "7     SVQLMERRMDKAGQYTDKGLRKCCEDGMRDIPMPYSCQRRARLITQ...          3642   \n",
      "8     TLQKKIEEIAAKYKHSVVKKCCYDGACVNNDETCEQRAARISLGPR...           573   \n",
      "...                                                 ...           ...   \n",
      "2846                                          CYIQNCPLG          2176   \n",
      "2864                                          CYIQNCPLG          2176   \n",
      "2882                                      XLYENKPRRPYIL          1579   \n",
      "2883  MRLLVLSSLLCILLLCFSIFSTEGKRRPAKAWSGRRTRLCCHRVPS...         10567   \n",
      "2884                                        LWVVPGALPQV         10568   \n",
      "\n",
      "                     Source                                  Receptor Sequence  \n",
      "1     Guide To Pharmacology  MEEPGAQCAPPPPAGSETWVPQANLSSAPSQNCSAKDYIYQDSISL...  \n",
      "4     Guide To Pharmacology  MSPLNQSAEGLPQEASNRSLNATETSEAWDPRTLQALKISLAVVLS...  \n",
      "6     Guide To Pharmacology  MASFSAETNSTDLLSQPWNEPPVILSMVILSLTFLLGLPGNGLVLW...  \n",
      "7     Guide To Pharmacology  MASFSAETNSTDLLSQPWNEPPVILSMVILSLTFLLGLPGNGLVLW...  \n",
      "8     Guide To Pharmacology  MASFSAETNSTDLLSQPWNEPPVILSMVILSLTFLLGLPGNGLVLW...  \n",
      "...                     ...                                                ...  \n",
      "2846  Guide To Pharmacology  MEGALAANWSAEAANASAAPPGAEGNRTAGPPRRNEALARVEVAVL...  \n",
      "2864  Guide To Pharmacology  MEGTPAANWSVELDLGSGVPPGEEGNRTAGPPQRNEALARVEVAVL...  \n",
      "2882  Guide To Pharmacology  MASPSLPGSDCSQIIDHSHVPEFEVATWIKITLILVYLIIFVMGLL...  \n",
      "2883  Guide To Pharmacology  MDPEETSVYLDYYYATSPNSDIRETHSHVPYTSVFLPVFYTAVFLT...  \n",
      "2884  Guide To Pharmacology  MDPEETSVYLDYYYATSPNSDIRETHSHVPYTSVFLPVFYTAVFLT...  \n",
      "\n",
      "[736 rows x 9 columns]\n"
     ]
    }
   ],
   "source": [
    "endo_data_unique = pd.read_csv(\"input_files/peptide_ligands_Arman_completed.csv\", index_col=0)\n",
    "print(endo_data_unique)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entries with X in sequences: 42\n"
     ]
    }
   ],
   "source": [
    "endo_data_unique[\"has_X\"] = False\n",
    "for index, row in endo_data_unique.iterrows():\n",
    "    if \"X\" in row[\"Receptor Sequence\"] or \"X\" in row[\"Sequence\"]:\n",
    "        endo_data_unique.loc[index, \"has_X\"] = True\n",
    "\n",
    "print(\"Entries with X in sequences:\", list(endo_data_unique[\"has_X\"].values).count(True))"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " >>> Complex rxfp1_human-1988 has a multi-chain ligand\n",
      "{'gtp_id': 1988, 'chainA_name': 'A chain', 'chainA_seq': 'PYVALFEKCCLIGCTKRSLAKYC', 'chainA_id': 'P04808', 'chainB_name': 'B chain', 'chainB_seq': 'VAAKWKDDVIKLCGRELVRAQIAICGMSTWS', 'chainB_id': 'P04808'}\n",
      " > Scripts for complex rxfp1_human-1988 generated \n",
      " >>> Complex rxfp1_human-1989 has a multi-chain ligand\n",
      "{'gtp_id': 1989, 'chainA_name': 'A chain', 'chainA_seq': 'QLYSALANKCCHVGCTKRSLARFC', 'chainA_id': 'P04090', 'chainB_name': 'B chain', 'chainB_seq': 'DSWMEEVIKLCGRELVRAQIAICGMSTWS', 'chainB_id': 'P04090'}\n",
      " > Scripts for complex rxfp1_human-1989 generated \n",
      " >>> Complex rxfp1_human-1990 has a multi-chain ligand\n",
      "{'gtp_id': 1990, 'chainA_name': 'A chain', 'chainA_seq': 'DVLAGLSSSCCKWGCSKSEISSLC', 'chainA_id': 'Q8WXF3', 'chainB_name': 'B chain', 'chainB_seq': 'RAAPYGVRLCGREFIRAVIFTCGGSRW', 'chainB_id': 'Q8WXF3'}\n",
      " > Scripts for complex rxfp1_human-1990 generated \n",
      " >>> Complex rl3r1_human-2000 has a multi-chain ligand\n",
      "{'gtp_id': 2000, 'chainA_name': 'Insulin-like peptide B chain', 'chainA_seq': 'KESVRLCGLEYIRTVIYICASSRW', 'chainA_id': 'Q9Y5Q6', 'chainB_name': 'Insuline-like B chain', 'chainB_seq': 'QDLQTLCCTDGCSMTDLSALC', 'chainB_id': 'Q9Y5Q6'}\n",
      " > Scripts for complex rl3r1_human-2000 generated \n",
      " >>> Complex rl3r2_human-2000 has a multi-chain ligand\n",
      "{'gtp_id': 2000, 'chainA_name': 'Insulin-like peptide B chain', 'chainA_seq': 'KESVRLCGLEYIRTVIYICASSRW', 'chainA_id': 'Q9Y5Q6', 'chainB_name': 'Insuline-like B chain', 'chainB_seq': 'QDLQTLCCTDGCSMTDLSALC', 'chainB_id': 'Q9Y5Q6'}\n",
      " > Scripts for complex rl3r2_human-2000 generated \n",
      " >>> Complex rl3r2_mouse-6519 has a multi-chain ligand\n",
      "{'gtp_id': 2000, 'chainA_name': 'Insulin-like peptide B chain', 'chainA_seq': 'KESVRLCGLEYIRTVIYICASSRW', 'chainA_id': 'Q9Y5Q6', 'chainB_name': 'Insuline-like B chain', 'chainB_seq': 'QDLQTLCCTDGCSMTDLSALC', 'chainB_id': 'Q9Y5Q6'}\n",
      " > Scripts for complex rl3r2_mouse-6519 generated \n"
     ]
    }
   ],
   "source": [
    "# build directories and make AF jobs\n",
    "af_dir = \"af_jobs_af23\"\n",
    "os.makedirs(af_dir, exist_ok=True)\n",
    "for i, row in endo_data_unique.iterrows():\n",
    "    # Make directories\n",
    "    rec = row[\"Receptor Name\"]\n",
    "    rec_dir = os.path.join(af_dir, rec)\n",
    "    os.makedirs(rec_dir, exist_ok=True)\n",
    "\n",
    "    lig = str(row[\"Alternate ID\"])\n",
    "    complex_dir = os.path.join(rec_dir, lig)\n",
    "    complex = \"%s-%s\" % (rec, lig)\n",
    "    os.makedirs(os.path.join(complex_dir), exist_ok=True)\n",
    "    if os.path.exists(os.path.join(complex_dir, \"%s.pdb\" % complex)):\n",
    "        continue\n",
    "\n",
    "    # save info about complex to json\n",
    "    row.to_json(os.path.join(complex_dir, \"%s.json\" % complex))\n",
    "\n",
    "    # prepare sequence file for AF\n",
    "    rec_record = SeqRecord(Seq(row[\"Receptor Sequence\"]),\n",
    "                           id=row[\"Receptor Name\"],\n",
    "                           description=row[\"Receptor Accession\"])\n",
    "\n",
    "    # check if multi-chain ligand\n",
    "    try:\n",
    "        multi_lig_entry = eval(row[\"Sequence\"])\n",
    "    except NameError:\n",
    "        # not a multi-chain ligand\n",
    "        lig_record = SeqRecord(Seq(row[\"Sequence\"]),\n",
    "                           id=row[\"Ligand Name\"],\n",
    "                           description=\"gtp:%i\" % row[\"Alternate ID\"])\n",
    "        records = [rec_record, lig_record]\n",
    "    else:\n",
    "        # multi-chain ligand\n",
    "        print(\" >>> Complex %s has a multi-chain ligand\" % complex)\n",
    "        print(multi_lig_entry)\n",
    "        json_file = os.path.join(complex_dir, \"multi_chain_ligand.json\")\n",
    "        json.dump(multi_lig_entry, open(json_file, \"w+\"), indent=4)\n",
    "\n",
    "        lig1_record = SeqRecord(Seq(str(multi_lig_entry[\"chainA_seq\"])),\n",
    "                   id=multi_lig_entry[\"chainA_name\"],\n",
    "                   description=\"| uniprot:%s | part of ligand '%s' with gtp:%s\" % (multi_lig_entry[\"chainA_id\"], row[\"Ligand Name\"], multi_lig_entry[\"gtp_id\"]))\n",
    "        lig2_record = SeqRecord(Seq(str(multi_lig_entry[\"chainB_seq\"])),\n",
    "                   id=multi_lig_entry[\"chainB_name\"],\n",
    "                   description=\"| uniprot:%s | part of ligand '%s' with gtp:%s\" % (multi_lig_entry[\"chainB_id\"], row[\"Ligand Name\"], multi_lig_entry[\"gtp_id\"]))\n",
    "        records = [rec_record, lig1_record, lig2_record]\n",
    "\n",
    "    # save sequence file\n",
    "    fastafile = os.path.join(complex_dir, \"%s.fasta\" % complex)\n",
    "    SeqIO.write(records, open(fastafile, \"w+\"), \"fasta\")\n",
    "\n",
    "    # make AF scripts\n",
    "    # first script to blast DBs and make MSAs (CPU only)\n",
    "    af_ilf_jobscript(jobname=complex,\n",
    "                     fastafile=os.path.basename(fastafile),\n",
    "                     jobscript_file=os.path.join(complex_dir, \"0_prepMSAs.run\"),\n",
    "                     use_precompute_MSAs=False, only_compute_MSAs=True,\n",
    "                     ncores=2)\n",
    "\n",
    "    do_amber = False if row[\"has_X\"] else True\n",
    "\n",
    "    # second script to build AF models (needs GPU)\n",
    "    af_ilf_jobscript(jobname=complex,\n",
    "                     fastafile=os.path.basename(fastafile),\n",
    "                     jobscript_file=os.path.join(complex_dir, \"1_buildmodels.run\"),\n",
    "                     use_precompute_MSAs=True, only_compute_MSAs=False,\n",
    "                     singleton_dependency=True,\n",
    "                     do_amber=do_amber,\n",
    "                     ncores=4)\n",
    "\n",
    "    print(\" > Scripts for complex %s generated \" % complex)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# If generated locally, move whole filesystem generated under ./af_jobs to the cluster. Then use go into ./af_jobs on the cluster and submit the jobs with \"WD=$PWD; for JOBDIR in */*; do cd $JOBDIR; cd $WD; sbatch < 0_prepMSAs.run; sbatch < 1_buildmodels.run; done\" or similar"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}