{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predict causal genes for cell state transitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will predict causal genes for cell state transitions using the global cell state manifold and a database of genetic perturbations.\n",
    "\n",
    "Let's begin by importing the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import scanpy as sc\n",
    "\n",
    "from scmg.model.causal_prediction import CausalGenePredictor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the required datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_ref = sc.read_h5ad('data/ref_global_cell_state_manifold.h5ad')\n",
    "\n",
    "adata_pert = sc.read_h5ad('data/pseudo_bulk_perturbation_database.h5ad')\n",
    "\n",
    "# Mask out the direct target genes\n",
    "for i in range(adata_pert.shape[0]):\n",
    "    pg = adata_pert.obs['perturbed_gene'].iloc[i]\n",
    "    \n",
    "    if pg in adata_pert.var_names:\n",
    "        adata_pert.X[i, adata_pert.var_names.get_loc(pg)] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, let's consider the transition from epiblast to nascent mesoderm. The cell state transition can be described as a vector of gene expression shifts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG2CAYAAACZEEfAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwmUlEQVR4nO3de1RVdf7/8dfhqtzDBCQFTU2lvGuKWWki6FhLv+qk5Sg6jLUc8ALlqDNeshtGmaZ5KafR5jtafa20X5g6RF5KSQ3FzNS0NCy5WAoIJijs3x99Pd/OoMWBAwe2z8daZy3OZ3/23u/Nacmrz/7sz7EYhmEIAADApFycXQAAAEBtIuwAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTc2rYeeKJJ2SxWGxe7du3t26/dOmS4uPj1aRJE/n4+GjEiBHKy8uzOUZ2draGDBkiLy8vBQUFafr06bpy5UpdXwoAAKin3JxdwO23364PP/zQ+t7N7f9KSkxM1KZNm7R+/Xr5+/srISFBw4cP165duyRJ5eXlGjJkiEJCQrR7927l5ORo3Lhxcnd317PPPlvn1wIAAOofizO/CPSJJ57Qxo0blZWVVWlbYWGhmjZtqnXr1mnkyJGSpKNHj6pDhw7KyMhQ7969tXnzZt1///06c+aMgoODJUkrV67UjBkzdPbsWXl4eNTl5QAAgHrI6SM7x48fV2hoqBo1aqTIyEglJycrLCxMmZmZunz5sqKioqx927dvr7CwMGvYycjIUMeOHa1BR5JiYmI0adIkHT58WF27dr3mOUtLS1VaWmp9X1FRoXPnzqlJkyayWCy1d7EAAMBhDMPQhQsXFBoaKheX68/McWrY6dWrl9asWaN27dopJydH8+fP1913360vvvhCubm58vDwUEBAgM0+wcHBys3NlSTl5ubaBJ2r269uu57k5GTNnz/fsRcDAACc4vTp02revPl1tzs17AwePNj6c6dOndSrVy+Fh4frf/7nf9S4ceNaO++sWbOUlJRkfV9YWKiwsDCdPn1afn5+tXZe4EZ3x7ytzi6hTn0xP8bZJQCmVlRUpBYtWsjX1/dX+zn9NtYvBQQE6LbbbtOJEyc0cOBAlZWVqaCgwGZ0Jy8vTyEhIZKkkJAQ7d271+YYV5/WutrnWjw9PeXp6Vmp3c/Pj7AD1CIXTy9nl1Cn+PcEqBu/NQWlXq2zU1xcrK+//lrNmjVT9+7d5e7urvT0dOv2Y8eOKTs7W5GRkZKkyMhIHTp0SPn5+dY+aWlp8vPzU0RERJ3XDwAA6h+njuw8/vjjeuCBBxQeHq4zZ85o3rx5cnV11UMPPSR/f3/FxcUpKSlJgYGB8vPz0+TJkxUZGanevXtLkqKjoxUREaGxY8cqJSVFubm5mj17tuLj4685cgMAAG48Tg073333nR566CH9+OOPatq0qfr27atPP/1UTZs2lSQtWrRILi4uGjFihEpLSxUTE6Ply5db93d1dVVqaqomTZqkyMhIeXt7KzY2Vk8++aSzLgkAANQzTl1np74oKiqSv7+/CgsLuccO1KKWMzc5u4Q6dWrBEGeXAJhaVf9+16s5OwAAAI5G2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZWb8LOggULZLFYNG3aNGvbpUuXFB8fryZNmsjHx0cjRoxQXl6ezX7Z2dkaMmSIvLy8FBQUpOnTp+vKlSt1XD0AAKiv6kXY2bdvn1555RV16tTJpj0xMVHvv/++1q9frx07dujMmTMaPny4dXt5ebmGDBmisrIy7d69W6+//rrWrFmjuXPn1vUlAACAesrpYae4uFhjxozRqlWrdNNNN1nbCwsL9dprr+nFF1/Ufffdp+7du2v16tXavXu3Pv30U0nSv//9b3355Zf617/+pS5dumjw4MF66qmntGzZMpWVlTnrkgAAQD3i9LATHx+vIUOGKCoqyqY9MzNTly9ftmlv3769wsLClJGRIUnKyMhQx44dFRwcbO0TExOjoqIiHT58+LrnLC0tVVFRkc0LAACYk5szT/7mm29q//792rdvX6Vtubm58vDwUEBAgE17cHCwcnNzrX1+GXSubr+67XqSk5M1f/78GlYPAAAaAqeN7Jw+fVpTp07V2rVr1ahRozo996xZs1RYWGh9nT59uk7PDwAA6o7Twk5mZqby8/PVrVs3ubm5yc3NTTt27NCSJUvk5uam4OBglZWVqaCgwGa/vLw8hYSESJJCQkIqPZ119f3VPtfi6ekpPz8/mxcAADAnp4WdAQMG6NChQ8rKyrK+evTooTFjxlh/dnd3V3p6unWfY8eOKTs7W5GRkZKkyMhIHTp0SPn5+dY+aWlp8vPzU0RERJ1fEwAAqH+cNmfH19dXd9xxh02bt7e3mjRpYm2Pi4tTUlKSAgMD5efnp8mTJysyMlK9e/eWJEVHRysiIkJjx45VSkqKcnNzNXv2bMXHx8vT07POrwkAANQ/Tp2g/FsWLVokFxcXjRgxQqWlpYqJidHy5cut211dXZWamqpJkyYpMjJS3t7eio2N1ZNPPunEqgEAQH1iMQzDcHYRzlZUVCR/f38VFhYyfweoRS1nbnJ2CXXq1IIhzi4BMLWq/v12+jo7AAAAtYmwAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATI2wAwAATM0hYaegoMARhwEAAHA4u8POc889p7feesv6/sEHH1STJk10yy236ODBgw4tDgAAoKbsDjsrV65UixYtJElpaWlKS0vT5s2bNXjwYE2fPt3hBQIAANSEm7075ObmWsNOamqqHnzwQUVHR6tly5bq1auXwwsEAACoCbtHdm666SadPn1akrRlyxZFRUVJkgzDUHl5uWOrAwAAqCG7R3aGDx+uhx9+WG3bttWPP/6owYMHS5IOHDigNm3aOLxAAACAmrA77CxatEgtW7bU6dOnlZKSIh8fH0lSTk6O/vznPzu8QAAAgJqwO+y4u7vr8ccfr9SemJjokIIAAAAcqVrr7Pz3f/+3+vbtq9DQUH377beSpMWLF+u9995zaHEAAAA1ZXfYWbFihZKSkjR48GAVFBRYJyUHBARo8eLFjq4PAACgRuwOO0uXLtWqVav0t7/9Ta6urtb2Hj166NChQw4tDgAAoKbsDjsnT55U165dK7V7enqqpKTEIUUBAAA4it1hp1WrVsrKyqrUvmXLFnXo0MERNQEAADiM3U9jJSUlKT4+XpcuXZJhGNq7d6/eeOMNJScn6+9//3tt1AgAAFBtdoedP/3pT2rcuLFmz56tixcv6uGHH1ZoaKheeukljR49ujZqBAAAqDa7w44kjRkzRmPGjNHFixdVXFysoKAgR9cFAADgENUKO1d5eXnJy8vLUbUAAAA4nN0TlPPy8jR27FiFhobKzc1Nrq6uNi8AAID6xO6wM378eO3fv19z5szR22+/rXfffdfmZY8VK1aoU6dO8vPzk5+fnyIjI7V582br9kuXLik+Pl5NmjSRj4+PRowYoby8PJtjZGdna8iQIfLy8lJQUJCmT5+uK1eu2HtZAADApOy+jfXJJ5/o448/VpcuXWp88ubNm2vBggVq27atDMPQ66+/rqFDh+rAgQO6/fbblZiYqE2bNmn9+vXy9/dXQkKChg8frl27dkmSysvLNWTIEIWEhGj37t3KycnRuHHj5O7urmeffbbG9QEAgIbPYhiGYc8OERERWrt27TUXFnSEwMBAPf/88xo5cqSaNm2qdevWaeTIkZKko0ePqkOHDsrIyFDv3r21efNm3X///Tpz5oyCg4MlSStXrtSMGTN09uxZeXh4VOmcRUVF8vf3V2Fhofz8/GrlugBILWducnYJderUgiHOLgEwtar+/bb7NtbixYs1c+ZMnTp1qib1VVJeXq4333xTJSUlioyMVGZmpi5fvqyoqChrn/bt2yssLEwZGRmSpIyMDHXs2NEadCQpJiZGRUVFOnz48HXPVVpaqqKiIpsXAAAwJ7tvY40aNUoXL15U69at5eXlJXd3d5vt586ds+t4hw4dUmRkpC5duiQfHx9t2LBBERERysrKkoeHhwICAmz6BwcHKzc3V5KUm5trE3Subr+67XqSk5M1f/58u+oEAAANk91hx9HfbN6uXTtlZWWpsLBQb7/9tmJjY7Vjxw6HnuM/zZo1S0lJSdb3RUVFatGiRa2eEwAAOIfdYSc2NtahBXh4eKhNmzaSpO7du2vfvn166aWXNGrUKJWVlamgoMBmdCcvL08hISGSpJCQEO3du9fmeFef1rra51o8PT3l6enp0OsAAAD1k91zdiTp66+/1uzZs/XQQw8pPz9fkrR58+ZfnSdTVRUVFSotLVX37t3l7u6u9PR067Zjx44pOztbkZGRkqTIyEgdOnTIWoMkpaWlyc/PTxERETWuBQAANHx2h50dO3aoY8eO2rNnj959910VFxdLkg4ePKh58+bZdaxZs2Zp586dOnXqlA4dOqRZs2Zp+/btGjNmjPz9/RUXF6ekpCRt27ZNmZmZmjBhgiIjI9W7d29JUnR0tCIiIjR27FgdPHhQW7du1ezZsxUfH8/IDQAAkFSN21gzZ87U008/raSkJPn6+lrb77vvPr388st2HSs/P1/jxo1TTk6O/P391alTJ23dulUDBw6UJC1atEguLi4aMWKESktLFRMTo+XLl1v3d3V1VWpqqiZNmqTIyEh5e3srNjZWTz75pL2XBQAATMrudXZ8fHx06NAhtWrVSr6+vjp48KBuvfVWnTp1Su3bt9elS5dqq9Zawzo7QN1gnR0AjlRr6+wEBAQoJyenUvuBAwd0yy232Hs4AACAWmV32Bk9erRmzJih3NxcWSwWVVRUaNeuXXr88cc1bty42qgRAACg2uyes/Pss88qPj5eLVq0UHl5uSIiIlReXq6HH35Ys2fPro0aAaBBqsltO26BAY5jd9jx8PDQqlWrNGfOHH3xxRcqLi5W165d1bZt29qoDwAAoEbsDjtXhYWFKSwszJG1AAAAOJzdYeeXX7PwSxaLRY0aNVKbNm00dOhQBQYG1rg4APXPjfZEFYCGz+6wc+DAAe3fv1/l5eVq166dJOmrr76Sq6ur2rdvr+XLl+uxxx7TJ598wirGAADA6ex+Gmvo0KGKiorSmTNnlJmZqczMTH333XcaOHCgHnroIX3//fe65557lJiYWBv1AgAA2MXuRQVvueUWpaWlVRq1OXz4sKKjo/X9999r//79io6O1g8//ODQYmsLiwoCVcdtrLrB01jAb6u1RQULCwttvnjzqrNnz6qoqEjSzwsPlpWV2XtoAAAAh6vWbaw//vGP2rBhg7777jt999132rBhg+Li4jRs2DBJ0t69e3Xbbbc5ulYAAAC72T1B+ZVXXlFiYqJGjx6tK1eu/HwQNzfFxsZq0aJFkqT27dvr73//u2MrBQAAqAa7w46Pj49WrVqlRYsW6ZtvvpEk3XrrrfLx8bH26dKli8MKBAAAqIlqLyro4+OjTp06ObIWAAAAh7N7zg4AAEBDQtgBAACmRtgBAACmVqWw061bN50/f16S9OSTT+rixYu1WhQAAICjVCnsHDlyRCUlJZKk+fPnq7i4uFaLAgAAcJQqPY3VpUsXTZgwQX379pVhGHrhhRdsHjX/pblz5zq0QAAAgJqoUthZs2aN5s2bp9TUVFksFm3evFlubpV3tVgshB0AAFCvVCnstGvXTm+++aYkycXFRenp6QoKCqrVwgAAABzB7kUFKyoqaqMOAACAWlGtFZS//vprLV68WEeOHJEkRUREaOrUqWrdurVDiwMAAKgpu9fZ2bp1qyIiIrR371516tRJnTp10p49e3T77bcrLS2tNmoEAACoNrtHdmbOnKnExEQtWLCgUvuMGTM0cOBAhxUHAABQU3aP7Bw5ckRxcXGV2v/4xz/qyy+/dEhRAAAAjmJ32GnatKmysrIqtWdlZfGEFgAAqHfsvo01ceJEPfLII/rmm2/Up08fSdKuXbv03HPPKSkpyeEFAgAA1ITdYWfOnDny9fXVwoULNWvWLElSaGionnjiCU2ZMsXhBQIAANSE3WHHYrEoMTFRiYmJunDhgiTJ19fX4YUBAAA4QrXW2bmKkAMAAOo7uycoAwAANCSEHQAAYGqEHQAAYGp2hZ3Lly9rwIABOn78eG3VAwAA4FB2hR13d3d9/vnntVULAACAw9l9G+sPf/iDXnvttdqoBQAAwOHsfvT8ypUr+sc//qEPP/xQ3bt3l7e3t832F1980WHFAQAA1JTdYeeLL75Qt27dJElfffWVzTaLxeKYqgAAABzE7rCzbdu22qgDAACgVlT70fMTJ05o69at+umnnyRJhmE4rCgAAABHsTvs/PjjjxowYIBuu+02/e53v1NOTo4kKS4uTo899pjDCwQAAKgJu8NOYmKi3N3dlZ2dLS8vL2v7qFGjtGXLFocWBwAAUFN2z9n597//ra1bt6p58+Y27W3bttW3337rsMIAAAAcwe6RnZKSEpsRnavOnTsnT09PhxQFAADgKHaHnbvvvlv//Oc/re8tFosqKiqUkpKi/v37O7Q4AACAmrL7NlZKSooGDBigzz77TGVlZfrLX/6iw4cP69y5c9q1a1dt1AgAAFBtdo/s3HHHHfrqq6/Ut29fDR06VCUlJRo+fLgOHDig1q1b10aNAAAA1Wb3yI4k+fv7629/+5ujawEAAHC4aoWd8+fP67XXXtORI0ckSREREZowYYICAwMdWhwAAEBN2X0ba+fOnWrZsqWWLFmi8+fP6/z581qyZIlatWqlnTt31kaNAAAA1Wb3yE58fLxGjRqlFStWyNXVVZJUXl6uP//5z4qPj9ehQ4ccXiQAAEB12T2yc+LECT322GPWoCNJrq6uSkpK0okTJxxaHAAAQE3ZHXa6detmnavzS0eOHFHnzp0dUhQAAICjVOk21ueff279ecqUKZo6dapOnDih3r17S5I+/fRTLVu2TAsWLKidKgEAAKrJYhiG8VudXFxcZLFY9FtdLRaLysvLHVZcXSkqKpK/v78KCwvl5+fn7HKAeq3lzE3OLuGGcGrBEGeXANR7Vf37XaWRnZMnTzqsMAAAgLpUpbATHh5e23UAAADUimotKnjmzBl98sknys/PV0VFhc22KVOmOKQwAAAAR7A77KxZs0aPPvqoPDw81KRJE1ksFus2i8VC2AEAAPWK3WFnzpw5mjt3rmbNmiUXF7ufXAcAAKhTdqeVixcvavTo0QQdAADQINidWOLi4rR+/XqHnDw5OVk9e/aUr6+vgoKCNGzYMB07dsymz6VLlxQfH68mTZrIx8dHI0aMUF5enk2f7OxsDRkyRF5eXgoKCtL06dN15coVh9QIAAAaNrtvYyUnJ+v+++/Xli1b1LFjR7m7u9tsf/HFF6t8rB07dig+Pl49e/bUlStX9Ne//lXR0dH68ssv5e3tLUlKTEzUpk2btH79evn7+yshIUHDhw/Xrl27JP38vVxDhgxRSEiIdu/erZycHI0bN07u7u569tln7b08AABgMlVaVPCXnn76ac2dO1ft2rVTcHBwpQnKH330UbWLOXv2rIKCgrRjxw7dc889KiwsVNOmTbVu3TqNHDlSknT06FF16NBBGRkZ6t27tzZv3qz7779fZ86cUXBwsCRp5cqVmjFjhs6ePSsPD4/fPC+LCgJVx6KCdYNFBYHf5tBFBX9p4cKF+sc//qHx48fXpL5rKiwslCQFBgZKkjIzM3X58mVFRUVZ+7Rv315hYWHWsJORkaGOHTtag44kxcTEaNKkSTp8+LC6du1a6TylpaUqLS21vi8qKnL4tQAAgPrB7jk7np6euuuuuxxeSEVFhaZNm6a77rpLd9xxhyQpNzdXHh4eCggIsOkbHBys3Nxca59fBp2r269uu5bk5GT5+/tbXy1atHDw1QAAgPrC7rAzdepULV261OGFxMfH64svvtCbb77p8GP/p1mzZqmwsND6On36dK2fEwAAOIfdt7H27t2rjz76SKmpqbr99tsrTVB+99137S4iISFBqamp2rlzp5o3b25tDwkJUVlZmQoKCmxGd/Ly8hQSEmLts3fvXpvjXX1a62qf/+Tp6SlPT0+76wQAAA2P3SM7AQEBGj58uO69917dfPPNNreD/P397TqWYRhKSEjQhg0b9NFHH6lVq1Y227t37y53d3elp6db244dO6bs7GxFRkZKkiIjI3Xo0CHl5+db+6SlpcnPz08RERH2Xh4AADAZu0d2Vq9e7bCTx8fHa926dXrvvffk6+trnWPj7++vxo0by9/fX3FxcUpKSlJgYKD8/Pw0efJkRUZGqnfv3pKk6OhoRUREaOzYsUpJSVFubq5mz56t+Ph4Rm8AAED1vgjUUVasWCFJ6tevn0376tWrrU97LVq0SC4uLhoxYoRKS0sVExOj5cuXW/u6uroqNTVVkyZNUmRkpLy9vRUbG6snn3yyri4DAADUY3avs9OqVSubtXX+0zfffFPjouoa6+wAVcc6O3WDdXaA31Zr6+xMmzbN5v3ly5d14MABbdmyRdOnT7e7UAAAgNpkd9iZOnXqNduXLVumzz77rMYFAQAAOJLDvrp88ODBeueddxx1OAAAAIdwWNh5++23rV/zAAAAUF/YfRura9euNhOUDcNQbm6uzp49a/OUFAAAQH1gd9gZNmyYzXsXFxc1bdpU/fr1U/v27R1VFwAAgEPYHXbmzZtXG3UAAADUCofN2QEAAKiPqjyy4+Li8quLCUqSxWLRlStXalwUAACAo1Q57GzYsOG62zIyMrRkyRJVVFQ4pCgAAABHqXLYGTp0aKW2Y8eOaebMmXr//fc1ZswYvo8KAADUO9Was3PmzBlNnDhRHTt21JUrV5SVlaXXX39d4eHhjq4PAACgRuwKO4WFhZoxY4batGmjw4cPKz09Xe+//77uuOOO2qoPAACgRqp8GyslJUXPPfecQkJC9MYbb1zzthYAAEB9YzEMw6hKRxcXFzVu3FhRUVFydXW9br93333XYcXVlap+RTwAqeXMTc4u4YZwasEQZ5cA1HtV/ftd5ZGdcePG/eaj5wAAAPVNlcPOmjVrarEMAACA2sEKygAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNQIOwAAwNTcnF0AAKCyljM3VXvfUwuGOLASoOFjZAcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJgaYQcAAJiaU8POzp079cADDyg0NFQWi0UbN2602W4YhubOnatmzZqpcePGioqK0vHjx236nDt3TmPGjJGfn58CAgIUFxen4uLiOrwKAABQnzk17JSUlKhz585atmzZNbenpKRoyZIlWrlypfbs2SNvb2/FxMTo0qVL1j5jxozR4cOHlZaWptTUVO3cuVOPPPJIXV0CAACo59ycefLBgwdr8ODB19xmGIYWL16s2bNna+jQoZKkf/7znwoODtbGjRs1evRoHTlyRFu2bNG+ffvUo0cPSdLSpUv1u9/9Ti+88IJCQ0Pr7FoAAED9VG/n7Jw8eVK5ubmKioqytvn7+6tXr17KyMiQJGVkZCggIMAadCQpKipKLi4u2rNnz3WPXVpaqqKiIpsXAAAwJ6eO7Pya3NxcSVJwcLBNe3BwsHVbbm6ugoKCbLa7ubkpMDDQ2udakpOTNX/+fAdXDDQcLWducnYJAFBn6u3ITm2aNWuWCgsLra/Tp087uyQAAFBL6m3YCQkJkSTl5eXZtOfl5Vm3hYSEKD8/32b7lStXdO7cOWufa/H09JSfn5/NCwAAmFO9DTutWrVSSEiI0tPTrW1FRUXas2ePIiMjJUmRkZEqKChQZmamtc9HH32kiooK9erVq85rBgAA9Y9T5+wUFxfrxIkT1vcnT55UVlaWAgMDFRYWpmnTpunpp59W27Zt1apVK82ZM0ehoaEaNmyYJKlDhw4aNGiQJk6cqJUrV+ry5ctKSEjQ6NGjeRILAABIcnLY+eyzz9S/f3/r+6SkJElSbGys1qxZo7/85S8qKSnRI488ooKCAvXt21dbtmxRo0aNrPusXbtWCQkJGjBggFxcXDRixAgtWbKkzq8FAADUTxbDMAxnF+FsRUVF8vf3V2FhIfN3cEPgaSxzO7VgiLNLAOpEVf9+19s5OwAAAI5A2AEAAKZG2AEAAKZG2AEAAKZG2AEAAKZWb78bCwBQPTV52o4nuWBGjOwAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTI+wAAABTc3N2AQCA+qPlzE3V3vfUgiEOrARwHEZ2AACAqRF2AACAqRF2AACAqTFnB2igajK3AgBuJIzsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAUyPsAAAAU2OdHQCAQ/C9WqivGNkBAACmRtgBAACmRtgBAACmRtgBAACmxgRlAIDTMbkZtYmRHQAAYGqEHQAAYGrcxgKcpCbD9gCAqmNkBwAAmBphBwAAmBphBwAAmBphBwAAmBoTlIEaYJIxANR/hB3c8AgswI2LxQxvDKYJO8uWLdPzzz+v3Nxcde7cWUuXLtWdd97p7LIAACZFUGo4TDFn56233lJSUpLmzZun/fv3q3PnzoqJiVF+fr6zSwMAAE5mMQzDcHYRNdWrVy/17NlTL7/8siSpoqJCLVq00OTJkzVz5szf3L+oqEj+/v4qLCyUn59fbZdrag3x/3S4jQWgIWFU6P9U9e93g7+NVVZWpszMTM2aNcva5uLioqioKGVkZFxzn9LSUpWWllrfFxYWSvr5l+Zod8zbWu19v5gf48BKqq4mNddEWOJ6p5wXABqShvhvZW39Pbv6d/u3xm0afNj54YcfVF5eruDgYJv24OBgHT169Jr7JCcna/78+ZXaW7RoUSs1Vpf/YmdXAABAzdX237MLFy7I39//utsbfNipjlmzZikpKcn6vqKiQufOnVOTJk1ksVicWFnVFBUVqUWLFjp9+jS33eoZPpv6jc+n/uKzqb/q82djGIYuXLig0NDQX+3X4MPOzTffLFdXV+Xl5dm05+XlKSQk5Jr7eHp6ytPT06YtICCgtkqsNX5+fvXuPzz8jM+mfuPzqb/4bOqv+vrZ/NqIzlUN/mksDw8Pde/eXenp6da2iooKpaenKzIy0omVAQCA+qDBj+xIUlJSkmJjY9WjRw/deeedWrx4sUpKSjRhwgRnlwYAAJzMFGFn1KhROnv2rObOnavc3Fx16dJFW7ZsqTRp2Sw8PT01b968Srfi4Hx8NvUbn0/9xWdTf5nhszHFOjsAAADX0+Dn7AAAAPwawg4AADA1wg4AADA1wg4AADA1wk4Dd+rUKcXFxalVq1Zq3LixWrdurXnz5qmsrMzZpUHSM888oz59+sjLy6tBLlxpJsuWLVPLli3VqFEj9erVS3v37nV2SZC0c+dOPfDAAwoNDZXFYtHGjRudXRL+V3Jysnr27ClfX18FBQVp2LBhOnbsmLPLqhbCTgN39OhRVVRU6JVXXtHhw4e1aNEirVy5Un/961+dXRr08xfV/v73v9ekSZOcXcoN7a233lJSUpLmzZun/fv3q3PnzoqJiVF+fr6zS7vhlZSUqHPnzlq2bJmzS8F/2LFjh+Lj4/Xpp58qLS1Nly9fVnR0tEpKSpxdmt149NyEnn/+ea1YsULffPONs0vB/1qzZo2mTZumgoICZ5dyQ+rVq5d69uypl19+WdLPq6y3aNFCkydP1syZM51cHa6yWCzasGGDhg0b5uxScA1nz55VUFCQduzYoXvuucfZ5diFkR0TKiwsVGBgoLPLAOqFsrIyZWZmKioqytrm4uKiqKgoZWRkOLEyoGEpLCyUpAb594WwYzInTpzQ0qVL9eijjzq7FKBe+OGHH1ReXl5pRfXg4GDl5uY6qSqgYamoqNC0adN011136Y477nB2OXYj7NRTM2fOlMVi+dXX0aNHbfb5/vvvNWjQIP3+97/XxIkTnVS5+VXnswGAhiw+Pl5ffPGF3nzzTWeXUi2m+G4sM3rsscc0fvz4X+1z6623Wn8+c+aM+vfvrz59+ujVV1+t5epubPZ+NnCum2++Wa6ursrLy7Npz8vLU0hIiJOqAhqOhIQEpaamaufOnWrevLmzy6kWwk491bRpUzVt2rRKfb///nv1799f3bt31+rVq+XiwoBdbbLns4HzeXh4qHv37kpPT7dOfK2oqFB6eroSEhKcWxxQjxmGocmTJ2vDhg3avn27WrVq5eySqo2w08B9//336tevn8LDw/XCCy/o7Nmz1m38X6vzZWdn69y5c8rOzlZ5ebmysrIkSW3atJGPj49zi7uBJCUlKTY2Vj169NCdd96pxYsXq6SkRBMmTHB2aTe84uJinThxwvr+5MmTysrKUmBgoMLCwpxYGeLj47Vu3Tq999578vX1tc5x8/f3V+PGjZ1cnX149LyBW7NmzXX/weajdb7x48fr9ddfr9S+bds29evXr+4LuoG9/PLLev7555Wbm6suXbpoyZIl6tWrl7PLuuFt375d/fv3r9QeGxurNWvW1H1BsLJYLNdsX7169W/eyq9vCDsAAMDUmNwBAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADAABMjbADwPTGjx9v/aoIZ7FYLNq4ceN1t2/fvl0Wi0UFBQXWto0bN6pNmzZydXXVtGnTar1GwKwIOwCsxo8ff81vcR80aJCzS6uRl156qd6vxtunTx/l5OTI39/f2vboo49q5MiROn36tJ566ql6EdqAhojvxgJgY9CgQVq9erVNm6enZ62es6ysTB4eHrV2/F8GiPrKw8PD5vvsiouLlZ+fr5iYGIWGhjqxMqDhY2QHgA1PT0+FhITYvG666SZJP99q8fDw0Mcff2ztn5KSoqCgIOXl5UmS+vXrp4SEBCUkJMjf318333yz5syZY/NdbS1bttRTTz2lcePGyc/PT4888ogk6ZNPPtHdd9+txo0bq0WLFpoyZYpKSkqs+y1fvlxt27ZVo0aNFBwcrJEjR1q3vf322+rYsaMaN26sJk2aKCoqyrrvf46IlJaWasqUKQoKClKjRo3Ut29f7du3z7r96i2l9PR09ejRQ15eXurTp4+OHTt23d9bWVmZEhIS1KxZMzVq1Ejh4eFKTk626fPDDz/ov/7rv+Tl5aW2bdvq//2//1fpnAUFBdq+fbt8fX0lSffdd58sFov69eun119/Xe+99551xG379u2//mEC+JkBAP8rNjbWGDp06K/2mT59uhEeHm4UFBQY+/fvNzw8PIz33nvPuv3ee+81fHx8jKlTpxpHjx41/vWvfxleXl7Gq6++au0THh5u+Pn5GS+88IJx4sQJ68vb29tYtGiR8dVXXxm7du0yunbtaowfP94wDMPYt2+f4erqaqxbt844deqUsX//fuOll14yDMMwzpw5Y7i5uRkvvviicfLkSePzzz83li1bZly4cOGa1zVlyhQjNDTU+OCDD4zDhw8bsbGxxk033WT8+OOPhmEYxrZt2wxJRq9evYzt27cbhw8fNu6++26jT58+1/29PP/880aLFi2MnTt3GqdOnTI+/vhjY926ddbtkozmzZsb69atM44fP25MmTLF8PHxqXTO8+fPG6WlpcaxY8cMScY777xj5OTkGIWFhcaDDz5oDBo0yMjJyTFycnKM0tLSKnyqAAg7AKxiY2MNV1dXw9vb2+b1zDPPWPuUlpYaXbp0MR588EEjIiLCmDhxos0x7r33XqNDhw5GRUWFtW3GjBlGhw4drO/Dw8ONYcOG2ewXFxdnPPLIIzZtH3/8seHi4mL89NNPxjvvvGP4+fkZRUVFlerOzMw0JBmnTp267nVdDTvFxcWGu7u7sXbtWuv2srIyIzQ01EhJSTEM4/+Cx4cffmjts2nTJkOS8dNPP13zHJMnTzbuu+8+m+v+JUnG7Nmzre+Li4sNScbmzZttznn+/HnDMAzj/PnzhiRj27Zt17wOAFXHnB0ANvr3768VK1bYtAUGBlp/9vDw0Nq1a9WpUyeFh4dr0aJFlY7Ru3dvWSwW6/vIyEgtXLhQ5eXlcnV1lST16NHDZp+DBw/q888/19q1a61thmGooqJCJ0+e1MCBAxUeHq5bb71VgwYN0qBBg6y3hDp37qwBAwaoY8eOiomJUXR0tEaOHGm9/fZLX3/9tS5fvqy77rrL2ubu7q4777xTR44csenbqVMn68/NmjWTJOXn5yssLKzSccePH6+BAweqXbt2GjRokO6//35FR0df93je3t7y8/NTfn5+pWMBcCzm7ACw4e3trTZt2ti8fhl2JGn37t2SpHPnzuncuXPVPs8vFRcX69FHH1VWVpb1dfDgQR0/flytW7eWr6+v9u/frzfeeEPNmjXT3Llz1blzZxUUFMjV1VVpaWnavHmzIiIitHTpUrVr104nT56s3i/hf7m7u1t/vhreKioqrtm3W7duOnnypJ566in99NNPevDBB23mFP3n8a4e83rHA+A4hB0Advn666+VmJioVatWqVevXoqNja30B3vPnj027z/99FO1bdvWOqpzLd26ddOXX35ZKWi1adPG+qSWm5uboqKilJKSos8//1ynTp3SRx99JOnn4HDXXXdp/vz5OnDggDw8PLRhw4ZK52ndurU8PDy0a9cua9vly5e1b98+RUREVPv3Ikl+fn4aNWqUVq1apbfeekvvvPNOtcPgtXh4eKi8vNxhxwNuFNzGAmCjtLRUubm5Nm1ubm66+eabVV5erj/84Q+KiYnRhAkTNGjQIHXs2FELFy7U9OnTrf2zs7OVlJSkRx99VPv379fSpUu1cOHCXz3vjBkz1Lt3byUkJOhPf/qTvL299eWXXyotLU0vv/yyUlNT9c033+iee+7RTTfdpA8++EAVFRVq166d9uzZo/T0dEVHRysoKEh79uzR2bNn1aFDh0rn8fb21qRJkzR9+nQFBgYqLCxMKSkpunjxouLi4qr9e3vxxRfVrFkzde3aVS4uLlq/fr1CQkIUEBBQ7WP+p5YtW2rr1q06duyYmjRpIn9//0qjRQAqI+wAsLFlyxbr/JSr2rVrp6NHj+qZZ57Rt99+q9TUVEk/z2N59dVX9dBDDyk6OlqdO3eWJI0bN04//fST7rzzTrm6umrq1KnWx8uvp1OnTtqxY4f+9re/6e6775ZhGGrdurVGjRolSQoICNC7776rJ554QpcuXVLbtm31xhtv6Pbbb9eRI0e0c+dOLV68WEVFRQoPD9fChQs1ePDga55rwYIFqqio0NixY3XhwgX16NFDW7duveYcn6ry9fVVSkqKjh8/LldXV/Xs2VMffPCBXFwcN4A+ceJEbd++XT169FBxcbG2bdumfv36Oez4gFlZDOMXi18AQA3169dPXbp00eLFi51dCgBIYs4OAAAwOcIOAAAwNW5jAQAAU2NkBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmBphBwAAmNr/B7gwwwn4LHZpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "adata = adata_ref[adata_ref.obs['cell_type'].isin(\n",
    "    ['Epiblast', 'Nascent mesoderm'])].copy()\n",
    "adata.X = adata.X.toarray()\n",
    "sc.pp.normalize_total(adata, target_sum=1e4)\n",
    "sc.pp.log1p(adata)\n",
    "\n",
    "source_mean_exp = adata.X[adata.obs['cell_type'] == 'Epiblast'].mean(axis=0)\n",
    "target_mean_exp = adata.X[adata.obs['cell_type'] == 'Nascent mesoderm'].mean(axis=0)\n",
    "\n",
    "exp_shift = target_mean_exp - source_mean_exp\n",
    "\n",
    "# Plot the distribution of expression shifts\n",
    "plt.hist(exp_shift, bins=30)\n",
    "plt.xlabel('Expression shift')\n",
    "plt.ylabel('Number of genes')\n",
    "plt.ylim(0, 500)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cell state transitions typically involve expression changes of hundreds to thousands of genes. Some of the genes are causal drivers for a cell state transition, while others are passengers that are dispensible for the transition. In order to predict which genes are causal, we make two assumptions:\n",
    "1. The expression shift of a causal gene should be large compared to its standard deviation in the global cell state manifold.\n",
    "2. The perturbation of a causal gene should be able to alter the expression of other genes that change through the cell state transition.\n",
    "\n",
    "Let's normalize the gene expression shifts by the global cell state manifold and match the expression shifts to the genetic perturbation database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 16080/16080 [00:02<00:00, 7757.79it/s]\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>condition</th>\n",
       "      <th>perturbed_gene</th>\n",
       "      <th>perturbation_sign</th>\n",
       "      <th>perturbed_gene_name</th>\n",
       "      <th>gene_shift_z</th>\n",
       "      <th>pert_sim</th>\n",
       "      <th>pert_match_score</th>\n",
       "      <th>causal_score</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pDS033</th>\n",
       "      <td>AdamsonWeissman2016_GSM2406681_10X010</td>\n",
       "      <td>ENSG00000106803</td>\n",
       "      <td>-1</td>\n",
       "      <td>SEC61B</td>\n",
       "      <td>0.594664</td>\n",
       "      <td>0.003110</td>\n",
       "      <td>-0.003110</td>\n",
       "      <td>-0.001850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052</th>\n",
       "      <td>AdamsonWeissman2016_GSM2406681_10X010</td>\n",
       "      <td>ENSG00000112249</td>\n",
       "      <td>-1</td>\n",
       "      <td>ASCC3</td>\n",
       "      <td>-0.281237</td>\n",
       "      <td>0.027262</td>\n",
       "      <td>-0.027262</td>\n",
       "      <td>0.007667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_pDS026</th>\n",
       "      <td>AdamsonWeissman2016_GSM2406681_10X010</td>\n",
       "      <td>ENSG00000205981</td>\n",
       "      <td>-1</td>\n",
       "      <td>DNAJC19</td>\n",
       "      <td>0.168256</td>\n",
       "      <td>-0.009787</td>\n",
       "      <td>0.009787</td>\n",
       "      <td>0.001647</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088</th>\n",
       "      <td>AdamsonWeissman2016_GSM2406681_10X010</td>\n",
       "      <td>ENSG00000113013</td>\n",
       "      <td>-1</td>\n",
       "      <td>HSPA9</td>\n",
       "      <td>-0.323847</td>\n",
       "      <td>-0.018247</td>\n",
       "      <td>0.018247</td>\n",
       "      <td>-0.005909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017</th>\n",
       "      <td>AdamsonWeissman2016_GSM2406681_10X010</td>\n",
       "      <td>ENSG00000044574</td>\n",
       "      <td>-1</td>\n",
       "      <td>HSPA5</td>\n",
       "      <td>-0.105252</td>\n",
       "      <td>0.008076</td>\n",
       "      <td>-0.008076</td>\n",
       "      <td>0.000850</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TianKampmann2021_CRISPRi_CTSD</th>\n",
       "      <td>TianKampmann2021_CRISPRi</td>\n",
       "      <td>ENSG00000117984</td>\n",
       "      <td>-1</td>\n",
       "      <td>CTSD</td>\n",
       "      <td>-0.169454</td>\n",
       "      <td>0.013613</td>\n",
       "      <td>-0.013613</td>\n",
       "      <td>0.002307</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TianKampmann2021_CRISPRi_RNF165</th>\n",
       "      <td>TianKampmann2021_CRISPRi</td>\n",
       "      <td>ENSG00000141622</td>\n",
       "      <td>-1</td>\n",
       "      <td>RNF165</td>\n",
       "      <td>-0.033471</td>\n",
       "      <td>0.003837</td>\n",
       "      <td>-0.003837</td>\n",
       "      <td>0.000128</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TianKampmann2021_CRISPRi_NDUFA2</th>\n",
       "      <td>TianKampmann2021_CRISPRi</td>\n",
       "      <td>ENSG00000131495</td>\n",
       "      <td>-1</td>\n",
       "      <td>NDUFA2</td>\n",
       "      <td>0.570147</td>\n",
       "      <td>0.015090</td>\n",
       "      <td>-0.015090</td>\n",
       "      <td>-0.008604</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TianKampmann2021_CRISPRi_KIF1B</th>\n",
       "      <td>TianKampmann2021_CRISPRi</td>\n",
       "      <td>ENSG00000054523</td>\n",
       "      <td>-1</td>\n",
       "      <td>KIF1B</td>\n",
       "      <td>-0.081678</td>\n",
       "      <td>0.005134</td>\n",
       "      <td>-0.005134</td>\n",
       "      <td>0.000419</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TianKampmann2021_CRISPRi_REEP2</th>\n",
       "      <td>TianKampmann2021_CRISPRi</td>\n",
       "      <td>ENSG00000132563</td>\n",
       "      <td>-1</td>\n",
       "      <td>REEP2</td>\n",
       "      <td>-0.031243</td>\n",
       "      <td>0.020592</td>\n",
       "      <td>-0.020592</td>\n",
       "      <td>0.000643</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>16080 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                                                condition  \\\n",
       "id                                                                                          \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...  AdamsonWeissman2016_GSM2406681_10X010   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052  AdamsonWeissman2016_GSM2406681_10X010   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...  AdamsonWeissman2016_GSM2406681_10X010   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088  AdamsonWeissman2016_GSM2406681_10X010   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017  AdamsonWeissman2016_GSM2406681_10X010   \n",
       "...                                                                                   ...   \n",
       "TianKampmann2021_CRISPRi_CTSD                                    TianKampmann2021_CRISPRi   \n",
       "TianKampmann2021_CRISPRi_RNF165                                  TianKampmann2021_CRISPRi   \n",
       "TianKampmann2021_CRISPRi_NDUFA2                                  TianKampmann2021_CRISPRi   \n",
       "TianKampmann2021_CRISPRi_KIF1B                                   TianKampmann2021_CRISPRi   \n",
       "TianKampmann2021_CRISPRi_REEP2                                   TianKampmann2021_CRISPRi   \n",
       "\n",
       "                                                     perturbed_gene  \\\n",
       "id                                                                    \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...  ENSG00000106803   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052  ENSG00000112249   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...  ENSG00000205981   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088  ENSG00000113013   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017  ENSG00000044574   \n",
       "...                                                             ...   \n",
       "TianKampmann2021_CRISPRi_CTSD                       ENSG00000117984   \n",
       "TianKampmann2021_CRISPRi_RNF165                     ENSG00000141622   \n",
       "TianKampmann2021_CRISPRi_NDUFA2                     ENSG00000131495   \n",
       "TianKampmann2021_CRISPRi_KIF1B                      ENSG00000054523   \n",
       "TianKampmann2021_CRISPRi_REEP2                      ENSG00000132563   \n",
       "\n",
       "                                                    perturbation_sign  \\\n",
       "id                                                                      \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...                 -1   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052                 -1   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...                 -1   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088                 -1   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017                 -1   \n",
       "...                                                               ...   \n",
       "TianKampmann2021_CRISPRi_CTSD                                      -1   \n",
       "TianKampmann2021_CRISPRi_RNF165                                    -1   \n",
       "TianKampmann2021_CRISPRi_NDUFA2                                    -1   \n",
       "TianKampmann2021_CRISPRi_KIF1B                                     -1   \n",
       "TianKampmann2021_CRISPRi_REEP2                                     -1   \n",
       "\n",
       "                                                   perturbed_gene_name  \\\n",
       "id                                                                       \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...              SEC61B   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052               ASCC3   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...             DNAJC19   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088               HSPA9   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017               HSPA5   \n",
       "...                                                                ...   \n",
       "TianKampmann2021_CRISPRi_CTSD                                     CTSD   \n",
       "TianKampmann2021_CRISPRi_RNF165                                 RNF165   \n",
       "TianKampmann2021_CRISPRi_NDUFA2                                 NDUFA2   \n",
       "TianKampmann2021_CRISPRi_KIF1B                                   KIF1B   \n",
       "TianKampmann2021_CRISPRi_REEP2                                   REEP2   \n",
       "\n",
       "                                                    gene_shift_z  pert_sim  \\\n",
       "id                                                                           \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...      0.594664  0.003110   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052     -0.281237  0.027262   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...      0.168256 -0.009787   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088     -0.323847 -0.018247   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017     -0.105252  0.008076   \n",
       "...                                                          ...       ...   \n",
       "TianKampmann2021_CRISPRi_CTSD                          -0.169454  0.013613   \n",
       "TianKampmann2021_CRISPRi_RNF165                        -0.033471  0.003837   \n",
       "TianKampmann2021_CRISPRi_NDUFA2                         0.570147  0.015090   \n",
       "TianKampmann2021_CRISPRi_KIF1B                         -0.081678  0.005134   \n",
       "TianKampmann2021_CRISPRi_REEP2                         -0.031243  0.020592   \n",
       "\n",
       "                                                    pert_match_score  \\\n",
       "id                                                                     \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...         -0.003110   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052         -0.027262   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...          0.009787   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088          0.018247   \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017         -0.008076   \n",
       "...                                                              ...   \n",
       "TianKampmann2021_CRISPRi_CTSD                              -0.013613   \n",
       "TianKampmann2021_CRISPRi_RNF165                            -0.003837   \n",
       "TianKampmann2021_CRISPRi_NDUFA2                            -0.015090   \n",
       "TianKampmann2021_CRISPRi_KIF1B                             -0.005134   \n",
       "TianKampmann2021_CRISPRi_REEP2                             -0.020592   \n",
       "\n",
       "                                                    causal_score  \n",
       "id                                                                \n",
       "AdamsonWeissman2016_GSM2406681_10X010_SEC61B_pD...     -0.001850  \n",
       "AdamsonWeissman2016_GSM2406681_10X010_ASCC3_pDS052      0.007667  \n",
       "AdamsonWeissman2016_GSM2406681_10X010_DNAJC19_p...      0.001647  \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA9_pDS088     -0.005909  \n",
       "AdamsonWeissman2016_GSM2406681_10X010_HSPA5_pDS017      0.000850  \n",
       "...                                                          ...  \n",
       "TianKampmann2021_CRISPRi_CTSD                           0.002307  \n",
       "TianKampmann2021_CRISPRi_RNF165                         0.000128  \n",
       "TianKampmann2021_CRISPRi_NDUFA2                        -0.008604  \n",
       "TianKampmann2021_CRISPRi_KIF1B                          0.000419  \n",
       "TianKampmann2021_CRISPRi_REEP2                          0.000643  \n",
       "\n",
       "[16080 rows x 8 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the standard deviation of the genes on the global cell state manifold\n",
    "sc.pp.normalize_total(adata_ref, target_sum=1e4)\n",
    "sc.pp.log1p(adata_ref)\n",
    "X_csc = adata_ref.X.tocsc()\n",
    "adata_ref.var['std'] = [np.std(X_csc.getcol(i).toarray()) for i in range(adata_ref.X.shape[1])]\n",
    "\n",
    "# Construct a causal gene predictor which will normalize the \n",
    "# gene expression shifts and match the shift to perturbations\n",
    "causal_gene_predictor = CausalGenePredictor(adata_pert, adata_ref.var['std'].values)\n",
    "\n",
    "# Calculate the match scores\n",
    "pert_match_df = causal_gene_predictor.calc_causal_scores(exp_shift)\n",
    "pert_match_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The perturbation of a gene can have different phenotypic effects under different experimental conditions. For causal gene prediction of cell state transition, let's only keep the experimental condition that have the best causal score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>condition</th>\n",
       "      <th>perturbed_gene</th>\n",
       "      <th>perturbation_sign</th>\n",
       "      <th>perturbed_gene_name</th>\n",
       "      <th>gene_shift_z</th>\n",
       "      <th>pert_sim</th>\n",
       "      <th>pert_match_score</th>\n",
       "      <th>causal_score</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_TBXT</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000164458</td>\n",
       "      <td>1</td>\n",
       "      <td>TBXT</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>0.198749</td>\n",
       "      <td>0.198749</td>\n",
       "      <td>0.993747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_EOMES</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000163508</td>\n",
       "      <td>1</td>\n",
       "      <td>EOMES</td>\n",
       "      <td>2.508312</td>\n",
       "      <td>0.155389</td>\n",
       "      <td>0.155389</td>\n",
       "      <td>0.389764</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PertOrg_Pertg00165</th>\n",
       "      <td>PertOrg_Pertg00165</td>\n",
       "      <td>ENSG00000142182</td>\n",
       "      <td>-1</td>\n",
       "      <td>DNMT3L</td>\n",
       "      <td>-4.887468</td>\n",
       "      <td>0.074330</td>\n",
       "      <td>-0.074330</td>\n",
       "      <td>0.363284</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PertOrg_Pertg01532</th>\n",
       "      <td>PertOrg_Pertg01532</td>\n",
       "      <td>ENSG00000088305</td>\n",
       "      <td>1</td>\n",
       "      <td>DNMT3B</td>\n",
       "      <td>-2.259588</td>\n",
       "      <td>-0.150434</td>\n",
       "      <td>-0.150434</td>\n",
       "      <td>0.339919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_SNAI1</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000124216</td>\n",
       "      <td>1</td>\n",
       "      <td>SNAI1</td>\n",
       "      <td>2.839346</td>\n",
       "      <td>0.107961</td>\n",
       "      <td>0.107961</td>\n",
       "      <td>0.306538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>knockTF_human_DataSet_01_45</th>\n",
       "      <td>knockTF_human_DataSet_01_45</td>\n",
       "      <td>ENSG00000147596</td>\n",
       "      <td>-1</td>\n",
       "      <td>PRDM14</td>\n",
       "      <td>-2.788014</td>\n",
       "      <td>0.106656</td>\n",
       "      <td>-0.106656</td>\n",
       "      <td>0.297358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>knockTF_human_DataSet_01_42</th>\n",
       "      <td>knockTF_human_DataSet_01_42</td>\n",
       "      <td>ENSG00000204531</td>\n",
       "      <td>-1</td>\n",
       "      <td>POU5F1</td>\n",
       "      <td>-1.897957</td>\n",
       "      <td>0.154678</td>\n",
       "      <td>-0.154678</td>\n",
       "      <td>0.293573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_EVX1</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000106038</td>\n",
       "      <td>1</td>\n",
       "      <td>EVX1</td>\n",
       "      <td>2.088400</td>\n",
       "      <td>0.130271</td>\n",
       "      <td>0.130271</td>\n",
       "      <td>0.272057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>knockTF_human_DataSet_01_94</th>\n",
       "      <td>knockTF_human_DataSet_01_94</td>\n",
       "      <td>ENSG00000111704</td>\n",
       "      <td>-1</td>\n",
       "      <td>NANOG</td>\n",
       "      <td>-2.193654</td>\n",
       "      <td>0.096152</td>\n",
       "      <td>-0.096152</td>\n",
       "      <td>0.210923</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PertOrg_Pertg01669</th>\n",
       "      <td>PertOrg_Pertg01669</td>\n",
       "      <td>ENSG00000121570</td>\n",
       "      <td>-1</td>\n",
       "      <td>DPPA4</td>\n",
       "      <td>-3.836496</td>\n",
       "      <td>0.051048</td>\n",
       "      <td>-0.051048</td>\n",
       "      <td>0.195845</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_GSC</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000133937</td>\n",
       "      <td>1</td>\n",
       "      <td>GSC</td>\n",
       "      <td>1.615573</td>\n",
       "      <td>0.112268</td>\n",
       "      <td>0.112268</td>\n",
       "      <td>0.181376</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_TBX6</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000149922</td>\n",
       "      <td>1</td>\n",
       "      <td>TBX6</td>\n",
       "      <td>2.763133</td>\n",
       "      <td>0.062253</td>\n",
       "      <td>0.062253</td>\n",
       "      <td>0.172014</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_FOXF1</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000103241</td>\n",
       "      <td>1</td>\n",
       "      <td>FOXF1</td>\n",
       "      <td>0.878242</td>\n",
       "      <td>0.185584</td>\n",
       "      <td>0.185584</td>\n",
       "      <td>0.162988</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PertOrg_Pertg09122</th>\n",
       "      <td>PertOrg_Pertg09122</td>\n",
       "      <td>ENSG00000163530</td>\n",
       "      <td>-1</td>\n",
       "      <td>DPPA2</td>\n",
       "      <td>-2.476860</td>\n",
       "      <td>0.060891</td>\n",
       "      <td>-0.060891</td>\n",
       "      <td>0.150819</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_ETV2</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000105672</td>\n",
       "      <td>1</td>\n",
       "      <td>ETV2</td>\n",
       "      <td>0.882748</td>\n",
       "      <td>0.161668</td>\n",
       "      <td>0.161668</td>\n",
       "      <td>0.142712</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hESC_TF_screen_CDX2</th>\n",
       "      <td>hESC_TF_screen</td>\n",
       "      <td>ENSG00000165556</td>\n",
       "      <td>1</td>\n",
       "      <td>CDX2</td>\n",
       "      <td>0.711745</td>\n",
       "      <td>0.199317</td>\n",
       "      <td>0.199317</td>\n",
       "      <td>0.141863</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PertOrg_Pertg09271</th>\n",
       "      <td>PertOrg_Pertg09271</td>\n",
       "      <td>ENSG00000138336</td>\n",
       "      <td>-1</td>\n",
       "      <td>TET1</td>\n",
       "      <td>-1.098587</td>\n",
       "      <td>0.114450</td>\n",
       "      <td>-0.114450</td>\n",
       "      <td>0.125733</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PertOrg_Pertg08127</th>\n",
       "      <td>PertOrg_Pertg08127</td>\n",
       "      <td>ENSG00000071655</td>\n",
       "      <td>-1</td>\n",
       "      <td>MBD3</td>\n",
       "      <td>1.182511</td>\n",
       "      <td>-0.106196</td>\n",
       "      <td>0.106196</td>\n",
       "      <td>0.125578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ReplogleWeissman2022_rpe1_XRCC6</th>\n",
       "      <td>ReplogleWeissman2022_rpe1</td>\n",
       "      <td>ENSG00000196419</td>\n",
       "      <td>-1</td>\n",
       "      <td>XRCC6</td>\n",
       "      <td>-1.479279</td>\n",
       "      <td>0.080169</td>\n",
       "      <td>-0.080169</td>\n",
       "      <td>0.118593</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ReplogleWeissman2022_rpe1_XRCC5</th>\n",
       "      <td>ReplogleWeissman2022_rpe1</td>\n",
       "      <td>ENSG00000079246</td>\n",
       "      <td>-1</td>\n",
       "      <td>XRCC5</td>\n",
       "      <td>-1.470233</td>\n",
       "      <td>0.078640</td>\n",
       "      <td>-0.078640</td>\n",
       "      <td>0.115620</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                   condition   perturbed_gene  \\\n",
       "id                                                                              \n",
       "hESC_TF_screen_TBXT                           hESC_TF_screen  ENSG00000164458   \n",
       "hESC_TF_screen_EOMES                          hESC_TF_screen  ENSG00000163508   \n",
       "PertOrg_Pertg00165                        PertOrg_Pertg00165  ENSG00000142182   \n",
       "PertOrg_Pertg01532                        PertOrg_Pertg01532  ENSG00000088305   \n",
       "hESC_TF_screen_SNAI1                          hESC_TF_screen  ENSG00000124216   \n",
       "knockTF_human_DataSet_01_45      knockTF_human_DataSet_01_45  ENSG00000147596   \n",
       "knockTF_human_DataSet_01_42      knockTF_human_DataSet_01_42  ENSG00000204531   \n",
       "hESC_TF_screen_EVX1                           hESC_TF_screen  ENSG00000106038   \n",
       "knockTF_human_DataSet_01_94      knockTF_human_DataSet_01_94  ENSG00000111704   \n",
       "PertOrg_Pertg01669                        PertOrg_Pertg01669  ENSG00000121570   \n",
       "hESC_TF_screen_GSC                            hESC_TF_screen  ENSG00000133937   \n",
       "hESC_TF_screen_TBX6                           hESC_TF_screen  ENSG00000149922   \n",
       "hESC_TF_screen_FOXF1                          hESC_TF_screen  ENSG00000103241   \n",
       "PertOrg_Pertg09122                        PertOrg_Pertg09122  ENSG00000163530   \n",
       "hESC_TF_screen_ETV2                           hESC_TF_screen  ENSG00000105672   \n",
       "hESC_TF_screen_CDX2                           hESC_TF_screen  ENSG00000165556   \n",
       "PertOrg_Pertg09271                        PertOrg_Pertg09271  ENSG00000138336   \n",
       "PertOrg_Pertg08127                        PertOrg_Pertg08127  ENSG00000071655   \n",
       "ReplogleWeissman2022_rpe1_XRCC6    ReplogleWeissman2022_rpe1  ENSG00000196419   \n",
       "ReplogleWeissman2022_rpe1_XRCC5    ReplogleWeissman2022_rpe1  ENSG00000079246   \n",
       "\n",
       "                                 perturbation_sign perturbed_gene_name  \\\n",
       "id                                                                       \n",
       "hESC_TF_screen_TBXT                              1                TBXT   \n",
       "hESC_TF_screen_EOMES                             1               EOMES   \n",
       "PertOrg_Pertg00165                              -1              DNMT3L   \n",
       "PertOrg_Pertg01532                               1              DNMT3B   \n",
       "hESC_TF_screen_SNAI1                             1               SNAI1   \n",
       "knockTF_human_DataSet_01_45                     -1              PRDM14   \n",
       "knockTF_human_DataSet_01_42                     -1              POU5F1   \n",
       "hESC_TF_screen_EVX1                              1                EVX1   \n",
       "knockTF_human_DataSet_01_94                     -1               NANOG   \n",
       "PertOrg_Pertg01669                              -1               DPPA4   \n",
       "hESC_TF_screen_GSC                               1                 GSC   \n",
       "hESC_TF_screen_TBX6                              1                TBX6   \n",
       "hESC_TF_screen_FOXF1                             1               FOXF1   \n",
       "PertOrg_Pertg09122                              -1               DPPA2   \n",
       "hESC_TF_screen_ETV2                              1                ETV2   \n",
       "hESC_TF_screen_CDX2                              1                CDX2   \n",
       "PertOrg_Pertg09271                              -1                TET1   \n",
       "PertOrg_Pertg08127                              -1                MBD3   \n",
       "ReplogleWeissman2022_rpe1_XRCC6                 -1               XRCC6   \n",
       "ReplogleWeissman2022_rpe1_XRCC5                 -1               XRCC5   \n",
       "\n",
       "                                 gene_shift_z  pert_sim  pert_match_score  \\\n",
       "id                                                                          \n",
       "hESC_TF_screen_TBXT                  5.000000  0.198749          0.198749   \n",
       "hESC_TF_screen_EOMES                 2.508312  0.155389          0.155389   \n",
       "PertOrg_Pertg00165                  -4.887468  0.074330         -0.074330   \n",
       "PertOrg_Pertg01532                  -2.259588 -0.150434         -0.150434   \n",
       "hESC_TF_screen_SNAI1                 2.839346  0.107961          0.107961   \n",
       "knockTF_human_DataSet_01_45         -2.788014  0.106656         -0.106656   \n",
       "knockTF_human_DataSet_01_42         -1.897957  0.154678         -0.154678   \n",
       "hESC_TF_screen_EVX1                  2.088400  0.130271          0.130271   \n",
       "knockTF_human_DataSet_01_94         -2.193654  0.096152         -0.096152   \n",
       "PertOrg_Pertg01669                  -3.836496  0.051048         -0.051048   \n",
       "hESC_TF_screen_GSC                   1.615573  0.112268          0.112268   \n",
       "hESC_TF_screen_TBX6                  2.763133  0.062253          0.062253   \n",
       "hESC_TF_screen_FOXF1                 0.878242  0.185584          0.185584   \n",
       "PertOrg_Pertg09122                  -2.476860  0.060891         -0.060891   \n",
       "hESC_TF_screen_ETV2                  0.882748  0.161668          0.161668   \n",
       "hESC_TF_screen_CDX2                  0.711745  0.199317          0.199317   \n",
       "PertOrg_Pertg09271                  -1.098587  0.114450         -0.114450   \n",
       "PertOrg_Pertg08127                   1.182511 -0.106196          0.106196   \n",
       "ReplogleWeissman2022_rpe1_XRCC6     -1.479279  0.080169         -0.080169   \n",
       "ReplogleWeissman2022_rpe1_XRCC5     -1.470233  0.078640         -0.078640   \n",
       "\n",
       "                                 causal_score  \n",
       "id                                             \n",
       "hESC_TF_screen_TBXT                  0.993747  \n",
       "hESC_TF_screen_EOMES                 0.389764  \n",
       "PertOrg_Pertg00165                   0.363284  \n",
       "PertOrg_Pertg01532                   0.339919  \n",
       "hESC_TF_screen_SNAI1                 0.306538  \n",
       "knockTF_human_DataSet_01_45          0.297358  \n",
       "knockTF_human_DataSet_01_42          0.293573  \n",
       "hESC_TF_screen_EVX1                  0.272057  \n",
       "knockTF_human_DataSet_01_94          0.210923  \n",
       "PertOrg_Pertg01669                   0.195845  \n",
       "hESC_TF_screen_GSC                   0.181376  \n",
       "hESC_TF_screen_TBX6                  0.172014  \n",
       "hESC_TF_screen_FOXF1                 0.162988  \n",
       "PertOrg_Pertg09122                   0.150819  \n",
       "hESC_TF_screen_ETV2                  0.142712  \n",
       "hESC_TF_screen_CDX2                  0.141863  \n",
       "PertOrg_Pertg09271                   0.125733  \n",
       "PertOrg_Pertg08127                   0.125578  \n",
       "ReplogleWeissman2022_rpe1_XRCC6      0.118593  \n",
       "ReplogleWeissman2022_rpe1_XRCC5      0.115620  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pert_match_df= pert_match_df.sort_values('causal_score', ascending=False)\n",
    "pert_match_df = pert_match_df.drop_duplicates('perturbed_gene', keep='first')\n",
    "pert_match_df.sort_values('causal_score', ascending=False)[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the normalized gene expression shifts and the perturbation match scores on a scatter plot. The genes along the diagonal and deviating from the point cloud around the origin are predicted to be causal for the cell state transition. Let's label the gene names of a few strong predictions on the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Perturbation match score')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAFzCAYAAADys0SZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABqT0lEQVR4nO3deVxUVf8H8M/MMAwDwrCIC4jgiqIIhAtqJa5YPi5ZPmaSa1aWWy6ppaKVS0Xlo5ktP5fMJR/L0KTcI8slFUXUFJdHwEAURUBhgGHm/P4Y73UuDDADd5hh+L5fL146d+7MPTPifOec7znfI2GMMRBCCCEiklq7AYQQQuwPBRdCCCGio+BCCCFEdBRcCCGEiI6CCyGEENFRcCGEECI6Ci6EEEJER8GFEEKI6Bys3QBbp9PpkJmZCVdXV0gkEms3hxBCaowxhgcPHsDHxwdSqWX6GBRcqpCZmQk/Pz9rN4MQQkR38+ZNNGvWzCLPTcGlCq6urgD0/whubm5Wbg2xVzt27MCIESOs3QxST+Tn58PPz4//fLMECi5V4IbC3NzcKLgQi3F2dqbfL1LrLDnUTwl9QgghoqPgQgghRHQUXAghhIiOggshhBDRUXAhhBAiOgouhBBCREfBhRBCiOgouBBCCBEdBRdCCCGio+BCCCFEdBRcCLEB/9xXo+eKw9h8Is3aTSFEFBRcCLEBqfcKkJGrxtqE69ZuCiGioMKVhNiAAC8X+MqVmBzZytpNIUQUFFwIsQHNPJQ4+kYfazeDENHQsBghhBDRUXAhhBAiOgouhNiIzSfSaMYYsRsUXAixEWsTrtOMMWK20NBQhIaGIigoCDKZjL89cuRIpKamCo61a9cOH3zwAf/YWbNm4Y033uBv//PPP2jatCn++OMP/jHNmzeHSqXib3/88ccmtYsS+oTYiMmRrbA24TrNGCNmSUpKAgCkpqYiNDSUv80dc3V15Y/l5eUhMDAQ/fv3BwC89957ePrpp3HgwAH069cPEyZMwKxZs/DUU0/xj9m4cSPi4uIQFxdnVrsouBBiI6Ij/BEd4W/tZhA7VlBQAMYYXF1dAQAuLi7YsGEDoqOj8dprr6GgoAAzZ84U5Vp1blhszZo1CAgIgJOTE7p164aTJ09WeO4333yDp556Ch4eHvDw8EC/fv0qPZ8QQuzNgwcPEBoaiuDgYLRo0QKvvvoqmjVrxt//5JNPIioqCjExMdi4cSOkUnHCQp0KLtu3b8fMmTMRExODM2fOICQkBFFRUbhz547R8xMSEjBq1Cj89ttvOH78OPz8/DBgwABkZGTUcssJIcQ6uGGx8+fP49atW9izZw9++eUX/v6HDx/i0KFD8Pb2xsWLF0W7bp0KLp9++ikmTZqE8ePHIygoCF9++SWcnZ2xfv16o+dv2bIFb7zxBp/I+r//+z/odDocOnSolltOCCHW5+npif79+ws+A2fPno1nnnkGP/74I958803cu3dPlGvVmeBSUlKCxMRE9OvXjz8mlUrRr18/HD9+3KTnKCwshEajgaenZ4XnFBcXIz8/X/BDSF1D05qJMcXFxTh69CjatGkDADh48CAOHTqEDz/8EBEREYiOjsbUqVNFuVadCS53796FVqtF48aNBccbN26MrKwsk55j7ty58PHxEQSospYvXw6VSsX/+Pn51ajdhFgDTWuuH0z5EsHlXLifkJAQTJw4EQAwffp0bNiwAc7OzgD0s8eSkpKwc+fOGret3swWW7FiBb7//nskJCTAycmpwvPmz58vmC2Rn59PAYbUOTStuX4w/BIRPa8PcnNzBfcHBARAq9WWexw3InPx4kW4ubnxxxUKBf7++2/BuePGjcO4cePMbludCS4NGzaETCbD7du3Bcdv376NJk2aVPrY2NhYrFixAgcPHkSnTp0qPVehUEChUNS4vYRYE01rrh9s+UtEnRkWc3R0RHh4uCARxSXnu3fvXuHjPvroI7z//vvYu3cvOnfuXBtNJYSQWhEd4Y+j8/rY5BeJOtNzAYCZM2di7Nix6Ny5M7p27YqVK1eioKAA48ePBwCMGTMGvr6+WL58OQDgww8/xKJFi7B161YEBATwuZkGDRqgQYMGVnsdhBBi7+pUcBk5ciSys7OxaNEiZGVlITQ0FHv37uWT/Onp6YIFQGvXrkVJSQleeOEFwfPExMRg8eLFtdl0Qiq0+UQajl+7C+2JNJv8BkpIdUgYY8zajbBl+fn5UKlUyMvLEyS+CBFLzxWH0a4wGZedO+HoPNowjFhebXyu1ZmcCyH2anJkK8hlUhQUl9K6FGI3KLgQYmXREf6QSSXIVWssti6FFlWS2kbBhRAbEODlAl93pcWmlNKiSlLb6lRCnxB71cxDiaNvWC7fYsvrIYh9ouBCSD1AiypJbaNhMUKsbNq2szh06TambTtr7aYQIhoKLoRYWXxyJtijPwmxFxRcCLGyQZ18IHn0JyH2goILIVbWtYUnFHIZuraoeJ8hQuoaCi6EWNnahOso0mhpmjCxKxRcCLEyWqFP7BEFF0JsgEarQ65ag9h9KdZuCiGioOBCiJXRcBixRxRcCLGyyZGtIJNIIAHwdFtvazfHZgUEBCAwMFCwH/z58+dRUlKCuXPnonXr1mjfvj2Cg4Px7bff8o9LTU2FRCLB0KFDBc8XExMDiUSCuLg4AMDixYvh7e0teP6lS5cCAG7evIkhQ4YgODgYwcHBCA0NxeHDh2vttddFtEKfECuLjvDH8QNSsGIgMe2+tZtj07Zv347Q0FDBsZdeegnFxcU4d+4cXFxckJqaimeeeQalpaWYOHEiAEClUuHKlSu4ffs2GjduDJ1Oh23btiE4OFjwXKNHj8bKlSvLXXfy5Mno27cvdu/eDQC4e/cuCgsLLfIa7QX1XAgRUXWqDw9Z/SeKNFoAgJeLo6WaZpeuXr2KuLg4fP3113BxcQGg7+F88sknWLJkieDc6OhobNq0CQBw8OBBhIWFwdPTtOnf//zzD3x9ffnbDRs2RPPmzUV6FfaJggshIqpO9eHkjDz+7xcz9X8Xu0S+vZTcHzlypGDY6tixY2jTpg28vLwE53Xv3h03b95EdnY2f2zs2LH8cNn69esxYcKEcs+/ZcsWwfNv374dADB37lxMnDgRPXv2xKxZs3DkyBELvkr7QMGFEBFNjmxldun8Tr4q/u8dfFToueIwYveliFoi315K7m/fvh1JSUn8j1KpNPmxzZo1Q7NmzbBnzx4kJiaif//+5c4ZPXq04PlHjhwJABg1ahTS09Mxa9YsAMDQoUPx8ccfi/Oi7BQFF0JEFB3hj6Pz+phVgfjfXfzgJJfhg2Edca+gBBm5agAQdX+X6gS9uiAsLAxXr17FvXv3BMePHz8OPz8/eHsLJ0iMHz8e48ePx4svvgip1LyPPw8PDwwfPhyffPIJ1q5di++++67G7bdnlNAnxMpi96UgXKvF0vi/oXCQwV0px+yoQEGA2nwijd+PpTql8+215H6bNm0wePBgvPrqq/juu+/g7OyM1NRUzJo1CwsXLix3/rBhw5Camoro6GizrrNnzx706dMHzs7OYIzh7NmzaNXKvgK12Ci4EGIjijQ6qDU6yCTl7zMc1rLHIFEZLrAWFJdi5MiRgqGwzz77DJs2bcKCBQsQHBwMR0dHyGQyzJkzx2hORaFQYO7cuRVea8uWLUhISOBv9+7dG5999hl+//13zJkzBw4ODmCMITAwEJ9//rmor9PeSBhjzNqNsGX5+flQqVTIy8uDm5ubtZtD7NCQ1X/CO/s0LjsHIyuvCFqmHxI7Ou/xzpQ17bnUZT1XHEZGrrrce0KqrzY+1yjnQoiVnX80WywztwhLhnY0mhupTi7HXthrvsje0bAYIVa0+UQauKEDB6nEbnMjNUHvSd1EPRdCrMhwarCLQv9dz17WpJD6jYILIVY0ObIV3JVyyGVSzI4KBPA4eR+z60K5AEOBh9QVFFwIsaLoCH/MjgqETPp4ili4vwcAQMvKV0y2l8WQxP5RcCHEyridKLmeCle8UiZBuSQ2JbdJXUEJfUKsLNzfAwWXHvdUJke2qnDaMSW3SV1BPRdCrKxsT6W+TDs2tj9L7969jS5ODAkJwc6dO7Fjxw6EhoZCo9EAALRaLSIiIrBhwwYA+oKUwcHBcHBwMFo6n9QeCi6EWFm4vwckAAZ18rH7gFJW2UKUU6ZM4QMF5/Tp07h16xYGDx6MESNGoF27dnw5/RUrVqBhw4YYP348ACA8PBz//e9/8dJLL9X6ayFCNCxGiJUlpt1HOwC7z2UCAFaNCrNug6xoyJAhmDx5MpKTk9GpUycA+t7ImDFjIJfLAQBffPEFQkJC4O/vj9WrV+Ps2bP840NCQgDA7KKURHz0L0CIlRkm5+OTM63YktpXdn+W0tJSvPzyy1i/fj0AoKioCNu2beN3lAQAT09PfPTRR3j11Vfx4YcfomnTptZqPqlEtYJLbm4u/u///g/z589HTk4OAODMmTPIyMgQtXGE1AfREf5o4uYEmUQ/NFafGNufZeLEidiyZQtKSkqwc+dOtG/fHu3btxc8Li4uDs2aNRP0WohtMXtYLDk5Gf369YNKpUJqaiomTZoET09P7Ny5E+np6fw2ooQQ03X0VeH67EHWboZNCAoKQuvWrfHzzz9j/fr1gl4LAHz//fe4evUqzpw5g/DwcDz//PN46qmnrNRaUhGzey4zZ87EuHHjcPXqVTg5OfHHn332Wdr6k5Bq+ue+mlbeG5g4cSKWLVuGkydP8rtBAsCtW7cwc+ZMfPvtt/D29saXX36JCRMmoLCw0IqtJcaYHVxOnTqF1157rdxxX19fZGVlidIoQuqb69kPkZGrRuy+FGs3xeK4Ejbc/iyGOZfffvsNgD4Xk5KSghEjRqBBgwb8YydNmoQpU6YgODgYgP5L7VNPPYV58+YBADZu3IhmzZphx44dWLx4MQ2dWZHZw2IKhQL5+fnljl+5cqXclqKEkKpN23YWGq0OAJCn1mDziTS7npLMlbBpO31ThfuzuLq64uHDh+WO79mzp9wxLvkPAOPGjcO4ceNEayupPrN7LkOGDMF7773HL2KSSCRIT0/H3Llz8fzzz4veQELs3c/nHs8QYyhfT8zeUAmb+sHs4PLJJ5/g4cOHaNSoEdRqNXr16oXWrVvD1dUVS5cutUQbCbFrTnIZAEAuk9SLD936UoGgvjN7WEylUuHAgQM4evQozp07h4cPH+KJJ55Av379LNE+QuxGRVsV9w9qjIJLl+EglcDLxRExuy7g5I2cer2YktR9ZgUXjUYDpVKJpKQk9OzZEz179rRUuwixO4bl8g2DC7dCX63RIfnRlsfxyZkUXEidZtawmFwuR/PmzaHVai3VHkLsVkW5hsmRrSCV6PdzcVc61MvFlMT+mJ1zeffdd/HOO+/wK/MJITWnYwwA8KCoFNeXD6JeSw3s3LkT4eHhCA0NRbt27dCnTx/odDpERkbCy8sLeXl5/LkvvPACNm7cKHh8TEwMZDIZ0tKEa44iIyMRFxcHAIiPj0d4eDgUCgVmzJhh4VdUN5mdc/n8889x7do1+Pj4wN/fHy4uLoL7z5w5I1rjCLEnZXeR5PIvsftSEP7oHOqx1MytW7fw6quvIjExEf7++qHHM2fOQPKoZ+jm5oYVK1Zg+fLlRh+v0+mwceNGREZGYsOGDVi8eLHR89q0aYP169djx44dRqdMk2oEl2HDhlmgGYTYP8NNwIxtVywB0LWFp/UaaAdu374NmUwGT8/H7+MTTzzB/33u3LlYtGgRpk6dCh+f8oH8wIEDaNy4MWJjYzFs2DAsWrTIaIXltm3bAgB++uknC7wK+2B2cImJibFEOwixe2V3keQCzX9P3QSy9WtcYvelmDRF13DmmeFz1ffpvZ06dcKTTz4Jf39/9OrVCz169MBLL70EX19fAECTJk3w2muvISYmBt988025x69btw4TJkxAWFgYvLy8cPDgQQwYMKC2X4ZdqHbJ/cTERGzevBmbN2+m8gqEmIlb6wGAnyFWGa5kCld7zLDnY6wXVBvKtskWSKVS/Pjjjzh27BgGDhyIo0ePokOHDrh27Rp/zpw5c7Bnzx5cvnxZ8Nh79+5h//79GDVqFABgwoQJWLduXa22356Y3XO5c+cOXnzxRSQkJMDd3R2AvgR/79698f3331MJGELMYBgQ5FIJZkcFVnie4TRmbmgt3N8DR65kw10pr/XFlxVNrbYF7dq1Q7t27fDaa69h4MCB2L17N3+fm5sb5s6di/nz50Mmk/HHv/vuO5SWlvIbjmm1Wty7dw/37t2Dl5dXrb+Gus7snsvUqVPx4MEDXLx4ETk5OcjJycGFCxeQn5+PadOmWaKNhNgtw4DgIJNU+CEd7u8BmUT/J/C455OYdh+5ag1cFA61/gFvi2VcMjIycPToUf72/fv3cePGDbRqVWb69+TJSEpKQmJiIn9s3bp1+OGHH5CamorU1FTcvHkTgwcPxubNm2ut/fbE7OCyd+9efPHFF4LNe4KCgrBmzRr8+uuvojaOEHsXHeEPyaO/F2l0FZ6XmHYfWqb/05A1P+BtqYwLN0S342Qq3nvvPbRt2xahoaF46qmnMHbsWAwdOlRwvkKhwHvvvYfU1FQAwMmTJ3Hnzp1ylUZGjx5tdGjs0KFDaNasGT799FOsW7cOzZo1E/SOCCBh7NEEexO5urrijz/+QGhoqOD42bNn0atXL6MVk+uy/Px8qFQq5OXlwc3NzdrNIXZoyOo/4Z19GkdK2yJmSAejH9YVlY4RU21cw1J6rjiMjFw1fN2VFVZaJo/Vxuea2T2XPn36YPr06cjMfFzJNSMjA2+99Rb69u0rauOMWbNmDQICAuDk5IRu3brh5MmTFZ578eJFPP/88wgICIBEIsHKlSst3j5CzHWvoAQAoNExQQ7GMGFuai+hJkl2a00MEIMtDtHVd2YHl88//xz5+fkICAhAq1at0KpVK7Ro0QL5+flYvXq1JdrI2759O2bOnImYmBicOXMGISEhiIqKwp07d4yeX1hYiJYtW2LFihVo0qSJRdtG6h5rz3birh/u7wG5TAqlXIqC4lKjM8JMVZMAUZc/oG1piI7omT0sBgCMMRw8eJCfyte+fftaqYrcrVs3dOnSBZ9//jkA/WpaPz8/TJ06ld+JriIBAQGYMWOG2aUaaFjMfll7KIW7vlIuQ0/JZRxlgVBrdHx7qjNMVZeHtkjtqY3PNbOnIgP6DcL69++P/v37i92eCpWUlCAxMRHz58/nj0mlUvTr1w/Hjx8X7TrFxcUoLi7mb9tbDok8Zrhi3prXz8hVgznqqyIDgJeLY7Wfs+xCTUKsxexhsWnTpmHVqlXljn/++ecWLeB29+5daLVaNG7cWHC8cePGyMrKEu06y5cvh0ql4n/8/PxEe25iW6w9lMKtVynrYqZ+UWVdzoEQYnZw+fHHH43u49KjRw/88MMPojTKmubPn4+8vDz+5+bNm9ZuErFjxgKHVCLB5hNpdToHQojZw2L37t2DSqUqd9zNzQ13794VpVHGNGzYEDKZDLdv3xYcv337tqjJeoVCAYVCIdrzEVKZyZGtsCDuguAYN2ussl4Vl1sJ9/dAYtp9UXIslK8hYjK759K6dWvs3bu33PFff/0VLVu2FKVRxjg6OiI8PByHDh3ij+l0Ohw6dAjdu3e32HUJEYux2WmGiygNebk4Cs6vqLZYfHJmuaGz6s6Co2E4Iiazey4zZ87ElClTkJ2djT599DNsDh06hE8++cTi60hmzpyJsWPHonPnzujatStWrlyJgoICjB8/HgAwZswY+Pr68ns1lJSU4O+//+b/npGRgaSkJDRo0ACtW7e2aFsJKauiWlzBviogW19brFTHwKAvZskVtOQ+7DNy1YjZpe/lGNYW43ouVV2nKtae4EDsi9nBZcKECSguLsbSpUvx/vvvA9BP8127di3GjBkjegMNjRw5EtnZ2Vi0aBGysrIQGhqKvXv38kn+9PR0wd4LmZmZCAt7vKNfbGwsYmNj0atXLyQkJFi0rYSUVdGH99U7D+At0dcWeya4KX4+lwnD9QE5BfrZixIAWqYvy++icKhw+Kq6QYJmmhExVWudCyc7OxtKpRINGjQQs002hda5EEsLmBePvo5XcaikDXzdlSgoLkWuWgMAUMql/BRlpVwGTxdH/v7qrM+hvAoBbLT8i1qtRmFhIQDA29sb9+7dw8qVK7F//37RG0dIfSA1SLpk5KpRUFIKAOjkq8K7g4L4nIzCQYqj8/pgdlRgtWeRUV6F1Bazg8vQoUOxadMmAPp9XLp27YpPPvkEQ4cOxdq1a0VvICF1TVUJ9bL3G44dKOVSlGr1By5m5iE6wh/vD+sIX3clv9eL4UZjhs9jSiKfpjeT2mJ2cDlz5gyeeuopAMAPP/yAJk2aIC0tDZs2bTK6uJKQ+qaq3kHZ+4N9H0/tL9Lo+NsdfPR/VrTYk3uemF0XMGT1n1gQd6HKXom1F46S+sPs4FJYWAhXV1cAwP79+zF8+HBIpVJEREQgLc12tjslxFqq6h1MjmwFd6WcL1KZnlPI3yeRPF6hz1VL5nok07adFfRMuI3DtEy4VTL1SogtqNY6l7i4ONy8eRP79u3DgAEDAOi3P6aENyFVi47wh4vCAblqDWL3pSDvUfIeAHRMHyxkksdBInZfCjJy1fj5nHBNy5Er2eWee0iID07eyEGr+fGYtu0sgIqDEyGWZHZwWbRoEWbPno2AgAB069aNX8C4f/9+wbRfQuoTw3wHFwxi96VUeF64vwfclXLkqTUwNl1zUCcfREf4Y8jqP/mZYw5SSYU9IqVcig+GdcSqUWGIT86ElgHxyfo9lypbcEmIpZi9zuWFF17Ak08+iVu3biEkJIQ/3rdvXzz33HOiNo6QusJYnqW4VIeeKw4Lpv1y5wHgey/GVujHJ2eiawtPwXCXi8JBMPV4dlSgYD3L2oTrOHkjB44OMhRptBjUyQeAfvgsK0+NJionZOUV8cNphFiS2T0XAGjSpAnCwsIECxa7du2Kdu3aidYwQmyVsVlZhnmWp9t6QyYBAFYu4Biex/3dSV7+v6GWATG7LgimKXOzxYzheku7z2VCrdHCx12JVaP0IwlHrmRDy4DM3CJoGZCYdr/G7wEhVanWfi6E1GfGyqsYrm5fm3AdWgYoHGTwdFEIhrHKroKPjvBHm3d+4f8nyiT6WWLJGXnQMv3CyZJSLT9MNm3bWcQnZ6KDjwrnM/LAgHLDb4b5GkNOjxZhUsKf1AYKLoSUUdUq9srKq2w+kYaC4lK4K+V8T8Ow52KsHphG9zjrMqiTD5+ol0C/cFKt0SIx7T42n0jD7nP6PIrhcBmg79Usjb/ED4cZtttw+IymIJPaQsGFkDKqKvxYWQ2utQnX+dIsgH5oS8uExSczc9V8j8NFIfwvWHbI6um23nwgMpaIlxuMm3m6OCIjV13uOahmGLGGauVcCLFnNVnFbvhYbniMG6YylmMpe42MXDU/NZlBn9j3cnF8VA2ZQSbRl4WRPYoppTrGT2kO9/eATAJK2BObUK3ClVevXsVvv/2GO3fuQKfTCe5btGiRaI2zBVS4klSX4fAa8Dg38nRbb37oq7mnM5Iz8vjClabwdVci3N+jXPVkd6XcaEFLKlZJyqqNzzWzh8W++eYbTJ48GQ0bNkSTJk0gkTzulkskErsLLoRUl+FwVM8Vh/n1KkeuZONBkQZaBuSWyZ0YI5XoF1fKpRI0cnPie0VlvxXmqTVwV8rL9Yaqu78LITVhdnD54IMPsHTpUsydO9cS7SHELk2ObIXYfSkoLtVWuHCyIjoGIxME9M8g4f+m/zNPrcHJGzmCngptAkasweycy/379zFixAhLtIWQOsfULYWjI/yRFDMAni4KkwKLXCZcWumicEB0hD/fC8nILQKgDyhDQnz4hZhcnsZwfY1hgDFsZ3W3Qxbr8cS+mR1cRowYQXu3EPII92Efuy/FpA9aLqmvNLJw0pBGKwxBXK+De7yvuxN/X3zy49yLBPp1MoaFMStqZ033drGFvWEowNkuk4bFDEvpt27dGgsXLsSJEycQHBwMuVwuOHfatGnitpAQCzI32V32fK5HUFBcalJeg8vDDFn9Z7m1KqbgHh+65PEXvA4+KtwrKOF3qLxXUAIXhYOgPcbaWdPhMlsYbqN8ku0yabZYixYtTHsyiQT/+9//atwoW0KzxexbzxWHkZGrrnTLYMOAwn2YySTAkqEd+Q80rmAloF+0WNkH3eYTaVgQd0FwzHC2mGEehVO2faFL9vMTBLj7uHaG+3vws9HKtsXeZo7Z2+upLbXxuVatqcj1CQUX+2bKh5NhAAr39+BXybsr5UiKGcCf137hXqg1WijlMlx6f2CF1+MWVhqqbCqyVKLfrdLH3QmZuUVwkkvRppErzmfkwUkuxbuDggRtNyVgkvqtNj7XaBElqddM2ZnRcHMvwz1U8tQawVh/kUYLAFBrtOVyAFxuIHZfCrQMRishG+OulPM9mYzcIjAAao0OyY/qinm6KAS9p9Al+5FTUGx0SnJNUX6DmMPs4PL888/jww8/LHf8o48+ollkxC4Zbu4FCD/wDZPZg0N8+L+vTbiOadvO8pt2ccNpxaU6yCT6HSeropRLMTsqEIM6+UAmgSCJz7mTX4Qhq/9Ei3nxWBB3AblqDdQaHVwUDuU2DaspW0jgk7rD7OBy5MgRPPvss+WOP/PMMzhy5IgojSLE1nCztGZHBSIpZgDeH9axXImYVaPC8IHBccNNu8L9PSCBvlejZfq1K1UpKdUhOsIfq0aF4fryQTg6ry9fs4yj0TG+F2Mo3N9DcH0xeh01KYtD6h+zcy5KpRJJSUkIDBTuLXH58mWEhYVBrVaL2kBro5wLqS6uPP6gTj5ITLvPbxImkwCuTg7IVZfyq+8ry7nIpRLEDOkAAHzlYy7/4iCVCKoqc7j8UNnrUx6GADaacwkODsb27dvLHf/+++8RFBQkSqMIsQdcj6NrC08UFJdCKZfBXSnHkqEd8aCoFIBpPRiNjiFm1wUsjLsAtUYLBiCnoAQqpbxcYJEAfL7F8Po5BSWQ4HFRS8qfEEszu/zLwoULMXz4cFy/fh19+ui/AR06dAjbtm3Djh07RG8gIbbOcAowVx7fcIKAYRl+rtdw8kYOP+vMFGVnlxVpdFBrhEVjlWVmjnHtKiguhfrRZAOuHD+tDyGWZnbPZfDgwYiLi8O1a9fwxhtvYNasWfjnn39w8OBBDBs2zAJNJMQ6yn67r+jbPvdBvftcptGEt+FsM+6xq0aFQSmXVbttwb4qwV4uwOMcTdl2AfrAY9hzofwJsbRqbRY2aNAgDBo0SOy2EGJTyn67N/Ztn9t5kps9ZmyL4egIf8TuS+H3XQH05fe53kRFjC2m5Fy98xA6g3SpTKLfxZJTdkdMru1cz4U2ECOWZnbPpWXLlrh3716547m5uWjZsqUojSLEFpT9dm/s2z4XNJzkMvi6KwWr9o0pLtXyU4arUlk6htvOWCbRF668vnwQVo0KK9eu4lItP2THtZ3yLaQ2mN1zSU1NhVZb/htXcXExMjIyRGkUIbag7Lf7yr7tKxyklc7Cmh0VyH/gVwc3q4wzOMQHq0aFCQKK4ew0TpFGxw+Nce3jVvBTvoVYksnBZffu3fzf9+3bB5VKxd/WarU4dOgQAgICRG0cIbaOG3LiejOGyf0Df99GkUbLL66sbmAB9IGlbIApO5GACyK7z2ViSIgPjlzJRnGpDgoHKd9j4c4Hyg/fmYrqeRFTmLzORSrVj6BJJBKUfYhcLkdAQAA++eQT/Otf/xK/lVZE61zqB7E+MLlegUzyeIZXZbkTjjnbHHPkFaxxAcAvtszIVcNdKYeLwoGvmlzTtS5Uu6zus6l1LjqdDjqdDs2bN8edO3f42zqdDsXFxUhJSbG7wELqD7FKm3Azw6QG9V0cpCbUejGBu1I40GAYWMqWhwn39+BzRAD4Xo0pM8SqysnQTDNiCrNzLjdu3LBEOwixqpruTWI45PSgSCNYl1JR78IUhr2eXHUplHIpSrUMGh0TDJN18FEJ9odJTLuPVaPCEB3hL+iVAY/roVXUQ6tqDQzNNCOmqNZU5IKCAvz+++9IT09HSUmJ4D7aLIzURVV9YBomyw2T6Jyl8Zeg1mj5HoJYfNyV8HJx5AMHt3Cy7JDY+TIbj4X7ewiCSri/B2J2XYCjgxRqja7SZL4tbAJG6j6za4udPXsWzz77LAoLC1FQUABPT0/cvXsXzs7OaNSoEW0WRuxSq/nxfG+EWzti+OHcYl58lXmVylSUcymbxDfGXSlHcakWao2O7+kYDof5uiuRlafmS/37PBrSKruSv2y+iRL39sumci6ct956C4MHD8b9+/ehVCpx4sQJpKWlITw8HLGxsZZoIyE1Ys66jorONZzem6vW8ENL3PnBvirIxEmtCFQVWAzzLEq5FE5yGZRyKQqKS+Hl4giZRN+L4dbEDA7xKbd/TUX5JiqxT2rC7OCSlJSEWbNmQSqVQiaTobi4GH5+fvjoo4/wzjvvWKKNhNSIOR+SFZ3LldN3V8oFG3Fx56fnFMLRofrlXKorI7eI38NF/6PvweSqNUjOyIOWAUeuZCMx7T6WDO1odEjPMEFvGFwpcU9qwuyci1wu56clN2rUCOnp6Wjfvj1UKhVu3rwpegMJqSlzcgiVncvlZbjhIkDfK8jKU6OguNRo4t5wSrIluCsdkKcurXBIjutNZeSqEbPrAoDyiXzDfJPhAsuqdugkpDJmB5ewsDCcOnUKbdq0Qa9evbBo0SLcvXsX3333HTp27GiJNhJiMmN5gqqS9WUfwwWQ0CX7AegXSgLg64Jx+Y3YfSlwUTg82vzL+Me7JQMLoJ9BZkgpl8HTxVFQoRkAYnZdgJaBT+RXlE+pKhBTHoaYyuxhsWXLlqFp06YAgKVLl8LDwwOTJ09GdnY2vv76a9EbSIg5qpMnMPYYrkw+l1/hSrdwQ1CAPvcS7u8Bd6VctLUsNeHrrkT/oMYAgK4tPPkFjmsTrmNQJx/BEFdF71N0hH+lPRbKwxBTmd1z6dy5M//3Ro0aYe/evaI2iJCaMGcIrLJyKJMjW/E9lcmRrbA0/m+jz3HkSjZcFA41Ku0iloxcNXIKivleFbeXS65ag6w8taCoZnWnG9M0ZWKqaq1zIaS2mDsMY84CP643UlBciqSYAeWuaTjdWF8iX1eulEueWgMXRe0n8ivC9aryi/S9LAn004+5ITEAgoDKMfV9pgWUxFRmD4vdu3cPb775JoKCgtCwYUN4enoKfggRkyWGYbgZUcWlxvdT4a65IO4Cpm07C0Cfd/F1V2JwiA/clXJ+oy8G/YwtW6Njj1f3c0N2Xi6O/GuLT9ZvbBa7LwU9VxxG7L4U/n02Nh2byvQTc5ndc3n55Zdx7do1TJw4EY0bN4ZEYv2xZmK/LDEMw33Auivl8HRRlHvuyZGtsCBOP7MqPjmTL6PCPZZL8C+MuwCGygtIWpOTXApPFwUyH1UNuJiZhyVDOwoqKRcUl/LvBZeT4d4fbmjN8BiV6SemMju4/PHHH/jzzz8REhJiifYQImCJYRjDgFVR7az/nrqJ5Iw8dPB5vLUE9wEbs+sCXJ3k/PCYrQUWw2CnzxfpS9Nw63AMKxlvPpEmyC0ZBlEu8BgGGGvmWmimWt1i9rBYu3btoFaLWz+J1L76PMxR1YwoALhXUCL4E9APKwH6/EVBSanRx1mbUi6Fi0L/nVHhIMPahOv8dspqjRYxu/RDfYb/9g8e5WcWxl3gp18fndeHHwrkPsyres8s/TtFM9XqFrODyxdffIF3330Xv//+O+7du4f8/HzBD6kb7OE/qikfZhWdY0pZeXelHDkFxQhdsh+bT6ThYubj4pClWoYhIT5GH1vbDAemudX5cn5qtLBXpWXg8y0xuy4gdl8KvxaHQVjaxpSAYsjSv1NUMaBuMTu4uLu7Iz8/H3369EGjRo3g4eEBDw8PuLu7w8PDo+onIDbBHv6jmvJhZuyczSfSsDDuAjJy1VgYV/6bPKD/YHVROPAf1gvjLgj2aAn2VeHIlWzLvDAzGRuU0+gYctUawWQDX3cn+Lor+TpjWgZ+RlknX1W50jbmsvTvlLnBjliX2TmX0aNHQy6XY+vWrZTQr8PsYUqpKXkAY+esTbjOfyAzAD+fywSDPkEPQLAWhJuuzKD/wOZ2X+y54rBNrG0xR0ZuEdyVcnRt4YmuLTz5VfsMQHpOodHp2ObkN+zhd4qIx+yS+87Ozjh79iwCAwMt1SabQiX37Q+XxOb2l+fKuQD6Iab3h+nLGHEfrlyAAfS1vB4UlaKJygmZuUU1KrNvqDrbHJuCG7rbfS6TP+aulCMpZgDfg2MGxzi0lbF9s8mS+507d6YClaTO4r6RP93WG54ujpgdFYh3BwXx9zPog4rhcBo39RjQ1/LSMn0vwLbmiJXn667EqlFhWDUqDO5Kebn7oyP88f6wjvB1VwpeI6AvyMmV6yekOswOLlOnTsX06dOxceNGJCYmIjk5WfBjaWvWrEFAQACcnJzQrVs3nDx5stLzd+zYgXbt2sHJyQnBwcH45ZdfLN5GYrvKLiLk1m0YfvjmFBTjVp5+RuSd/CKcvJFjrebWSGauGtO2ncW0bWeR9yjJr5Tr/8tz+aWK8hiJaff55L81ZhTW59mM9sLs4DJy5EhcunQJEyZMQJcuXRAaGoqwsDD+T0vavn07Zs6ciZiYGJw5cwYhISGIiorCnTt3jJ5/7NgxjBo1ChMnTsTZs2cxbNgwDBs2DBcuXLBoO4nt4pLOhoUcuQ8wpVwGCfQzrrilKxodEwwp1SUM+uEwLqekYwyeLgrBjLCKTI5sxSf9rTGj0B5mM9Z3Zudc0tIq/ybh72+5hF63bt3QpUsXfP755wAAnU4HPz8/TJ06FfPmzSt3/siRI1FQUIA9e/bwxyIiIhAaGoovv/zSpGtSzkVctbUQzth1Kro2l1+w9N4rlREz5yKT4FEw0d+WSvS5pA4+KqTnFAJAuW2azXm/agMtmLSs2vhcM3u2mCWDR2VKSkqQmJiI+fPn88ekUin69euH48ePG33M8ePHMXPmTMGxqKgoxMXFVXid4uJiFBcX87dp7Y64aquMiLHrVHRtLmlfXKqFwkFW4cZfdYWW6RdTcpMUdEyff0nPKUSuWgNfd6VgT5dwfw/EJ2cK9nsBrDv7i2ae1X11piry3bt3odVq0bhxY8Hxxo0b4/Lly0Yfk5WVZfT8rKysCq+zfPlyLFmypNzxHTt2wNnZuRotJ4aGuamRqilAgJsLtm3bVqvXqejaMgARyIZGogO0gEwugda8Dn2NNZQWoK/jVVGf081NjvyiR9OlCwG5TAqZi4R//cev3UU7jRaFl4BIub53E+jmVuG/yz/31Ui9V4AALxc081CK2lZSuwoLCy1+jToTXGrL/PnzBb2d/Px8+Pn5YcSIETQsZmNMGTrZfCINnz86p2kLIC7hOro/0Qqjypz/3sJf+W/61ihEWd1hMQkAH3cno5WZZRrA1UnOL5RUKeWY3V8/K+zzhOvwcvfDxUx9/bR7BSVVDkH1XHEYGQVq+MqVOPoGTU+uy/Lz8/HKK69Y9Bp1Jrg0bNgQMpkMt2/fFhy/ffs2mjRpYvQxTZo0Met8AFAoFFAoFDVvsA2zl/FsU4bYDCv8PijSQMv0W/6evJHDbwMcHeEPhYOMDy6ljwKLTAI4OjweXgJsrwJyZSX/uRX43AyxXLUGS+P/RpFG9+hx+hlx9wpKTFrLYgvFK0ndYfZsMWtxdHREeHg4Dh06xB/T6XQ4dOgQunfvbvQx3bt3F5wPAAcOHKjw/PrCXmbimFJuhDsHeJysN6yvxb0HT7f1hgT6XMXgEP1MsiVDO+LdQUGCacq2FFiM4aYbyw22XS4p1UHxqCKy+lFg4cgkMPr+GZsKTOVXiDmqHVxKSkrwzz//ID09XfBjSTNnzsQ333yDb7/9FpcuXcLkyZNRUFCA8ePHAwDGjBkjSPhPnz4de/fuxSeffILLly9j8eLFOH36NKZMmWLRdto6e6grBpj2Ycedw1X4HfIocJTdUz4x7T4Y9JWEj1zJxp38IiyM0/dwuCrDdUHMkA4o0ugEQZDbNqDsQkp3pRxLhuqrERgGks0n0hCz64JdfAEh1mP2/5qrV69iwoQJOHbsmOA4YwwSiQRarfHd/cQwcuRIZGdnY9GiRcjKykJoaCj27t3LJ+3T09MhlT6Olz169MDWrVuxYMECvPPOO2jTpg3i4uLQsWNHi7WxLqiPM3HKvubNJ9KQmHYfJ2/k8DPF5DJJuXphu89l8sNKdQG3yRmgz8cMDvHhZ4IZvo5Ovip+OwGuvE3svhRER/hjbcJ1aFnFvRrymL0MMVuC2etcevbsCQcHB8ybNw9NmzYtV7jS3jYRq+11LvTLWnOmvIehS/bbVOFJMde58Mn7qEB+CFQm0U9JZo/ulz5a0+OulCPvUWFO5aOdK7ldKul3sGp1tQabTa5zSUpKQmJiItq1a2eJ9tR7tJ1szXBDOmXXbBjevzbhOopLy/ewpRKAsUd7zksAzaMkDbcXfV2hUjrgQZEGJ2/kINzfA1l5anTwUeHqnYco0mjhJJdBrdFC9uh7IYO+l6JwkPFJflM/KOv7lyGa5FAxs/v7QUFBuHv3riXaQmA/+RBr4YZ0JAAKikvL1aaK3ZeCjFw11BodlHIZ3JVyfrjIzUkOH3clNDoGnUHOgvu2X1dwxTV/PpfJ1wi7mJkHtUYLH3cl3h3Uns87AY9zL4Y7T5rKXiaHVBdNcqiY2T2XDz/8EG+//TaWLVuG4OBgyOXCJCGtBamZ+pgPqa5p284iPjkTgzr5oGsLT361OaAPLLlqDWJ2CfdoMeTp4oij8/oI9pHnHs8NDXHf5OtSz4UrY8Og35qZ67kYrmWJjvDn96ThVuwD5d+nqnom9M2dVMTsnAuXMC+ba6mNhL41UG0x29VqfjyfeG6iUgrGvg2HxwzHww0DyeyoQJy8kYP45Ew4OuiHisqeu2T3xVqZfixmzsVwLQ43pOeulMNF4cAHAS4QV5VbqSynUN+HxOoym8y5/Pbbb5ZoByFmG9TJp1zPhfvwNKwnZvitumzPkAtARY8CS9kdK219XYsxmjJDekq5fkM0bkaYi8LB5NxKZT0Tyg+Sypjdc6lvqOdin7hv3QBDRm4ROvmqsHvqk/x9sftSarWApZg9F6VcCoWDTDAbTvkoie9uMIvMWI/DnN4I9VzqLpvsuQBAbm4u1q1bh0uXLgEAOnTogAkTJkClUonaOEIsxXCKLgB+zQd3X9lpyrZW9sUYuVSCq8ue5T/0XRQyg9IwjO+ZGfbqAGGexZTeiGFQqUvTb0ntMnu22OnTp9GqVSt89tlnyMnJQU5ODj799FO0atUKZ86csUQbCRFd2U3Dwv09ELpkP9ov/BU5BcWC8imAfqjJ1meMcTXRuBlxhjXHijQ6PrBUtgLflNmK9X2GGDGN2cHlrbfewpAhQ5CamoqdO3di586duHHjBv71r39hxowZFmgiIeYxZ4vcri08cXReHySm3UeuWgO1RvdoJ0r9B7VSLoUEgFwmgZNcCqkNR5hgXxV6rjiM4tLHhTYleJzUj9l1ge91cBMhwv09BO+VKVNrabo8MUW1ei5z586Fg8PjETUHBwe8/fbbOH36tKiNI6Q6TPlmXfacyZGt+DUv7ko536NROMjAoF9Qabj9sS1KzylERq4aRRotlHIplHIpVEo5Bof4CLYs5oLDkqEd+enW5vRCaG0HMYXZORc3Nzekp6eXW6F/8+ZNuLq6itYwQqrLlLUXZc+paH3RtG1nsftcpsXaKiau6gBXgDPvUU/syJVsLBnaUZB8L5tnoV4IEZvZs8WmTZuGn376CbGxsejRowcA4OjRo5gzZw6ef/55rFy50hLttBqaLVZ3GS6yXDUqzKTHlJ0Bxa3zMKSUy9A/qDGOXMkWbUZZTWeLGc4Q40bumMF9VDOMGLLJ2WKxsbGQSCQYM2YMSktLAQByuRyTJ0/GihUrRG8gsV+WnsrKVQOOT86sNLgYtoMbLuNW9nu5OCIjVw25TMLXGisp1eLA37eh1tT+gmHpowKUZSkcZHi6rbdgQShXwLK4VIuMXDUyc9VgMF5zjRCxmZ1zcXR0xH/+8x/cv38fSUlJSEpKQk5ODj777DO738GRiMvSs44GddLnGgZ18qk0yW/YjsmRrQT5iYuZeQAeF7EE9PdZI7AAxgMLh6sjpnCQwtddifeH6euFFT3aSdNJLqVEPKk11d4FydnZGcHBwWK2hdQzlq5LtWpUGN9j4Ya3jH1rL9sOVyc5f/zkjRz8fC5TUFvMXSm3iXL93NobCfQ7aR65ks2vxgfAl7bhqh6/OyiIeiyk1piUcxk+fDg2btwINzc3DB8+vNJzd+7cKVrjbAHlXOyDqUNwxmpplc27KOVSlGqZqIsqa5Jz4eqGZeSqBdsDcD0wmQRYMrQjBRbCs5mci0ql4gtVurm5lStaSYitq6raNBd8uKrIhr2pyZGt+N0aAf2CRFuakfx0W2++tlpOQQmfbxnUyYcS+MRqqLZYFajnUj9wvROuF8DNrDKcYQXo8zBeLo44n5EnaoAxpecil0rg8KhejWHPqWwlZ6r3RapSG59rZif0+/Tpg9zc3HLH8/Pz0acP1Rki1WfOynoxrmH493B/D8gkQHGpDhm5asQnZwr+5HI1R+f1wb2CEqv0XDQ6/UJOhYMMjdycAJTf576mCxxr49+A1A9mB5eEhASUlJSUO15UVIQ//vhDlEaR+qk2alYZXsPw79xMq6JHlYMHdfKBUi6DlunXtZQdJpPV0siwUi4rdyxXrcGd/CIo5VJ+8oFYqG4YEYvJwSU5ORnJyckAgL///pu/nZycjLNnz2LdunXw9fW1WEOJ/auNmlWG1+B6K+H+HnzAYABcFA7o2sKTn25cdtpxdIQ/lgztyG+PbEiM2mPcU3TyVeHdQe3LFdEE9L2YklIdv9sm19Mo2/OYtu0sWs2Px7RtZ026NtUNI2IxOecilUr5RL6xhyiVSqxevRoTJkwQt4VWRjkX+1VRnqVsAh+AYIdLLqfBbTTGKbsy3hzGci7cNVvMiy/3nFyVAG6hKHcu95q4GWJcG2US4PryQdVomX2gXJSQzcwWA4AbN26AMYaWLVvi5MmT8Pb25u9zdHREo0aNIJOV78ITYqu49S0FxaWCnRk3n0hD3qPAwpVOmRzZSrB18sK4C+AmTUolgEyqX8Ev1miZXCpBTkEx2i/8tdwam9lRgfwHZNkdOA2D3tqE64LdOusz2jWz9pkcXPz99f8gOp2uijMJsQ7D6cRVTcE1/CYLoNzfjS087LniMN9TYQC4DrybkxwPijT8cTHoGINa8/jZJADeH9aR348ldMl+FJdq+bIvZTf+Mnw9iWn30bWFZ7nXXZ8+ZC29YJeUV+2pyH///TfS09PLJfeHDBkiSsNsRV0bFquvHx4ABENChkNFHGM1xMqeU/Y8w/eQK4QplQh3pTSs91XdHSu5YbFOviqk5xSiuFQLtebxFzl3pRxJMQMEr7MsY68ldMl+5Ko1/OONLRIl9Y9NTkX+3//+h5CQEHTs2BGDBg3CsGHDMGzYMDz33HN47rnnLNFGYob6PNun7O6SZb+llq0hVvYcLhkOwOh03iNXsqFlgINMIpgtZhhLri57Fu7K6s/gSs8pxIMifal87hISALOjAgWvs+w1TE3CU8Ke1Bazg8v06dPRokUL3LlzB87Ozrh48SKOHDmCzp07IyEhwQJNJOaozx8e3BqPVaPCjAaHqt4bUwOzwkH2aLZY+RxjzxWH8XRbb/i6KyGvxnzlPLWGT8APDvHhC1Byr4XrVc2OCkQnXxUA/awyw1I1hmtUZkcFwtddyQcn2uiL1Bazh8UaNmyIw4cPo1OnTlCpVDh58iQCAwNx+PBhzJo1C2fPmjblsa6oa8NixDTGhocqysNweY7YfSkA9OVWEtPuo6C4tMICluYOjxnOFqusFljZGW6GQ3c05EVMZZPDYlqtlt9xsmHDhsjM1O/S5+/vj5SUFHFbR4iFGOvFcN/qAf1+84a9GMOpyUeuZCMjV43iUm2Fs8OqW9TSXSnnAws3TDdt21m+R8K1m9ujhQt4Fb0mQqzF7JL7HTt2xLlz59CiRQt069YNH330ERwdHfH111+jZcuWlmgjIaKrqJCl4XTjsqVVAP3qeG5Ro8JBhncHBWFB3AXR2wY8HqbLylPzU4u5Ia3QJfsFCf/KXhMh1mB2z2XBggX8dOT33nsPN27cwFNPPYVffvkFq1atEr2BxH7UhbpVaxOuQ8v0SXTD0iqGCfVSHRPkMSylsgkKZXMphNgaUaoi5+TkwMPDwy5L8VPORTx1ISfA5V24fIphW7mpyIM6+WDVqDBBL6emuJzLB8PK51rMWb9DiClsMudi6ObNm7h58yY8PT3tMrAQcdlSTqCiXhSXd+F6BoZtXTUqDNeXD6owsHAlwCQAP5PLXMaCBjc8ZlihuabqQi+S1G1mB5fS0lIsXLgQKpUKAQEBCAgIgEqlwoIFC6DRWH/rV2K7xJ4Ga+4HpGGCnEvYx+5LqTTIVJWXMaR7VEF5cIgPzmfkVft1lVXV+p3qqM/roUjtMDuhP3XqVOzcuRMfffQRunfvDgA4fvw4Fi9ejHv37mHt2rWiN5IQY4zVi6qsQkHZBDm3DMXwOUypcMDlZWQS/W6PR65k8zPJ1BotjlzJFnW/F0sk6qkcCrE0s4PL1q1b8f333+OZZ57hj3Xq1Al+fn4YNWoUBRdSa4x9QFZWoJA739jukoZ/r6rAoeF1uXOmbTuL3ef00/KLS7VGH2cqwzU1hkUqxUQzy4ilmZ3Qb9SoEX7//Xe0b99ecPzSpUt4+umnkZ2dLWoDrY0S+nVLTWurVVXQ0tjzlk24c/vYm4NL6Pu6KwWLM02d/FCfa8oR89XG55rZweW9997D5cuXsWHDBigUCgBAcXExJk6ciDZt2iAmJsYiDbUWCi71l+HsNgCCmW7GimByq+ZzCorLrUGpSl/HqziibQuNlgnqhpnac6kLM/GI7bCZ/VyGDx8uuH3w4EE0a9YMISEhAIBz586hpKQEffv2Fb+FhFhJ2WE3Y8NnsftS+JX6xaU65KrVUMqlkEnAV0+WwLRS/I1cFQAk1ep9UA6F2BqTgotKJZxW+fzzzwtu+/n5idciO0HDFHWfYV6i7Gwyw43GHvdSGD+spdbo4Oqk/+9VXKpFkUZXZYDJyC1C6orHu0Wa8ztEORRia0wKLhs2bACg39745s2b8Pb2hlKptGjD6jra+c6+lP335H42n0jDwrgLgsDBFbYM9/fgtyGWQLjvS3WuSUhdYtY6F8YYWrdujX/++cdS7bEbtrRgsL4Sc6FgRf+e0RH+eH9Yx0d5GQkyctU4ciUbR+f14fd/AfTDYubWsqTfIVKXmZ3Q79ChA9atW4eIiAhLtcmmUEK/7iq7C2NtXY/bkpirpKyUS1FSquN7MA5GyvFzs8UMh8UIsRSbLP+yYsUKzJkzBxcuiFsJlpC6bnZUIGQSfS8lZtcFPN3WG+5KORQOMgzq5AN3pRwMj8vxly2YVI29xQixWWYvohwzZgwKCwsREhICR0fHcrmXnJwc0RpHSE3MjgoUZQaVqYl17j6uNMyRK9l4UKTfWTI+OROODsLvciqlHMWlWn5CgGEVZkLqOrODy8qVKy3QDELEJ9YMKlMT61wQGtTJh9+pksu5aBkEa1+4PWGKDI4193SucVsJsRVmB5exY8daoh2E8AICAuDk5IQLFy7AwUH/K9q5c2fExsYiMjISgP738KeffsKtW7fg4uJi1mP37NmDJUuWID8/HxqNBgMGDMBHH33Ejz0zxrBmzRp89dVX0Gg0KGIOyNU54aUZsyptNxeEAPALLQ0LXMqlEjRyc8LkyFZYsvtiuS2Sk8sUu6Tp7KQuMzvnkp6eXukPIWIoLi7GunXrjN6Xn5+Pn3/+GSEhIdixY4dZj927dy9ee+01fP3110hJScGVK1cgl8vxr3/9C9zcloULF2Lr1q349ddfcfnyZaSmXMCBrWvRXFZ5peOys7uiI/wxqJMPf7+DTMKvnjdlG+TYfSnltjImpK4wu+cSEBBQ6d4tWm3NivYRAgCLFy/Gu+++i5dffhnOzsLhom3btqFfv34YNWoUPv30U4wbN87kx37wwQd49913ERYWBgBwcHDAJ598gpYtW+KdNd/jt3tuOPPxx0g+dw7NmjXjH9elSxd06dKl0jYbG4ZLTLsvuN1+4a8VlobxdXeq9PkJqUvM7rmcPXsWZ86c4X/++usvfPnll2jbtq3Rb5GEVEdISAh69+6Nzz77rNx969atw4QJE/Cvf/0LV69eRUpKismPPXPmDL9VBMfR0RHh4eHY9ksC0q6nQCeRo127dqK8jsmRreCulPOzxgwDS9mvaFl5RYLbpmxlTJt+EVtlds+FqydmqHPnzvDx8cHHH39crg4ZIdX1/vvvo2vXrnj99df5Y+fPn8etW7cwYMAASKVSREdHY/369fjwww+rfGxVurf0wkXmhHz54+9carUa3bt3R0lJCXx8fHDw4MFql2UxLMuvlMv4yslckDEcQiv72IrQKn5iq2q0zbGhwMBAnDp1SqynIwQBAQF46aWX8MEHH/DH1q1bhwcPHqBly5YICAjAtm3bsGnTJpSWllb5WAB44okncPz4ccGxkpISJCYm4pXn+uLYh2Og1ZTwvaEfz92By4uf4rkpi3D37l0Apu/iOG3bWbSaH49p284CeDxEpu/FPP6vJ5UA/do3xqpRYea8PQBoFT+xXWYHl/z8fMFPXl4eLl++jAULFqBNmzaWaCMA/fqZ0aNHw83NDe7u7pg4cSIePnxY6WO+/vprREZGws3NDRKJBLm5uRZrH9ELCAhAYGAgQkNDERQUhDVr1iA1NRUymQyhoaEICQlBSEgI4uPj+cds3LgRKpUKYWFhaN++PTIzM/Hll19CrVZjwYIF2Lx5M5KSkjBs2DBs3rwZJ06cQGpqKjZs2IDMzEwwxhAfH4/U1FRkZWXhySefRGhoKP/YzMxM/lrz58/HW2+9BVdXVwD6bbtnzZqFgIAA9OnTBw0aNMDMmTPxyiuvICMjgw8kO0/+j38ObqiroLi00uEorq5YfHIm/ziufL/hTDFHBxn+ua+u1vst9tbRhIjF7ODi7u4ODw8P/sfT0xNBQUE4fvy4RXehHD16NC5evIgDBw5gz549OHLkCF599dVKH1NYWIiBAwfinXfesVi7SHnbt29HUlISfv31V7zzzjvIz8+Hq6srkpKScO7cOSxduhSjRo0STP7o3bs3Zq2Ng+fYNXD1aIhLly5h5MiRaNiwIaZNmwatVgtPT0+oVCo+H7Ju3Tp07twZ7dq1w7p16+Dm5gYPDw8sX74cAPjH3rp1i79OSkoKevXqhcLCQgQGBqJt27YoLi5GfHw8P1Fl6dKleOGFFzBw4EDcWDsJOd/Pg/zyPv55oyP84aJwQK5aU2nvZVAnH8gkQAcfFXquOAxAP0WZW8nPUWu0uJ5d+RclQuoas2uL/f7774LbUqkU3t7eaN26Nb+uQGyXLl1CUFAQTp06hc6dOwPQTyl99tln8c8//8DHx6fSxyckJKB37964f/8+3N3dzbo21RYzT0BAAOLi4hAaGgoA6Nq1K/7973/jgw8+4HuOarUazs7OuHPnDry9vbFx40bExcUhO2Iav+HVntfC4Ovri1OnTqFDhw4ICAjAnDlz8MsvvyA+Ph55eXkIDw/HqFGj8ODBA8Hi3oSEBMyYMQNJSUmCtl28eBGTJ0/Ghg0bEB4eXqOerOHuk0eu6HdfrWhjL2MbeZWtptzX8Sp6P/sc9UBIrbDJ2mISiQQ9e/ZEr1690KtXLzz11FP8N8kjR46I3kAAOH78ONzd3fnAAgD9+vWDVCrFX3/9Jeq1iouLyw39keo5f/48Ll++jKFDhwqO//DDD+jTpw+8vb0Fxw3zBx4eHmjTpg0uXrzI39+zZ0+kpqYiMzMT27Ztw4gRIyCTyUxqi0ajwaRJk/DVV1+Z/JjKcMNRiWn3kavW8L0YY7O3wv09IJPo/zR8vEopLPdSVQ6HkLrE7K5G7969cevWLTRq1EhwPC8vD71797bIOpesrKxy13NwcICnpyeysrJEvdby5cuxZMkSUZ+zvhk5ciSUSiWcnZ2xfv16yOVyPHjwAKGhocjJycHdu3dx+PDhco8rOzvKWKf65Zdf5ns6W7ZswZYtW0xq05IlSzB8+HC0b98eqamp1X5tHMOeS0GxfjJBuL8HvyLfcPZWYtp9aFn5NS9c7bNwfw9or/2PkvLErpjdc2GMGV1Eee/ePUEZDlPMmzcPEomk0p/Lly+b28QamT9/PvLy8vifmzdv1ur17QGXczl27BheeOEFAOBzLmlpaZg3bx5efPFFFBUVVfgc9+/fx7Vr19CxY0fB8TFjxmDVqlVwcnIyawLJ77//jtWrVyMgIABPPvkk8vPzERAQgOzs7Gq9Ri7Rn5h2H0kxA5AUM4APIhJAkOw3NqPLcDrzqlFheLJ1QxoSI3bF5J4Lt35FIpFg3LhxUCgU/H1arRbJycno0aOHWRefNWtWudXVZbVs2RJNmjTBnTt3BMdLS0uRk5ODJk2amHXNqigUCsFrI+KSSCRYuHAhdu/ejbVr1+Ktt94qd052djZeffVV9OvXD0FBQYL7fHx8sHz5crMXOf7xxx/831NTUxEaGlqjHoyxPesNtz7mhskMd600VHY68/Frd6E9kUYBhtgNk4OLSqUCoO+5uLq6CkrtOzo6IiIiApMmTTLr4t7e3uXG3Y3p3r07cnNzkZiYiPDwcADA4cOHodPp0K1bN7OuSSyD+ybODRFVRiKR4JNPPsHIkSPx2muvAQB+++03hIWFQa1WQ6FQ4LnnnsPcuXONPn78+PFGjxcWFvKzv/Ly8tCsWTO8/PLL/CwvMRkLGIZbH3OBp6IFl1wgCvf3wMK4C+jjqEXsvhQKLsRumDVbjDGGCRMmYPXq1WjQoIEl21XOM888g9u3b+PLL7+ERqPB+PHj0blzZ2zduhUAkJGRgb59+2LTpk3o2rUrAH2uJisrC6dPn8akSZNw5MgRuLq6onnz5vD09DTpujRbzDTGZkSRqt8X7v6+jldxjAXi0vvPWKGVpL6xudlijDFs2bJFsG6gtmzZsgXt2rVD37598eyzz+LJJ5/E119/zd+v0WiQkpKCwsJC/tiXX36JsLAwvkf19NNPIywsDLt376719ts7WiluXFXvi+EMMoVDzWexEWIrzF7n0qFDB6xbtw4RERGWapNNoZ4LEZvhUBmXe+nreBUu7XtVqwQMIeayuZ4LAKxYsQJz5szBhQsXLNEeQuxGRRWLuYASuy8FBcWlfOHKslOVCanLzA4uY8aMwcmTJxESEgKlUglPT0/BDyF1kSVK1xsrcLn5RBoKikvh/mgBJVdjTC6T0pAisStmL6I0LLNBiL2wROn6stOVDUu+uCvlmB0VyC+6LNUa30CMkLrK7OAyduxYS7SDEKsytm6lpspOV16bcB2szP0A+IBDU5GJPanWfi7Xr1/HggULMGrUKH5x46+//iqoA0VIXVIbpevD/T0ggX6jMG53SWM1xgixB2YHl99//x3BwcH466+/sHPnTn5PlXPnziEmJkb0BhJiLxLT7oMB8HRx5Bdb9lxxGM09nSEB8HTbqhcUE1JXmB1c5s2bhw8++AAHDhyAo6Mjf7xPnz44ceKEqI0jxJ6UXfPC5XkuZuaBgWaLEftids7l/Pnz/Kp4Q40aNeK3gSWElFc2BxPu74GsPDU6+KjglCuj2WLErlRrJ0pjK/TPnj0LX19fURpFSH3AVVFOzyms+mRC6hizg8uLL76IuXPnIisrCxKJBDqdDkePHsXs2bMxZswYS7SRELvEDZMBQJFGS5uFEbtidnBZtmwZ2rVrBz8/Pzx8+BBBQUF4+umn0aNHDyxYsMASbSTELnEz1GZHBcJJTsNixL6YnXNxdHTEN998g0WLFuH8+fN4+PAhwsLCzNq4iRDy2MkbOSjWaHHyRg6tcyF2w+TgotPp8PHHH2P37t0oKSlB3759ERMTI9jXhRBivvjkTETK9X9S4UpiL0weFlu6dCneeecdNGjQAL6+vvjPf/6DN99805JtI6ReGNTJB5JHfxJiL0wOLps2bcIXX3yBffv2IS4uDj///DO2bNkCnY5qIhFSE11beEIhl6FrCyr8SuyHycElPT0dzz77LH+7X79+kEgkyMzMtEjDCKkv1iZcp9lixO6YHFxKS0vh5OQkOCaXy6HRaERvFCH1CVdzzHBXSkLqOpODC2MM48aNw/Dhw/mfoqIivP7664JjhFRXQEAAAgMDERISgtatW2Po0KE4duwYAGDjxo2QSCT47rvv+PP37NmDyMhI/rZEIkFoaKjgOTds2ACJRIKVK1fiyy+/RGhoKEJDQ+Hp6QlfX1/+9m+//YZ3330XwcHB/LHvv/+ef55x48bx5wcHB+Ppp5/G5cuXq3xNpuwTw9Uco/IvxJ6YHFzGjh2LRo0aQaVS8T/R0dHw8fERHCOkJrZv345z587h2rVrGDt2LJ599ln89ddfAAB/f38sWrQIJSUlFT7ewcEBiYmJ/O3169ejc+fOAIDXX38dSUlJSEpKwpAhQzBnzhz+du/evTFnzhycP38eSUlJiI+Px6uvviooacSdf/78eTz77LNYuHBhla/H2IZhZU2ObEXrXIjdMXkq8oYNGyzZDkLKGT58OE6ePInY2FgMGjQIoaGhkMlkWLNmDd566y2jjxk/fjzWr1+P8PBwXLlyBRqNBh06dDDpeu7u7vzfHz58CMaY0QkrjDHk5+fDw6PqYSxT9omJjvCH7EZDjKI1LsSOVGs/F0JqS7du3QT7BC1btgyzZ89GmzZtMG3aNJw6dQpr1qzh79+/fz++/vprtG7dGv369UOvXr34+xYvXowZM2YInn/jxo0YNmwYACAhIQGOjo5wcnJC+/bt4eXlBVdXVwBAVlYW3nrrLSiVSjg6OiI2NhZdunQBAEybNg0BAQGQSCRISkoSPH9t7BNDiC2i4EJsGmNMcDswMBDOzs7o1asXVq1aheDgYLzzzjs4d+4cACAsLAwTJ05ETEwMSktLsW7dOty/b3ouIygoCEVFRTh79ixcXV1RWPi4qKSPjw/UajU0Gg2++uorfPXVVwCAF154AX/++Sf8/asXQDafSMOf1+5WmpchpK6h4EJs2qlTp9CxY0fBMXd3d/z444+4desWnJycEBgYiAMHDgAAZsyYgfHjx2PmzJno27cvJkyYgAsXLph93ZCQEPj6+iIhIcHo/SNHjkRiYiKys7Px9NNPo1mzZmZfg0NTkYk9Mru2GCG1ZdeuXVi7di327duHS5cu8cdlMhmee+45LFu2DF5eXrhy5QoOHTrE39+tWzcsWLAA/fv3R0pKCr7++muTrpeamorr16/jiSeeQGlpKW7cuIHPPvuMvz87O5ufjZafnw8vLy94eXnV+HVOjmyF4wcuUkKf2BUKLsTqNp9Iw9qE6ygoLsXIkSPh5OSEgoICBAUF4ZdffkG3bt0EwQUA/vjjD9y8eRM5OTlYv349Dh8+LLh/+vTpAICUlBT+mEQiMXp97vi2bdvg6+vLJ/E9PT1x/vx5tG/fHsDjITrGGBo2bIgNGzZAKqXOPyHGUHAhVsdN1207fROOzutj9Jxx48Zh3Lhx/O0dO3YI1rTk5uaiS5cughlfAHD8+HEMHjwYM2bMwBdffIFr164B0CfyASA2NhaNGjUCAOzbt0/w2OXLl+OPP/7Av//9b8ybNw9ZWVnlEvZiWJtwHe0eDYtR4p/YC/raRayu7N7yFeEWJBYUl5a7b9SoUbh37x4+/PBD/tjhw4exfv16zJkzBwDQp08fHDx4EOnp6QD0Q1tbtmzBgAEDAAC3bt3iey0PHjzAnj17EBZm+SrFtM6F2CPquRCrK7u3fEW4Hs6DovLBxcXFBQkJCZg1axZatGgBBwcHNG3aFLt370anTp0AAO3atcPq1asxfPhwlJaWQqfTITo6Gs8//zwA4Mcff8TatWvh4OCA0tJSjBgxAuPHj6+yXa+99hri4+ORlZWFqKgouLq68j0kU18/rXMh9kbCys71JAL5+flQqVTIy8uDm5ubtZtTr3G5mcmRrexi+Mjw9chuHMOoUaOs3SRST9TG5xr1XEidYWoPp64wLA0zxX5eFiEAKOdCiNWYmmsipC6ingshVmLYE9t245iVW0OIuKjnQgghRHQUXAghhIiOggshhBDRUXAhhBAiOgouhBBCREfBhRBCiOgouBBCCBEdBRdCCCGio+BCCCFEdBRcCCGEiI6CCyGEENFRcCGEECI6Ci6EEEJER8GFEEKI6Ci4EEIIER0FF0Js0OYTaei54jA2n0izdlMIqRYKLoRYSWUBxHALZELqIgouhFhJZQGEtkAmdR1tc0yIlUyObIW1Cdf1AeTGbcF9hlsgE1IXUXAhxEoMA8i2G8es3BpCxFVnhsVycnIwevRouLm5wd3dHRMnTsTDhw8rPX/q1KkIDAyEUqlE8+bNMW3aNOTl5dViqwkhpH6qM8Fl9OjRuHjxIg4cOIA9e/bgyJEjePXVVys8PzMzE5mZmYiNjcWFCxewceNG7N27FxMnTqzFVhNStc0n0vDntbs0M4zYFQljjFm7EVW5dOkSgoKCcOrUKXTu3BkAsHfvXjz77LP4559/4OPjY9Lz7NixA9HR0SgoKICDg2kjgvn5+VCpVMjLy4Obm1u1XwMhFem54jDaFSbjsnMnHJ3Xx9rNIfVAbXyu1Ymey/Hjx+Hu7s4HFgDo168fpFIp/vrrL5Ofh3sjKwssxcXFyM/PF/wQYgncVORwfw84yWU0M4zYlToRXLKystCoUSPBMQcHB3h6eiIrK8uk57h79y7ef//9SofSAGD58uVQqVT8j5+fX7XbTUhluKnIiWn38WTrhjQ7jNgVqwaXefPmQSKRVPpz+fLlGl8nPz8fgwYNQlBQEBYvXlzpufPnz0deXh7/c/PmzRpfnxBjaC0LsWdWnYo8a9YsjBs3rtJzWrZsiSZNmuDOnTuC46WlpcjJyUGTJk0qffyDBw8wcOBAuLq64qeffoJcLq/0fIVCAYVCYVL7CakJmopM7JlVg4u3tze8vb2rPK979+7Izc1FYmIiwsPDAQCHDx+GTqdDt27dKnxcfn4+oqKioFAosHv3bjg5OYnWdkIIIRWrEzmX9u3bY+DAgZg0aRJOnjyJo0ePYsqUKXjxxRf5mWIZGRlo164dTp48CUAfWAYMGICCggKsW7cO+fn5yMrKQlZWFrRarTVfDiGE2L06s0J/y5YtmDJlCvr27QupVIrnn38eq1at4u/XaDRISUlBYWEhAODMmTP8TLLWrVsLnuvGjRsICAiotbYTQkh9U2eCi6enJ7Zu3Vrh/QEBATBcshMZGYk6sISHEELsUp0YFiOEEFK3UHAhhBAiOgouhBBCREfBhRBCiOgouBBCCBEdBRdCCCGio+BCCCFEdBRcCLEy2iyM2CMKLoRY2dqE6yjSaLE24bq1m0KIaCi4EGJlkyNb0WZhxO7UmfIvhNir6Ah/yG40xCjaLIzYEeq5EEIIER0FF0IIIaKj4EIIIUR0FFwIIYSIjoILIYQQ0VFwIYQQIjoKLoQQQkRHwYUQQojoKLgQQggRHQUXQgghoqPyL1VgjAEA8vPzrdwSYs8KCwvpd4zUGu53jft8swQJs+Sz24F//vkHfn5+1m4GIYSI7ubNm2jWrJlFnpuCSxV0Oh0yMzPh6uoKiURi7eYI5Ofnw8/PDzdv3oSbm5u1m2MR9BrtA71G28IYw4MHD+Dj4wOp1DLZERoWq4JUKrVYZBeLm5ubzf8y1xS9RvtAr9F2qFQqiz4/JfQJIYSIjoILIYQQ0VFwqcMUCgViYmKgUCis3RSLoddoH+g11j+U0CeEECI66rkQQggRHQUXQgghoqPgQgghRHQUXAghhIiOgosdKi4uRmhoKCQSCZKSkqzdHNGkpqZi4sSJaNGiBZRKJVq1aoWYmBiUlJRYu2k1smbNGgQEBMDJyQndunXDyZMnrd0k0SxfvhxdunSBq6srGjVqhGHDhiElJcXazbKoFStWQCKRYMaMGdZuilVRcLFDb7/9Nnx8fKzdDNFdvnwZOp0OX331FS5evIjPPvsMX375Jd555x1rN63atm/fjpkzZyImJgZnzpxBSEgIoqKicOfOHWs3TRS///473nzzTZw4cQIHDhyARqPBgAEDUFBQYO2mWcSpU6fw1VdfoVOnTtZuivUxYld++eUX1q5dO3bx4kUGgJ09e9baTbKojz76iLVo0cLazai2rl27sjfffJO/rdVqmY+PD1u+fLkVW2U5d+7cYQDY77//bu2miO7BgwesTZs27MCBA6xXr15s+vTp1m6SVVHPxY7cvn0bkyZNwnfffQdnZ2drN6dW5OXlwdPT09rNqJaSkhIkJiaiX79+/DGpVIp+/frh+PHjVmyZ5eTl5QFAnf03q8ybb76JQYMGCf496zMqXGknGGMYN24cXn/9dXTu3BmpqanWbpLFXbt2DatXr0ZsbKy1m1Itd+/ehVarRePGjQXHGzdujMuXL1upVZaj0+kwY8YM9OzZEx07drR2c0T1/fff48yZMzh16pS1m2IzqOdi4+bNmweJRFLpz+XLl7F69Wo8ePAA8+fPt3aTzWbqazSUkZGBgQMHYsSIEZg0aZKVWk7M8eabb+LChQv4/vvvrd0UUd28eRPTp0/Hli1b4OTkZO3m2Awq/2LjsrOzce/evUrPadmyJf7973/j559/Fuw5o9VqIZPJMHr0aHz77beWbmq1mfoaHR0dAQCZmZmIjIxEREQENm7caLH9KCytpKQEzs7O+OGHHzBs2DD++NixY5Gbm4tdu3ZZr3EimzJlCnbt2oUjR46gRYsW1m6OqOLi4vDcc89BJpPxx7RaLSQSCaRSKYqLiwX31RcUXOxEenq6YJvczMxMREVF4YcffkC3bt1sfk8aU2VkZKB3794IDw/H5s2b6/x/2m7duqFr165YvXo1AP3QUfPmzTFlyhTMmzfPyq2rOcYYpk6dip9++gkJCQlo06aNtZskugcPHiAtLU1wbPz48WjXrh3mzp1rd0OApqKci51o3ry54HaDBg0AAK1atbKrwBIZGQl/f3/ExsYiOzubv69JkyZWbFn1zZw5E2PHjkXnzp3RtWtXrFy5EgUFBRg/fry1myaKN998E1u3bsWuXbvg6uqKrKwsAPqNqpRKpZVbJw5XV9dyAcTFxQVeXl71NrAAFFxIHXLgwAFcu3YN165dKxcw62oHfOTIkcjOzsaiRYuQlZWF0NBQ7N27t1ySv65au3YtACAyMlJwfMOGDRg3blztN4jUGhoWI4QQIrq6mQklhBBi0yi4EEIIER0FF0IIIaKj4EIIIUR0FFwIIYSIjoILIYQQ0VFwIYQQIjoKLoRYUEJCAiQSCXJzc63WhnHjxglqlxkTEBCAlStX8rezsrLQv39/uLi4wN3d3exrvvzyy1i2bJnZj7Nle/fuRWhoKHQ6nbWbUidQcKlHsrKyMH36dLRu3RpOTk5o3LgxevbsibVr16KwsNDazbNLPXr0wK1bt6BSqazdlEqdOnUKr776Kn/7s88+w61bt5CUlIQrV66YFSTPnTuHX375BdOmTbNgi2vfwIEDIZfLsWXLFms3pU6g4FJP/O9//0NYWBj279+PZcuW4ezZszh+/Djefvtt7NmzBwcPHrR2E62ipKTEos/v6OiIJk2aCKpV2yJvb2/BBnPXr19HeHg42rRpg0aNGpn1XKtXr8aIESP4+na2TqvVmtwbGTduHFatWmXhFtkJq+2BSWpVVFQUa9asGXv48KHR+3U6Hf/3+/fvs4kTJ7KGDRsyV1dX1rt3b5aUlMTfHxMTw0JCQtimTZuYv78/c3NzYyNHjmT5+fn8OVqtli1btowFBAQwJycn1qlTJ7Zjx45K21hUVMRmzZrFfHx8mLOzM+vatSv77bffGGOMqdVqFhQUxCZNmsSff+3aNdagQQO2bt06xhhjGzZsYCqViv3000+sdevWTKFQsAEDBrD09PRybf/mm29YQEAAk0gkJr3mpKQkFhkZyRo0aMBcXV3ZE088wU6dOsUYYyw1NZX961//Yu7u7szZ2ZkFBQWx+Ph4xhhjv/32GwPA7t+/zz/XDz/8wIKCgpijoyPz9/dnsbGxgvfB39+fLV26lI0fP541aNCA+fn5sa+++qrS927Hjh2sY8eOzMnJiXl6erK+ffvy/9Zjx45lQ4cOZR9//DFr0qQJ8/T0ZG+88QYrKSkRXPOzzz7j/w6A/xk7dqzgNnfMmNLSUqZSqdiePXv4Y9x7YOpzcI/p0qULc3Z2ZiqVivXo0YOlpqby9+/evZt17tyZKRQK5uXlxYYNG8bfl5OTw15++WXm7u7OlEolGzhwILty5Qp/P/d7smvXLta+fXsmk8nYjRs3Kv3946SlpTEA7Nq1a5X+exDGKLjUA3fv3mUSicTkfdn79evHBg8ezE6dOsWuXLnCZs2axby8vNi9e/cYY/oP6AYNGrDhw4ez8+fPsyNHjrAmTZqwd955h3+ODz74gLVr147t3buXXb9+nW3YsIEpFAqWkJBQ4XVfeeUV1qNHD3bkyBF27do19vHHHzOFQsF/MJw9e5Y5OjqyuLg4VlpayiIiIthzzz3HP37Dhg1MLpezzp07s2PHjrHTp0+zrl27sh49evDnxMTEMBcXFzZw4EB25swZdu7cOZNec4cOHVh0dDS7dOkSu3LlCvvvf//LB59Bgwax/v37s+TkZHb9+nX2888/83vElw0up0+fZlKplL333nssJSWFbdiwgSmVSrZhwwa+jf7+/szT05OtWbOGXb16lS1fvpxJpVJ2+fJlo+9bZmYmc3BwYJ9++im7ceMGS05OZmvWrGEPHjxgjOmDi5ubG3v99dfZpUuX2M8//8ycnZ3Z119/LbgmF1zu3LnDBg4cyP7973+zW7dusdzcXPbjjz8yACwlJYU/ZsyZM2cYAJaVlcUfKy4uZrdu3eJ/Dh8+zJycnPgvBWVpNBqmUqnY7Nmz2bVr19jff//NNm7cyNLS0hhjjO3Zs4fJZDK2aNEi9vfff7OkpCS2bNky/vFDhgxh7du3Z0eOHGFJSUksKiqKtW7dmg+m3O9Jjx492NGjR9nly5dZQUFBlb9/nMaNGwv+vYhxFFzqgRMnTjAAbOfOnYLjXl5ezMXFhbm4uLC3336bMcbYH3/8wdzc3FhRUZHg3FatWvHfnmNiYpizs7OgpzJnzhzWrVs3xpi+B+Ls7MyOHTsmeI6JEyeyUaNGGW1jWloak8lkLCMjQ3C8b9++bP78+fztjz76iDVs2JBNmTKFNW3alN29e5e/b8OGDQwAO3HiBH/s0qVLDAD766+/+LbL5XJ2584d/hxTXrOrqyvbuHGj0bYHBwezxYsXG72vbHB56aWXWP/+/QXnzJkzhwUFBfG3/f39WXR0NH9bp9OxRo0asbVr1xq9RmJiIgMg+GZvaOzYsczf35+Vlpbyx0aMGMFGjhwpuCYXXBhjbOjQoYKehbEemDE//fQTk8lkgp6wobt377KWLVuyN954o8LnuHfvHgNQ4ReR7t27s9GjRxu978qVKwwAO3r0qOCaSqWS/fe//2WMPf49MeyZmvr7xxhjYWFhFf57k8co51KPnTx5EklJSejQoQOKi4sB6JOxDx8+hJeXFxo0aMD/3LhxA9evX+cfGxAQAFdXV/5206ZNcefOHQD6ve0LCwvRv39/wXNs2rRJ8ByGzp8/D61Wi7Zt2woe8/vvvwseM2vWLLRt2xaff/451q9fDy8vL8HzODg4oEuXLvztdu3awd3dHZcuXeKP+fv7w9vbm79tymueOXMmXnnlFfTr1w8rVqwQtGnatGn44IMP0LNnT8TExCA5ObnC9/zSpUvo2bOn4FjPnj1x9epVaLVa/linTp34v0skEjRp0oR/f8sKCQlB3759ERwcjBEjRuCbb77B/fv3Bed06NBBsLGa4b+XmNRqNRQKhdEck0ajwfPPPw9/f3/85z//AaDf5M7wPV+2bBk8PT0xbtw4REVFYfDgwfjPf/6DW7du8c+TlJSEvn37Gr3+pUuX4ODggG7duvHHvLy8EBgYKPgdcHR0FLzHpv7+AYBSqaQJMCag/VzqgdatW0MikSAlJUVwvGXLlgAg2LTp4cOHaNq0KRISEso9j+GUVLlcLrhPIpHwSdGHDx8CAOLj4+Hr6ys4T6FQGG3jw4cPIZPJkJiYWG53ScPE8J07d3DlyhXIZDJcvXoVAwcONPp8lXFxcSl37ape8+LFi/HSSy8hPj4ev/76K2JiYvD999/jueeewyuvvIKoqCjEx8dj//79WL58OT755BNMnTrV7LZxKnt/y5LJZDhw4ACOHTuG/fv3Y/Xq1Xj33Xfx119/8VsKm/N8NdGwYUMUFhaipKSE35aaM3nyZNy8eRMnT56Eg4P+o8fHxwdJSUn8OZ6engD0+71MmzYNe/fuxfbt27FgwQIcOHAAERERomwyplQqBQHQ1N8/AMjJyRF8OSHGUc+lHvDy8kL//v3x+eefo6CgoNJzn3jiCWRlZcHBwQGtW7cW/DRs2NCk6wUFBUGhUCA9Pb3cc/j5+Rl9TFhYGLRaLe7cuVPuMYa7TE6YMAHBwcH49ttvMXfuXMG3UQAoLS3F6dOn+dspKSnIzc1F+/bta/ya27Zti7feegv79+/H8OHDsWHDBv4+Pz8/vP7669i5cydmzZqFb775xui12rdvj6NHjwqOHT16FG3btq3Rls0SiQQ9e/bEkiVLcPbsWTg6OuKnn36q9vOVxQUKw96VMaGhoQCAv//+W3D8008/xX//+1/s2rVL0Nss+55zwQXQ/07Mnz8fx44dQ8eOHbF161YA+l7doUOHjF6/ffv2KC0txV9//cUfu3fvHlJSUhAUFFRhu039/SsqKsL169cRFhZW6ftAKLjUG1988QVKS0vRuXNnbN++HZcuXUJKSgo2b96My5cv8x9s/fr1Q/fu3TFs2DDs378fqampOHbsGN59913Bh3ZlXF1dMXv2bLz11lv49ttvcf36dZw5cwarV6/Gt99+a/Qxbdu2xejRozFmzBjs3LkTN27cwMmTJ7F8+XLEx8cDANasWYPjx4/j22+/xejRozFs2DCMHj1aMJ1YLpdj6tSp+Ouvv5CYmIhx48YhIiICXbt2rbC9Vb1mtVqNKVOmICEhAWlpaTh69ChOnTrFB6wZM2Zg3759uHHjBs6cOYPffvutwmA2a9YsHDp0CO+//z6uXLmCb7/9Fp9//jlmz55t0ntrzF9//YVly5bh9OnTSE9Px86dO5GdnV1pQDWXv78/JBIJ9uzZg+zsbL53Wpa3tzeeeOIJ/Pnnn/yxgwcP4u2338bHH3+Mhg0bIisrC1lZWcjLyzP6HDdu3MD8+fNx/PhxpKWlYf/+/bh69Sr/emJiYrBt2zbExMTg0qVLOH/+PD788EMAQJs2bTB06FBMmjQJf/75J86dO4fo6Gj4+vpi6NChFb4+U37/AODEiRNQKBTo3r272e9hvWPtpA+pPZmZmWzKlCmsRYsWTC6XswYNGrCuXbuyjz/+mBUUFPDn5efns6lTpzIfHx8ml8uZn58fGz16ND+ll5vOa+izzz5j/v7+/G2dTsdWrlzJAgMDmVwuZ97e3iwqKoqfRWVMSUkJW7RoEQsICGByuZw1bdqUPffccyw5OZldunSJKZVKtnXrVv78+/fvMz8/P34yAjfF9Mcff2QtW7ZkCoWC9evXj59lVFHbq3rNxcXF7MUXX2R+fn7M0dGR+fj4sClTpjC1Ws0YY2zKlCmsVatWTKFQMG9vb/byyy/zEw0qm4osl8tZ8+bN2ccffyxoS9nkOmOMhYSEsJiYGKPv299//82ioqKYt7c3UygUrG3btmz16tX8/dxUZEPTp09nvXr1qvCaZRP6jDH23nvvsSZNmjCJRFLpNOIvvviCRURE8LdjYmLMmoqclZXFhg0bxpo2bcpP1160aBHTarX8OT/++CMLDQ1ljo6OrGHDhmz48OH8fdxUZJVKxZRKJYuKijI6Fbmsyn7/OK+++ip77bXXKnzt5DHa5pjYjY0bN2LGjBlWLbVC9En9wMBAbN++3a6+4d+9exeBgYE4ffo0n8siFaNhMUKIqJRKJTZt2oS7d+9auymiSk1NxRdffEGBxUQ0W4wQIrrIyEhrN0F0nTt3RufOna3djDqDhsUIIYSIjobFCCGEiI6CCyGEENFRcCGEECI6Ci6EEEJER8GFEEKI6Ci4EEIIER0FF0IIIaKj4EIIIUR0FFwIIYSI7v8BC668gSlmENkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(4, 4), dpi=100)\n",
    "\n",
    "ax.axhline(c='grey', lw=0.5)\n",
    "ax.axvline(c='grey', lw=0.5)\n",
    "ax.scatter(pert_match_df['gene_shift_z'], pert_match_df['pert_match_score'], s=1, rasterized=True)\n",
    "\n",
    "label_df = pert_match_df[pert_match_df['perturbed_gene_name'].isin(\n",
    "    ['EOMES', 'TBXT', 'EVX1', 'SNAI1', 'POU5F1', 'NANOG', 'PRDM14', 'DNMT3B'])]\n",
    "for i, row in label_df.iterrows():\n",
    "    plt.text(row['gene_shift_z'], row['pert_match_score'], row['perturbed_gene_name'], \n",
    "             fontsize=8)\n",
    "    \n",
    "ax.set_xlabel('Gene expression shift (z-score)')\n",
    "ax.set_ylabel('Perturbation match score')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we can see that the EOMES, TBXT, EVX1 and SNAI1 are predicted to be the positive causal genes for the cell state transition. Indeed, these are transcription factors known to be crucial for primitive streak and mesoderm development. The transcription factors POU5F1 and NANOG, which are responsible for maintaining the pluripotency state, are predicted to be negative causal genes."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "scmg",
   "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.10.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}