examples/preload-labels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preload Labels for 3KB Windows\n",
    "\n",
    "In this notebook we're using existing peak calls for **Encode e11.5's face and hindbrain** dataset and search for differentially expressed peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import bbi\n",
    "import json\n",
    "import numpy as np\n",
    "import os\n",
    "import sqlite3\n",
    "import sys\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Configurations**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "MAX_PRELOADED_LABELS = 100\n",
    "NUM_SEARCHES_TO_BE_PRELOADED = 10\n",
    "CLEAR_DB = True\n",
    "\n",
    "###########################################\n",
    "# Only change if you know what you're doing\n",
    "###########################################\n",
    "\n",
    "base = \"../\"\n",
    "settings_filepath = \"config-user-study-encode-e11-5-face-hindbrain.json\"\n",
    "window_size = 3000\n",
    "step_size = 1500\n",
    "resolution = 25\n",
    "# 1395142003 is the absolute offset of chr10\n",
    "target_from = 1395142003 + 57039000\n",
    "target_to = 1395142003 + 57042000\n",
    "\n",
    "assert target_to - target_from == window_size\n",
    "\n",
    "# Minimum value to consider a peak annotation a peak for differential accessible peak annotations\n",
    "min_peak_val_diff = 0.75 \n",
    "# Minimum value to consider a peak annotation a peak for equally accessible peak annotations\n",
    "min_peak_val_same = 1\n",
    "\n",
    "with open(os.path.join(base, settings_filepath), \"r\") as f:\n",
    "    settings = json.load(f)\n",
    "\n",
    "signal_face = \"data/ENCFF373NJX.bigWig\"\n",
    "signal_hindbrain = \"data/ENCFF943PHW.bigWig\"\n",
    "\n",
    "narrow_peaks_face = \"data/ENCFF545ITR.bigBed\"\n",
    "narrow_peaks_hindbrain = \"data/ENCFF007GMX.bigBed\"\n",
    "\n",
    "broad_peaks_face = \"data/ENCFF285BLZ.bigBed\"\n",
    "broad_peaks_hindbrain = \"data/ENCFF007GMX.bigBed\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ignore warnings as they just pollute the output\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# Enable importing modules from the parent directory\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "module_path = os.path.abspath(os.path.join('../experiments'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)\n",
    "module_path = os.path.abspath(os.path.join('../server'))\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract windows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Extracted 87129 windows from chr10 with a max value of 1.0.\n",
      "Extracted 87129 windows from chr10 with a max value of 1.0.\n"
     ]
    }
   ],
   "source": [
    "from server.bigwig import chunk\n",
    "\n",
    "windows_face = chunk(\n",
    "    signal_face,\n",
    "    window_size,\n",
    "    resolution,\n",
    "    window_size // settings['step_freq'],\n",
    "    settings['chroms'],\n",
    "    verbose=True,\n",
    ")\n",
    "\n",
    "windows_hindbrain = chunk(\n",
    "    signal_hindbrain,\n",
    "    window_size,\n",
    "    resolution,\n",
    "    window_size // settings['step_freq'],\n",
    "    settings['chroms'],\n",
    "    verbose=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Get the max signal per window**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_signal_face = np.max(windows_face, axis=1)\n",
    "max_signal_hindbrain = np.max(windows_hindbrain, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find differentially accessible peaks much faster\n",
    "\n",
    "`chunk_beds_binary()` extracts only a binary value per window: `1` if a window contains an annotation, i.e., a peak, or `0` if not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chunk_beds(bigbed):\n",
    "    bins = 11\n",
    "\n",
    "    chrom_sizes = bbi.chromsizes(bigbed)\n",
    "    chrom_size = chrom_sizes[settings['chroms'][0]]\n",
    "    num_total_windows = np.ceil((chrom_size - window_size) / step_size).astype(int) + 1\n",
    "\n",
    "    num_windows = np.ceil((chrom_size - window_size) / step_size).astype(int) + 1\n",
    "    start_pos = np.arange(0, step_size * num_total_windows, step_size)\n",
    "    end_pos = np.arange(window_size, step_size * num_total_windows + window_size, step_size)\n",
    "\n",
    "    return bbi.stackup(\n",
    "        bigbed,\n",
    "        settings['chroms'] * num_total_windows,\n",
    "        start_pos,\n",
    "        end_pos,\n",
    "        bins=bins,\n",
    "        missing=0,\n",
    "        oob=0,\n",
    "    ).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks_face = chunk_beds(narrow_peaks_face)\n",
    "peaks_hindbrain = chunk_beds(narrow_peaks_hindbrain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Face peaks: 11008\n",
      "Hindbrain peaks: 8347\n",
      "Diff peaks: 1762\n",
      "Same peaks: 5130\n"
     ]
    }
   ],
   "source": [
    "print('Face peaks: {}'.format(np.sum(np.max(peaks_face[:,2:9], axis=1))))\n",
    "print('Hindbrain peaks: {}'.format(np.sum(np.max(peaks_hindbrain[:,2:9], axis=1))))\n",
    "\n",
    "diff_peaks = (\n",
    "    (\n",
    "        np.max(peaks_face[:,2:9], axis=1) + np.max(peaks_hindbrain[:,2:9], axis=1) == 1\n",
    "    ) & (\n",
    "        np.abs(np.sum(peaks_face[:,2:9], axis=1) - np.sum(peaks_hindbrain[:,2:9], axis=1)) > 2\n",
    "    )\n",
    ")\n",
    "print('Diff peaks: {}'.format(np.sum(diff_peaks)))\n",
    "\n",
    "same_peaks = (\n",
    "    np.max(peaks_face[:,2:9], axis=1) + np.max(peaks_hindbrain[:,2:9], axis=1)\n",
    ") == 2\n",
    "print('Same peaks: {}'.format(np.sum(same_peaks)))\n",
    "\n",
    "diff_peaks_win_ids = np.where(diff_peaks)[0]\n",
    "same_peaks_win_ids = np.where(same_peaks)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diff peaks with max val >= 0.75: 55\n",
      "Same peaks with max val >= 1: 1201\n"
     ]
    }
   ],
   "source": [
    "diff_peaks_with_max = diff_peaks & ((max_signal_face >= min_peak_val_diff) | (max_signal_hindbrain >= min_peak_val_diff))\n",
    "diff_peaks_with_max_ids = np.where(diff_peaks_with_max)[0]\n",
    "\n",
    "print('Diff peaks with max val >= {}: {}'.format(min_peak_val_diff, np.sum(diff_peaks_with_max)))\n",
    "\n",
    "same_peaks_with_max = same_peaks & ((max_signal_face >= min_peak_val_same) | (max_signal_hindbrain >= min_peak_val_same))\n",
    "same_peaks_with_max_ids = np.where(same_peaks_with_max)[0]\n",
    "\n",
    "print('Same peaks with max val >= {}: {}'.format(min_peak_val_same, np.sum(same_peaks_with_max)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preload Search DB with some Labels\n",
    "\n",
    "Preload at most `MAX_PRELOADED_LABELS` positive and negative differentially accessible peaks. We are limiting the number to not overrepresent negative examples as there seem to be many more peaks that are equally accessible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "db6f614cd1774f999e54ba1ad336dc01",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Checkbox(value=False, description='Clear DB (Make sure you know what you do!)')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets.widgets import Checkbox\n",
    "\n",
    "clear_db = Checkbox(value=False, description='Clear DB (Make sure you know what you do!)')\n",
    "clear_db"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "from server.config import Config\n",
    "from server.database import DB\n",
    "\n",
    "db_path = os.path.join(base, settings[\"db_path\"])\n",
    "\n",
    "if os.path.exists(db_path) and not clear_db.value:\n",
    "    print('Database already exist. Check above to delete!')\n",
    "else:\n",
    "    os.remove(db_path)\n",
    "    DB(db_path=db_path, clear=True)\n",
    "\n",
    "    with sqlite3.connect(db_path) as db:\n",
    "        for search_id in range(1, NUM_SEARCHES_TO_BE_PRELOADED + 1):\n",
    "            db.execute(\n",
    "                \"\"\"\n",
    "                    INSERT INTO\n",
    "                        search(id, target_from, target_to, config)\n",
    "                    VALUES\n",
    "                        (?, ?, ?, ?);\n",
    "                \"\"\",\n",
    "                (int(search_id), int(target_from), int(target_to), json.dumps(settings)),\n",
    "            )\n",
    "\n",
    "            for window_idx in np.random.choice(\n",
    "                diff_peaks_with_max_ids,\n",
    "                np.min((diff_peaks_with_max_ids.size, MAX_PRELOADED_LABELS)),\n",
    "                replace=False\n",
    "            ):\n",
    "                db.execute(\n",
    "                    \"\"\"\n",
    "                        INSERT INTO\n",
    "                            classification(search_id, window_id, is_positive)\n",
    "                        VALUES\n",
    "                            (?, ?, ?);\n",
    "                    \"\"\",\n",
    "                    (int(search_id), int(window_idx), 1),\n",
    "                )\n",
    "\n",
    "\n",
    "            for window_idx in np.random.choice(\n",
    "                same_peaks_with_max_ids,\n",
    "                np.min((same_peaks_with_max_ids.size, MAX_PRELOADED_LABELS)),\n",
    "                replace=False\n",
    "            ):\n",
    "                db.execute(\n",
    "                    \"\"\"\n",
    "                        INSERT INTO\n",
    "                            classification(search_id, window_id, is_positive)\n",
    "                        VALUES\n",
    "                            (?, ?, ?);\n",
    "                    \"\"\",\n",
    "                    (int(search_id), int(window_idx), -1),\n",
    "                )\n",
    "\n",
    "            db.commit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Make sure to start the server first!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import time\n",
    "\n",
    "for search_id in range(NUM_SEARCHES_TO_BE_PRELOADED, 0, -1):\n",
    "    r = requests.post(\n",
    "        url = f'http://localhost:5000/api/v1/classifier/?s={search_id}'\n",
    "    )\n",
    "    time.sleep(5)\n",
    "    r = requests.post(\n",
    "        url = f'http://localhost:5000/api/v1/progress/?s={search_id}&u=1'\n",
    "    )\n",
    "    time.sleep(5)"
   ]
  },
  {
   "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": 4
}