norman_scETM_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b8270f66-7a27-4263-aa97-25abeb2cf34d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "529f64f3-c83e-481d-872c-df89b9dd063b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from utils import (\n",
    "    auprc,\n",
    "    factor_enrichment_gsea,\n",
    "    get_gprofiler,\n",
    "    read_aws_csv,\n",
    "    read_aws_h5ad,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00d0db8a-5c3d-4e9f-bce0-46c90135cab4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in scETM results\n",
    "adata = read_aws_h5ad(\"s3://pert-spectra/scETM_checkpoints/scetm_norman/fold_0.h5ad\")\n",
    "# read in pertspectra results to retrieve gene labels\n",
    "ref_adata = read_aws_h5ad(\n",
    "    \"s3://pert-spectra/PertSpectra_checkpoints/pertspectra_norman/fold_0.h5ad\"\n",
    ")\n",
    "adata.var_names = ref_adata.var_names"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b7bf773-5e0c-43b8-8df2-4ce3376be959",
   "metadata": {},
   "source": [
    "# Reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b4b092b-2120-48c6-a492-84656e8b7169",
   "metadata": {},
   "source": [
    "## Spearman Coefficient\n",
    "- Spearman correlation between predicted and observed expression (on DE genes if available)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e2885bd3-5433-47d6-88d5-47be83df5b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_folds = 5\n",
    "model_adatas = []\n",
    "for n in range(0, n_folds):\n",
    "    # new adata\n",
    "    adata_n = adata.copy()\n",
    "    # load model from checkpoint\n",
    "    s3_dir = \"s3://pert-spectra/scETM_checkpoints/\"\n",
    "    experiment_name = \"scetm_norman/\"\n",
    "    model_name = f\"fold_{n}.h5ad\"\n",
    "    m_adata = read_aws_h5ad(s3_dir + experiment_name + model_name)\n",
    "    model_adatas.append(m_adata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cfe07156-f95b-481e-a1ac-495e69f7e954",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "(18547, 4990)\n",
      "test\n",
      "1\n",
      "(20156, 4990)\n",
      "test\n",
      "2\n",
      "(21410, 4990)\n",
      "test\n",
      "3\n",
      "(20621, 4990)\n",
      "test\n",
      "4\n",
      "(20752, 4990)\n",
      "test\n"
     ]
    }
   ],
   "source": [
    "# iterate through all models and get the losses and correlations\n",
    "# take the mean loss and correlation for the test set\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "test_corr_singles = []\n",
    "test_corr_combos = []\n",
    "\n",
    "\n",
    "for n in range(n_folds):\n",
    "    adata_n = model_adatas[n]\n",
    "    loss_weights = np.ones(adata_n.shape[0])\n",
    "\n",
    "    for pert in adata_n.obs[\"perturbation_name\"].unique():\n",
    "        hold_idx = [\n",
    "            i for i, x in enumerate(adata_n.obs[\"perturbation_name\"]) if x == pert\n",
    "        ]\n",
    "        recon = adata_n[hold_idx].uns[\"recon\"]\n",
    "        # correlation\n",
    "        mean_reconstruction = recon.mean(axis=0)\n",
    "        mean_observed = np.squeeze(np.array(adata_n[hold_idx].X.mean(axis=0)))\n",
    "        if (\"+\" in pert) and (\"ctrl\" not in pert):\n",
    "            test_corr_combos.append(\n",
    "                [pert, spearmanr(mean_reconstruction, mean_observed)[0], n]\n",
    "            )\n",
    "        else:\n",
    "            test_corr_singles.append(\n",
    "                [pert, spearmanr(mean_reconstruction, mean_observed)[0], n]\n",
    "            )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "72f27f1d-a5a5-46f8-a9a4-efcd5f3119fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAHFCAYAAACQIN+BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLSElEQVR4nO3deZyN9f//8edhdrNYhlmYjZgZJFsJyU6hqD6VRJbysX6kRfERM/QxQnyUQipStvZSP1kSolFJUTS0YRQTI1kyhpl5//7oO+fjmBlmzDlzZq553G+3c7s513mf9/W63uea6zxdy7lsxhgjAAAAWEIFdxcAAAAA5yHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAAC7micPftt99q4MCBiomJkY+Pj/z9/dW0aVNNnz5df/zxh7NrLJaNGzfKZrNp48aNRX7v999/r8TERO3fvz/PawMGDFB0dHSx6yuq/fv3y2az6emnn3Zan4cOHVJiYqJ27NjhtD6R1yuvvCKbzeawPrVr107t2rUrcl9JSUl67733nFZbrtz165VXXnF631fK2X9r7dq1U8OGDZ3WX1lRGj/b0ioxMVE2m83dZZQLrv7+yc0Ab731ltP6LI3Z4GJFDncvvviimjVrpm3btmnMmDFavXq13n33Xd15552aP3++7r//flfU6Rbff/+9Jk2alO8HOGHCBL377rslX5QLHDp0SJMmTSLcucHcuXM1d+7cIr/PVeEOgPTAAw9o69at7i6jXCiL3z9lIRt4FKXx1q1bNWzYMHXu3FnvvfeevL297a917txZjzzyiFavXu2Uws6cOSM/P78807Ozs5WVleUwb3eoU6eOW+ePkuPKda5+/fpO7xNlQ2nZlpV2GRkZ8vX1LdF51qpVS7Vq1SrReZak0rDu5dZQVvsvSGnJBkXac5eUlCSbzaYFCxbku1J4eXnp1ltvtT/PycnR9OnTFRcXJ29vb9WoUUP33Xeffv31V4f35R4m+fTTT9WqVSv5+flp0KBB9sMI06dP13/+8x/FxMTI29tbGzZskCR99dVXuvXWW1W1alX5+PioSZMmeuONNy67HF999ZV69+6t6Oho+fr6Kjo6Wvfcc48OHDhgb/PKK6/ozjvvlCS1b99eNpvN4ZBGfrtez549q3HjxikmJkZeXl6qWbOmRowYoT///NOhXXR0tHr06KHVq1eradOm8vX1VVxcnBYuXHjZ2i8c2ylTpigyMlI+Pj5q3ry51q9fn6fdjz/+qD59+qhGjRry9vZWfHy8nn/+efvrGzdu1LXXXitJGjhwoH05ExMT9f/+3/+TzWbTtm3b7O3ffvtt2Ww2de/e3WE+jRo10h133GF/bozR3Llz1bhxY/n6+qpKlSr6xz/+oV9++SVPjR9//LE6duyowMBA+fn5qXXr1nmWJfcwye7du3XPPfcoKChIISEhGjRokE6cOHHZ8cpdxzZv3qzrr79evr6+qlmzpiZMmKDs7Gx7O2etc59//rlat24tHx8fhYeHa9y4cTp//ny+dV18WDYzM1OTJ09WfHy8fHx8VK1aNbVv317JycmSJJvNpr/++kuLFy+2f14X9pGWlqYhQ4aoVq1a8vLyUkxMjCZNmpRnQ3fo0CHdddddCggIUFBQkO6++26lpaVddiwl6ejRoxo+fLjq168vf39/1ahRQx06dNDmzZsd2l14GsGsWbMUExMjf39/tWzZUp9//nmefl955RXFxsba19VXX321UPVIhd/e5LrcuiBJ8+bN0zXXXCN/f38FBAQoLi5O//73vx3aFGa8C1qv3njjDXl5eWnChAl56tuzZ49sNpueffbZIs1LKt5nW5j5GGPUrVs3VatWTampqfb3nTlzRg0aNFB8fLz++usvSf/72/3mm290++23KzAwUEFBQerbt6+OHj3qMN/cbeM777yjJk2ayMfHR5MmTSrSsl/uMztz5oweffRR+2lFVatWVfPmzbV8+XJ7m/wOyxb1+2zbtm1q06aN/Pz8VLt2bT311FPKycm57NjbbDaNHDlSL7zwgurVqydvb2/Vr19fK1asKPLnJF1+m+buGgr6/skdy/xOW7n4+7cwy3j27Fk9/PDDCg0Nla+vr9q2batvvvnGoV8rZAM7U0hZWVnGz8/PtGjRorBvMf/85z+NJDNy5EizevVqM3/+fFO9enUTERFhjh49am/Xtm1bU7VqVRMREWHmzJljNmzYYDZt2mT27dtnJJmaNWua9u3bm7feesusXbvW7Nu3z3zyySfGy8vLtGnTxrz++utm9erVZsCAAUaSWbRokb3vDRs2GElmw4YN9mlvvvmmmThxonn33XfNpk2bzIoVK0zbtm1N9erV7XUdOXLEJCUlGUnm+eefN1u3bjVbt241R44cMcYY079/fxMVFWXvMycnx3Tt2tV4eHiYCRMmmLVr15qnn37aVKpUyTRp0sScPXvW3jYqKsrUqlXL1K9f37z66qtmzZo15s477zSSzKZNmy45prljEhERYW644Qbz9ttvmzfffNNce+21xtPT0yQnJ9vb7t692wQFBZmrr77avPrqq2bt2rXmkUceMRUqVDCJiYnGGGNOnDhhFi1aZCSZJ554wr6cBw8eNKdOnTKenp4mKSnJ3ufQoUONr6+vqVSpkjl37pwxxpjff//d2Gw2M3fuXHu7wYMHG09PT/PII4+Y1atXm2XLlpm4uDgTEhJi0tLS7O1ee+01Y7PZTK9evcw777xjPvjgA9OjRw9TsWJF8/HHH9vbJSQkGEkmNjbWTJw40axbt87MmjXLeHt7m4EDB15yzIz5ex2rVq2aCQ8PN88++6xZs2aNGTVqlJFkRowYkWd8i7PO7d692/j5+Zn69eub5cuXm/fff9907drVREZGGklm3759DnW1bdvW/vz8+fOmffv2xsPDwzz66KNm1apVZuXKlebf//63Wb58uTHGmK1btxpfX1/TrVs3++e1e/duY4wxhw8fNhERESYqKsq88MIL5uOPPzZPPvmk8fb2NgMGDLDP58yZMyY+Pt4EBQWZOXPm2Mcjt8YLlyc/e/bsMcOGDTMrVqwwGzduNB9++KG5//77TYUKFRz+1nLHMzo62tx0003mvffeM++99565+uqrTZUqVcyff/5pb5u7Hvbs2dN88MEHZsmSJeaqq66yL8/lFGV7U5h1Yfny5UaS+de//mXWrl1rPv74YzN//nwzatQoe5vCjvel1qvbbrvNREREmOzsbIfleeyxx4yXl5dJT08v0c+2sPNJT083tWrVMi1atLBvC/r37298fX3Nt99+a2+X+7cbFRVlxowZY9asWWNmzZpl3zbmvteYv7eNYWFhpnbt2mbhwoVmw4YN5ssvvyx0TYX5zIYMGWL8/PzMrFmzzIYNG8yHH35onnrqKTNnzpw8NV+oqOtX3bp1zfz58826devM8OHDjSSzePHiS469Mca+fc/dfqxcudLcdNNNRpJ58803i/w5XWrdc3cNO3fuLPD7J3csL9w+5rr4+/dSy5ibASIiIvJsWwIDA83PP/9s76csZ4M8n2FhG6alpRlJpnfv3oVqn5KSYiSZ4cOHO0z/4osvjCTz73//2z6tbdu2RpJZv369Q9vcD6xOnToOGwBjjImLizNNmjQx58+fd5jeo0cPExYWZt9Q5hfuLpaVlWVOnz5tKlWqZJ555hn79DfffLPA9178Aa5evdpIMtOnT3do9/rrrxtJZsGCBfZpUVFRxsfHxxw4cMA+LSMjw1StWtUMGTKkwDqN+d+YhIeHm4yMDPv0kydPmqpVq5pOnTrZp3Xt2tXUqlXLnDhxwqGPkSNHGh8fH/PHH38YY4zZtm1bgRv9G264wXTo0MH+/KqrrjJjxowxFSpUsK9sS5cuNZLMDz/8YIz5O3xIMjNnznTo6+DBg8bX19c89thjxhhj/vrrL1O1alVzyy23OLTLzs4211xzjbnuuuvs03I3theP7/Dhw42Pj4/Jycm55LjlrmPvv/++w/TBgwebChUq2D8LZ6xzd999t/H19XUIsVlZWSYuLu6y4e7VV181ksyLL754yeWpVKmS6d+/f57pQ4YMMf7+/g7rljHGPP3000aSPQTOmzevwPEoTAC4WFZWljl//rzp2LGjue222+zTc8fz6quvNllZWfbpX375pZFkD6zZ2dkmPDzcNG3a1OGz3L9/v/H09LxsuLuS7c3l1oWRI0eaypUrX3K+hR3vS61XK1euNJLM2rVr7dOysrJMeHi4ueOOO4o8r+J+toWdjzHGbNmyxXh4eJjRo0ebhQsXGknmpZdecnhf7t/uQw895DA9d7uxZMkS+7SoqChTsWJFs3fv3iuqqTCfWcOGDU2vXr0u2ebicHcl69cXX3zh0LZ+/fqma9eul5yvMX8Hq4K2H1dddZV9mjPWvdJQw6W+f4oa7vLrPzcDFLRteeCBBwoch7KUDS7msp9Cyd0dOmDAAIfp1113neLj4/McdqtSpYo6dOiQb1+33nqrPD097c9/+ukn7dmzR/fee68kKSsry/7o1q2bDh8+rL179xZY2+nTp/X444/rqquukoeHhzw8POTv76+//vpLKSkpV7K4+uSTTyTlXd4777xTlSpVyrO8jRs3VmRkpP25j4+P6tWr57D791Juv/12+fj42J8HBATolltu0aeffqrs7GydPXtW69ev12233SY/P788Y3T27Nl8D4tdrGPHjvrss8+UkZGhAwcO6KefflLv3r3VuHFjrVu3TtLfh1UjIyNVt25dSdKHH34om82mvn37Osw3NDRU11xzjf3K5eTkZP3xxx/q37+/Q7ucnBzddNNN2rZtm/3QTq4LD/tLfx8OPnv2rI4cOXLZZQkICMjz/j59+ignJ0effvppnvlc6Tq3YcMGdezYUSEhIfb3V6xYUXffffdla/zoo4/k4+OjQYMGXbZtfj788EO1b99e4eHhDjXefPPNkqRNmzbZayxoPApr/vz5atq0qXx8fOTh4SFPT0+tX78+37+h7t27q2LFivbnjRo1kiT7+r53714dOnRIffr0cTgcFhUVpVatWl22lqJubwqzLlx33XX6888/dc899+j9999Xenp6nvkWdrxzXbxeSdLNN9+s0NBQLVq0yD5tzZo1OnTokMN6UFKfbVGWqXXr1poyZYpmz56tYcOGqW/fvgVeVJf7t5PrrrvukoeHR57Dg40aNVK9evWuqKbCfGbXXXedPvroI40dO1YbN25URkbGZcekqOtXaGiorrvuujzLVdjte0Hbj59++sl+GNgZ615pr6GoLtV/QduWC9c/K2SDXIUOd8HBwfLz89O+ffsK1f7YsWOSpLCwsDyvhYeH21/PlV+7gl77/fffJUmPPvqoPD09HR7Dhw+XpHz/qHP16dNHzz33nB544AGtWbNGX375pbZt26bq1asX6g89P8eOHZOHh4eqV6/uMN1msyk0NDTP8larVi1PH97e3oWef2hoaL7Tzp07p9OnT+vYsWPKysrSnDlz8oxRt27dJF16jHJ16tRJmZmZ2rJli9atW6fg4GA1adJEnTp10scffyxJWr9+vTp16mR/z++//y5jjEJCQvLM+/PPP7fPN/dz/Mc//pGn3bRp02SMyfPTOhePW+65n4UZtws3VBeOmaTLro9FWeeOHTtW4OdzOUePHlV4eLgqVLiy/3f9/vvv+uCDD/LU2KBBgzw1Xmo8LmfWrFkaNmyYWrRoobfffluff/65tm3bpptuuinfz+Jyn1vu+F/puBV1e1OYdaFfv35auHChDhw4oDvuuEM1atRQixYt7P+pkQo/3rnyq8/Dw0P9+vXTu+++az8H55VXXlFYWJi6du1a5HkV97Mt6jLde++98vLyUmZmpsaMGVNgvxfP38PDQ9WqVSvUd0FhayrMZ/bss8/q8ccf13vvvaf27duratWq6tWrl3788ccCay/q+uWq7fuFtThj3SvtNRTVpfovaHku/OyskA1yFfpq2YoVK6pjx4766KOP9Ouvv172SqLcAg8fPpyn7aFDhxQcHOww7VK/KXTxa7nvHTdunG6//fZ83xMbG5vv9BMnTujDDz9UQkKCxo4da5+emZlZrN/oq1atmrKysnT06FGHD9EYo7S0NPtJo86S38nRaWlp8vLykr+/vzw9PVWxYkX169dPI0aMyLePmJiYy86nRYsW8vf318cff6z9+/erY8eOstls6tixo2bOnKlt27YpNTXVIdwFBwfLZrNp8+bN+V54kzst93OcM2eOrr/++nznn9+X1JXKDWgXyh3Hi/+girPOVatWrcDP53KqV6+uLVu2KCcn54oCXnBwsBo1aqQpU6bk+3p4eLi9xi+//PKKapSkJUuWqF27dpo3b57D9FOnThWxYtnrKWj+hampqNubwq4LAwcO1MCBA/XXX3/p008/VUJCgnr06KEffvhBUVFRhR7vXAVt5wYOHKgZM2ZoxYoVuvvuu7Vy5UqNHj3aYW9nSX22RVmm7Oxs3XvvvapSpYq8vb11//3367PPPpOXl1e+869Zs6b9eVZWlo4dO3bZv72i1nS5z6xSpUqaNGmSJk2apN9//92+F++WW27Rnj178u2/qOtXcV3q7yC3Fmete6W5Bh8fn3wvmCtox8Sl+i9oeXKXxSrZIFeRfgpl3LhxWrVqlQYPHqz3338/zx/w+fPntXr1at1yyy32Q6xLlixxKH7btm1KSUnR+PHjr7jo2NhY1a1bVzt37lRSUlKR3muz2WSMyRM6XnrppTxXyhVlr1DHjh01ffp0LVmyRA899JB9+ttvv62//vpLHTt2LFKdl/POO+9oxowZ9kOzp06d0gcffKA2bdqoYsWK8vPzU/v27fXNN9+oUaNG+W5sc11qOT09PXXjjTdq3bp1OnjwoJ566ilJUps2beTh4aEnnnjCHvZy9ejRQ0899ZR+++033XXXXQXOt3Xr1qpcubK+//57jRw58orGoShOnTqllStXOhyuWrZsmSpUqKAbb7zxku8tyjrXvn17rVy5Ur///rs9nGZnZ+v111+/bI0333yzli9frldeeeWSh2YL+p9cjx49tGrVKtWpU0dVqlS5ZI1vvPFGvuNRGDabLc/f0LfffqutW7cqIiKiUH1cKDY2VmFhYVq+fLkefvhh+0b6wIEDSk5OzvMlcbGibm+Kui5UqlRJN998s86dO6devXpp9+7dioqKKvR4X058fLxatGihRYsWKTs7W5mZmRo4cKBDm5L6bIuyTAkJCdq8ebPWrl2rSpUq6cYbb9SYMWP0zDPP5Gm7dOlSNWvWzP78jTfeUFZWVqF+xPtKxrmgz+xCISEhGjBggHbu3KnZs2cX+BNcrvw+y8/69evz3X7UqVPHHi6dte65u4ZLff9ER0frzTffVGZmpr3dsWPHlJycrMDAwCLNp6Bty3333SfJOtkgV5HCXcuWLTVv3jwNHz5czZo107Bhw9SgQQOdP39e33zzjRYsWKCGDRvqlltuUWxsrP75z39qzpw5qlChgm6++Wbt379fEyZMUEREhMNCXokXXnhBN998s7p27aoBAwaoZs2a+uOPP5SSkqKvv/5ab775Zr7vCwwM1I033qgZM2YoODhY0dHR2rRpk15++WVVrlzZoW3ur9gvWLBAAQEB8vHxUUxMTL67TTt37qyuXbvq8ccf18mTJ9W6dWt9++23SkhIUJMmTdSvX79iLe/FKlasqM6dO+vhhx9WTk6Opk2bppMnT9p/NkCSnnnmGd1www1q06aNhg0bpujoaJ06dUo//fSTPvjgA/u5AHXq1JGvr6+WLl2q+Ph4+fv7Kzw83P6F2rFjRz3yyCOSZN9D5+vrq1atWmnt2rVq1KiRatSoYZ9v69at9c9//lMDBw7UV199pRtvvFGVKlXS4cOHtWXLFl199dUaNmyY/P39NWfOHPXv319//PGH/vGPf6hGjRo6evSodu7cqaNHj+bZM1Qc1apV07Bhw5Samqp69epp1apVevHFFzVs2DCHcxwKUth17oknntDKlSvVoUMHTZw4UX5+fnr++efznD+Yn3vuuUeLFi3S0KFDtXfvXrVv3145OTn64osvFB8fr969e0uSrr76am3cuFEffPCBwsLCFBAQoNjYWE2ePFnr1q1Tq1atNGrUKMXGxurs2bPav3+/Vq1apfnz56tWrVq677779N///lf33XefpkyZorp162rVqlVas2ZNocayR48eevLJJ5WQkKC2bdtq7969mjx5smJiYq7ot6UqVKigJ598Ug888IBuu+02DR48WH/++acSExMLdTixqNubwqwLgwcPlq+vr1q3bq2wsDClpaVp6tSpCgoKsn/BF3a8C2PQoEEaMmSIDh06pFatWuU5+lBSn21h57Nu3TpNnTpVEyZMsH9BTZ06VY8++qjatWun2267zaHfd955Rx4eHurcubN2796tCRMm6JprrrnkfwCLWlNhPrMWLVqoR48eatSokapUqaKUlBS99tpratmyZb7BTir6+lVcwcHB6tChgyZMmKBKlSpp7ty52rNnj8NPkThz3XNnDZf6/unXr59eeOEF9e3bV4MHD9axY8c0ffr0Igc7STpy5Ih923LixAklJCTIx8dH48aNk2SdbGBXpMsv/s+OHTtM//79TWRkpPHy8rJf0jtx4kT75cDG/H0F3LRp00y9evWMp6enCQ4ONn379rVf5pyrbdu2pkGDBnnmk3sFzIwZM/KtY+fOneauu+4yNWrUMJ6eniY0NNR06NDBzJ8/394mv6tlf/31V3PHHXeYKlWqmICAAHPTTTeZXbt2maioqDxXIM6ePdvExMSYihUrOlzRc/EVMcb8fVXL448/bqKiooynp6cJCwszw4YNM8ePH3doFxUVZbp3755neQq6Mii/MZk2bZqZNGmSqVWrlvHy8jJNmjQxa9asybf9oEGDTM2aNY2np6epXr26adWqlfnPf/7j0G758uUmLi7OeHp6GkkmISHB/trOnTuNJFO3bl2H90yZMsVIMg8//HC+tS5cuNC0aNHCVKpUyfj6+po6deqY++67z3z11VcO7TZt2mS6d+9uqlatajw9PU3NmjVN9+7dHS65z7167cKfHDDmfz+fcanL+o353zq2ceNG07x5c+Pt7W3CwsLMv//9b4erX52xzhljzGeffWauv/564+3tbUJDQ82YMWPMggULLnu1rDF/r0cTJ040devWNV5eXqZatWqmQ4cODj9zs2PHDtO6dWvj5+dnJDn0cfToUTNq1CgTExNjPD09TdWqVU2zZs3M+PHjzenTp+3tcv8O/P39TUBAgLnjjjtMcnJyoa6ozMzMNI8++qipWbOm8fHxMU2bNjXvvfdegVex5TeeF69nxhjz0ksv2Ze7Xr16ZuHChfn+reWnqNuby60LixcvNu3btzchISHGy8vLhIeHm7vuusvhZz6MKdx4X269MubvnyXy9fW95NXSJfHZFmY+hw4dMjVq1DAdOnRw+AmXnJwcc8stt5jKlSvb1/Pcv93t27ebW265xV7TPffcY37//XeH+Ra0bSzsshfmMxs7dqxp3ry5qVKlivH29ja1a9c2Dz30kP0nZy6s+ULF/T4r7Hqs//tJnrlz55o6deoYT09PExcXZ5YuXXpFY1KYdc/dNVzq+2fx4sUmPj7e+Pj4mPr165vXX3+9SNuZ3Azw2muvmVGjRpnq1asbb29v06ZNmzzfRWU5G1zMZowxromNQOnRrl07paena9euXe4uBShXEhMTNWnSJB09etTp56ZZkc1m04gRI/Tcc8+V6xpQPC77KRQAAACUPMIdAACAhXBYFgAAwELYcwcAAGAhhDsAAAALIdwBAABYSJF+xNjqcnJydOjQIQUEBBT5NikAAMA9jDE6depUse7NbSWEuwscOnToim6dBAAA3O/gwYPFujOHVRDuLhAQECDp75XjSm5vAgAASt7JkycVERFh/x4v7wh3F8g9FBsYGEi4AwCgjOGUqr9xYBoAAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCEe7i4AAOAcqampSk9Pd0pfwcHBioyMdEpfAEoW4Q4ALCA1NVVxcfHKyDjjlP58ff20Z08KAQ8ogwh3AGAB6enpysg4oxaDEhQYFl2svk4e3q8vFk5Seno64Q4ogwh3AGAhgWHRqhoZ6+4yALgRF1QAAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACzEw90FAEB5lpqaqvT09GL3k5KS4oRqAFgB4Q4A3CQ1NVVxcfHKyDjjtD7PZ55zWl8AyibCHQC4SXp6ujIyzqjFoAQFhkUXq6/D323VrpULlJWV5ZziAJRZhDsAcLPAsGhVjYwtVh8nD+93TjEAyjwuqAAAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFeLi7AABA6ZSSkuKUfoKDgxUZGemUvgBcHuEOAOAg48QxSTb17dvXKf35+vppz54UAh5QQgh3AAAH58+ckmTUuM/jqh4TV6y+Th7ery8WTlJ6ejrhDighhDsAQL78a0SqamSsu8sAUERl4oKKrKwsPfHEE4qJiZGvr69q166tyZMnKycnx97GGKPExESFh4fL19dX7dq10+7du91YNQAAQMkrE+Fu2rRpmj9/vp577jmlpKRo+vTpmjFjhubMmWNvM336dM2aNUvPPfectm3bptDQUHXu3FmnTp1yY+UAAAAlq0yEu61bt6pnz57q3r27oqOj9Y9//ENdunTRV199JenvvXazZ8/W+PHjdfvtt6thw4ZavHixzpw5o2XLlrm5egAAgJJTJsLdDTfcoPXr1+uHH36QJO3cuVNbtmxRt27dJEn79u1TWlqaunTpYn+Pt7e32rZtq+Tk5AL7zczM1MmTJx0eAAAAZVmZuKDi8ccf14kTJxQXF6eKFSsqOztbU6ZM0T333CNJSktLkySFhIQ4vC8kJEQHDhwosN+pU6dq0qRJriscAACghJWJPXevv/66lixZomXLlunrr7/W4sWL9fTTT2vx4sUO7Ww2m8NzY0yeaRcaN26cTpw4YX8cPHjQJfUDAACUlDKx527MmDEaO3asevfuLUm6+uqrdeDAAU2dOlX9+/dXaGiopL/34IWFhdnfd+TIkTx78y7k7e0tb29v1xYPAABQgsrEnrszZ86oQgXHUitWrGj/KZSYmBiFhoZq3bp19tfPnTunTZs2qVWrViVaKwAAgDuViT13t9xyi6ZMmaLIyEg1aNBA33zzjWbNmqVBgwZJ+vtw7OjRo5WUlKS6deuqbt26SkpKkp+fn/r06ePm6gEAAEpOmQh3c+bM0YQJEzR8+HAdOXJE4eHhGjJkiCZOnGhv89hjjykjI0PDhw/X8ePH1aJFC61du1YBAQFurBwAAKBklYlwFxAQoNmzZ2v27NkFtrHZbEpMTFRiYmKJ1QUAAFDalIlz7gAAAFA4hDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIV4uLsAAChrUlNTlZ6eXux+UlJSnFANADgi3AFAEaSmpiouLl4ZGWec1uf5zHNO6wsACHcAUATp6enKyDijFoMSFBgWXay+Dn+3VbtWLlBWVpZzigMAEe4A4IoEhkWramRssfo4eXi/c4oBgAtwQQUAAICFEO4AAAAshHAHAABgIYQ7AAAACykz4e63335T3759Va1aNfn5+alx48bavn27/XVjjBITExUeHi5fX1+1a9dOu3fvdmPFAAAAJa9MhLvjx4+rdevW8vT01EcffaTvv/9eM2fOVOXKle1tpk+frlmzZum5557Ttm3bFBoaqs6dO+vUqVPuKxwAAKCElYmfQpk2bZoiIiK0aNEi+7To6Gj7v40xmj17tsaPH6/bb79dkrR48WKFhIRo2bJlGjJkSEmXDAAA4BZlYs/dypUr1bx5c915552qUaOGmjRpohdffNH++r59+5SWlqYuXbrYp3l7e6tt27ZKTk4usN/MzEydPHnS4QEAAFCWlYlw98svv2jevHmqW7eu1qxZo6FDh2rUqFF69dVXJUlpaWmSpJCQEIf3hYSE2F/Lz9SpUxUUFGR/REREuG4hAAAASkCZCHc5OTlq2rSpkpKS1KRJEw0ZMkSDBw/WvHnzHNrZbDaH58aYPNMuNG7cOJ04ccL+OHjwoEvqBwAAKCllItyFhYWpfv36DtPi4+OVmpoqSQoNDZWkPHvpjhw5kmdv3oW8vb0VGBjo8AAAACjLykS4a926tfbu3esw7YcfflBUVJQkKSYmRqGhoVq3bp399XPnzmnTpk1q1apVidYKAADgTmXiatmHHnpIrVq1UlJSku666y59+eWXWrBggRYsWCDp78Oxo0ePVlJSkurWrau6desqKSlJfn5+6tOnj5urBwAAKDllItxde+21evfddzVu3DhNnjxZMTExmj17tu699157m8cee0wZGRkaPny4jh8/rhYtWmjt2rUKCAhwY+UAAAAlq0yEO0nq0aOHevToUeDrNptNiYmJSkxMLLmiAAAASpkycc4dAAAACodwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCXBruateurWPHjuWZ/ueff6p27dqunDUAAEC55NJwt3//fmVnZ+eZnpmZqd9++82VswYAACiXPFzR6cqVK+3/XrNmjYKCguzPs7OztX79ekVHR7ti1gAAAOWaS8Jdr169JEk2m039+/d3eM3T01PR0dGaOXOmK2YNAABQrrkk3OXk5EiSYmJitG3bNgUHB7tiNgAAALiIS8Jdrn379rmyewAAAFzEpeFOktavX6/169fryJEj9j16uRYuXOjq2QMAAJQrLg13kyZN0uTJk9W8eXOFhYXJZrO5cnYAAADlnkvD3fz58/XKK6+oX79+rpwNAAAA/o9Lf+fu3LlzatWqlStnAQAAgAu4NNw98MADWrZsmStnAQAAgAu49LDs2bNntWDBAn388cdq1KiRPD09HV6fNWuWK2cPAABQ7rg03H377bdq3LixJGnXrl0Or3FxBQAAgPO5NNxt2LDBld0DAADgIi495w4AAAAly6V77tq3b3/Jw6+ffPKJK2cPAABQ7rg03OWeb5fr/Pnz2rFjh3bt2qX+/fu7ctYAAADlkkvD3X//+998pycmJur06dOunDUAAEC55JZz7vr27ct9ZQEAAFzALeFu69at8vHxccesAQAALM2lh2Vvv/12h+fGGB0+fFhfffWVJkyY4MpZAwAAlEsuDXdBQUEOzytUqKDY2FhNnjxZXbp0ceWsAQAAyiWXhrtFixa5snsAAABcxKXhLtf27duVkpIim82m+vXrq0mTJiUxWwAAgHLHpeHuyJEj6t27tzZu3KjKlSvLGKMTJ06offv2WrFihapXr+7K2QMAAJQ7Lr1a9l//+pdOnjyp3bt3648//tDx48e1a9cunTx5UqNGjXLlrAEAAMoll+65W716tT7++GPFx8fbp9WvX1/PP/88F1QAAAC4gEv33OXk5MjT0zPPdE9PT+Xk5Lhy1gAAAOWSS/fcdejQQQ8++KCWL1+u8PBwSdJvv/2mhx56SB07dnTlrAEAFpOamqr09HSn9BUcHKzIyEin9AWUNi4Nd88995x69uyp6OhoRUREyGazKTU1VVdffbWWLFniylkDACwkNTVVcXHxysg445T+fH39tGdPCgEPluTScBcREaGvv/5a69at0549e2SMUf369dWpUydXzhYAYDHp6enKyDijFoMSFBgWXay+Th7ery8WTlJ6ejrhDpbkknD3ySefaOTIkfr8888VGBiozp07q3PnzpKkEydOqEGDBpo/f77atGnjitkDcDMOn8FVAsOiVTUy1t1lAKWaS8Ld7NmzNXjwYAUGBuZ5LSgoSEOGDNGsWbMId4AFcfgMANzLJeFu586dmjZtWoGvd+nSRU8//bQrZg3AzTh8BgDu5ZJw9/vvv+f7Eyj2mXp46OjRo66YNYBSgsNnAOAeLvmdu5o1a+q7774r8PVvv/1WYWFhrpg1AABAueaScNetWzdNnDhRZ8+ezfNaRkaGEhIS1KNHD1fMGgAAoFxzyWHZJ554Qu+8847q1aunkSNHKjY2VjabTSkpKXr++eeVnZ2t8ePHu2LWAAAA5ZpLwl1ISIiSk5M1bNgwjRs3TsYYSZLNZlPXrl01d+5chYSEuGLWAAAA5ZrLfsQ4KipKq1at0vHjx/XTTz/JGKO6deuqSpUqrpolAKCUSklJcev7gfLEpXeokKQqVaro2muvdfVsAAClUMaJY5Js6tu3r1P6O595zin9AFbm8nAHACi/zp85JcmocZ/HVT0m7or7OfzdVu1auUBZWVnOKw6wKMIdAMDl/GtEFut3D08e3u+8YgCLc8lPoQAAAMA9CHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAspk+Fu6tSpstlsGj16tH2aMUaJiYkKDw+Xr6+v2rVrp927d7uvSAAAADcoc+Fu27ZtWrBggRo1auQwffr06Zo1a5aee+45bdu2TaGhoercubNOnTrlpkoBAABKXpkKd6dPn9a9996rF198UVWqVLFPN8Zo9uzZGj9+vG6//XY1bNhQixcv1pkzZ7Rs2TI3VgwAAFCyylS4GzFihLp3765OnTo5TN+3b5/S0tLUpUsX+zRvb2+1bdtWycnJJV0mAACA25SZe8uuWLFCX3/9tbZt25bntbS0NElSSEiIw/SQkBAdOHCgwD4zMzOVmZlpf37y5EknVQsAAOAeZWLP3cGDB/Xggw9qyZIl8vHxKbCdzWZzeG6MyTPtQlOnTlVQUJD9ERER4bSaAQAA3KFMhLvt27fryJEjatasmTw8POTh4aFNmzbp2WeflYeHh32PXe4evFxHjhzJszfvQuPGjdOJEyfsj4MHD7p0OQAAAFytTByW7dixo7777juHaQMHDlRcXJwef/xx1a5dW6GhoVq3bp2aNGkiSTp37pw2bdqkadOmFdivt7e3vL29XVo7AABASSoT4S4gIEANGzZ0mFapUiVVq1bNPn306NFKSkpS3bp1VbduXSUlJcnPz099+vRxR8kAAABuUSbCXWE89thjysjI0PDhw3X8+HG1aNFCa9euVUBAgLtLAwAAKDFlNtxt3LjR4bnNZlNiYqISExPdUg8AAEBpUGbDHYDyIyUlxSn9BAcHKzIy0il9AUBpRbgDUGplnDgmyaa+ffs6pT9fXz/t2ZNCwANgaYQ7AKXW+TOnJBk17vO4qsfEFauvk4f364uFk5Senk64A2BphDsApZ5/jUhVjYx1dxkAUCaUiR8xBgAAQOEQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFsLv3AEoV4p7KzNn3QoNAFyFcAegXHD2rczOZ55zSj8A4GyEOwDlgrNuZXb4u63atXKBsrKynFccADgR4Q5AuVLcW5mdPLzfecUAgAtwQQUAAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAsxMPdBQAoHVJTU5Wenl7sflJSUpxQDQDgShHuACg1NVVxcfHKyDjjtD7PZ55zWl8AgMIj3AFQenq6MjLOqMWgBAWGRRerr8PfbdWulQuUlZXlnOIAAEVCuANgFxgWraqRscXq4+Th/c4pBgBwRbigAgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQ7i0LACiXUlJSnNJPcHCwIiMjndIX4AyEOwBAuZJx4pgkm/r27euU/nx9/bRnTwoBD6UG4Q4AUK6cP3NKklHjPo+rekxcsfo6eXi/vlg4Senp6YQ7lBqEOwBAueRfI1JVI2PdXQbgdFxQAQAAYCGEOwAAAAsh3AEAAFhImQh3U6dO1bXXXquAgADVqFFDvXr10t69ex3aGGOUmJio8PBw+fr6ql27dtq9e7ebKgYAAHCPMhHuNm3apBEjRujzzz/XunXrlJWVpS5duuivv/6yt5k+fbpmzZql5557Ttu2bVNoaKg6d+6sU6dOubFyAACAklUmrpZdvXq1w/NFixapRo0a2r59u2688UYZYzR79myNHz9et99+uyRp8eLFCgkJ0bJlyzRkyBB3lA0AAFDiysSeu4udOHFCklS1alVJ0r59+5SWlqYuXbrY23h7e6tt27ZKTk52S40AAADuUCb23F3IGKOHH35YN9xwgxo2bChJSktLkySFhIQ4tA0JCdGBAwcK7CszM1OZmZn25ydPnnRBxQAAACWnzO25GzlypL799lstX748z2s2m83huTEmz7QLTZ06VUFBQfZHRESE0+sFAAAoSWUq3P3rX//SypUrtWHDBtWqVcs+PTQ0VNL/9uDlOnLkSJ69eRcaN26cTpw4YX8cPHjQNYUDAACUkDIR7owxGjlypN555x198skniomJcXg9JiZGoaGhWrdunX3auXPntGnTJrVq1arAfr29vRUYGOjwAAAAKMvKxDl3I0aM0LJly/T+++8rICDAvocuKChIvr6+stlsGj16tJKSklS3bl3VrVtXSUlJ8vPzU58+fdxcPQAAQMkpE+Fu3rx5kqR27do5TF+0aJEGDBggSXrssceUkZGh4cOH6/jx42rRooXWrl2rgICAEq4WAADAfcpEuDPGXLaNzWZTYmKiEhMTXV8QAABAKVUmzrkDAABA4RDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALMTD3QUAZUVqaqrS09Od0ldwcLAiIyOd0hcAABci3AGFkJqaqri4eGVknHFKf76+ftqzJ4WABwBwOsIdUAjp6enKyDijFoMSFBgWXay+Th7ery8WTlJ6ejrhDgDgdIQ7oAgCw6JVNTLW3WUAAFAgwh1QhjnrPMCUlBQnVAMAKA0Id0AZ5ezzACXpfOY5p/UFAHAPwh1QRjnzPMDD323VrpULlJWV5ZziAABuQ7gDyjhnnAd48vB+5xQDAHA7fsQYAADAQgh3AAAAFkK4AwAAsBDOuYNTcGsuAOWZs35OiO0fnIFwh2Lj1lwAyquME8ck2dS3b1+n9Mf2D85AuEOxcWsuAOXV+TOnJBk17vO4qsfEFasvtn9wFsIdnIZbcwEor/xrRLL9Q6lBuCvnnHGuXGm+dVVpvj1XcfsszeMOAHAfwl055uxz5UrbratK6+25nH2OTmkbdwCAexHuyjFnnStXWm9dVVpvz+Wsc3RK67gDANyLcIdinytX2m9dVVpvz1Xcc3RK+7gDANyDHzEGAACwEMIdAACAhRDuAAAALIRz7koIt+cCAJRFfH+VPYS7EsDtuQAAZRHfX2UT4a4EcHsuAEBZxPdX2US4K0HcngsAUBbx/VW2cEEFAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWwu/cAQBQiqSkpBS7D27zVb4R7gAAKAUyThyTZFPfvn2L3Re3+SrfCHcAAJQC58+ckmTUuM/jqh4Td8X9cJsvEO4AAChF/GtEcqsvFIvlwt3cuXM1Y8YMHT58WA0aNNDs2bPVpk0bd5cFAECJcsa5e87oAyXPUuHu9ddf1+jRozV37ly1bt1aL7zwgm6++WZ9//337JoGAJQLzjx3L9f5zHNO6wuuZ6lwN2vWLN1///164IEHJEmzZ8/WmjVrNG/ePE2dOtXN1QEA4HrOOndPkg5/t1W7Vi5QVlaWc4pDibBMuDt37py2b9+usWPHOkzv0qWLkpOT3VQVAADu4Yxz904e3u+cYlCiLBPu0tPTlZ2drZCQEIfpISEhSktLy/c9mZmZyszMtD8/ceKEJOnkyZNOre306dOSpD8O7FVWZkax+jqZlipJ2r59u73fK7V3716n1HXy8AFJ0onffpSnh61YNZXG5ZOcvIxO6qs01lQe+iqNNZWHvkpjTaW1r9JYk/S/7fvp06ed/j2b258xxqn9llnGIn777TcjySQnJztM/89//mNiY2PzfU9CQoKRxIMHDx48ePCwwOPgwYMlETlKPcvsuQsODlbFihXz7KU7cuRInr15ucaNG6eHH37Y/jwnJ0d//PGHqlWrJputeP9DcZaTJ08qIiJCBw8eVGBgoLvLKXMYv+Jh/IqH8Ssexq94ytP4GWN06tQphYeHu7uUUsEy4c7Ly0vNmjXTunXrdNttt9mnr1u3Tj179sz3Pd7e3vL29naYVrlyZVeWecUCAwMt/8fpSoxf8TB+xcP4FQ/jVzzlZfyCgoLcXUKpYZlwJ0kPP/yw+vXrp+bNm6tly5ZasGCBUlNTNXToUHeXBgAAUCIsFe7uvvtuHTt2TJMnT9bhw4fVsGFDrVq1SlFRUe4uDQAAoERYKtxJ0vDhwzV8+HB3l+E03t7eSkhIyHP4GIXD+BUP41c8jF/xMH7Fw/iVXzZjuG4YAADAKiq4uwAAAAA4D+EOAADAQgh3AAAAFkK4AwAAsBDCnRvMnTtXMTEx8vHxUbNmzbR58+ZLtl+6dKmuueYa+fn5KSwsTAMHDtSxY8cc2syePVuxsbHy9fVVRESEHnroIZ09e9aVi+E2RR2/559/XvHx8fL19VVsbKxeffXVPG3efvtt1a9fX97e3qpfv77effddV5Xvds4evxdffFFt2rRRlSpVVKVKFXXq1ElffvmlKxfBrVyx/uVasWKFbDabevXq5eSqSw9XjN+ff/6pESNGKCwsTD4+PoqPj9eqVatctQhu5YrxK0/fH+WGm29/Vu6sWLHCeHp6mhdffNF8//335sEHHzSVKlUyBw4cyLf95s2bTYUKFcwzzzxjfvnlF7N582bToEED06tXL3ubJUuWGG9vb7N06VKzb98+s2bNGhMWFmZGjx5dUotVYoo6fnPnzjUBAQFmxYoV5ueffzbLly83/v7+ZuXKlfY2ycnJpmLFiiYpKcmkpKSYpKQk4+HhYT7//POSWqwS44rx69Onj3n++efNN998Y1JSUszAgQNNUFCQ+fXXX0tqsUqMK8Yv1/79+03NmjVNmzZtTM+ePV28JO7hivHLzMw0zZs3N926dTNbtmwx+/fvN5s3bzY7duwoqcUqMa4Yv/L0/VGeEO5K2HXXXWeGDh3qMC0uLs6MHTs23/YzZswwtWvXdpj27LPPmlq1atmfjxgxwnTo0MGhzcMPP2xuuOEGJ1VdehR1/Fq2bGkeffRRh2kPPvigad26tf35XXfdZW666SaHNl27djW9e/d2UtWlhyvG72JZWVkmICDALF68uPgFlzKuGr+srCzTunVr89JLL5n+/ftbNty5YvzmzZtnateubc6dO+f8gksZV4xfefr+KE84LFuCzp07p+3bt6tLly4O07t06aLk5OR839OqVSv9+uuvWrVqlYwx+v333/XWW2+pe/fu9jY33HCDtm/fbj8U9ssvv2jVqlUObazgSsYvMzNTPj4+DtN8fX315Zdf6vz585KkrVu35umza9euBfZZVrlq/C525swZnT9/XlWrVnVO4aWEK8dv8uTJql69uu6//37nF15KuGr8Vq5cqZYtW2rEiBEKCQlRw4YNlZSUpOzsbNcsiJu4avzKy/dHeUO4K0Hp6enKzs5WSEiIw/SQkBClpaXl+55WrVpp6dKluvvuu+Xl5aXQ0FBVrlxZc+bMsbfp3bu3nnzySd1www3y9PRUnTp11L59e40dO9aly1PSrmT8unbtqpdeeknbt2+XMUZfffWVFi5cqPPnzys9PV2SlJaWVqQ+yypXjd/Fxo4dq5o1a6pTp05OXwZ3ctX4ffbZZ3r55Zf14osvunwZ3MlV4/fLL7/orbfeUnZ2tlatWqUnnnhCM2fO1JQpU1y+TCXJVeNXXr4/yhvCnRvYbDaH58aYPNNyff/99xo1apQmTpyo7du3a/Xq1dq3b5+GDh1qb7Nx40ZNmTJFc+fO1ddff6133nlHH374oZ588kmXLoe7FGX8JkyYoJtvvlnXX3+9PD091bNnTw0YMECSVLFixSvqs6xzxfjlmj59upYvX6533nknzx4Dq3Dm+J06dUp9+/bViy++qODgYFeXXio4e/3LyclRjRo1tGDBAjVr1ky9e/fW+PHjNW/ePJcuh7s4e/zK2/dHeUG4K0HBwcGqWLFinv9lHTlyJM//xnJNnTpVrVu31pgxY9SoUSN17dpVc+fO1cKFC3X48GFJf/8B9+vXTw888ICuvvpq3XbbbUpKStLUqVOVk5Pj8uUqKVcyfr6+vlq4cKHOnDmj/fv3KzU1VdHR0QoICLB/mYaGhhapz7LKVeOX6+mnn1ZSUpLWrl2rRo0auWw53MUV4/fzzz9r//79uuWWW+Th4SEPDw+9+uqrWrlypTw8PPTzzz+XxKKVCFetf2FhYapXr57Dfzbi4+OVlpamc+fOuW6BSpirxq+8fH+UN4S7EuTl5aVmzZpp3bp1DtPXrVunVq1a5fueM2fOqEIFx48pdyNm/u+2wAW1MX9fMOOs8t3uSsYvl6enp2rVqqWKFStqxYoV6tGjh33MWrZsmafPtWvXXrbPssZV4ydJM2bM0JNPPqnVq1erefPmLqnf3VwxfnFxcfruu++0Y8cO++PWW29V+/bttWPHDkVERLhykUqUq9a/1q1b66effnIIIj/88IPCwsLk5eXl/AVxE1eNX3n5/ih3SvoKjvIu91L2l19+2Xz//fdm9OjRplKlSmb//v3GGGPGjh1r+vXrZ2+/aNEi4+HhYebOnWt+/vlns2XLFtO8eXNz3XXX2dskJCSYgIAAs3z5cvPLL7+YtWvXmjp16pi77rqrxJfP1Yo6fnv37jWvvfaa+eGHH8wXX3xh7r77blO1alWzb98+e5vPPvvMVKxY0Tz11FMmJSXFPPXUU5b/KRRnjt+0adOMl5eXeeutt8zhw4ftj1OnTpX04rmcK8bvYla+WtYV45eammr8/f3NyJEjzd69e82HH35oatSoYf7zn/+U9OK5nCvGrzx9f5QnhDs3eP75501UVJTx8vIyTZs2NZs2bbK/1r9/f9O2bVuH9s8++6ypX7++8fX1NWFhYebee+91+A2x8+fPm8TERFOnTh3j4+NjIiIizPDhw83x48dLaIlKVlHG7/vvvzeNGzc2vr6+JjAw0PTs2dPs2bMnT59vvvmmiY2NNZ6eniYuLs68/fbbJbEobuHs8YuKijKS8jwSEhJKaIlKlivWvwtZOdwZ45rxS05ONi1atDDe3t6mdu3aZsqUKSYrK6skFqfEOXv8ytv3R3lhM4b9rgAAAFbBOXcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4A3CgxMVGNGzcuNf0AKPsId0AZduTIEQ0ZMkSRkZHy9vZWaGiounbtqq1bt7q7tFLr7bffVrt27RQUFCR/f381atRIkydP1h9//OHu0grNZrPpvffec5j26KOPav369e4pCECpQrgDyrA77rhDO3fu1OLFi/XDDz9o5cqVateuXakIKufPn3d3CXmMHz9ed999t6699lp99NFH2rVrl2bOnKmdO3fqtddeu+J+81vWkl5+f39/VatWrUTnCaCUcvf9zwBcmePHjxtJZuPGjZdsJ8nMnTvX3HTTTcbHx8dER0ebN954w6HNr7/+au666y5TuXJlU7VqVXPrrbc63Fz8yy+/NJ06dTLVqlUzgYGB5sYbbzTbt2/PM5958+aZW2+91fj5+ZmJEyeahIQEc80115iXX37ZREREmEqVKpmhQ4earKwsM23aNBMSEmKqV6+e5ybvM2fONA0bNjR+fn6mVq1aZtiwYebUqVP21xctWmSCgoLM6tWrTVxcnKlUqZLp2rWrOXToUIHj8MUXXxhJZvbs2QWOZ665c+ea2rVrG09PT1OvXj3z6quvFmlZY2JijM1mMzk5OebPP/80gwcPNtWrVzcBAQGmffv2ZseOHfa+ct9X2LG++F6+UVFR+faTnZ1tJk2aZGrWrGm8vLzMNddcYz766CP76/v27TOSzNtvv23atWtnfH19TaNGjUxycnKBYwigbGDPHVBG+fv7y9/fX++9954yMzMv2XbChAn2vXx9+/bVPffco5SUFEnSmTNn1L59e/n7++vTTz/Vli1b5O/vr5tuuknnzp2TJJ06dUr9+/fX5s2b9fnnn6tu3brq1q2bTp065TCfhIQE9ezZU999950GDRokSfr555/10UcfafXq1Vq+fLkWLlyo7t2769dff9WmTZs0bdo0PfHEE/r888/t/VSoUEHPPvusdu3apcWLF+uTTz7RY4895jCvM2fO6Omnn9Zrr72mTz/9VKmpqXr00UcLHIOlS5fK399fw4cPz/f1ypUrS5LeffddPfjgg3rkkUe0a9cuDRkyRAMHDtSGDRsuu6w//fST3njjDb399tvasWOHJKl79+5KS0vTqlWrtH37djVt2lQdO3YscO/q5cZ627ZtkqRFixbp8OHD9ucXe+aZZzRz5kw9/fTT+vbbb9W1a1fdeuut+vHHHx3ajR8/Xo8++qh27NihevXq6Z577lFWVlaB4wigDHB3ugRw5d566y1TpUoV4+PjY1q1amXGjRtndu7c6dBGkhk6dKjDtBYtWphhw4YZY4x5+eWXTWxsrMnJybG/npmZaXx9fc2aNWvynW9WVpYJCAgwH3zwgcN8Ro8e7dAuISHB+Pn5mZMnT9qnde3a1URHR5vs7Gz7tNjYWDN16tQCl/ONN94w1apVsz9ftGiRkWR++ukn+7Tnn3/ehISEFNjHzTffbBo1alTg67latWplBg8e7DDtzjvvNN26dbM/L2hZPT09zZEjR+zT1q9fbwIDA83Zs2cd2tapU8e88MIL9vdduMftYgWN9bvvvptn/hf2Ex4ebqZMmeLQ5tprrzXDhw83xvxvz91LL71kf3337t1GkklJSSmwHgClH3vugDLsjjvu0KFDh7Ry5Up17dpVGzduVNOmTfXKK684tGvZsmWe57l77rZv366ffvpJAQEB9r2BVatW1dmzZ/Xzzz9L+vvCjaFDh6pevXoKCgpSUFCQTp8+rdTUVId+mzdvnqfG6OhoBQQE2J+HhISofv36qlChgsO0I0eO2J9v2LBBnTt3Vs2aNRUQEKD77rtPx44d019//WVv4+fnpzp16tifh4WFOfRxMWOMbDZbga/nSklJUevWrR2mtW7d2j5el1rWqKgoVa9e3f58+/btOn36tKpVq2YfW39/f+3bt88+thcr7FhfysmTJ3Xo0KFCLUejRo3s/w4LC7PXAKDs8nB3AQCKx8fHR507d1bnzp01ceJEPfDAA0pISNCAAQMu+b7coJOTk6NmzZpp6dKledrkBpUBAwbo6NGjmj17tqKiouTt7a2WLVvaD9vmqlSpUp4+PD0988w3v2k5OTmSpAMHDqhbt24aOnSonnzySVWtWlVbtmzR/fff73CRQn59GGMKXN569eppy5YtOn/+fJ73XuziEJhfMMxvWS+elpOTo7CwMG3cuDFP29zDwBcr7FgXRmGW48KxuHCdAFB2secOsJj69es77OGS5HA+W+7zuLg4SVLTpk31448/qkaNGrrqqqscHkFBQZKkzZs3a9SoUerWrZsaNGggb29vpaenu6T+r776SllZWZo5c6auv/561atXT4cOHSp2v3369NHp06c1d+7cfF//888/JUnx8fHasmWLw2vJycmKj48v8jybNm2qtLQ0eXh45Bnb4ODgfN9TmLH29PRUdnZ2gfMNDAxUeHi405YDQNlCuAPKqGPHjqlDhw5asmSJvv32W+3bt09vvvmmpk+frp49ezq0ffPNN7Vw4UL98MMPSkhI0JdffqmRI0dKku69914FBwerZ8+e2rx5s/bt26dNmzbpwQcf1K+//ipJuuqqq/Taa68pJSVFX3zxhe699175+vq6ZLnq1KmjrKwszZkzR7/88otee+01zZ8/v9j9tmjRQo899pgeeeQRPfbYY9q6dasOHDig9evX684779TixYslSWPGjNErr7yi+fPn68cff9SsWbP0zjvvXPJijYJ06tRJLVu2VK9evbRmzRrt379fycnJeuKJJ/TVV1/l+57CjHV0dLTWr1+vtLQ0HT9+PN9+xowZo2nTpun111/X3r17NXbsWO3YsUMPPvhgkZcDQNlCuAPKKH9/f7Vo0UL//e9/deONN6phw4aaMGGCBg8erOeee86h7aRJk7RixQo1atRIixcv1tKlS1W/fn1Jf5+79umnnyoyMlK333674uPjNWjQIGVkZCgwMFCStHDhQh0/flxNmjRRv379NGrUKNWoUcMly9W4cWPNmjVL06ZNU8OGDbV06VJNnTrVKX1PmzZNy5Yt0xdffKGuXbuqQYMGevjhh9WoUSP1799fktSrVy8988wzmjFjhho0aKAXXnhBixYtUrt27Yo8P5vNplWrVunGG2/UoEGDVK9ePfXu3Vv79+9XSEhIvu8pzFjPnDlT69atU0REhJo0aZJvP6NGjdIjjzyiRx55RFdffbVWr16tlStXqm7dukVeDgBli81c6iQVAGWezWbTu+++q169erm7FABACWDPHQAAgIUQ7gAAACyEn0IBLI4zLwCgfGHPHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIX8f+twJi5lWttKAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot train correlation distribution vs test correlation distribution\n",
    "import seaborn as sns\n",
    "\n",
    "test_corr = np.array(test_corr_singles + test_corr_combos)[:, 1].astype(float)\n",
    "sns.histplot(test_corr, label=\"test_corr\")\n",
    "plt.title(\"Correlation between predicted and observed expression per perturbation\")\n",
    "plt.xlabel(\"Spearman Correlation\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e4dcca93-bc7e-4bae-b434-50f2a2fb52ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\n",
    "    \"figures/reconstruction_spearmans/scETM_norman_spearman_correlations.csv\",\n",
    "    test_corr,\n",
    "    delimiter=\",\",\n",
    "    fmt=\"%.2f\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64fc0b35-f87e-4d8a-87f1-e6080a3ee83f",
   "metadata": {},
   "source": [
    "# Perturbation Embedding Heatmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c36ebe2b-2835-44ed-9f71-a52b522125b1",
   "metadata": {},
   "source": [
    "## Hierarchical Clustering + Enrichment\n",
    "- Perform enrichment tests (gprofiler) on hierarchical clustering of perturbation embeddings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aaf8dc3-7e5b-4748-8685-fc8560bf7fc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "from scipy.cluster.hierarchy import dendrogram, fcluster, linkage\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "\n",
    "# aggregate cell embeddings to perturbation embeddings\n",
    "cell_emb = adata.uns[\"cell_emb\"] @ adata.uns[\"topics\"]\n",
    "perts = []\n",
    "pert_emb = []\n",
    "for i in adata.obs[\"perturbation_name\"].unique():\n",
    "    if \"+\" not in i and i != \"nan\":\n",
    "        perts.append(i)\n",
    "        pert_emb.append(cell_emb[adata.obs[\"perturbation_name\"] == i].mean(axis=0))\n",
    "pert_emb = np.array(pert_emb)\n",
    "pert_emb_df = pd.DataFrame(pert_emb, index=perts)\n",
    "\n",
    "# Compute the pairwise distances\n",
    "df = pert_emb_df.drop(index=[\"control\"])\n",
    "\n",
    "distance = \"euclidean\"\n",
    "distances = pdist(df.values, metric=distance)\n",
    "\n",
    "# Convert the distances into a square distance matrix\n",
    "distance_matrix = pd.DataFrame(squareform(distances), index=df.index, columns=df.index)\n",
    "\n",
    "# Compute the linkage matrix using the condensed distance matrix\n",
    "linkage_matrix = linkage(distance_matrix, method=\"ward\")\n",
    "\n",
    "# Plot the dendrogram (optional)\n",
    "dendro = dendrogram(linkage_matrix, no_plot=True)\n",
    "\n",
    "# Assign clusters\n",
    "max_d = 1.2  # Max distance for flat clusters\n",
    "clusters = fcluster(linkage_matrix, max_d, criterion=\"distance\")\n",
    "\n",
    "# Create a color palette for clusters\n",
    "palette = sns.color_palette(\"husl\", len(np.unique(clusters)))\n",
    "\n",
    "# Map each cluster id to a color\n",
    "cluster_colors = [palette[i - 1] for i in clusters]\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "clustermap = sns.clustermap(\n",
    "    distance_matrix, cmap=\"viridis_r\", row_colors=cluster_colors\n",
    ")\n",
    "clustermap.fig.suptitle(f\"Pairwise {distance} distance of perturbation latent vectors\")\n",
    "# clustermap.ax_row_dendrogram.set_visible(False)\n",
    "clustermap.ax_col_dendrogram.set_visible(False)\n",
    "# clustermap.cax.set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "c655f7dc-7d39-43dd-9af5-6ca0faf6d8eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract clusters and perform gprofiler\n",
    "index_names = list(distance_matrix.index)\n",
    "cluster_df = pd.DataFrame({\"Pert\": index_names, \"Cluster\": clusters})\n",
    "cluster_process = {}\n",
    "for cluster_id in np.unique(clusters):\n",
    "    points_in_cluster = cluster_df[cluster_df[\"Cluster\"] == cluster_id][\"Pert\"].tolist()\n",
    "    gprofiler_in = pd.DataFrame(points_in_cluster, columns=[\"gene_symbol\"])\n",
    "    res = get_gprofiler(gprofiler_in)\n",
    "    cluster_process[cluster_id] = res[res[\"p_value\"] <= 0.05]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e41dff85-043e-46d2-8737-4db7d7bf969c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create descriptions for each cluster (for demonstration purpose)\n",
    "cluster_descriptions = {\n",
    "    1: \"Protein tyrosine phosphatase activity\",\n",
    "    2: \"\",\n",
    "    3: \"protein serine/threonine kinase activity\",\n",
    "    4: \"DNA-binding transcription factor activity, RNA polymerase II transcription\",\n",
    "    5: \"MAPK activity, regulation of IRE1-mediated unfolded protein\",\n",
    "    6: \"morphogenesis, DNA-binding transcription factor activity, transcription by RNA polymerase II\",\n",
    "    7: \"negative regulation of ERBB signaling pathway\",\n",
    "    8: \"MAPK activity, procentriole replication\",\n",
    "    9: \"vesicle lumen\",\n",
    "    10: \"collagen trimer\",\n",
    "    11: \"DNA-binding transcription factor activity, RNA polymerase II transcription\",\n",
    "    12: \"\",\n",
    "    13: \"\",\n",
    "    14: \"C/EBP complex\",\n",
    "}\n",
    "\n",
    "# Create a consistent color palette\n",
    "unique_clusters = sorted(np.unique(clusters))\n",
    "palette = sns.color_palette(\"husl\", len(unique_clusters))\n",
    "cluster_colors = [palette[i - 1] for i in clusters]\n",
    "cluster_colors = {cid: palette[i] for i, cid in enumerate(unique_clusters)}\n",
    "\n",
    "# Create clustermap using consistent colors\n",
    "row_colors = [cluster_colors[cid] for cid in clusters]\n",
    "clustermap = sns.clustermap(\n",
    "    distance_matrix, cmap=\"viridis_r\", row_colors=[palette[i - 1] for i in clusters]\n",
    ")\n",
    "clustermap.savefig(\n",
    "    \"figures/figure_pngs/scETM_norman_clustermap.png\", dpi=600, bbox_inches=\"tight\"\n",
    ")\n",
    "\n",
    "# Prepare data for the table\n",
    "cluster_data = {\n",
    "    \"Cluster ID\": unique_clusters,\n",
    "    \"Description\": [cluster_descriptions[cid] for cid in unique_clusters],\n",
    "}\n",
    "\n",
    "cluster_df = pd.DataFrame(cluster_data)\n",
    "\n",
    "# Plot the table with colors matching the clustermap\n",
    "fig, ax = plt.subplots(figsize=(18, 4))\n",
    "ax.axis(\"off\")\n",
    "ax.axis(\"tight\")\n",
    "\n",
    "# Create table while applying the color to the row background\n",
    "table = ax.table(\n",
    "    cellText=cluster_df.values,\n",
    "    colLabels=[\"Cluster ID\", \"Description\"],\n",
    "    cellColours=[[cluster_colors[cid]] * 2 for cid in cluster_df[\"Cluster ID\"]],\n",
    "    cellLoc=\"center\",\n",
    "    loc=\"center\",\n",
    ")\n",
    "\n",
    "# Customize table appearance\n",
    "table.auto_set_font_size(False)\n",
    "table.set_fontsize(10)\n",
    "table.scale(1, 1.5)\n",
    "\n",
    "plt.savefig(\n",
    "    \"figures/figure_pngs/scETM_norman_clustermap_descriptions.png\",\n",
    "    dpi=600,\n",
    "    bbox_inches=\"tight\",\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "a53c2cec-1c3d-465d-a627-47bcc7aa4da7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save grpofiler results\n",
    "with open(\n",
    "    \"figures/pert_embedding_cluster_gprofiler/scETM_norman_pert_emb_gprofiler.pickle\",\n",
    "    \"wb\",\n",
    ") as handle:\n",
    "    pickle.dump(cluster_process, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ad3121b-7c81-435e-9d09-ab431200706a",
   "metadata": {},
   "source": [
    "## PR Curve+AUC\n",
    "- AUCPR using prior graph as binary label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "4509c613-b436-4143-8246-b9d83125b4bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               AUC\n",
      "StringDB  0.726058\n",
      "CORUM     0.946866\n",
      "                F1  Threshold\n",
      "StringDB  0.820312        1.0\n",
      "CORUM     0.972222        1.0\n"
     ]
    }
   ],
   "source": [
    "aucpr_df, f1_df, pr_dict = auprc(distance_matrix)\n",
    "print(aucpr_df)\n",
    "print(f1_df)\n",
    "aucpr_df.to_csv(\"figures/pert_embedding_recall/scETM_norman_aucpr.csv\")\n",
    "f1_df.to_csv(\"figures/pert_embedding_recall/scETM_norman_f1.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fdff66a-8f20-4268-8520-d33a0fa4969a",
   "metadata": {},
   "source": [
    "# Interpretability of our learned latent space"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ec54f80-8a0f-41f0-8c55-edd2b1fbfb3c",
   "metadata": {},
   "source": [
    "## Factor Enrichment\n",
    "- GSEA on log transformed latent-by-gene factors\n",
    "- Then associate each perturbation to its top latent factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35eb8169-15b1-4b5b-b926-d5b7238f372a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "factor_to_go = factor_enrichment_gsea(adata, np.abs(adata.uns[\"gene_emb\"]), fdr=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "814189ce-53f3-4716-834a-a4f0ea6efcf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter and add description to processes\n",
    "go_df = read_aws_csv(\"s3://pert-spectra/references/GO_to_Description.txt\")\n",
    "go_df.set_index(\"Term\", inplace=True)\n",
    "go_dict = go_df.to_dict()[\"Description\"]\n",
    "\n",
    "filtered_factor_to_go = {}\n",
    "for i in factor_to_go:\n",
    "    proc = factor_to_go[i]\n",
    "    proc[\"descr\"] = [go_dict[x] for x in proc[\"GO_ID\"]]\n",
    "    filtered_factor_to_go[i] = proc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d1e2990-edee-430e-957f-328342271a42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save latent enrichment results\n",
    "with open(\n",
    "    \"figures/factor_enrichments/scETM_norman_factor_enrichment.pickle\", \"wb\"\n",
    ") as handle:\n",
    "    pickle.dump(filtered_factor_to_go, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pertspectra",
   "language": "python",
   "name": "pertspectra"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}