notebooks/rastermap_interactive.ipynb

{
  "cells": [
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "ZdQZj3mhDJ1o"
      },
      "source": [
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/MouseLand/rastermap/blob/main/notebooks/rastermap_interactive.ipynb)\n",
        "\n",
        "# Rastermap sorting of 34k neurons\n",
        "\n",
        "We will use a spontaneous activity recording from [Syeda et al, 2023](https://www.biorxiv.org/content/10.1101/2022.11.03.515121v1.abstract). We recorded 34,086 neurons from mouse sensorimotor cortex for 2+ hours using two-photon calcium imaging at a rate of 3.2Hz. FYI to make the download of the dataset faster, we are analyzing only the first half of the recording. During the recording, the mouse was free to run on an air floating ball, and we recorded the mouse face with a camera at a rate of 50Hz and tracked keypoints on the mouse face.\n",
        "\n",
        "This notebook includes an **interactive** plotting section to explore the spatial relationships among neurons in the dataset."
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "aLdd4ERXPIm9"
      },
      "source": [
        "If you are using colab, to ensure a smooth and efficient experience while running the interactive plot, it's recommended to clear your browser's cache before executing the code. Here's how you can do it for different browsers:\n",
        "\n",
        "**Chrome**: Settings > Privacy and security > Clear browsing data > Check \"Cached images and files\" > Click \"clear data\"\n",
        "\n",
        "**Safari**: Settings > Advanced > check \"Show Develop menu in menu bar\" > Go back to the menu bar > Develop > Empty Caches\n",
        "\n",
        "**Firefox**: Settings > Privacy & security > Clear data > check \"Cached Web Content\" > Clear"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "V-y9iv8rDJ1q"
      },
      "source": [
        "First we will install the required packages, if not already installed. If on google colab, it will require you to click the \"RESTART RUNTIME\" button because we are updating numpy. Also, select the GPU runtime to make the interactive plotting faster:\n",
        "**Runtime > Change runtime type > Hardware accelerator = GPU**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nlFEhwGCDJ1r"
      },
      "outputs": [],
      "source": [
        "!pip install numpy>=1.24 # (required for google colab)\n",
        "!pip install rastermap\n",
        "!pip install matplotlib"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "zh382SYZDJ1r"
      },
      "source": [
        "### Load data and import libraries\n",
        "\n",
        "If not already downloaded, the following cell will automatically download the processed data stored [here](https://osf.io/8xg7n)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Kbm3fPaDDJ1s"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "# importing rastermap\n",
        "# (this will be slow the first time since it is compiling the numba functions)\n",
        "from rastermap import Rastermap, utils\n",
        "from scipy.stats import zscore\n",
        "\n",
        "# download spontaneous activity\n",
        "filename = utils.download_data(data_type=\"spont2\")\n",
        "\n",
        "dat = np.load(filename)\n",
        "\n",
        "# spks is neurons by time\n",
        "# (each timepoint is 313 ms)\n",
        "spks = dat[\"spks\"]\n",
        "n_neurons, n_time = spks.shape\n",
        "print(f\"{n_neurons} neurons by {n_time} timepoints\")\n",
        "\n",
        "# zscore activity (each neuron activity trace is then mean 0 and standard-deviation 1)\n",
        "spks = zscore(spks, axis=1)\n",
        "\n",
        "# XY position of each neuron in the recording\n",
        "xpos, ypos = dat[\"xpos\"], dat[\"ypos\"]\n",
        "\n",
        "# for your own data, you will need \"spks\" and \"xpos\" and \"ypos\""
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "k6ocUQLlDJ1s"
      },
      "source": [
        "### Run Rastermap\n",
        "\n",
        "Let's sort the single neurons with Rastermap, with clustering and upsampling:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JiCXd65nDJ1s"
      },
      "outputs": [],
      "source": [
        "model = Rastermap(n_clusters=100, # number of clusters to compute\n",
        "                  n_PCs=128, # number of PCs to use\n",
        "                  locality=0.75, # locality in sorting to find sequences (this is a value from 0-1)\n",
        "                  time_lag_window=5, # use future timepoints to compute correlation\n",
        "                  grid_upsample=10, # default value, 10 is good for large recordings\n",
        "                ).fit(spks)\n",
        "y = model.embedding # neurons x 1\n",
        "isort = model.isort"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "P5zRt0duDJ1t"
      },
      "source": [
        "Let's create superneurons from Rastermap -- we sort the data and then sum over neighboring neurons:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3BQ46VXgDJ1t"
      },
      "outputs": [],
      "source": [
        "nbin = 200 # number of neurons to bin over\n",
        "sn = utils.bin1d(spks[isort], bin_size=nbin, axis=0) # bin over neuron axis"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "31RN-4zhUAEc"
      },
      "source": [
        "### Interactive Visualization"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "8aaYzoar5cpu"
      },
      "source": [
        "Use the Rastermap sorting to visualize neural activity of all neurons and show the positions of selected neurons. GPU is required for fast rendering."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "m8aeDprVTZ_x"
      },
      "outputs": [],
      "source": [
        "!pip install dash\n",
        "!pip install plotly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "m-jbQW87QRNZ"
      },
      "outputs": [],
      "source": [
        "# @title Interactive Plot - press play and wait around 5sec\n",
        "import dash\n",
        "import plotly.graph_objs as go\n",
        "from dash.dependencies import Input, Output\n",
        "from dash import callback, dcc, html, State\n",
        "\n",
        "matrix = sn\n",
        "NN, NT = matrix.shape\n",
        "\n",
        "# get indices of neurons corresponding to each superneuron\n",
        "NN_all = spks.shape[0]\n",
        "indices_bin = isort[:NN_all // nbin * nbin].reshape((NN_all // nbin, nbin))\n",
        "\n",
        "# decide how many time points to show per frame\n",
        "NT_show = min(1000, NT-1)\n",
        "nmin, nmax = 0, NN\n",
        "tmin, tmax = 0, NT_show\n",
        "\n",
        "# swap the xpos and ypos\n",
        "xpos_plot = ypos\n",
        "ypos_plot = xpos\n",
        "\n",
        "# visualize first frame\n",
        "fig = go.Figure(\n",
        "    data=[\n",
        "        go.Heatmap(\n",
        "            x=np.arange(tmin, tmax).tolist(),\n",
        "            z=matrix[:, tmin:tmax],\n",
        "            colorscale=\"Greys\",\n",
        "            zmin=0,\n",
        "            zmax=0.8,\n",
        "        )\n",
        "    ]\n",
        ")\n",
        "\n",
        "# initialize the positions of the selecting bar\n",
        "x0 = tmin\n",
        "x1 = tmin+NT_show\n",
        "y0 = nmin+int(NN/10) * 8\n",
        "y1 = nmin+int(NN/10) * 9\n",
        "\n",
        "# visualize neurons with their positions\n",
        "color_values = np.ones(len(xpos_plot)) * 0.1\n",
        "size_values = np.ones(len(xpos_plot)) * 5\n",
        "\n",
        "neuron_fig = go.Figure(\n",
        "    data=[\n",
        "        go.Scattergl(x=xpos_plot, y=ypos_plot, mode='markers',\n",
        "                    marker=dict(\n",
        "                  size=size_values,\n",
        "                  color=color_values,\n",
        "                  colorscale='Purples',\n",
        "                  cmin=0,\n",
        "                  cmax=1,\n",
        "                  )\n",
        "                  )\n",
        "    ]\n",
        ")\n",
        "\n",
        "neuron_fig.update_layout(\n",
        "    width=500,\n",
        "    height=500,\n",
        "    yaxis={\"title\": 'y position'},\n",
        "    xaxis={\"title\": 'x position'},\n",
        "    template='simple_white',\n",
        "    margin=dict(l=10, r=0, t=100, b=0),\n",
        ")\n",
        "\n",
        "# define the dash app layout\n",
        "app = dash.Dash(__name__)\n",
        "app.layout = html.Div(\n",
        "    style={'display': 'flex', 'flex-direction': 'row', \"padding\": \"0\", \"margin\": \"0\"},\n",
        "    children=[\n",
        "        html.Div(\n",
        "            style={\"width\": \"60%\", \"display\": \"flex\", \"flex-direction\": \"column\", \"padding\": \"0\", \"margin\": \"0\"},\n",
        "            children=[\n",
        "                html.Div(\n",
        "                    style={'display': 'flex', 'flex-direction': 'row', \"padding\": \"10\", \"margin\": \"10\"},\n",
        "                    children=[\n",
        "                        html.H2(\"Rastermap\", style={'margin': '0'}),\n",
        "                        html.Div(dcc.Input(id='input-on-submit', type='text', placeholder=\"{}\".format(NT_show)), style={'margin-left': '10px', 'margin-top': '50px'}),\n",
        "                        html.Button('Submit', id='submit-val', n_clicks=0, style={'height': '20px', 'margin-left': '10px', 'margin-top': '50px'}),\n",
        "                        html.Div(id='button-output', children=f'number of time points to show: {NT_show}',\n",
        "                                 style = {'margin-left': '10px', 'margin-top': '50px'})\n",
        "                    ]\n",
        "                ),\n",
        "                dcc.Graph(id=\"matrix-plot\",\n",
        "                          figure=fig,\n",
        "                          config={\n",
        "                              'edits': {\n",
        "                                  'shapePosition': True\n",
        "                              }\n",
        "                          },\n",
        "                          style={'margin-top': '0px'}\n",
        "                          ),\n",
        "                dcc.Slider(0, 1, step=1/NT, id='slider-time',\n",
        "                           marks={i*NT/(NT-NT_show): '{}'.format(int(i*NT)) for i in (np.arange(NT+1, step=int(NT/10))/NT).tolist()})\n",
        "            ]\n",
        "        ),\n",
        "        html.Div(\n",
        "            style={\"width\": \"40%\", \"display\": \"inline-block\", \"padding\": \"0\", \"margin\": \"0\"},\n",
        "            children=[\n",
        "                html.H2(\"Neuron locations\"),\n",
        "                dcc.Graph(id=\"neuron-plot\")\n",
        "            ]\n",
        "        )\n",
        "    ]\n",
        ")\n",
        "\n",
        "# call back for slider to change time points to show\n",
        "@app.callback(\n",
        "    Output(\"matrix-plot\", 'figure'),\n",
        "    Input('slider-time', 'value'),\n",
        "    )\n",
        "def update_output(tvalue):\n",
        "    if tvalue is not None:\n",
        "        tmin = int((NT-NT_show)*tvalue)\n",
        "        tmax = tmin + NT_show\n",
        "    else:\n",
        "        tmin, tmax = 0, NT_show\n",
        "\n",
        "    fig = go.Figure(\n",
        "        data=[\n",
        "            go.Heatmap(\n",
        "                x=np.arange(tmin, tmax).tolist(),\n",
        "                z=matrix[:, tmin:tmax],\n",
        "                colorscale=\"Greys\",\n",
        "                zmin=0,\n",
        "                zmax=0.8,\n",
        "            )\n",
        "        ]\n",
        "    )\n",
        "\n",
        "    fig.add_shape(type=\"rect\",\n",
        "        xref=\"x\", yref=\"y\",\n",
        "        x0=tmin, y0=y0,\n",
        "        x1=tmin+x1-x0, y1=y1,\n",
        "        line=dict(\n",
        "            color=\"grey\",\n",
        "            width=3,\n",
        "        ),\n",
        "        fillcolor=\"grey\",\n",
        "        opacity=0.5,\n",
        "        xanchor=tmin,\n",
        "    )\n",
        "\n",
        "    fig.update_layout(\n",
        "        width=800,\n",
        "        height=500,\n",
        "        yaxis={\"title\": 'Neuron'},\n",
        "        xaxis={\"title\": 'Time'},\n",
        "        margin=dict(l=0, r=0, t=50, b=0),\n",
        "    )\n",
        "    return fig\n",
        "\n",
        "# call back for moving the selecting bar to select neurons\n",
        "@app.callback(\n",
        "    Output(\"neuron-plot\", \"figure\"),\n",
        "    Input(\"matrix-plot\", \"relayoutData\"))\n",
        "def update_matrix_plot(relayout_data):\n",
        "    global x0, y0, x1, y1\n",
        "    color_values = np.ones(len(xpos_plot)) * 0.1\n",
        "    size_values = np.ones(len(xpos_plot)) * 5\n",
        "    if relayout_data is not None:\n",
        "        x0, y0 = int(relayout_data[\"shapes[0].x0\"]), int(relayout_data[\"shapes[0].y0\"])\n",
        "        x1, y1 = int(relayout_data[\"shapes[0].x1\"]), int(relayout_data[\"shapes[0].y1\"])\n",
        "        neuron_range = np.arange(y0, y1)\n",
        "        neuron_range = indices_bin[neuron_range].reshape(-1)\n",
        "        color_values[neuron_range] = 1\n",
        "        size_values[neuron_range] = 5\n",
        "\n",
        "    neuron_fig['data'][0]['marker']['color'] = color_values\n",
        "    neuron_fig['data'][0]['marker']['size'] = size_values\n",
        "\n",
        "    return neuron_fig\n",
        "\n",
        "# call back for updating number of time points to show per frame\n",
        "@callback(\n",
        "    Output('button-output', 'children'),\n",
        "    Output('slider-time', 'marks'),\n",
        "    Input('submit-val', 'n_clicks'),\n",
        "    State('input-on-submit', 'value')\n",
        ")\n",
        "def update_timepoints_show(n_clicks, value):\n",
        "    global NT_show, x1\n",
        "    if value is not None:\n",
        "        NT_show = int(value)\n",
        "        x1 = x0+NT_show\n",
        "    else:\n",
        "        value = NT_show\n",
        "    text_to_show = f'number of time points to show: {value}'\n",
        "    new_marks = {i*NT/(NT-NT_show): '{}'.format(int(i*NT)) for i in (np.arange(NT+1, step=int(NT/10))/NT).tolist()}\n",
        "    return text_to_show, new_marks\n",
        "\n",
        "# run the app\n",
        "app.run_server(jupyter_mode='inline')"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "1TJTYl4XDJ1w"
      },
      "source": [
        "### Settings\n",
        "\n",
        "You can see all the rastermap settings with `Rastermap?`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "792xHut7DJ1w"
      },
      "outputs": [],
      "source": [
        "Rastermap?"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {
        "id": "_L1yX4l1DJ1w"
      },
      "source": [
        "### Outputs\n",
        "\n",
        "All the attributes assigned to the Rastermap `model` are listed with `Rastermap.fit?`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QnyOGi2ADJ1w"
      },
      "outputs": [],
      "source": [
        "Rastermap.fit?"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "gpuType": "T4",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.9.17"
    },
    "orig_nbformat": 4,
    "vscode": {
      "interpreter": {
        "hash": "998540cc2fc2836a46e99cd3ca3c37c375205941b23fd1eb4b203c48f2be758f"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}