b2ar_bias_analysis-HTBC-Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "In this iPython notebook, we will featurize MOR ligand binding simulation by pairwise distances between the ligand and different receptor residues. We will then perform tICA and prospectively build an MSM. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# changing matplotlib the default style\n",
    "matplotlib.style.use('ggplot')\n",
    "#matplotlib.rcParams[\"figure.facecolor\"] = \"white\"\n",
    "#matplotlib.rcPar|ams[\"savefig.transparent\"] = \"True\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/enf/anaconda3/lib/python3.5/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect\n",
      "because the backend has already been chosen;\n",
      "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "we are operating on biox3\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "\n",
    "from PDB_Order_Fixer import PDB_Order_Fixer\n",
    "import mdtraj as md\n",
    "import os\n",
    "import numpy as np\n",
    "import h5py\n",
    "\n",
    "import datetime\n",
    "import glob\n",
    "import copy\n",
    "from functools import partial \n",
    "import operator\n",
    "import time\n",
    "\n",
    "import random \n",
    "import subprocess\n",
    "from subprocess import Popen\n",
    "import sys\n",
    "from custom_clusterer import *\n",
    "from custom_tica import *\n",
    "from custom_featurizer import *\n",
    "from pdb_editing import *\n",
    "from analysis import *\n",
    "from io_functions import *\n",
    "#from topology_fixing import *\n",
    "from subsampling import *\n",
    "from conversions import *\n",
    "from custom_msm import *\n",
    "from grids import *\n",
    "from docking_analysis import *\n",
    "\n",
    "from scipy import stats\n",
    "import os\n",
    "from efficacy_scripts import *\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from detect_intermediates import *\n",
    "from interpret_tICs import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from msmbuilder.utils import verbosedump, verboseload\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from b2ar_feature_types import *\n",
    "#from b2ar_feature_types import *\n",
    "from get_variable_names import *\n",
    "from b2ar_tica_config import *\n",
    "from residue import Residue, Atom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn import tree\n",
    "from sklearn.tree import DecisionTreeClassifier, export_graphviz\n",
    "from sklearn.preprocessing import scale\n",
    "from random import shuffle\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import roc_curve, auc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "ori_feature_name = copy.deepcopy(feature_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#schemes = [\"closest-heavy\", \"CA\"]\n",
    "#feature_name = \"%s-CA\" %ori_feature_name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rho = 0.01\n",
    "rho_string = \"_rho0pt01\"\n",
    "n_clusters = 25\n",
    "n_samples = 1\n",
    "lag_time=5\n",
    "precision = \"XP\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A\n",
      "/home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A\n"
     ]
    }
   ],
   "source": [
    "(active_ref_dir, inactive_ref_dir, simulation_ref_dir, scripts_dir,\n",
    "          ligand_dir, agonist_dir, inverse_agonist_dir, biased_agonist_dir, ref_receptors_dir, whole_trajectory_pnas,\n",
    "          sasa_file) = get_base_files(base)\n",
    "\n",
    "tica_dir = get_tica_dir(base, is_sparse, lag_time, n_components, feature_name, \n",
    "                                 wolf_string, shrinkage_string, rho_string)\n",
    "ori_tica_dir = copy.deepcopy(tica_dir)\n",
    "features_dir = get_features_dir(base, feature_name)\n",
    "\n",
    "landmarks_dir = get_landmarks_dir(tica_dir)\n",
    "analysis_dir = get_analysis_dir(tica_dir, n_clusters, sampling_method)\n",
    "gmm_dir = get_gmm_dir(tica_dir)\n",
    "rf_dirdir = get_rf_dir(tica_dir)\n",
    "\n",
    "\n",
    "ref_tica_dir, ref_tica_coords = get_ref_tica_dirs(tica_dir)\n",
    "\n",
    "graph_file = get_graph_file(tica_dir, msm_lag_time, n_clusters)\n",
    "\n",
    "pnas_titles =  [\"tm6_tm3_dist\", \"rmsd_npxxy_inactive\", \"rmsd_npxxy_active\", \"rmsd_connector_inactive\", \"rmsd_connector_active\"]\n",
    "pnas_features_dir = analysis_dir\n",
    "\n",
    "\n",
    "(clusterer_dir, msm_model_dir, macrostate_dir, features_known, model_dir, projected_features_dir,\n",
    "         projection_operator_dir, ktica_fit_model_filename, ktica_projected_data_filename, nystroem_data_filename,\n",
    "         mutual_information_csv, pearson_csv) = get_tica_files(base, tica_dir, n_clusters, msm_lag_time, n_macrostates)\n",
    "\n",
    "(standardized_features_dir, feature_residues_csv, feature_residues_pkl,\n",
    "          contact_csv, ref_features_dir) = get_feature_files(features_dir)\n",
    "\n",
    "(kmeans_csv, tica_coords_csv, features_csv, active_rmsd_dir, inactive_rmsd_dir, active_pnas_dir, inactive_pnas_joined, active_pnas_joined,\n",
    "        clusters_map_file, ktica_clusters_map_file, analysis_file, combined_file, docking_summary, docking_joined, docking_z_scores_csv,\n",
    "        aggregate_docking, aggregate_docking_joined, docking_pnas_joined, aggregate_docking_pnas, aggregate_docking_pnas_joined, docking_multiple_ligands,\n",
    "        docking_distances_file, docking_pdf, mmgbsa_docking_distances, pnas_coords, mmgbsa_dir, mmgbsa_csv, mmgbsa_pdf, aggregate_mmgbsa,\n",
    "        aggregate_mmgbsa_joined, aggregate_mmgbsa_pnas_joined, mmgbsa_z_scores_csv, active_clusters_csv, intermediate_clusters_csv,\n",
    "        inactive_clusters_csv, pnas_clusters_averages, tica_clusters_averages, tica_classes_csv, tica_samples_csv, subgraph_save_base,\n",
    "        degree_save_base, degree_map_csv, degree_z_map_csv, aggregate_docking_pnas_degree_z_joined, tic_residue_csv, feature_coefs_csv,\n",
    "        duplicated_feature_coefs_csv) = get_analysis_files(analysis_dir, n_clusters, tica_dir, tica_dir, sampling_method, n_samples, precision,\n",
    "                                                           msm_lag_time)\n",
    "\n",
    "(inactive_pnas_distances_dir, active_pnas_distances_dir, active_pnas_all_distances_dir,\n",
    "          inactive_pnas_distances_new_csv, active_pnas_distances_new_csv, active_pnas_joined, active_pnas_means, pnas_coords_dir,\n",
    "          pnas_coords_csv, pnas_all_coords_csv, pnas_coords_hexbin_dir, pnas_coords_co_crystallized_docking_dir,\n",
    "          pnas_coords_active_colors_dir, user_defined_features_file, reaction_coordinates_trajs_file) = get_pnas_files(whole_trajectory_pnas, pnas_features_dir)\n",
    "\n",
    "features_dir = get_features_dir(base, feature_name)\n",
    "\n",
    "\n",
    "\n",
    "graph_file = get_graph_file(tica_dir, msm_lag_time, n_clusters)\n",
    "(scripts_dir, pymol_fixpdb_dir) = get_script_dir(scripts_dir)\n",
    "(save_dir, reimaged_dir, mae_dir, combined_reimaged_dir, grid_dir, docking_dir) = get_docking_dirs(tica_dir, n_clusters, n_components, n_samples, sampling_method, precision)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import multiprocessing as mp\n",
    "mp.cpu_count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from imp import reload\n",
    "import analysis\n",
    "reload(analysis)\n",
    "from analysis import *\n",
    "\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "\n",
    "reference_docking_dir = \"/home/enf/b2ar_analysis/reference_docking/docking_SP\"\n",
    "ligands_dir = \"/home/enf/htbc/sdfs\"\n",
    "reference_docking_df, reference_poses_df = analyze_docking_results_in_dir(reference_docking_dir, ligands_dir, write_to_disk=False, redo=False)\n",
    "#reference_docking_df = analyze_docking_results_in_dir(reference_docking_dir, ligands_dir, write_to_disk=False, redo=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['2rh1_grid', '3p0g_grid']"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reference_docking_df.columns.values.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = reference_docking_df.sort(columns='3p0g_grid', inplace=False, ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df[\"mean\"] = df.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>2rh1_grid</th>\n",
       "      <th>3p0g_grid</th>\n",
       "      <th>mean</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>_46365</th>\n",
       "      <td>11.15</td>\n",
       "      <td>10.35</td>\n",
       "      <td>10.750</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_43235</th>\n",
       "      <td>10.87</td>\n",
       "      <td>9.70</td>\n",
       "      <td>10.285</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_62901</th>\n",
       "      <td>10.67</td>\n",
       "      <td>8.85</td>\n",
       "      <td>9.760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_119833</th>\n",
       "      <td>8.53</td>\n",
       "      <td>10.85</td>\n",
       "      <td>9.690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_103698</th>\n",
       "      <td>10.82</td>\n",
       "      <td>8.49</td>\n",
       "      <td>9.655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_42622</th>\n",
       "      <td>10.83</td>\n",
       "      <td>8.19</td>\n",
       "      <td>9.510</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_103874</th>\n",
       "      <td>10.10</td>\n",
       "      <td>8.91</td>\n",
       "      <td>9.505</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_124773</th>\n",
       "      <td>8.52</td>\n",
       "      <td>10.49</td>\n",
       "      <td>9.505</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_124756</th>\n",
       "      <td>9.39</td>\n",
       "      <td>9.60</td>\n",
       "      <td>9.495</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>_120850</th>\n",
       "      <td>9.38</td>\n",
       "      <td>9.51</td>\n",
       "      <td>9.445</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         2rh1_grid  3p0g_grid    mean\n",
       "_46365       11.15      10.35  10.750\n",
       "_43235       10.87       9.70  10.285\n",
       "_62901       10.67       8.85   9.760\n",
       "_119833       8.53      10.85   9.690\n",
       "_103698      10.82       8.49   9.655\n",
       "_42622       10.83       8.19   9.510\n",
       "_103874      10.10       8.91   9.505\n",
       "_124773       8.52      10.49   9.505\n",
       "_124756       9.39       9.60   9.495\n",
       "_120850       9.38       9.51   9.445"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = df.sort(\"mean\", ascending=False, inplace=False)\n",
    "df.iloc[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "338\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<AsyncMapResult: chdir>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from ipyparallel import Client\n",
    "rc = Client()\n",
    "print(len(rc.ids))\n",
    "dview = rc[:]\n",
    "dview.map(os.chdir, ['/home/enf/b2ar_analysis/conformation']*len(rc.ids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import grids\n",
    "reload(grids)\n",
    "from grids import *\n",
    "\n",
    "grid_center = \"64.4, 16.9, 11.99\"\n",
    "\n",
    "\n",
    "#reimaged_dir = samples_dir\n",
    "#mae_dir = reimaged_dir\n",
    "#remove_ter(reimaged_dir)\n",
    "#reorder(reimaged_dir)\n",
    "\n",
    "\n",
    "#agonist_ligands = [a for a in agonist_ligands if \"TA\" not in a]\n",
    "\n",
    "#pprep(mae_dir, ref = active_ref_dir, chosen_receptors = chosen_receptors, worker_pool=None, parallel=True)\n",
    "#generate_grids(mae_dir, grid_center, grid_dir, remove_lig = \"BIA\", chosen_receptors = chosen_receptors, worker_pool=None, outer_box=25.)\n",
    "\n",
    "precision = \"XP\"\n",
    "htbc_dir = \"/home/enf/htbc/sdfs\"\n",
    "docking_dir = \"/home/enf/htbc/b2ar/docking_xp\"\n",
    "grid_dir = \"/home/enf/htbc/b2ar/grids\"\n",
    "\n",
    "if not os.path.exists(docking_dir):\n",
    "    os.makedirs(docking_dir)\n",
    "if not os.path.exists(grid_dir):\n",
    "    os.makedirs(grid_dir)\n",
    "\n",
    "chosen_ligands = df.index.values.tolist()[:1000]\n",
    "#dock_ligands_and_receptors(grid_dir, docking_dir, htbc_dir, precision = precision, ext = \"-out.maegz\", chosen_ligands=chosen_ligands, chosen_receptors = False, parallel = False, grid_ext = \".zip\", worker_pool=dview)\n",
    "\n",
    "#dock_ligands_and_receptors(grid_dir, docking_dir,  biased_agonist_dir, precision = precision, ext = \"-out.maegz\", chosen_ligands = biased_ligands, chosen_receptors = chosen_receptors, parallel = None, grid_ext = \".grd\", worker_pool=dview)\n",
    "#dock_ligands_and_receptors(grid_dir, docking_dir, agonist_dir, precision = precision, ext = \"-out.maegz\", chosen_ligands = agonist_ligands, chosen_receptors = chosen_receptors, parallel = None, grid_ext = \".grd\", worker_pool=dview)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 729.62933227  289.65120394]\n"
     ]
    }
   ],
   "source": [
    "tica_object = compat_verboseload(projection_operator_dir)\n",
    "print(tica_object.timescales_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Interpreting tIC 1\n",
      "feature_importances_df.shape\n",
      "(6, 5)\n",
      "residue_importances_df.shape\n",
      "(11, 2)\n",
      "       feature_name   res_i   res_j  importance           feature\n",
      "3  Phe108 To Gln170  Phe108  Gln170    0.278917  PHE108 to GLN170\n",
      "5  Thr274 To Cys327  Thr274  Cys327   -0.216950  THR274 to CYS327\n",
      "0   Asn51 To Ser319   Asn51  Ser319    0.209493   ASN51 to SER319\n",
      "4   Asp79 To Tyr326   Asp79  Tyr326    0.195303   ASP79 to TYR326\n",
      "2   Leu75 To Pro323   Leu75  Pro323   -0.172376   LEU75 to PRO323\n",
      "1   Asn51 To Gly320   Asn51  Gly320   -0.163539   ASN51 to GLY320\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Gln170  Gln170    0.278917\n",
      "Phe108  Phe108    0.278917\n",
      "Cys327  Cys327   -0.216950\n",
      "Thr274  Thr274   -0.216950\n",
      "Ser319  Ser319    0.209493\n",
      "Tyr326  Tyr326    0.195303\n",
      "Asp79    Asp79    0.195303\n",
      "Asn51    Asn51    0.190841\n",
      "Pro323  Pro323   -0.172376\n",
      "Leu75    Leu75   -0.172376\n",
      "Gly320  Gly320   -0.163539\n",
      "Using dark_background\n",
      "Interpreting tIC 2\n",
      "feature_importances_df.shape\n",
      "(7, 5)\n",
      "residue_importances_df.shape\n",
      "(12, 2)\n",
      "       feature_name   res_i   res_j  importance           feature\n",
      "6  Tyr219 To Leu272  Tyr219  Leu272    0.442840  TYR219 to LEU272\n",
      "2   Asp79 To Tyr326   Asp79  Tyr326    0.325170   ASP79 to TYR326\n",
      "3  Ile127 To Leu275  Ile127  Leu275    0.305118  ILE127 to LEU275\n",
      "4  Tyr141 To Glu268  Tyr141  Glu268    0.272796  TYR141 to GLU268\n",
      "0   Asp79 To Asn318   Asp79  Asn318    0.257021   ASP79 to ASN318\n",
      "5  Tyr219 To Ala271  Tyr219  Ala271   -0.205418  TYR219 to ALA271\n",
      "1   His93 To Cys190   His93  Cys190    0.193080   HIS93 to CYS190\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Leu272  Leu272    0.442840\n",
      "Tyr219  Tyr219    0.410427\n",
      "Tyr326  Tyr326    0.325170\n",
      "Asp79    Asp79    0.321762\n",
      "Leu275  Leu275    0.305118\n",
      "Ile127  Ile127    0.305118\n",
      "Glu268  Glu268    0.272796\n",
      "Tyr141  Tyr141    0.272796\n",
      "Asn318  Asn318    0.257021\n",
      "Ala271  Ala271   -0.205418\n",
      "Cys190  Cys190    0.193080\n",
      "His93    His93    0.193080\n",
      "Using dark_background\n"
     ]
    }
   ],
   "source": [
    "tic_components_dir = tica_dir\n",
    "important_contact_features = interpret_tIC_components(projection_operator_dir, tic_components_dir, feature_residues_pkl, n_tica_components=5, percentile=95)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "feature_names = compat_verboseload(feature_residues_pkl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-04.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-06.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-10.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-16.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-12.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-08.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-14.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-18.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-22.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-20.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-13.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-17.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-09.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-07.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-05.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-24.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-19.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-21.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-23.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-11.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-05.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-07.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-04.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-06.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-09.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-15.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/C-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/C-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/B-08.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-04.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/C-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-06.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-05.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-08.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-10.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/A-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/C-04.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-05.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-07.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-04.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/F-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/D-09.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/C-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/G-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-06.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/H-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/H-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/E-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/H-04.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/F-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/I-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/I-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/J-00.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/J-02.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/G-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/J-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/H-01.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/H-05.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/H-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/I-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/J-03.dataset\n",
      "loading /home/enf/b2ar_analysis/featuresall_residues_2rh1_3sn6_under_cutoff6A/I-01.dataset\n"
     ]
    }
   ],
   "source": [
    "tic_subsampled_features_file = \"%s/features_subsampled.pkl\" % tica_dir\n",
    "subsampled_features_dir = os.path.join(tica_dir, \"subsampled_features\")\n",
    "if not os.path.exists(subsampled_features_dir): os.makedirs(subsampled_features_dir)\n",
    "important_contact_features_pruned, important_contact_features_indices = find_non_zero_features(important_contact_features[0], feature_names)\n",
    "subsample_features(features_dir, important_contact_features_indices, important_contact_features_pruned, tic_subsampled_features_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tica_coords = compat_verboseload(projected_features_dir)\n",
    "pnas_coords = compat_verboseload(pnas_coords_dir)\n",
    "for pnas_coord in pnas_coords: pnas_coord[:,0]*=7.14\n",
    "tica_names = [\"tIC.%d\" %i for i in range(1,n_components+1)]\n",
    "pnas_names = [\"tm6_tm3_dist\", \"rmsd_npxxy_inactive\", \"rmsd_npxxy_active\", \"rmsd_connector_inactive\", \"rmsd_connector_active\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import plots\n",
    "reload(plots)\n",
    "from plots import *\n",
    "#plot_histograms(projected_features_dir, analysis_dir, \"tICA histogram\", titles=[\"tIC.%d\" %i for i in range(1,n_components+1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lag_time = 25\n",
    "msm_model_dir = \"%s/msm_lag_time%d.h5\" % (tica_dir, lag_time)\n",
    "#build_msm(clusterer_tICs_1_2_3_filename, lag_time=lag_time, msm_model_dir=msm_model_dir)\n",
    "msm_object = compat_verboseload(msm_model_dir)\n",
    "prior_counts = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#compute_aggregate_scores(docking_multiple_ligands, inverse_agonists = inverse_ligands, summary = aggregate_docking, z_scores_csv = docking_z_scores_csv)\n",
    "#aggregate_docking_joined_map = convert_csv_to_joined_map(aggregate_docking, aggregate_docking_joined)[0]\n",
    "#aggregate_docking_means = calc_mean(aggregate_docking_joined_map)\n",
    "#write_map_to_csv(aggregate_docking_joined, aggregate_docking_means, [\"cluster\", \"mean_aggregate_docking_z_score\"])\n",
    "#r['do.analysis'](tica_dir, analysis_dir, pnas_coords_csv, tica_coords_csv, features_dir, docking_multiple_ligands)\n",
    "#tics_vs_docking_file = \"%s/tICA_vs_docking_carazolol.pdf\" % analysis_dir\n",
    "#plot_tICs_vs_docking(docking_multiple_ligands, tica_coords_csv, tics_vs_docking_file, chosen_ligand=\"s-carazolol\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "projection_operator_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Interpreting tIC 1\n",
      "feature_importances_df.shape\n",
      "(5, 5)\n",
      "residue_importances_df.shape\n",
      "(10, 2)\n",
      "       feature_name   res_i   res_j  importance           feature\n",
      "0   Ile72 To Thr123   Ile72  Thr123   -1.180972   ILE72 to THR123\n",
      "1   Leu75 To Pro323   Leu75  Pro323    1.069516   LEU75 to PRO323\n",
      "3  Phe332 To Ala335  Phe332  Ala335    0.763328  PHE332 to ALA335\n",
      "4  Ile135 To Gln224  Ile135  Gln224   -0.662528  ILE135 to GLN224\n",
      "2  Phe104 To Glu188  Phe104  Glu188    0.654189  PHE104 to GLU188\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Thr123  Thr123   -1.180972\n",
      "Ile72    Ile72   -1.180972\n",
      "Pro323  Pro323    1.069516\n",
      "Leu75    Leu75    1.069516\n",
      "Ala335  Ala335    0.763328\n",
      "Phe332  Phe332    0.763328\n",
      "Gln224  Gln224   -0.662528\n",
      "Ile135  Ile135   -0.662528\n",
      "Glu188  Glu188    0.654189\n",
      "Phe104  Phe104    0.654189\n",
      "Using dark_background\n",
      "Interpreting tIC 2\n",
      "feature_importances_df.shape\n",
      "(10, 5)\n",
      "residue_importances_df.shape\n",
      "(17, 2)\n",
      "       feature_name   res_i   res_j  importance           feature\n",
      "7  Ala128 To Val216  Ala128  Val216   -1.457099  ALA128 to VAL216\n",
      "0   Leu75 To Pro323   Leu75  Pro323    1.288608   LEU75 to PRO323\n",
      "8  Ala128 To Met279  Ala128  Met279    1.076218  ALA128 to MET279\n",
      "4  Phe332 To Ala335  Phe332  Ala335    1.027000  PHE332 to ALA335\n",
      "1    Phe89 To Trp99   Phe89   Trp99    0.888000    PHE89 to TRP99\n",
      "9  Val129 To Val222  Val129  Val222    0.848733  VAL129 to VAL222\n",
      "5    Phe61 To Thr66   Phe61   Thr66    0.705307    PHE61 to THR66\n",
      "2  Val213 To Val216  Val213  Val216   -0.637757  VAL213 to VAL216\n",
      "3  Phe332 To Ile334  Phe332  Ile334    0.581121  PHE332 to ILE334\n",
      "6   His93 To Lys305   His93  Lys305    0.417727   HIS93 to LYS305\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Leu75    Leu75    1.288608\n",
      "Pro323  Pro323    1.288608\n",
      "Met279  Met279    1.076218\n",
      "Ala335  Ala335    1.027000\n",
      "Phe332  Phe332    1.004706\n",
      "Ala128  Ala128    0.949553\n",
      "Trp99    Trp99    0.888000\n",
      "Phe89    Phe89    0.888000\n",
      "Val129  Val129    0.848733\n",
      "Val222  Val222    0.848733\n",
      "Thr66    Thr66    0.705307\n",
      "Phe61    Phe61    0.705307\n",
      "Val216  Val216   -0.678724\n",
      "Val213  Val213   -0.637757\n",
      "Ile334  Ile334    0.581121\n",
      "Lys305  Lys305    0.417727\n",
      "His93    His93    0.417727\n",
      "Using dark_background\n",
      "Interpreting tIC 3\n",
      "feature_importances_df.shape\n",
      "(7, 5)\n",
      "residue_importances_df.shape\n",
      "(14, 2)\n",
      "       feature_name   res_i   res_j  importance           feature\n",
      "6  Gly280 To Thr283  Gly280  Thr283    2.363325  GLY280 to THR283\n",
      "2  Val114 To Thr164  Val114  Thr164    1.655296  VAL114 to THR164\n",
      "5  Ser203 To Asn293  Ser203  Asn293    1.251909  SER203 to ASN293\n",
      "4  Ala198 To Ile201  Ala198  Ile201    0.924236  ALA198 to ILE201\n",
      "3  Tyr185 To Thr195  Tyr185  Thr195    0.679447  TYR185 to THR195\n",
      "0    Phe89 To Trp99   Phe89   Trp99   -0.438659    PHE89 to TRP99\n",
      "1  Phe104 To Glu188  Phe104  Glu188    0.405791  PHE104 to GLU188\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Thr283  Thr283    2.363325\n",
      "Gly280  Gly280    2.363325\n",
      "Thr164  Thr164    1.655296\n",
      "Val114  Val114    1.655296\n",
      "Asn293  Asn293    1.251909\n",
      "Ser203  Ser203    1.251909\n",
      "Ile201  Ile201    0.924236\n",
      "Ala198  Ala198    0.924236\n",
      "Thr195  Thr195    0.679447\n",
      "Tyr185  Tyr185    0.679447\n",
      "Trp99    Trp99   -0.438659\n",
      "Phe89    Phe89   -0.438659\n",
      "Glu188  Glu188    0.405791\n",
      "Phe104  Phe104    0.405791\n",
      "Using dark_background\n",
      "Interpreting tIC 4\n",
      "feature_importances_df.shape\n",
      "(12, 5)\n",
      "residue_importances_df.shape\n",
      "(24, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "2     Phe89 To Trp99   Phe89   Trp99    2.068668    PHE89 to TRP99\n",
      "6   Glu107 To Thr110  Glu107  Thr110    1.394653  GLU107 to THR110\n",
      "0    Leu75 To Pro323   Leu75  Pro323    1.000663   LEU75 to PRO323\n",
      "8   Leu167 To Tyr199  Leu167  Tyr199   -0.951380  LEU167 to TYR199\n",
      "1    Val87 To Tyr316   Val87  Tyr316   -0.933353   VAL87 to TYR316\n",
      "7   Met156 To Ser161  Met156  Ser161    0.914508  MET156 to SER161\n",
      "5   Phe104 To Glu188  Phe104  Glu188   -0.818826  PHE104 to GLU188\n",
      "9   Ala198 To Ile201  Ala198  Ile201    0.766072  ALA198 to ILE201\n",
      "11  Asp130 To Gln142  Asp130  Gln142   -0.755161  ASP130 to GLN142\n",
      "3    Ala91 To Trp109   Ala91  Trp109   -0.731690   ALA91 to TRP109\n",
      "4   Asn103 To Tyr185  Asn103  Tyr185   -0.643264  ASN103 to TYR185\n",
      "10   His93 To Lys305   His93  Lys305   -0.425386   HIS93 to LYS305\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Phe89    Phe89    2.068668\n",
      "Trp99    Trp99    2.068668\n",
      "Thr110  Thr110    1.394653\n",
      "Glu107  Glu107    1.394653\n",
      "Leu75    Leu75    1.000663\n",
      "Pro323  Pro323    1.000663\n",
      "Tyr199  Tyr199   -0.951380\n",
      "Leu167  Leu167   -0.951380\n",
      "Val87    Val87   -0.933353\n",
      "Tyr316  Tyr316   -0.933353\n",
      "Ser161  Ser161    0.914508\n",
      "Met156  Met156    0.914508\n",
      "Phe104  Phe104   -0.818826\n",
      "Glu188  Glu188   -0.818826\n",
      "Ala198  Ala198    0.766072\n",
      "Ile201  Ile201    0.766072\n",
      "Asp130  Asp130   -0.755161\n",
      "Gln142  Gln142   -0.755161\n",
      "Trp109  Trp109   -0.731690\n",
      "Ala91    Ala91   -0.731690\n",
      "Tyr185  Tyr185   -0.643264\n",
      "Asn103  Asn103   -0.643264\n",
      "His93    His93   -0.425386\n",
      "Lys305  Lys305   -0.425386\n",
      "Using dark_background\n",
      "Interpreting tIC 5\n",
      "feature_importances_df.shape\n",
      "(17, 5)\n",
      "residue_importances_df.shape\n",
      "(30, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "9   Tyr174 To Thr195  Tyr174  Thr195    1.866706  TYR174 to THR195\n",
      "10  Val213 To Val216  Val213  Val216   -1.358702  VAL213 to VAL216\n",
      "0     Phe89 To Trp99   Phe89   Trp99   -1.345337    PHE89 to TRP99\n",
      "5   Met156 To Ser161  Met156  Ser161    1.216413  MET156 to SER161\n",
      "8   Ile169 To Tyr199  Ile169  Tyr199    1.134999  ILE169 to TYR199\n",
      "7   Ile169 To Tyr174  Ile169  Tyr174    1.069908  ILE169 to TYR174\n",
      "6   Leu167 To Gln170  Leu167  Gln170    1.067885  LEU167 to GLN170\n",
      "11  Gly280 To Thr283  Gly280  Thr283   -1.030235  GLY280 to THR283\n",
      "2   Glu107 To Thr110  Glu107  Thr110    0.959543  GLU107 to THR110\n",
      "3   Trp109 To Tyr316  Trp109  Tyr316   -0.886761  TRP109 to TYR316\n",
      "14   Val48 To Ser319   Val48  Ser319   -0.851106   VAL48 to SER319\n",
      "13  Phe332 To Ile334  Phe332  Ile334   -0.765788  PHE332 to ILE334\n",
      "16  Ile135 To His269  Ile135  His269    0.764262  ILE135 to HIS269\n",
      "1    Phe89 To Trp109   Phe89  Trp109   -0.640965   PHE89 to TRP109\n",
      "4   Ser111 To Val114  Ser111  Val114    0.628804  SER111 to VAL114\n",
      "12  Phe289 To Leu311  Phe289  Leu311   -0.510120  PHE289 to LEU311\n",
      "15  Ala134 To Ser143  Ala134  Ser143    0.362678  ALA134 to SER143\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Thr195  Thr195    1.866706\n",
      "Tyr174  Tyr174    1.826866\n",
      "Val213  Val213   -1.358702\n",
      "Val216  Val216   -1.358702\n",
      "Trp99    Trp99   -1.345337\n",
      "Ser161  Ser161    1.216413\n",
      "Met156  Met156    1.216413\n",
      "Tyr199  Tyr199    1.134999\n",
      "Ile169  Ile169    1.131745\n",
      "Gln170  Gln170    1.067885\n",
      "Leu167  Leu167    1.067885\n",
      "Gly280  Gly280   -1.030235\n",
      "Thr283  Thr283   -1.030235\n",
      "Thr110  Thr110    0.959543\n",
      "Glu107  Glu107    0.959543\n",
      "Tyr316  Tyr316   -0.886761\n",
      "Val48    Val48   -0.851106\n",
      "Ser319  Ser319   -0.851106\n",
      "Phe332  Phe332   -0.765788\n",
      "Ile334  Ile334   -0.765788\n",
      "His269  His269    0.764262\n",
      "Ile135  Ile135    0.764262\n",
      "Phe89    Phe89   -0.676184\n",
      "Trp109  Trp109   -0.653255\n",
      "Val114  Val114    0.628804\n",
      "Ser111  Ser111    0.628804\n",
      "Phe289  Phe289   -0.510120\n",
      "Leu311  Leu311   -0.510120\n",
      "Ala134  Ala134    0.362678\n",
      "Ser143  Ser143    0.362678\n",
      "Using dark_background\n",
      "Interpreting tIC 6\n",
      "feature_importances_df.shape\n",
      "(13, 5)\n",
      "residue_importances_df.shape\n",
      "(24, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "3     Phe89 To Trp99   Phe89   Trp99    2.171456    PHE89 to TRP99\n",
      "10  Tyr326 To Cys327  Tyr326  Cys327   -1.294144  TYR326 to CYS327\n",
      "11  His178 To Asp300  His178  Asp300   -1.268916  HIS178 to ASP300\n",
      "0    Asn51 To Ser319   Asn51  Ser319    1.086656   ASN51 to SER319\n",
      "4    Ala91 To Trp109   Ala91  Trp109   -1.014983   ALA91 to TRP109\n",
      "2     Ala76 To Ala78   Ala76   Ala78    0.962534    ALA76 to ALA78\n",
      "7   Phe104 To Glu188  Phe104  Glu188    0.836899  PHE104 to GLU188\n",
      "5   Asn103 To Ala186  Asn103  Ala186    0.782822  ASN103 to ALA186\n",
      "9   Phe289 To Leu311  Phe289  Leu311    0.753307  PHE289 to LEU311\n",
      "12  Ala181 To Thr195  Ala181  Thr195   -0.702196  ALA181 to THR195\n",
      "6   Asn103 To Glu188  Asn103  Glu188   -0.659318  ASN103 to GLU188\n",
      "1    Ile72 To Thr123   Ile72  Thr123    0.584199   ILE72 to THR123\n",
      "8   Val213 To Val218  Val213  Val218    0.518289  VAL213 to VAL218\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Phe89    Phe89    2.171456\n",
      "Trp99    Trp99    2.171456\n",
      "Cys327  Cys327   -1.294144\n",
      "Tyr326  Tyr326   -1.294144\n",
      "Asp300  Asp300   -1.268916\n",
      "His178  His178   -1.268916\n",
      "Asn51    Asn51    1.086656\n",
      "Ser319  Ser319    1.086656\n",
      "Trp109  Trp109   -1.014983\n",
      "Ala91    Ala91   -1.014983\n",
      "Ala76    Ala76    0.962534\n",
      "Ala78    Ala78    0.962534\n",
      "Phe104  Phe104    0.836899\n",
      "Ala186  Ala186    0.782822\n",
      "Glu188  Glu188    0.762089\n",
      "Phe289  Phe289    0.753307\n",
      "Leu311  Leu311    0.753307\n",
      "Asn103  Asn103    0.710715\n",
      "Ala181  Ala181   -0.702196\n",
      "Thr195  Thr195   -0.702196\n",
      "Thr123  Thr123    0.584199\n",
      "Ile72    Ile72    0.584199\n",
      "Val213  Val213    0.518289\n",
      "Val218  Val218    0.518289\n",
      "Using dark_background\n",
      "Interpreting tIC 7\n",
      "feature_importances_df.shape\n",
      "(19, 5)\n",
      "residue_importances_df.shape\n",
      "(29, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "13  Arg175 To Ala176  Arg175  Ala176    4.251978  ARG175 to ALA176\n",
      "14  Arg175 To Thr177  Arg175  Thr177   -3.520309  ARG175 to THR177\n",
      "11  Trp173 To Arg175  Trp173  Arg175    2.597660  TRP173 to ARG175\n",
      "12  Tyr174 To Thr195  Tyr174  Thr195    2.512260  TYR174 to THR195\n",
      "2     Ala76 To Ala78   Ala76   Ala78    2.020504    ALA76 to ALA78\n",
      "7   Ile169 To Tyr199  Ile169  Tyr199    1.995631  ILE169 to TYR199\n",
      "8   Gln170 To Tyr174  Gln170  Tyr174    1.943937  GLN170 to TYR174\n",
      "10  His172 To Arg175  His172  Arg175   -1.898820  HIS172 to ARG175\n",
      "16  Asp331 To Ile334  Asp331  Ile334   -1.624373  ASP331 to ILE334\n",
      "17  Phe332 To Ala335  Phe332  Ala335    1.566055  PHE332 to ALA335\n",
      "4   Glu107 To Thr110  Glu107  Thr110   -1.464964  GLU107 to THR110\n",
      "0     Leu75 To Asp79   Leu75   Asp79   -1.168642    LEU75 to ASP79\n",
      "5   Met156 To Ser161  Met156  Ser161   -1.140607  MET156 to SER161\n",
      "6   Ile169 To Arg175  Ile169  Arg175    1.130920  ILE169 to ARG175\n",
      "15  Ser204 To Val206  Ser204  Val206    1.090576  SER204 to VAL206\n",
      "9   Gln170 To Arg175  Gln170  Arg175    0.953853  GLN170 to ARG175\n",
      "1    Leu75 To Trp158   Leu75  Trp158    0.793085   LEU75 to TRP158\n",
      "18   Met98 To Phe101   Met98  Phe101   -0.452538   MET98 to PHE101\n",
      "3    Trp99 To Cys191   Trp99  Cys191   -0.224606   TRP99 to CYS191\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Ala176  Ala176    4.251978\n",
      "Arg175  Arg175    3.838398\n",
      "Thr177  Thr177   -3.520309\n",
      "Trp173  Trp173    2.597660\n",
      "Thr195  Thr195    2.512260\n",
      "Tyr174  Tyr174    2.483844\n",
      "Ala76    Ala76    2.020504\n",
      "Ala78    Ala78    2.020504\n",
      "Tyr199  Tyr199    1.995631\n",
      "Ile169  Ile169    1.952395\n",
      "His172  His172   -1.898820\n",
      "Gln170  Gln170    1.894433\n",
      "Asp331  Asp331   -1.624373\n",
      "Ile334  Ile334   -1.624373\n",
      "Phe332  Phe332    1.566055\n",
      "Ala335  Ala335    1.566055\n",
      "Thr110  Thr110   -1.464964\n",
      "Glu107  Glu107   -1.464964\n",
      "Asp79    Asp79   -1.168642\n",
      "Ser161  Ser161   -1.140607\n",
      "Met156  Met156   -1.140607\n",
      "Ser204  Ser204    1.090576\n",
      "Val206  Val206    1.090576\n",
      "Trp158  Trp158    0.793085\n",
      "Leu75    Leu75    0.694999\n",
      "Met98    Met98   -0.452538\n",
      "Phe101  Phe101   -0.452538\n",
      "Trp99    Trp99   -0.224606\n",
      "Cys191  Cys191   -0.224606\n",
      "Using dark_background\n",
      "Interpreting tIC 8\n",
      "feature_importances_df.shape\n",
      "(13, 5)\n",
      "residue_importances_df.shape\n",
      "(22, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "5   Phe133 To Arg221  Phe133  Arg221    6.547490  PHE133 to ARG221\n",
      "9   Lys270 To Thr274  Lys270  Thr274   -5.607832  LYS270 to THR274\n",
      "2   Tyr132 To Phe217  Tyr132  Phe217    4.469084  TYR132 to PHE217\n",
      "12   Val87 To Ser319   Val87  Ser319   -4.365522   VAL87 to SER319\n",
      "4   Phe133 To Ala134  Phe133  Ala134    4.130300  PHE133 to ALA134\n",
      "11   Ala85 To Cys116   Ala85  Cys116   -3.809418   ALA85 to CYS116\n",
      "3   Tyr132 To Glu225  Tyr132  Glu225    3.188799  TYR132 to GLU225\n",
      "8   Lys267 To Lys270  Lys267  Lys270    2.111741  LYS267 to LYS270\n",
      "0   Cys125 To Leu212  Cys125  Leu212    2.007287  CYS125 to LEU212\n",
      "7   Tyr174 To Thr195  Tyr174  Thr195   -0.957343  TYR174 to THR195\n",
      "6   Met156 To Ser161  Met156  Ser161   -0.650240  MET156 to SER161\n",
      "10  Met279 To Thr281  Met279  Thr281    0.626873  MET279 to THR281\n",
      "1   Ile127 To Tyr132  Ile127  Tyr132    0.544145  ILE127 to TYR132\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Arg221  Arg221    6.547490\n",
      "Phe133  Phe133    6.426631\n",
      "Thr274  Thr274   -5.607832\n",
      "Phe217  Phe217    4.469084\n",
      "Ser319  Ser319   -4.365522\n",
      "Val87    Val87   -4.365522\n",
      "Tyr132  Tyr132    4.341056\n",
      "Ala134  Ala134    4.130300\n",
      "Cys116  Cys116   -3.809418\n",
      "Ala85    Ala85   -3.809418\n",
      "Glu225  Glu225    3.188799\n",
      "Lys267  Lys267    2.111741\n",
      "Leu212  Leu212    2.007287\n",
      "Cys125  Cys125    2.007287\n",
      "Lys270  Lys270    1.725762\n",
      "Tyr174  Tyr174   -0.957343\n",
      "Thr195  Thr195   -0.957343\n",
      "Met156  Met156   -0.650240\n",
      "Ser161  Ser161   -0.650240\n",
      "Met279  Met279    0.626873\n",
      "Thr281  Thr281    0.626873\n",
      "Ile127  Ile127    0.544145\n",
      "Using dark_background\n",
      "Interpreting tIC 9\n",
      "feature_importances_df.shape\n",
      "(18, 5)\n",
      "residue_importances_df.shape\n",
      "(28, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "5   Ser137 To Glu225  Ser137  Glu225    3.516679  SER137 to GLU225\n",
      "7   Pro138 To Gln142  Pro138  Gln142    3.173817  PRO138 to GLN142\n",
      "8   Leu163 To Leu167  Leu163  Leu167    2.470449  LEU163 to LEU167\n",
      "1     Thr66 To Thr73   Thr66   Thr73    2.333143    THR66 to THR73\n",
      "6   Pro138 To Phe139  Pro138  Phe139    2.267625  PRO138 to PHE139\n",
      "3   Phe133 To Arg221  Phe133  Arg221   -2.131505  PHE133 to ARG221\n",
      "9   Leu163 To Pro168  Leu163  Pro168   -2.031260  LEU163 to PRO168\n",
      "4   Thr136 To Arg221  Thr136  Arg221   -1.848268  THR136 to ARG221\n",
      "13  Pro168 To His172  Pro168  His172    1.558619  PRO168 to HIS172\n",
      "10  Thr164 To Leu167  Thr164  Leu167    1.556199  THR164 to LEU167\n",
      "12  Pro168 To Met171  Pro168  Met171   -1.425854  PRO168 to MET171\n",
      "14  Pro168 To Arg175  Pro168  Arg175   -0.955627  PRO168 to ARG175\n",
      "15  Tyr174 To Thr195  Tyr174  Thr195   -0.947528  TYR174 to THR195\n",
      "11  Pro168 To Ile169  Pro168  Ile169   -0.929003  PRO168 to ILE169\n",
      "2   Cys106 To Phe193  Cys106  Phe193    0.797698  CYS106 to PHE193\n",
      "16  His269 To Lys270  His269  Lys270    0.631766  HIS269 to LYS270\n",
      "0    Leu64 To Ala335   Leu64  Ala335   -0.615229   LEU64 to ALA335\n",
      "17  His178 To Asp300  His178  Asp300   -0.538523  HIS178 to ASP300\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Ser137  Ser137    3.516679\n",
      "Glu225  Glu225    3.516679\n",
      "Gln142  Gln142    3.173817\n",
      "Pro138  Pro138    3.128508\n",
      "Leu167  Leu167    2.424737\n",
      "Thr73    Thr73    2.333143\n",
      "Thr66    Thr66    2.333143\n",
      "Phe139  Phe139    2.267625\n",
      "Leu163  Leu163    2.245364\n",
      "Phe133  Phe133   -2.131505\n",
      "Arg221  Arg221   -1.862430\n",
      "Thr136  Thr136   -1.848268\n",
      "His172  His172    1.558619\n",
      "Thr164  Thr164    1.556199\n",
      "Met171  Met171   -1.425854\n",
      "Pro168  Pro168    1.061095\n",
      "Arg175  Arg175   -0.955627\n",
      "Thr195  Thr195   -0.947528\n",
      "Tyr174  Tyr174   -0.947528\n",
      "Ile169  Ile169   -0.929003\n",
      "Phe193  Phe193    0.797698\n",
      "Cys106  Cys106    0.797698\n",
      "Lys270  Lys270    0.631766\n",
      "His269  His269    0.631766\n",
      "Ala335  Ala335   -0.615229\n",
      "Leu64    Leu64   -0.615229\n",
      "His178  His178   -0.538523\n",
      "Asp300  Asp300   -0.538523\n",
      "Using dark_background\n",
      "Interpreting tIC 10\n",
      "feature_importances_df.shape\n",
      "(25, 5)\n",
      "residue_importances_df.shape\n",
      "(37, 2)\n",
      "        feature_name   res_i   res_j  importance           feature\n",
      "13  Lys270 To Thr274  Lys270  Thr274    3.794095  LYS270 to THR274\n",
      "7   Phe133 To Arg221  Phe133  Arg221    2.389175  PHE133 to ARG221\n",
      "12  Glu268 To Ala271  Glu268  Ala271   -2.087100  GLU268 to ALA271\n",
      "24  Tyr141 To Ala226  Tyr141  Ala226    1.823428  TYR141 to ALA226\n",
      "18   Val87 To Ser319   Val87  Ser319   -1.777070   VAL87 to SER319\n",
      "10  Lys267 To Lys270  Lys267  Lys270   -1.745172  LYS267 to LYS270\n",
      "6   Phe133 To Ile135  Phe133  Ile135    1.502533  PHE133 to ILE135\n",
      "14  Ala271 To Leu272  Ala271  Leu272   -1.435527  ALA271 to LEU272\n",
      "11  Glu268 To Lys270  Glu268  Lys270   -1.316713  GLU268 to LYS270\n",
      "9   Lys267 To His269  Lys267  His269   -1.205864  LYS267 to HIS269\n",
      "19  Val129 To Met279  Val129  Met279   -1.193216  VAL129 to MET279\n",
      "4   Cys125 To Leu212  Cys125  Leu212    1.000602  CYS125 to LEU212\n",
      "17   Thr68 To Ser143   Thr68  Ser143   -0.966338   THR68 to SER143\n",
      "23  Lys140 To Glu225  Lys140  Glu225   -0.945357  LYS140 to GLU225\n",
      "5   Tyr132 To Phe217  Tyr132  Phe217    0.927413  TYR132 to PHE217\n",
      "8   Val213 To Phe217  Val213  Phe217    0.895219  VAL213 to PHE217\n",
      "15  Leu272 To Ile277  Leu272  Ile277   -0.877006  LEU272 to ILE277\n",
      "1     Phe89 To Trp99   Phe89   Trp99    0.827264    PHE89 to TRP99\n",
      "16  Thr274 To Ile278  Thr274  Ile278    0.783670  THR274 to ILE278\n",
      "20  Arg131 To Met279  Arg131  Met279    0.739430  ARG131 to MET279\n",
      "22  Thr136 To Leu144  Thr136  Leu144   -0.719756  THR136 to LEU144\n",
      "3   Ile121 To Leu212  Ile121  Leu212   -0.490657  ILE121 to LEU212\n",
      "21  Ile135 To Phe139  Ile135  Phe139    0.456708  ILE135 to PHE139\n",
      "0     Thr68 To Ile72   Thr68   Ile72    0.441253    THR68 to ILE72\n",
      "2   Asn103 To Ala186  Asn103  Ala186    0.262933  ASN103 to ALA186\n",
      "Using dark_background\n",
      "       residue  importance\n",
      "Thr274  Thr274    3.643574\n",
      "Lys270  Lys270    3.283014\n",
      "Arg221  Arg221    2.389175\n",
      "Phe133  Phe133    2.344842\n",
      "Ala226  Ala226    1.823428\n",
      "Tyr141  Tyr141    1.823428\n",
      "Ser319  Ser319   -1.777070\n",
      "Val87    Val87   -1.777070\n",
      "Ala271  Ala271   -1.468105\n",
      "Ile135  Ile135    1.450242\n",
      "Glu268  Glu268   -1.355232\n",
      "Lys267  Lys267   -1.232829\n",
      "His269  His269   -1.205864\n",
      "Val129  Val129   -1.193216\n",
      "Cys125  Cys125    1.000602\n",
      "Ser143  Ser143   -0.966338\n",
      "Glu225  Glu225   -0.945357\n",
      "Lys140  Lys140   -0.945357\n",
      "Tyr132  Tyr132    0.927413\n",
      "Leu212  Leu212    0.926039\n",
      "Phe217  Phe217    0.925804\n",
      "Leu272  Leu272   -0.904932\n",
      "Val213  Val213    0.895219\n",
      "Ile277  Ile277   -0.877006\n",
      "Trp99    Trp99    0.827264\n",
      "Phe89    Phe89    0.827264\n",
      "Ile278  Ile278    0.783670\n",
      "Arg131  Arg131    0.739430\n",
      "Leu144  Leu144   -0.719756\n",
      "Thr136  Thr136   -0.719756\n",
      "Met279  Met279    0.642798\n",
      "Ile121  Ile121   -0.490657\n",
      "Phe139  Phe139    0.456708\n",
      "Ile72    Ile72    0.441253\n",
      "Thr68    Thr68    0.370873\n",
      "Asn103  Asn103    0.262933\n",
      "Ala186  Ala186    0.262933\n",
      "Using dark_background\n"
     ]
    }
   ],
   "source": [
    "alt_pp_tica_dir = \"/home/enf/b2ar_analysis/sparse-tICA_t5_n_components25all_residues_2rh1_3sn6_under_cutoff6A_regularization_wolf_autoShrinkage_rho0pt01\"\n",
    "alt_projection_operator_dir = \"%s/phi_psi_chi2_allprot_tica_coords.h5\" %alt_pp_tica_dir\n",
    "tic_subsampled_features_file = \"%s/features_subsampled.pkl\" % alt_pp_tica_dir\n",
    "\n",
    "alt_important_contact_features = interpret_tIC_components(alt_projection_operator_dir, alt_pp_tica_dir, feature_residues_pkl, n_tica_components=10, percentile=95)\n",
    "\n",
    "subsampled_features_dir = os.path.join(alt_pp_tica_dir, \"subsampled_features\")\n",
    "if not os.path.exists(subsampled_features_dir): os.makedirs(subsampled_features_dir)\n",
    "alt_important_contact_features_pruned, alt_important_contact_features_indices = find_non_zero_features(alt_important_contact_features[0], feature_names)\n",
    "if not os.path.exists(tic_subsampled_features_file):\n",
    "    subsample_features(features_dir, alt_important_contact_features_indices, alt_important_contact_features_pruned, tic_subsampled_features_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading /home/enf/b2ar_analysis/sparse-tICA_t5_n_components25all_residues_2rh1_3sn6_under_cutoff6A_regularization_wolf_autoShrinkage_rho0pt01/features_subsampled.pkl\n"
     ]
    }
   ],
   "source": [
    "top_features = load_file(tic_subsampled_features_file)\n",
    "top_features = [t*10. for t in top_features]\n",
    "\n",
    "user_defined_coords = compat_verboseload(user_defined_features_file)\n",
    "user_defined_names = sorted(feature_name_residues_dict.keys())\n",
    "user_defined_dfs = [pd.DataFrame(t, columns=user_defined_names) for t in user_defined_coords]\n",
    "\n",
    "tica_dfs = [pd.DataFrame(t, columns=[\"tIC.%d\" %i for i in range(1,n_components+1)]) for t in tica_coords]\n",
    "\n",
    "all_feature_dfs = [pd.concat([top_features[i], user_defined_dfs[i], tica_dfs[i]], axis=1) for i in range(0, len(top_features))]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'clusterer' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-30-5a0f06aedd18>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mall_traj_features_np\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvalues\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mf\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mall_feature_dfs\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mcluster_features_averages\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcalculate_cluster_averages_per_feature\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mclusterer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mall_traj_features_np\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mcluster_features_averages\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcluster_features_averages\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcolumns\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mall_feature_dfs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"cluster%d\"\u001b[0m \u001b[1;33m%\u001b[0m\u001b[0mi\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mn_clusters\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'clusterer' is not defined"
     ]
    }
   ],
   "source": [
    "all_traj_features_np = [f.values for f in all_feature_dfs]\n",
    "cluster_features_averages = calculate_cluster_averages_per_feature(clusterer, all_traj_features_np)\n",
    "cluster_features_averages = pd.DataFrame(cluster_features_averages, columns=all_feature_dfs[0].columns, index=[\"cluster%d\" %i for i in range(0,n_clusters)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import custom_clusterer\n",
    "reload(custom_clusterer)\n",
    "from custom_clusterer import *\n",
    "feature_name = \"Asn148-Leu266_ca_dist\"\n",
    "find_snapshots_within_feature_range(all_feature_dfs, feature_name, [31., 32.], \n",
    "                                    get_trajectory_files(traj_dir, traj_ext), analysis_dir,\n",
    "                                    \"%s_31_32\" %feature_name, 5, lig_name=\"BIA\", \n",
    "                                    structure=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cluster_features_averages.loc[(cluster_features_averages[\"Asn148-Leu266_ca_dist\"] > 35.) & (cluster_features_averages[\"Asn148-Leu266_ca_dist\"] < 40.)][\"Asn148-Leu266_ca_dist\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "all_features_onehot, names_onehot = multi_onehot_trajectories([t.values for t in all_feature_dfs], all_feature_dfs[0].columns.values.tolist(), subsample=100)\n",
    "all_features_onehot = [pd.DataFrame(t, columns=names_onehot) for t in all_features_onehot]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[n for n in top_features[0].columns.values.tolist() if \"TYR326\" in n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from imp import reload\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "\n",
    "clusterer, cluster_averages, active_clusters, inactive_clusters, biased_ligands, agonist_ligands, inverse_ligands, all_ligands, c, feature_residues, tica_coords, user_defined_coords, pp_n_components, apriori_dfs, tica_dfs, cluster_pnas_averages, cluster_tica_averages, cluster_tica_pnas, top_features, clusters_map, tica_resampled_file, projected_features, num_trajs, features_eq, all_traj_features, samples_indices_file, samples_dir, samples_tica_avg_df, samples_pnas_avg_df, samples_features_avg_df, samples_normalized_features_avg_df, feature_names, feature_strings, samples_pnas_tica, reference_docking = initialize_analysis(clusterer_dir, user_defined_coords, user_defined_names, biased_agonist_dir, agonist_dir, inverse_agonist_dir, docking_dir, precision, docking_multiple_ligands, aggregate_docking, feature_residues_pkl, n_components, all_feature_dfs,\n",
    "                        lag_time, n_clusters, projected_features_dir, traj_dir, traj_ext, tica_dir,\n",
    "                        prior_counts, msm_object, analysis_dir, n_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "analysis_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import msm_resampled\n",
    "reload(msm_resampled)\n",
    "from msm_resampled import *\n",
    "\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "\n",
    "docking_cluster_averages = c\n",
    "all_ligands_dir = \"/home/enf/b2ar_analysis/all_ligands\"\n",
    "ligands = get_ligands(all_ligands_dir)\n",
    "apo_populations, df_agg, aggregate_docking_msm, docking_normalized, ddg_scaled, deltas_tica, delta_delta_g, lig_features_eq, new_populations, bi_msm, num_trajs, features, null_features, classes, agonists, antagonists, labels, X, N, C, y = compute_docking_dg(docking_cluster_averages, msm_object, samples_tica_avg_df, samples_pnas_avg_df, samples_normalized_features_avg_df, important_contact_features, traj_dir, traj_ext, tica_dir, ligands, reference_docking, clusters_map, all_features_onehot, analysis_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lig_continuous_features = msm_reweighted_features_per_ligand(all_feature_dfs, new_populations, bi_msm, 10000, clusters_map, num_trajs, apo_populations, save_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import detect_intermediates\n",
    "reload(detect_intermediates)\n",
    "from detect_intermediates import *\n",
    "\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "#compare_feature_to_apo(lig_continuous_features, [\"r_isopreterenol\", \"nebivolol\", \"3p0g_lig\", \"r_epinephrine\", \"s-carvedilol\"], \"\")\n",
    "ligands = [\"r_isopreterenol\", \"3p0g_lig\", \"r_epinephrine\", 'Ici118551', \"s-carazolol\", \"salbutamol\", \"salmeterol\"]\n",
    "compare_feature_to_apo(lig_continuous_features, ligands, \"s-carazolol\", 'Asn148-Leu266_ca_dist')\n",
    "plot_overall_kde(lig_continuous_features, ligands, 'Asn148-Leu266_ca_dist')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "X_binarized = np.zeros((len(lig_features_eq.keys()), lig_features_eq[lig_features_eq.keys()[0]].shape[1]))\n",
    "for i, lig in enumerate(lig_features_eq.keys()):\n",
    "    print(i)\n",
    "    x = lig_features_eq[lig]\n",
    "    X_binarized[i,:] = x.astype(bool).sum(axis=0).values\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_onehot_df = pd.DataFrame(X_binarized, index=lig_features_eq.keys(), columns=names_onehot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "[n for n in names_onehot if \"TYR326\" in n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_onehot_df['rmsd_npxxy_active < 0.334651'].sort(inplace=False).plot(kind='barh')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#X_all_features = np.zeros((len(lig_features_eq.keys()), lig_features_eq[lig_features_eq.keys()[0]].shape[1]))\n",
    "#for i, lig in enumerate(lig_features_eq.keys()):\n",
    "#    x = lig_features_eq[lig].mean().values\n",
    "#    X_all_features[i,:] = x\n",
    "#X_all_features_df = pd.DataFrame(X_all_features, index=lig_features_eq.keys(), columns=lig_features_eq[lig].columns)\n",
    "#X_all_features_df = standardize_df(X_all_features_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "salt_bridge = []\n",
    "for ligand in common_ligands:\n",
    "    salt_bridge.append(lig_features_eq[ligand].loc[lig_features_eq[ligand][\"Glu268-Arg328_dist\"] < 5.].shape[0])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "salt_bridge = []\n",
    "for ligand in common_ligands:\n",
    "    salt_bridge.append(lig_features_eq[ligand].loc[lig_features_eq[ligand][\"Asn148-Leu266_dist\"] > 37.5].shape[0])\n",
    "pd.DataFrame(salt_bridge, index=common_ligands, columns=[\"dist\"]).sort(columns=\"dist\", inplace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.hist(lig_features_eq[\"r_isopreterenol\"][\"Asn148-Leu266_ca_dist\"].values, bins=50)\n",
    "plt.show()\n",
    "plt.hist(lig_features_eq[\"s-carazolol\"][\"Asn148-Leu266_ca_dist\"].values, bins=50)\n",
    "plt.show()\n",
    "plt.hist(lig_features_eq[\"practolol\"][\"Asn148-Leu266_ca_dist\"].values, bins=50)\n",
    "plt.show()\n",
    "plt.hist(lig_features_eq[\"practolol\"][\"tm6_tm3_dist\"].values, bins=50)\n",
    "plt.show()\n",
    "plt.hist(lig_features_eq[\"r_isopreterenol\"][\"tm6_tm3_dist\"].values, bins=50)\n",
    "plt.show()\n",
    "plt.hist(lig_features_eq[\"s-carvedilol\"][\"tm6_tm3_dist\"].values, bins=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(standardize_df(X_onehot_df.loc[common_ligands][X_onehot_df.columns.values[-13]]).values, bret.loc[common_ligands][\"B2AR-Gprotein, Mean\"])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(compute_pearson_matrix(standardize_df(X_onehot_df.loc[common_ligands].transpose()).values.T,bret[\"B2AR-Gprotein, Mean\"].subtract(bret[\"B2AR-Arrestin, Mean\"]).loc[common_agonists].values.reshape((-1,1))), index=X_onehot_df.columns, columns=[\"correlation\"]).sort(\"correlation\", inplace=False)\n",
    "print(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "list(feature_name_residues_dict.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pd.set_option('display.max_rows', None)\n",
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cluster_pnas_averages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples_features_avg_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bret = pd.read_excel(\"/home/enf/b2ar_analysis/bret_bias_study2.xlsx\", header=0).dropna().set_index(\"EvanName\")\n",
    "#bret = pd.read_csv(\"/home/enf/b2ar_analysis/bias_analysis/bret_bias_study.csv\", header=0).dropna().set_index(\"EvanName\")\n",
    "common_ligands = [n for n in bret.index.values if n in delta_delta_g.columns.values]\n",
    "bret[\"B2AR-Arrestin, Mean\"].sort(inplace=False)#.subtract(bret[\"B2AR-Arrestin, Mean\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "#arrestin_antagonists = [\"s-carvedilol\", \"nebivolol\"]\n",
    "#non_arrestin_antagonists = [n for n in antagonists if n not in arrestin_antagonists and n not in [\"Carvedilol\"]]\n",
    "#y = np.array([1. for i in arrestin_antagonists] + [0. for i in non_arrestin_antagonists]).reshape((-1,1))\n",
    "\n",
    "\n",
    "total_activity = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands])\n",
    "#common_agonists = arrestin_antagonists + non_arrestin_antagonists\n",
    "#biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\"]\n",
    "#non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "#y = np.array([1. for i in biased_ligands] + [0. for i in non_biased_ligands]).reshape((-1,1))\n",
    "#common_agonists = biased_ligands + non_biased_ligands\n",
    "\n",
    "common_agonists = total_activity.loc[total_activity > 0.2].index.values\n",
    "y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].subtract(bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists]).values.reshape((-1,1))\n",
    "\n",
    "top_clusters = delta_delta_g.index.values\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"nebivolol\").index.values[:10].tolist() + delta_delta_g.sort(\"3p0g_lig\").index.values[:10].tolist()))\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"N-Cyclopentylbutanephrine\", inplace=False).index.values[:4].tolist() + delta_delta_g.sort(\"procaterol\", inplace=False).index.values[:4].tolist()))\n",
    "#agonists_df = [a for a in agonists if a in delta_delta_g.columns.values]\n",
    "#common_agonists = agonists_df + antagonists\n",
    "#y = np.array([1. for i in agonists_df] + [0. for i in antagonists]).reshape((-1,1))\n",
    "#common_agonists = common_ligands\n",
    "y_arr = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "y_gpr = bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "#y_ori = y_arr - y_gpr\n",
    "#y_ori = y_arr\n",
    "#y = y_arr\n",
    "#y = multi_binarizer(y_gpr, [0.2, 0.8])\n",
    "y = multi_binarizer(y_ori, [-0.2])\n",
    "print(y)\n",
    "\n",
    "C = null_features.loc[common_agonists].values\n",
    "#X = delta_delta_g.loc[top_clusters][common_agonists].values.T\n",
    "X = np.hstack([delta_delta_g.loc[top_clusters][common_agonists].values.T, C*-1.0])\n",
    "#X_scaled = ddg_scaled.loc[top_clusters][common_agonists].values.T\n",
    "X_scaled = np.hstack([ddg_scaled.loc[top_clusters][common_agonists].values.T, C*-1.0])\n",
    "#D_scaled = docking_normalized.loc[top_clusters][common_agonists].values.T\n",
    "D_scaled = np.hstack([docking_normalized.loc[top_clusters][common_agonists].values.T, C*-1.0])\n",
    "X_diff = np.zeros((X.shape[0], X.shape[1]**2/2))\n",
    "k=0\n",
    "for i in range(0,X.shape[1]):\n",
    "    for j in range(i+1,X.shape[1]):\n",
    "        X_diff[:,k] = X[:,i] - X[:,j]\n",
    "        k+=1\n",
    "X_diff = np.hstack([X_scaled, X_diff, C*-1.0])\n",
    "\n",
    "all_features_df = pd.concat([ddg_scaled.loc[top_clusters], null_features.multiply(-1.0).transpose(), X_onehot_df.transpose()], axis=0)[common_agonists].transpose()\n",
    "all_features = all_features_df.values\n",
    "features = [C*-1.0, X_scaled]\n",
    "features_y = [C*-1.0, X_scaled]\n",
    "feature_names = [\"Crystal Structures\", \"Docking ddG Scaled\"]#, \"Docking ddG and Observables\"]\n",
    "\n",
    "n_trials=1000\n",
    "a_vs_g_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=True, model=\"logistic_cv\")\n",
    "with open(\"%s/a_vs_g_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(a_vs_g_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs, f)\n",
    "\n",
    "\n",
    "analyze_multiclass_experiment(a_vs_g_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs, \n",
    "                              [\"Crystal Structures\", \"MSM Docking\"],\n",
    "                              top_clusters.tolist() + null_features.columns.values.tolist(), common_agonists, analysis_dir,\n",
    "                              [\"Arrestin vs. G Protein\", \"G Protein Agonists\"], X_scaled, \n",
    "                              exp_title=\"Arrestin vs. G Protein Two Class\", coef_name=\"Logistic Coefficient\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_features_df.loc[\"r_isopreterenol\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "analysis_dir = \"%s_%dsamples\" %(analysis_dir, n_samples)\n",
    "print(analysis_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(analysis_dir):\n",
    "    os.makedirs(analysis_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_trials = 1000\n",
    "\n",
    "y = multi_binarizer(y_gpr, [0.33, 0.66])\n",
    "gprot_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/gprot_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f, f)\n",
    "\n",
    "y = multi_binarizer(y_arr, [0.33, 0.66])\n",
    "arr_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/arr_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f, f)\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "    \n",
    "    \n",
    "y = multi_binarizer(y_gpr, [0.2, 0.8])\n",
    "gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f, f)\n",
    "\n",
    "y = multi_binarizer(y_arr, [0.2, 0.8])\n",
    "arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f, f)\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))   \n",
    "\n",
    "y = multi_binarizer(y_gpr, [0.5])\n",
    "gprot_results_t1000_single0pt5_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/gprot_results_t1000_single0pt5_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_single0pt5_split0pt6_logistic_2f, f)\n",
    "\n",
    "y = multi_binarizer(y_arr, [0.5])\n",
    "arr_results_t1000_single0pt5_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/arr_results_t1000_single0pt5_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_single0pt5_split0pt6_logistic_2f, f)\n",
    "    \n",
    "print(np.median(np.array(gprot_results_t1000_single0pt5_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_single0pt5_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "\n",
    "    \n",
    "y = multi_binarizer(y_gpr, [0.2])\n",
    "gprot_results_t1000_single0pt2_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/gprot_results_t1000_single0pt2_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_single0pt2_split0pt6_logistic_2f, f)\n",
    "\n",
    "y = multi_binarizer(y_arr, [0.2])\n",
    "arr_results_t1000_single0pt2_split0pt6_logistic_2f = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False, model=\"logistic_cv\")\n",
    "with open(\"%s/arr_results_t1000_single0pt2_split0pt6_logistic_2f.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_single0pt2_split0pt6_logistic_2f, f)\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_single0pt2_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_single0pt2_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "\n",
    "        \n",
    "gprot_results_t1000_0pt8_ridge = do_regression_experiment(features, y_gpr, feature_names, n_trials, .8, regularize=False, model=\"RidgeCV\")\n",
    "with open(\"%s/gprot_results_t1000_0pt8_ridge.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_0pt8_ridge, f)\n",
    "\n",
    "arr_results_t1000_0pt8_ridge = do_regression_experiment(features, y_arr, feature_names, n_trials, .8, regularize=False, model=\"RidgeCV\")\n",
    "with open(\"%s/arr_results_t1000_0pt8_ridge.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_0pt8_ridge, f)\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_0pt8_ridge[\"test_r2s\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_0pt8_ridge[\"test_r2s\"]), axis=0))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "n_trials = 1000\n",
    "\n",
    "y = multi_binarizer(y_gpr, [0.2, 0.8])\n",
    "gprot_results_t1000_multi0pt2_split0pt8_logistic_cv_2f_obs = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=True, model=\"logistic_cv\")\n",
    "with open(\"%s/gprot_results_t1000_multi0pt2_split0pt8_logistic_cv_2f_obs.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_multi0pt2_split0pt8_logistic_cv_2f_obs, f)\n",
    "\n",
    "y = multi_binarizer(y_arr, [0.2, 0.8])\n",
    "arr_results_t1000_multi0pt2_split0pt8_logistic_cv_2f_obs = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=True, model=\"logistic_cv\")\n",
    "with open(\"%s/arr_results_t1000_multi0pt2_split0pt8_logistic_cv_2f_obs.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_multi0pt2_split0pt8_logistic_cv_2f_obs, f)\n",
    "\n",
    "\n",
    "y = multi_binarizer(y_gpr, [0.2])\n",
    "gprot_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=True, model=\"logistic_cv\")\n",
    "with open(\"%s/gprot_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs, f)\n",
    "\n",
    "y = multi_binarizer(y_arr, [0.2])\n",
    "arr_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=True, model=\"logistic_cv\")\n",
    "with open(\"%s/arr_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs, f)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(a_vs_g_results_t1000_single0pt2_split0pt8_logistic_cv_2f_obs[\"test_roc_aucs\"]), axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arr_results_t100_single0pt2_split0pt9_rfr_2f_obs[\"feature_importances\"][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "analyze_multiclass_experiment(gprot_results_t100_single0pt2_split0pt9_logistic_cv_2f_obs, \n",
    "                              [\"Crystal Structures\", \"MSM Docking\"],\n",
    "                              top_clusters.tolist() + null_features.columns.values.tolist(), common_agonists, analysis_dir,\n",
    "                              [\"GProtein Antagonists\", \"GProtein Agonists\"], X_scaled, \n",
    "                              exp_title=\"GProtein Class\", coef_name=\"Logistic Coefficient\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_features_df.columns.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arr_results_t100_single0pt2_split0pt9_rfr_2f_obs[\"feature_importances\"][0][2].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_features_df.columns.values.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.median(np.array(gprot_results_t1000_0pt8_ridge[\"test_r2s\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_0pt8_ridge[\"test_r2s\"]), axis=0))\n",
    "print(\"/n\")\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_multi0pt33_0pt66_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(\"/n\")\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(\"/n\")\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_single0pt5_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_single0pt5_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(\"/n\")\n",
    "\n",
    "print(np.median(np.array(gprot_results_t1000_single0pt2_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "print(np.median(np.array(arr_results_t1000_single0pt2_split0pt6_logistic_2f[\"test_roc_aucs\"]), axis=0))\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import plots\n",
    "reload(plots)\n",
    "from plots import *\n",
    "\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "\n",
    "\n",
    "feature_names = top_clusters.tolist() +  [\"Inactive Crystal\", \"Active Crystal\", \"Crystal Difference\"] \n",
    "feature_names = [s.replace(\"cluster\", \"MSM State \") for s in feature_names]\n",
    "print(feature_names)\n",
    "\n",
    "#with open(\"%s/gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f.pkl\" %analysis_dir, \"rb\") as f:\n",
    "#    gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f = pickle.load(f)\n",
    "#analyze_multiclass_experiment(gprot_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f, [\"Crystal Structures\", \"MSM Docking\"], feature_names, common_agonists, analysis_dir, [\"G Protein Antagonists\", \"G Protein Partial Agonists\", \"G Protein Full Agonists\"], X_scaled, exp_title=\"GProtein Three Clas 1000 Trials\", coef_name=\"Logistic Coefficient\")\n",
    "\n",
    "\n",
    "#with open(\"%s/arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f.pkl\" %analysis_dir, \"rb\") as f:\n",
    "#    arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f = pickle.load(f)\n",
    "#analyze_multiclass_experiment(arr_results_t1000_multi0pt2_0pt8_split0pt6_logistic_2f, [\"Crystal Structures\", \"MSM Docking\"], feature_names, common_agonists, analysis_dir, [\"Arrestin Antagonists\", \"Arrestin Partial Agonists\", \"Arrestin Full Agonists\"], X_scaled, exp_title=\"Arrestin Three Class 1000 Trials\", coef_name=\"Logistic Coefficient\")\n",
    "\n",
    "#plt.clf()\n",
    "#analyze_multiclass_experiment(gprot_results_t100_single0pt2_split0pt6_logistic_2f, [\"Crystal Structures\", \"MSM Docking\"], feature_names, common_agonists, analysis_dir, [\"G Protein Antagonists\", \"G Protein Agonists\"], X_scaled, exp_title=\"Two Class\", coef_name=\"Logistic Coefficient\")\n",
    "with open(\"%s/arr_results_t1000_single0pt2_split0pt6_logistic_2f.pkl\" %analysis_dir, \"rb\") as f:\n",
    "    arr_results_t1000_single0pt2_split0pt6_logistic_2f = pickle.load(f)\n",
    "analyze_multiclass_experiment(arr_results_t1000_single0pt2_split0pt6_logistic_2f, [\"Crystal Structures\", \"MSM Docking\"],feature_names, common_agonists, analysis_dir, [\"Arrestin Antagonists\", \"Arrestin Agonists\"], X_scaled, exp_title=\"Arrestin Two Class\", coef_name=\"Logistic Coefficient\")\n",
    "\n",
    "with open(\"%s/gprot_results_t1000_single0pt2_split0pt6_logistic_2f.pkl\" %analysis_dir, \"rb\") as f:\n",
    "    gprot_results_t1000_single0pt2_split0pt6_logistic_2f = pickle.load(f)\n",
    "analyze_multiclass_experiment(gprot_results_t1000_single0pt2_split0pt6_logistic_2f, [\"Crystal Structures\", \"MSM Docking\"],feature_names, common_agonists, analysis_dir, [\"G Protein Antagonists\", \"G Protein Agonists\"], X_scaled, exp_title=\"G Protein Two Class\", coef_name=\"Logistic Coefficient\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "analysis_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_scaled[:,25]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ddg_scaled[\"salbutamol\"].sort(inplace=False\n",
    "                             )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arr_results_t100_multi_0pt6_logistic = arr_results_t100_multi_0pt_6_rfr \n",
    "\n",
    "= arr_results_t1000_multi_0pt_6_rfr "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "save_file = \"%s/arr_results_t1000_multi_0pt6_logistic.pkl\" %(analysis_dir)\n",
    "with open(save_file, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_multi_0pt6_logistic, f)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "save_file = \"%s/gprot_results_t1000_multi_0pt6_logistic.pkl\" %(analysis_dir)\n",
    "with open(save_file, \"wb\") as f:\n",
    "    pickle.dump(gprot_results_t1000_multi_0pt6_logistic, f)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(arr_results_t1000_0pt5_md3_ridge[\"test_r2s\"]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(arr_results_t1000_0pt5_md3_ridge[\"kendall_pvalues\"]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.median(np.array(arr_results_t1000_0pt5_md3_ridge[\"kendall_pvalues\"]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(arr_results_t1000_multi_0pt6_logistic[\"test_roc_aucs\"]), axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(y_gpr, bins=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "#arrestin_antagonists = [\"s-carvedilol\", \"nebivolol\"]\n",
    "#non_arrestin_antagonists = [n for n in antagonists if n not in arrestin_antagonists and n not in [\"Carvedilol\"]]\n",
    "#y = np.array([1. for i in arrestin_antagonists] + [0. for i in non_arrestin_antagonists]).reshape((-1,1))\n",
    "\n",
    "\n",
    "total_activity = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands])\n",
    "#common_agonists = arrestin_antagonists + non_arrestin_antagonists\n",
    "#biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\"]\n",
    "#non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "#y = np.array([1. for i in biased_ligands] + [0. for i in non_biased_ligands]).reshape((-1,1))\n",
    "#common_agonists = biased_ligands + non_biased_ligands\n",
    "\n",
    "#common_agonists = total_activity.loc[total_activity > 0.2].index.values\n",
    "#y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].subtract(bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists]).values.reshape((-1,1))\n",
    "\n",
    "top_clusters = delta_delta_g.index.values\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"nebivolol\").index.values[:10].tolist() + delta_delta_g.sort(\"3p0g_lig\").index.values[:10].tolist()))\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"N-Cyclopentylbutanephrine\", inplace=False).index.values[:4].tolist() + delta_delta_g.sort(\"procaterol\", inplace=False).index.values[:4].tolist()))\n",
    "#agonists_df = [a for a in agonists if a in delta_delta_g.columns.values]\n",
    "#common_agonists = agonists_df + antagonists\n",
    "#y = np.array([1. for i in agonists_df] + [0. for i in antagonists]).reshape((-1,1))\n",
    "common_agonists = common_ligands\n",
    "y_arr = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "y_gpr = bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "#y_ori = y_arr - y_gpr\n",
    "#y_ori = y_arr\n",
    "y = multi_binarizer(y_arr, [0.33, 0.66])\n",
    "#y = binarize(y_gpr, threshold=0.2)\n",
    "\n",
    "C = null_features.loc[common_agonists].values\n",
    "#X = delta_delta_g.loc[top_clusters][common_agonists].values.T\n",
    "X = np.hstack([delta_delta_g.loc[top_clusters][common_agonists].values.T, C])\n",
    "#X_scaled = ddg_scaled.loc[top_clusters][common_agonists].values.T\n",
    "X_scaled = np.hstack([ddg_scaled.loc[top_clusters][common_agonists].values.T, C])\n",
    "#D_scaled = docking_normalized.loc[top_clusters][common_agonists].values.T\n",
    "D_scaled = np.hstack([docking_normalized.loc[top_clusters][common_agonists].values.T, C])\n",
    "X_diff = np.zeros((X.shape[0], X.shape[1]**2/2))\n",
    "k=0\n",
    "for i in range(0,X.shape[1]):\n",
    "    for j in range(i+1,X.shape[1]):\n",
    "        X_diff[:,k] = X[:,i] - X[:,j]\n",
    "        k+=1\n",
    "X_diff = np.hstack([X_scaled, X_diff, C])\n",
    "\n",
    "features = [C, D_scaled, X, X_scaled, X_diff]\n",
    "features_y = [C, D_scaled, X, X_scaled, X_diff]\n",
    "feature_names = [\"Crystal Structures\", \"Normalized Docking\", \"Docking ddG\", \"Docking ddg Scaled\", \"Docking ddG Differences\"]\n",
    "\n",
    "n_trials = 1000\n",
    "test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=False)\n",
    "#arr_results_t100_0pt9_md3 = do_regression_experiment(features, y, feature_names, n_trials, 0.9, regularize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "save_file = \"%s/arr_results_t1000_0pt9_md3.pkl\" %analysis_dir\n",
    "with open(save_file, \"wb\") as f:\n",
    "    pickle.dump(arr_results_t1000_0pt9_md3, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.median(np.array(arr_results_t1000_0pt9_md3[\"test_r2s\"]), axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "#arrestin_antagonists = [\"s-carvedilol\", \"nebivolol\"]\n",
    "#non_arrestin_antagonists = [n for n in antagonists if n not in arrestin_antagonists and n not in [\"Carvedilol\"]]\n",
    "#y = np.array([1. for i in arrestin_antagonists] + [0. for i in non_arrestin_antagonists]).reshape((-1,1))\n",
    "\n",
    "\n",
    "total_activity = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands])\n",
    "#common_agonists = arrestin_antagonists + non_arrestin_antagonists\n",
    "#biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\"]\n",
    "#non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "#y = np.array([1. for i in biased_ligands] + [0. for i in non_biased_ligands]).reshape((-1,1))\n",
    "#common_agonists = biased_ligands + non_biased_ligands\n",
    "\n",
    "common_agonists = total_activity.loc[total_activity > 0.2].index.values\n",
    "#y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].subtract(bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists]).values.reshape((-1,1))\n",
    "\n",
    "top_clusters = delta_delta_g.index.values\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"nebivolol\").index.values[:10].tolist() + delta_delta_g.sort(\"3p0g_lig\").index.values[:10].tolist()))\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"N-Cyclopentylbutanephrine\", inplace=False).index.values[:4].tolist() + delta_delta_g.sort(\"procaterol\", inplace=False).index.values[:4].tolist()))\n",
    "#agonists_df = [a for a in agonists if a in delta_delta_g.columns.values]\n",
    "#common_agonists = agonists_df + antagonists\n",
    "#y = np.array([1. for i in agonists_df] + [0. for i in antagonists]).reshape((-1,1))\n",
    "#common_agonists = common_ligands\n",
    "y_arr = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "y_gpr = bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "y_ori = y_arr - y_gpr\n",
    "#y_ori = y_arr\n",
    "#y = y_arr\n",
    "y = binarize(y_ori, threshold=-0.2)\n",
    "\n",
    "C = null_features.loc[common_agonists].values\n",
    "#X = delta_delta_g.loc[top_clusters][common_agonists].values.T\n",
    "X = np.hstack([delta_delta_g.loc[top_clusters][common_agonists].values.T, C])\n",
    "#X_scaled = ddg_scaled.loc[top_clusters][common_agonists].values.T\n",
    "X_scaled = np.hstack([ddg_scaled.loc[top_clusters][common_agonists].values.T, C])\n",
    "#D_scaled = docking_normalized.loc[top_clusters][common_agonists].values.T\n",
    "D_scaled = np.hstack([docking_normalized.loc[top_clusters][common_agonists].values.T, C])\n",
    "X_diff = np.zeros((X.shape[0], X.shape[1]**2/2))\n",
    "k=0\n",
    "for i in range(0,X.shape[1]):\n",
    "    for j in range(i+1,X.shape[1]):\n",
    "        X_diff[:,k] = X[:,i] - X[:,j]\n",
    "        k+=1\n",
    "X_diff = np.hstack([X_scaled, X_diff, C])\n",
    "\n",
    "features = [C, D_scaled, X, X_scaled, X_diff]\n",
    "features_y = [C, D_scaled, X, X_scaled]\n",
    "feature_names = [\"Crystal Structures\", \"Normalized Docking\", \"Docking ddG\", \"Docking ddg Scaled\", \"Docking ddG Differences\"]\n",
    "\n",
    "n_trials = 1000\n",
    "#arr_classification_results_t1000_0pt6 = do_classification_experiment(features, y, feature_names, n_trials, 0.6, regularize=False)\n",
    "arr_vs_gprot_classification_results_t1000_0pt8 = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=False)\n",
    "#arr_results_1000 = do_regression_experiment(features, y, feature_names, n_trials, 0.8, regularize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open(\"%s/arr_vs_gprot_classification_results_t100_0pt8.pkl\" %analysis_dir, \"wb\") as f:\n",
    "    pickle.dump(arr_vs_gprot_classification_results_t100_0pt8, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.median(np.array(arr_vs_gprot_classification_results_t100_0pt8[\"test_aucs\"]), axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.hist(y_ori,bins=50)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "save_file = \"%s/arr_classification_results_t1000_0pt6.pkl\" %analysis_dir\n",
    "with open(save_file, \"wb\") as f:\n",
    "    pickle.dump(arr_classification_results_t1000_0pt6, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(arr_classification_results_t1000_0pt6[\"test_log_aucs\"]),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "#arrestin_antagonists = [\"s-carvedilol\", \"nebivolol\"]\n",
    "#non_arrestin_antagonists = [n for n in antagonists if n not in arrestin_antagonists and n not in [\"Carvedilol\"]]\n",
    "#y = np.array([1. for i in arrestin_antagonists] + [0. for i in non_arrestin_antagonists]).reshape((-1,1))\n",
    "\n",
    "\n",
    "total_activity = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands])\n",
    "#common_agonists = arrestin_antagonists + non_arrestin_antagonists\n",
    "#biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\"]\n",
    "#non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "#y = np.array([1. for i in biased_ligands] + [0. for i in non_biased_ligands]).reshape((-1,1))\n",
    "#common_agonists = biased_ligands + non_biased_ligands\n",
    "\n",
    "#common_agonists = total_activity.loc[total_activity > 0.2].index.values\n",
    "#y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].subtract(bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists]).values.reshape((-1,1))\n",
    "\n",
    "top_clusters = delta_delta_g.index.values\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"nebivolol\").index.values[:10].tolist() + delta_delta_g.sort(\"3p0g_lig\").index.values[:10].tolist()))\n",
    "#top_clusters = list(set(delta_delta_g.sort(\"N-Cyclopentylbutanephrine\", inplace=False).index.values[:4].tolist() + delta_delta_g.sort(\"procaterol\", inplace=False).index.values[:4].tolist()))\n",
    "#agonists_df = [a for a in agonists if a in delta_delta_g.columns.values]\n",
    "#common_agonists = agonists_df + antagonists\n",
    "#y = np.array([1. for i in agonists_df] + [0. for i in antagonists]).reshape((-1,1))\n",
    "common_agonists = common_ligands\n",
    "y_arr = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "y_gpr = bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "#y_ori = y_arr - y_gpr\n",
    "y_ori = y_arr\n",
    "#y = y_arr\n",
    "y = binarize(y_arr, threshold=0.2)\n",
    "\n",
    "C = null_features.loc[common_agonists].values\n",
    "#X = delta_delta_g.loc[top_clusters][common_agonists].values.T\n",
    "X = np.hstack([delta_delta_g.loc[top_clusters][common_agonists].values.T, C])\n",
    "#X_scaled = ddg_scaled.loc[top_clusters][common_agonists].values.T\n",
    "X_scaled = np.hstack([ddg_scaled.loc[top_clusters][common_agonists].values.T, C])\n",
    "#D_scaled = docking_normalized.loc[top_clusters][common_agonists].values.T\n",
    "D_scaled = np.hstack([docking_normalized.loc[top_clusters][common_agonists].values.T, C])\n",
    "X_diff = np.zeros((X.shape[0], X.shape[1]**2/2))\n",
    "k=0\n",
    "for i in range(0,X.shape[1]):\n",
    "    for j in range(i+1,X.shape[1]):\n",
    "        X_diff[:,k] = X[:,i] - X[:,j]\n",
    "        k+=1\n",
    "X_diff = np.hstack([X_scaled, X_diff, C])\n",
    "\n",
    "features = [C, D_scaled, X, X_scaled, X_diff]\n",
    "features_y = [C, D_scaled, X, X_scaled]\n",
    "feature_names = [\"Crystal Structures\", \"Normalized Docking\", \"Docking ddG\", \"Docking ddg Scaled\", \"Docking ddG Differences\"]\n",
    "\n",
    "n_trials = 1000\n",
    "test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances = do_classification_experiment(features, y, feature_names, n_trials, 0.8, regularize=False)\n",
    "#arr_auc_results = do_regression_experiment(features, y, feature_names, n_trials, 0.8, regularize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(test_log_aucs),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(gprot_results[\"test_r2s\"]),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gprot_results[\"test_r2s\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "features_y = copy.deepcopy(features)\n",
    "features_y.append(y)\n",
    "train_test_arrays = train_test_split(*features_y, train_size=0.8)\n",
    "print([a.shape for a in train_test_arrays])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "train_test_arrays = train_test_split(*features_y, train_size=0.8)\n",
    "print([a.shape for a in train_test_arrays])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gprot_results = copy.deepcopy([test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances])\n",
    "#import pickle\n",
    "gprot_file = \"%s/gprot_results_0pt2_XP_no-regularization.pkl\"\n",
    "with open(\"%s/gprot_results_0pt2_XP_no-regularization.pkl\", \"wb\") as f:\n",
    "    pickle.dump(gprot_results, f)\n",
    "with open(gprot_file) as f:\n",
    "    test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#arrestin_results = copy.deepcopy([test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances])\n",
    "#import pickle\n",
    "arrestin_file = \"%s/arrestin_results_0pt2_XP_no-regularization.pkl\"\n",
    "#with open(\"%s/arrestin_results_0pt2_XP_no-regularization.pkl\", \"wb\") as f:\n",
    "#    pickle.dump(arrestin_results, f)\n",
    "with open(arrestin_file) as f:\n",
    "    test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#arrestin_rfr_results = copy.deepcopy([test_r2s, rfr_feature_importances])\n",
    "arrestin_rfr_file = \"%s/arrestin_results_rfr_XP_no-regularization.pkl\"\n",
    "#with open(\"%s/arrestin_results_rfr_XP_no-regularization.pkl\", \"wb\") as f:\n",
    "#    pickle.dump(arrestin_rfr_results, f)\n",
    "with open(arrestin_rfr_file) as f:\n",
    "    test_r2s, rfr_feature_importances = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#gprot_rfr_results = copy.deepcopy([test_r2s, rfr_feature_importances])\n",
    "gprot_rfr_file = \"%s/gprot_results_rfr_XP_no-regularization.pkl\"\n",
    "#with open(\"%s/gprot_results_rfr_XP_no-regularization.pkl\", \"wb\") as f:\n",
    "#    pickle.dump(gprot_rfr_results, f)\n",
    "with open(gprot_rfr_file) as f:\n",
    "    test_r2s, rfr_feature_importances = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test_aucs = np.array(test_aucs)\n",
    "plt.clf()\n",
    "plt.boxplot(test_aucs[:,2]-test_aucs[:,0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test_log_aucs = np.array(test_log_aucs)\n",
    "np.percentile(test_aucs[:,2]-test_aucs[:,0], 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sklearn\n",
    "null_preds = binarize(C[:,2])\n",
    "sklearn.metrics.roc_auc_score(binarize(y, 0.2).ravel(), null_preds.ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test_aucs = np.array(test_aucs)\n",
    "print(np.median(test_aucs, axis=0))\n",
    "print(np.median(test_aucs[:,2] - test_aucs[:,0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_successes = len(np.where(test_aucs[:,2]-test_aucs[:,0] > 0.)[0])\n",
    "nobs = test_aucs.shape[0]\n",
    "statsmodels.stats.proportion.proportion_confint(count=n_successes, nobs=nobs, alpha=0.01, method='wilson')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(np.array(test_r2s),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arrestin_vs_gprot_results = copy.deepcopy([test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances])\n",
    "import pickle\n",
    "with open(\"%s/arrestin_vs_gprot_results_-0pt2_XP.pkl\", \"wb\") as f:\n",
    "    pickle.dump(arrestin_vs_gprot_results, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arrestin_vs_gprot_rfr_results = copy.deepcopy([test_r2s, rfr_feature_importances])\n",
    "import pickle\n",
    "with open(\"%s/arrestin_vs_gprot_rfr_results_XP.pkl\", \"wb\") as f:\n",
    "    pickle.dump(arrestin_vs_gprot_rfr_results, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances = arrestin_vs_gprot_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#gprotein_results = copy.deepcopy([test_accuracies, test_aucs, test_log_aucs, C_test_aucs, C_test_log_aucs, feature_importances])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "results_dict = gprot_results_t1000_multi_0pt6_logistic\n",
    "import efficacy_scripts\n",
    "reload(efficacy_scripts)\n",
    "from efficacy_scripts import *\n",
    "importances_df, results_df = analyze_classification_experiment(results_dict[\"test_roc_aucs\"],\n",
    "                                               results_dict[\"feature_importances\"],\n",
    "                                               feature_names,\n",
    "                                               X_scaled, y, pd.concat([ddg_scaled, null_features.transpose()], axis=0),\n",
    "                                               top_clusters.tolist() + null_features.columns.values.tolist(), \n",
    "                                               common_agonists, \n",
    "                                                \"Predicting G Protein Three Class, 1000 trials, 60:40 Split, XP\", analysis_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "feature_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "importances_df, results_df = analyze_experiment(test_aucs, test_log_aucs, feature_importances, feature_names,\n",
    "                        X_scaled, y, pd.concat([ddg_scaled, null_features.transpose()], axis=0) , top_clusters.tolist() + null_features.columns.values.tolist(), common_agonists, \"Predicting Arrestin, w Crystal Features, 0.2, SP\", analysis_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "importances_df, results_df = analyze_regression_experiment(test_r2s, rfr_feature_importances, feature_names,\n",
    "                        X, y, pd.concat([delta_delta_g, null_features.transpose()], axis=0) , top_clusters.tolist() + null_features.columns.values.tolist(), common_agonists, \"Predicting Arrestin, RFR, w Crystal Features, SP\", analysis_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_top_measurable_features(samples_normalized_features_avg_df, cluster_name):\n",
    "    import re\n",
    "    top_features_cluster = []\n",
    "    #top_features_cluster = samples_normalized_features_avg_df.loc[\"cluster_name\"].loc[samples_normalized_features_avg_df.loc[\"cluster_name\"].abs() > .75].index.values\n",
    "    #print(top_features_cluster)\n",
    "    #print(len(top_features_cluster))\n",
    "    [top_features_cluster.append(pair) for pair in samples_normalized_features_avg_df.loc[cluster_name].abs().sort(inplace=False, ascending=False).index.values[:10]]\n",
    "    all_features = []\n",
    "    features = []\n",
    "    for f in top_features_cluster:\n",
    "        fs = f.split(\",\")\n",
    "        for i in range(0, len(fs)):\n",
    "            res = int(re.findall(r'\\d+', fs[i])[0])\n",
    "            all_features.append(res)\n",
    "            if \"TRP\" in fs[i] or \"CYS\" in fs[i] or \"TYR\" in fs[i] or \"LYS\" in fs[i]:\n",
    "                features.append(res)\n",
    "            \n",
    "    top_features_cluster = sorted(list(set(features)))\n",
    "    #print(sorted(list(set([int(re.findall(r'\\d+', r)[0]) for r in all_features]))))\n",
    "    #print(top_features_cluster)\n",
    "    return top_features_cluster, all_features\n",
    "a, b = get_top_measurable_features(samples_normalized_features_averages_df, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ddg_scaled[[\"norepinephrine\", \"r_epinephrine\", \"ethylnorepinephrine\", \"nebivolol\", \"s-carvedilol\", \"s-carazolol\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.mean(test_aucs, axis=0))\n",
    "print(np.median(test_aucs, axis=0))\n",
    "print(np.mean(test_log_aucs, axis=0))\n",
    "print(np.median(test_log_aucs, axis=0))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(ddg_scaled.loc[\"cluster2\"][common_agonists], y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_ori = bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "y = binarize(y_ori, threshold=0.05)\n",
    "\n",
    "X = delta_delta_g.loc[top_clusters][common_agonists].values.T\n",
    "X_scaled = ddg_scaled.loc[top_clusters][common_agonists].values.T\n",
    "C = null_features.loc[common_agonists].values\n",
    "D_scaled = docking_normalized.loc[top_clusters][common_agonists].values.T\n",
    "\n",
    "X_train = D_scaled\n",
    "y_train = y\n",
    "f = np.zeros(X_train.shape[1])\n",
    "\n",
    "rfr = RandomForestClassifier(n_estimators=1000, max_features='sqrt', n_jobs=-1, oob_score=True)\n",
    "rfr.fit(X_train, y_train)\n",
    "\n",
    "#top_indices = np.argsort(rfr.feature_importances_*-1.)[:min(20, X_train.shape[1])]\n",
    "#rfr = RandomForestClassifier(n_estimators=100, max_features='sqrt', n_jobs=-1, oob_score=True)\n",
    "#X_train = X_train[:, top_indices]\n",
    "#rfr.fit(X_train, y_train)\n",
    "#f[top_indices] = rfr.feature_importances_\n",
    "#y_pred = rfr.predict(X_train)\n",
    "#y_score = rfr.predict_proba(X_train)\n",
    "top_indices=range(0,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_scaled.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#G Protein, Agonist Results\n",
    "test_drugs = secret_compounds + [\"nebivolol\", \"s-carvedilol\", \"xamoterol\", \"3p0g_lig\", \"isoetharine\", \"ethylnorepinephrine\", \"N-Cyclopentylbutanephrine\", \"ta-2005\", \"procaterol\"]\n",
    "X_test = docking_normalized.transpose().loc[test_drugs].values[:, top_indices]\n",
    "print(X_test)\n",
    "pd.DataFrame(rfr.predict_proba(X_test), index=test_drugs, columns=[\"P(Antagonist\", \"P(Agonist)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ARRESTIN, Agonist Results\n",
    "X_test = docking_normalized.transpose().loc[test_drugs].values[:, top_indices]\n",
    "X_test.shape\n",
    "pd.DataFrame(rfr.predict_proba(X_test), index=test_drugs, columns=[\"P(Antagonist\", \"P(Agonist)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#G Protein, Full Agonist Results\n",
    "test_drugs = secret_compounds + [\"nebivolol\", \"s-carvedilol\", \"xamoterol\", \"3p0g_lig\", \"isoetharine\", \"ethylnorepinephrine\", \"N-Cyclopentylbutanephrine\"]\n",
    "X_test = docking_normalized.transpose().loc[test_drugs].values[:, top_indices]\n",
    "pd.DataFrame(rfr.predict_proba(X_test), index=test_drugs, columns=[\"P(Antagonist\", \"P(Agonist)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ARRESTIN, Agonist Results\n",
    "X_test = docking_normalized.transpose().loc[test_drugs].values[:, top_indices]\n",
    "X_test.shape\n",
    "pd.DataFrame(rfr.predict_proba(X_test), index=test_drugs, columns=[\"P(Antagonist\", \"P(Agonist)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "common_agonists = total_activity.loc[total_activity > 0.2].index.values\n",
    "plt.scatter(docking_normalized.loc[\"cluster21\"][common_agonists], bret[\"B2AR-Gprotein, Mean\"].subtract(bret[\"B2AR-Arrestin, Mean\"])[common_agonists])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "len(common_ligands)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ARRESTIN, Partial Agonist Results\n",
    "X_test = docking_normalized.transpose().loc[test_drugs].values[:, top_indices]\n",
    "X_test.shape\n",
    "pd.DataFrame(rfr.predict_proba(X_test), index=test_drugs, columns=[\"P(Antagonist\", \"P(Agonist)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#ARRESTIN, Full Agonist Results\n",
    "X_test = docking_normalized.transpose().loc[test_drugs].values[:, top_indices]\n",
    "X_test.shape\n",
    "pd.DataFrame(rfr.predict_proba(X_test), index=test_drugs, columns=[\"P(Antagonist\", \"P(Agonist)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top_indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ddg_scaled.transpose().loc[secret_compounds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_clustermap(docking_normalized[secret_compounds].iloc[top_indices].transpose(), save_file=\"%s/mehrdad_clustermap.pdf\" %(save_dir), method='average', z_score=None)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.median(test_aucs, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(ddg_scaled.loc[\"cluster11\"][common_agonists], y)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arrestin_top = [16, 80, 43, 21, 84, 38, 44, 6, 13, 99]\n",
    "gprot_top = [44, 6, 83, 4, 76, 99, 62, 92, 39, 80]\n",
    "\n",
    "arrestin_only = sorted(list(set(arrestin_top).difference(set(gprot_top))))\n",
    "print(arrestin_only)\n",
    "gprot_only = sorted(list(set(gprot_top).difference(set(arrestin_top))))\n",
    "print(gprot_only)\n",
    "both = sorted(list(set(arrestin_top).intersection(set(gprot_top))))\n",
    "print(both)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples_pnas_tica.loc[importances_df.index.values.tolist()[:5]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "importances_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sklearn\n",
    "reload(sklearn)\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.preprocessing import binarize\n",
    "from sklearn import cross_validation\n",
    "\n",
    "test_accuracies = []\n",
    "test_aucs = []\n",
    "test_log_aucs = []\n",
    "C_test_aucs = []\n",
    "C_test_log_aucs = []\n",
    "n_trials = 10\n",
    "feature_importances = []\n",
    "total_activity = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands])\n",
    "common_agonists = total_activity.loc[total_activity > 0.1].index.values\n",
    "\n",
    "#biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\"]\n",
    "#biased_ligands += [\"nebivolol\", \"s-carvedilol\"]\n",
    "#non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "#non_biased_ligands += [\"s-carazolol\", \"Ici215001\", \"bisoprolol\", \"timolol\", \"s-atenolol\"]\n",
    "#non_biased_ligands = [n for n in df.columns.values.tolist() if n not in biased_ligands and \"Carvedilol\" not in n]\n",
    "#common_agonists = biased_ligands + non_biased_ligands\n",
    "#common_agonists = common_ligands\n",
    "#top_clusters = [\"cluster80\", \"cluster62\", \"cluster11\", \"cluster21\", \"cluster16\", \"cluster43\", \"cluster38\"]\n",
    "#differences = np.zeros((len(top_clusters), len(top_clusters)))\n",
    "#for i, cluster in enumerate(top_clusters):\n",
    "#    for j in range(i, len(top_clusters)):\n",
    "#        differences[i][j] = \n",
    "\n",
    "#top_clusters = importances_df.index.values\n",
    "top_clusters = delta_delta_g.index.values\n",
    "\n",
    "y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].subtract(bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists]).values.reshape((-1,1))\n",
    "#y_ori = np.vstack([y_ori, np.ones(3).reshape((-1,1))])\n",
    "#common_agonists = common_agonists.tolist() + biased_ligands\n",
    "\n",
    "X = delta_delta_g.loc[top_clusters][common_agonists].values.T\n",
    "X_scaled = ddg_scaled.loc[top_clusters][common_agonists].values.T\n",
    "C = null_features.loc[common_agonists].values\n",
    "D_scaled = docking_normalized.loc[top_clusters][common_agonists].values.T\n",
    "\n",
    "#y = np.array([1. for n in biased_ligands] + [0. for n in non_biased_ligands]).reshape((-1,1))\n",
    "#print(y_ori)\n",
    "#y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "#plt.hist(y_ori, bins=25)\n",
    "y = binarize(y_ori, threshold=-0.2)\n",
    "\n",
    "features = [C, X_scaled, D_scaled]\n",
    "features_y = [C, X_scaled, D_scaled, y]\n",
    "feature_names = [\"Crystal Structures\", \"MSM ddG\", \"Docking\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "xt = ddg_scaled[biased_ligands].values.T\n",
    "xt_preds = []\n",
    "\n",
    "\n",
    "for j in range(0,n_trials):\n",
    "    print(j)\n",
    "    aucs = []\n",
    "    log_aucs = []\n",
    "    train_test_arrays = sklearn.cross_validation.train_test_split(*features_y, train_size=0.8, stratify=y) \n",
    "    y_train = train_test_arrays[2*len(features)]\n",
    "    y_test = train_test_arrays[2*len(features) + 1]\n",
    "    feature_importance = []\n",
    "    \n",
    "    for i in range(0, len(features)):\n",
    "        X_train = train_test_arrays[2*i]\n",
    "        X_test = train_test_arrays[2*i+1]\n",
    "\n",
    "        sc = StandardScaler()\n",
    "        sc.fit(X_train)\n",
    "        X_train = sc.transform(X_train)\n",
    "        X_test = sc.transform(X_test)\n",
    "\n",
    "        rfr = RandomForestClassifier(n_estimators=100, max_features='sqrt', max_depth=3, n_jobs=-1, oob_score=True)\n",
    "        rfr.fit(X_train, y_train)\n",
    "        #top_indices = np.argsort(rfr.feature_importances_*-1.)[:min(10, X.shape[1])]\n",
    "        feature_importance.append(rfr.feature_importances_)\n",
    "        #rfr = RandomForestClassifier(n_estimators=10, max_features=None, n_jobs=-1, oob_score=True)\n",
    "        #X_train = X_train[:, top_indices]\n",
    "        #X_test = X_test[:, top_indices]\n",
    "        #rfr.fit(X_train, y_train)\n",
    "        #f = np.zeros(X.shape[1])\n",
    "        #f[top_indices] = rfr.feature_importances_\n",
    "        #feature_importance.append(f)\n",
    "        \n",
    "        if i == 1:\n",
    "            xt_preds.append(rfr.predict(xt))\n",
    "        \n",
    "        y_pred = rfr.predict(X_test)\n",
    "        y_score = rfr.predict_proba(X_test)\n",
    "        auc, logauc = compute_auc(y_test, y_score)\n",
    "        aucs.append(auc)\n",
    "        log_aucs.append(logauc)  \n",
    "    feature_importances.append(feature_importance)\n",
    "    test_aucs.append(aucs)\n",
    "    test_log_aucs.append(log_aucs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\"]\n",
    "\n",
    "non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "\n",
    "ddg_scaled.loc[importances_df.index.values[:5]][biased_ligands + non_biased_ligands]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ddg_scaled.sort(\"procaterol\", inplace=False).iloc[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(ddg_scaled.loc[\"cluster36\"][common_agonists], y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import sklearn\n",
    "reload(sklearn)\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.preprocessing import binarize\n",
    "from sklearn import cross_validation\n",
    "\n",
    "test_accuracies = []\n",
    "test_aucs = []\n",
    "test_log_aucs = []\n",
    "C_test_aucs = []\n",
    "C_test_log_aucs = []\n",
    "n_trials = 1000\n",
    "feature_importances = []\n",
    "reg = []\n",
    "total_activity = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands])\n",
    "common_agonists = total_activity.loc[total_activity > 0.3].index.values\n",
    "\n",
    "#biased_ligands = [\"ethylnorepinephrine\", \"isoetharine\", \"N-Cyclopentylbutanephrine\", \"3p0g_lig\"]\n",
    "#biased_ligands = [\"nebivolol\", \"s-carvedilol\"]\n",
    "#non_biased_ligands =  [\"r_isopreterenol\", \"r_epinephrine\", \"norepinephrine\", \"zinterol\", \"orciprenaline\", \"epinine\", \"terbutaline\", \"fenoterol\", \"procaterol\", \"formoterol\", \"salbutamol\", \"salmeterol\"]\n",
    "#non_biased_ligands = [\"s-carazolol\", \"Ici215001\", \"bisoprolol\", \"timolol\", \"s-atenolol\"]\n",
    "#non_biased_ligands = [n for n in df.columns.values.tolist() if n not in biased_ligands and \"Carvedilol\" not in n]\n",
    "#common_agonists = biased_ligands + non_biased_ligands\n",
    "\n",
    "X = delta_delta_g[common_agonists].values.T\n",
    "X_scaled = ddg_scaled[common_agonists].values.T\n",
    "C = null_features.loc[common_agonists].values\n",
    "\n",
    "#y = np.array([1. for n in biased_ligands] + [0. for n in non_biased_ligands]).reshape((-1,1))\n",
    "#y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].divide(bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].add(bret[\"B2AR-Gprotein, Mean\"].loc[common_agonists])).values.reshape((-1,1))\n",
    "y_ori = bret[\"B2AR-Arrestin, Mean\"].loc[common_agonists].values.reshape((-1,1))\n",
    "plt.hist(y_ori, bins=25)\n",
    "y = binarize(y_ori, threshold=0.2) \n",
    "\n",
    "features = [C, X, X_scaled, docking_normalized[common_agonists].values.T]\n",
    "features_y = [C, X, X_scaled, docking_normalized[common_agonists].values.T, y]\n",
    "feature_names = [\"Crystal Structures\", \"MSM States\", \"Normalized MSM States\", \"Normalized Docking\"]\n",
    "  \n",
    "\n",
    "for j in range(0,n_trials):\n",
    "    print(j)\n",
    "    aucs = []\n",
    "    log_aucs = []\n",
    "    train_test_arrays = sklearn.cross_validation.train_test_split(*features_y, train_size=0.8, stratify=y) \n",
    "    y_train = train_test_arrays[2*len(features)]\n",
    "    y_test = train_test_arrays[2*len(features) + 1]\n",
    "    feature_importance = []\n",
    "    r = []\n",
    "    \n",
    "    for i in range(0, len(features)):\n",
    "        X_train = train_test_arrays[2*i]\n",
    "        X_test = train_test_arrays[2*i+1]\n",
    "\n",
    "        sc = StandardScaler()\n",
    "        sc.fit(X_train)\n",
    "        X_train = sc.transform(X_train)\n",
    "        X_test = sc.transform(X_test)\n",
    "        \n",
    "        cs = np.logspace(-3., 20.)\n",
    "        rfr = LogisticRegressionCV(Cs=cs, penalty='l2')\n",
    "        rfr.fit(X_train, y_train)\n",
    "        feature_importance.append(rfr.coef_)\n",
    "        y_pred = rfr.predict(X_test)\n",
    "        y_score = rfr.predict_proba(X_test)\n",
    "        auc, logauc = compute_auc(y_test, y_score)\n",
    "        aucs.append(auc)\n",
    "        log_aucs.append(logauc)  \n",
    "        r.append(rfr.C_)\n",
    "    reg.append(r)\n",
    "    feature_importances.append(feature_importance)\n",
    "    test_aucs.append(aucs)\n",
    "    test_log_aucs.append(log_aucs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(docking_normalized.loc[\"cluster80\"][common_ligands], -1.0*bret.loc[common_ligands][\"B2AR-Gprotein, Mean\"].subtract(bret.loc[common_ligands][\"B2AR-Arrestin, Mean\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "auc_df = pd.DataFrame(np.array(test_aucs), columns=feature_names)\n",
    "auc_df.plot(kind='box')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "auc_df.median(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "normalized_docking_importances = [f[1] for f in feature_importances]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "importances_df = make_importances_df(normalized_docking_importances, top_clusters)\n",
    "importances_df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import binarize\n",
    "X = ddg_scaled[common_ligands].values.T\n",
    "y = bret[\"B2AR-Gprotein, Mean\"].loc[common_ligands].values.reshape((-1,1))\n",
    "y = binarize(y, threshold=0.5)\n",
    "print(np.shape(y))\n",
    "from sklearn.svm import l1_min_c\n",
    "from sklearn import linear_model\n",
    "\n",
    "#cs = l1_min_c(X, y, loss='log') * np.logspace(0, 3)\n",
    "cs = np.logspace(-3., 20.)\n",
    "print(\"Computing regularization path ...\")\n",
    "clf = linear_model.LogisticRegression(C=1.0, penalty='l2', tol=1e-6)\n",
    "coefs_ = []\n",
    "for c in cs:\n",
    "    clf.set_params(C=c)\n",
    "    clf.fit(X, y)\n",
    "    coefs_.append(clf.coef_.ravel().copy())\n",
    "\n",
    "coefs_ = pd.DataFrame(np.array(coefs_), columns=ddg_scaled.index.values, index=np.log10(cs))\n",
    "#coefs_[list(set(inactive_clusters.tolist()).intersection(set(importances_df.iloc[10:20].index.values.tolist())))].plot()\n",
    "coefs_[importances_df.index.values[:5]].plot()\n",
    "#plt.plot(np.log10(cs), coefs_)\n",
    "#ymin, ymax = plt.ylim()\n",
    "###plt.xlabel('log(C)')\n",
    "#plt.ylabel('Coefficients')\n",
    "##plt.title('Logistic Regression Path')\n",
    "#plt.axis('tight')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples_pnas_tica.loc[importances_df.index.values[:10]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(deltas_tica.loc[importances_df.index.values[:10]][\"tIC.6\"], coefs_[importances_df.index.values[:10]].values[49])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_clustermap(ddg_scaled[common_agonists.tolist()].loc[importances_df.index.values.tolist()[:5]].transpose(), save_file=\"%s/msm_n-clusters%d_lag-time%d_tICs%d.pdf\" %(tica_dir, n_clusters, msm_lag_time, n_components), method='average')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples_pnas_tica.loc[samples_pnas_tica[\"tm6_tm3_dist\"] < 18.0].loc[importances_df.iloc[0:5].index].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ddg_scaled[\"nebivolol\"].subtract(ddg_scaled[\"s-carazolol\"]).sort(inplace=False).iloc[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ddg_scaled[\"s-carvedilol\"].subtract(ddg_scaled[\"s-carazolol\"]).sort(inplace=False).iloc[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples_pnas_tica.loc[[\"cluster74\", \"cluster69\", \"cluster13\", \"cluster12\", \"cluster66\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import plots\n",
    "reload(plots)\n",
    "from plots import *\n",
    "#plot_importances_barh(importances_df.values, importances_df.index.values, \"MSM State Importance in Arrestin Prediction\", \"Feature Importance\", \"MSM State\", \"%s/arrestin_0pt5_classification_rfr.pdf\" %(tica_dir), n_features=50)\n",
    "importances_df.iloc[0:25].plot(kind='barh')\n",
    "plt.xlabel(\"Feature Importance\")\n",
    "plt.ylabel(\"MSM State\")\n",
    "plt.title(\"Importance of MSM States in Predicting Arrestin Activity\")\n",
    "plt.savefig(\"%s/msm_%dstates_arrestin_0pt2_agonists_only_classification_rfr.pdf\" %(tica_dir, n_clusters))\n",
    "#plt.clf()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = docking_normalized[common_agonists].values.T\n",
    "y = bret[\"B2AR-Arrestin, Mean\"].subtract(bret[\"B2AR-Gprotein, Mean\"]).loc[common_agonists].values.reshape((-1,1))\n",
    "y = binarize(y, threshold=-0.2)\n",
    "print(y)\n",
    "from sklearn import linear_model\n",
    "\n",
    "cs = np.logspace(-3., 200.)\n",
    "clf = linear_model.LogisticRegression(C=1.0, penalty='l2', tol=1e-6)\n",
    "coefs_ = []\n",
    "for c in cs:\n",
    "    clf.set_params(C=c)\n",
    "    clf.fit(X, y)\n",
    "    coefs_.append(clf.coef_.ravel().copy())\n",
    "\n",
    "coefs_ = pd.DataFrame(np.array(coefs_), columns=ddg_scaled.index.values, index=np.log10(cs))\n",
    "coefs_[importances_df.iloc[0:10].loc[samples_pnas_tica[\"tm6_tm3_dist\"] < 9.].index].plot()\n",
    "plt.xlabel(\"Log Regularization Parameter\")\n",
    "plt.ylabel(\"Coefficient for Arrestin Activity\")\n",
    "plt.title(\"Logistic Regression Coefficient in Predicting Arrestin Activity\")\n",
    "plt.savefig(\"%s/msm_%dstates_arrestin_0pt2_agonists_only_classification_logistic.pdf\" %(tica_dir, n_clusters))\n",
    "#plt.clf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.shape(coefs_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.median(np.nan_to_num(test_aucs)))\n",
    "print(np.median(np.nan_to_num(C_test_aucs)))\n",
    "print(np.median(np.nan_to_num(test_log_aucs)))\n",
    "print(np.median(np.nan_to_num(C_test_log_aucs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lr = LinearRegression()\n",
    "states = importances_df.index.values.tolist()\n",
    "model = lr.fit(X, y_ori)\n",
    "pd.DataFrame(model.coef_.T, index=delta_delta_g.index, columns=[\"importance\"]).loc[states]#.sort(\"importance\", inplace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model.coef_.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import seaborn\n",
    "reload(seaborn)\n",
    "import seaborn as sns\n",
    "plt.style.use('ggplot')\n",
    "plt.figure(figsize=(5, 5))\n",
    "sns.set_style(\"darkgrid\")\n",
    "g = (auc_df\n",
    "    .pipe((sns.boxplot, 'data'), orient='v', showfliers=True))\n",
    "g.set_xticklabels(auc_df.columns.values, rotation=90)\n",
    "sns.despine()\n",
    "plt.title(\"AUC for Arrestin Prediction\")\n",
    "plt.ylabel(\"Frequency AUCs over Random Splits\")\n",
    "plt.xlabel(\"Featurization\")\n",
    "plt.show()\n",
    "plt.savefig(\"%s/auc_arrestin_prediction_all_ligands_0pt2_cutoff.pdf\" %tica_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "corr_matrix = compute_pearson_matrix(delta_delta_g[common_agonists].values.T, y)\n",
    "corr_df = pd.DataFrame(model.coef_.T, index=delta_delta_g.index.values, columns=[\"Correlation\"]).sort(\"Correlation\",inplace=False)\n",
    "#corr_df.loc[[\"cluster80\", \"cluster16\", \"cluster43\", \"cluster44\"]].plot(kind='barh')\n",
    "corr_df.loc[importances_df.index.values[:20]].sort(\"Correlation\", inplace=False).plot(kind='barh')#, figsize=(5,20))\n",
    "plt.xlabel(\"Pearsson Correlation with Arrestin Activity\")\n",
    "plt.ylabel(\"MSM State\")\n",
    "plt.title(\"Correlation of MSM States with Arrestin Activity\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samples_pnas_tica.loc[corr_df.loc[importances_df.index.values[:20]].sort(\"Correlation\", inplace=False).index.values]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = copy.deepcopy(aggregate_docking_msm)\n",
    "df[df.columns.values] = scale(df.values)\n",
    "plt.scatter(df[common_ligands].loc[\"cluster13\"].values, bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].values.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import sklearn\n",
    "reload(sklearn)\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.preprocessing import binarize\n",
    "from sklearn import cross_validation\n",
    "\n",
    "test_accuracies = []\n",
    "test_aucs = []\n",
    "test_log_aucs = []\n",
    "C_test_aucs = []\n",
    "C_test_log_aucs = []\n",
    "n_trials = 100\n",
    "feature_importances = []\n",
    "\n",
    "for j in range(0,n_trials):\n",
    "    print(j)\n",
    "\n",
    "    X = delta_delta_g[common_ligands].values.T\n",
    "    C = null_features.loc[common_ligands].values\n",
    "    y = bret[\"B2AR-Arrestin, Mean\"].loc[common_ligands].values.reshape((-1,1))\n",
    "    y = binarize(y, threshold=0.5)\n",
    "\n",
    "    X_train, X_test, y_train, y_test, C_train, C_test = sklearn.cross_validation.train_test_split(X, y, C, train_size=0.8, stratify=y)\n",
    "    \n",
    "    sc = StandardScaler()\n",
    "    sc.fit(X_train)\n",
    "    X_train = sc.transform(X_train)\n",
    "    X_test = sc.transform(X_test)\n",
    "\n",
    "    sc = StandardScaler()\n",
    "    sc.fit(C_train)\n",
    "    C_train = sc.transform(C_train)\n",
    "    C_test = sc.transform(C_test)\n",
    " \n",
    "    rfr = RandomForestClassifier(n_estimators=100, max_features='sqrt', max_depth=2, n_jobs=-1, oob_score=True)\n",
    "    rfr.fit(X_train, y_train)\n",
    "    feature_importances.append(rfr.feature_importances_)\n",
    "    y_pred = rfr.predict(X_test)\n",
    "    test_accuracies.append(np.sqrt(np.mean(np.square(y_test-y_pred.reshape((-1,1))))))\n",
    "    y_score = rfr.predict_proba(X_test)\n",
    "    auc, logauc = compute_auc(y_test, y_score)\n",
    "    test_aucs.append(auc)\n",
    "    test_log_aucs.append(logauc)\n",
    "    \n",
    "    rfr = RandomForestClassifier(n_estimators=100, max_features='sqrt', max_depth=2, n_jobs=-1, oob_score=True)\n",
    "    rfr.fit(C_train, y_train)\n",
    "    C_y_pred = rfr.predict(C_test)\n",
    "    y_score = rfr.predict_proba(C_test)\n",
    "    auc, logauc = compute_auc(y_test, y_score)\n",
    "    C_test_aucs.append(auc)\n",
    "    C_test_log_aucs.append(logauc)\n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import seaborn\n",
    "reload(seaborn)\n",
    "import seaborn as sns\n",
    "auc_df = pd.DataFrame(np.vstack([test_aucs, C_test_aucs]).T, columns=[\"MSM States\", \"Crystal Structures\"])\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.figure(figsize=(5, 5))\n",
    "sns.set_style(\"darkgrid\")\n",
    "g = (auc_df\n",
    "    .pipe((sns.boxplot, 'data'), orient='v', showfliers=True))\n",
    "#g.set_xticklabels(experiments.columns.values, rotation=90)\n",
    "sns.despine()\n",
    "plt.title(\"AUC for G Protein Prediction\")\n",
    "plt.ylabel(\"Frequency AUCs over Random Splits\")\n",
    "plt.xlabel(\"Featurization\")\n",
    "plt.show()\n",
    "plt.savefig(\"%s/msm_n-states%d_auc_gprot_prediction_cutoff0pt5.pdf\" %(tica_dir, n_clusters))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "importances_df = make_importances_df(feature_importances, delta_delta_g.index.values.tolist())\n",
    "importances_df\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sklearn\n",
    "from sklearn.linear_model import LinearRegression, Lasso\n",
    "model = Lasso(alpha=0.0001)                                \n",
    "model.fit(ddg_scaled[common_agonists].values.T, bret.loc[common_agonists][\"B2AR-Arrestin, Mean\"].subtract(bret.loc[common_agonists][\"B2AR-Gprotein, Mean\"].values))\n",
    "pd.DataFrame(model.coef_, index=ddg_scaled.index, columns=[\"importance\"]).sort(\"importance\", inplace=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(np.median(np.nan_to_num(test_aucs)))\n",
    "print(np.median(np.nan_to_num(C_test_aucs)))\n",
    "print(np.median(np.nan_to_num(test_log_aucs)))\n",
    "print(np.median(np.nan_to_num(C_test_log_aucs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import plots\n",
    "reload(plots)\n",
    "from plots import *\n",
    "#plot_importances_barh(importances_df.values, importances_df.index.values, \"MSM State Importance in Arrestin Prediction\", \"Feature Importance\", \"MSM State\", \"%s/arrestin_0pt5_classification_rfr.pdf\" %(tica_dir), n_features=50)\n",
    "importances_df.iloc[0:25].plot(kind='barh')\n",
    "plt.xlabel(\"Feature Importance\")\n",
    "plt.ylabel(\"MSM State\")\n",
    "plt.title(\"Importance of MSM States in Predicting G Protein Activity\")\n",
    "#plt.savefig(\"%s/msm_%dstates_gprot_0pt5_classification_rfr.pdf\" %(tica_dir, n_clusters))\n",
    "#plt.clf()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "corr_matrix = compute_pearson_matrix(ddg_scaled[common_ligands].values.T, y)\n",
    "corr_df = pd.DataFrame(model.coef_, index=ddg_scaled.index.values, columns=[\"Correlation\"])\n",
    "corr_df.loc[importances_df.index.values[:10]].sort(\"Correlation\", inplace=False).plot(kind='barh')\n",
    "plt.xlabel(\"Pearsson Correlation with G Protein Activity\")\n",
    "plt.ylabel(\"MSM State\")\n",
    "plt.title(\"Correlation of MSM States with G Protein Activity\")\n",
    "#plt.savefig(\"%s/msm_%dstates_gprot_0pt5_classification_correlations.pdf\" %(tica_dir, n_clusters))\n",
    "#plt.clf()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pd.set_option('display.max_rows', 500)\n",
    "pd.set_option('display.max_columns', 500)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}