[0ad989]: / train_neural_network / agg_fig_data.ipynb

Download this file

373 lines (372 with data), 23.4 kB

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "78ad58c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import anndata as ad\n",
    "import os\n",
    "import scipy\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sys\n",
    "import multivelo as mv\n",
    "import scanpy as sc\n",
    "import scvelo as scv\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sys.path.append(\"/..\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "923af518-b096-4e63-94e5-06f08c944fee",
   "metadata": {},
   "source": [
    "## Read in the Appropriate Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51361505-5c31-4f11-9212-9a5cfab0f6d1",
   "metadata": {},
   "source": [
    "Uncomment the top cell to aggragate data for training. Uncomment the bottom cell to aggregate data for validation. (You will need to supply your own AnnData object of validation data.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ebbb4e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# training data\n",
    "\n",
    "# data used to generate figures in paper\n",
    "fig3 = sc.read_h5ad(\"../Examples/multivelo_result_fig3.h5ad\")\n",
    "fig4 = sc.read_h5ad(\"../Examples/multivelo_result_fig4.h5ad\")\n",
    "fig5 = sc.read_h5ad(\"../Examples/multivelo_result_fig5.h5ad\")\n",
    "fig6 = sc.read_h5ad(\"../Examples/multivelo_result_fig6.h5ad\")\n",
    "\n",
    "figs = [fig3, fig4, fig5, fig6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "13dca1bb-4ab9-4168-a696-fa5e3bb305d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# vali data\n",
    "# You will need to supply your own validation data in order to run this code!\n",
    "# val_data = sc.read_h5ad()\n",
    "# figs = [val_data]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "46b3767a-63fa-43c7-ad58-4af3b67357fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = len(figs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34e2fb82-18e3-4fd0-922a-a0187ae68b2b",
   "metadata": {},
   "source": [
    "## Define Appropriate Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dc5347e3-21bf-4d87-94fa-ae96f264ace5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# graph a set of rate parameters\n",
    "def graph_params(alpha_c, alpha, beta, gamma):\n",
    "\n",
    "    # the names of the parameters for axis labels\n",
    "    names = [\"alpha_c\", \"alpha\", \"beta\", \"gamma\"]\n",
    "\n",
    "    # create a subplot\n",
    "    fig, axs = plt.subplots(4, 3, figsize=(12, 12))\n",
    "\n",
    "    # axis boundaries for each rate parameter\n",
    "    lims = [0.2, 2, 1.25, 2]\n",
    "\n",
    "    # assemble the rate parameters into a matrix\n",
    "    mvln = np.array([alpha_c, alpha, beta, gamma])\n",
    "\n",
    "    # for each combination of rate parameters...\n",
    "    for i in range(4):\n",
    "        col = 0\n",
    "        for j in range(4):\n",
    "\n",
    "            # if we're about to plot a rate parameter against\n",
    "            # itself, then skip this iteration\n",
    "            if i == j:\n",
    "                continue\n",
    "\n",
    "            # graph the rate parameters\n",
    "            h = axs[i][col].hist2d(np.ravel(mvln[i,:]), np.ravel(mvln[j,:]),\n",
    "                                   # range=[[0, lims[i]],\n",
    "                                   #        [0, lims[j]]],\n",
    "                                   bins=10,\n",
    "                                   cmap=\"Greens\")\n",
    "            axs[i][col].set_xlabel(names[i])\n",
    "            axs[i][col].set_ylabel(names[j])\n",
    "            fig.colorbar(h[3], ax=axs[i][col])\n",
    "            col += 1\n",
    "    \n",
    "    # # fig.colorbar(axs)\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "936b562d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# a function for filtering out genes that fall below a certain fit_likelihood threshold\n",
    "def filter_likelihood(adata, thresh):\n",
    "\n",
    "    # print out the shape of the AnnData object before filtering\n",
    "    print(\"Before shape:\", adata.shape)\n",
    "\n",
    "    # get the fit likelihood of each gene\n",
    "    fig_likelihood = adata.var['fit_likelihood']\n",
    "\n",
    "    # draw a histogram of the fit likelihood\n",
    "    plt.hist(fig_likelihood, bins=40, range=(0, 0.2))\n",
    "\n",
    "    # do the filtering\n",
    "    filtered = fig_likelihood > thresh\n",
    "    return_val = adata[:, filtered]\n",
    "\n",
    "    # print out the shape of the AnnData object after filtering\n",
    "    print(\"After shape:\", return_val.shape)\n",
    "    print()\n",
    "    \n",
    "    return return_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a591b8fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# a function for choosing a subset of the data to use\n",
    "# like if we wanted to only save enough for a validation set, for example\n",
    "def subset_data(adata, n=None):\n",
    "\n",
    "    # if a subset size isn't specified just\n",
    "    # return the original dataset\n",
    "    if n is None:\n",
    "        return adata\n",
    "\n",
    "    total_genes = adata.shape[1]\n",
    "    \n",
    "    full_data = range(total_genes)\n",
    "\n",
    "    # make a random choice of indices\n",
    "    choice = np.random.choice(total_genes, size=n, replace=False)\n",
    "\n",
    "    # subset the AnnData object with our random set\n",
    "    subset_adata = adata[:,choice]\n",
    "    \n",
    "    return subset_adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "eef193fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function responsible for saving the relevant data\n",
    "def write_files(outfile, adatas, graph=False):\n",
    "\n",
    "        # print the name of the file we're saving to\n",
    "        print(outfile)\n",
    "\n",
    "        alpha_c = np.array([])\n",
    "        alpha = np.array([])\n",
    "        beta = np.array([])\n",
    "        gamma = np.array([])\n",
    "\n",
    "        # assemble the rate parameters of each AnnData object\n",
    "        for adata in adatas:\n",
    "        \n",
    "            alpha_c = np.concatenate((alpha_c, np.array(adata.var['fit_alpha_c'])))\n",
    "            alpha = np.concatenate((alpha, np.array(adata.var['fit_alpha'])))\n",
    "            beta = np.concatenate((beta, np.array(adata.var['fit_beta'])))\n",
    "            gamma = np.concatenate((gamma, np.array(adata.var['fit_gamma'])))\n",
    "\n",
    "        # graph results if the user specifies it\n",
    "        if graph:\n",
    "            graph_results(alpha_c, alpha, beta, gamma)\n",
    "\n",
    "        # save all of the data\n",
    "        np.savetxt(outfile + \"/alpha_c.txt\", alpha_c)\n",
    "        np.savetxt(outfile + \"/alpha.txt\", alpha)\n",
    "        np.savetxt(outfile + \"/beta.txt\", beta)\n",
    "        np.savetxt(outfile + \"/gamma.txt\", gamma)\n",
    "\n",
    "        # print the total number of genes saved\n",
    "        print(\"Number of genes:\", alpha_c.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "780c9b89-5cf8-406d-adde-227cbd2af61b",
   "metadata": {},
   "source": [
    "## Process and Save Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e966528d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before shape: (3365, 865)\n",
      "After shape: (3365, 665)\n",
      "\n",
      "Before shape: (6436, 960)\n",
      "After shape: (6436, 771)\n",
      "\n",
      "Before shape: (11605, 936)\n",
      "After shape: (11605, 655)\n",
      "\n",
      "Before shape: (4693, 747)\n",
      "After shape: (4693, 507)\n",
      "\n",
      "The number of remaining genes is: 26099\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAk5UlEQVR4nO3de3BU9f3/8VcCJEFMNgRMNlsTiIqCGrygxIi3SkZARqEyRZRWUApegoqIQjpy0aqJaJHBIqijgOMFsVOx1YrFIF5DkAiKihFpKqhsUGk2XCRc8vn+4Y/z65oo2XBO9rPr8zGzU3LOZ8++3z2b7MvPuWyCMcYIAADAIonRLgAAAODHCCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOu0j3YBrdHY2Kivv/5aqampSkhIiHY5AACgBYwx2rFjhwKBgBITf36OJCYDytdff62cnJxolwEAAFphy5YtOvroo392TEwGlNTUVEk/NJiWlhblagAAQEvU19crJyfH+Rz/OTEZUA4e1klLSyOgAAAQY1pyegYnyQIAAOsQUAAAgHUIKAAAwDoRB5Q333xTl1xyiQKBgBISErR06dKw9cYYTZs2TdnZ2erYsaOKioq0cePGsDHbt2/XyJEjlZaWpvT0dI0ZM0Y7d+48rEYAAED8iDig7Nq1S6eccormzp3b7PqZM2dqzpw5mj9/viorK9WpUycNGDBAe/bsccaMHDlSH3/8sZYvX66XXnpJb775psaNG9f6LgAAQFxJMMaYVj85IUEvvPCChg4dKumH2ZNAIKBbb71VkyZNkiSFQiFlZWVp4cKFGjFihDZs2KATTzxR7733ns444wxJ0rJly3TxxRfryy+/VCAQOOTr1tfXy+fzKRQKcRUPAAAxIpLPb1fPQampqVEwGFRRUZGzzOfzqaCgQBUVFZKkiooKpaenO+FEkoqKipSYmKjKyspmt9vQ0KD6+vqwBwAAiF+uBpRgMChJysrKCluelZXlrAsGg8rMzAxb3759e2VkZDhjfqy0tFQ+n895cBdZAADiW0xcxVNSUqJQKOQ8tmzZEu2SAACAh1wNKH6/X5JUW1sbtry2ttZZ5/f7tW3btrD1+/fv1/bt250xP5acnOzcNZa7xwIAEP9cDSh5eXny+/0qLy93ltXX16uyslKFhYWSpMLCQtXV1amqqsoZs2LFCjU2NqqgoMDNcgAAQIyK+Lt4du7cqc8//9z5uaamRuvWrVNGRoZyc3M1YcIE3X333erRo4fy8vI0depUBQIB50qfXr16aeDAgRo7dqzmz5+vffv2afz48RoxYkSLruABAADxL+KAsmbNGv361792fp44caIkadSoUVq4cKFuv/127dq1S+PGjVNdXZ3OOeccLVu2TCkpKc5znn76aY0fP179+/dXYmKihg0bpjlz5rjQDgAAiAeHdR+UaOE+KAAAxJ5IPr8jnkEBrDfDF+H4kDd1AABaLSYuMwYAAL8sBBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWKe92xs8cOCAZsyYoaeeekrBYFCBQECjR4/WHXfcoYSEBEmSMUbTp0/XY489prq6OvXr10/z5s1Tjx493C4H8WKGL9oVAADakOszKPfdd5/mzZunv/zlL9qwYYPuu+8+zZw5Uw899JAzZubMmZozZ47mz5+vyspKderUSQMGDNCePXvcLgcAAMQg12dQ3n33XQ0ZMkSDBw+WJHXv3l3PPvusVq9eLemH2ZPZs2frjjvu0JAhQyRJTz75pLKysrR06VKNGDHC7ZIAAECMcX0G5eyzz1Z5ebk+++wzSdIHH3ygt99+W4MGDZIk1dTUKBgMqqioyHmOz+dTQUGBKioqmt1mQ0OD6uvrwx4AACB+uT6DMmXKFNXX16tnz55q166dDhw4oHvuuUcjR46UJAWDQUlSVlZW2POysrKcdT9WWlqqO++80+1SgR9Ecn7LjJB3dQAAHK7PoCxZskRPP/20nnnmGb3//vtatGiRHnjgAS1atKjV2ywpKVEoFHIeW7ZscbFiAABgG9dnUG677TZNmTLFOZckPz9fX3zxhUpLSzVq1Cj5/X5JUm1trbKzs53n1dbW6tRTT212m8nJyUpOTna7VAAAYCnXZ1B2796txMTwzbZr106NjY2SpLy8PPn9fpWXlzvr6+vrVVlZqcLCQrfLAQAAMcj1GZRLLrlE99xzj3Jzc3XSSSdp7dq1mjVrlq655hpJUkJCgiZMmKC7775bPXr0UF5enqZOnapAIKChQ4e6XQ4AAIhBrgeUhx56SFOnTtUNN9ygbdu2KRAI6Nprr9W0adOcMbfffrt27dqlcePGqa6uTuecc46WLVumlJQUt8sBAAAxKMEYY6JdRKTq6+vl8/kUCoWUlpYW7XLQFmy5kyxX8QBAq0Xy+e36DAoQ17gkGQDaBF8WCAAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADW4Vb3iB5bvl8HAGAdZlAAAIB1CCgAAMA6HOJpQ/mL8ls8dv2o9R5W4iEO2wAAXMAMCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOtwH5TBFcm8TAADQMsygAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFjHk4Dy1Vdf6Xe/+526dOmijh07Kj8/X2vWrHHWG2M0bdo0ZWdnq2PHjioqKtLGjRu9KAUAAMQg1wPKf//7X/Xr108dOnTQK6+8ok8++UR//vOf1blzZ2fMzJkzNWfOHM2fP1+VlZXq1KmTBgwYoD179rhdDgAAiEHt3d7gfffdp5ycHC1YsMBZlpeX5/zbGKPZs2frjjvu0JAhQyRJTz75pLKysrR06VKNGDHC7ZIAAECMcX0G5e9//7vOOOMM/fa3v1VmZqZOO+00PfbYY876mpoaBYNBFRUVOct8Pp8KCgpUUVHR7DYbGhpUX18f9gAAAPHL9YDy73//W/PmzVOPHj306quv6vrrr9dNN92kRYsWSZKCwaAkKSsrK+x5WVlZzrofKy0tlc/ncx45OTlulw0AACziekBpbGzU6aefrnvvvVennXaaxo0bp7Fjx2r+/Pmt3mZJSYlCoZDz2LJli4sVAwAA27geULKzs3XiiSeGLevVq5c2b94sSfL7/ZKk2trasDG1tbXOuh9LTk5WWlpa2AMAAMQv1wNKv379VF1dHbbss88+U7du3ST9cMKs3+9XeXm5s76+vl6VlZUqLCx0uxwAABCDXL+K55ZbbtHZZ5+te++9V8OHD9fq1av16KOP6tFHH5UkJSQkaMKECbr77rvVo0cP5eXlaerUqQoEAho6dKjb5QAAgBjkekA588wz9cILL6ikpER33XWX8vLyNHv2bI0cOdIZc/vtt2vXrl0aN26c6urqdM4552jZsmVKSUlxuxwAABCDEowxJtpFRKq+vl4+n0+hUCjq56PkL8r3ZLvrR633ZLuem+GLdgX2mBGKdgUAYJVIPr/5Lh4AAGAdAgoAALAOAQUAAFiHgAIAAKzj+lU8AP6fSE4Y5oRaAAjDDAoAALAOAQUAAFiHQzxtaEnp/pYPHuVdHQAA2I4ZFAAAYB0CCgAAsA4BBQAAWIdzUCy1oWeviMb3+nSDR5UAAND2mEEBAADWIaAAAADrcIgHcSc/Lzei8etrNntUCQCgtZhBAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh/ugIGoiuV8J9yoBgF8WZlAAAIB1CCgAAMA6HOLBz8pflB/R+PUe1QEA+GVhBgUAAFiHgAIAAKxDQAEAANbhHJQ4Ecm5IutHeXemSCSXDtuCy50BwD7MoAAAAOsQUAAAgHU4xPMLtKFnr5YPLuEtAgBoe8ygAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh2tIm9F9ysstHpsawRW7vwRLSve3eOxwLmEGAPwEZlAAAIB1CCgAAMA6BBQAAGAdTgI4TJGcc4HWi8VvSQYAtB4zKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArENAAQAA1uEy4zjB5c4AgHjCDAoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsI7nAaWsrEwJCQmaMGGCs2zPnj0qLi5Wly5ddOSRR2rYsGGqra31uhQAABAjPL3V/XvvvadHHnlEvXv3Dlt+yy236OWXX9bzzz8vn8+n8ePH67LLLtM777zjZTmAvWb4Ihwf8qYOALCEZzMoO3fu1MiRI/XYY4+pc+fOzvJQKKTHH39cs2bN0oUXXqg+ffpowYIFevfdd7Vq1SqvygEAADHEs4BSXFyswYMHq6ioKGx5VVWV9u3bF7a8Z8+eys3NVUVFhVflAACAGOLJIZ7Fixfr/fff13vvvddkXTAYVFJSktLT08OWZ2VlKRgMNru9hoYGNTQ0OD/X19e7Wi8AALCL6zMoW7Zs0c0336ynn35aKSkprmyztLRUPp/PeeTk5LiyXQAAYCfXA0pVVZW2bdum008/Xe3bt1f79u31xhtvaM6cOWrfvr2ysrK0d+9e1dXVhT2vtrZWfr+/2W2WlJQoFAo5jy1btrhdNgAAsIjrh3j69++v9evXhy27+uqr1bNnT02ePFk5OTnq0KGDysvLNWzYMElSdXW1Nm/erMLCwma3mZycrOTkZLdLBQAAlnI9oKSmpurkk08OW9apUyd16dLFWT5mzBhNnDhRGRkZSktL04033qjCwkKdddZZbpfTKq8snRTtEgAA+EXz9D4oP+XBBx9UYmKihg0bpoaGBg0YMEAPP/xwNEoBAAAWapOAsnLlyrCfU1JSNHfuXM2dO7ctXh4AAMQYvosHAABYJyqHeBA7lpTuj3YJVsnPy23x2PU1mz2sBADiGzMoAADAOgQUAABgHQ7xIGq8Onw0vMSOtzWHgwCg9ZhBAQAA1iGgAAAA6xBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYJ320S4AiLYlpftbPHZ4iSW/MjN8EYwNeVcHAHiEGRQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOtYcs0k4J5ILhsGANiJGRQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA63ugc8EtEt90d4VwcAxCJmUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA4BBQAAWIeAAgAArMNlxoAF8vNyIxq/vmazR5UAgB2YQQEAANYhoAAAAOsQUAAAgHUIKAAAwDoEFAAAYB0CCgAAsA6XGQMRiOgbii2Rvyi/xWN3bCiLaNv/KRscaTkA0CLMoAAAAOsQUAAAgHVcDyilpaU688wzlZqaqszMTA0dOlTV1dVhY/bs2aPi4mJ16dJFRx55pIYNG6ba2lq3SwEAADHK9YDyxhtvqLi4WKtWrdLy5cu1b98+XXTRRdq1a5cz5pZbbtE//vEPPf/883rjjTf09ddf67LLLnO7FAAAEKNcP0l22bJlYT8vXLhQmZmZqqqq0nnnnadQKKTHH39czzzzjC688EJJ0oIFC9SrVy+tWrVKZ511ltslAQCAGOP5OSihUEiSlJGRIUmqqqrSvn37VFRU5Izp2bOncnNzVVFR0ew2GhoaVF9fH/YAAADxy9PLjBsbGzVhwgT169dPJ598siQpGAwqKSlJ6enpYWOzsrIUDAab3U5paanuvPNOL0sFYsqGxYGWDy7xrg4A8IqnMyjFxcX66KOPtHjx4sPaTklJiUKhkPPYsmWLSxUCAAAbeTaDMn78eL300kt68803dfTRRzvL/X6/9u7dq7q6urBZlNraWvn9/ma3lZycrOTkZK9KBQAAlnF9BsUYo/Hjx+uFF17QihUrlJeXF7a+T58+6tChg8rLy51l1dXV2rx5swoLC90uBwAAxCDXZ1CKi4v1zDPP6MUXX1RqaqpzXonP51PHjh3l8/k0ZswYTZw4URkZGUpLS9ONN96owsJCruABPBDJ7fkHDfWuDgCIhOsBZd68eZKkCy64IGz5ggULNHr0aEnSgw8+qMTERA0bNkwNDQ0aMGCAHn74YbdLAQAAMcr1gGKMOeSYlJQUzZ07V3PnznX75QEAQBzgu3gAAIB1CCgAAMA6BBQAAGAdAgoAALCOp7e6BxDfuk95ucVjU3tNafHY9aPWt6YcAHGEGRQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOtwmTEA6+Qvym/xWC5JBuITMygAAMA6BBQAAGAdAgoAALAOAQUAAFiHgAIAAKxDQAEAANYhoAAAAOtwHxQAjleWTvJs28N72fHnpvuUl1s89j9lg6O+XbQd9qFdmEEBAADWIaAAAADr2DHnCiDuLSnd3+Kxw0v40wT80jGDAgAArENAAQAA1iGgAAAA63CgF0BMy1+UH+Ezyjyp45eAy3DRlphBAQAA1iGgAAAA63CIB7BAJJfg4vCk9poSwWgOUwDRwgwKAACwDgEFAABYh4ACAACswzkoANDGIrlcV/Lukl0u0YbNmEEBAADWIaAAAADrEFAAAIB1OAcFAH5CZOdotPz8jFeWToqskDi/bbyXt9D3ah/Ce8ygAAAA6xBQAACAdTjEA8A6Xt76f3iJN3/2IrqF/tLIts23COOXiBkUAABgHQIKAACwDgEFAABYh3NQAOAnRHIujFfntkgRnt+i+D4HJdKvCUjtFf06OC+odZhBAQAA1iGgAAAA63CIB8AvipeXMHslkpq7b2j5oQevDn94KbLDXfhftnyLdksxgwIAAKxDQAEAANYhoAAAAOtwDgoAuMCWc1si+abk4b0i+wiI6FuYubTWEdk3KkvrR633ZNuvLI3wPco5KAAAAOEIKAAAwDoEFAAAYB3OQQEAuC6y8y7KWjwy0nN9vPwKAq9Ees5KvGIGBQAAWIeAAgAArBN7c18AAOtFdigmgsuXPRTRZdQRiPRy7kjYcnm7F5hBAQAA1iGgAAAA60Q1oMydO1fdu3dXSkqKCgoKtHr16miWAwAALBG1c1Cee+45TZw4UfPnz1dBQYFmz56tAQMGqLq6WpmZmdEqCwAAV8XzeSJeitoMyqxZszR27FhdffXVOvHEEzV//nwdccQReuKJJ6JVEgAAsERUZlD27t2rqqoqlZSUOMsSExNVVFSkioqKJuMbGhrU0NDg/BwKhSRJ9fX1ntS388ABT7YLADY58H1CRONj8W9jJD3GYn9e8uIz9uA2jTGHHBuVgPLtt9/qwIEDysrKCluelZWlTz/9tMn40tJS3XnnnU2W5+TkeFYjAMS96yMb3tebKrwVQY8x2Z+XfD7PNr1jxw75DrH9mLgPSklJiSZOnOj83NjYqO3bt6tLly5KSIjsvwAOpb6+Xjk5OdqyZYvS0tJc3bYN6C/2xXuP9Bf74r3HeO9P8q5HY4x27NihQCBwyLFRCShdu3ZVu3btVFtbG7a8trZWfr+/yfjk5GQlJyeHLUtPT/eyRKWlpcXtG0+iv3gQ7z3SX+yL9x7jvT/Jmx4PNXNyUFROkk1KSlKfPn1UXl7uLGtsbFR5ebkKCwujURIAALBI1A7xTJw4UaNGjdIZZ5yhvn37avbs2dq1a5euvvrqaJUEAAAsEbWAcvnll+ubb77RtGnTFAwGdeqpp2rZsmVNTpxta8nJyZo+fXqTQ0rxgv5iX7z3SH+xL957jPf+JDt6TDAtudYHAACgDfFdPAAAwDoEFAAAYB0CCgAAsA4BBQAAWCfuAsrcuXPVvXt3paSkqKCgQKtXr/7Z8c8//7x69uyplJQU5efn65///GfYemOMpk2bpuzsbHXs2FFFRUXauHFj2Jjt27dr5MiRSktLU3p6usaMGaOdO3e63pvkbn/79u3T5MmTlZ+fr06dOikQCOiqq67S119/HbaN7t27KyEhIexRVlbmSX+S+/tw9OjRTeofOHBg2JhY3YeSmvR28HH//fc7Y9pyH0bS38cff6xhw4Y59c2ePbtV29yzZ4+Ki4vVpUsXHXnkkRo2bFiTG0G6ye0eS0tLdeaZZyo1NVWZmZkaOnSoqqurw8ZccMEFTfbhdddd53Zrktzvb8aMGU1q79mzZ9iYttyHbvfX3O9XQkKCiouLnTFtuf+kyHp87LHHdO6556pz587q3LmzioqKmoyPymehiSOLFy82SUlJ5oknnjAff/yxGTt2rElPTze1tbXNjn/nnXdMu3btzMyZM80nn3xi7rjjDtOhQwezfv16Z0xZWZnx+Xxm6dKl5oMPPjCXXnqpycvLM99//70zZuDAgeaUU04xq1atMm+99ZY57rjjzBVXXGF9f3V1daaoqMg899xz5tNPPzUVFRWmb9++pk+fPmHb6datm7nrrrvM1q1bncfOnTtd78+LHo0xZtSoUWbgwIFh9W/fvj1sO7G6D40xYX1t3brVPPHEEyYhIcFs2rTJGdNW+zDS/lavXm0mTZpknn32WeP3+82DDz7Yqm1ed911Jicnx5SXl5s1a9aYs846y5x99tmu9+dVjwMGDDALFiwwH330kVm3bp25+OKLTW5ubtg+Ov/8883YsWPD9mEoFIqJ/qZPn25OOumksNq/+eabsDFttQ+96G/btm1hvS1fvtxIMq+//rozpq32X2t6vPLKK83cuXPN2rVrzYYNG8zo0aONz+czX375pTMmGp+FcRVQ+vbta4qLi52fDxw4YAKBgCktLW12/PDhw83gwYPDlhUUFJhrr73WGGNMY2Oj8fv95v7773fW19XVmeTkZPPss88aY4z55JNPjCTz3nvvOWNeeeUVk5CQYL766ivXejPG/f6as3r1aiPJfPHFF86ybt26NftL6QUvehw1apQZMmTIT75mvO3DIUOGmAsvvDBsWVvtw0j7+18/VeOhtllXV2c6dOhgnn/+eWfMhg0bjCRTUVFxGN00z4sef2zbtm1GknnjjTecZeeff765+eabW1NyRLzob/r06eaUU075yee15T5si/138803m2OPPdY0NjY6y9pq/xlzeD0aY8z+/ftNamqqWbRokTEmep+FcXOIZ+/evaqqqlJRUZGzLDExUUVFRaqoqGj2ORUVFWHjJWnAgAHO+JqaGgWDwbAxPp9PBQUFzpiKigqlp6frjDPOcMYUFRUpMTFRlZWVVvfXnFAopISEhCbfdVRWVqYuXbrotNNO0/3336/9+/e3vpmf4GWPK1euVGZmpk444QRdf/31+u6778K2ES/7sLa2Vi+//LLGjBnTZJ3X+7A1/bmxzaqqKu3bty9sTM+ePZWbm9vq1z2cetwQCoUkSRkZGWHLn376aXXt2lUnn3yySkpKtHv3btdeU/K2v40bNyoQCOiYY47RyJEjtXnzZmddW+3Dtth/e/fu1VNPPaVrrrmmyZfZer3/Dr7+4fa4e/du7du3z3n/ReuzMCa+zbglvv32Wx04cKDJnWizsrL06aefNvucYDDY7PhgMOisP7js58ZkZmaGrW/fvr0yMjKcMW7wor8f27NnjyZPnqwrrrgi7MuhbrrpJp1++unKyMjQu+++q5KSEm3dulWzZs06zK7CedXjwIEDddlllykvL0+bNm3SH//4Rw0aNEgVFRVq165dXO3DRYsWKTU1VZdddlnY8rbYh63pz41tBoNBJSUlNQnVP/f/U2t50eOPNTY2asKECerXr59OPvlkZ/mVV16pbt26KRAI6MMPP9TkyZNVXV2tv/3tb668ruRdfwUFBVq4cKFOOOEEbd26VXfeeafOPfdcffTRR0pNTW2zfdgW+2/p0qWqq6vT6NGjw5a3xf6T3Olx8uTJCgQCTiCJ1mdh3AQUHJ59+/Zp+PDhMsZo3rx5YesmTpzo/Lt3795KSkrStddeq9LS0pi41fOIESOcf+fn56t379469thjtXLlSvXv3z+KlbnviSee0MiRI5WSkhK2PNb34S9JcXGxPvroI7399tthy8eNG+f8Oz8/X9nZ2erfv782bdqkY489tq3LjMigQYOcf/fu3VsFBQXq1q2blixZ0uxsXyx7/PHHNWjQIAUCgbDlsbL/ysrKtHjxYq1cubLJ35G2FjeHeLp27ap27do1Oeu7trZWfr+/2ef4/f6fHX/wfw81Ztu2bWHr9+/fr+3bt//k67aGF/0ddDCcfPHFF1q+fPkhv1q7oKBA+/fv13/+85/IG/kZXvb4v4455hh17dpVn3/+ubONWN+HkvTWW2+purpaf/jDHw5Zixf7sDX9ubFNv9+vvXv3qq6uzrXXPZx6Dsf48eP10ksv6fXXX9fRRx/9s2MLCgokyXkfu8Hr/g5KT0/X8ccfH/Y72Bb70Ov+vvjiC7322mst/h2U3N1/0uH1+MADD6isrEz/+te/1Lt3b2d5tD4L4yagJCUlqU+fPiovL3eWNTY2qry8XIWFhc0+p7CwMGy8JC1fvtwZn5eXJ7/fHzamvr5elZWVzpjCwkLV1dWpqqrKGbNixQo1NjY6b0Bb+5P+fzjZuHGjXnvtNXXp0uWQtaxbt06JiYlNpvMOl1c9/tiXX36p7777TtnZ2c42YnkfHvT444+rT58+OuWUUw5Zixf7sDX9ubHNPn36qEOHDmFjqqurtXnz5la/7uHU0xrGGI0fP14vvPCCVqxYoby8vEM+Z926dZLkvI/d4FV/P7Zz505t2rTJqb2t9qHX/S1YsECZmZkaPHjwIcd6sf+k1vc4c+ZM/elPf9KyZcvCziORovhZ2KpTay21ePFik5ycbBYuXGg++eQTM27cOJOenm6CwaAxxpjf//73ZsqUKc74d955x7Rv39488MADZsOGDWb69OnNXmacnp5uXnzxRfPhhx+aIUOGNHtp1WmnnWYqKyvN22+/bXr06OHZJapu9rd3715z6aWXmqOPPtqsW7cu7PK3hoYGY4wx7777rnnwwQfNunXrzKZNm8xTTz1ljjrqKHPVVVe53p8XPe7YscNMmjTJVFRUmJqaGvPaa6+Z008/3fTo0cPs2bPH2U6s7sODQqGQOeKII8y8efOavGZb7sNI+2toaDBr1641a9euNdnZ2WbSpElm7dq1ZuPGjS3epjE/XKKam5trVqxYYdasWWMKCwtNYWGh6/151eP1119vfD6fWblyZdjv4e7du40xxnz++efmrrvuMmvWrDE1NTXmxRdfNMccc4w577zzYqK/W2+91axcudLU1NSYd955xxQVFZmuXbuabdu2OWPaah960Z8xP1wpk5ubayZPntzkNdty/7Wmx7KyMpOUlGT++te/hr3/duzYETamrT8L4yqgGGPMQw89ZHJzc01SUpLp27evWbVqlbPu/PPPN6NGjQobv2TJEnP88cebpKQkc9JJJ5mXX345bH1jY6OZOnWqycrKMsnJyaZ///6muro6bMx3331nrrjiCnPkkUeatLQ0c/XVV4ftWFv7q6mpMZKafRy8fr+qqsoUFBQYn89nUlJSTK9evcy9994b9uFuc4+7d+82F110kTnqqKNMhw4dTLdu3czYsWPDPtyMid19eNAjjzxiOnbsaOrq6pqsa+t9GEl/P/UePP/881u8TWOM+f77780NN9xgOnfubI444gjzm9/8xmzdutWT/rzo8ad+DxcsWGCMMWbz5s3mvPPOMxkZGSY5Odkcd9xx5rbbbvPsPhpu93f55Zeb7Oxsk5SUZH71q1+Zyy+/3Hz++edhr9mW+9CL9+irr75qJDX5fDCm7fefMZH12K1bt2Z7nD59ujMmGp+FCcYY07q5FwAAAG/EzTkoAAAgfhBQAACAdQgoAADAOgQUAABgHQIKAACwDgEFAABYh4ACAACsQ0ABAADWIaAAAADrEFAAAIB1CCgAAMA6BBQAAGCd/wMPclQBKkxoegAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# filter out genes that have a fit likelihood lower\n",
    "# than this specified value\n",
    "\n",
    "# used on fig data:\n",
    "# likelihood_thresh = 0.0\n",
    "\n",
    "# used on hspc data:\n",
    "likelihood_thresh = 0.025\n",
    "\n",
    "# keep a running total of the final number of genes\n",
    "new_gene_num = 0\n",
    "\n",
    "# For each AnnData object we're pulling data from,\n",
    "# filter out the bad data\n",
    "for i in range(N):\n",
    "    \n",
    "    figs[i] = filter_likelihood(figs[i], likelihood_thresh)\n",
    "    new_gene_num += figs[i].shape[0]\n",
    "\n",
    "print(\"The number of remaining genes is:\", new_gene_num)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "13590e02",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_figs = []\n",
    "\n",
    "# Set subset_size to number of genes per\n",
    "# AnnData object you want if you're generating\n",
    "# validation data\n",
    "# subset_size = 30\n",
    "\n",
    "subset_size = None\n",
    "\n",
    "# For each AnnData object we're pulling data from,\n",
    "# pull out a subset:\n",
    "for i in range(N):\n",
    "    \n",
    "    test_fig = subset_data(figs[i],  subset_size)\n",
    "    test_figs.append(test_fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b86e54b4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./data/rates/fig_rates_train\n",
      "Number of genes: 2598\n"
     ]
    }
   ],
   "source": [
    "# write the data to the specified file:\n",
    "# write_files(\"./data/rates/val_rates\", test_figs)\n",
    "write_files(\"./data/rates/fig_rates_train\", test_figs)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}