ccbm_model_quality.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating figures and analysis for assessing the quality of the cbmlib models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download Smiles structure representations from ChEMBL\n",
    "... You could also simply use your browser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve: ftp://ftp.ebi.ac.uk/pub/databases/chembl/ChEMBLdb/releases/chembl_25/chembl_25_chemreps.txt.gz\n",
    "from ftplib import FTP\n",
    "\n",
    "ftp = FTP('ftp.ebi.ac.uk')     # connect to host, default port\n",
    "ftp.login()\n",
    "ftp.cwd('pub/databases/chembl/ChEMBLdb/releases/chembl_25')\n",
    "with open('chembl_25_chemreps.txt.gz', 'wb') as fout:\n",
    "    ftp.retrbinary('RETR chembl_25_chemreps.txt.gz', fout.write)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decompress data\n",
    "import gzip\n",
    "with  gzip.open(\"chembl_25_chemreps.txt.gz\",\"rt\") as fin:\n",
    "    with open(\"chembl_25_chemreps.txt\",\"wt\") as fout:\n",
    "        for line in fin:\n",
    "            fout.write(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wash the Smiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ccbmlib.preprocessing import wash_molecules,export_washed\n",
    "import rdkit.Chem as Chem\n",
    "import os\n",
    "\n",
    "mol_suppl = Chem.SmilesMolSupplier(\"chembl_25_chemreps.txt\",delimiter=\"\\t\",smilesColumn=1,nameColumn=0)\n",
    "if not os.path.exists(\"chembl25.smi\"):\n",
    "    washed = wash_molecules(mol_suppl)\n",
    "    export_washed(washed, \"chembl25.smi\", \"chembl25_duplicates.txt\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn on logging as genrating staitstics can be very time consuming\n",
    "\n",
    "import logging\n",
    "logging.basicConfig()\n",
    "logging.getLogger().setLevel(logging.DEBUG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import ccbmlib.models as ccbm\n",
    "\n",
    "# Make sure a data folder is specified!\n",
    "data_folder = \"ccbm_models\"\n",
    "os.makedirs(data_folder, exist_ok=True)\n",
    "ccbm.set_data_folder(data_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup different fingerprints\n",
    "fp_names = [ \"atom_pairs\", \"avalon\", \"maccs\", \"morgan\", \"morgan\", \"torsions\", \n",
    "        \"hashed_atom_pairs\",\"hashed_morgan\",\"hashed_morgan\",\"hashed_torsions\", \"rdkit\"]\n",
    "titles = [\"Atom pairs\",\"Avalon\",\"MACCS\",\"Morgan, radius 1\",\"Morgan, radius 2\",\"Topological torsions\",\"Atom pairs, hashed\",\n",
    "         \"Morgan, radius 1, hashed\",\"Morgan, radius 2, hashed\", \"Topological torsions, hashed\",\"RDKit\" ]\n",
    "pars = [{} for _ in range(11)]\n",
    "pars[3] = pars[7] = {\"radius\":1}\n",
    "pars[4] = pars[8] = {\"radius\":2}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Only a sample ...\n",
    "Generating statistics for 1.7 mio compounds is time consuming. You might want to try the code on a small sample first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def skip_some(seq_like,step):\n",
    "    for i,x in enumerate(seq_like):\n",
    "        if i%step==0:\n",
    "            yield x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"chembl25.smi\") as fin:\n",
    "    with open(\"chembl25_sample.smi\",\"w\") as fout:\n",
    "        for line in skip_some(fin,100):\n",
    "            fout.write(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "db_name = \"chembl25_sample\"\n",
    "filename = \"chembl25_sample.smi\"\n",
    "\n",
    "#uncomment to run on full data set\n",
    "# db_name = \"chembl25\"\n",
    "# filename = \"chembl25.smi\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats=[]\n",
    "for fp_name,fp_par in zip(fp_names,pars):\n",
    "    pwstats = ccbm.get_feature_statistics(db_name, fp_name, fp_par, filename)\n",
    "    stats.append(pwstats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.WARN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code for Model evaluation\n",
    "\n",
    "For random sampling the code will tr to keep all fingerprints in memory.\n",
    "This can be very taxing on low memory machines especially for fingerprints with a large number of features present per fingerpint\n",
    "like atom pairs or the rdkit fingerprint. \n",
    "\n",
    "Modify the code below for reading fingerprints using skip_some to reduce the memory load. The code below illustrates how to do it\n",
    "for the rdkit fingerprint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fingerprint reader\n",
    "def fingerprint_source(fname):\n",
    "    with ccbm.auto_open(fname) as f:\n",
    "        next(f)\n",
    "        for line in f:\n",
    "            if line.startswith(\"#\"):\n",
    "                continue\n",
    "            entries = line.strip(\"\\n\").split(\"\\t\")\n",
    "            fp = set(map(int,entries[1].split()))\n",
    "            yield fp\n",
    "\n",
    "            # Calculate Tanimoto coefficient\n",
    "def tc(A,B):\n",
    "    a=len(A)\n",
    "    b=len(B)\n",
    "    c=len(A.intersection(B))\n",
    "    if c>0:\n",
    "        return c/(a+b-c)\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "# Generate random Tc's\n",
    "# if given one reference fingerprint stays fixed\n",
    "def random_tcs(fps,n=100000,ref=None):\n",
    "    tcs=[]\n",
    "    if ref:\n",
    "        for _ in range(n):\n",
    "            a = random.choice(fps)\n",
    "            tcs.append(tc(ref,a))\n",
    "    else:\n",
    "        for _ in range(n):\n",
    "            a = random.choice(fps)\n",
    "            b = random.choice(fps)\n",
    "            tcs.append(tc(a,b))\n",
    "    tcs.sort()\n",
    "    return tcs\n",
    "\n",
    "# Retrieve empirical CDF from sample of Tcs\n",
    "# using 'num_data_pts' data points\n",
    "def get_cdf_from_sample(tcs,num_data_pts=1000):\n",
    "    tc_x = []\n",
    "    hist_y = []\n",
    "    n=len(tcs)\n",
    "    ct = 0\n",
    "    next_idx = 0\n",
    "    for i,v in enumerate(tcs):\n",
    "        if i>=next_idx:\n",
    "            tc_x.append(v)\n",
    "            hist_y.append(i/n)\n",
    "            ct +=1\n",
    "            next_idx = n*ct/num_data_pts\n",
    "    tc_x.append(tcs[-1])\n",
    "    hist_y.append(1)\n",
    "    tc_x.append(1)\n",
    "    hist_y.append(1)\n",
    "    return tc_x, hist_y\n",
    "\n",
    "# Retrieve modeled CDF at values given in 'tcs_x'\n",
    "def get_cdf_from_model(tcs_x,cnd):\n",
    "    return [cnd.cdf(x) for x in tcs_x]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1.inset_locator import (inset_axes, InsetPosition,\n",
    "                                                  mark_inset)\n",
    "\n",
    "dpi = 300\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n=1000000\n",
    "\n",
    "outf = open(db_name+\"-KS-results.txt\",\"w\")\n",
    "for fp_name,par,title in zip(fp_names[:],pars[:],titles):\n",
    "    title0=fp_name\n",
    "    if \"radius\" in par:\n",
    "        title0+=str(par[\"radius\"])\n",
    "    img_filename=db_name+\"-global_model-cdf-\"+title0+\".png\"\n",
    "    fp_file=ccbm.get_full_filename(ccbm.get_fingerprints(db_name,fp_name,par))\n",
    "    fps=[]\n",
    "    \n",
    "    # Modify by demand using skip_some to reduce memory load\n",
    "    fps =list(fingerprint_source(fp_file)) if fp_name!=\"rdkit\" else list(skip_some(fingerprint_source(fp_file),3)) \n",
    "\n",
    "    pw = ccbm.get_feature_statistics(db_name,fp_name,par)\n",
    "    tcs=random_tcs(fps,n)\n",
    "    cnd=pw.get_tc_distribution()\n",
    "    tc_x,hist_y = get_cdf_from_sample(tcs)\n",
    "    cnd_y = get_cdf_from_model(tc_x,cnd)\n",
    "    \n",
    "    deltas=[abs(a-b) for a,b in zip(hist_y,cnd_y)]\n",
    "    for i,v in enumerate(hist_y):\n",
    "        if v>0.9:\n",
    "            break\n",
    "    KS_all=max(deltas)\n",
    "    deltas=deltas[i:]\n",
    "    KS_90=max(deltas)\n",
    "\n",
    "    print(title,\"    KS_all: {} KS_90: {}\".format(KS_all,KS_90),file=outf)\n",
    "    print(title,\"    KS_all: {} KS_90: {}\".format(KS_all,KS_90))\n",
    "\n",
    "    fig, ax1 = plt.subplots(figsize=(8*0.75,6*0.75))\n",
    "\n",
    "    # These are in unitless percentages of the figure size. (0,0 is bottom left)\n",
    "    left, bottom, width, height = [0.53, 0.18, 0.35, 0.35]\n",
    "    ax2 = fig.add_axes([left, bottom, width, height])\n",
    "    mark_inset(ax1, ax2, loc1=1, loc2=3, fc=\"none\", ec='0.5')\n",
    "    axs=[ax1,ax2]\n",
    "    plt.sca(axs[0])\n",
    "    plt.plot(tc_x,hist_y,'#4080c0')\n",
    "    plt.plot(tc_x,cnd_y,'#d03020')\n",
    "    plt.gca().set_xlim(0,1)\n",
    "    plt.sca(axs[1])\n",
    "    \n",
    "    plt.gca().patch.set_facecolor('w')\n",
    "    plt.gca().patch.set_alpha(0)\n",
    "    plt.plot(tc_x[i:-2],hist_y[i:-2],'#4080c0')\n",
    "    plt.plot(tc_x[i:-2],cnd_y[i:-2],'#d03020')\n",
    "\n",
    "    plt.sca(ax1)\n",
    "    plt.title(title,fontsize=18)\n",
    "    plt.xlabel(\"Tc\",fontsize=16)\n",
    "    plt.ylabel(\"Significance\",fontsize=16)\n",
    "    plt.gca().tick_params(axis=\"both\",which=\"major\",labelsize=14)\n",
    "\n",
    "    yticks = [\"{:.3}\".format(x) for x in ax2.get_yticks()]\n",
    "    ax2.set_yticklabels(yticks, backgroundcolor='w')\n",
    "    plt.savefig(img_filename,dpi=dpi)\n",
    "    plt.show()\n",
    "    #break\n",
    "\n",
    "outf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate figures for conditional distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import numpy as np\n",
    "\n",
    "logger = logging.getLogger(\"FpDistModels\")\n",
    "logger.setLevel(logging.WARN)\n",
    "\n",
    "n=100000\n",
    "\n",
    "for fp_name,par,title in zip(fp_names,pars,titles):\n",
    "    title0=fp_name\n",
    "    if \"radius\" in par:\n",
    "        title0+=str(par[\"radius\"])\n",
    "    img_filename=db_name+\"-conditional_model-sig2sig-\"+title0+\".png\"\n",
    "\n",
    "    fp_file=ccbm.get_full_filename(ccbm.get_fingerprints(db_name,fp_name,par))\n",
    "    fps=[]\n",
    "    cnds=[]\n",
    "    \n",
    "    # Modify by demand using skip_some to reduce memory load\n",
    "    fps = list(fingerprint_source(fp_file)) if fp_name!=\"rdkit\" else list(skip_some(fingerprint_source(fp_file),3))\n",
    "\n",
    "    pw = ccbm.get_feature_statistics(db_name,fp_name,par)\n",
    "    cnds = []\n",
    "    for _ in range(100):\n",
    "        ref = random.choice(fps)\n",
    "        tcs=random_tcs(fps,n,ref)\n",
    "        cnd=pw.get_tc_distribution(ref)\n",
    "        tc_x,hist_y = get_cdf_from_sample(tcs)\n",
    "        cnd_y = get_cdf_from_model(tc_x,cnd)\n",
    "        cnds.append(cnd_y)\n",
    "    bp=[]\n",
    "    for v in zip (*cnds):\n",
    "        np.percentile(v,[0,5,25,50,75,95,100])\n",
    "        bp.append(np.percentile(v,[0,5,25,50,75,95,100]))\n",
    "    bp=list(zip(*bp))\n",
    "    \n",
    "    fig, ax1 = plt.subplots(figsize=(8*0.75,6*0.75))\n",
    "\n",
    "    #These are in unitless percentages of the figure size. (0,0 is bottom left)\n",
    "    left, bottom, width, height = [0.58, 0.18, 0.3, 0.3]\n",
    "    ax2 = fig.add_axes([left, bottom, width, height])\n",
    "    mark_inset(ax1, ax2, loc1=2, loc2=4, fc=\"none\", ec='0.5')\n",
    "    axs=[ax1,ax2]\n",
    "    plt.sca(axs[0])\n",
    "    plt.title(title,fontsize=18)\n",
    "    colors = ['#c0c0c0','#808080','#404040','#000000','#404040','#808080','#c0c0c0']\n",
    "    #plt.fill_between(hist_y,bp[0],bp[-1],color='#e0e0e0')\n",
    "    plt.fill_between(hist_y,bp[1],bp[-2],color='#c0c0c0')\n",
    "    plt.fill_between(hist_y,bp[2],bp[-3],color='#808080')\n",
    "    plt.plot(hist_y,hist_y,'#00d060')\n",
    "    plt.plot(hist_y,bp[3],'#000000')\n",
    "    \n",
    "    plt.sca(axs[1])\n",
    "    for i,v in enumerate(hist_y):\n",
    "        if v>=0.9: break\n",
    "    #plt.fill_between(hist_y,bp[0],bp[-1],color='#e0e0e0')\n",
    "    plt.fill_between(hist_y[i:],bp[1][i:],bp[-2][i:],color='#c0c0c0')\n",
    "    plt.fill_between(hist_y[i:],bp[2][i:],bp[-3][i:],color='#808080')\n",
    "    plt.plot(hist_y[i:],hist_y[i:],'#00d060')\n",
    "    plt.plot(hist_y[i:],bp[3][i:],'#000000')\n",
    "    \n",
    "    plt.sca(axs[0])\n",
    "    plt.xlabel(\"Empirical significance\",fontsize=16)\n",
    "    plt.ylabel(\"Modeled significance\",fontsize=16)\n",
    "    plt.gca().tick_params(axis=\"both\",which=\"major\",labelsize=14)\n",
    "    plt.savefig(img_filename,dpi=dpi)\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some simple fingerprint statistics\n",
    "- Average number of features per fingerprint\n",
    "- standard deviation of the number of features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "for fp_name,par,title in zip(fp_names[:],pars[:],titles):\n",
    "    pw = ccbm.get_feature_statistics(db_name,fp_name,par)  \n",
    "    print(title,\"{:8.1f} {:8.1f}\".format(pw.marginal.avg_no_of_features,sqrt(pw.marginal.no_of_features_variance)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relating Tanimoto coefficients of different fingerprints to each other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse(f,s):\n",
    "    tl = -5\n",
    "    th = 5\n",
    "    sl = f(tl)\n",
    "    sh = f(th)\n",
    "    if s <= sl:\n",
    "        return tl\n",
    "    if s >= sh:\n",
    "        return th\n",
    "    while th - tl > 1e-6:\n",
    "        tm = (th+tl)/2\n",
    "        sm = f(tm)\n",
    "        if sm < s:\n",
    "            tl = tm\n",
    "            sl = sm\n",
    "        elif sm > s:\n",
    "            th = tm\n",
    "            sh = sm\n",
    "        else:\n",
    "            return tm\n",
    "    alpha = (s-sl)/(sh-sl)\n",
    "    return (1-alpha)*tl + alpha*th\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "img_filename=db_name+\"-maccs_tc_to_tc.png\"\n",
    "fig, ax1 = plt.subplots(figsize=(8*0.75,6*0.75))\n",
    "plt.sca(ax1)\n",
    "ref=2\n",
    "sel = [1,2,3,4,5,10]\n",
    "sel.remove(ref)\n",
    "fps = [list(zip(fp_names,pars,stats))[i] for i in sel]\n",
    "modelx = stats[ref].get_tc_distribution()\n",
    "for i in sel:\n",
    "    x = np.linspace(0,1,1001)\n",
    "    modely = stats[i].get_tc_distribution()\n",
    "    pvalues = [modelx.cdf(t) for t in x]\n",
    "    y = [inverse(modely.cdf,p) for p in pvalues]\n",
    "    #x,y = zip(*[(x,y) for x,y,p in zip(x,y,pvalues) if p>=0.01 and p<=0.999])\n",
    "    plt.plot(x,y,label=titles[i])\n",
    "plt.legend(fontsize=14)\n",
    "plt.xlabel(\"MACCS Tc\",fontsize=16)\n",
    "plt.ylabel(\"Tc\",fontsize=16)\n",
    "plt.gca().tick_params(axis=\"both\",which=\"major\",labelsize=14)\n",
    "xupp=modelx.icdf(0.99)\n",
    "print(xupp)\n",
    "r=plt.Rectangle((xupp,0),1-xupp,1,fill=True,color=\"w\",zorder=2.5,alpha=0.5)\n",
    "plt.gca().add_patch(r)\n",
    "plt.plot([xupp,xupp],[0,1],\"--\",color=\"#808080\")\n",
    "plt.gca().set_xlim(0,1)\n",
    "plt.gca().set_ylim(0,1)\n",
    "plt.savefig(img_filename,dpi=dpi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Corresponding Tc values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Corresponding Tc\")\n",
    "sig=stats[ref].get_cached_tc_distribution().cdf(0.6)\n",
    "for i in sel+[ref]:\n",
    "    print(titles[i],\"Tc: %.2f\"%inverse(stats[i].get_tc_distribution().cdf,sig))"
   ]
  },
  {
   "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.4"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}