OldBindingDB_Investigation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from rdkit.Chem import AllChem as Chem\n",
    "import pandas as pd\n",
    "from Bio import pairwise2\n",
    "import numpy as np\n",
    "from functools import partial\n",
    "from multiprocessing import Pool, cpu_count\n",
    "import pickle,re\n",
    "# from rdkit.Chem import rdFMCS\n",
    "# from rdkit.Chem import MolToInchiKey\n",
    "from glob import glob\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import permutations,combinations,combinations_with_replacement\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/anm329/anaconda3/envs/RDK/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3071: DtypeWarning: Columns (15,19,20,51,57,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n"
     ]
    }
   ],
   "source": [
    "ligand_info = pd.read_csv('Ligand_info.tsv', sep='\\t')\n",
    "rec_pdbs = pd.read_csv('cri_bindb_old.txt', sep=' ',header=None,usecols=[0,2])\n",
    "rec_pdbs.columns = ['Rec','Seq']\n",
    "rec_pdbs.drop_duplicates(subset='Rec',inplace=True)\n",
    "rec_pdbs.set_index(['Rec'],inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "print(os.path.isfile('simi_mat/matrix.pickle'))\n",
    "with open('simi_mat/matrix.pickle', 'rb') as file:\n",
    "        (distanceMatrix, target_names,ligandsim) = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcClusterGroups(dists, target_names, t):\n",
    "    '''dists is a distance matrix (full) for target_names'''\n",
    "    assigned = set()\n",
    "    groups = []\n",
    "    for i in range(dists.shape[0]):\n",
    "        if i not in assigned:\n",
    "            group = assignGroup(dists, t, set([i]),target_names)\n",
    "            groups.append(group)\n",
    "            assigned.update(group)\n",
    "    return [set(target_names[i] for i in g) for g in groups]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def assignGroup(dists, t, explore, names):\n",
    "    '''group targets that are less than t away from each other and what's in explore'''\n",
    "    group = set(explore)\n",
    "    while explore:\n",
    "        frontier = set()\n",
    "        for i in explore:\n",
    "            for j in range(dists.shape[1]):\n",
    "                if j not in group:\n",
    "                    #add to the group if protein is close by threshold t (these are distances - default 0.5)\n",
    "                    #also add if the ligands are more similar (not distance) than ligandt and \n",
    "                    #the protein is closer than t2 (default 0.8 - meaning more than 20% similar)\n",
    "                    if dists[i][j] < t:\n",
    "                        group.add(j)\n",
    "                        frontier.add(j)                \n",
    "                                        \n",
    "        explore = frontier\n",
    "    return group\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cUTDM2(needed, options):\n",
    "    '''compute distance between target pair'''\n",
    "    mindist = 1.0\n",
    "    best_opt=[]\n",
    "    substring = False\n",
    "    for idx,opt in enumerate(options):\n",
    "        score = pairwise2.align.globalxx(needed, opt, score_only=True)\n",
    "        length = max(len(needed), len(opt))\n",
    "        distance = (length-score)/length\n",
    "        if distance < mindist:\n",
    "            mindist = distance\n",
    "            best_opt = [idx]\n",
    "            substring = needed in opt\n",
    "        elif distance == mindist:\n",
    "            best_opt.append(idx)\n",
    "    return best_opt, mindist, substring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "similarity = 0.9\n",
    "thresh = 1-similarity \n",
    "cluster_groups = calcClusterGroups(distanceMatrix, target_names, thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for group in cluster_groups:\n",
    "final_df = pd.DataFrame(columns=['Receptor','mol','file','BindingDB MonomerID','Inchi Key'])\n",
    "seen_recs = set()\n",
    "for idx,group in enumerate(cluster_groups):\n",
    "    lig_mols = []\n",
    "    for rec in group:\n",
    "        for lig in glob('separated_sets/{}/*.sdf'.format(rec)):\n",
    "            mol = Chem.MolFromMolFile(lig,sanitize=False)\n",
    "            filename= lig.split('/')[-1]\n",
    "            if mol is None:\n",
    "                print(filename)\n",
    "            lig_id = re.findall(r'\\d+',mol.GetProp('_Name'))[0]\n",
    "            lig_mols.append([rec,mol,filename,lig_id,MolToInchiKey(mol)])\n",
    "    group_df = pd.DataFrame(data=lig_mols,columns=['Receptor','mol','file','BindingDB MonomerID','Inchi Key'])\n",
    "    no_dup_group = group_df.drop_duplicates(subset='Inchi Key')\n",
    "    final_df = final_df.append(no_dup_group,ignore_index=True)\n",
    "    \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('{} unique receptors'.format(len(final_df['Receptor'].unique())))\n",
    "experimental_df = final_df.groupby('Receptor').filter(lambda x: len(x) > 1) # make sure receptors have more than 1 ligand\n",
    "print('{} unique receptors with greater than 1 ligand'.format(len(experimental_df['Receptor'].unique())))\n",
    "experimental_df.to_csv('best_experimental_df_475recs.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','Inchi Key'],header=False,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experimental_df = pd.read_csv('best_experimental_df_475recs.txt', sep=' ', names=['Receptor','file','BindingDB MonomerID','Inchi Key'],header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(experimental_df.groupby('Receptor').count()['file'],bins=np.arange(1,max(experimental_df.groupby('Receptor').count()['file']),1),align='mid')\n",
    "plt.axvline(x=np.mean(experimental_df.groupby('Receptor').count()['file']),label='Average',color='k',linestyle='--')\n",
    "plt.xlabel('Ligands in Series')\n",
    "plt.xlim([1.5,max(experimental_df.groupby('Receptor').count()['file'])])\n",
    "plt.ylabel('Frequency')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IC50=dict()\n",
    "bad_rows=[]\n",
    "experimental_df.assign(IC50=np.nan)\n",
    "for idx, row in experimental_df.iterrows():\n",
    "    base_info = ligand_info[ligand_info['BindingDB MonomerID'] == int(row['BindingDB MonomerID'])]\n",
    "    pertinent_info = base_info[base_info['BindingDB Target Chain  Sequence'].notna()]\n",
    "    key='{}:{}'.format(row['Receptor'],row['BindingDB MonomerID'])\n",
    "    if len(pertinent_info) == 1:\n",
    "        IC50[key] = pertinent_info['IC50 (nM)'].item()\n",
    "        experimental_df.at[idx,'IC50'] = pertinent_info['IC50 (nM)'].item()\n",
    "    elif len(pertinent_info) > 1:\n",
    "        sequences = pertinent_info['BindingDB Target Chain  Sequence'].values.tolist()\n",
    "        idx_good, dist, sub = cUTDM2(rec_pdbs.loc[row['Receptor'],'Seq'],sequences)\n",
    "        IC50[key] = pertinent_info.iloc[idx_good,:]['IC50 (nM)']\n",
    "        if len(idx_good) == 1:\n",
    "            experimental_df.at[idx,'IC50'] = pertinent_info.iloc[idx_good]['IC50 (nM)'].item()\n",
    "            experimental_df.at[idx,'distance'] = dist\n",
    "        else:\n",
    "            experimental_df.at[idx,'distance'] = dist\n",
    "            experimental_df.at[idx,'IC50'] = pertinent_info.iloc[idx_good]['IC50 (nM)'].dropna().values.tolist()\n",
    "    else:\n",
    "        bad_rows.append(row)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def isFloat(val):\n",
    "    try:\n",
    "        float(val)\n",
    "        return True\n",
    "    except:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pullOutNumeric(invar):\n",
    "    if isinstance(invar,str):\n",
    "        if invar.strip(' ><').isnumeric() or isFloat(invar):\n",
    "            return float(invar.strip(' ><'))\n",
    "        else:\n",
    "            return np.nan\n",
    "    elif isinstance(invar, (int, float)):\n",
    "            return float(invar)\n",
    "    else:\n",
    "        return [float(val.strip(' ><')) for val in invar if val.strip(' ><').isnumeric() or isFloat(val.strip(' ><'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experimental_df['IC50 nums'] = experimental_df['IC50'].apply(pullOutNumeric)\n",
    "experimental_df['IC50 avg'] = experimental_df['IC50 nums'].apply(lambda x: np.nanmean(np.asarray(x, dtype=np.float32)))\n",
    "experimental_df['IC50 var'] = experimental_df['IC50 nums'].apply(lambda x: np.nanvar(np.asarray(x, dtype=np.float32)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with_sdf_data = experimental_df.copy()\n",
    "for idx, row in with_sdf_data.iterrows():\n",
    "    if not np.isnan(row['IC50 avg']):\n",
    "        continue\n",
    "    try:\n",
    "        with_sdf_data.at[idx, 'IC50 avg'] = BindingDBAffinities[(BindingDBAffinities['Receptor'] == row['Receptor']) & (BindingDBAffinities['BindingDB MonomerID'] == row['BindingDB MonomerID'])]['IC50']\n",
    "    except:\n",
    "        print(row['Receptor'],row['BindingDB MonomerID'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdf_bef_lig_check = with_sdf_data[with_sdf_data['IC50 avg'].notna()]\n",
    "sdf_ic50_vals = sdf_bef_lig_check.groupby('Receptor').filter(lambda x: len(x) > 1) # make sure receptors have more than 1 ligand\n",
    "sdf_ic50_vals['pIC50'] = sdf_ic50_vals['IC50 avg'].apply(lambda x: -np.log10(x)) #note that this is the negative log10 of the value\n",
    "print('{} unique receptors with greater than 1 ligand'.format(len(sdf_ic50_vals['Receptor'].unique())))\n",
    "sdf_ic50_vals.sort_values(by=['Receptor','file'],inplace=True)\n",
    "sdf_ic50_vals.to_csv('sdf_final_ic50_vals.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','IC50 avg', 'IC50 var','pIC50'],header=False,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "before_lig_check = experimental_df[experimental_df['IC50 avg'].notna()]\n",
    "ic50_vals = before_lig_check.groupby('Receptor').filter(lambda x: len(x) > 1) # make sure receptors have more than 1 ligand\n",
    "ic50_vals['pIC50'] = ic50_vals['IC50 avg'].apply(lambda x: -np.log10(x)) #note that this is the negative log10 of the value\n",
    "print('{} unique receptors with greater than 1 ligand'.format(len(ic50_vals['Receptor'].unique())))\n",
    "ic50_vals.sort_values(by=['Receptor','file'],inplace=True)\n",
    "ic50_vals.to_csv('final_ic50_vals.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','IC50 avg', 'IC50 var','pIC50'],header=False,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(ic50_vals.groupby('Receptor').count()['mol'],bins=np.arange(1,max(ic50_vals.groupby('Receptor').count()['mol']),1),align='mid')\n",
    "plt.hist(sdf_ic50_vals.groupby('Receptor').count()['mol'],bins=np.arange(1,max(sdf_ic50_vals.groupby('Receptor').count()['mol']),1),align='mid',label='with_sdf',alpha=0.2)\n",
    "plt.axvline(x=np.mean(ic50_vals.groupby('Receptor').count()['mol']),label='Average',color='k',linestyle='--')\n",
    "plt.xlabel('Ligands in Series')\n",
    "plt.xlim([1.5,max(experimental_df.groupby('Receptor').count()['mol'])])\n",
    "plt.ylabel('Frequency')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(ic50_vals[['Receptor','IC50 avg']].groupby('Receptor').max()['IC50 avg']-ic50_vals[['Receptor','IC50 avg']].groupby('Receptor').min()['IC50 avg'],bins=np.arange(0,10,1),align='mid')\n",
    "plt.hist(sdf_ic50_vals[['Receptor','IC50 avg']].groupby('Receptor').max()['IC50 avg']-sdf_ic50_vals[['Receptor','IC50 avg']].groupby('Receptor').min()['IC50 avg'],bins=np.arange(0,10,1),align='mid',alpha=0.2,label='with sdf')\n",
    "plt.xlabel('Range of IC50 in Series')\n",
    "plt.xlim([0,10])\n",
    "plt.ylabel('Frequency')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using just the affinity values provided in the downloaded sdfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "important_recs = experimental_df['Receptor'].unique().tolist()\n",
    "\n",
    "BindingDBAffinities = pd.DataFrame(columns=['Receptor','BindingDB MonomerID','IC50'])\n",
    "for file in glob('separated_sets/sdf_files/*.sdf'):\n",
    "    if '3D.sdf' in file:\n",
    "        continue\n",
    "    liginfo = []\n",
    "    receptor = re.findall(r'/(....)_Validation_Affinities.sdf',file)[0]\n",
    "    if receptor not in important_recs:\n",
    "        continue\n",
    "    mono_id,ic=\"\",0\n",
    "    bdbid,affdata=False,False\n",
    "    with open(file) as mol_file:\n",
    "        for line in mol_file:\n",
    "            if 'BindingDB monomerid' in line:\n",
    "                bdbid=True\n",
    "                continue\n",
    "            if bdbid:\n",
    "                mono_id = line.strip()\n",
    "                bdbid=False\n",
    "                continue\n",
    "            if 'Enzymologic: Ki nM' in line:\n",
    "                affdata=True\n",
    "                continue\n",
    "            if affdata:\n",
    "                ic = float(line)\n",
    "                affdata=False\n",
    "                continue\n",
    "            if 'Enzymologic' in line:\n",
    "                print(line)\n",
    "            if '$$$$' in line:\n",
    "                liginfo.append([receptor,mono_id,ic])\n",
    "                mono_id,ic=\"\",0\n",
    "    receptor_df = pd.DataFrame(data=liginfo,columns=['Receptor','BindingDB MonomerID','IC50'])\n",
    "    BindingDBAffinities = BindingDBAffinities.append(receptor_df,ignore_index=True)\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_of_diffs = []\n",
    "for rec in important_recs:\n",
    "    sdf_info = BindingDBAffinities[BindingDBAffinities['Receptor'] == rec]\n",
    "    tsv_info = experimental_df[(experimental_df['Receptor'] == rec)]\n",
    "    tsv_info.head()\n",
    "    for lig in sdf_info['BindingDB MonomerID'].unique().tolist():\n",
    "        sdf_val=sdf_info[sdf_info['BindingDB MonomerID']==lig]['IC50'].item()\n",
    "        tsv_stuff = tsv_info[tsv_info['BindingDB MonomerID']==lig]['IC50 avg']\n",
    "        if len(tsv_stuff) != 1:\n",
    "            print(rec)\n",
    "            print(tsv_stuff.head())\n",
    "            continue\n",
    "        tsv_val = tsv_stuff.item()\n",
    "        diff = np.abs(sdf_val-tsv_val)\n",
    "        hist_of_diffs.append(diff)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.scatter(range(len(hist_of_diffs)),hist_of_diffs)\n",
    "arraydif = np.array(hist_of_diffs)\n",
    "print(len(arraydif))\n",
    "smul = arraydif[arraydif>5]\n",
    "zeros = arraydif[arraydif<1e-5]\n",
    "plt.scatter(range(len(smul)),smul)\n",
    "print(len(zeros))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see above, there is about half that are the same between the two data types (tsv and sdf). Though the half that are not identical between the two data sources have very major discrepancies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "only_sdf_data = experimental_df.copy()\n",
    "for idx, row in only_sdf_data.iterrows():\n",
    "    only_sdf_data.at[idx,'IC50 sdf'] = BindingDBAffinities[(BindingDBAffinities['Receptor'] == row['Receptor']) & (BindingDBAffinities['BindingDB MonomerID'] == str(row['BindingDB MonomerID']))]['IC50'].item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "only_bef_lig_check = only_sdf_data[only_sdf_data['IC50 sdf'].notna()]\n",
    "sdf_only_ic50_vals = only_bef_lig_check.groupby('Receptor').filter(lambda x: len(x) > 1) # make sure receptors have more than 1 ligand\n",
    "sdf_only_ic50_vals['pIC50'] = sdf_only_ic50_vals['IC50 sdf'].apply(lambda x: -np.log10(x)) #note that this is the negative log10 of the value\n",
    "print('{} unique receptors with greater than 1 ligand'.format(len(sdf_only_ic50_vals['Receptor'].unique())))\n",
    "sdf_only_ic50_vals.sort_values(by=['Receptor','file'],inplace=True)\n",
    "sdf_only_ic50_vals.to_csv('sdf_only_final_ic50_vals.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','IC50 sdf','pIC50'],header=False,index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make the Model Training Files (Types File)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#stuff to make the training files\n",
    "import requests\n",
    "from bs4 import BeautifulSoup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def searchRCSBStructPage(s):\n",
    "    ligs = []\n",
    "    ligand_bigp = s.find(id=\"smallMoleculespanel\")\n",
    "    ligand_smallt = ligand_bigp.find(id=\"LigandsMainTable\")\n",
    "    for row in ligand_smallt.find_all(id=re.compile('ligand_row')):\n",
    "        if 'UNL' in row['id']:\n",
    "            continue\n",
    "        data = row.find_all('td')[2]\n",
    "        ligs.append(data.find_all('br')[1].next_sibling)\n",
    "    return ligs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getReferenceLigs(receptors):\n",
    "    lig_dict = dict()\n",
    "    for receptor in receptors:\n",
    "        possible_ligs =[] \n",
    "        response = requests.get('https://www.rcsb.org/structure/{}'.format(receptor))\n",
    "        soupyness = BeautifulSoup(response.content, 'html.parser')\n",
    "        try:\n",
    "            lig_dict[receptor] = searchRCSBStructPage(soupyness)\n",
    "        except:\n",
    "            try:\n",
    "                newlig= soupyness.find(id='note_obsoletedBy')\n",
    "                newreceptor = newlig.find(href=re.compile('/structure/')).text\n",
    "                response = requests.get('https://www.rcsb.org/structure/{}'.format(newreceptor))\n",
    "                soupyness = BeautifulSoup(response.content, 'html.parser')\n",
    "                lig_dict[receptor] = searchRCSBStructPage(soupyness)\n",
    "            except:\n",
    "                print(receptor)\n",
    "    return lig_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "recs = sdf_only_ic50_vals['Receptor'].unique()\n",
    "reference_dict = getReferenceLigs(recs)\n",
    "train_groups = sdf_only_ic50_vals.groupby('Receptor')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdf_only_ic50_vals['Reference'] = sdf_only_ic50_vals.apply(lambda x: x['Inchi Key'] in reference_dict[x.Receptor],axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdf_only_ic50_vals.sort_values(by=['Receptor','Reference'],ascending=[True,False],ignore_index=True,inplace=True)\n",
    "sdf_only_ic50_vals.to_csv('sdf_only_final_ic50_vals.txt', sep=' ', columns=['Receptor','file','BindingDB MonomerID','IC50 sdf','pIC50','Reference'],header=False,index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Start here if you want to use the IC<sub>50</sub> values from the sdfs downloaded from the BindingDB page and the 475 receptors (maximal amount I have ever gotten)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdf_only_ic50_vals = pd.read_csv('sdf_only_final_ic50_vals.txt', sep=' ', header=None)\n",
    "sdf_only_ic50_vals.columns = ['Receptor','file','BindingDB MonomerID','IC50 sdf','pIC50','Reference']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### maybe using combinations instead of permutations will be better (reduce it to minimum trainable set)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'combinations' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-593329a3d3af>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mrec\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgroup\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mgrouped\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m         \u001b[0mgroupinfo\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m         \u001b[0;32mfor\u001b[0m \u001b[0mrow1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrow2\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcombinations\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m##switched to combinations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m             \u001b[0mcls\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrow1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'pIC50'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrow2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'pIC50'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m             \u001b[0mdif\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrow1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'pIC50'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgroup\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mrow2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'pIC50'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'combinations' is not defined"
     ]
    }
   ],
   "source": [
    "possible_datas = [sdf_only_ic50_vals]\n",
    "trainfile_names = ['sdf_only_old_bdb_comb.types']\n",
    "for i,d in enumerate(possible_datas):\n",
    "    trainfile = pd.DataFrame(columns=['class','reg','rec','lig1','lig2'])\n",
    "    grouped = d.groupby('Receptor')\n",
    "    for rec, group in grouped:\n",
    "        groupinfo=[]\n",
    "        for row1,row2 in list(combinations(group.index,2)): ##switched to combinations\n",
    "            cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))\n",
    "            dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])\n",
    "            receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])\n",
    "            lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])\n",
    "            lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])\n",
    "            if dif == np.inf or dif == -np.inf or math.isnan(dif):\n",
    "                print(lig1,lig2,dif)\n",
    "                continue\n",
    "            groupinfo.append([cls,dif,receptor,lig1,lig2])\n",
    "        group_df = pd.DataFrame(data=groupinfo,columns=['class','reg','rec','lig1','lig2'])\n",
    "        trainfile = trainfile.append(group_df,ignore_index=True)\n",
    "    trainfile.to_csv(trainfile_names[i],header=False,index=False,sep=' ',float_format='%.4f')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the same train/test split as the DeltaDelta paper\n",
    "Using the reference ligand and increasing amounts of additional ligands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "random=45\n",
    "possible_datas = [sdf_only_ic50_vals]\n",
    "trfile_names = ['sdfo_train_papersplit_rand{}_allpos'] ##note the allpos in the name because all affinity values are positive (may not make a difference)\n",
    "ttfile_names = ['sdfo_test_papersplit_rand{}_allpos'] \n",
    "trainfile_names = [str_val.format(random) for str_val in trfile_names]\n",
    "testfile_names = [str_val.format(random) for str_val in ttfile_names]\n",
    "for i,d in enumerate(possible_datas):\n",
    "    for addnl in range(1,7):\n",
    "        trainfile = pd.DataFrame(columns=['class','reg','rec','lig1','lig2'])\n",
    "        testfile = pd.DataFrame(columns=['class','reg','rec','lig1','lig2'])\n",
    "        grouped = d.groupby('Receptor')\n",
    "        for rec, group in grouped:\n",
    "            traingroupinfo=[]\n",
    "            testgroupinfo=[]\n",
    "            if not group['Reference'].sum(): ## have to make sure there is one reference structure per group (arbitrarily assigning the reference to the first row if it isn't already assigned)\n",
    "                group.at[group.index[0],'Reference'] = True\n",
    "            if len(group)-1 < addnl:\n",
    "                train_group = pd.concat([group[(group['Reference'] == True)],group[(group['Reference'] == False)].sample(n=len(group)-1,random_state=random)])\n",
    "            else:\n",
    "                train_group = pd.concat([group[(group['Reference'] == True)],group[group['Reference'] == False].sample(n=addnl,random_state=random)])\n",
    "            ## Training set first\n",
    "            for row1,row2 in list(combinations(train_group.index,2)): ##switched to combinations\n",
    "                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])\n",
    "                if dif == np.inf or dif == -np.inf or math.isnan(dif):\n",
    "                    continue\n",
    "#                 if dif < 0: ## Lets see if all positive DDG values helps training, so if regression value < 0 then swap order\n",
    "#                     dif = np.abs(dif)\n",
    "#                     tmp = row2\n",
    "#                     row2 = row1\n",
    "#                     row1 = tmp\n",
    "                try:\n",
    "                    cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))\n",
    "                except:\n",
    "                    print(train_group)\n",
    "                    break\n",
    "                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])\n",
    "                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])\n",
    "                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])\n",
    "                traingroupinfo.append([cls,dif,receptor,lig1,lig2])\n",
    "            # Test set. Using idea that can only compare to values that are in the training set\n",
    "            for row1,row2 in list(combinations(group.index,2)):\n",
    "                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])\n",
    "                if dif == np.inf or dif == -np.inf or math.isnan(dif):\n",
    "                    continue\n",
    "#                 if dif < 0: ## Lets see if all positive DDG values helps training, so if regression value < 0 then swap order\n",
    "#                     dif = np.abs(dif)\n",
    "#                     tmp = row2\n",
    "#                     row2 = row1\n",
    "#                     row1 = tmp\n",
    "                if not (bool(row1 in train_group.index) ^ bool(row2 in train_group.index)): ##supposed to be xor\n",
    "                    continue\n",
    "                cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))\n",
    "                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])\n",
    "                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])\n",
    "                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])\n",
    "                testgroupinfo.append([cls,dif,receptor,lig1,lig2])\n",
    "            group_df = pd.DataFrame(data=traingroupinfo,columns=['class','reg','rec','lig1','lig2'])\n",
    "            trainfile = trainfile.append(group_df,ignore_index=True)\n",
    "            tgroup_df = pd.DataFrame(data=testgroupinfo,columns=['class','reg','rec','lig1','lig2'])\n",
    "            testfile = testfile.append(tgroup_df,ignore_index=True)\n",
    "        trainfile.to_csv('{}_{}.types'.format(trainfile_names[i],str(addnl)),header=False,index=False,sep=' ',float_format='%.4f')\n",
    "        testfile.to_csv('{}_{}.types'.format(testfile_names[i],str(addnl)),header=False,index=False,sep=' ',float_format='%.4f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Making 5 different train/test splits based on the random seeds\n",
    "# Additionally making two new columns before receptor for the absolute binding affinity of each ligand\n",
    "# Use permutations of the ligands rather than combinations (order matters)\n",
    "randos=[0,1,2,3,4]\n",
    "trainfile_names = [f'sdfo_train_papersplit_rand{rand}_p' for rand in randos]\n",
    "validfile_names = [f'sdfo_valid_papersplit_rand{rand}_p' for rand in randos]\n",
    "testfile_names = [f'sdfo_testv_papersplit_rand{rand}_p' for rand in randos]\n",
    "for i,random in enumerate(randos):\n",
    "    for addnl in range(1,7):\n",
    "        trainfile = pd.DataFrame(columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "        testfile = pd.DataFrame(columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "        validfile = pd.DataFrame(columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "        grouped = sdf_only_ic50_vals.groupby('Receptor')\n",
    "        for rec, group in grouped:\n",
    "            traingroupinfo=[]\n",
    "            testgroupinfo=[]\n",
    "            validgroupinfo=[]\n",
    "            if not group['Reference'].sum(): ## have to make sure there is one reference structure per group (arbitrarily assigning the reference to the first row if it isn't already assigned)\n",
    "                group.at[group.index[0],'Reference'] = True\n",
    "            if len(group)-1 < addnl:\n",
    "                train_group = pd.concat([group[(group['Reference'] == True)],group[(group['Reference'] == False)].sample(n=len(group)-1,random_state=random)])\n",
    "            else:\n",
    "                train_group = pd.concat([group[(group['Reference'] == True)],group[group['Reference'] == False].sample(n=addnl,random_state=random)])\n",
    "            not_train_group = set(group.index)-set(train_group.index)\n",
    "            if len(not_train_group) > 2:\n",
    "                valid,_ = train_test_split(list(set(group.index)-set(train_group.index)),test_size=0.6,random_state=random)\n",
    "            else:\n",
    "                valid = []\n",
    "            ## Training set first\n",
    "            for row1,row2 in list(permutations(train_group.index,2)): ##switched to combinations\n",
    "                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])\n",
    "                if dif == np.inf or dif == -np.inf or math.isnan(dif):\n",
    "                    continue\n",
    "                try:\n",
    "                    cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))\n",
    "                except:\n",
    "                    print(train_group)\n",
    "                    break\n",
    "                dg_lig1 = float(group.loc[row1,'pIC50'])\n",
    "                dg_lig2 = float(group.loc[row2,'pIC50'])\n",
    "                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])\n",
    "                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])\n",
    "                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])\n",
    "                traingroupinfo.append([cls,dif,dg_lig1,dg_lig2,receptor,lig1,lig2])\n",
    "            # Test set. Using idea that can only compare to values that are in the training set (i.e. 'Known' values)\n",
    "            for row1,row2 in list(permutations(group.index,2)):\n",
    "                if not (bool(row1 in train_group.index) ^ bool(row2 in train_group.index)): ##supposed to be xor\n",
    "                    continue\n",
    "                dif = float(group.loc[row1,'pIC50']) - float(group.loc[row2,'pIC50'])\n",
    "                if dif == np.inf or dif == -np.inf or math.isnan(dif):\n",
    "                    continue\n",
    "                dg_lig1 = float(group.loc[row1,'pIC50'])\n",
    "                dg_lig2 = float(group.loc[row2,'pIC50'])\n",
    "                cls = int(float(group.loc[row1,'pIC50']) > float(group.loc[row2,'pIC50']))\n",
    "                receptor = '{0}/{0}_R_0.gninatypes'.format(group.loc[row1,'Receptor'])\n",
    "                lig1 = '{0}/{1}_0.gninatypes'.format(group.loc[row1,'Receptor'],group.loc[row1,'file'].split('.')[0])\n",
    "                lig2 = '{0}/{1}_0.gninatypes'.format(group.loc[row2,'Receptor'],group.loc[row2,'file'].split('.')[0])\n",
    "                if row1 in valid or row2 in valid:\n",
    "                    validgroupinfo.append([cls,dif,dg_lig1,dg_lig2,receptor,lig1,lig2])\n",
    "                else:\n",
    "                    testgroupinfo.append([cls,dif,dg_lig1,dg_lig2,receptor,lig1,lig2])\n",
    "            group_df = pd.DataFrame(data=traingroupinfo,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "            trainfile = trainfile.append(group_df,ignore_index=True)\n",
    "            tgroup_df = pd.DataFrame(data=testgroupinfo,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "            vgroup_df = pd.DataFrame(data=validgroupinfo,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "            testfile = testfile.append(tgroup_df,ignore_index=True)\n",
    "            validfile = validfile.append(vgroup_df,ignore_index=True)\n",
    "        trainfile.to_csv(f'CV_Sets/{trainfile_names[i]}_{addnl}.types',header=False,index=False,sep=' ',float_format='%.4f')\n",
    "        testfile.to_csv(f'CV_Sets/{testfile_names[i]}_{addnl}.types',header=False,index=False,sep=' ',float_format='%.4f')\n",
    "        validfile.to_csv(f'CV_Sets/{validfile_names[i]}_{addnl}.types',header=False,index=False,sep=' ',float_format='%.4f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "group_rec = sdf_only_ic50_vals.groupby('Receptor')\n",
    "full_types_list=[]\n",
    "for rec_name, group in group_rec:\n",
    "    for idx1, idx2 in list(combinations(group.index,2)):\n",
    "        assert group.loc[idx1,'Receptor'] == group.loc[idx2,'Receptor']\n",
    "        regression = float(group.loc[idx1,'pIC50']) - float(group.loc[idx2,'pIC50'])\n",
    "        if regression == np.inf or regression == -np.inf or math.isnan(regression):\n",
    "                    continue\n",
    "#         if regression < 0: ## Lets see if all positive DDG values helps training, so if regression value < 0 then swap order\n",
    "#             regression = np.abs(regression)\n",
    "#             tmp = idx2\n",
    "#             idx2 = idx1\n",
    "#             idx1 = tmp\n",
    "        try:\n",
    "            classification = int(float(group.loc[idx1,'pIC50']) > float(group.loc[idx2,'pIC50']))\n",
    "        except:\n",
    "            print(train_group)\n",
    "            break\n",
    "        \n",
    "        receptor = '{0}/{0}_R.pdb'.format(group.loc[idx1,'Receptor'])\n",
    "        lig1 = '{0}/{1}'.format(group.loc[idx1,'Receptor'],group.loc[idx1,'file'])\n",
    "        lig2 = '{0}/{1}'.format(group.loc[idx2,'Receptor'],group.loc[idx2,'file'])\n",
    "        full_types_list.append([classification, regression, receptor,lig1,lig2])\n",
    "full_types_df = pd.DataFrame(full_types_list,columns=['Classification','∆pIC50','Rec','Lig1','Lig2'])\n",
    "full_types_df.to_csv('all_bdb_types.types',sep=' ',header=False,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uses permutations\n",
    "# Outputs cls, reg, absaff1, absaff2, rec, lig1, lig2\n",
    "group_rec = sdf_only_ic50_vals.groupby('Receptor')\n",
    "full_types_list=[]\n",
    "for rec_name, group in group_rec:\n",
    "    for idx1, idx2 in list(permutations(group.index,2)):\n",
    "        assert group.loc[idx1,'Receptor'] == group.loc[idx2,'Receptor']\n",
    "        regression = float(group.loc[idx1,'pIC50']) - float(group.loc[idx2,'pIC50'])\n",
    "        if regression == np.inf or regression == -np.inf or math.isnan(regression):\n",
    "                    continue\n",
    "        try:\n",
    "            classification = int(float(group.loc[idx1,'pIC50']) > float(group.loc[idx2,'pIC50']))\n",
    "        except:\n",
    "            print(train_group)\n",
    "            break\n",
    "        dg_lig1 = float(group.loc[idx1,'pIC50'])\n",
    "        dg_lig2 = float(group.loc[idx2,'pIC50'])\n",
    "        receptor = '{0}/{0}_R.pdb'.format(group.loc[idx1,'Receptor'])\n",
    "        lig1 = '{0}/{1}'.format(group.loc[idx1,'Receptor'],group.loc[idx1,'file'])\n",
    "        lig2 = '{0}/{1}'.format(group.loc[idx2,'Receptor'],group.loc[idx2,'file'])\n",
    "        full_types_list.append([classification, regression, dg_lig1, dg_lig2, receptor,lig1,lig2])\n",
    "full_types_df_best = pd.DataFrame(full_types_list,columns=['class','reg','dg_lig1','dg_lig2','rec','lig1','lig2'])\n",
    "full_types_df_best.to_csv('all_bdb_types_perm_mult.types',sep=' ',header=False,index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}