norman_GSFA_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b8270f66-7a27-4263-aa97-25abeb2cf34d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
    "    perturbation_signal_recovery,\n",
    "    read_aws_csv,\n",
    "    read_aws_h5ad,\n",
    "    read_aws_npz,\n",
    "    read_aws_pickle,\n",
    "    retrieve_stringdb_neighbors,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "39762015-fa94-4ae8-9038-405f7c56a03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read learned parameters\n",
    "Z = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/norman_gsfa_outputs/Z.csv\"\n",
    ").to_numpy()\n",
    "W = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/norman_gsfa_outputs/W.csv\"\n",
    ").to_numpy()\n",
    "F = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/norman_gsfa_outputs/F.csv\"\n",
    ").to_numpy()\n",
    "beta = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/norman_gsfa_outputs/beta.csv\"\n",
    ").to_numpy()\n",
    "lsfr = read_aws_csv(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/norman_gsfa_outputs/lsfr.csv\"\n",
    ").to_numpy()\n",
    "# read gene and perturbation information\n",
    "gene_labels = read_aws_csv(\"s3://pert-spectra/gsfa_checkpoints/norman_top_genes.csv\")[\n",
    "    \"x\"\n",
    "].to_numpy()\n",
    "pert_labels = read_aws_npz(\"s3://pert-spectra/gsfa_checkpoints/norman_G_labels.npz\")[\n",
    "    \"arr_0\"\n",
    "]\n",
    "pert_labels = [x.split(\"_\")[1] for x in pert_labels[:-1]]\n",
    "pert_labels = pert_labels + [\"ctrl\"]\n",
    "# read inputs and preprocessed input\n",
    "preprocessed_Y = read_aws_npz(\n",
    "    \"s3://pert-spectra/gsfa_checkpoints/norman_GSFA_preprocessed.npz\"\n",
    ")[\"array1\"]\n",
    "G = read_aws_npz(\"s3://pert-spectra/gsfa_checkpoints/norman_GSFA_inputs.npz\")[\"array2\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "601bbfb0-17ed-4694-850c-fcbf1fef961b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in adata of raw data for reference\n",
    "adata = read_aws_h5ad(\"path to raw Norman adata here\")\n",
    "adata = adata[:, gene_labels - 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7888facf-9e1c-4b41-9f48-2e0d194a36af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset to powered perturbations\n",
    "obs_df = pd.DataFrame(adata.obs[\"perturbation_name\"])\n",
    "category_counts = obs_df[\"perturbation_name\"].value_counts()\n",
    "filtered_categories = category_counts[category_counts >= 50].index\n",
    "adata = adata[adata.obs[\"perturbation_name\"].isin(filtered_categories)]"
   ]
  },
  {
   "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": 7,
   "id": "e2885bd3-5433-47d6-88d5-47be83df5b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for GSFA, measure reconstruction on training\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "reconstruction = Z @ W.T\n",
    "# rebuild pert labels\n",
    "G_labels = []\n",
    "for row in G:\n",
    "    pert_idx = np.where(row == 1)[0].tolist()\n",
    "    perts = [pert_labels[i] for i in pert_idx]\n",
    "    G_labels.append(\"+\".join(perts))\n",
    "test_corr_singles = []\n",
    "test_corr_combos = []\n",
    "\n",
    "for pert in set(G_labels):\n",
    "    hold_idx = [i for i, x in enumerate(G_labels) if x == pert]\n",
    "    if not hold_idx:\n",
    "        continue\n",
    "    recon = reconstruction[hold_idx]\n",
    "    # correlation\n",
    "    mean_reconstruction = recon.mean(axis=0)\n",
    "    mean_observed = preprocessed_Y[hold_idx].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",
    "    else:\n",
    "        test_corr_singles.append(\n",
    "            [pert, spearmanr(mean_reconstruction, mean_observed)[0]]\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bc94b654-4864-485b-94e7-c1cb1f6f8baf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAncAAAHFCAYAAACQIN+BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABI9UlEQVR4nO3dd3xUVf7/8fdAehUIpEBIIiWhKkURWHpvAro2FCnKIsgigqCImOCuICiIiojsSrFQZEVEF4GI1A0lIEUwWHaBoBBCEAglBJKc3x/+Ml+GJJBAJkMur+fjMY8Hc+bMPZ9z7507b+7MndiMMUYAAACwhDKuLgAAAADFh3AHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAs5LrC3Z49ezRgwABFRUXJy8tLfn5+atiwoaZMmaLff/+9uGu8IevWrZPNZtO6deuK/NwffvhBcXFxOnjwYJ7H+vfvr8jIyBuur6gOHjwom82mN954o9iWeeTIEcXFxWnXrl3FtkzkNW/ePNlsNof9qXXr1mrdunWRlzVx4kQtW7as2GrLlbt/zZs3r9iXfb2K+7XWunVr1a1bt9iWV1rcjNv2ZhUXFyebzebqMm4Jzn7/yc0A//rXv4ptmTdjNrhSkcPdP/7xDzVq1EiJiYkaPXq0Vq5cqc8//1wPPPCAZs2apSeeeMIZdbrEDz/8oAkTJuS7AcePH6/PP/+85ItygiNHjmjChAmEOxeYOXOmZs6cWeTnOSvcAZCefPJJbd682dVl3BJK4/tPacgGbkXpvHnzZg0ZMkQdOnTQsmXL5OnpaX+sQ4cOGjVqlFauXFkshZ0/f14+Pj552rOzs5WVleUwtitUq1bNpeOj5Dhzn6tdu3axLxOlw81yLLvZZWRkyNvbu0THrFKliqpUqVKiY5akm2Hfy62htC6/IDdLNijSmbuJEyfKZrNp9uzZ+e4UHh4euvfee+33c3JyNGXKFMXExMjT01OVKlXS448/rl9//dXhebkfk2zYsEHNmjWTj4+PBg4caP8YYcqUKfr73/+uqKgoeXp6au3atZKk7du3695771X58uXl5eWlBg0a6NNPP73mPLZv366HH35YkZGR8vb2VmRkpB555BEdOnTI3mfevHl64IEHJElt2rSRzWZz+Egjv1OvFy5c0NixYxUVFSUPDw9VrlxZTz/9tE6dOuXQLzIyUt27d9fKlSvVsGFDeXt7KyYmRnPmzLlm7Zev21dffVVVq1aVl5eXGjdurDVr1uTp9/PPP6tPnz6qVKmSPD09VatWLb377rv2x9etW6e77rpLkjRgwAD7POPi4vTvf/9bNptNiYmJ9v6fffaZbDabunXr5jBO/fr1df/999vvG2M0c+ZM3XnnnfL29la5cuX05z//Wf/73//y1PjNN9+oXbt2CggIkI+Pj5o3b55nLrkfk+zbt0+PPPKIAgMDFRwcrIEDB+r06dPXXF+5+9jGjRt1zz33yNvbW5UrV9b48eOVnZ1t71dc+9yWLVvUvHlzeXl5KSwsTGPHjtWlS5fyrevKj2UzMzP1yiuvqFatWvLy8lKFChXUpk0bJSQkSJJsNpvOnTun+fPn27fX5ctISUnR4MGDVaVKFXl4eCgqKkoTJkzIc6A7cuSIHnzwQfn7+yswMFAPPfSQUlJSrrkuJen48eMaOnSoateuLT8/P1WqVElt27bVxo0bHfpd/jWCadOmKSoqSn5+fmratKm2bNmSZ7nz5s1TdHS0fV/98MMPC1WPVPjjTa5r7QuS9N577+mOO+6Qn5+f/P39FRMToxdffNGhT2HWd0H71aeffioPDw+NHz8+T3379++XzWbT22+/XaSxpBvbtoUZxxijrl27qkKFCkpOTrY/7/z586pTp45q1aqlc+fOSfq/1+7OnTt13333KSAgQIGBgXrsscd0/Phxh3Fzj41Lly5VgwYN5OXlpQkTJhRp7tfaZufPn9dzzz1n/1pR+fLl1bhxYy1cuNDeJ7+PZYv6fpaYmKgWLVrIx8dHt99+u1577TXl5ORcc93bbDYNGzZM77//vmrWrClPT0/Vrl1bixYtKvJ2kq59THN1DQW9/+Suy/y+tnLl+29h5njhwgWNHDlSISEh8vb2VqtWrbRz506H5VohG9iZQsrKyjI+Pj6mSZMmhX2K+ctf/mIkmWHDhpmVK1eaWbNmmYoVK5rw8HBz/Phxe79WrVqZ8uXLm/DwcPPOO++YtWvXmvXr15sDBw4YSaZy5cqmTZs25l//+pdZvXq1OXDggPn222+Nh4eHadGihVm8eLFZuXKl6d+/v5Fk5s6da1/22rVrjSSzdu1ae9uSJUvMyy+/bD7//HOzfv16s2jRItOqVStTsWJFe12pqalm4sSJRpJ59913zebNm83mzZtNamqqMcaYfv36mYiICPsyc3JyTKdOnYybm5sZP368Wb16tXnjjTeMr6+vadCggblw4YK9b0REhKlSpYqpXbu2+fDDD82qVavMAw88YCSZ9evXX3Wd5q6T8PBw86c//cl89tlnZsmSJeauu+4y7u7uJiEhwd533759JjAw0NSrV898+OGHZvXq1WbUqFGmTJkyJi4uzhhjzOnTp83cuXONJPPSSy/Z53n48GFz5swZ4+7ubiZOnGhf5lNPPWW8vb2Nr6+vuXjxojHGmGPHjhmbzWZmzpxp7zdo0CDj7u5uRo0aZVauXGkWLFhgYmJiTHBwsElJSbH3++ijj4zNZjO9evUyS5cuNV9++aXp3r27KVu2rPnmm2/s/WJjY40kEx0dbV5++WUTHx9vpk2bZjw9Pc2AAQOuus6M+WMfq1ChggkLCzNvv/22WbVqlRk+fLiRZJ5++uk86/dG9rl9+/YZHx8fU7t2bbNw4ULzxRdfmE6dOpmqVasaSebAgQMOdbVq1cp+/9KlS6ZNmzbGzc3NPPfcc2bFihVm+fLl5sUXXzQLFy40xhizefNm4+3tbbp27WrfXvv27TPGGHP06FETHh5uIiIizPvvv2+++eYb87e//c14enqa/v3728c5f/68qVWrlgkMDDTvvPOOfX3k1nj5fPKzf/9+M2TIELNo0SKzbt0689VXX5knnnjClClTxuG1lrs+IyMjTefOnc2yZcvMsmXLTL169Uy5cuXMqVOn7H1z98OePXuaL7/80nz88cemevXq9vlcS1GON4XZFxYuXGgkmb/+9a9m9erV5ptvvjGzZs0yw4cPt/cp7Pq+2n7Vu3dvEx4ebrKzsx3mM2bMGOPh4WHS0tJKdNsWdpy0tDRTpUoV06RJE/uxoF+/fsbb29vs2bPH3i/3tRsREWFGjx5tVq1aZaZNm2Y/NuY+15g/jo2hoaHm9ttvN3PmzDFr164127ZtK3RNhdlmgwcPNj4+PmbatGlm7dq15quvvjKvvfaaeeedd/LUfLmi7l81atQws2bNMvHx8Wbo0KFGkpk/f/5V170xxn58zz1+LF++3HTu3NlIMkuWLCnydrravufqGnbv3l3g+0/uurz8+Jjryvffq80xNwOEh4fnObYEBASY//73v/bllOZskGcbFrZjSkqKkWQefvjhQvVPSkoykszQoUMd2rdu3WokmRdffNHe1qpVKyPJrFmzxqFv7garVq2awwHAGGNiYmJMgwYNzKVLlxzau3fvbkJDQ+0HyvzC3ZWysrLM2bNnja+vr3nrrbfs7UuWLCnwuVduwJUrVxpJZsqUKQ79Fi9ebCSZ2bNn29siIiKMl5eXOXTokL0tIyPDlC9f3gwePLjAOo35v3USFhZmMjIy7O3p6emmfPnypn379va2Tp06mSpVqpjTp087LGPYsGHGy8vL/P7778YYYxITEws86P/pT38ybdu2td+vXr26GT16tClTpox9Z/vkk0+MJPPTTz8ZY/4IH5LM1KlTHZZ1+PBh4+3tbcaMGWOMMebcuXOmfPnypkePHg79srOzzR133GHuvvtue1vuwfbK9Tt06FDj5eVlcnJyrrrecvexL774wqF90KBBpkyZMvZtURz73EMPPWS8vb0dQmxWVpaJiYm5Zrj78MMPjSTzj3/846rz8fX1Nf369cvTPnjwYOPn5+ewbxljzBtvvGEk2UPge++9V+D6KEwAuFJWVpa5dOmSadeunendu7e9PXd91qtXz2RlZdnbt23bZiTZA2t2drYJCwszDRs2dNiWBw8eNO7u7tcMd9dzvLnWvjBs2DBz2223XXXcwq7vq+1Xy5cvN5LM6tWr7W1ZWVkmLCzM3H///UUe60a3bWHHMcaYTZs2GTc3NzNixAgzZ84cI8n885//dHhe7mv32WefdWjPPW58/PHH9raIiAhTtmxZ8+OPP15XTYXZZnXr1jW9evW6ap8rw9317F9bt2516Fu7dm3TqVOnq45rzB/BqqDjR/Xq1e1txbHv3Qw1XO39p6jhLr/l52aAgo4tTz75ZIHroTRlgys57adQck+H9u/f36H97rvvVq1atfJ87FauXDm1bds232Xde++9cnd3t9//5ZdftH//fj366KOSpKysLPuta9euOnr0qH788ccCazt79qyef/55Va9eXW5ubnJzc5Ofn5/OnTunpKSk65muvv32W0l55/vAAw/I19c3z3zvvPNOVa1a1X7fy8tLNWvWdDj9ezX33XefvLy87Pf9/f3Vo0cPbdiwQdnZ2bpw4YLWrFmj3r17y8fHJ886unDhQr4fi12pXbt2+s9//qOMjAwdOnRIv/zyix5++GHdeeedio+Pl/THx6pVq1ZVjRo1JElfffWVbDabHnvsMYdxQ0JCdMcdd9ivXE5ISNDvv/+ufv36OfTLyclR586dlZiYaP9oJ9flH/tLf3wcfOHCBaWmpl5zLv7+/nme36dPH+Xk5GjDhg15xrnefW7t2rVq166dgoOD7c8vW7asHnrooWvW+PXXX8vLy0sDBw68Zt/8fPXVV2rTpo3CwsIcauzSpYskaf369fYaC1ofhTVr1iw1bNhQXl5ecnNzk7u7u9asWZPva6hbt24qW7as/X79+vUlyb6///jjjzpy5Ij69Onj8HFYRESEmjVrds1ainq8Kcy+cPfdd+vUqVN65JFH9MUXXygtLS3PuIVd37mu3K8kqUuXLgoJCdHcuXPtbatWrdKRI0cc9oOS2rZFmVPz5s316quvavr06RoyZIgee+yxAi+qy33t5HrwwQfl5uaW5+PB+vXrq2bNmtdVU2G22d13362vv/5aL7zwgtatW6eMjIxrrpOi7l8hISG6++6788yrsMf3go4fv/zyi/1j4OLY9272Gorqassv6Nhy+f5nhWyQq9DhLigoSD4+Pjpw4ECh+p84cUKSFBoamuexsLAw++O58utX0GPHjh2TJD333HNyd3d3uA0dOlSS8n1R5+rTp49mzJihJ598UqtWrdK2bduUmJioihUrFuqFnp8TJ07Izc1NFStWdGi32WwKCQnJM98KFSrkWYanp2ehxw8JCcm37eLFizp79qxOnDihrKwsvfPOO3nWUdeuXSVdfR3lat++vTIzM7Vp0ybFx8crKChIDRo0UPv27fXNN99IktasWaP27dvbn3Ps2DEZYxQcHJxn7C1bttjHzd2Of/7zn/P0mzx5sowxeX5a58r1lvvdz8Kst8sPVJevM0nX3B+Lss+dOHGiwO1zLcePH1dYWJjKlLm+/3cdO3ZMX375ZZ4a69Spk6fGq62Pa5k2bZqGDBmiJk2a6LPPPtOWLVuUmJiozp0757strrXdctf/9a63oh5vCrMv9O3bV3PmzNGhQ4d0//33q1KlSmrSpIn9PzVS4dd3rvzqc3NzU9++ffX555/bv4Mzb948hYaGqlOnTkUe60a3bVHn9Oijj8rDw0OZmZkaPXp0gcu9cnw3NzdVqFChUO8Fha2pMNvs7bff1vPPP69ly5apTZs2Kl++vHr16qWff/65wNqLun856/h+eS3Fse/d7DUU1dWWX9B8Lt92VsgGuQp9tWzZsmXVrl07ff311/r111+veSVRboFHjx7N0/fIkSMKCgpyaLvabwpd+Vjuc8eOHav77rsv3+dER0fn23769Gl99dVXio2N1QsvvGBvz8zMvKHf6KtQoYKysrJ0/Phxh41ojFFKSor9S6PFJb8vR6ekpMjDw0N+fn5yd3dX2bJl1bdvXz399NP5LiMqKuqa4zRp0kR+fn765ptvdPDgQbVr1042m03t2rXT1KlTlZiYqOTkZIdwFxQUJJvNpo0bN+Z74U1uW+52fOedd3TPPffkO35+b1LXKzegXS53PV75grqRfa5ChQoFbp9rqVixojZt2qScnJzrCnhBQUGqX7++Xn311XwfDwsLs9e4bdu266pRkj7++GO1bt1a7733nkP7mTNnilix7PUUNH5hairq8aaw+8KAAQM0YMAAnTt3Ths2bFBsbKy6d++un376SREREYVe37kKOs4NGDBAr7/+uhYtWqSHHnpIy5cv14gRIxzOdpbUti3KnLKzs/Xoo4+qXLly8vT01BNPPKH//Oc/8vDwyHf8ypUr2+9nZWXpxIkT13ztFbWma20zX19fTZgwQRMmTNCxY8fsZ/F69Oih/fv357v8ou5fN+pqr4PcWopr37uZa/Dy8sr3grmCTkxcbfkFzSd3LlbJBrmK9FMoY8eO1YoVKzRo0CB98cUXeV7Aly5d0sqVK9WjRw/7R6wff/yxQ/GJiYlKSkrSuHHjrrvo6Oho1ahRQ7t379bEiROL9FybzSZjTJ7Q8c9//jPPlXJFOSvUrl07TZkyRR9//LGeffZZe/tnn32mc+fOqV27dkWq81qWLl2q119/3f7R7JkzZ/Tll1+qRYsWKlu2rHx8fNSmTRvt3LlT9evXz/dgm+tq83R3d1fLli0VHx+vw4cP67XXXpMktWjRQm5ubnrppZfsYS9X9+7d9dprr+m3337Tgw8+WOC4zZs312233aYffvhBw4YNu671UBRnzpzR8uXLHT6uWrBggcqUKaOWLVte9blF2efatGmj5cuX69ixY/Zwmp2drcWLF1+zxi5dumjhwoWaN2/eVT+aLeh/ct27d9eKFStUrVo1lStX7qo1fvrpp/muj8Kw2Wx5XkN79uzR5s2bFR4eXqhlXC46OlqhoaFauHChRo4caT9IHzp0SAkJCXneJK5U1ONNUfcFX19fdenSRRcvXlSvXr20b98+RUREFHp9X0utWrXUpEkTzZ07V9nZ2crMzNSAAQMc+pTUti3KnGJjY7Vx40atXr1avr6+atmypUaPHq233norT99PPvlEjRo1st//9NNPlZWVVagf8b6e9VzQNrtccHCw+vfvr927d2v69OkF/gSXM9/P8rNmzZp8jx/VqlWzh8vi2vdcXcPV3n8iIyO1ZMkSZWZm2vudOHFCCQkJCggIKNI4BR1bHn/8cUnWyQa5ihTumjZtqvfee09Dhw5Vo0aNNGTIENWpU0eXLl3Szp07NXv2bNWtW1c9evRQdHS0/vKXv+idd95RmTJl1KVLFx08eFDjx49XeHi4wySvx/vvv68uXbqoU6dO6t+/vypXrqzff/9dSUlJ+u6777RkyZJ8nxcQEKCWLVvq9ddfV1BQkCIjI7V+/Xp98MEHuu222xz65v6K/ezZs+Xv7y8vLy9FRUXle9q0Q4cO6tSpk55//nmlp6erefPm2rNnj2JjY9WgQQP17dv3huZ7pbJly6pDhw4aOXKkcnJyNHnyZKWnp9t/NkCS3nrrLf3pT39SixYtNGTIEEVGRurMmTP65Zdf9OWXX9q/C1CtWjV5e3vrk08+Ua1ateTn56ewsDD7G2q7du00atQoSbKfofP29lazZs20evVq1a9fX5UqVbKP27x5c/3lL3/RgAEDtH37drVs2VK+vr46evSoNm3apHr16mnIkCHy8/PTO++8o379+un333/Xn//8Z1WqVEnHjx/X7t27dfz48Txnhm5EhQoVNGTIECUnJ6tmzZpasWKF/vGPf2jIkCEO33EoSGH3uZdeeknLly9X27Zt9fLLL8vHx0fvvvtunu8P5ueRRx7R3Llz9dRTT+nHH39UmzZtlJOTo61bt6pWrVp6+OGHJUn16tXTunXr9OWXXyo0NFT+/v6Kjo7WK6+8ovj4eDVr1kzDhw9XdHS0Lly4oIMHD2rFihWaNWuWqlSposcff1xvvvmmHn/8cb366quqUaOGVqxYoVWrVhVqXXbv3l1/+9vfFBsbq1atWunHH3/UK6+8oqioqOv6bakyZcrob3/7m5588kn17t1bgwYN0qlTpxQXF1eojxOLerwpzL4waNAgeXt7q3nz5goNDVVKSoomTZqkwMBA+xt8Ydd3YQwcOFCDBw/WkSNH1KxZszyfPpTUti3sOPHx8Zo0aZLGjx9vf4OaNGmSnnvuObVu3Vq9e/d2WO7SpUvl5uamDh06aN++fRo/frzuuOOOq/4HsKg1FWabNWnSRN27d1f9+vVVrlw5JSUl6aOPPlLTpk3zDXZS0fevGxUUFKS2bdtq/Pjx8vX11cyZM7V//36HnyIpzn3PlTVc7f2nb9++ev/99/XYY49p0KBBOnHihKZMmVLkYCdJqamp9mPL6dOnFRsbKy8vL40dO1aSdbKBXZEuv/j/du3aZfr162eqVq1qPDw87Jf0vvzyy/bLgY354wq4yZMnm5o1axp3d3cTFBRkHnvsMftlzrlatWpl6tSpk2ec3CtgXn/99Xzr2L17t3nwwQdNpUqVjLu7uwkJCTFt27Y1s2bNsvfJ72rZX3/91dx///2mXLlyxt/f33Tu3Nns3bvXRERE5LkCcfr06SYqKsqULVvW4YqeK6+IMeaPq1qef/55ExERYdzd3U1oaKgZMmSIOXnypEO/iIgI061btzzzKejKoPzWyeTJk82ECRNMlSpVjIeHh2nQoIFZtWpVvv0HDhxoKleubNzd3U3FihVNs2bNzN///neHfgsXLjQxMTHG3d3dSDKxsbH2x3bv3m0kmRo1ajg859VXXzWSzMiRI/Otdc6cOaZJkybG19fXeHt7m2rVqpnHH3/cbN++3aHf+vXrTbdu3Uz58uWNu7u7qVy5sunWrZvDJfe5V69d/pMDxvzfz2dc7bJ+Y/5vH1u3bp1p3Lix8fT0NKGhoebFF190uPq1OPY5Y4z5z3/+Y+655x7j6elpQkJCzOjRo83s2bOvebWsMX/sRy+//LKpUaOG8fDwMBUqVDBt27Z1+JmbXbt2mebNmxsfHx8jyWEZx48fN8OHDzdRUVHG3d3dlC9f3jRq1MiMGzfOnD171t4v93Xg5+dn/P39zf33328SEhIKdUVlZmamee6550zlypWNl5eXadiwoVm2bFmBV7Hltz6v3M+MMeaf//ynfd41a9Y0c+bMyfe1lp+iHm+utS/Mnz/ftGnTxgQHBxsPDw8TFhZmHnzwQYef+TCmcOv7WvuVMX/8LJG3t/dVr5YuiW1bmHGOHDliKlWqZNq2bevwEy45OTmmR48e5rbbbrPv57mv3R07dpgePXrYa3rkkUfMsWPHHMYt6NhY2LkXZpu98MILpnHjxqZcuXLG09PT3H777ebZZ5+1/+TM5TVf7kbfzwq7H+v//yTPzJkzTbVq1Yy7u7uJiYkxn3zyyXWtk8Lse66u4WrvP/Pnzze1atUyXl5epnbt2mbx4sVFOs7kZoCPPvrIDB8+3FSsWNF4enqaFi1a5HkvKs3Z4Eo2Y4xxTmwEbh6tW7dWWlqa9u7d6+pSgFtKXFycJkyYoOPHjxf7d9OsyGaz6emnn9aMGTNu6RpwY5z2UygAAAAoeYQ7AAAAC+FjWQAAAAvhzB0AAICFEO4AAAAshHAHAABgIUX6EWOry8nJ0ZEjR+Tv71/kP5MCAABcwxijM2fO3NDf5rYSwt1ljhw5cl1/OgkAALje4cOHb+gvc1gF4e4y/v7+kv7YOa7nz5sAAICSl56ervDwcPv7+K2OcHeZ3I9iAwICCHcAAJQyfKXqD3wwDQAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEDdXFwAAsJbk5GSlpaW5uowiCQoKUtWqVV1dBlAsCHcAgGKTnJysmJhaysg47+pSisTb20f79ycR8GAJhDsAQLFJS0tTRsZ5NRkYq4DQSFeXUyjpRw9q65wJSktLI9zBEgh3AIBiFxAaqfJVo11dBnBL4oIKAAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIaUy3E2aNEk2m00jRoywtxljFBcXp7CwMHl7e6t169bat2+f64oEAABwgVIX7hITEzV79mzVr1/foX3KlCmaNm2aZsyYocTERIWEhKhDhw46c+aMiyoFAAAoeaUq3J09e1aPPvqo/vGPf6hcuXL2dmOMpk+frnHjxum+++5T3bp1NX/+fJ0/f14LFixwYcUAAAAlq1SFu6efflrdunVT+/btHdoPHDiglJQUdezY0d7m6empVq1aKSEhoaTLBAAAcBk3VxdQWIsWLdJ3332nxMTEPI+lpKRIkoKDgx3ag4ODdejQoQKXmZmZqczMTPv99PT0YqoWAADANUrFmbvDhw/rmWee0ccffywvL68C+9lsNof7xpg8bZebNGmSAgMD7bfw8PBiqxkAAMAVSkW427Fjh1JTU9WoUSO5ubnJzc1N69ev19tvvy03Nzf7GbvcM3i5UlNT85zNu9zYsWN1+vRp++3w4cNOnQcAAICzlYqPZdu1a6fvv//eoW3AgAGKiYnR888/r9tvv10hISGKj49XgwYNJEkXL17U+vXrNXny5AKX6+npKU9PT6fWDgAAUJJKRbjz9/dX3bp1Hdp8fX1VoUIFe/uIESM0ceJE1ahRQzVq1NDEiRPl4+OjPn36uKJkAAAAlygV4a4wxowZo4yMDA0dOlQnT55UkyZNtHr1avn7+7u6NAAAgBJTasPdunXrHO7bbDbFxcUpLi7OJfUAAADcDErFBRUAAAAoHMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwkFL7FyoA4FaQnJystLQ0V5dRaElJSa4uAbjlEe4A4CaVnJysmJhaysg47+pSiuxS5kVXlwDcsgh3AHCTSktLU0bGeTUZGKuA0EhXl1MoR7/frL3LZysrK8vVpQC3LMIdANzkAkIjVb5qtKvLKJT0owddXQJwy+OCCgAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAsh3AEAAFgI4Q4AAMBCCHcAAAAW4ubqAgAAuBkkJSW5uoQiCQoKUtWqVV1dBm5ChDsAwC0t4/QJSTY99thjri6lSLy9fbR/fxIBD3kQ7gAAt7RL589IMrqzz/OqGBXj6nIKJf3oQW2dM0FpaWmEO+RBuAMAQJJfpaoqXzXa1WUAN4wLKgAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIWUinD33nvvqX79+goICFBAQICaNm2qr7/+2v64MUZxcXEKCwuTt7e3WrdurX379rmwYgAAANcoFeGuSpUqeu2117R9+3Zt375dbdu2Vc+ePe0BbsqUKZo2bZpmzJihxMREhYSEqEOHDjpz5oyLKwcAAChZpSLc9ejRQ127dlXNmjVVs2ZNvfrqq/Lz89OWLVtkjNH06dM1btw43Xfffapbt67mz5+v8+fPa8GCBa4uHQAAoESVinB3uezsbC1atEjnzp1T06ZNdeDAAaWkpKhjx472Pp6enmrVqpUSEhJcWCkAAEDJc3N1AYX1/fffq2nTprpw4YL8/Pz0+eefq3bt2vYAFxwc7NA/ODhYhw4duuoyMzMzlZmZab+fnp5e/IUDuGkkJycrLS3N1WUUWlJSkqtLAFAKlZpwFx0drV27dunUqVP67LPP1K9fP61fv97+uM1mc+hvjMnTdqVJkyZpwoQJTqkXwM0lOTlZMTG1lJFx3tWlFNmlzIuuLgFAKVJqwp2Hh4eqV68uSWrcuLESExP11ltv6fnnn5ckpaSkKDQ01N4/NTU1z9m8K40dO1YjR460309PT1d4eLgTqgfgamlpacrIOK8mA2MVEBrp6nIK5ej3m7V3+WxlZWW5uhQApUipCXdXMsYoMzNTUVFRCgkJUXx8vBo0aCBJunjxotavX6/JkydfdRmenp7y9PQsiXIB3CQCQiNVvmq0q8solPSjB11dAoBSqFSEuxdffFFdunRReHi4zpw5o0WLFmndunVauXKlbDabRowYoYkTJ6pGjRqqUaOGJk6cKB8fH/Xp08fVpQMAAJSoUhHujh07pr59++ro0aMKDAxU/fr1tXLlSnXo0EGSNGbMGGVkZGjo0KE6efKkmjRpotWrV8vf39/FlQMAAJSsUhHuPvjgg6s+brPZFBcXp7i4uJIpCAAA4CZV6n7nDgAAAAUrFWfuANxcStvvxUn8ZhyAWwfhDkCRlObfi5P4zTgA1ke4A1AkpfH34iR+Mw7ArYNwB+C6lKbfi5P4zTgAtw4uqAAAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFODXc3X777Tpx4kSe9lOnTun222935tAAAAC3JKeGu4MHDyo7OztPe2Zmpn777TdnDg0AAHBLcnPGQpcvX27/96pVqxQYGGi/n52drTVr1igyMtIZQwMAANzSnBLuevXqJUmy2Wzq16+fw2Pu7u6KjIzU1KlTnTE0AADALc0p4S4nJ0eSFBUVpcTERAUFBTljGAAAAFzBKeEu14EDB5y5eAAAAFzBqeFOktasWaM1a9YoNTXVfkYv15w5c5w9PAAAwC3FqeFuwoQJeuWVV9S4cWOFhobKZrM5czgAAIBbnlPD3axZszRv3jz17dvXmcMAAADg/3Pq79xdvHhRzZo1c+YQAAAAuIxTz9w9+eSTWrBggcaPH+/MYYBSLTk5WWlpaa4uo9CSkpJcXQIA4CqcGu4uXLig2bNn65tvvlH9+vXl7u7u8Pi0adOcOTxw00tOTlZMTC1lZJx3dSlFdinzoqtLAADkw6nhbs+ePbrzzjslSXv37nV4jIsrACktLU0ZGefVZGCsAkIjXV1OoRz9frP2Lp+trKwsV5cCAMiHU8Pd2rVrnbl4wDICQiNVvmq0q8solPSjB11dAgDgKpx6QQUAAABKllPP3LVp0+aqH79+++23zhweAADgluPUcJf7fbtcly5d0q5du7R3717169fPmUMDAADckpwa7t5888182+Pi4nT27FlnDg0AAHBLcsl37h577DH+riwAAIATuCTcbd68WV5eXq4YGgAAwNKc+rHsfffd53DfGKOjR49q+/bt/NUKAAAAJ3BquAsMDHS4X6ZMGUVHR+uVV15Rx44dnTk0AADALcmp4W7u3LnOXDwAAACu4NRwl2vHjh1KSkqSzWZT7dq11aBBg5IYFgAA4Jbj1HCXmpqqhx9+WOvWrdNtt90mY4xOnz6tNm3aaNGiRapYsaIzhwcAALjlOPVq2b/+9a9KT0/Xvn379Pvvv+vkyZPau3ev0tPTNXz4cGcODQAAcEty6pm7lStX6ptvvlGtWrXsbbVr19a7777LBRVwiuTkZKWlpbm6jEJLSkpydQkAAItxarjLycmRu7t7nnZ3d3fl5OQ4c2jcgpKTkxUTU0sZGeddXUqRXcq86OoSAAAW4dRw17ZtWz3zzDNauHChwsLCJEm//fabnn32WbVr186ZQ+MWlJaWpoyM82oyMFYBoZGuLqdQjn6/WXuXz1ZWVparSwEAWIRTw92MGTPUs2dPRUZGKjw8XDabTcnJyapXr54+/vhjZw6NW1hAaKTKV412dRmFkn70oKtLAABYjFPDXXh4uL777jvFx8dr//79Msaodu3aat++vTOHBQAAuGU55WrZb7/9VrVr11Z6erokqUOHDvrrX/+q4cOH66677lKdOnW0ceNGZwwNAABwS3NKuJs+fboGDRqkgICAPI8FBgZq8ODBmjZtmjOGBgAAuKU5Jdzt3r1bnTt3LvDxjh07aseOHc4YGgAA4JbmlHB37NixfH8CJZebm5uOHz/ujKEBAABuaU4Jd5UrV9b3339f4ON79uxRaGioM4YGAAC4pTkl3HXt2lUvv/yyLly4kOexjIwMxcbGqnv37s4YGgAA4JbmlJ9Ceemll7R06VLVrFlTw4YNU3R0tGw2m5KSkvTuu+8qOztb48aNc8bQAAAAtzSnnLkLDg5WQkKC6tatq7Fjx6p3797q1auXXnzxRdWtW1f/+c9/FBwcXOjlTZo0SXfddZf8/f1VqVIl9erVSz/++KNDH2OM4uLiFBYWJm9vb7Vu3Vr79u0r7qkBAADc1JwS7iQpIiJCK1asUFpamrZu3aotW7YoLS1NK1asUGRkZJGWtX79ej399NPasmWL4uPjlZWVpY4dO+rcuXP2PlOmTNG0adM0Y8YMJSYmKiQkRB06dNCZM2eKeWYAAAA3L6f+hQpJKleunO66664bWsbKlSsd7s+dO1eVKlXSjh071LJlSxljNH36dI0bN0733XefJGn+/PkKDg7WggULNHjw4BsaHwAAoLRw2pk7Zzp9+rQkqXz58pKkAwcOKCUlRR07drT38fT0VKtWrZSQkOCSGgEAAFzB6WfuipsxRiNHjtSf/vQn1a1bV5KUkpIiSXm+xxccHKxDhw4VuKzMzExlZmba7+f+uTQAAIDSqtSduRs2bJj27NmjhQsX5nnMZrM53DfG5Gm73KRJkxQYGGi/hYeHF3u9AAAAJalUhbu//vWvWr58udauXasqVarY20NCQiT93xm8XKmpqVe9Knfs2LE6ffq0/Xb48GHnFA4AAFBCSkW4M8Zo2LBhWrp0qb799ltFRUU5PB4VFaWQkBDFx8fb2y5evKj169erWbNmBS7X09NTAQEBDjcAAIDSrFR85+7pp5/WggUL9MUXX8jf399+hi4wMFDe3t6y2WwaMWKEJk6cqBo1aqhGjRqaOHGifHx81KdPHxdXDwAAUHJKRbh77733JEmtW7d2aJ87d6769+8vSRozZowyMjI0dOhQnTx5Uk2aNNHq1avl7+9fwtUCAAC4TqkId8aYa/ax2WyKi4tTXFyc8wsCAAC4SZWK79wBAACgcAh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQtxcXQAAALg+SUlJri6hSIKCglS1alVXl2F5hDsAAEqZjNMnJNn02GOPubqUIvH29tH+/UkEPCcj3AEAUMpcOn9GktGdfZ5XxagYV5dTKOlHD2rrnAlKS0sj3DkZ4Q4AgFLKr1JVla8a7eoycJPhggoAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhhDsAAAALIdwBAABYCOEOAADAQgh3AAAAFkK4AwAAsBDCHQAAgIUQ7gAAACyEcAcAAGAhpSbcbdiwQT169FBYWJhsNpuWLVvm8LgxRnFxcQoLC5O3t7dat26tffv2uaZYAAAAFyk14e7cuXO64447NGPGjHwfnzJliqZNm6YZM2YoMTFRISEh6tChg86cOVPClQIAALiOm6sLKKwuXbqoS5cu+T5mjNH06dM1btw43XfffZKk+fPnKzg4WAsWLNDgwYNLslQAAACXKTVn7q7mwIEDSklJUceOHe1tnp6eatWqlRISElxYGQAAQMkqNWfuriYlJUWSFBwc7NAeHBysQ4cOFfi8zMxMZWZm2u+np6c7p0AAAIASYokzd7lsNpvDfWNMnrbLTZo0SYGBgfZbeHi4s0sEAABwKkuEu5CQEEn/dwYvV2pqap6zeZcbO3asTp8+bb8dPnzYqXUCAAA4myXCXVRUlEJCQhQfH29vu3jxotavX69mzZoV+DxPT08FBAQ43AAAAEqzUvOdu7Nnz+qXX36x3z9w4IB27dql8uXLq2rVqhoxYoQmTpyoGjVqqEaNGpo4caJ8fHzUp08fF1YNAABQskpNuNu+fbvatGljvz9y5EhJUr9+/TRv3jyNGTNGGRkZGjp0qE6ePKkmTZpo9erV8vf3d1XJAAAAJa7UhLvWrVvLGFPg4zabTXFxcYqLiyu5ogAAAG4ylvjOHQAAAP5AuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAACyHcAQAAWAjhDgAAwEIIdwAAABbi5uoCbhXJyclKS0tzdRlFkpmZKU9PT1eXUWhJSUmuLgEAAJcj3JWA5ORkxcTUUkbGeVeXUjQ2m2SMq6soskuZF11dAgAALkO4KwFpaWnKyDivJgNjFRAa6epyCuXo95u1d/ls3dnneVWMinF1OYWSW3NWVparSwEAwGUIdyUoIDRS5atGu7qMQkk/elCS5FepaqmrGQCAWxkXVAAAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwAAAAuxXLibOXOmoqKi5OXlpUaNGmnjxo2uLgkAAKDEWCrcLV68WCNGjNC4ceO0c+dOtWjRQl26dFFycrKrSwMAACgRlgp306ZN0xNPPKEnn3xStWrV0vTp0xUeHq733nvP1aUBAACUCMuEu4sXL2rHjh3q2LGjQ3vHjh2VkJDgoqoAAABKlpurCyguaWlpys7OVnBwsEN7cHCwUlJS8n1OZmamMjMz7fdPnz4tSUpPTy/W2s6ePStJ+v3Qj8rKzCjWZTtL+tFDkqTTv/0sdzebi6spHGouGaWxZql01k3NJYOaS0Z6yh9fkTp79myxv8/mLs8YU6zLLbWMRfz2229GkklISHBo//vf/26io6PzfU5sbKyRxI0bN27cuHGzwO3w4cMlETluepY5cxcUFKSyZcvmOUuXmpqa52xerrFjx2rkyJH2+zk5Ofr9999VoUIF2Ww37/+E0tPTFR4ersOHDysgIMDV5ZS4W3n+zJ25M/dbB3Mv/NyNMTpz5ozCwsJKoLqbn2XCnYeHhxo1aqT4+Hj17t3b3h4fH6+ePXvm+xxPT095eno6tN12223OLLNYBQQE3HIv+MvdyvNn7sz9VsPcmfu1BAYGOrma0sMy4U6SRo4cqb59+6px48Zq2rSpZs+ereTkZD311FOuLg0AAKBEWCrcPfTQQzpx4oReeeUVHT16VHXr1tWKFSsUERHh6tIAAABKhKXCnSQNHTpUQ4cOdXUZTuXp6anY2Ng8HynfKm7l+TN35n6rYe7MHUVnM4brhgEAAKzCMj9iDAAAAMIdAACApRDuAAAALIRwBwAAYCGEu1Li5MmT6tu3rwIDAxUYGKi+ffvq1KlTV31O//79ZbPZHG733HNPyRRcjK5n7pcbPHiwbDabpk+f7rQaneV65h4XF6eYmBj5+vqqXLlyat++vbZu3VoyBRejos790qVLev7551WvXj35+voqLCxMjz/+uI4cOVJyRReT69nuS5cuVadOnRQUFCSbzaZdu3aVSK3FYebMmYqKipKXl5caNWqkjRs3XrX/+vXr1ahRI3l5een222/XrFmzSqjS4leUuR89elR9+vRRdHS0ypQpoxEjRpRcoU5QlLkvXbpUHTp0UMWKFRUQEKCmTZtq1apVJVht6UK4KyX69OmjXbt2aeXKlVq5cqV27dqlvn37XvN5nTt31tGjR+23FStWlEC1xet65y5Jy5Yt09atW0vtn6S5nrnXrFlTM2bM0Pfff69NmzYpMjJSHTt21PHjx0uo6uJR1LmfP39e3333ncaPH6/vvvtOS5cu1U8//aR77723BKsuHtez3c+dO6fmzZvrtddeK6Eqi8fixYs1YsQIjRs3Tjt37lSLFi3UpUsXJScn59v/wIED6tq1q1q0aKGdO3fqxRdf1PDhw/XZZ5+VcOU3rqhzz8zMVMWKFTVu3DjdcccdJVxt8Srq3Dds2KAOHTpoxYoV2rFjh9q0aaMePXpo586dJVx5KeHiv22LQvjhhx+MJLNlyxZ72+bNm40ks3///gKf169fP9OzZ88SqNB5rnfuxhjz66+/msqVK5u9e/eaiIgI8+abbzq52uJ1I3O/3OnTp40k88033zijTKcorrlv27bNSDKHDh1yRplOcaNzP3DggJFkdu7c6cQqi8/dd99tnnrqKYe2mJgY88ILL+Tbf8yYMSYmJsahbfDgweaee+5xWo3OUtS5X65Vq1bmmWeecVJlzncjc89Vu3ZtM2HChOIuzRI4c1cKbN68WYGBgWrSpIm97Z577lFgYKASEhKu+tx169apUqVKqlmzpgYNGqTU1FRnl1usrnfuOTk56tu3r0aPHq06deqURKnF7ka2e66LFy9q9uzZCgwMLFX/0y+OuUvS6dOnZbPZStXfjC6uuZcGFy9e1I4dO9SxY0eH9o4dOxY4182bN+fp36lTJ23fvl2XLl1yWq3F7XrmbhXFMfecnBydOXNG5cuXd0aJpR7hrhRISUlRpUqV8rRXqlRJKSkpBT6vS5cu+uSTT/Ttt99q6tSpSkxMVNu2bZWZmenMcovV9c598uTJcnNz0/Dhw51ZnlNd79wl6auvvpKfn5+8vLz05ptvKj4+XkFBQc4qtdjdyNxzXbhwQS+88IL69OlTqv7oenHMvbRIS0tTdna2goODHdqDg4MLnGtKSkq+/bOyspSWlua0Wovb9czdKopj7lOnTtW5c+f04IMPOqPEUo9w50JxcXF5Lni48rZ9+3ZJks1my/N8Y0y+7bkeeughdevWTXXr1lWPHj309ddf66efftK///1vp82psJw59x07duitt97SvHnzrrp+XMXZ212S2rRpo127dikhIUGdO3fWgw8+eFOctS2JuUt/XFzx8MMPKycnRzNnziz2eVyPkpp7aXTlvK411/z659deGhR17lZyvXNfuHCh4uLitHjx4nz/IwQL/m3Z0mTYsGF6+OGHr9onMjJSe/bs0bFjx/I8dvz48Tz/87ma0NBQRURE6Oeffy5yrcXNmXPfuHGjUlNTVbVqVXtbdna2Ro0apenTp+vgwYM3VPuNKont7uvrq+rVq6t69eq65557VKNGDX3wwQcaO3bsDdV+o0pi7pcuXdKDDz6oAwcO6Ntvv71pztqV9Ou9NAgKClLZsmXznK1JTU0tcK4hISH59ndzc1OFChWcVmtxu565W8WNzH3x4sV64okntGTJErVv396ZZZZqhDsXCgoKKtRHZU2bNtXp06e1bds23X333ZKkrVu36vTp02rWrFmhxztx4oQOHz6s0NDQ6665uDhz7n379s3zou/UqZP69u2rAQMG3HjxN6ikt7v0x/+Ib4aP450999xg9/PPP2vt2rU31Zu9K7b7zc7Dw0ONGjVSfHy8evfubW+Pj49Xz549831O06ZN9eWXXzq0rV69Wo0bN5a7u7tT6y1O1zN3q7jeuS9cuFADBw7UwoUL1a1bt5IotfRy1ZUcKJrOnTub+vXrm82bN5vNmzebevXqme7duzv0iY6ONkuXLjXGGHPmzBkzatQok5CQYA4cOGDWrl1rmjZtaipXrmzS09NdMYXrVtS556c0Xi1rTNHnfvbsWTN27FizefNmc/DgQbNjxw7zxBNPGE9PT7N3715XTOG6FXXuly5dMvfee6+pUqWK2bVrlzl69Kj9lpmZ6YopXLfr2edPnDhhdu7caf79738bSWbRokVm586d5ujRoyVdfpEsWrTIuLu7mw8++MD88MMPZsSIEcbX19ccPHjQGGPMCy+8YPr27Wvv/7///c/4+PiYZ5991vzwww/mgw8+MO7u7uZf//qXq6Zw3Yo6d2OM2blzp9m5c6dp1KiR6dOnj9m5c6fZt2+fK8q/IUWd+4IFC4ybm5t59913HV7bp06dctUUbmqEu1LixIkT5tFHHzX+/v7G39/fPProo+bkyZMOfSSZuXPnGmOMOX/+vOnYsaOpWLGicXd3N1WrVjX9+vUzycnJJV/8DSrq3PNTWsNdUeeekZFhevfubcLCwoyHh4cJDQ019957r9m2bVvJF3+Dijr33J8Aye+2du3aEq//RlzPPj937tx85x4bG1uitV+Pd99910RERBgPDw/TsGFDs379evtj/fr1M61atXLov27dOtOgQQPj4eFhIiMjzXvvvVfCFRefos49v20cERFRskUXk6LMvVWrVvnOvV+/fiVfeClgM+b/fxMVAAAApR5XywIAAFgI4Q4AAMBCCHcAAAAWQrgDAACwEMIdAACAhRDuAAAALIRwBwAAYCGEOwBwobi4ON155503zXIAlH6EO6AUS01N1eDBg1W1alV5enoqJCREnTp10ubNm11d2k3rs88+U+vWrRUYGCg/Pz/Vr19fr7zyin7//XdXl1ZoNptNy5Ytc2h77rnntGbNGtcUBOCmQrgDSrH7779fu3fv1vz58/XTTz9p+fLlat269U0RVC5duuTqEvIYN26cHnroId111136+uuvtXfvXk2dOlW7d+/WRx99dN3LzW+uJT1/Pz8/VahQoUTHBHCTcvXfPwNwfU6ePGkkmXXr1l21nyQzc+ZM07lzZ+Pl5WUiIyPNp59+6tDn119/NQ8++KC57bbbTPny5c29995rDhw4YH9827Ztpn379qZChQomICDAtGzZ0uzYsSPPOO+995659957jY+Pj3n55ZdNbGysueOOO8wHH3xgwsPDja+vr3nqqadMVlaWmTx5sgkODjYVK1Y0f//73x2WNXXqVFO3bl3j4+NjqlSpYoYMGWLOnDljf3zu3LkmMDDQrFy50sTExBhfX1/TqVMnc+TIkQLXw9atW40kM3369ALXZ66ZM2ea22+/3bi7u5uaNWuaDz/8sEhzjYqKMjabzeTk5JhTp06ZQYMGmYoVKxp/f3/Tpk0bs2vXLvuycp9X2HUdERGR798VvXI52dnZZsKECaZy5crGw8PD3HHHHebrr7+2P577t3g/++wz07p1a+Pt7W3q169vEhISClyHAEoHztwBpZSfn5/8/Py0bNkyZWZmXrXv+PHj7Wf5HnvsMT3yyCNKSkqSJJ0/f15t2rSRn5+fNmzYoE2bNsnPz0+dO3fWxYsXJUlnzpxRv379tHHjRm3ZskU1atRQ165ddebMGYdxYmNj1bNnT33//fcaOHCgJOm///2vvv76a61cuVILFy7UnDlz1K1bN/36669av369Jk+erJdeeklbtmyxL6dMmTJ6++23tXfvXs2fP1/ffvutxowZ4zDW+fPn9cYbb+ijjz7Shg0blJycrOeee67AdfDJJ5/Iz89PQ4cOzffx2267TZL0+eef65lnntGoUaO0d+9eDR48WAMGDNDatWuvOddffvlFn376qT777DPt2rVLktStWzelpKRoxYoV2rFjhxo2bKh27doVeHb1Wus6MTFRkjR37lwdPXrUfv9Kb731lqZOnao33nhDe/bsUadOnXTvvffq559/dug3btw4Pffcc9q1a5dq1qypRx55RFlZWQWuRwClgKvTJYDr969//cuUK1fOeHl5mWbNmpmxY8ea3bt3O/SRZJ566imHtiZNmpghQ4YYY4z54IMPTHR0tMnJybE/npmZaby9vc2qVavyHTcrK8v4+/ubL7/80mGcESNGOPSLjY01Pj4+Jj093d7WqVMnExkZabKzs+1t0dHRZtKkSQXO89NPPzUVKlSw3587d66RZH755Rd727vvvmuCg4MLXEaXLl1M/fr1C3w8V7NmzcygQYMc2h544AHTtWtX+/2C5uru7m5SU1PtbWvWrDEBAQHmwoULDn2rVatm3n//ffvzLj/jdqWC1vXnn3+eZ/zLlxMWFmZeffVVhz533XWXGTp0qDHm/87c/fOf/7Q/vm/fPiPJJCUlFVgPgJsfZ+6AUuz+++/XkSNHtHz5cnXq1Enr1q1Tw4YNNW/ePId+TZs2zXM/98zdjh079Msvv8jf399+NrB8+fK6cOGC/vvf/0r648KNp556SjVr1lRgYKACAwN19uxZJScnOyy3cePGeWqMjIyUv7+//X5wcLBq166tMmXKOLSlpqba769du1YdOnRQ5cqV5e/vr8cff1wnTpzQuXPn7H18fHxUrVo1+/3Q0FCHZVzJGCObzVbg47mSkpLUvHlzh7bmzZvb19fV5hoREaGKFSva7+/YsUNnz55VhQoV7OvWz89PBw4csK/bKxV2XV9Nenq6jhw5Uqh51K9f3/7v0NBQew0ASi83VxcA4MZ4eXmpQ4cO6tChg15++WU9+eSTio2NVf/+/a/6vNygk5OTo0aNGumTTz7J0yc3qPTv31/Hjx/X9OnTFRERIU9PTzVt2tT+sW0uX1/fPMtwd3fPM25+bTk5OZKkQ4cOqWvXrnrqqaf0t7/9TeXLl9emTZv0xBNPOFykkN8yjDEFzrdmzZratGmTLl26lOe5V7oyBOYXDPOb65VtOTk5Cg0N1bp16/L0zf0Y+EqFXdeFUZh5XL4uLt8nAJRenLkDLKZ27doOZ7gkOXyfLfd+TEyMJKlhw4b6+eefValSJVWvXt3hFhgYKEnauHGjhg8frq5du6pOnTry9PRUWlqaU+rfvn27srKyNHXqVN1zzz2qWbOmjhw5csPL7dOnj86ePauZM2fm+/ipU6ckSbVq1dKmTZscHktISFCtWrWKPGbDhg2VkpIiNze3POs2KCgo3+cUZl27u7srOzu7wHEDAgIUFhZWbPMAULoQ7oBS6sSJE2rbtq0+/vhj7dmzRwcOHNCSJUs0ZcoU9ezZ06HvkiVLNGfOHP3000+KjY3Vtm3bNGzYMEnSo48+qqCgIPXs2VMbN27UgQMHtH79ej3zzDP69ddfJUnVq1fXRx99pKSkJG3dulWPPvqovL29nTKvatWqKSsrS++8847+97//6aOPPtKsWbNueLlNmjTRmDFjNGrUKI0ZM0abN2/WoUOHtGbNGj3wwAOaP3++JGn06NGaN2+eZs2apZ9//lnTpk3T0qVLr3qxRkHat2+vpk2bqlevXlq1apUOHjyohIQEvfTSS9q+fXu+zynMuo6MjNSaNWuUkpKikydP5ruc0aNHa/LkyVq8eLF+/PFHvfDCC9q1a5eeeeaZIs8DQOlCuANKKT8/PzVp0kRvvvmmWrZsqbp162r8+PEaNGiQZsyY4dB3woQJWrRokerXr6/58+frk08+Ue3atSX98d21DRs2qGrVqrrvvvtUq1YtDRw4UBkZGQoICJAkzZkzRydPnlSDBg3Ut29fDR8+XJUqVXLKvO68805NmzZNkydPVt26dfXJJ59o0qRJxbLsyZMna8GCBdq6das6deqkOnXqaOTIkapfv7769esnSerVq5feeustvf7666pTp47ef/99zZ07V61bty7yeDabTStWrFDLli01cOBA1axZUw8//LAOHjyo4ODgfJ9TmHU9depUxcfHKzw8XA0aNMh3OcOHD9eoUaM0atQo1atXTytXrtTy5ctVo0aNIs8DQOliM1f7kgqAUs9ms+nzzz9Xr169XF0KAKAEcOYOAADAQgh3AAAAFsJPoQAWxzcvAODWwpk7AAAACyHcAQAAWAjhDgAAwEIIdwAAABZCuAMAALAQwh0AAICFEO4AAAAshHAHAABgIYQ7AAAAC/l/YwcTQOJoy+sAAAAASUVORK5CYII=",
      "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": 91,
   "id": "e7033192-d47d-4568-aa9a-63a7aac918c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savetxt(\n",
    "    \"figures/reconstruction_spearmans/GSFA_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 pandas as pd\n",
    "import seaborn as sns\n",
    "from scipy.cluster.hierarchy import dendrogram, fcluster, linkage\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "\n",
    "# Put into df\n",
    "pert_emb_df = pd.DataFrame(beta[:-1], index=pert_labels)\n",
    "pert_emb_df = pert_emb_df.drop(columns=[0]).astype(\"float\")\n",
    "\n",
    "# Compute the pairwise distances\n",
    "df = pert_emb_df.drop(index=[\"ctrl\"])\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 = 10  # 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": 10,
   "id": "d187a4f5-b9d8-4dfd-a98b-c759446de3e6",
   "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": "a7dedf5b-de4c-4142-bbd4-1a2e4a4a4ec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create descriptions for each cluster (for demonstration purpose)\n",
    "cluster_descriptions = {\n",
    "    1: \"\",\n",
    "    2: \"Embryonic development\",\n",
    "    3: \"Transcription cis-regulatory region binding, metabolic process\",\n",
    "    4: \"protein serine/threonine kinase activity\",\n",
    "    5: \"\",\n",
    "    6: \"DNA-binding transcription factor activity, RNA polymerase II transcription regulatory region\",\n",
    "    7: \"\",\n",
    "    8: \"DNA binding, negative regulation of biosynthetic process\",\n",
    "    9: \"\",\n",
    "    10: \"chromatin\",\n",
    "    11: \"epidermal growth factor\",\n",
    "    12: \"collagen type II trimer\",\n",
    "    13: \"\",\n",
    "    14: \"\",\n",
    "    15: \"C/EBP complex\",\n",
    "    16: \"\",\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/GSFA_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/GSFA_norman_clustermap_descriptions.png\",\n",
    "    dpi=600,\n",
    "    bbox_inches=\"tight\",\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "id": "23b3a2e6-f20f-4626-bcc9-32e9e4ef0d54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save grpofiler results\n",
    "with open(\n",
    "    \"figures/pert_embedding_cluster_gprofiler/GSFA_norman_pert_emb_gprofiler.pickle\",\n",
    "    \"wb\",\n",
    ") as handle:\n",
    "    pickle.dump(cluster_process, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc06a81d-1fbe-4b46-82f2-0464e55727e7",
   "metadata": {},
   "source": [
    "## PR Curve+AUC\n",
    "- AUCPR using prior graph as binary label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "id": "73b5c140-05f1-486b-9d28-7592ed1556bb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "               AUC\n",
      "StringDB  0.732776\n",
      "CORUM     0.946837\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/GSFA_norman_aucpr.csv\")\n",
    "f1_df.to_csv(\"figures/pert_embedding_recall/GSFA_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, W.T, fdr=5e-2)"
   ]
  },
  {
   "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": 105,
   "id": "ff945301-46b0-4beb-9f4e-9bc17ede92d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\n",
    "    \"figures/factor_enrichments/GSFA_norman_factor_enrichment.pickle\", \"wb\"\n",
    ") as handle:\n",
    "    pickle.dump(filtered_factor_to_go, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e5f583b-5c2e-45ea-a210-33b0f6ab2b3e",
   "metadata": {},
   "source": [
    "## Overlap with prior knowledge/ground truth (stringdb)\n",
    "- Group A: a set of GO terms associated with a perturbation (either drivers from msigdb, or from literature) and its neighbors in stringdb\n",
    "- Group B: a set of GO terms from the interpretability analysis\n",
    "- Hypergeometric test on the overlap of the two groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3f70b845-17f9-429e-9bdf-6d36d6fb0a44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load precomputed factor_to_go dict if available\n",
    "filtered_factor_to_go = read_aws_pickle(\n",
    "    \"s3://pert-spectra/figures/factor_enrichments/GSFA_norman_factor_enrichment.pickle\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0a0ca673-3177-412f-937d-f0ec50c3824f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get neighbors for each perturbation\n",
    "pert_neighbors = retrieve_stringdb_neighbors(pert_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "28d8b858-9326-489e-be3f-00a756eb3fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct group A - known processes for each perturbation\n",
    "gene_sets = read_aws_pickle(\"s3://pert-spectra/references/GO_to_Gene.pickle\")\n",
    "# BP only go terms\n",
    "go_reference = read_aws_csv(\"s3://pert-spectra/references/GO_terms.txt.gz\", zipped=True)\n",
    "go_bp = go_reference[go_reference[\"go_category\"] == \"biological_process\"]\n",
    "go_bp_ids = set(go_bp[\"go_id\"].values)\n",
    "filtered_go_terms = {key: gene_sets[key] for key in go_bp_ids if key in gene_sets}\n",
    "\n",
    "# GO terms per perturbation AND its neighbors in stringdb\n",
    "pert_to_go = {key: set() for key in df.index}\n",
    "for goterm in filtered_go_terms:\n",
    "    for pert in df.index:\n",
    "        if pert in filtered_go_terms[goterm] and set(\n",
    "            filtered_go_terms[goterm]\n",
    "        ).intersection(pert_neighbors[pert]):\n",
    "            pert_to_go[pert].add(goterm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b3b81d9d-7264-4dac-b14b-1f35c84479c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# construct set B - model identified processes for each perturbation\n",
    "n = 10  # number of top factors to get processes from\n",
    "model_pert_to_go = {}\n",
    "for pert in pert_emb_df.index:\n",
    "    if pert in [\"ctrl\", \"intergenic\", \"basal\"]:\n",
    "        continue\n",
    "    # get top factors\n",
    "    pert_emb_df = pd.DataFrame(beta[:-1], index=pert_labels)\n",
    "    pert_loading = pert_emb_df.loc[pert].to_numpy()[1:].astype(float)\n",
    "    ctrl_loading = pert_emb_df.loc[\"ctrl\"].to_numpy()[1:].astype(float)\n",
    "    delta_loading = np.abs(np.log(np.abs(pert_loading)) - np.log(np.abs(ctrl_loading)))\n",
    "    top_n_factors = np.argpartition(np.array(delta_loading), -n)[-n:]\n",
    "    # get processes\n",
    "    model_processes = set()\n",
    "    for f in top_n_factors:\n",
    "        proc = filtered_factor_to_go[f]\n",
    "        model_processes = model_processes.union(set(proc[\"Term\"]))\n",
    "    model_pert_to_go[pert] = model_processes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f3a8ecfa-1184-4a30-b787-1ce5e56a5562",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overlap for AHR: 5 out of 12 in researchDB\n",
      "P-value for AHR: 2.1319858454369287e-07\n",
      "Overlap for ARID1A: 5 out of 16 in researchDB\n",
      "P-value for ARID1A: 2.2083684093482343e-06\n",
      "Overlap for BAK1: 5 out of 43 in researchDB\n",
      "P-value for BAK1: 0.00046615592180678417\n",
      "Overlap for BCL2L11: 3 out of 29 in researchDB\n",
      "P-value for BCL2L11: 0.008310750164771933\n",
      "Overlap for CBL: 5 out of 28 in researchDB\n",
      "P-value for CBL: 0.00014285653454198518\n",
      "Overlap for CDKN1A: 3 out of 37 in researchDB\n",
      "P-value for CDKN1A: 0.0090102759706699\n",
      "Overlap for CDKN1B: 2 out of 38 in researchDB\n",
      "P-value for CDKN1B: 0.08621078058192475\n",
      "Overlap for CDKN1C: 3 out of 15 in researchDB\n",
      "P-value for CDKN1C: 0.0008590529302742615\n",
      "Overlap for CEBPA: 10 out of 34 in researchDB\n",
      "P-value for CEBPA: 9.605689713987963e-11\n",
      "Overlap for CEBPB: 8 out of 33 in researchDB\n",
      "P-value for CEBPB: 3.5589724354327065e-08\n",
      "Overlap for CEBPE: 5 out of 9 in researchDB\n",
      "P-value for CEBPE: 9.287161274422125e-08\n",
      "Overlap for CITED1: 6 out of 14 in researchDB\n",
      "P-value for CITED1: 1.0848735073457748e-08\n",
      "Overlap for COL1A1: 3 out of 29 in researchDB\n",
      "P-value for COL1A1: 0.005941572402352278\n",
      "Overlap for COL2A1: 0 out of 12 in researchDB\n",
      "P-value for COL2A1: 1.0\n",
      "Overlap for DUSP9: 0 out of 5 in researchDB\n",
      "P-value for DUSP9: 1.0\n",
      "Overlap for EGR1: 9 out of 18 in researchDB\n",
      "P-value for EGR1: 6.076427546891436e-13\n",
      "Overlap for ETS2: 4 out of 5 in researchDB\n",
      "P-value for ETS2: 1.2643534179980056e-07\n",
      "Overlap for FOSB: 3 out of 11 in researchDB\n",
      "P-value for FOSB: 0.0005638568648125962\n",
      "Overlap for FOXA1: 7 out of 15 in researchDB\n",
      "P-value for FOXA1: 1.0193815225213786e-09\n",
      "Overlap for FOXL2: 4 out of 7 in researchDB\n",
      "P-value for FOXL2: 1.5792583093509072e-06\n",
      "Overlap for FOXO4: 4 out of 8 in researchDB\n",
      "P-value for FOXO4: 3.8416764921953e-06\n",
      "Overlap for HK2: 1 out of 13 in researchDB\n",
      "P-value for HK2: 0.15293459959857644\n",
      "Overlap for HNF4A: 6 out of 18 in researchDB\n",
      "P-value for HNF4A: 2.0064461479781388e-07\n",
      "Overlap for IKZF3: 2 out of 6 in researchDB\n",
      "P-value for IKZF3: 0.0033188033960969143\n",
      "Overlap for IRF1: 4 out of 19 in researchDB\n",
      "P-value for IRF1: 0.00011774780147504422\n",
      "Overlap for JUN: 9 out of 38 in researchDB\n",
      "P-value for JUN: 6.800257732829314e-10\n",
      "Overlap for KIF18B: 4 out of 5 in researchDB\n",
      "P-value for KIF18B: 1.02966073190364e-07\n",
      "Overlap for KIF2C: 1 out of 8 in researchDB\n",
      "P-value for KIF2C: 0.11717918583009142\n",
      "Overlap for KMT2A: 4 out of 14 in researchDB\n",
      "P-value for KMT2A: 5.737318595565771e-05\n",
      "Overlap for LHX1: 3 out of 15 in researchDB\n",
      "P-value for LHX1: 0.001550750199766002\n",
      "Overlap for MAP2K3: 5 out of 18 in researchDB\n",
      "P-value for MAP2K3: 6.727417522597989e-06\n",
      "Overlap for MAP2K6: 2 out of 21 in researchDB\n",
      "P-value for MAP2K6: 0.027671721378544775\n",
      "Overlap for MAPK1: 9 out of 55 in researchDB\n",
      "P-value for MAPK1: 1.2457826721768169e-08\n",
      "Overlap for MEIS1: 2 out of 11 in researchDB\n",
      "P-value for MEIS1: 0.011117126714865002\n",
      "Overlap for NCL: 3 out of 8 in researchDB\n",
      "P-value for NCL: 0.0002326125897755479\n",
      "Overlap for PLK4: 1 out of 5 in researchDB\n",
      "P-value for PLK4: 0.059561646093364874\n",
      "Overlap for POU3F2: 2 out of 9 in researchDB\n",
      "P-value for POU3F2: 0.008122844993866516\n",
      "Overlap for PTPN1: 1 out of 14 in researchDB\n",
      "P-value for PTPN1: 0.1896342601320953\n",
      "Overlap for PTPN12: 0 out of 5 in researchDB\n",
      "P-value for PTPN12: 1.0\n",
      "Overlap for S1PR2: 0 out of 5 in researchDB\n",
      "P-value for S1PR2: 1.0\n",
      "Overlap for SGK1: 7 out of 17 in researchDB\n",
      "P-value for SGK1: 4.12440746718372e-09\n",
      "Overlap for SLC38A2: 1 out of 11 in researchDB\n",
      "P-value for SLC38A2: 0.12020964276782159\n",
      "Overlap for SNAI1: 3 out of 13 in researchDB\n",
      "P-value for SNAI1: 0.0005804369395604471\n",
      "Overlap for SPI1: 12 out of 25 in researchDB\n",
      "P-value for SPI1: 8.848408979799951e-16\n",
      "Overlap for STIL: 1 out of 13 in researchDB\n",
      "P-value for STIL: 0.15204982937479228\n",
      "Overlap for TBX2: 4 out of 8 in researchDB\n",
      "P-value for TBX2: 1.9410436163852128e-06\n",
      "Overlap for TBX3: 4 out of 14 in researchDB\n",
      "P-value for TBX3: 2.0954324972493857e-05\n",
      "Overlap for TGFBR2: 3 out of 55 in researchDB\n",
      "P-value for TGFBR2: 0.020165641355063494\n",
      "Overlap for TP73: 4 out of 17 in researchDB\n",
      "P-value for TP73: 5.875369519782811e-05\n"
     ]
    }
   ],
   "source": [
    "pvals = perturbation_signal_recovery(\n",
    "    pert_to_go,\n",
    "    model_pert_to_go,\n",
    "    list(filtered_go_terms.keys()),\n",
    "    list(pert_emb_df.index),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "712834c7-e823-4117-a2f6-dc704c391d7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save as csv for visualization\n",
    "pd.DataFrame.from_dict(data=pvals, orient=\"index\").to_csv(\n",
    "    \"figures/process_recovery_hypergeo_pvals/GSFA_norman_hypergeo_neighbors_recovery_pvalues.csv\",\n",
    "    header=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde32c38-f16b-4f39-9dcd-5e5f480b3720",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}