Average_distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import itertools\n",
    "from itertools import combinations_with_replacement\n",
    "import networkx as nx\n",
    "from itertools import product\n",
    "\n",
    "\n",
    "def read_gene_file(gene_file):\n",
    "    \"\"\"\n",
    "    Reads a list genes from an external file.\n",
    "\n",
    "    * The genes must be provided as a table in rows. The first column will\n",
    "    be the identifier, the rest are the values.\n",
    "\n",
    "    * Lines that start with '#' will be ignored\n",
    "    \"\"\"\n",
    "    no_genes=0\n",
    "    genes_set ={}\n",
    "    fp=open(gene_file, 'r', encoding='utf-8', errors='ignore')\n",
    "    for line in fp:\n",
    "    \n",
    "        # lines starting with '#' will be ignored\n",
    "        if line[0]=='#':\n",
    "            continue\n",
    "        # the first column in the line will be interpreted as a seed\n",
    "        # gene:\n",
    "        line_data = line.strip().split('\\t')\n",
    "        disease      = line_data[0]\n",
    "        genes      = set(line_data[1:])\n",
    "        no_genes+=len(genes)\n",
    "        genes_set[disease]=genes\n",
    "    fp.close()\n",
    "\n",
    "    return genes_set \n",
    "\n",
    "\n",
    "# funnction to give you the averate shortest dist between two sets of nodes \n",
    "def av_short_dist_comb(G,a, b):\n",
    "    set1=set(G.nodes())&set(a)\n",
    "    set2=set(G.nodes())&set(b)\n",
    "    c = list(itertools.product(set1, set2))\n",
    "    sp=[]\n",
    "    for what in c:\n",
    "        one_path=nx.shortest_path_length(G, what[0],  what[1])\n",
    "        sp.append(one_path)\n",
    "    totsp = sum(sp)\n",
    "    average_shortest_dist= totsp / len(sp)\n",
    "    return average_shortest_dist\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'mydata/mono_AIM_AIF_IDs.txt'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m--------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                        Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-47b97208ce36>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m#this will be your dictionary with the pathways and their assocaited nodes\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mgenelist_A\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mread_gene_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'mydata/mono_AIM_AIF_IDs.txt'\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      7\u001b[0m \u001b[0mgenelist_B\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mread_gene_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'mydata/mono_AIM_AIF_IDs.txt'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-1-a746922a5292>\u001b[0m in \u001b[0;36mread_gene_file\u001b[0;34m(gene_file)\u001b[0m\n\u001b[1;32m     18\u001b[0m     \u001b[0mno_genes\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0mgenes_set\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---> 20\u001b[0;31m     \u001b[0mfp\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgene_file\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mencoding\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'utf-8'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'ignore'\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     21\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mfp\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'mydata/mono_AIM_AIF_IDs.txt'"
     ]
    }
   ],
   "source": [
    "# calling the function on your favorite dictionaries with a set of nodes\n",
    "\n",
    "\n",
    "\n",
    "#this will be your dictionary with the pathways and their assocaited nodes \n",
    "genelist_A=read_gene_file('mydata/mono_AIM_AIF_IDs.txt')\n",
    "genelist_B=read_gene_file('mydata/mono_AIM_AIF_IDs.txt')\n",
    "\n",
    "# this is the network that you are using\n",
    "G0=nx.read_edgelist('mydata/Supplementary_File_1_interactome_edgelist_cleaned.txt')\n",
    "#taking the lcc of the network (to make sure that there are always shortest paths available)\n",
    "G = (G0.subgraph(c) for c in nx.connected_components(G0))\n",
    "G = list(G)[0] # the lcc of the network that you are using\n",
    "\n",
    "\n",
    "# distance calculations starting\n",
    "t0 = time.time()\n",
    "mycombinations= list(itertools.product(set(genelist_A.keys()), set(genelist_B.keys())))\n",
    "\n",
    "\n",
    "distance_dict={} # this dicitonary will store all your distance pairs\n",
    "for pathwaytubple in mycombinations:\n",
    "    genes_A = set(genelist_A[pathwaytubple[0]]) & set(G.nodes())\n",
    "    genes_B = set(genelist_B[pathwaytubple[1]]) & set(G.nodes())\n",
    "    if (len(genes_A)>0) &(len(genes_B)>0):\n",
    "        mydist=av_short_dist_comb(G,genes_A,genes_B)\n",
    "        mykey = pathwaytubple[0], pathwaytubple[1]\n",
    "        distance_dict[mykey]=mydist            \n",
    "print('computing time - group dist calculation: %.2f (sec)' %(float(time.time()-t0)))\n",
    "\n",
    "#print(distance_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# significance calculation\n",
    "#\n",
    "\n",
    "import time \n",
    "t0 = time.time()\n",
    "print('randomisation')\n",
    "import random as rand\n",
    "import statistics as stat\n",
    "\n",
    "rand_res_dict={}\n",
    "usage_of_dict=0\n",
    "number=1000\n",
    "rand_resu={}\n",
    "for groupa, genesa in genelist_A.items():\n",
    "    for groupb, genesb in genelist_B.items():\n",
    "        genes_A = set(genesa) & set(G.nodes())\n",
    "        genes_B = set(genesb) & set(G.nodes())\n",
    "        mykey = groupa, groupb\n",
    "        if (len(genes_A)>0)  &(len(genes_B)>0):\n",
    "            rand_list=[]\n",
    "            rand_res_key=len(genes_A), len(genes_B)\n",
    "            #print(rand_res_key)\n",
    "            if rand_res_key in rand_res_dict.keys():\n",
    "                rand_list=rand_res_dict[rand_res_key]\n",
    "                usage_of_dict+=1\n",
    "                print(usage_of_dict)\n",
    "            else:\n",
    "                while len(rand_list)<number:\n",
    "\n",
    "                    random_nodes_1=set((rand.sample(G.nodes(),len(genes_A))))\n",
    "                    random_nodes_2=set((rand.sample(G.nodes(),len(genes_B))))\n",
    "                    d_AB_rand = distance_min(G,random_nodes_1,random_nodes_2)\n",
    "                    rand_list.append(d_AB_rand)\n",
    "                print('done randomisation for %s, %s' %(groupa,groupb ))\n",
    "                rand_res_dict[rand_res_key]=rand_list\n",
    "            rand_resu[mykey]=rand_list\n",
    "            \n",
    "\n",
    "print('computing time - min distance significance with %sX permut and rand res dict %.2f (sec)' %( number, float(time.time()-t0)))\n",
    "print('dict used: %s times' %(usage_of_dict))\n",
    "\n",
    "zscores={}\n",
    "for my_key, results in rand_resu.items():\n",
    "    nums = list(map(int, results))\n",
    "    try:\n",
    "        devi=stat.stdev(nums)\n",
    "        if devi==0:\n",
    "            devi=0.000000000001\n",
    "    except:\n",
    "        print(\"stdev not possible\")\n",
    "\n",
    "    my_Z=(min_distance_dict[my_key]-(stat.mean(nums)))/(devi)\n",
    "    zscores[my_key]=my_Z"
   ]
  }
 ],
 "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
}