14.rwr-proximities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute gene-pair proximities using random walk with restart\n",
    "\n",
    "Export a table where the most proximal genes for each gene are exported.\n",
    "\n",
    "Opened [this issue](https://github.com/sknetwork-team/scikit-network/issues/347) to see whether a more efficient implementation exists.\n",
    "\n",
    "The resulting dataset `combined-score-top-proximities.tsv.xz` can be large. Users might want to filter the dataset using a command-line script for convenience, like:\n",
    "\n",
    "```shell\n",
    "# filter combined-score-top-proximities.tsv.xz with max_proximal_genes set to 50\n",
    "curl --location --silent \\\n",
    "  https://github.com/related-sciences/string-protein-network/raw/2ac8e35e7983fa2d762142e7f287ade3bc8cfb63/data/proximities/combined-score-top-proximities.tsv.xz \\\n",
    "  | awk 'NR==1 || ($6 <= 50)' \\\n",
    "  | gzip --no-name >| combined-score-top-proximities-filtered.tsv.gz\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.sparse\n",
    "import sknetwork as skn\n",
    "import tqdm.notebook\n",
    "\n",
    "from utils import get_protein_info_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "max_proximal_genes = 100\n",
    "proximity_threshold = 1.0\n",
    "path = \"data/proximities/combined-score-top-proximities.tsv.xz\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>protein_external_id</th>\n",
       "      <th>preferred_name</th>\n",
       "      <th>protein_size</th>\n",
       "      <th>annotation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>9606.ENSP00000000233</td>\n",
       "      <td>ARF5</td>\n",
       "      <td>180</td>\n",
       "      <td>ADP-ribosylation factor 5; GTP-binding protein...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>9606.ENSP00000000412</td>\n",
       "      <td>M6PR</td>\n",
       "      <td>277</td>\n",
       "      <td>Cation-dependent mannose-6-phosphate receptor;...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   index   protein_external_id preferred_name  protein_size  \\\n",
       "0      0  9606.ENSP00000000233           ARF5           180   \n",
       "1      1  9606.ENSP00000000412           M6PR           277   \n",
       "\n",
       "                                          annotation  \n",
       "0  ADP-ribosylation factor 5; GTP-binding protein...  \n",
       "1  Cation-dependent mannose-6-phosphate receptor;...  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_df = get_protein_info_df()\n",
    "gene_df.head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix = scipy.sparse.load_npz(\"data/score-matrices/combined_score.sparse.npz\")\n",
    "matrix = matrix / 1000  # rescale STRING scores to be from 0 to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_rwr_scores(seeds):\n",
    "    pagerank = skn.ranking.PageRank(damping_factor=0.85)\n",
    "    np.random.seed(0)  # set seed in case there is some non-determinism\n",
    "    scores = pagerank.fit_transform(adjacency=matrix, seeds=seeds)\n",
    "    scores /= scores.mean()  # rescale scores to have a mean of 1\n",
    "    return scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19566"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nrows = matrix.shape[0]\n",
    "nrows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.15830465, 1.25733971, 1.98393592, ..., 1.        , 1.        ,\n",
       "       0.41587283])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "null_scores = get_rwr_scores(seeds=np.ones(nrows))\n",
    "null_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "828f3aced40d4dbaa162b679b6c46c72",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=19566.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "CPU times: user 4h 38min 3s, sys: 32min 37s, total: 5h 10min 40s\n",
      "Wall time: 5h 10min 33s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "iterrows = tqdm.notebook.tqdm(gene_df.iterrows(), total=len(gene_df))\n",
    "kwargs_chunk = dict(mode=\"a\", header=False)\n",
    "for i, row in iterrows:\n",
    "    scores = get_rwr_scores(seeds={i: 1.0})\n",
    "    df = pd.DataFrame(dict(\n",
    "        gene_a=row.preferred_name,\n",
    "        gene_b=gene_df.preferred_name,\n",
    "        score=(1000 * np.asarray(matrix[i].todense()).flatten()).astype(int),\n",
    "        proximity=scores,\n",
    "        proximity_corrected=scores - null_scores,\n",
    "    ))\n",
    "    df = (\n",
    "        df\n",
    "        .query(\"proximity_corrected >= @proximity_threshold and gene_a != gene_b\")\n",
    "        .sort_values([\"proximity_corrected\"], ascending=False)\n",
    "    )\n",
    "    df = df.iloc[:max_proximal_genes, :]\n",
    "    df[\"rank\"] = df.proximity_corrected.rank(method=\"first\", ascending=False).astype(int)\n",
    "    kwargs = kwargs_chunk if i else {}\n",
    "    df.to_csv(path, sep=\"\\t\", index=False, **kwargs, float_format=\"%.3g\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>gene_a</th>\n",
       "      <th>gene_b</th>\n",
       "      <th>score</th>\n",
       "      <th>proximity</th>\n",
       "      <th>proximity_corrected</th>\n",
       "      <th>rank</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1517</th>\n",
       "      <td>OR6Q1</td>\n",
       "      <td>GNGT1</td>\n",
       "      <td>900</td>\n",
       "      <td>141.113393</td>\n",
       "      <td>134.955583</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19077</th>\n",
       "      <td>OR6Q1</td>\n",
       "      <td>GNB1</td>\n",
       "      <td>900</td>\n",
       "      <td>142.026442</td>\n",
       "      <td>134.765930</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7976</th>\n",
       "      <td>OR6Q1</td>\n",
       "      <td>GNAL</td>\n",
       "      <td>900</td>\n",
       "      <td>139.234503</td>\n",
       "      <td>134.187622</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9549</th>\n",
       "      <td>OR6Q1</td>\n",
       "      <td>RTP2</td>\n",
       "      <td>902</td>\n",
       "      <td>136.466782</td>\n",
       "      <td>132.641310</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6176</th>\n",
       "      <td>OR6Q1</td>\n",
       "      <td>RTP1</td>\n",
       "      <td>900</td>\n",
       "      <td>136.308976</td>\n",
       "      <td>132.444731</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      gene_a gene_b  score   proximity  proximity_corrected  rank\n",
       "1517   OR6Q1  GNGT1    900  141.113393           134.955583     1\n",
       "19077  OR6Q1   GNB1    900  142.026442           134.765930     2\n",
       "7976   OR6Q1   GNAL    900  139.234503           134.187622     3\n",
       "9549   OR6Q1   RTP2    902  136.466782           132.641310     4\n",
       "6176   OR6Q1   RTP1    900  136.308976           132.444731     5"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head(5)"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}