point-data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Labeled Point Data\n",
    "\n",
    "In this notebook we're going to implement a custom tile generator to serve custom data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HOST = 'localhost'\n",
    "PORT = 8002\n",
    "\n",
    "# If using Binder, set this to the current host: https://hub.{domain}/user-{username}-{hash}\n",
    "LOCATION = ''\n",
    "\n",
    "if LOCATION:\n",
    "    # Use the jupyterhub proxy address to reach our server\n",
    "    API_URL = LOCATION.replace('/lab', '') + f'/proxy/{PORT}/api/v1'\n",
    "else:\n",
    "    API_URL = f'http://{HOST}:{PORT}/api/v1'\n",
    "        \n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data\n",
    "\n",
    "In this example we show a dataset on image-based spatial transcriptomics by [Wang et al., 2018.](https://doi.org/10.1038/s41598-018-22297-7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('data/spatial_gene_expression.csv.gz')\n",
    "\n",
    "print('Number of points: {}'.format(len(df)))\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Custom tileset API\n",
    "\n",
    "To generate and serve custom tiles we need a `tileset_info()` and a `tiles()` function as demonstrated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from higlass import Tileset\n",
    "from clodius.tiles.format import format_dense_tile\n",
    "from clodius.tiles.utils import tile_bounds\n",
    "\n",
    "\n",
    "def dfdensity(df, max_zoom=5, **kwargs):\n",
    "    \"\"\"\n",
    "    Create an \"on-the-fly\" point density tileset from a dataframe of (x, y) points.\n",
    "    To be rendered as a heatmap track.\n",
    "    \n",
    "    \"\"\"\n",
    "    data = df.reindex(columns=['x', 'y']).values\n",
    "\n",
    "    tsinfo = {\n",
    "        'tile_size': 256,\n",
    "        'min_pos': [df['x'].min(), df['y'].min()],\n",
    "        'max_pos': [df['x'].max(), df['y'].max()],\n",
    "        'max_width': max(\n",
    "            df['x'].max() - df['x'].min(),\n",
    "            df['y'].max() - df['y'].min()\n",
    "        ),\n",
    "        'max_zoom': max_zoom,\n",
    "        'mirror_tiles': 'false'\n",
    "    }\n",
    "    \n",
    "    def tileset_info():\n",
    "        return tsinfo\n",
    "            \n",
    "    def _get_tile(z, x, y):\n",
    "        extent = tile_bounds(tsinfo, z, x, y)\n",
    "        \n",
    "        # get all the points within the extent\n",
    "        points = data[\n",
    "            (data[:, 0] > extent[0]) &\n",
    "            (data[:, 0] < extent[2]) &\n",
    "            (data[:, 1] > extent[1]) &\n",
    "            (data[:, 1] < extent[3])\n",
    "        ]\n",
    "        \n",
    "        # Generate a 2D histogram\n",
    "        hist, _, _ = np.histogram2d(\n",
    "            points[:, 0],\n",
    "            points[:, 1],\n",
    "            bins=256\n",
    "        )\n",
    "        \n",
    "        # Set empty bins to `nan` to make them transparent\n",
    "        hist[hist == 0.] = np.nan\n",
    "        \n",
    "        return hist.T\n",
    "    \n",
    "    def tiles(tile_ids):\n",
    "        tiles = []\n",
    "        \n",
    "        for tile_id in tile_ids:\n",
    "            # decompose the tile zoom and location\n",
    "            _, z, x, y = tile_id.split('.')\n",
    "            \n",
    "            # generate the tile\n",
    "            data = _get_tile(int(z), int(x), int(y))\n",
    "            \n",
    "            # format the tile response\n",
    "            tiles.append((tile_id, format_dense_tile(data)))\n",
    "    \n",
    "        return tiles\n",
    "\n",
    "    return Tileset(\n",
    "        tileset_info=tileset_info,\n",
    "        tiles=tiles,\n",
    "        **kwargs\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the data as usual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from higlass import Track, View, display\n",
    "\n",
    "all_genes_tileset = dfdensity(df, name='Wang et al., 2018. All Genes')\n",
    "\n",
    "track_config = {\n",
    "    'track_type': 'heatmap',\n",
    "    'position': 'center',\n",
    "    'height': 600,\n",
    "    'options': {\n",
    "        'colorRange': ['#ffbb33', '#e5001c', 'black'],\n",
    "        'backgroundColor': 'white',\n",
    "    },\n",
    "    'server': API_URL,\n",
    "}\n",
    "\n",
    "widget, server, viewconf = display([\n",
    "        View([Track(tileset=all_genes_tileset, **track_config)]),\n",
    "    ],\n",
    "    host=HOST, \n",
    "    server_port=PORT\n",
    ")\n",
    "\n",
    "widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we implemented a custom tileset API we can add any kind of customization like picking which channels (i.e., genes) we actually want to visualize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Important change ##########\n",
    "tilesets = {\n",
    "    'all_genes': dfdensity(\n",
    "        df, \n",
    "        name='Wang et al., 2018. All Genes'\n",
    "    ),\n",
    "    'eno1': dfdensity(\n",
    "        df[df['gene'] == 'ENO1'], \n",
    "        name='Wang et al., 2018. ENO1 Only'\n",
    "    ),\n",
    "    'notch2': dfdensity(\n",
    "        df[df['gene'] == 'NOTCH2'], \n",
    "        name='Wang et al., 2018. NOTCH2 Only'\n",
    "    ),\n",
    "}\n",
    "\n",
    "tracks = {\n",
    "    gene: Track(tileset=ts, **track_config) for gene, ts in tilesets.items()\n",
    "}\n",
    "##############################\n",
    "\n",
    "view_options = {\n",
    "    'y': 0,\n",
    "    'width': 4,\n",
    "}\n",
    "\n",
    "widget, server, viewconf = display([\n",
    "    View(x=0, tracks=[tracks['all_genes']], **view_options),\n",
    "    View(x=4, tracks=[tracks['eno1']], **view_options),\n",
    "    View(x=8, tracks=[tracks['notch2']], **view_options)],\n",
    "    host=HOST, \n",
    "    server_port=PORT,\n",
    ")\n",
    "\n",
    "widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plasma vs nuclear transcribed RNAs**\n",
    "\n",
    "In the following closeup we can see that ENO1 is primarily transcriped in the cell body (outside the nucleus) while NOTCH2 is transcribed within the nucleus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "\n",
    "closeup_view_options = deepcopy(view_options)\n",
    "closeup_view_options['initialXDomain'] = [-333, -209]\n",
    "closeup_view_options['initialYDomain'] = [-202, -30]\n",
    "\n",
    "closeup_tracks = {\n",
    "    gene: track.change_options(maxZoom=2) for gene, track in tracks.items()\n",
    "}\n",
    "\n",
    "widget, server, viewconf = display([\n",
    "    View(x=0, tracks=[closeup_tracks['all_genes']], **closeup_view_options),\n",
    "    View(x=4, tracks=[closeup_tracks['eno1']], **closeup_view_options),\n",
    "    View(x=8, tracks=[closeup_tracks['notch2']], **closeup_view_options)],\n",
    "    host=HOST,\n",
    "    server_port=PORT\n",
    ")\n",
    "\n",
    "widget"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tileset Response from Server\n",
    "\n",
    "Most HiGlass tracks expect the data to be base64 encoded numpy arrays (`dense`) with a the data type (`dtype`), length `size`, and precalculated min and max values (`min_value` and `max_value`) for efficient value scaling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "response = server.tiles(tilesets['all_genes'].uuid, 0, 0, 0)\n",
    "print(\n",
    "    f'dtype: {response[\"dtype\"]}', \n",
    "    f'tile_size: {response[\"size\"]}',\n",
    "    f'min_value: {response[\"min_value\"]}' ,\n",
    "    f'max_value: {response[\"max_value\"]}' ,\n",
    "    f'dense: {response[\"dense\"][:20]}...',\n",
    "    sep='\\n'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View config\n",
    "\n",
    "To ultimately visualize the data, HiGlass' viewer requires a view configuration telling HiGlass the what, how, and where."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewconf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "☝️ This dictionary ultimately tells the JavaScript library what dataset to visualize where and how."
   ]
  },
  {
   "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
}