2.1.TileDetect.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b25e5e4e-ea76-44ad-aad3-24712d6c3313",
   "metadata": {},
   "source": [
    "# Tile Detection\n",
    "\n",
    "As we don't completely know which tiles belong to which Nova-ST chips, we need to detect this using the data itself.\n",
    "\n",
    "We do this by checking a small subset of the reads from the spatial data vs. the subset files we saved from the HDMI analysis earlier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9868271-849d-49e7-a1d5-1bd08f624114",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pickle\n",
    "from collections import defaultdict\n",
    "\n",
    "import pysam\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa1f185b-dd1e-40d6-beae-d38635c5a3ed",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tiles_dir = 'HDMI_Tiles_Data'\n",
    "subsets = defaultdict(dict)\n",
    "for lane in range(1, 5):\n",
    "    for surface in range(1, 3):\n",
    "        for swath in range(1, 7):\n",
    "            with open(f'{tiles_dir}/{lane}_{surface}{swath}_subset_barcodes.pickle', 'rb') as fh:\n",
    "                subsets[f'{lane}_{surface}{swath}'] = pickle.load(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb7ea9be-a89c-437d-ad6e-8650c94af014",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = \"\" # Demultiplexed spatial libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a89543d1-231c-49a1-9fa0-f187400e0f0f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples = [\n",
    "    # Put the names of the samples here, this should match the demultiplexing sheet and will be appended with _R1_001.fastq.gz later\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f94af9e0-449b-47e0-9fe0-4b138aa93890",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples_data = defaultdict(set)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fb01c4c-d14c-4996-aabf-cf646d96916a",
   "metadata": {},
   "source": [
    "We load the first 100k barcodes from each sample to profile this small subset of reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "701b238b-876f-455e-ae4c-1cf3be8eafde",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for sample in samples:\n",
    "    print(f'On {sample}:')\n",
    "    with pysam.FastqFile(f'{data_dir}/{sample}_R1_001.fastq.gz', persist=False) as fh:\n",
    "        for en, entry in enumerate(fh):\n",
    "            seq = entry.sequence[:32]\n",
    "            samples_data[sample].add(seq)\n",
    "            \n",
    "            if (en + 1) % 100000 == 0:\n",
    "                print(f'    Finished {en + 1} reads')\n",
    "                break\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eaffa6c-7418-4db7-9ee5-79382f76685e",
   "metadata": {},
   "source": [
    "Plot the barcode pattern from the spatial library and ensure it looks correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b46a4ab-556a-44c1-9370-25e087785bc1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for sample in samples:\n",
    "    base_pairs = defaultdict(lambda: defaultdict(lambda: 0))\n",
    "    for seq in samples_data[sample]:\n",
    "        for n, base in enumerate(seq):\n",
    "            base_pairs[n][base] += 1\n",
    "        if (en) % 1000000 == 0:\n",
    "            print(f'Finished {en} reads')\n",
    "            break \n",
    "\n",
    "    bp_df = pd.DataFrame(base_pairs, index=['A', 'C', 'G', 'T', 'N'])\n",
    "\n",
    "    bp_df = bp_df.drop('N')\n",
    "    bp_df.columns = [int(x) + 1 for x in bp_df.columns]\n",
    "\n",
    "    bp_df/bp_df.sum()\n",
    "\n",
    "\n",
    "\n",
    "    (bp_df/bp_df.sum()).T.plot(kind='line', figsize=(18, 12))\n",
    "    plt.grid(True, which='major')\n",
    "    plt.title(sample)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "824a8c5d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Helper functions to convert tiles to plotted columns and back again\n",
    "\n",
    "def get_x(lane, surface, swath):\n",
    "    # Return x coordinate of the tile, swaths are grouped, each lane has 2 surfaces, each surface has 6 swaths, add a gap of 1 between every surface\n",
    "    return (lane - 1) * 14 + (surface - 1) * 7 + swath    \n",
    "\n",
    "def reverse_get_x(x):\n",
    "    # Return lane, surface, swath of the tile\n",
    "    # 1 -> 1, 1, 1\n",
    "    # 2 -> 1, 1, 2\n",
    "    # 7 -> None\n",
    "    # 8 -> 1, 2, 1\n",
    "    # 44 -> 4, 1, 2\n",
    "\n",
    "    lane = (x - 1) // 14 + 1\n",
    "    surface = ((x - 1) % 14) // 7 + 1\n",
    "    swath = ((x - 1) % 14) % 7 + 1\n",
    "    return lane, surface, swath\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eafaf9df-c8ea-46da-b800-a341580e1e95",
   "metadata": {},
   "source": [
    "Here, we count how many times we see a barcode from each tile (subset) in the reads we took from the library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54d92390-326d-434e-a949-ffa52760761b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "lanes = 4\n",
    "surfaces = 2\n",
    "swaths = 6\n",
    "tiles = 78\n",
    "counts_dfs = {}\n",
    "for s, data in samples_data.items():\n",
    "    df = pd.DataFrame(\n",
    "        data=np.nan,\n",
    "        columns=list(range((lanes * surfaces * swaths) + lanes + surfaces + 1)),\n",
    "        index=[x for x in range(1, tiles + 1)],\n",
    "    )\n",
    "    col = 0\n",
    "    for lane in range(1, lanes + 1):\n",
    "        for surface in range(1, surfaces + 1):\n",
    "            for swath in range(1, swaths + 1):\n",
    "                for tile, tile_bcs in subsets[f\"{lane}_{surface}{swath}\"].items():\n",
    "                    x = get_x(lane, surface, swath) - 1\n",
    "                    y = int(tile)\n",
    "                    df.loc[y, x] = len(tile_bcs & data)\n",
    "\n",
    "    counts_dfs[s] = df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16dc5cd-2da8-4644-b4aa-eecfebc1b28f",
   "metadata": {},
   "source": [
    "Higlight every tile where one or more reads was identified to visually inspect where the tissue is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0b934b1-dad5-4c73-af47-f60e6c2f7d64",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, len(samples) ,figsize=((20 * len(samples)), 78), facecolor='white')\n",
    "current_cmap = matplotlib.cm.get_cmap()\n",
    "current_cmap.set_bad(color='black')\n",
    "for ax, sample, df in zip(axs, counts_dfs.keys(), counts_dfs.values()):\n",
    "    ax.imshow(df, vmax=1, vmin=0, )\n",
    "    ax.set_title(sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2832b44c-bf40-449a-b2dc-e64742884681",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def extract_whitelist(start, end, path, lane=1, surface=1, swaths=[1, 2, 3, 4, 5, 6], trim_to_31=True, sample=None):\n",
    "    final_path = f\"{path}/whitelist\"\n",
    "    if sample:\n",
    "        final_path += f\"_{sample}\"\n",
    "    final_path += f\"_{lane}_{surface}_tiles_{start}_to_{end}\"\n",
    "    if trim_to_31:\n",
    "        final_path += \"_31bp\"\n",
    "    final_path += \".tsv\"\n",
    "\n",
    "    with open(final_path, \"w\") as whitelist_fh:\n",
    "        for tile_y in range(start, end + 1):\n",
    "            for swath in swaths:\n",
    "                tile = f\"{surface}{swath}{tile_y}\"\n",
    "                with open(f\"{tiles_dir}/{lane}_{tile}_barcodes.pickle\", \"rb\") as tile_fh:\n",
    "                    tile_data = pickle.load(tile_fh)\n",
    "                for bc in tile_data.keys():\n",
    "                    if trim_to_31:\n",
    "                        whitelist_fh.write(f\"{bc[:31]}\\n\")\n",
    "                    else:\n",
    "                        whitelist_fh.write(f\"{bc}\\n\")\n",
    "        return final_path"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51fe9d82-9e8c-409c-8fb7-fe20bc198294",
   "metadata": {},
   "source": [
    "Detect the surface for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b31428c2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sample_surface_detections = {}\n",
    "for sample, df in counts_dfs.items():\n",
    "    surface_counts = defaultdict(int)\n",
    "    for ix, val in df.sum().items():\n",
    "        lane, surface, swath = reverse_get_x(ix+1)\n",
    "        surface_counts[f\"{lane}_{surface}\"] += val\n",
    "\n",
    "    max_lane, max_surface = max(surface_counts, key=surface_counts.get).split(\"_\")[0:2]\n",
    "    sample_surface_detections[sample] = (max_lane, max_surface)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "316d78d5-cdeb-415f-a3de-a1d628bdca3b",
   "metadata": {},
   "source": [
    "Save identify the exact tiles, create and save the whitelist for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d61f1531",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "path = \"output\"\n",
    "\n",
    "for sample, (lane, surface) in sample_surface_detections.items():\n",
    "    # Detect edges tile\n",
    "    THRESHOLD = 2\n",
    "\n",
    "    x = get_x(int(lane), int(surface), 1) - 1\n",
    "    selected_counts = counts_dfs[sample].loc[:, x : x + 5]\n",
    "    selected_counts = selected_counts.loc[:, (selected_counts != 0).any(axis=0)]\n",
    "    selected_counts[selected_counts < THRESHOLD] = 0\n",
    "    first = int(selected_counts.ne(0).idxmax().median()) - 3\n",
    "    last = int(selected_counts.iloc[::-1].ne(0).idxmax().median()) + 3\n",
    "    print(f\"{sample} - {lane}_{surface} - {first} to {last}\")\n",
    "    plt.figure()\n",
    "    cumsum = counts_dfs[sample].loc[:, x : x + 5].sum(axis=1).cumsum()\n",
    "    plt.plot(cumsum)\n",
    "    plt.title(f\"{sample} - {lane}_{surface} - {first} to {last}\")\n",
    "    plt.vlines(first, 0, cumsum.max(), colors=\"red\")\n",
    "    plt.vlines(last, 0, cumsum.max(), colors=\"red\")\n",
    "\n",
    "    extract_whitelist(first, last, path, lane=lane, surface=surface, trim_to_31=True, sample=sample)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e273dc3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.9.5 - BiomedBundle",
   "language": "Python",
   "name": "biomedbundle"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}