2.4.gefCreation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d4fc2f8-4478-45ab-9beb-fc94e7bfed40",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pysam\n",
    "import collections as c\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pickle\n",
    "import numpy as np\n",
    "import cv2 as cv\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import time\n",
    "import os\n",
    "import time\n",
    "import datetime\n",
    "from multiprocessing import Pool\n",
    "import PIL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59c13b9b-24dd-458f-9654-9053fa39c885",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import coo_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c22e5da0-d2bf-41d8-b552-a0b563d27389",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "TILE_DIR = 'HDMI_Tiles_Data'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ac35bc8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tiles = os.listdir(TILE_DIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ed7cddc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tiles = ['_'.join(x.split('_')[:2]) for x in tiles if not 'subset' in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2136d72",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Measured based on a single circle in both H&E and in tile coord space\n",
    "measured_size_px = 2323\n",
    "measured_size_nm = 79763\n",
    "nm_per_coord = measured_size_nm / measured_size_px\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "776e9a34",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "circle_resize_scale = 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a37a7725",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "with open(f'{TILE_DIR}/circle_info.pickle', 'rb') as fh:\n",
    "    ci_info = pickle.load(fh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "936a14aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_final_centroid(data, circle_resize_scale = 25):\n",
    "    final_coords_x = []\n",
    "    final_coords_y = []\n",
    "    if data['circles1'].shape[1] == 8:\n",
    "        final_coords_x.append(int(data['centroid1'][0]) * circle_resize_scale)\n",
    "        final_coords_y.append(int(data['centroid1'][1]) * circle_resize_scale)\n",
    "    if data['circles2'].shape[1] == 8:\n",
    "        final_coords_x.append(int(data['centroid2'][0]) * circle_resize_scale)\n",
    "        final_coords_y.append(int(data['centroid2'][1]) * circle_resize_scale)\n",
    "    if len(final_coords_x) == 0:\n",
    "        raise ValueError('No centroids available for use')\n",
    "\n",
    "    return (np.average(final_coords_y).astype(int), np.average(final_coords_x).astype(int))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d42a412",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# These numbers are based on read measurements from the H&E image and EM images\n",
    "nm_dist_between_circles_y = 310781\n",
    "coords_dist_between_circles_y = int(nm_dist_between_circles_y / nm_per_coord)\n",
    "\n",
    "nm_dist_between_centroids_y = nm_dist_between_circles_y * 4\n",
    "coords_dist_between_centroids_y = int(nm_dist_between_centroids_y / nm_per_coord)\n",
    "\n",
    "nm_dist_between_centroids_x = 1155750\n",
    "coords_dist_between_centroids_x = int(nm_dist_between_centroids_x / nm_per_coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9705c7e4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_dir = \"output\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddf5cbd8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples = ! ls {data_dir}/whitelist*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25f64386",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples = [x.split('/')[-1].split('.')[0] for x in samples]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fde1d67",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples = [x[10:-5] for x in samples]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f66311e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import scipy.io as sio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d63393b9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import gzip\n",
    "from tqdm import tqdm, trange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf93b398",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e26f1dbc-432e-4a3a-99ba-b9ec12644564",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_folder = 'output'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dc203bd-d356-4c28-9a26-4b1b93cd6cc1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfa67e28",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for sample in samples:\n",
    "    _split = sample.split(\"_\")\n",
    "    sample_name = \"_\".join(_split[:-7])\n",
    "    lane, surface, _, start_tile_y, _, end_tile_y = _split[-6:]\n",
    "    lane = int(lane)\n",
    "    surface = int(surface)\n",
    "    start_tile_y = int(start_tile_y)\n",
    "    end_tile_y = int(end_tile_y)\n",
    "\n",
    "    if os.path.exists(\n",
    "        f\"{data_folder}/bc_coords_{lane}_{surface}_tiles_{start_tile_y}_to_{end_tile_y}.pickle\"\n",
    "    ):\n",
    "        with open(\n",
    "            f\"{data_folder}/bc_coords_{lane}_{surface}_tiles_{start_tile_y}_to_{end_tile_y}.pickle\",\n",
    "            \"rb\",\n",
    "        ) as fh:\n",
    "            bc_map = pickle.load(fh)\n",
    "    else:\n",
    "        start_tile_id = f\"{lane}_{surface}{1}{start_tile_y:02}\"\n",
    "        with open(f\"{TILE_DIR}/{start_tile_id}_barcodes.pickle\", \"rb\") as tile_fh:\n",
    "            tile_1_data = pickle.load(tile_fh)\n",
    "        tile_1_ci = ci_info[start_tile_id]\n",
    "        tile_1_cx = int(\n",
    "            ((tile_1_ci[\"centroid1\"][1] * circle_resize_scale) + (tile_1_ci[\"centroid2\"][1]) * circle_resize_scale) / 2\n",
    "        )\n",
    "        tile_1_cy = int(\n",
    "            ((tile_1_ci[\"centroid1\"][0] * circle_resize_scale) + (tile_1_ci[\"centroid2\"][0]) * circle_resize_scale) / 2\n",
    "        )\n",
    "\n",
    "        bcs = list(tile_1_data.keys())\n",
    "        xs = [int(x[0]) for x in tile_1_data.values()]\n",
    "        ys = [int(x[1]) for x in tile_1_data.values()]\n",
    "\n",
    "        for _y_tile in range(start_tile_y, end_tile_y + 1):\n",
    "            for x_tile in range(6):\n",
    "                if _y_tile == 0 and x_tile == 0:\n",
    "                    continue\n",
    "                x_tile = x_tile + 1\n",
    "                y_tile = _y_tile + 1\n",
    "\n",
    "                cur_tile_id = f\"{lane}_{surface}{x_tile}{y_tile:02}\"\n",
    "                with open(f\"{TILE_DIR}/{cur_tile_id}_barcodes.pickle\", \"rb\") as tile_fh:\n",
    "                    cur_tile_data = pickle.load(tile_fh)\n",
    "                cur_tile_ci = ci_info[cur_tile_id]\n",
    "                # TODO: This can have multiple nones\n",
    "                if cur_tile_ci == None:\n",
    "                    # try to average the centroids from the previous and next tiles\n",
    "                    cur_tile_ci = {}\n",
    "                    next_tile_id = cur_tile_id\n",
    "                    next_y_tile = y_tile\n",
    "                    while next_tile_id not in ci_info or ci_info[next_tile_id] == None:\n",
    "                        next_y_tile += 1\n",
    "                        next_tile_id = f\"{lane}_{surface}{x_tile}{next_y_tile:02}\"\n",
    "                    prev_tile_id = cur_tile_id\n",
    "                    prev_y_tile = y_tile\n",
    "                    while prev_tile_id not in ci_info or ci_info[prev_tile_id] == None:\n",
    "                        prev_y_tile -= 1\n",
    "                        prev_tile_id = f\"{lane}_{surface}{x_tile}{prev_y_tile:02}\"\n",
    "\n",
    "                    cur_tile_ci[\"centroid1\"] = (\n",
    "                        (ci_info[next_tile_id][\"centroid1\"][0] + ci_info[prev_tile_id][\"centroid1\"][0]) / 2,\n",
    "                        (ci_info[next_tile_id][\"centroid1\"][1] + ci_info[prev_tile_id][\"centroid1\"][1]) / 2,\n",
    "                    )\n",
    "                    cur_tile_ci[\"centroid2\"] = (\n",
    "                        (ci_info[next_tile_id][\"centroid2\"][0] + ci_info[prev_tile_id][\"centroid2\"][0]) / 2,\n",
    "                        (ci_info[next_tile_id][\"centroid2\"][1] + ci_info[prev_tile_id][\"centroid2\"][1]) / 2,\n",
    "                    )\n",
    "\n",
    "                cur_tile_cx = int(\n",
    "                    (\n",
    "                        (cur_tile_ci[\"centroid1\"][1] * circle_resize_scale)\n",
    "                        + (cur_tile_ci[\"centroid2\"][1]) * circle_resize_scale\n",
    "                    )\n",
    "                    / 2\n",
    "                )\n",
    "                cur_tile_cy = int(\n",
    "                    (\n",
    "                        (cur_tile_ci[\"centroid1\"][0] * circle_resize_scale)\n",
    "                        + (cur_tile_ci[\"centroid2\"][0]) * circle_resize_scale\n",
    "                    )\n",
    "                    / 2\n",
    "                )\n",
    "                y_offset = cur_tile_cy - tile_1_cy\n",
    "                x_offset = cur_tile_cx - tile_1_cx\n",
    "                bcs += list(cur_tile_data.keys())\n",
    "                if x_tile % 2 == 0 and surface == 2:\n",
    "                    for x, y in cur_tile_data.values():\n",
    "                        xs.append(int(x) - x_offset + ((x_tile - 1) * coords_dist_between_centroids_x))\n",
    "                        ys.append(\n",
    "                            (int(y) - y_offset + ((y_tile - start_tile_y - 1) * coords_dist_between_centroids_y))\n",
    "                            + coords_dist_between_circles_y\n",
    "                        )  # Even swaths need moving 1 circle down on thick surface\n",
    "                elif x_tile % 2 == 1 and surface == 1:\n",
    "                    for x, y in cur_tile_data.values():\n",
    "                        xs.append(int(x) - x_offset + ((x_tile - 1) * coords_dist_between_centroids_x))\n",
    "                        ys.append(\n",
    "                            (int(y) - y_offset + ((y_tile - start_tile_y - 1) * coords_dist_between_centroids_y))\n",
    "                            + coords_dist_between_circles_y\n",
    "                        )  # Odd swaths need moving 1 circle down on thin surface\n",
    "                else:\n",
    "                    for x, y in cur_tile_data.values():\n",
    "                        xs.append(int(x) - x_offset + ((x_tile - 1) * coords_dist_between_centroids_x))\n",
    "                        ys.append(int(y) - y_offset + ((y_tile - start_tile_y - 1) * coords_dist_between_centroids_y))\n",
    "        bc_map = {x[0][:31]: x[1] for x in zip(bcs, zip(xs, ys))}\n",
    "\n",
    "        with open(\n",
    "            f\"{data_folder}/bc_coords_{lane}_{surface}_tiles_{start_tile_y}_to_{end_tile_y}.pickle\",\n",
    "            \"wb\",\n",
    "        ) as fh:\n",
    "            pickle.dump(bc_map, fh)\n",
    "\n",
    "    hdmis = pd.read_csv(\n",
    "        f\"{data_folder}/{sample_name}/{sample_name}_Solo.out/GeneFull/raw/barcodes.tsv\",\n",
    "        header=None,\n",
    "    )\n",
    "    mtx = sio.mmread(\n",
    "        f\"{data_folder}/{sample_name}/{sample_name}_Solo.out/GeneFull/raw/matrix.mtx\"\n",
    "    )\n",
    "    with open(\n",
    "        f\"{data_folder}/{sample_name}/{sample_name}_Solo.out/GeneFull/raw/features.tsv\",\n",
    "        \"r\",\n",
    "    ) as fh:\n",
    "        gene_list = [x.split(\"\\t\")[1] for x in fh.readlines()]\n",
    "\n",
    "    mtx_csr = mtx.tocsr()\n",
    "    hdmi_bcs = list(hdmis[0].values)\n",
    "\n",
    "    # Gem creation\n",
    "\n",
    "    header = f\"\"\"#FileFormat=GEMv0.1\n",
    "#SortedBy=None\n",
    "#BinSize=1\n",
    "#STOmicsChip={sample}\n",
    "#OffsetX=0\n",
    "#OffsetY=0\n",
    "geneID\\tx\\ty\\tMIDCount\\tExonCount\\n\"\"\"\n",
    "\n",
    "    with gzip.open(\n",
    "        f\"{data_folder}/{sample_name}/{sample_name}.gem.gz\",\n",
    "        mode=\"wt\",\n",
    "    ) as fh:\n",
    "        fh.write(header)\n",
    "        # Loop over all genes\n",
    "        for gene_idx in trange(len(gene_list)):\n",
    "            gene = gene_list[gene_idx]\n",
    "            gene_data = mtx_csr.getrow(gene_idx)\n",
    "            for cell_idx, gene_count in zip(gene_data.indices, gene_data.data):\n",
    "                try:\n",
    "                    cell_bc = hdmi_bcs[cell_idx]\n",
    "                    x, y = bc_map[cell_bc]\n",
    "                    fh.write(f\"{gene}\\t{x}\\t{y}\\t{gene_count}\\t{gene_count}\\n\")\n",
    "                except (KeyError, IndexError):\n",
    "                    continue\n",
    "    print(\"Finished writing GEM file\")\n",
    "\n",
    "    del mtx_csr\n",
    "\n",
    "    max_x = max([x[0] for x in bc_map.values()])\n",
    "    max_y = max([x[1] for x in bc_map.values()])\n",
    "\n",
    "    gene_sums = mtx.sum(axis=1)\n",
    "    hdmi_sums = mtx.sum(axis=0)\n",
    "    hdmi_mask = hdmi_sums > 0\n",
    "    hdmi_stats = hdmis[np.ravel(hdmi_mask)]\n",
    "    hdmi_stats = hdmi_stats.reindex(hdmi_stats[0])\n",
    "    hdmi_stats = hdmi_stats.drop(0, axis=1)\n",
    "    hdmi_stats[\"Total_Counts\"] = np.ravel(hdmi_sums[hdmi_mask])\n",
    "    genes_detected = (mtx > 0).sum(axis=0)\n",
    "    hdmi_stats[\"Total_Genes\"] = np.ravel(genes_detected[hdmi_mask])\n",
    "\n",
    "    hdmi_stats.plot(kind=\"box\")\n",
    "    hdmi_stats.plot(kind=\"scatter\", x=\"Total_Counts\", y=\"Total_Genes\")\n",
    "\n",
    "    plt.xlim(0, 200)\n",
    "    plt.ylim(0, 200)\n",
    "\n",
    "    mtx_csc = mtx.tocsc()\n",
    "    filtered_mtx = mtx_csc[:, np.ravel(hdmi_mask)]\n",
    "    filtered_mtx_csr = filtered_mtx.tocsr()\n",
    "\n",
    "    nbin_to_500nm = 500 / nm_per_coord\n",
    "    bin_size = int(200 * nbin_to_500nm)\n",
    "\n",
    "    gene_sums = pd.DataFrame(index=gene_list, data=np.ravel(filtered_mtx.sum(axis=1)), columns=[\"Sum\"])\n",
    "    gene_sums.sort_values(by=\"Sum\", ascending=False)\n",
    "\n",
    "    tot_counts = np.ravel(mtx_csc.sum(axis=0))\n",
    "\n",
    "    filtered_bcs = [x[0] for x in zip(hdmis[0], np.ravel(hdmi_mask)) if x[1]]\n",
    "\n",
    "    tile_coo_genes = coo_array((max_x + 1, max_y + 1), dtype=np.int8).toarray()\n",
    "    hdmi_sums = filtered_mtx.sum(axis=0)\n",
    "\n",
    "    del mtx\n",
    "    del mtx_csc\n",
    "    del filtered_mtx\n",
    "    del filtered_mtx_csr\n",
    "\n",
    "    hdmi_sums = hdmi_sums.astype(np.uint16)\n",
    "    missing_bcs = set()\n",
    "\n",
    "    for n, bc in enumerate(filtered_bcs):\n",
    "        if bc in bc_map:\n",
    "            x, y = bc_map[bc]\n",
    "            tile_coo_genes[x, y] = hdmi_sums[0, n]\n",
    "        else:\n",
    "            missing_bcs.add(bc)\n",
    "\n",
    "    del bc_map\n",
    "\n",
    "    print(f\"Missing BCs: {len(missing_bcs)}\")\n",
    "    del missing_bcs\n",
    "\n",
    "    # Bin the data in 25x25 bins\n",
    "    bin_size = 25\n",
    "\n",
    "    # function to bin 2d matrix\n",
    "    def bin2d(a, K):\n",
    "        m_bins = a.shape[0] // K\n",
    "        n_bins = a.shape[1] // K\n",
    "        return a[: m_bins * K, : n_bins * K].reshape(m_bins, K, n_bins, K).sum(3).sum(1)\n",
    "\n",
    "    # Coo to np\n",
    "    tile_coo_genes_np = np.array(tile_coo_genes)\n",
    "\n",
    "    # Bin\n",
    "    gene_binned = bin2d(tile_coo_genes_np, bin_size)\n",
    "    # add min and then rescale to have max at 255\n",
    "    if gene_binned.min() < 0:\n",
    "        gene_binned = gene_binned + abs(gene_binned.min())\n",
    "\n",
    "    # gene_binned = (gene_binned / gene_binned.max()) * 255\n",
    "    gene_binned = (gene_binned / np.percentile(gene_binned, 0.999995)) * 255\n",
    "\n",
    "    gene_binned = np.clip(gene_binned, 0, 255)\n",
    "\n",
    "    PIL.Image.fromarray(gene_binned.astype(np.uint8), \"L\").save(\n",
    "        f\"{data_folder}/{sample_name}/GeneFull_HDMI_Locs_Global_HDMIs.png\"\n",
    "    )"
   ]
  }
 ],
 "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
}