AutoCore_definition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################################################################\n",
    "#\n",
    "#\n",
    "#\n",
    "# THIS SCRIPT EXPANDS A GIVEN GENE SET TO DEFINE THE FULL SUBGRAPH THAT CONTAINS\n",
    "#\n",
    "# ALL NODES WITHIN THAT GENE SET BASED ON A RANDOM WALK.\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "#\n",
    "##########################################################################\n",
    "\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import networkx as nx\n",
    "import itertools as it\n",
    "import random as rd\n",
    "import pymysql as mysql\n",
    "import pickle \n",
    "#from fa2 import ForceAtlas2\n",
    "import os.path\n",
    "from networkx.utils import pairwise, not_implemented_for\n",
    "from itertools import combinations, chain\n",
    "from collections import defaultdict\n",
    "import time\n",
    "from sklearn.preprocessing import normalize\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "#########################################################################\n",
    "#\n",
    "#          CALCULATES THE RANDOM WALK OPERATOR \n",
    "#             REQUIRES ADJACENCY MATRIX\n",
    "#\n",
    "#########################################################################\n",
    "\n",
    "def rnd_walk_matrix(A, r, a, num_nodes):\n",
    "\n",
    "    num = 1*num_nodes\n",
    "    n = num_nodes\n",
    "    factor = float((1-a)/n)\n",
    "    E = np.multiply(factor,np.ones([n,n])) \n",
    "    A_tele = np.multiply(a,A) + E  \n",
    "\n",
    "#     print(A_tele)\n",
    "    M = normalize(A_tele, norm='l1', axis=0)                                 # column wise normalized MArkov matrix\n",
    "\n",
    "    del A\n",
    "    del A_tele\n",
    "    del E\n",
    "\n",
    "    U = np.identity(n,dtype=int) \n",
    "    H = (1-r)*M\n",
    "    H1 = np.subtract(U,H)\n",
    "    del U\n",
    "    del M\n",
    "    del H    \n",
    "\n",
    "    W = r*np.linalg.inv(H1)   \n",
    "\n",
    "    return W\n",
    "\n",
    "\n",
    "\n",
    "#########################################################################\n",
    "#\n",
    "#          CONSTRUCTING STEINER TREE \n",
    "#\n",
    "#########################################################################\n",
    "\n",
    "def metric_closure(G, weight='weight'):\n",
    "    \"\"\"  Return the metric closure of a graph.\n",
    "\n",
    "    The metric closure of a graph *G* is the complete graph in which each edge\n",
    "    is weighted by the shortest path distance between the nodes in *G* .\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    G : NetworkX graph\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    NetworkX graph\n",
    "        Metric closure of the graph `G`.\n",
    "\n",
    "    \"\"\"\n",
    "    M = nx.Graph()\n",
    "\n",
    "    Gnodes = set(G)\n",
    "\n",
    "    # check for connected graph while processing first node\n",
    "    all_paths_iter = nx.all_pairs_dijkstra(G, weight=weight)\n",
    "    u, (distance, path) = next(all_paths_iter)\n",
    "    if Gnodes - set(distance):\n",
    "        msg = \"G is not a connected graph. metric_closure is not defined.\"\n",
    "        raise nx.NetworkXError(msg)\n",
    "    Gnodes.remove(u)\n",
    "    for v in Gnodes:\n",
    "        M.add_edge(u, v, distance=distance[v], path=path[v])\n",
    "\n",
    "    # first node done -- now process the rest\n",
    "    for u, (distance, path) in all_paths_iter:\n",
    "        Gnodes.remove(u)\n",
    "        for v in Gnodes:\n",
    "            M.add_edge(u, v, distance=distance[v], path=path[v])\n",
    "\n",
    "    return M\n",
    "\n",
    "def steiner_tree(G, terminal_nodes, weight='weight'):\n",
    "    \"\"\" Return an approximation to the minimum Steiner tree of a graph.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    G : NetworkX graph\n",
    "\n",
    "    terminal_nodes : list\n",
    "         A list of terminal nodes for which minimum steiner tree is\n",
    "         to be found.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    NetworkX graph\n",
    "        Approximation to the minimum steiner tree of `G` induced by\n",
    "        `terminal_nodes` .\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Steiner tree can be approximated by computing the minimum spanning\n",
    "    tree of the subgraph of the metric closure of the graph induced by the\n",
    "    terminal nodes, where the metric closure of *G* is the complete graph in\n",
    "    which each edge is weighted by the shortest path distance between the\n",
    "    nodes in *G* .\n",
    "    This algorithm produces a tree whose weight is within a (2 - (2 / t))\n",
    "    factor of the weight of the optimal Steiner tree where *t* is number of\n",
    "    terminal nodes.\n",
    "\n",
    "    \"\"\"\n",
    "    # M is the subgraph of the metric closure induced by the terminal nodes of\n",
    "    # G.\n",
    "    M = metric_closure(G, weight=weight)\n",
    "    # Use the 'distance' attribute of each edge provided by the metric closure\n",
    "    # graph.\n",
    "    H = M.subgraph(terminal_nodes)\n",
    "    mst_edges = nx.minimum_spanning_edges(H, weight='distance', data=True)\n",
    "    # Create an iterator over each edge in each shortest path; repeats are okay\n",
    "    edges = chain.from_iterable(pairwise(d['path']) for u, v, d in mst_edges)\n",
    "    T = G.edge_subgraph(edges)\n",
    "    return T\n",
    "\n",
    "\n",
    "\n",
    "#########################################################################\n",
    "#########################################################################\n",
    "\n",
    "def read_network_from_file(network_file):\n",
    "    \"\"\"\n",
    "    Reads a network from an external file.\n",
    "\n",
    "    * The edgelist must be provided as a tab-separated table. The\n",
    "    first two columns of the table will be interpreted as an\n",
    "    interaction gene1 <==> gene2\n",
    "\n",
    "    * Lines that start with '#' will be ignored\n",
    "    \"\"\"\n",
    "\n",
    "    G = nx.Graph()\n",
    "    for line in open(network_file,'r'):\n",
    "        # lines starting with '#' will be ignored\n",
    "        if line[0]=='#':\n",
    "            continue\n",
    "        # The first two columns in the line will be interpreted as an\n",
    "        # interaction gene1 <=> gene2\n",
    "        line_data   = line.strip().split('\\t')\n",
    "        node1 = line_data[0]\n",
    "        node2 = line_data[1]\n",
    "        G.add_edge(node1,node2)\n",
    "\n",
    "    print (\"\\n> done loading network:\")\n",
    "    print (\"> network contains %s nodes and %s links\" %(G.number_of_nodes(),\n",
    "                                                       G.number_of_edges()))\n",
    "    \n",
    "    return G\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "> done loading network:\n",
      "> network contains 18853 nodes and 483037 links\n"
     ]
    }
   ],
   "source": [
    "# LOAD PPI \n",
    "\n",
    "\n",
    "G_ppi= read_network_from_file('mydata/Supplementary_File_1_interactome_edgelist_cleaned.txt')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RANDOM WALK MATRIX OF FULL NETWORK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inversion done\n",
      "rnd walk computing time: 646.21\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# ██████╗  █████╗ ███╗   ██╗██████╗  ██████╗ ███╗   ███╗    ██╗    ██╗ █████╗ ██╗     ██╗  ██╗\n",
    "# ██╔══██╗██╔══██╗████╗  ██║██╔══██╗██╔═══██╗████╗ ████║    ██║    ██║██╔══██╗██║     ██║ ██╔╝\n",
    "# ██████╔╝███████║██╔██╗ ██║██║  ██║██║   ██║██╔████╔██║    ██║ █╗ ██║███████║██║     █████╔╝ \n",
    "# ██╔══██╗██╔══██║██║╚██╗██║██║  ██║██║   ██║██║╚██╔╝██║    ██║███╗██║██╔══██║██║     ██╔═██╗ \n",
    "# ██║  ██║██║  ██║██║ ╚████║██████╔╝╚██████╔╝██║ ╚═╝ ██║    ╚███╔███╔╝██║  ██║███████╗██║  ██╗\n",
    "# ╚═╝  ╚═╝╚═╝  ╚═╝╚═╝  ╚═══╝╚═════╝  ╚═════╝ ╚═╝     ╚═╝     ╚══╝╚══╝ ╚═╝  ╚═╝╚══════╝╚═╝  ╚═╝\n",
    "                                                                                            \n",
    "\n",
    "\n",
    "\n",
    "t0 = time.time()\n",
    "\n",
    "# RANDOM WALK PARAMETERS\n",
    "r = .9\n",
    "alpha = 1.\n",
    "# infile = sys.argv[4]\n",
    "\n",
    "\n",
    "\n",
    "# ###########################################################################\n",
    "#                                                                           #\n",
    "#    CALCULATES MATRIX INVERSION FOR GIVEN PARAMETERS                       #\n",
    "#                                                                           #\n",
    "                                                                            #\n",
    "num_nodes = G_ppi.number_of_nodes()                                             #\n",
    "A = nx.adjacency_matrix(G_ppi,sorted(G_ppi.nodes()))                                #\n",
    "                                                                            #\n",
    "#################################################################           #\n",
    "#                                                                           #\n",
    "#    INVERT MARKOV MATRIX & GENERATE RW MATRIX                              #\n",
    "#                                                                           #\n",
    "###################                                                         #\n",
    "                                                                            #\n",
    "W = rnd_walk_matrix(A, r, alpha, num_nodes)                                 #\n",
    "                                                                            #\n",
    "print('Inversion done')                                                     #\n",
    "#################################################################           #\n",
    "#                                                                           #\n",
    "#    GENERATE DICT FOR NODE-LABELS TO INTEGERS                              #\n",
    "#                                                                           #\n",
    "###################                                                         #\n",
    "                                                                            #\n",
    "d_idx_entz = {}                                                             #\n",
    "cc = 0                                                                      #\n",
    "for entz in sorted(G_ppi.nodes()):                                              #\n",
    "    d_idx_entz[cc] = entz                                                   #\n",
    "    cc += 1                                                                 #\n",
    "d_entz_idx = dict((y,x) for x,y in d_idx_entz.items())                      #\n",
    "print('rnd walk computing time: %.2f' %float(time.time()-t0))               #\n",
    "                                                                            #\n",
    "#############################################################################\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LOAD SEED SETS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('mydata/mono_AIM_AIF_IDs.txt','r')\n",
    "lines = f.readlines()\n",
    "f.close()\n",
    "mono_all = lines[0].split('\\n')[0].split('\\t')[1:-1]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RND WALK FROM SEEDs UNTIL GRAPH IS CONNECTED "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of initial set: 186 seed genes\n",
      "# nodes: 399\n",
      "# edges: 5244\n"
     ]
    }
   ],
   "source": [
    "seedset_all = mono_all\n",
    "\n",
    "# CHECK PPI OVERLAP\n",
    "\n",
    "seedset_ppi = set(seedset_all) & set(G_ppi.nodes())\n",
    "\n",
    "\n",
    "# PERFORM RANDOM WALK FROM SEED GENES\n",
    "nodeset = []\n",
    "for s_ent in seedset_ppi:\n",
    "    nodeset.append(d_entz_idx[s_ent])\n",
    "\n",
    "normprob = len(nodeset)\n",
    "print('Size of initial set: %s seed genes' %normprob)\n",
    "p0 = np.zeros(G_ppi.number_of_nodes())\n",
    "# generate start vector\n",
    "for n in range(len(p0)):\n",
    "    if n in nodeset:\n",
    "        p0[n] = 1.\n",
    "        \n",
    "pinf =  np.array(W.dot(p0))\n",
    "   \n",
    "# DICT WITH NODE ID AND PVIS \n",
    "d_n_p = {}\n",
    "i = 0\n",
    "for x in pinf:\n",
    "#     print(i,x)\n",
    "    d_n_p[d_idx_entz[i]] = x/normprob\n",
    "    i += 1    \n",
    "       \n",
    "# # write file \n",
    "# outfile = 'ranked_genes_pvis.txt'\n",
    "# f_out = open(outfile,'w')\n",
    "# f_out.write('entrez_id,visit probab,cumulative\\n')\n",
    "\n",
    "pcum = 0\n",
    "l_rndwnodes = []\n",
    "cc = 0\n",
    "for node, pvis in sorted(d_n_p.items(), key = lambda x: x[1], reverse = True):\n",
    "    pcum += pvis\n",
    "#     print(node,pvis,pcum)\n",
    "    l_rndwnodes.append(node)\n",
    "    G_test = G_ppi.subgraph(l_rndwnodes)\n",
    "#     print(nx.is_connected(G_test))\n",
    "    \n",
    "    if nx.is_connected(G_test) and cc>len(seedset_ppi):\n",
    "        break\n",
    "#     del G_test\n",
    "\n",
    "    cc += 1\n",
    "print('# nodes: %s' %G_test.number_of_nodes())\n",
    "print('# edges: %s' %G_test.number_of_edges())\n",
    "\n",
    "nx.write_edgelist(G_test,'mydata/AutoCore.elist')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'd_gene_terms_all' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-1320a41a201e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# position for nodes without any functional features\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mposG\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mnode\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md_gene_terms_all\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkeys\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0mn_index\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0md_gene_idx_feat\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnode\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mxpos\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mX_tsne2D\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn_index\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'd_gene_terms_all' is not defined"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}