FP_Pheno_Networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.spatial.distance import squareform, pdist,jaccard\n",
    "import sys,  datetime, os\n",
    "import requests\n",
    "import time\n",
    "from scipy.stats.stats import pearsonr\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import metrics\n",
    "import seaborn as sns\n",
    "from scipy import stats, integrate\n",
    "from sklearn.metrics import pairwise_distances\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load FPS and Resouces Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fps = []\n",
    "files = []\n",
    "\n",
    "for filename in os.listdir('RDKit_fps/'):\n",
    "    fps.append(filename[:-24])\n",
    "    df = pd.read_table('RDKit_fps/'+ filename)\n",
    "    df[df.columns[0]] = df[df.columns[0]].astype(str)\n",
    "    files.append(df.set_index(df.columns[0]))\n",
    "\n",
    "fpsfiles_dict=dict(zip(fps, files)) \n",
    "del fpsfiles_dict['']\n",
    "fps.remove('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['AtomPair', 'Avalon']"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fps = sorted(fps)\n",
    "fps[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = []\n",
    "files = []\n",
    "\n",
    "for filename in os.listdir('All_bmat/'):\n",
    "    names.append(filename[:-12])\n",
    "    df = pd.read_table('All_bmat/'+ filename)\n",
    "    files.append(df.set_index(df.columns[0]))\n",
    "\n",
    "namefiles_dict=dict(zip(names, files))   \n",
    "del namefiles_dict['']\n",
    "names.remove('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['DrugBank_Targets', 'PharmagKB_SE']"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "names[:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20163"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "name_pcid = pd.read_csv('Input/L1000_Drugs_metadata.csv')\n",
    "\n",
    "name_pcid_dict = {}\n",
    "for index, row in name_pcid.iterrows():\n",
    "    name_pcid_dict[row.loc['pubchem_cid']] = row.loc['pert_iname']\n",
    "len(name_pcid_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions for creating the network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pairwise_dis(df, metric):\n",
    "    array_matrix = metrics.pairwise_distances(df, metric = metric)\n",
    "    return array_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_adjcency_mat(X, metric='euclidean'):\n",
    "    '''X is a samples by features matrix.'''\n",
    "    pdist = pairwise_distances(X, metric=metric)\n",
    "    adj_mat = 1 - pdist / pdist.max()\n",
    "    # remove 1's on the diagonal\n",
    "    adj_mat -= np.eye(X.shape[0])\n",
    "    return adj_mat\n",
    "\n",
    "def create_graph_by_threshold(adj_mat, percentile):\n",
    "    '''This function convert an adjacency matrix to a Graph object \n",
    "    by keeping the drug-drug connections in top percentile.\n",
    "    '''\n",
    "    triu_idx = np.tril_indices(adj_mat.shape[0], 1)\n",
    "    thresold = np.percentile(adj_mat[triu_idx], percentile)\n",
    "    adj_mat_ = adj_mat.copy()\n",
    "    adj_mat_[adj_mat<thresold] = 0\n",
    "    G = nx.from_numpy_matrix(adj_mat_)\n",
    "    return G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a network for each drug with each resource"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "count = 0\n",
    "\n",
    "for name in names:\n",
    "    for fp in fps:\n",
    "        # Load datasets\n",
    "        X1 = fpsfiles_dict[fp] # Dataset1, np.array, e.g. drugs by fingerprints\n",
    "        X2 = namefiles_dict[name]\n",
    "        X2= X2.T # Dataset2, np.array, e.g. drugs by genes\n",
    "        # # make sure the rows of X1 and X2 are the drugs in the same order\n",
    "\n",
    "        shared_drugs = sorted(list(set(X1.index) & set(X2.index)))\n",
    "        for i in shared_drugs:\n",
    "            if i not in list(name_pcid_dict.keys()):\n",
    "                shared_drugs.remove(i)\n",
    "        X1 = X1.loc[shared_drugs]\n",
    "        X2 = X2.loc[shared_drugs]\n",
    "        \n",
    "        # distance metric\n",
    "        A1 = compute_adjcency_mat(X1,'manhattan')\n",
    "        A2 = compute_adjcency_mat(X2, 'cosine')\n",
    "        # Convert to Graphs\n",
    "        percentile = 99.7 # try adjust this param \n",
    "        G1 = create_graph_by_threshold(A1, percentile)\n",
    "#         print(G1.number_of_nodes(), G1.number_of_edges())\n",
    "\n",
    "        G2 = create_graph_by_threshold(A2, percentile)\n",
    "#         print(G2.number_of_nodes(), G2.number_of_edges())\n",
    "\n",
    "        G12 = nx.intersection(G1, G2)\n",
    "        nodes_wo_edges = [n for n, k in G12.degree() if k == 0]\n",
    "        G12.remove_nodes_from(nodes_wo_edges)\n",
    "#         print(G12.number_of_nodes(), G12.number_of_edges())\n",
    "    \n",
    "        G12 = nx.relabel_nodes(G12, dict(zip(range(len(shared_drugs)), shared_drugs)))\n",
    "        G12 = nx.relabel_nodes(G12, name_pcid_dict)\n",
    "        # write to a .gml file for visualization in Cytoscape\n",
    "        filename = 'Networks/' + fp + '_' + name + '_intersection_network.gml'\n",
    "        nx.write_gml(G12, filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create networkswith just L1000 signature data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "count = 0\n",
    "\n",
    "for fp in fps:\n",
    "\n",
    "    # Load datasets\n",
    "    X1 = fpsfiles_dict[fp] # Dataset1, np.array, e.g. drugs by fingerprints\n",
    "    X2 = namefiles_dict['L1000_sig_new']\n",
    "    X2= X2.T # Dataset2, np.array, e.g. drugs by genes\n",
    "    # # make sure the rows of X1 and X2 are the drugs in the same order\n",
    "#         drugs = [...] # a list of drug IDs/names in X1 and X2\n",
    "\n",
    "    shared_drugs = sorted(list(set(X1.index) & set(X2.index)))\n",
    "    for i in shared_drugs:\n",
    "        if i not in list(name_pcid_dict.keys()):\n",
    "            shared_drugs.remove(i)\n",
    "    X1 = X1.loc[shared_drugs]\n",
    "    X2 = X2.loc[shared_drugs]\n",
    "\n",
    "    # distance metric\n",
    "    A1 = compute_adjcency_mat(X1,'manhattan')\n",
    "    A2 = compute_adjcency_mat(X2, 'cosine')\n",
    "    # Convert to Graphs\n",
    "    percentile = 90 # try adjust this param \n",
    "    G1 = create_graph_by_threshold(A1, percentile)\n",
    "#   print(G1.number_of_nodes(), G1.number_of_edges())\n",
    "\n",
    "    G2 = create_graph_by_threshold(A2, percentile)\n",
    "#   print(G2.number_of_nodes(), G2.number_of_edges())\n",
    "\n",
    "    G12 = nx.intersection(G1, G2)\n",
    "    nodes_wo_edges = [n for n, k in G12.degree() if k == 0]\n",
    "    G12.remove_nodes_from(nodes_wo_edges)\n",
    "    print(G12.number_of_nodes(), G12.number_of_edges())\n",
    "\n",
    "\n",
    "    # label the nodes with drug names/ID\n",
    "    G12 = nx.relabel_nodes(G12, dict(zip(range(len(shared_drugs)), shared_drugs)))\n",
    "    G12 = nx.relabel_nodes(G12, name_pcid_dict)\n",
    "    # write to a .gml file for visualization in Cytoscape\n",
    "    filename = 'Networks/' + fp + '_' + 'L1000_sig_new' + '_intersection_network.gml'\n",
    "    nx.write_gml(G12, filename)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}