{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "intro",
   "metadata": {},
   "source": [
    "# Calling hits against genes that cannot respond\n",
    "\n",
    "A genome-wide knockout screen asks which genes change the cell when they are removed. Most do not, so the work is telling the genes that moved the cell from the ones that only look like they did. A per-feature test does not settle the question on its own: with thousands of control wells behind it, almost any gene differs from the controls on some feature, and the number that reaches a fixed p-value says more about the test's power than about biology.\n",
    "\n",
    "PERISCOPE calibrates the call against genes a cell line does not express. DepMap reports expression for the cell line, and a gene it records at zero TPM cannot change anything when it is knocked out, so the scores those genes reach are an empirical null. A real knockout counts as a hit when it scores past what the non-expressed genes reach. This case study builds the per-gene score, reads the null from DepMap, and calls hits with {func}`~mantispy.tl.empirical_fdr` on the JUMP CRISPR screen, which was run in U2OS cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "imports",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:44:22.946587Z",
     "iopub.status.busy": "2026-10-21T00:44:22.946060Z",
     "iopub.status.idle": "2026-10-21T00:45:22.295242Z",
     "shell.execute_reply": "2026-10-21T00:45:22.293094Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import fisher_exact\n",
    "\n",
    "import mantispy as mt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "load",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:45:22.299505Z",
     "iopub.status.busy": "2026-10-21T00:45:22.298935Z",
     "iopub.status.idle": "2026-10-21T00:45:23.980306Z",
     "shell.execute_reply": "2026-10-21T00:45:23.979177Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(51185, 595) {'treatment': 43138, 'negcon': 7478, 'poscon': 569}\n"
     ]
    }
   ],
   "source": [
    "wells = mt.ds.jump_crispr()\n",
    "kinds = wells.obs[\"Metadata_Control_Type\"].astype(str).value_counts().to_dict()\n",
    "print(wells.shape, kinds)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "score-intro",
   "metadata": {},
   "source": [
    "## A score that rises with the phenotype\n",
    "\n",
    "The score is how many features separate a gene's wells from the negative controls. For each feature it is a Mann-Whitney test of the gene's wells against the control wells, and the gene's score is the count of features that come out significant. The per-feature Mann-Whitney test is the one PERISCOPE uses, here on JUMP's Cell Painting features rather than PERISCOPE's own. The count is unbounded, so a strong knockout keeps climbing rather than saturating, which matters for the calibration below. The gene is the unit, with its several guides as replicate wells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "score",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:45:23.983729Z",
     "iopub.status.busy": "2026-10-21T00:45:23.983481Z",
     "iopub.status.idle": "2026-10-21T00:49:56.788412Z",
     "shell.execute_reply": "2026-10-21T00:49:56.786578Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7974 genes scored, median 5, max 431\n"
     ]
    }
   ],
   "source": [
    "# For each feature, a Mann-Whitney test of every gene's wells against the negative controls, via\n",
    "# effect_size; the gene's score is the count of features significant at p<0.001, the test PERISCOPE uses.\n",
    "mt.tl.effect_size(wells, groupby=\"Metadata_Perturbation\", reference=\"negcon\", pvalues=True, exact=False)\n",
    "table = wells.uns[\"mantispy\"][\"effect\"]\n",
    "treatment_genes = (\n",
    "    wells.obs.loc[wells.obs[\"Metadata_Control_Type\"].astype(str) == \"treatment\", \"Metadata_Perturbation\"]\n",
    "    .astype(str)\n",
    "    .unique()\n",
    ")\n",
    "counts = table.assign(hit=table[\"pvalue\"] < 0.001).groupby(\"group\", observed=True)[\"hit\"].sum()\n",
    "counts = counts[counts.index.astype(str).isin(treatment_genes)]\n",
    "genes = counts.index.to_numpy().astype(str)\n",
    "score = counts.to_numpy().astype(float)\n",
    "print(f\"{len(genes)} genes scored, median {np.median(score):.0f}, max {score.max():.0f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "controls-intro",
   "metadata": {},
   "source": [
    "## The genes U2OS cannot express\n",
    "\n",
    "The null is the genes the screened cell line does not express, read from DepMap. {func}`~mantispy.ds.depmap_expression` fetches the expression matrix (`OmicsExpressionProteinCodingGenesTPMLogp1.csv`) from DepMap's own public release (CC BY 4.0) on first use, and you point {func}`~mantispy.io.unexpressed_genes` at it. DepMap keys its cell lines by a model id, so the screen's line, U2OS, is named by its id rather than by a looked-up name. A gene the model reports at zero TPM cannot respond to a knockout, so these genes are the empirical null."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "controls",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:49:56.792762Z",
     "iopub.status.busy": "2026-10-21T00:49:56.792447Z",
     "iopub.status.idle": "2026-10-21T00:50:08.572266Z",
     "shell.execute_reply": "2026-10-21T00:50:08.571185Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2137 genes at zero TPM in U2OS, 559 of them screened here\n"
     ]
    }
   ],
   "source": [
    "# DepMap's expression matrix, fetched from DepMap's own public release on first use.\n",
    "# U2OS, the line this screen was run in, is model ACH-000364 in DepMap's own key.\n",
    "controls = mt.io.unexpressed_genes(mt.ds.depmap_expression(), \"ACH-000364\")\n",
    "present = sorted(set(genes) & controls)\n",
    "print(f\"{len(controls)} genes at zero TPM in U2OS, {len(present)} of them screened here\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "scores-fig",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:08.575780Z",
     "iopub.status.busy": "2026-10-21T00:50:08.575583Z",
     "iopub.status.idle": "2026-10-21T00:50:09.967556Z",
     "shell.execute_reply": "2026-10-21T00:50:09.966255Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAE2CAYAAACukdFuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPhlJREFUeJzt3XlYVGX/BvB7WBzZQUVBAVHckTBcADXFLXFL00xIsVQ0SV8rscIyS1vcs1JTy71UEjHxxX0XEXIBMcQdEQJElF0QFM7vj5+c1wkGZobBOcb9uS6umnOe55nvc0bm5iwzRyYIggAiIiKSJD1dF0BERETKMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwA10XIFVlZWVIS0uDmZkZZDKZrsshIqJ/EUEQkJ+fj6ZNm0JPr+p9Zga1EmlpabC3t9d1GURE9C+WkpICOzu7KtswqJUwMzMDnm5Ec3NzXZdDRET/Inl5ebC3txezpioMaiXKD3ebm5szqImIqFaocmqVF5MRERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiojqkfv36ePToka7LUJu3tzcSEhJ0XUalnjx5gq5duyIrK6tWxucXnhDRv9qBAwdqZVxvb+9aGbc6/fv3R1BQEPr376+T59eFU6dOwdDQEB06dNDKeJ06dUJcXJz4eMmSJZg1a1a160JDQ/Hpp58iOTkZ3bt3x+bNm2FnZwcDAwO88847WLZsGb755hut1PgsBvVzUNUbha5+2YmIXhQbNmzAmDFjtDpmbGwsOnXqpPK61NRUTJgwAbt27UL37t2xdOlSTJ8+Hbt37wYAvPnmm+jUqRO+/vprrd/IiYe+iYhq0d9//43XX38dVlZWkMlk8Pf3BwBcv34dXl5eMDMzQ8eOHRX+oPfy8sKXX36J9u3bw9LSEnPmzAEALF26FEePHsWAAQMgk8kQHh4OLy8vfP7552jfvj169eoFANi7dy86dOgAMzMz9O3bFzdv3lSp1pCQEDg6OsLW1hZLly6FgcH/9uXOnDmDbt26wcTEBM7Ozti/f7+4TiaTYfny5WjatCmaNWsmhldV/UpKSjB27FiYm5tDJpPB0dFRaV3Hjx+Hp6enSs9XWxISEuDo6Ij+/fvD2NgY7777Lv773/8iPz8fAGBtbQ0TE5NaOTzPoCYiqkW+vr5wcHBAYmIiBEHAunXrxOV9+/ZFeno6vvnmG/j4+CAzM1Psd/bsWZw8eRKxsbH45ZdfcOvWLcyaNQv9+vXD4cOHIQgChg4dCgCIjo7GiRMncOrUKdy9exdjx47FsmXLkJaWBg8PD4wbN67aOtPT0zFlyhRs3LgR165dQ3JysrguKysLAQEBWLt2LbKzs7Fp0yZMnjwZhYWFYpvs7Gxcv34dv/zyC2bMmFFtv0OHDiE9PR0pKSkQBAFJSUmV1lVaWoqUlBQ0a9ZMYXllzwcAHTt2hEwmq/SnoKBAbOfl5QULCwsMHjwYt27dUhi7snXt27fH7du3cfz4cRQWFuKXX35BWVkZUlJSxH729vZITEysdluri0FNRFRLcnNzce7cOSxbtgxWVlbi8uzsbCQkJGDOnDkwNTXF8OHD8fLLLyM6OlpsM2PGDDRu3BgtWrSAm5tblQEwY8YMNGnSBAAQFRWFbt26YdCgQTAzM8O8efMQExMj7vkpExUVhV69eqFPnz4wNzfHl19+Ka47ffo0Ll26BDc3N8jlcnTr1g2pqakKNX366acwNTXF4MGDcf/+fTx+/LjKfq1bt0ZCQgI+++wz7Ny5UyH0VVHZ8wFAfHw8BEGo9MfU1BQAcPHiReTk5CAxMRGtWrXC2LFjxXGVrbOzs8PPP/8Mf39/WFtbo7CwEGZmZjA0NBT7CoKg9cPeYFATEdUeQRDUavvsm3z9+vXF/9fX10dpaSmg5G5Lz94qUZ3nVLVWQRDg5eVVIfg6duxYab16enooLS2tsl/btm1x9uxZtG7dGlu2bMFLL72EoqKiCs+tr68Pe3t7pKamKiyv7Pmgxh41ADRs2BBBQUGIiYmpMP/K1vn6+uLWrVt4+PAhfHx88OTJE4VD9qmpqWjRooVK21sdDGoiolpiaWmJrl274uOPP0ZOTo643MrKCh06dMDXX3+NgoIC/Pe//8XFixfh7u5e7ZgWFhZV7l13794dZ8+exf79+1FQUIB58+bBzc2t2vsed+/eHadOncLJkyeRn5+P+fPni+t69OiB+Ph4rF+/Xq093+r6OTg44P3338fKlSuRm5urcOj/WX369FE42lAVVfaoyz148AALFy5E586dK/wBVNW6mzdvYurUqZgxY4a4R33//n0UFBRo7cr0ZzGoiYhq0bZt23Dz5k04ODgoXEy2fft2HDt2DDY2NggKCkJwcDCsra2rHW/KlCmYP38+9PT0EB4eXmG9jY0Ntm7disDAQNjY2ODMmTP47bffqh3X1tYWa9aswfjx49GmTRs0aNAAFhYWAIBGjRohPDwcmzZtgrW1NWQyWYXQq0xV/dasWQOZTAY9PT14eHjgww8/hIODQ6XjTJo0CcHBwdU+nyquXr0q7mG3atUKSUlJ2LZtW7XrAKBLly6QyWTo1q0bOnfurPDHzI4dO/D222/XyqFvmaDpcZJ/uby8PFhYWCA3Nxfm5uY1GosfzyKiF0leXh4WL16M2NhY7N27V9flAE/fK7/77rta2WOtqdLSUnh4eODAgQNo2LChSn3UyRh+jpqIiAAAb7zxBkJDQ2FiYgIPDw/8/PPPui5JVFtfXKMN+vr6OHfuXK2Nz6AmIiIAwM6dO3VdAlWC56iJiIgkjEFNREQkYQxqIiIiCdN5UIeGhqJ3795o164dRo8ejWvXrtW4T1FREb7++mt0794d7du3x8CBA7Fr165anAUREVHt0GlQh4WFwcfHB6NHj8bWrVtRr1499OrVC/fv369Rn//85z/4+eef8fnnn2PXrl3w9vbG6NGjK/3MIRERkZTpNKjnz58PPz8/TJ8+HZ07d8amTZtQWlqKNWvW1KjP0aNHMXHiRAwaNAjt27fHhx9+CGdnZxw9evQ5zYyISH3vvfceYmJidF3GcxcaGqrSl7JIUVBQkEpHgmtCZx/Pys/PR0xMDD7++GNxmaGhIfr164cTJ06It3XTpE/fvn1x5MgRzJw5E+bm5rh48SISExOxYMGC5zQ7Rdv0lL+I3uAXnhDVps+u1M5pr2/aj6zxGDNmzICfnx+6du0KPL2VYl5enhaqe3GUlpZizpw5iIyM1Mp4Pj4+Fe7EtXz5coXbZCYnJ8PX1xeTJ0/GO++8AwBYsWIFtm7dCgCYNm0a/Pz8xPbff/+9wjej9ezZE0uXLgUAvPrqq5gzZw5CQkK0Un9ldBbUf//9N/D06+6eZWNjg0uXLtWoz9q1a+Hn5wdra2tYWVkhNzcXq1evxpAhQ5TWU1xcjOLiYvFxXftlIaLnLyEhAbm5ubouQ6f27dsHZ2dnNGjQQCvjzZ49W7y5R3JyMvz9/fHSSy+J6wVBQEBAAGQymZgpADB8+HB07doVCxYsqHADkKSkJHh6emLMmDHA0+9qL9enTx9MmjQJ6enpsLW11coc/klnh77LysoAQOHG5Hi6h1x+FxRN+8yePRtRUVHYuXMnjh07hgULFmD69Ok4fvy40noWLFgACwsL8cfe3r5G8yMiAoCCggIEBQVhwIAB8Pf3x+3bt4Gn3wF+7tw5zJgxAx4eHuIeZUpKCiZMmABvb2+Fi2BLSkqwcOFCeHt7Y9SoUQqn8t577z3s2rULEyZMQEBAQKV15OTkYPr06Rg0aBBWr16N9957T9zBqW7ssLAwtWs6dOgQhg8fDk9PT3h4eCA9Pb3Sug4cOIBevXqp9HyqcHV1hYeHBzw8PHD58mX4+PjAxMREXL9mzRp4enoq3PkLT28Q4uHhofT71ps3by6O27ZtW3G5TCZDjx49cPjwYbXqVIfO9qgbNWoEPL1DybPu37+vdEOp0icnJwffffcdNm/ejGHDhgEAOnTogNOnT+Obb75Bnz59Kh179uzZmDlzpvg4Ly+PYU1ENebn5wc9PT188sknOHr0KLy8vHD16lX06dMHbdu2FQ99t2vXDgCwbNkyzJ8/Hw8fPsTEiRPRo0cPNGnSBFOmTIG1tTWCgoKQmZmJKVOmIDw8HO3bt0dCQgIiIiLw9ddfo3Xr1pXWMXbsWFhZWSEwMBDh4eFYv3493nrrLeDpjT6qGvvMmTNq1dSiRQu89dZbWLVqFZo3bw4ASveYr1y5gkGDBomPq3q+ZcuWKT3EvGzZMvTo0UN8XFZWhk2bNmHHjh3istu3byM4OBhHjx7F9OnT1XodN2zYgD179qBdu3b4+OOPFW5v2bJlSyQkJKg1njp0FtRNmjSBg4MDzpw5g9dee01cHhkZqfQQtSp9Hj9+jLKyMlhaWir0tbCwUPoXHQDI5XLI5XItzIyI6P8VFhZi7969yMrKgqmpKfr374+DBw8iMjIS/fv3h7m5Odq3bw8PDw+xz9y5c8X3tw0bNiAhIQEWFhbYunUrOnfujIiICODpnvqpU6fQvn17AMCnn36K4cOHV1rHw4cPcfz4cWRnZ0Mul6Nfv34ICwsDADx69KjasdWtqUWLFtDX18fdu3fRoUMHuLi4KN1GxcXFFd57K3u+Jk2aYNSoUQph/KzyP3TKHTx4EGZmZuKtQwVBwNSpU7Fy5coKR2Wr8+GHH8LHxweFhYUICQlB7969ceXKFRgbGwNP741dm6dLdfpd3++99x6WLFmCcePGwdnZGatXr0ZSUhImT54stgkKCkJsbCwOHjyoUh9ra2u4urqKf11ZWVkhNjYWoaGh+PDDD3U2VyKqe/Lz8yGXyxVuCdm4ceMqz0s/u+cpl8tRXFyMvLw86OvrY/ny5Qq3UXx2r66q86P5+fkwNTUVA1Emk4lHKFUZW92a6tevj0OHDuHnn38Wb2d56NAh2NnZVajN1ta2wn2oK3s+PL06XNU96nXr1mHSpEni4/PnzyM6OlrMitu3b6NevXoQBAGff/650m2Hp4e9y48M9O3bF05OToiLixMvULt37x6aNm1a5Rg1odOg/uijj5CamoouXbpALpejfv362LZtG5ydncU29+/fVzixr0qfkJAQTJ06Fba2tjA3NxcPn8yePfu5z5GI6q4mTZrAwsIC4eHhGDp0KJKSkhAdHY1Vq1YBAIyNjVFQUFDtOI0bN4adnR3i4+MVdmTUqcPAwAAnT55E7969cfXqVfH8tKZjV9fP1dVVnOfw4cNx6NAhTJw4sUK7Hj16ICYmRjwMXxVV96jv3buHgwcPYu3ateKy9u3bizt8ALB48WI0btxYped91uXLl5GRkaEQzOfPn8eSJUvUGkcdOg1qPT09/Pjjj1i0aBFycnLQpEkT6OkpXt+2aNEihauxVenTunVrHD16FCUlJcjOzkbjxo1r5WbeRETVWb9+PcaNG4cmTZogJSUFn332GZycnAAAgwYNgr+/PxYuXIhly5ZVOc62bdvg4+ODRYsWiXvDmzZtqnDItzIymQyrVq3Ca6+9BgcHB9SrVw/NmzeHkZFRjcZW1s/IyEi8Qjo3NxdlZWVYvXp1pWP4+PigT58+4sedquLo6Kiwp6/M5s2bMXjwYLEmADA1NVU4xVD+h0b5a7Fnzx58++23SExMRP369bF79278/PPPeOmllzBkyBA8ePAAhYWFSEpKwldffSXuYd+9exeZmZkKY2ubTBAEodZGf4Gpc1Pv6ow/9IPSdVtefb9GYxOR9BUVFeHWrVuws7OrcP3M33//jbS0NLRp0wZpaWmws7MT33OuXbsGGxsbWFhYAE8/c5yYmCheUOvi4gITExMkJCQo9FMmLy8PaWlpMDQ0hLu7O27evCnWo+rYqtSkr6+PixcvQiaTwcLCAm3atKmwQ/Wsd999F8OGDcPQoUOrfT5VXL9+Hebm5hU+yvus8kPfzZo1AwBkZmbi1q1bCm2cnZ1hZmaGmJgYlJSUwNjYGC1btlQ4lTFnzhy0bNmy0qMFVVEnYxjUSjCoiejfZOfOnVi6dCmKi4tx69YtzJ8/Hx988IGuywKevt8+ePAALVq00HUpavvrr7/g7Oxc5R8ilVEnY3R66JuIiJ6PV155BXZ2dqhXrx5atmxZYc9el8zNzWu8Q6QrVV3Rri0MaiKiOqBJkyZo0qSJrssgDej8NpdERESkHIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCSMQU1ERCRhBrou4NKlS1i7di0yMjLg4uKC999/H5aWljXu8+TJE2zZsgXHjh2DsbExJk2aBHd391qeDRERkXbpdI/67NmzcHd3R1lZGYYNG4aDBw+iZ8+eKCoqqlGfR48eoW/fvli6dCl69+6N3r174+OPP0ZMTMxzmhkREZF2yARBEHT15P369YOpqSnCwsIAADk5ObCzs8PChQsxffp0jfvMnTsXq1evxpUrV9CoUSMAQFlZGR4+fAgzMzOVasvLy4OFhQVyc3Nhbm5eo3mOP/SD0nVbXn2/RmMTEdGLR52M0dke9aNHj3Dy5EmMHDlSXGZpaYl+/frhwIEDNeqzceNG+Pn5iSENAHp6eiqHNBERkVToLKiTk5NRWloKOzs7heX29va4ffu2xn2ysrLw999/o1OnTli4cCFGjx6N999/H+fPn6+ynuLiYuTl5Sn8EBER6ZrOgrqkpAQAYGxsrLDc2NhYXKdJn8LCQgBAUFAQ0tLSMHr0aBgYGMDDwwN79uxRWs+CBQtgYWEh/tjb29dwhkRERDWns6u+y6/SzsrKUlj+4MEDpVd9q9Kn/L8eHh748ccfAQBvvvkm0tPTsWDBArz22muVjj179mzMnDlTfJyXl8ewJiIindPZHnWzZs3QsGFDxMXFKSy/ePEiXF1dNe5jamqKVq1awcnJSaFNy5Ytce/ePaX1yOVymJubK/wQERHpms6CWiaTwc/PD+vWrcODBw8AAIcPH0ZMTAzGjx8vtluxYgWmTZumVp+JEydiz5494nnmgoIChIWFoWfPns95lkRERDWj0y88+eqrrxAXF4e2bduiTZs2iI2Nxfz589GrVy+xTVxcHKKjo9XqM2vWLFy8eBFOTk5wcXHB5cuX0a5dO3z33XfPfY5EREQ1odPPUZe7dOkSMjIy4OzsjKZNm1ZYl5OToxDE1fUpd+PGDdy5cwcODg5o06aNWjXxc9RERFRb1MkYnX+FKAC89NJLaq+rqk+51q1bo3Xr1jWqjYiISJd4Uw4iIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYRpFNSPHz/WfiVERERUgUZBbWdnh1mzZiEhIUH7FREREZFIo6D+8ssvcfLkSTg7O8PDwwO//PIL8vPztV8dERFRHadRUAcEBODcuXO4dOkSPD098emnn8LGxgZvv/02Tp48CQl8KykREdG/Qo0uJnNxccHy5cuRlpaGb775BsHBwfDy8kKbNm2watUqPHnyRHuVEhER1UE1+q7v4uJihIWFYcOGDTh8+DC6du0Kf39/3Lt3D99++y1Onz6N7du3a69aIiKiOkajoL548SI2bNiArVu3AgDGjRuHpUuXomPHjmKb8ePHo1WrVtqrlIiIqA7SKKjd3Nzg5eWFlStXYuTIkZDL5RXa2NnZ4Y033tBGjURERHWWRueo582bh2PHjsHX17dCSK9Zs0b8/99++63mFRIREdVhGgX13Llzla4LCAioST1ERET0DK1+hWhKSgoaNGigzSGJiIjqNLXOUXfq1KnS/weAsrIyJCUlYdiwYdqrjoiIqI5TK6jHjRsHAIiLixP/v5yhoSEcHR0xZMgQ7VZIRERUh6kV1LNmzQIANGrUCO+8805t1URERERPaXSOmiFNRET0fKi8R13+ZSbx8fEKX2xSmfj4+JpXRkRERKoHtb+/f6X/T0RERLVH5aD+4IMPKv1/IiIiqj0anaPOzs7Ghg0bxMfbt29Hly5dMGbMGDx48ECb9REREdVpGgX1J598AplMBgDIzMyEv78/+vTpg3v37olXhhMREVHNaXRTjj179mDBggUAgP3796Nnz55YsmQJUlJS0LVrV23XSEREVGdptEddVFQEQRAAAMeOHcOAAQMAAJaWligsLNRuhURERHWYRkHt7u6OwMBAbNmyBTt37sTQoUMBAOfPn+ceNRERkRZpFNQrV65EYmIi5syZg7lz56Jdu3YAgEWLFuGzzz7Tdo1ERER1lkbnqNu0aYOIiIgKyw8cOKCNmoiIiOgprd7mkoiIiLRLo6BOT0+Hn58fWrZsCUtLywo/REREpB0aHfqeNGkS7t+/j9mzZ8PKykr7VRERERGgaVBHRETg6tWraNasmfYrIiIiIpFGh76tra0hl8u1Xw0REREp0Ciox48fj6+++gqlpaXar4iIiIhEGh36Dg8Px4ULFxAcHAxHR0fxe7/LRUdHa6s+IiKiOk2joB4xYgRGjBih/WqIiIhIgUZBPWfOHO1XQkRERBXU6AtPSktLkZKSor1qiIiISIFGQV1YWIipU6fCxMQEDg4O4nI/Pz/89ddf2qyPiIioTtMoqD///HPExcXh0KFDCsvfeOMNfPnll2qNlZmZiRUrVmDOnDn4/fffVbqSXJ0+58+fx6xZs7Bnzx616iIiIpICjYI6JCQEW7ZsQa9evRSW9+jRA4cPH1Z5nFu3bsHFxQW7d+9GaWkpZs+ejcGDB1cZvOr0ycnJgY+PD9avX49jx46pOUsiIiLd0yio7927h6ZNmwKAwkezHj9+rNZnqz/55BO0atUKhw8fxoIFC3Ds2DEcP34cwcHBWukzefJkjBs3Ds2bN1d7jkRERFKgUVB37NgRR44cAf4R1GvWrEHnzp1VGuPJkyfYu3cvxo4dCz29/y/D0dERvXv3xu7du2vcZ82aNfj777/x+eefazJFIiIiSdDo41nz5s3D2LFjMX36dADAqlWrcODAAezfvx8HDx5UaYw7d+7g0aNHcHJyUlju5OSEqKioGvWJj4/H3LlzERUVBX19fZXqKS4uRnFxsfg4Ly9PpX5ERES1SaM96iFDhiA0NBRRUVEwNTVFUFAQ8vLycPDgQfTr10+lMQoLCwEA5ubmCsstLCzw8OFDjfsUFhZizJgxWLx4cYVAr8qCBQtgYWEh/tjb26vcl4iIqLZotEcNAP369VM5lCtjamoKPL3g61k5OTkwMzPTuM9vv/2G1NRUxMfHY9asWQCAtLQ0REREYNasWVi8eLF42PxZs2fPxsyZM8XHeXl5DGsiItK5Gn3hSU04ODjAxMQE165dU1h+9epVtG/fXuM+3bp1w5w5c2BjYyP+GBoawsTEBDY2NhW+l7ycXC6Hubm5wg8REZGuqbxH3apVK5UHvXnzZrVt9PX1MXLkSGzevBlTp06FXC5HfHw8Tp8+jZ07d4rtQkJCkJycjMDAQJX6dOrUCZ06dVJ4rt9++w1ubm7iHjYREdGLQuWgDgoKEv//1q1bWLp0KUaNGoWuXbsCAM6dO4fQ0FC1wnDRokV45ZVX4O7uDjc3N4SHh8PX1xevv/662ObgwYOIjo5GYGCgyn2IiIj+LVQOan9/f/H/Bw4ciA0bNsDPz0+hzZYtW7B161aVn9zW1haXLl1CeHg4MjIyMH78eHh5eSm0efPNN9GzZ0+1+vzTzJkzeb6ZiIheSDJBEAR1O1laWiI5ObnCedy8vDw0b94c2dnZ2qxRJ/Ly8mBhYYHc3Nwan68ef+gHpeu2vPp+jcYmIqIXjzoZo9HFZEZGRuIXnjzryJEjMDIy0mRIIiIiqoRGH8+aNWsWxo0bh3feeQddu3aFIAg4f/48Nm3ahG+//Vb7VRIREdVRGgV1YGAgnJycsHz5coSEhAAAOnTogO3bt2P48OHarpGIiKjO0vgLT0aMGIERI0ZotxoiIiJSoLMvPCEiIqLqMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBjUREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBjUREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBjUREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBjUREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIMdF3A48ePceLECWRkZMDFxQWurq5a6ZOamoqzZ8/CwMAAXbp0ga2tbS3NgIiIqPbodI/6wYMH6Nq1KwICAhAaGopevXph+vTpNeojCAL8/PzQo0cP/Prrr1i1ahVatWqFFStWPIcZERERaZdO96hnz56Nx48fIy4uDiYmJjh37hzc3d0xZMgQDBo0SKM+giBg4MCB2LRpE/T19QEAmzZtwqRJkzBkyBC0bNnyOc+SiIhIczrboy4rK8Pvv/+OiRMnwsTEBADQtWtXeHh4YPv27Rr30dPTw7hx48SQBoDBgwejrKwMV69efS5zIyIi0had7VGnpKQgLy8PHTp0UFju7OyMCxcuaK0PAOzbtw96enpwcXFR2qa4uBjFxcXi47y8PDVmQ0REVDt0tkddHoRWVlYKyxs0aKA0JDXpc/36dcycORMffPAB7O3tldazYMECWFhYiD9VtSUiInpedLZHbWRkBADIz89XWJ6fny+uq2mfpKQk9O/fH/3798fixYurrGf27NmYOXOm+DgvL++5hPWBAwcqXe7t7V3rz01ERNKns6B2cHCAoaEhkpKSFJbfvn0brVq1qnGfO3fuwMvLC926dcO2bdsUzllXRi6XQy6XazwfTW3Tu1bpcm8wqImISIeHvuvVq4eBAwciODgYgiAAANLT03H8+HEMGzZMbHf69Gns2rVLrT7Jycnw8vJCly5dEBwcDAMDnX9cnIiISCMyoTzxdCAhIQHdu3dH37594eHhgc2bN8PS0hInTpyAoaEhAMDf3x/R0dGIj49Xqc+jR4/QsWNH5OfnY968eQoh/corr6Bt27Yq1ZaXlwcLCwvk5ubC3Ny8RvMcf+gHtftsefX9Gj0nERFJlzoZo9NdzQ4dOuDSpUvYtGkT7ty5g//85z+YMGGCGNJ4Gq7NmjVTuc+TJ0/g5eUFADh//rzC87Vr107loCYiIpICne5RSxn3qImIqLaokzG8KQcREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBjUREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBjUREZGEMaiJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIMdF0AVW78oR8qXb7l1fefey1ERKQ73KMmIiKSMAY1ERGRhDGoiYiIJIxBTUREJGEMaiIiIgnjVd8vmM+u7Kp0+TftRz73WoiIqPZxj5qIiEjCGNREREQSxqAmIiKSMJ6j/pc4cOCA0nXe3t7PtRYiItIe7lETERFJGIOaiIhIwnjo+wWTkpKiZE3b51wJERE9D9yjJiIikjDuUf9LbNO7pnSdN3gxGRHRi4p71ERERBLGoCYiIpIwHvquA8Yf+qHS5W+VVX4BGj93TUQkHQzqOkzZeW2e0yYikg5JBHVKSgoyMjLQpk0bmJuba62PJuMSv+WMiEhKdHqO+tGjRxg1ahTatm0LPz8/2NjYYMWKFTXuo8m49D/b9K4p/SEioudLp3vU8+bNw9mzZ3Hr1i3Y2tpi9+7deP3119GtWze4u7tr3EeTcUk1yu6HrQzvk01EVDM6DeqNGzciICAAtra2AIARI0agY8eO2Lhxo9JAVaWPJuNS7agq2NUNcW2ORUT0otBZUKelpSEjIwOdO3dWWN6tWzfExsZq3EeTcQGguLgYxcXF4uPc3FwAQF5engazU1Ty8FGNx5CKW1dvqNXezq6Z0nWzzv2mhYr+n7LXaf61PZUun9v2Na09NxGRusrfswRBqLatzoI6KysLANCwYUOF5Q0bNhTXadJHk3EBYMGCBZg3b16F5fb29irOiHRpWS23JyKqDfn5+bCwsKiyjc6C2tDQEHh64dezioqKUK9ePY37aDIuAMyePRszZ84UH5eVlSErKwsNGzaETCZTc3b/k5eXB3t7e6SkpPDKcy3jtq0d3K61g9u19ryI21YQBOTn56Np06bVttVZUNvb20NPTw+pqakKy1NTU+Hg4KBxH03GBQC5XA65XK6wzNLSUu15KWNubv7C/AN60XDb1g5u19rB7Vp7XrRtW92edDmdfTzL2NgY3bt3x549/zuH+PDhQxw5cgQDBgwQl928eVM8t6xKH1XHJSIiehHIBFXOZNeSkydPYsCAAQgMDISnpydWrFiB27dv4+LFizA1NQUA+Pv7Izo6GvHx8Sr3UaXN85KXlwcLCwvk5ua+UH/pvQi4bWsHt2vt4HatPf/2bavTLzzp3bs3jh8/jjt37uCHH36As7MzTp8+rRCmrVu3hpubm1p9VGnzvMjlcnzxxRcVDqtTzXHb1g5u19rB7Vp7/u3bVqd71ERERFQ13uaSiIhIwhjUREREEsagJiIikjAGdS27fv06YmJiFL6elKqWl5eHyMhIpKWlKW1z//59nDt3DhkZGTVqU5dkZ2cjNjYW9+/fV9qmpKQEsbGxuHr1ao3a1CVlZWW4du0aEhISKnzR0rNu3ryJCxcuoKioqEZt6pqysjJERkYiLi6u0vVZWVk4d+4c0tPTlY6hShtJE6hWpKenC926dRMsLS2Fli1bCg0aNBDCw8N1XZakJSUlCVOmTBFsbGwEAwMDYfny5ZW2++STTwS5XC506NBBkMvlwnvvvSeUlZWp3aauiI+PFwYNGiQ0aNBAePnllwUTExNh1KhRQn5+vkK7I0eOCNbW1kKLFi2Ehg0bCq6urkJycrLabeqSdevWCQ4ODkKHDh2EVq1aCZaWlsJPP/2k0ObBgwdCz549BXNzc6F169aChYWFEBISonabuuqbb74R9PT0hM6dO1dYN2/ePIXf83feeUd48uSJ2m2kjkFdS4YOHSq4u7sLhYWFgiAIwldffSWYmpoKGRkZui5Nsg4fPiysXr1ayM/PFxo2bFhpUAcHBwtyuVw4e/asIAiC8NdffwkmJibC2rVr1WpTl/zxxx/Cvn37xMdpaWmCo6OjEBAQIC7Lzs4WrKyshE8//VQQBEEoLi4WXnnlFaFPnz5qtalrli9fLmRmZoqP169fLwAQrl+/Li7z8fERXF1dxT+Mli9fLsjlcuHOnTtqtamLzpw5Izg6OgpvvfVWhaAODw8XDAwMhFOnTgmCIAjXr18XLC0thWXLlqnV5kXAoK4Fd+/eFWQymbBz505xWWFhoWBiYiKsWLFCp7W9KJQF9auvviqMGDFCYdm4ceMEd3d3tdrUdR988IHg4uIiPt6wYYNQr149ITc3V1wWHh4uABBu376tcpu6LjExUQAgnDhxQhAEQcjNzRUMDQ2FTZs2iW0eP34sNGjQQFiwYIHKbeqi7OxsoUWLFsKxY8eEadOmVQjqkSNHCv3791dYNnXqVMHZ2VmtNi8CnqOuBXFxcRAEQeFWm0ZGRnB2dq7yVptUvdjY2EpvYXrx4kXxdnGqtKnrzp8/j1atWomPY2Nj0bp1a4VvderWrZu4TtU2ddHdu3dx+vRp7N69GxMnTsTgwYPRs2dPAEB8fDweP36s8O/RwMAAnTp1EreZKm3qIn9/f7zxxhvo06dPpeuV/Z5fuXJFvCZIlTYvAp3dlOPfTNNbbVL1yu9o9qyGDRuiuLgYhYWFMDExUalNXbZy5UpER0cjMjJSXFbZNmvQoIG4TtU2ddGff/6JJUuWIC0tDSUlJVi9ejX09fWBWr4177/Z6tWrkZiYiG3btilto+z3vKysDDk5OWjSpIlKbV4E3KOuBZreapOqZ2hoWOl2BaBwq9Pq2tRVO3bswMyZM7Fu3TpxbxhKtln546q26z/b1EXDhw/H6dOnkZiYiCVLlmD48OGIiooCavnWvP9WqampCAwMREBAAM6ePYvTp08jPT0dBQUFOH36NHJzc4E69l7AoK4FzZs3B57+g3tWdbfapOo1b9680u1qY2MjvuGp0qYu2rlzJ/z8/LBmzRq8/fbbCuuUbTMA4r9ZVdrUdb6+vmjWrBkOHjwIqPhewPcLRUVFRXBzc8PmzZsRFBSEoKAgREVFITU1FUFBQUhMTASq+PdoZmYGKysrldu8CBjUtaBTp05o1KiRwq02r1y5ghs3bvBWmzU0YMAA7N27F2VlZeKyPXv2KGxXVdrUNaGhoRg7dix++uknTJw4scL6AQMGIDU1FTExMeKysLAwmJubw93dXeU2dUlxcTFKS0sVluXk5OD+/fvi4da2bdvC3t5e4b3gzp07uHjxovjvUZU2dUmrVq1w+vRphZ+RI0eibdu2OH36NF5++WXg6b/H/fv348mTJ2LfsLAw9O/fX3ysSpsXgq6vZvu3Wrt2rSCXy4Xvv/9eCAkJEZydnYW+ffvquixJKygoECIiIoSIiAjBwsJCmDFjhhARESEkJCSIbZKTk4WGDRsKvr6+wp49e4SJEycKZmZmwtWrV9VqU5ccOHBAMDQ0FCZNmiRu34iICCE6Olqh3bBhw4Q2bdoIv//+u7Bq1SrByMhI+O6779RuU1fcvHlT6NKli/DTTz8Jhw4dErZs2SJ07txZcHJyErKzs8V2v/32m2BgYCAsXrxYCA0NFdzc3AQPDw+htLRUrTZ1WWVXfd+9e1ewsbERRo4cKezZs0cICAgQjIyMhIsXL6rV5kXAu2fVol27duHXX39FQUEBXnnlFQQGBtb5C5mqcvPmTbzzzjsVlvfq1QvffvutQrvFixfj1q1baN68OQIDA+Hs7FxhrOra1BWrVq3C9u3bKyy3tLREeHi4+PjRo0dYvnw5jh8/DiMjI/j4+MDX11ehjypt6pIbN27gp59+QkJCAqysrODp6Ql/f/8Kv+fh4eHYuHEjcnNz4enpiY8++qjCfZNVaVNXLV++HJcvX8a6desUlt+5cweLFi3CtWvXYGdnhw8++EDc41anjdQxqImIiCSM56iJiIgkjEFNREQkYQxqIiIiCWNQExERSRiDmoiISMIY1ERERBLGoCYiIpIwBrVECYKAyMhI7NixAwkJCbouR6eOHj2K+Ph4XZdBdfy1yM3NxcGDBxEcHCzeGEJqSkpKEBwcjJycHF2XIgnJycn4448/dF1GjTGoJWrYsGGYMGECQkNDce3aNa2Ne+PGDezbt09r42nbiRMncOnSJYVlX3zxBXbu3FkrY/8baWuetflavGjS0tLQunVrLFiwALt370Z+fr6uS8KTJ08QHByMBw8eiMvy8vLg6+uLpKQkndamqcrmVBNnzpzB5MmTtTKWLvF+1BKUmpqKvXv34vr162jdurVWxz548CBWrlyJwYMHa3Vcbfn666/RpUsXvPTSSy/U2FKirXnWle2lij/++APW1tY4ceKErksRPXr0CL6+voiKiqpwz+UX1b9xTtrAoJaYxMREhISEAE/3aC5cuIChQ4fC1NQUAHD79m1cunQJjRo1gpubG4yMjMS+qampiIiIAAAYGxujXbt2aNOmjbg+OTkZMTExyM/PR3BwMADAzc0Njx49QlZWFry8vMS2t2/fxl9//YXXXnsNAHDt2jUkJSWhT58+iIqKQkZGBt544w3o6emhrKwM58+fR1paGpycnODi4lJhXlevXsX169dhb28PV1dX6OlVPJhTPu6VK1fE+oYNGyauz8nJQWxsLADAw8NDYe5//vknbt++DZlMhoYNG4p3MKtu7Mq+e72qWqub69GjR9GkSRPY2tri4sWLEAQBvXr1Urj3bXW1PjuOjY0N/vzzT5iamqJ37941mqcqr5M2XgtVtlNlcnNzERERAUtLS7z88su4ePEijIyM4Obmpvb2t7Oz01p9EREROHHiBEpLSxEcHAxLS0t4e3srfY0A4N69e/jzzz9hYGCA7t27w8LCQhyvpKQEu3btgre3N7KyspCQkABbW1t07twZAHD58mXcvHkT7dq1Q9u2bZXWVX5I98iRI0hKSoKVlZU4BgAkJSXh8uXLsLW1VdiG5QoKChAVFYUnT57A1dUVTZs2reLV+X+FhYWIjo5GSUkJPDw8YGlpqbBe1Xnn5ORUWltlc+rTp4/Y7+7du0hISICrqyucnJzE7XX9+nU0btwYHh4e0NfXr3IOWVlZuHDhAvT09NClSxeFGqWKQS0xd+7cwfHjxwEA+/fvR7169dCnTx8YGxtj2rRpCAkJgYeHBzIyMnD//n388ccf6NSpEwAgPT0du3fvBp7+EkZGRuL111/Hhg0bgKdBHhcXh/z8fLGdhYUFTp48ifPnzysE9cmTJzFnzhwxqP/73//ihx9+gLW1NczNzWFjY4PXX38d6enpeO2111BQUIB27dohNjYWzs7O2LVrl/jmOGnSJISFhaFHjx5IT09H/fr1ERYWVuF+sOfOnUNmZiZKS0vF+vr16wcAOHbsGH777Te0a9dOPBXw559/okGDBgCAmJgYnDx5UtwOMTEx+Pnnn8UbRigb+59BXVWtqamp1c71iy++gEwmw+3bt+Hi4oKEhASYmpri+PHjaNy4sUq1PjtOSkoKXFxc0KNHD/Tu3Vvjeebk5FRbu7ZeC1W20z9duHABAwcOhI2NDWxtbXH79m0YGRmhd+/e4hu5qttfT08PaWlpWqsvKioKV65cwYMHD7B7927Y29vD29tb6Wu0YcMGTJ8+HV27dkVRURGuX7+O7du3Y9CgQcAzh6d79+6Ne/fuoUWLFjh27BjefvttPHnyBNHR0XBwcMCxY8ewfPlyBAQEVLrN9u/fDzz9XY2Pj4ejo6MY1HPnzhWPyEVERGDkyJHi+wAA7Nu3D+PHj0f79u1hZmaGM2fO4LPPPsNHH31U6XMBwN69ezF+/Hg0bdoUzZo1w82bN7F+/XrxjxNV5z1s2DCltSmbk6+vL4YMGYLr16/D1dUV5ubmaNGiBd566y0cOHAAnp6euHz5MqysrHDw4EHY2Ngo3WY+Pj7o1KkTjIyMcOPGDaxcuVKsUbJ0ffsuqig2NlYAIGRmZorLfvzxR6Fdu3YKt8+bM2eO4OrqqnSc9PR0wdraWggPDxeXrVixQmjbtq1Cu08++UTo16+fwrKNGzcKzZo1Ex8vWbJEACBs375doV3fvn2FgIAAoaysTBAEQSgqKhK6dOkifPHFF4Lw9FaAAITExESxT3R0tJCamlppzf369RM++eQThWU9evQQmjZtKmRkZAiCIAjFxcWCk5OT8O233yqde2hoqGBpaSkUFBRUOfazqqu1urmW12pkZCRcu3ZNEARBePjwoeDm5iZMmTJFrVp79OghNGjQQOl2UneeqtT+T5q+Fpo8l7u7u+Dr6yv22b9/vwBAmDZtmlrj1lZ9X3zxhdCjR48K2+Kfr1FKSopgZGQkrFu3TlwWFBQk2Nraiq9RZmamAEDw9fUVb2P566+/CgCECRMmiHWtXr1aaNCggdKa8vPzBQBCVFSUuKx87DFjxohjnz9/XgAg3i42LS1NMDU1Ffbt2yf2i4+PF4yMjISYmJhKn6t8Xl9//bW47MGDB0JERITa866qtqrmNHToUOHx48fi8nXr1gnm5ubCrVu3BOHpbXI7d+4s+Pn5iW22b98uNGzYUHzcs2dPYe7cueLjnJwc4ejRo0q3sVTwYrIXxMaNG+Hq6oojR44gJCQEO3bsgIWFBeLi4pCdnS22KykpwZkzZxAaGooTJ07Azs4OZ8+e1UoN1tbW8PHxER8nJyfj2LFjcHJyQmhoKEJCQrBnzx60bNlSPCogl8shk8kQFxcn9nN3d1fpMNuzRo0aJe6R1qtXD56enhUussvMzMTRo0fx+++/o6ioCDk5Obh+/brKz1FVrarMtdzw4cPFUw7GxsaYPn06duzYoXato0aNqnQ7qTtPdWpXRVWvhSbPlZKSgj///BOBgYGQyWQAAG9vb3Ts2FGjOWi7vuq2xbOvUVhYGCwsLDBx4kRx2aeffoqMjAzxSEi5yZMni6dVPD09xWXl28DT0xNZWVm4f/++2nVNmTJFHLtz584wNTUVt0FISAiMjY3x8OFDhISEICQkRDwMrewcfEhICExMTBAUFCQua9CgAXr27Kn2vKuqrSoBAQEwMPjfQeDg4GD4+PigZcuWAAATExN8+OGH2LFjB8rKyiodo3wvurCwEHh6RLFv377VPreu8dD3CyIpKQkGBgYVrrgdM2YMioqKYGVlhQsXLuC1116DmZkZWrduDRMTE2RmZuLevXtaqcHW1rZCTQAQGRmJc+fOictlMhm6dOkCALCzs8NPP/2EgIAA/Oc//0GfPn3w9ttvi4dRVVV+2LKcXC5HQUGB+Pj777/HnDlz4OrqChsbGxgaGkImk6k196pqVWWu5RwdHRUet2jRAjk5OcjLy4O5ubnKtf5ze2s6T3VqV0VVr4Umz5WSkgJUst2efazOuNquryr/fI3u3LkDR0dHMWwBwMzMDNbW1rhz545C22dP/cjlcqXLHj16pHZdlW2D8nGSkpJQVlZW4b2ka9euSg8ZJycno2XLlkrP/6oz76pqq0pl27p///4Ky5ycnFBcXIyMjIxKf3+WL1+OyZMno3HjxvD09MSwYcMwZcoU1K9fv9rn1yUG9QvC3Nwc3t7emD9/vtI2s2fPxtChQ7F27VpxmZeXF6q75Xj5BWHPquwX59lfwvKaAGDOnDmVXqxSburUqXj33Xfx119/ISwsDIMGDcLOnTvF8981lZ+fj8DAQISHh4vnmnJycvD7779XO3dVa3VwcABUmCsAhSMc5Y/lcjlMTU3VqvWf21vTear6OmmDJs9V/sadk5OjcKVvdnY2mjdvrvG42qqvKv98jRo1aoSsrCyFZYIgICcnp8IFg7pibm4OIyMj8SJBVVhaWlb5kannMW9VtnVWVhb09PQq/DFQztnZGWfOnEFGRgaOHTuG+fPn48SJE9i1a5dWaqwtPPT9gvD29saWLVvEQzblUlNTxf+/e/euwlWiKSkpFQ57m5qaVgjhZs2aITExUeHNXpXDgC4uLmjWrBnWrFlTYV1aWhrw9BenqKgIMpkML730Ej7//HN07doV0dHRlY5ZWX3VyczMRFlZmcLcK/usb3VjV1WrKnMtt2/fPpSUlIiPd+3aBU9PT+jp6alcqzbnqU7tVY2jCk2ey8nJCTY2NggLCxOXpaenK+zxajoHbdSnjp49e+LGjRv466+/xGV79uxBWVkZunXrVuPxn2VkZAR9fX21XyNvb2+kpKRg7969CssLCwuVflHKq6++ilu3biEyMlJheWZmJqDFeaszp549eyIsLAylpaXispCQEHTp0kU8GvFP5e+XTZo0ga+vL2bOnKn0vUhKuEf9gvjqq69w4sQJdOvWDRMnTkT9+vURFRWFu3fv4vDhwwCAESNGYOHCheL5nx9//LHCP9jOnTsjJSUFCxcuhKOjI9zc3DBy5EgEBQXBz88Pffv2xcmTJ3HixAkYGhpWWZO+vj7Wr1+P119/HQ8ePMDAgQORlZWF8PBwvPHGG/jggw9w584djBkzBqNHj0br1q1x+fJlxMbGYtmyZZWO2aVLF6xbtw4uLi4wMTFR+EiQMo6OjnB1dYWfnx8mTpyIGzduYOPGjRX+Aq9s7Gev+q6qVlXmWq6oqAgDBgzAW2+9hfPnz2PHjh3iHz6q1qrteapae01fC3W2UzlDQ0N89dVXmDZtGrKysmBra4vVq1fD1NRUnJsm42qrPnV0794dPj4+GDx4MAIDA1FUVIRFixbho48+Eo8OaIu+vj46deqEZcuWITU1FY0aNVL4eJYyHh4emDlzJkaPHo3p06ejbdu2uHnzJkJDQ/HHH39U+MgVnp4vf++99zBkyBDMmDEDzZo1w969ezF06FBMmTJFa/NWZ05BQUHYvn07Xn31VYwZMwbnzp3D9u3bceTIEaXjl5/T9vT0RElJCb7//nuMHj1a5fp0hUEtQVZWVhgzZoxCyFpbWyMmJga//vorzp49C2NjYwwePBhvvvmm2ObLL7+Eo6MjIiMjUb9+faxZswYJCQkwMzMT27i4uCAsLAz79u1DXFwcLCwsMGjQIJw9exbr16/HuXPn0K9fP0yYMEH8PDcAtGvXDt7e3hVqHThwIC5fvowtW7YgMjISdnZ2WLZsGdzd3QEAL7/8Mo4fP47Nmzfj5MmTsLW1xblz5+Ds7Fzp3GfNmgUzMzNERkaisLAQ/fr1Q//+/RUuLAKAbt26iX916+np4ejRo1i5ciVOnTqFZs2a4cyZM5g7d67CearKxn42qKurtbq5lps2bRo6dOiAU6dOoV69eoiMjBTPf6paa2Vzrsk8Va29pq+FOtvpWf7+/rCxscHu3btRVFSElStXYt68eQr/dlUZt7bq++eYyp4LAH799Vds3boVJ0+ehIGBATZt2oQRI0aI6+VyOcaMGaNwPtrY2BhjxowRD83j6SHqMWPGwNjYWGldoaGhWLNmDfbv3w87Ozt07969wtgAMHLkSIXAXLZsGYYMGYI9e/YgMjIS7dq1w6lTp5SeowaAVatWYdCgQdi3bx/u37+Pd999F0OGDKnRvCurTdU5NWjQQPx44unTp9GkSROcP39e4TVp3rw5Ro4cKT4+fvw4tm/fjjNnzkBfXx+LFi3CqFGjlM5ZKmSCuifxiEipnj17on///vjyyy91XcoLJTs7W+GNODMzE05OTti4ceML8UZKVJu4R01EOhcbG4tvv/0WI0aMwOPHj7FmzRp06NBBpcPtRP92DGoiLVJ2OJSq1rdvX+jr6yMsLAxFRUUIDAzEhAkTqr1Ogqgu4KFvIiIiCePHs4iIiCSMQU1ERCRhDGoiIiIJY1ATERFJGIOaiIhIwhjUREREEsagJiIikjAGNRERkYQxqImIiCTs/wBaurHxlUqa2QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x320 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.linspace(0, score.max(), 60)\n",
    "is_control = np.isin(genes, list(controls))\n",
    "plt.figure(figsize=(5, 3.2))\n",
    "plt.hist(\n",
    "    score[is_control], bins=bins, density=True, alpha=0.6, color=\"#888\", label=f\"control genes (n={is_control.sum()})\"\n",
    ")\n",
    "plt.hist(\n",
    "    score[~is_control], bins=bins, density=True, alpha=0.6, color=\"#2b7\", label=f\"other genes (n={(~is_control).sum()})\"\n",
    ")\n",
    "plt.xlabel(\"features that separate the gene from the controls\")\n",
    "plt.ylabel(\"density\")\n",
    "plt.legend(frameon=False, fontsize=8)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "scores-read",
   "metadata": {},
   "source": [
    "Most genes sit near zero whether or not they are expressed: a typical knockout separates on a handful of features, and so does a typical non-expressed gene. The hits are the long tail on the right, and the control genes almost never reach it. That gap is what the calibration turns into a rate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdr-intro",
   "metadata": {},
   "source": [
    "## Calibrating the hit list\n",
    "\n",
    "{func}`~mantispy.tl.empirical_fdr` takes the per-gene score and the control set and writes two numbers per gene. It wants one row per gene, so the score goes on a small object keyed by gene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fdr",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:09.971116Z",
     "iopub.status.busy": "2026-10-21T00:50:09.970915Z",
     "iopub.status.idle": "2026-10-21T00:50:10.477509Z",
     "shell.execute_reply": "2026-10-21T00:50:10.476400Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hits at p<=0.01: 331\n",
      "hits at q<=0.01: 102\n"
     ]
    }
   ],
   "source": [
    "import anndata as ad\n",
    "\n",
    "per_gene = ad.AnnData(score.reshape(-1, 1), obs=pd.DataFrame({\"gene\": genes, \"score\": score}, index=genes))\n",
    "mt.tl.empirical_fdr(per_gene, control_genes=controls, score=\"score\", group=\"gene\", alpha=0.01)\n",
    "mt.tl.empirical_fdr(\n",
    "    per_gene, control_genes=controls, score=\"score\", group=\"gene\", alpha=0.01, criterion=\"q\", key_added=\"efdr_q\"\n",
    ")\n",
    "\n",
    "p = per_gene.obs[\"empirical_fdr_pvalue\"]\n",
    "q = per_gene.obs[\"empirical_fdr_qvalue\"]\n",
    "print(f\"hits at p<=0.01: {int((p <= 0.01).sum())}\")\n",
    "print(f\"hits at q<=0.01: {int((q <= 0.01).sum())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "hit-list",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:10.480587Z",
     "iopub.status.busy": "2026-10-21T00:50:10.480280Z",
     "iopub.status.idle": "2026-10-21T00:50:10.536323Z",
     "shell.execute_reply": "2026-10-21T00:50:10.534696Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          score pvalue qvalue\n",
      "ITGAV  431.0000 0.0018 0.0018\n",
      "MDM2   422.0000 0.0018 0.0018\n",
      "PSMD7  384.0000 0.0018 0.0018\n",
      "NXF1   373.0000 0.0018 0.0018\n",
      "POLR2A 372.0000 0.0018 0.0018\n",
      "HSPA5  367.0000 0.0018 0.0018\n",
      "LRR1   360.0000 0.0018 0.0018\n",
      "KIF23  344.0000 0.0018 0.0018\n",
      "TSG101 342.0000 0.0018 0.0018\n",
      "SUPT6H 331.0000 0.0018 0.0018\n",
      "RAN    330.0000 0.0018 0.0018\n",
      "EIF4A3 330.0000 0.0018 0.0018\n"
     ]
    }
   ],
   "source": [
    "hits = per_gene.obs.loc[per_gene.obs[\"empirical_fdr\"], [\"score\", \"empirical_fdr_pvalue\", \"empirical_fdr_qvalue\"]]\n",
    "print(\n",
    "    hits.sort_values(\"score\", ascending=False)\n",
    "    .head(12)\n",
    "    .to_string(header=[\"score\", \"pvalue\", \"qvalue\"], float_format=lambda v: f\"{v:.4f}\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hit-list-read",
   "metadata": {},
   "source": [
    "The two top scores are ITGAV, an integrin, and MDM2, a regulator of p53. Most of the rest of the top twelve are machinery a cell cannot run without: RNA polymerase (POLR2A), the proteasome (PSMD7), mRNA export and ribosome biogenesis (NXF1, EIF4A3, RAN), the chaperone HSPA5, and cell-division genes (KIF23, TSG101). These are the genes to follow up."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "two-numbers",
   "metadata": {},
   "source": [
    "Two numbers come out, and they answer different questions.\n",
    "\n",
    "The p-value is the share of control genes that reach a gene's score, with one added to the count and the control total so a gene past every control does not get an exact zero. A cutoff on it fixes the rate at which a control gene is called, the way PERISCOPE sets its threshold, and here it calls 331 genes.\n",
    "\n",
    "The q-value is a target-decoy false discovery rate. At a gene's score it weighs the share of control genes reaching it against the share of all tested genes reaching it, so it estimates how much of the hit list is false. It accounts for the thousands of genes tested, which the p-value does not, so it is stricter: 102 genes at `q <= 0.01`. It treats every tested gene as a potential null, so it is a conservative bound rather than an exact rate. Quote the q-value as the honest rate, and read the p-value as the PERISCOPE-style cutoff. Every q-hit is also a p-hit.\n",
    "\n",
    "The `empirical_fdr` call above used `criterion=\"p\"` by default, so `obs[\"empirical_fdr\"]` holds the p-value hit list. Pass `criterion=\"q\"` for the honest list, or read `obs[\"empirical_fdr_qvalue\"]` directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "hits-fig",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:10.539590Z",
     "iopub.status.busy": "2026-10-21T00:50:10.539259Z",
     "iopub.status.idle": "2026-10-21T00:50:10.631360Z",
     "shell.execute_reply": "2026-10-21T00:50:10.630394Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAE2CAYAAACukdFuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAT81JREFUeJzt3XdYFFfbBvB7aUsHBURUqo2gKGKvqNg7SYyiscUejfqaoqap0URN9IsxMXaNscSY2LHFKEYxsYOx+4qgdOm97p7vD3VfV9qyLOwC9++6uC535pyZZwbh4Zw5c45ECCFAREREOklP2wEQERFR8ZioiYiIdBgTNRERkQ5joiYiItJhTNREREQ6jImaiIhIhzFRExER6TAmaiIiIh3GRE1ERKTDDLQdgEwmw/Hjx3Hv3j2MGDECjo6OhcpkZWXh6NGjiIuLg6enJ3x8fCqsDBERkS6RaHMK0d9//x0ffvghHB0dcf78eQQGBqJ79+5KZWJiYtCtWzeYmZmhVatWOHr0KPr06YOdO3dqvAwREZGu0WqL2s7ODmfPnoW+vn6RLWkAmDdvHiwsLPDPP/9AKpXi1q1baNmyJd58800MGzZMo2WIiIh0jVafUfv4+MDZ2bnY/TKZDPv378e4ceMglUoBAM2bN0eXLl2wd+9ejZYhIiLSRVp/Rl2SiIgIZGZmomnTpkrb3d3dcfnyZY2WKUpubi5yc3MVn+VyOZKSkmBjYwOJRKKRayQioppHCIH09HTUq1cPenolt5l1OlGnp6cDAKytrZW2W1tbK/ZpqkxRli1bhsWLF2vkWoiIiF4VERGBBg0alFhGpxO1qakpACAtLU1pe2pqKszMzDRapigLFizA3Llzlco7OTkhIiIClpaW5bw6IiKqqdLS0uDo6AgLC4tSy+p0onZ2doZUKkVoaKjS9tDQUDRu3FijZYoilUoVz7RfZmlpyURNRETlpspjVJ2e8MTAwAADBw7Erl27IJfLAQDh4eH466+/4Ofnp9EyREREukir71Hfvn0bx48fR2pqKpYuXYpp06ahYcOG6NSpEzp16gQ8b/V26tQJnp6eaNeuHfbs2YPGjRvj2LFj0NfX12iZ0qSlpcHKygqpqalsURMRkdrKkk+0mqivXr2KPXv2FNrep08f9OnTR/E5Pj4ev/zyi2JGseHDhxdKrpoqUxImaiIi0oQqk6irGiZqIiLShLLkE51+Rk1ERFTTMVETERHpMCZqIiIiHabT71FTzXXx4kWsXr26yMGGZZWVlYU333wThw8fhoGBbv2XX716NQBgzpw52g6lUm3duhVCCEycOFHboRApJGUWYNLPj/HnvbQSy73eqhZ+Gu9SaXGxRU1quXjxIkaOHKn2/tLk5OQgNjZW7fov+/777+Hj41MhSfrChQt4++231a6fkpKClJQUjcakLatWrcIPP/ygUtm33noLy5cvR15eXoXHRaSKe7E5aL/8Hg6EpCA9R17iV3aevFJj063mBVUZXl5eitagOvsr08aNG3H27NkKOXZ2drbG/qCo6pKTk2FsbKxSWXNzc7Rt2xaHDh3C8OHDKzw2opKcuJWKkZvDkJotg7ONEXZMcEE9a6Niy5sZVW4bly1qDRNCIDNXVqFfqrxRd/ToUfj6+qJ58+aYNGkSEhISAAB//fUXRo8ejQ8++AAtW7aEv78/oqOjMWHCBHh4eGDWrFmK2dv++usvjBkzBvPnz4eXlxf8/f0RHx8PAAgJCVF0174o98knn6Bly5a4evWq0n4ACAgIQK9evdCwYUO4uLggOTkZd+7cgYuLC1xdXdGuXTuVW2MlxfWq0NBQFBQUKK13XlQsAPDTTz+hc+fOaNOmDb766ivIZDKl83322Wfw8vLCG2+8gadPnwIAxo0bhwsXLsDFxQXDhw8v8l7k5ubik08+QatWrdCtWzfs379fpevMzs7GrFmz0KJFC0yYMAGffvopNm3aBDxfunXVqlXo0KEDvL29sWTJEsX3bcWKFfjqq68wZswYeHp6Yu7cuYp9JdV71YoVK7B06VKMHj0aHh4ewPMuaxcXFzRs2BC+vr44c+YM8HxOhB9++AGrVq2Ci4sLVq1aVeq5unXrhpMnT6p0L4gqghACq/+Mw8AfHiI1W4Yujcxxeb47uja2QEM7abFfda0MKz1QUlFqaqoAIFJTU4stk5FTIDDlaoV+ZeQUlBjnP//8Izp16iSCgoLE3bt3xfz588XIkSOFEEIcP35cGBgYiE2bNolbt26J7t27C0dHR7F3714REhIiPD09xYEDBxRl9fT0xP/93/+JmzdvikmTJonBgwcLIYQIDAwUPj4+SuXWrl0r/vvf/4qcnByl/UFBQcLW1lbs3btXPHz4UISFhQmZTCZyc3NFWFiYCAsLExcvXhTe3t4iKCio0PFfVVJcRZXt2rWr4nNxsZw6dUrUq1dPnDx5Uly8eFF4enqKb7/9VumerVu3Tty+fVuMHj1azJ49WwghxM6dO0Xnzp1FWFiYiImJKfJefPzxx6Jr167i2rVr4sCBA8La2lpcu3ZNCCHEwoULxcKFC4uM/cMPPxS9e/cWwcHBYu/evcLMzEx88803QgghVqxYISZOnChu3LghgoODRf/+/cX69euFEELMmzdP1KtXTxw7dkxcv35d6XtaUr1XzZs3T9jb24ujR4+KsLAwIZ7/DISFhYnQ0FAREBAgnJycRGZmpsjJyREzZ84U77//vggLCxPJycmlnuvV7w1RZcrNl4lJP4crfq9O+ClM5OTJKu38quSTF9j1XQ398ssvuHfvHkaPHg08b0W9PPF7u3btMGnSJOD5s8KzZ88quh8HDx6M27dvY9iwYQAAT09P/Oc//wEArFmzBrVq1UJ+fn6hc3p7e+Pdd98tMp5ff/0Vs2bNKtTFKZfL8f333yMwMBDJyclISkrC9evX0blz51Kvsbi4DA2V/9ItKChQejZdXCwHDhzAzJkzFTPiLVmyBN9++62iV6Bdu3aYNm0aAGDixIn46quvAAD29vYwNjaGi8v/Bpa8ei8OHjyIzZs3w9vbG97e3vjrr79w+PBheHt7K8Xw4hh2dna4cuUKAgICsHv3bnh5ecHLywtHjhxRlN21axcSEhLw559/AgAyMzNha2uLqVOnAgDGjx+P/v37AwD8/PwU39PS6r1qzJgxGDBggOJzUlISPv74Y9y4cQOZmZmIjY1FaGgoPD09YWVlpXQvSjuXoaFhkf+XiCpafHo+3lj/COcfZkBPAqx8swHm+NZRaYEMbWCi1jBTIz1krPGq8HOUJDMzE9OmTcPkyZMV216eKvXlZdUMDQ0LfX75l+fLy4AaGxtDIpEU+cvVxsam2Hiys7OLnHnnu+++w71797BhwwbY2tpi0aJFyM3NLfHaSovr1UTt6Oio9Ay5uFhycnJgbm6u+GxhYYHs7Gylzy+UlmBevRelHfuFF8/RX3yviqr3QmZmJtavXw9PT88i70lx8ZZWr7RrmT59Ory8vPDxxx/D3Nwcffv2LfZ7Vtq5YmJi4OTkVOy5qeYqkAmsOBmLsw/SK+T4d2JyEJ2SD0tjPeyZ7Ib+za0q5DyawkStYRKJBGZS1ecPrwi+vr5YtGgRpk+fXuqC5KW5cuUKrl+/Dm9vb2zduhXu7u6K9b1V1b17d6xYsQIjR46Evb29Ynt0dDRatGiBtm3bIiIiAmfOnFH6pa6JuJo3b46nT58iLS0NlpaWxcbSqVMnbNiwAePGjYOpqSnWrVunUsve3Nxc8Yy7OJ06dcK6deuwdu1aJCYmYs+ePVi5cmWhci+3ygGgffv22LRpE1asWIGYmBgcOnRI0cL39fXFzz//jO3bt5fp+6FuvReio6MxZ84cNG/eHKdPn1ZaOtbc3FxprEBp57p8+TJ8fHzKHANVbzGp+Rix8VlrtyI1tJPiyIyGeM3BpELPowlM1NWQv78/7ty5g9deew1WVlYwMDDAgAED8OOPP5b5WF5eXpg5cyYePHgAIyMjHDp0qMzHGDVqFIKDg+Hq6gpbW1vo6ekhODgYkydPRu/evbFlyxZYWFigefPmGo9LX18fo0aNwoEDBzBu3LhiYxk/fjzOnz+PevXqQV9fH+3atcOGDRtUiiMvLw/29vbo1q1bke8Fr1ixAm+++aaie37y5MkYMmRIqcdetmwZBg0ahI0bN8LOzg7NmzdXtEhXrFiBiRMnws7ODra2tpBIJPjiiy8wduzYEo+pbr0XPvroI7z++uuwsLBAkyZN0KRJE8W+gQMHokePHvjtt98we/bsEs8ll8tx7NgxLFy4UKXzUs0QeD8d/psfIS6tABbGelgypB7sLDQ/cMtIX4I+HpawNNFuo0pllfLUvJooy8N/XSCTyURERIQICwsTT58+FUIIkZWVJWJjYxVl0tPTRUJCguJzcnKySE5OFuL5YJ++ffsKIYRISkoScrlcUS47O1vExMQUecxX97+Qn58vHj9+rBjAJYQQBQUFIj4+XnGOF/e2qPovlBRXUZ48eSI6duyoVK6oWF5cS1pamlL9V68vJydHKTa5XC5iYmJETExMkffiheTkZJGbm1to24v7XZzExEQRGxsrmjRpIi5cuKC0LycnR4SHh4uwsDBF3ElJSSIlJUVRJiUlpdA5iqr3qleP80J2drbi+xQdHa10TXl5eSIiIkLpfEWd65dffhHvv/9+iddNNYdMJhdfHYsWelOfDezyXHxbPIjN1nZYFaos+YSrZ5VBTVs968SJE1i9ejVOnDih7VCUqBNXbGws7OzsyrSsqbbdvHkTgwcPhkwmQ2JiImbOnImvv/5a22FpREJCAiwsLCCVSrUdCmlZcmYBxm4LR8DNVADAuI42+HGUU6ljcao6LnNZQWpaos7OzkZaWprSs1xdoKtxaVpeXh6io6Ohp6cHBweHQgPliKq6a48z8eaGRwhPzIPUQIK1/k54p7ONzo6+1qSy5BM+o6ZimZiYwMRE9wZa6GpcmmZkZFRogBlRZbofm4Npux7j70eZFXL8vIJn7UQ3WyP8PrUhWjmVfYBjTcBETURESuRyge/OPMXHB6OQk1+xna7DvKyxbZwzrE2ZjorDO0NERAqh8bmY8FO44vWoPh6W+L/hDWBVASOkjfQlqGPJRzqlYaImIiLI5QLrz8Xjw31RyMqTw0yqh1VvNsCUrrY14pmxLmOiJiKq4R4n5mLiz49x+t6zmcC6NzHH1nEucLXlqHxdUL3Hv5NOO3z4MPbu3avtMMosKSkJixYtqtBznDhxAseOHavQcxDh+RKPnl/cwel76TAxlGDNCEec/k8TJmkdwhZ1NXTr1i0cOnQIn3zySaWeNyAgAKmpqYrFQEpz584dZGRk4K233qrw2DTpq6++Ukx1GhgYWGgGs6VLlyIiIgIbNmyARCKBqakp3Nzc8MYbb8Dd3b1QPSsrK/j4+MDf31/RxdiuXTt07doVffr0UVpUhEiTjtxIwZsbHyGvQKBTQzP8NM4Fje1VW1OcKg9b1NVQbGwsAgMDK/289+7dQ3BwcKWftzJlZ2djz549itW3QkND8eDBAwwbNkzxVatWLcX2oUOHom3btnjy5Am8vb0Vyfnlei1atMBHH32E7777TnGe2rVro0WLFkorZhFp0oHgZLyx4VmSfsPbGmffb8okraP4p7qGCSEgz8mq0HPoGZuWOLhj+fLluHXrFkaOHAl3d3dMmDAB8+bNg0QigZ2dHfz9/dGxY0cAwI0bN3D8+HF4eHjg6NGjmDt3Lpo2bYr9+/fj5MmTaNasGRwdHZGfn69o+V65cgV79uxBVlYW+vfvjyFDhiAyMhK7du1CdnY2IiMj0bdvX0yYMKFQbAcPHsTx48fRrFkzyGQypX1FHfflfb/++iuioqJQu3ZtrF27FgDwxx9/4NChQzA2NsaECRPQvHlz5OXlYerUqVi3bh2MjZ/94jly5AgSExMxfvz4Iu/Zy3G9er0vO3/+fKEFQOrVq4eRI0cWKvvq9oEDB2LkyJGKJP/yfiEEjhw5olh0A88XtTh8+DD8/PyK/V4TqeP3a8nw3/wIBXJgRJta2PGOKwz1OWBMV7FFrWHynCyEdDWv0K/S/hDo2LEj7OzsMGzYMHTt2hVWVlYYNmwYhg4dioYNG8Lf3x/3798Hni81uHz5cuzbtw9du3aFjY0Ndu/ejTlz5qBly5aIjY3FlClT8O+//wIA/vzzT7z77rto1KgRvL29sXDhQuzatQsWFhZo1qwZGjdujGHDhqFly5aF4tq7dy/ee+89tGjRAjExMVi9erViX3HHBYCTJ0+if//+sLa2xtChQ9G3b18AwL59+zB+/Hg0bdoU5ubm6Ny5s2KRjuzsbKXn34sXL1Z0O5cW18vX+6qbN28qLUQBAMHBwRg5ciRGjhyJ9957r9jvy5AhQ2BmZobLly8X2hcdHV1odqLXXnut2DiI1LXnShJGPk/So9vVxk4maZ3HFnU15OPjg3/++UepNaevr4/AwEAkJydDKpXi7NmzaNq0KQDA3t4e27dvV5TdsmUL1qxZg2HDhgEAIiMjFftWrVoFMzMz/PXXXwAAU1NTHDhwAKNHj4aXlxdiY2MV5504cSIyMzNhZGSEn3/+GRs3bsR3332H119/HQAQFRWl0nG//fZbfPPNN4Va6OvWrcPKlSsxatQoAEBKSgq2bduGZcuWYcaMGViwYAHGjh2L69evQyKRoEOHDkXer5LielVmZmahJRvt7e0V96qktZ0BwM7ODmlpacBLCT4pKQmhoaE4cOCAUllTU1Okp1fMerxUM+28mIhxP4VDLoDxHW2weawz9PWYpHUdE7WG6Rmbwut8xa6jqmdctmn2tm/fjq+//hrTpk2Dra0t0tLSlBKAq6urUvn4+HilqStdXV0V3dTR0dHw8/NTap06ODgUed6BAwciLy9PsRDG06dPlc6l6nFjY2MVf1S87NXjNWzYEDdu3AAAdO3aFRkZGbh9+zbWr1+PWbNmFXt/SorrVXZ2dggPD1faVlzX96tyc3MRHh6O+vXrK+YqHzZsGCIjI/H999/j1Wn3ExMTq/185lR5fvo7Ae/8/BhCAJO62GLDaCfoMUlXCUzUGiaRSKBvUnKrqqJJpVLk5uYqPl+5cgXvvPMO3nvvPeTn5+Obb74psb67uzsCAwPh5eUFuVyOs2fPomvXrsDz9ZeTkpKKTEyvnvdFC/UFDw8PnDp1Cq1atYJcLseZM2fg4+NT6nFbtmyJQ4cOoVOnTkUer2PHjhBC4I8//kCPHj0U+999912sWrUK586dww8//FDs9ZYU16vat2+PzZs3l3D3ilZQUICPPvoI9vb2aNeuHe7evauU4OvVq4dJkybh8uXLivEHISEhirEEROWx9UICJu14lqSndbPFWn8m6aqEiboaatasGW7fvo2hQ4eiVatWGDRoEEaMGIFz584pnuGW5OOPP0bv3r1x/PhxxMfHw8DAQFHnyy+/RJ8+feDl5YWmTZtCIpFgxIgR8PPzQ/v27fH5558jNjYWAwYMKNRV/emnn6Jnz544deoU4uPjlZY4LOm4X3zxBXx9fREYGAg3NzfY2Nhg7dq1WLhwIXr16oWzZ88iJSUFQghMnTpVcczRo0dj3rx5mDVrVonXXFJcr/L29lYMmGvQoEGJ9/FF13ZKSgpu3boFKysrHDx4sMhVsEaNGoXvv/8e27dvVwx4O3z4sNJzfCJ1XAnPxNSdz5L0zO52WDPSkTONVTFM1NVQ7dq1cevWLVy+fBkWFhbw9fXF1atXcf/+fbRs2VLxnBrPW7ILFixQqu/l5YU7d+7gypUrcHd3x6xZs9CoUSMAQIMGDXDjxg1cuXIFkZGRkMvlim7pdu3a4eLFi7h161ah7nQAaN68ueK4TZs2RW5uLgoKCko9rrOzM27duoVLly4hNjZWsXLWa6+9hrt37+LixYswNjZGx44dlZKgVCqFg4MDpk2bVuL9ejWuvXv3IiOj+McXH3zwAdavX4+lS5eiZ8+eRXbL9+zZE+bm5pBIJDAxMYGbmxuaN2+utP/Vejt37sTjx4+B54PWjI2N4e3tXWLsRCXJyJFh9JYwFMiB4a1rMUlXUVyPugxqynrUcXFxmD17NuRyOe7evQupVIoLFy6U2NLUNVu2bMGPP/6Idu3aYd26dWWqu3z5cmRkZGDp0qVF7pfL5YqR6BXl1q1bsLCwgLOzc4Wdg6q/yTseY3NQAhrUMsS/n3mglhnbZrqiLPmEiboMakqizsjIQEBAAPT09FC/fn20b9++ys2OdfXqVcTHx6N3795ljv3u3bsoKChQzD5GVBUdCE7G6+sfQSIBzvynCbo3tdB2SPQSJuoKUlMSNRFVbdEpeWjxxR0kZsrwUR97rHij5PEUVPnKkk844QkRUTUilwuM/ykciZkyeDuZYsnQetoOicqpSvRnJicn4/z580hOToaTkxO6deumeDf3hby8PJw+fRpxcXHw9PRE69atCx1HlTJERFXZmjNPcerus5Wwdk10hZEB22NVnc4n6mPHjmHEiBHw9vaGi4sLgoKCIJVKERgYqJgMIj4+Hj179kReXh5atGiB//znPxgxYgTWr1+vOI4qZYiIqrJ/I7Mw78CzmfVWDXeEe10uslEd6HyiXrBgAV5//XXFFJeZmZlo2LAh1q1bp1gT+MXrRcHBwTA1NcW1a9fQtm1bDB48GAMHDlS5DBFRVZWTL8foLWHIKxAY5GmFad1stR0SaYjO94lIJBLY2v7vP5ypqSnMzc2hp/csdLlcjr1792LChAmKOZhbt26NTp06Yc+ePSqXISKqyubvj8Kt6BzUsTDAlrHOfF+6GtH5FvXatWsxadIkyOVyODs74/Tp02jcuLFi7uaIiAikp6fDw8NDqZ6HhweuXr2qcpmi5ObmKk2J+WIxBSIiTXqSlIdv/4zD7stJyMqTq3WMjNxn9baNc0Edy8Kz31HVpfOJ2tzcHFZWVrh8+TISExPx4MEDdOrUSfFu7IvkaW1trVSvVq1ain2qlCnKsmXLsHjxYo1fExERnj9T/uaPOPxyJQky9fKzkjm+dTDA00oToZEO0elELZPJMGTIEAwbNgzfffcdACA7OxutWrXCxx9/jDVr1iimk3x1OcC0tDRFN7cqZYqyYMECzJ07V6m8o6OjBq+QiGoaIQQC76fj65NxOHnnfw2Fnk0t8H5vezRVcwCYob4EjrXYkq6OdDpRR0dH48mTJ0qDvUxMTNCjRw/8/fffwPN5oI2MjBAWFqZUNywsTDE/tSpliiKVSqvUtJlEVPH+Cc3AH3fUewwmF8DRm6m49iQLAKAneTYH94d97NHaWbur7pHu0ulE7eDgAKlUiuDgYPTp00exPSQkRLFesqGhIfr164c9e/Zg8uTJkEgkiI6ORmBgIDZs2KByGSKi0py8nYqBPzwsdze1iaEE73S2xdxe9nCzY2OASqbzU4iuXLkSn376KSZNmgQ3NzecOnUKQUFBCAoKQsuWLQEA9+7dQ6dOndClSxd06NABP//8M+rUqYMzZ84onmWrUqY0nEKUqOa6GZWNzl/fQ3qOHN2bmKv9jrKzjRSTutjC1lyn20lUwardXN9XrlzByZMnkZiYCGdnZ/j7+ysmO3khMjISP//8s2LWsbFjxxZag1iVMiVhoiaqmWJS89F+2V1EJOfDp4k5/pjdmDN+UblUu0StK5ioiWqezFwZuq96gKuPs9DEXop/5rmjNpeLpHLiohxERBogkwu8vTUcVx9nwcZMH8fea8wkTZWOiZqIqBjz9kfhYEgKjAwkOPRuIzTkwC/SAiZqIqIirP8rHqtOxQEAfhrngs6NzLUdEtVQTNRERK84cSsVM/c8AQAsGVIP/u1qazskqsH4sIWIapT0HBmWHo3B5fDMYstcfZwFmRwY19EGnwyoW6nxEb2KiZqIaoxjN1MxbddjRCTnl1rWp4k5Nr7txFWoSOuYqImo2otPz8ecvZHYfTkJAOBqa4RP+jvAXFr00z9TIz30bWbJd6VJJzBRE1G1JYTA7stJmP1rBBIzZdCTPFth6osh9WAm1dd2eEQqYaImomrpcWIupu9+guO3ni2g0aKBCTaPcUZbFy5+QVULEzURVTl3Y7Lx/u+RuPY4q9gyKdky5BUIGBlI8PlAB3zUty4M9fm8maoeJmoiqjKy8uRYejQGK0/FIV9W+uzHnRuaYfNYF7UX0CDSBUzURFQlHL2Zipm/PEF4Yh4AYJCnFT4d6AAzo6IHfEkNJGhUR8pR21TlMVETkU6LSMrD7F8jcCAkBQDgWMsQa0Y6YWhLKyZhqhGYqIlIJ8nkAqtPP8XCI9HIzJXDQA/4Ty97fD7QAebGHLFNNQcTNRHpHCEEpu16gs1BCQCALo3MsW6UE5rXN9F2aESVjomaiHTO1yfjsDkoAXoS4MdRTpjcxRZ6euzmppqJiZqIdMpv15Ix/0AUAOC7EY6Y2s1O2yERaRXnxyMinXHxUQbGbgsDAMzuWQcze9TRdkhEWsdETUQ6ISwhF0N+DEVOvsDgFlZYNbyBtkMi0glM1ESkdSlZBRjw/UPEpxeglaMJdk90hT6fSRMBTNREpG15BXK8sf4R7sXmoL61IY7MaMTXr4hewkRNRFojhMD03U9w5n46zKV6ODqzEerXMtJ2WEQ6haO+iUgrMnJkWHosBlsvJEJPAvw62Q0tHU21HRaRzmGiJqJK9Sg+F2vPPsWWC4lIzZYBz1/DGuBppe3QiHSSyom6TZs2Kh/06tWr6sZDRNWQEAJn7qVjTeBTHPk3FeL5wleN6kgxv29dTOxiq+0QiXSWyol65MiRin9HRkZizZo16NOnD9q2bQsAuHLlCk6ePInZs2dXTKREpJP+G5eDkMjsYvfHpuZj/bl43InJUWzr62GJWT3roF8zS844RlQKlRP1Bx98oPj34MGD8eOPP2LatGlKZdatW4fjx49rNkIi0kn5MoEvj8Xgy2MxKJCXXt5MqofxHW0ws0cdrg9NVAYSIUTpq6+/onbt2ggPD4elpaXS9tTUVLi5uSExMVGTMeqMtLQ0WFlZITU1tdC1E9UkN6OyMW5bGIIjnrWkvZ1MYWFc9EskhvoSDGxuhQmdbWFlwteuiFDGfKLWYDJ9fX1cuHAB/fv3V9r+999/w8CA49OIqqsCmcDXJ2OxKCAG+TKB2mb6+NHfCSPa1tZ2aETVllpZ9b333sPIkSMxY8YMtG3bFkIIXL16FWvXrsWHH36o+SiJSOvuxmRj/E/huByeBQAY0tIKG0Y7o66VobZDI6rW1Or6BoCffvoJ3377Le7fvw8AaNq0KebOnYtx48ZpOkadwa5vqqoycmS4/iQLeTK1ftxx7XEWFh6JRm6BgJWJPtaMcMSYDrUhkXAgGJE6ypJP1E7UL8jlz0aR6OlV/0nOmKipqpDLBW5EZuOPO2k4eScNQQ8zkK9mkn5ZXw9LbB7rjAacPYyoXCr8GfXLMjMzYWFhUd7DEFE5JWQU4PitVPxxJw1/3EnD0/QCpf2OtQxRy0y9H3mpgQRTutphYmcbtqKJKplaP7V5eXn4/PPPsWXLFiQkJOBFo3zKlCl4//330bRpU03HSURFkMsFAu+nY+P5BBwISVFqNZtJ9dCzqQX6eFiir4clGtWRMskSVUFqJeolS5bg+PHj2LJlC4YOHarY3rt3byxevBi7d+/WZIwAgJSUFFy8eBGmpqbo2LEjDA2VB7DI5XJcunQJcXFxaN68ORo1alToGKqUIaoKnqbl46d/ErEpKAEPn+YqtrdsYIL+za3Q18MSnRqawcig+j+SIqru1HpG7eLigqNHj6JZs2aQSCSKFnVcXBzc3d2RnJys0SDXrl2L+fPnw9vbG6ampkhISMDBgwdRv3594Pn72/369cOTJ0/g4eGBv//+G++99x6WL1+uOIYqZUrDZ9SkTXK5wJnnreeDL7WeLY318HZ7G0zuagsvLmpBVCVU+DPqmJgYuLq6AoBSV5pEIkFOTk4JNcvu8OHDmDVrFgICAhTvbd+5cwfZ2f+bsvCTTz5BUlIS7ty5AysrKwQFBaFr167o3bs3fH19VS5DpItiU5+3ns/H41FCnmJ7e1czTOlqixFtasFMyolEiKortRK1u7s7zp07h379+ikl6u3bt8PLy0uT8eGrr77CsGHDlCZX8fDwUPxbCIFdu3ZhwYIFsLJ6tvpOly5d0K5dO+zcuRO+vr4qlSHSJXK5wJ/30rHxfDwOhaQopui0NNbDmA42mNzFlktCEtUQaiXqTz/9FGPGjMG8efMAALt27cKJEyewe/duHDp0SGPBZWdn48qVK5gwYQLCw8Nx48YN1KtXD97e3tDXf9aCiIyMREpKCpo3b65U19PTEyEhISqXKUpubi5yc//3/C8tLU1j10ZUlNjUfGz7OwGbghIQ9lLruYOrGaZ2s8Xw1mw9E9U0aiXq4cOHw8DAAEuXLgUAjBkzBp6envj9998xaNAgjQWXmJgIuVyOP/74A8uXL4enpydu3LgBKysrHD16FI6OjkhNTQWezz/+MhsbG6SkpADPn0+XVqYoy5Ytw+LFizV2PUTFuRyWiTVnnmLvtWTFs2crE32MaV8bk7vaokUDtp6Jaiq136P28/ODn58fZDIZhBAVMse3sfGzFXbu37+PO3fuwMTEBDk5OejUqRPmzp2L3377DVKpFACQkZGhVDcjI0NRX5UyRVmwYAHmzp2r+JyWlgZHR0cNXiHVZHkFcvx+PQVrzjzFpbBMxfYOrmaY5mOH4a1rwdSIo7aJarpyZ9cXXdAVwdbWFlZWVujfvz9MTEyA58l74MCB2LVrFwDAyckJBgYGePLkiVLdx48fw83NTeUyRZFKpYokT6QpcWn52HAuHuv+ikds2rNJSYwMJBjZphbe61EHbVzMtB0iEekQlRN1ly5dVD5oUFCQuvEUMmjQINy7d09p2927dxUtW6lUCl9fX/z222945513AAAJCQk4c+YMvv32W5XLEFWGTefjMXNPBPIKnnVvO1gZYrqPHaZ0tYW9JRe3IKLCVH6PetGiRSoftCxlSxMWFob27dtj8ODB6Ny5My5duoTt27fj5MmT8PHxAQCEhISgS5cu8PPzQ8eOHbF582ZIJBL8/fffihaxKmVKw/eoqTz+71Qc3v89Enj+atXsnnXwhrc1JyUhqoEqdVGOyhAdHY0NGzbg8ePHcHJywtixYwvNKvbgwQNs3rwZcXFx8PT0xPTp02FmZlbmMiVhoiZ1CCGw5GgMFh6JAQDM62uPZX71OZ0nUQ1W7RK1rmCiprISQuCjfVFYeSoOALB0aD183L8ukzRRDVchM5Np6xk1UVUllwvM+OUJ1p9LAACsfqsBZvvaazssIqpiVE7UvXr1qthIiKqRApnAOz+HY8fFJEgkwKa3nTGxi622wyKiKohd32XArm9SRW6+HKO2hGF/cAr09YAdE1zh3662CjWJqKao8EU5iKgwIQT+epCBzw5HI+hhBowMJNg72Q1Dvay1HRoRVWFqJ+q4uDicPHkST548QUFBgdI+Tb6eRaTrZHKBA8Ep+PqPWFwJzwIAmBrp4eD0hujtwZ4XIioftRL1uXPnMHjwYDg6OuL27dto3bo17t+/j4yMDLRt25aJmmqE7Dw5tv+TiJWn4hAa/2zxFmNDCSZ0ssUHve3hZsdZ7Yio/NRK1PPmzcMXX3yB2bNnQyKR4OrVq0hNTcW4ceNKnJKTqKq4FJaJm1HZxe6PSMrDunPxiE9/1ptUy1QfM3vUwczudqjDGcaISIPUGkxmbm6O6OhoWFpaQl9fH1lZWZBKpXjy5Anat2+PmJiYiolWyziYrPq7G5ONj/ZFIeBmqkrlnWob4f3e9ninkw3Mjbn8JBGppsIHk2VmZioObG9vj/DwcDRt2hTm5uYlLhtJpKtiU/OxKCAam4MSIJMD+npAL3dLSA2LnpjESF8Cv1a1MLx1LRjqc/ISIqo45R713bt3b8yePRuTJk3C7t270aZNG81ERlQJMnNl+L8/n2LFyVhk5soBAMO8rLHcrz6a1i1+CVQiosqiVqJet26d4t9ff/01Jk2ahGnTpsHd3R2bNm3SZHxEansUn4uI5Lxi99+LzcEXR2MQnZIPAGjrYopVbzZA18YWlRglEVHJOOFJGfAZtW6TywUuh2fi0I1UHL6RgjsxOSrVc7ExwjK/+nirdS3o6bEbm4gqXqVMeCKTyaCvr1/qNqKKlJ0nx+l7aTh8IxVH/k1BbNr/3unX1wMa2UmLTb5SAwnebm+Dmd3tIDXkUpNEpJvUStQ7d+7EiRMnsHPnTqXtY8eOxcCBAzFq1ChNxUdUSHx6Po7eTMXhG6k4eScNWXlyxT4LYz30b2aFoV7W6N/MErXMOPkeEVVtanV9u7u74/Dhw2jSpInS9gcPHmDYsGG4c+eOJmPUGez61p4HcTk4FJKCw/+m4u/QDMhf+l/boJYhhra0xpCW1ujexBxGBmwdE5Fuq/Cu7/DwcNjZ2RXabmtri0ePHqlzSKJCkjIL8P2Zp/jlShLux+Uq7WvlaIIhz5NzK0cTru9MRNWWWonaw8MDe/fuxdSpU5W2//rrr3B3d9dUbFRDxafn4//+fIofAp8i4/krUwZ6QI+mFork7FTbSNthEhFVCrUS9YIFC/D222/jxo0b6NatG4QQOHfuHLZs2YJdu3ZpPkqqEeLS8rHyjzj8+Fe84rlzywYm+KC3PQa3tIaVCQcqElHNo1aiHj58OABg6dKlWL9+PQDA09MTu3btUuwjUlV0Sh6+PhmHDefjkZP/7OFzaydTfD7IAYNbWLFbm4hqtHK/R11QUACJRFIjXsviYDLNSs+R4ctjMVh9+ilyC579N+zgaobPBzmgXzNLJmgiqrYq5T1qxQEM+PoLlY1cLrDzUhLmH4hCTOqzWcG6NDLHwkEO8HW3YIImInoJsyxVqivhmZi1JwIXwzIBAA3tpFj9VgMM9GQXNxFRUZioqVLEpeVjwYEobPs7EQBgLtXDpwMcMMe3DmcFIyIqARM1VagCmcD3gU+x6Eg00nKejeQe26E2lvnVRz1rvmJFRFQatZoyGRkZ2Lt3r+JzQEAAevXqhalTpyI9PV2T8VEVdiMiCx1X3MPc3yKRliNHG2dT/DOvKbZPcGWSJiJSkVqJ+uOPP0ZCQgIAICkpCf7+/nB1dUVwcDA+/PBDTcdIVUxOvhyfHIxCm6/u4urjLFib6mPTGGdcmu+ODm7m2g6PiKhKUavre9++fQgJCQEAnDx5Em3atMGmTZvw6NEjdOnSRfFuNdU85/+bjsk7Hium/HzD2xrfj3SCg5WhtkMjIqqS1ErUaWlpiteyzpw5g759+wIA6tSpw67vGiotW4b5B6Kw7q94AEBdSwP8OMoJfq1qaTs0IqIqTa1E3bp1a3z66afo27cvfv31V1y4cAEAEBISglatWmk6RtJx159kYeiPDxGZ/Oyd6EldbPHNG/VhbcqxikRE5aXWb9LvvvsOb7/9Nnbs2IFZs2bB09MTALB8+XLMmzdP0zGSDkvOLIDfulBEJufDzdYIm8Y4o6c7Z20jItKUck8h+rLMzEyYmZlp6nA6h1OIKhNC4I31j3AgJAUN7aS49slrXDiDiEgFZckn5Z5p4uVn0tU5SVNha8/G40BICgz1JdgzyZVJmoioAqiVqPPy8jB//nzY2dkp/SUwZcoU3L9/X5PxkY4KfpKF93+PBAB880Z9tHHhH2lERBVBrUS9ZMkSHD9+HFu2bFHa3rt3byxevFhTsZGOSs+R4a1Nj5BXIDCkpRVm9ayj7ZCIiKottQaT7dixA0ePHkWzZs2Utnfr1g1TpkzRVGyFxMfH49KlS3Bzc4OHh0eh/Tdv3kRcXBw8PDxQr169Io+hShkqnhAC03Y9wcOnuXCsZYht41y4mAYRUQVSq0UdExMDV1dXAFD6JS2RSJCTk6O56F4il8sxYsQI+Pn5YePGjUr7MjIy4Ovrix49euCzzz5Dw4YNsXTp0jKXodJt+zsRuy8nQV8P+GWSG2qb8RUsIqKKpNZvWXd3d5w7dw79+vVTStTbt2+Hl5eXJuNTWLp0KWrVqlWoFQ8An332GcLCwnD//n3Y2Njg1KlT6NOnD7p164Zu3bqpXIZKdjs6GzN/eQIAWDKkHjo34nSgREQVTa0W9aeffooxY8Zg5cqVAIBdu3ZhzJgxmD9/Pj755BNNx4igoCBs2bKlUEsaz7tid+zYgUmTJsHGxgZ4/qzc29sbP//8s8plqGRZeXKM2PQI2fkCvV+zwLy+dbUdEhFRjaBWi3r48OEwMDBQdB2PGTMGnp6e+P333zFo0CCNBpiUlIS3334bW7ZsUSTZl0VFRSExMREtW7ZU2u7l5YUbN26oXKYoubm5yM3NVXxOS0vTwBVVPZm5Mkzd+QS3o3NQ19IAO95xhZ4en0sTEVUGtR8w+vn5wc/PDzKZDEIIxdzfmjZx4kS88cYb6NWrV5H7U1JSAAC1a9dW2m5jY6PYp0qZoixbtqxGj2KXyQV++jsRnx+JRnRKPiQSYOc7rrC35AIbRESVpdwTnujr61dYkv79998RGBiITp06ISAgAAEBAUhLS0NYWBgCAgIghICR0bN1jbOyspTqZmVlKfapUqYoCxYsQGpqquIrIiKiAq5SN528nYpWS+9i0o7HiE7Jh4uNEfZPawjf1zgjGxFRZVIrwyYmJmLRokW4cOECkpKSCu0PDw/XRGwwMTFBly5dsG3bNqVz37x5E+vXr8eAAQPg5OQEfX19REZGKtWNiIhQjExXpUxRpFIppFKpRq6lqrgRkYUP90Xi1N1nM85Zm+rj0wEOmNndDlLDcv9dR0REZaRWop48eTLu37+P8ePHo1atilvGcODAgRg4cKDSNi8vL3Tv3h2rV68GABgbG8PHxwf79+/HuHHjgOdd3adPn8by5ctVLlPTpWQV4D97I7H9YiKEAIwMJJjZ3Q6fDHDgK1hERFqk1m/g06dPIyQkpMTWaGVatmwZunXrhunTp6Njx45Yv349XFxcMHHixDKVqamEEBj3UzgO30gFAIxoUwtfDasPN7ua1ZtARKSL1OrLrFWrFqytrTUfjQq6detW6F3qdu3a4fLly5BIJDh06BD69OmDoKAgmJiYlKlMTbXvegoO30iFob4EZ99vgj2T3ZikiYh0hFrLXH700UcwMjKqcTN7VcdlLpMzC+Cx6DZi0wrw2UAHfDGE06oSEVW0suQTtbq+r1y5grNnz2LPnj1o2LBhobmeT5w4oc5hSQs+2h+F2LQCNLWX4uP+nMSEiEjXqJWoW7dujdatW2s+GqpUZ++nY3NQAgBg0xhnGHNUNxGRzlErUb+YOpSqruw8OabsfAwAmNrNFl0bW2g7JCIiKkK5m1Dp6emaiYQq1dJjMfjv01w4WBlixesNtB0OEREVQ61EnZeXh/nz58POzk7pIfiUKVNw//59TcZHFeDfyCx8fTIWALDW3xFWJvraDomIiIqhVqJesmQJjh8/ji1btiht7927d42eG7sqkMkFJu14jAI58Hora/i1qrgJa4iIqPzUeka9Y8cOHD16tND7zN26dcOUKVM0FRtVgO/PPMWV8CxYmejj+5GO2g6HiIhKoVaLOiYmRjEr2cuvZkkkEuTk5GguOtKo8IRcfHo4GgDw9ev1Uc+6+AVJiIhIN6iVqN3d3XHu3DnglUS9fft2eHl5aS460qgZvzxBZq4cXRuZY1IXW22HQ0REKlCr6/vTTz/FmDFjMG/ePADArl27cOLECezevRuHDh3SdIykAeEJuTh2Kw36esDGMc7Q05OoUIuIiLRNrUQ9fPhwGBgYKKYQHTNmDDw9PfH7779j0KBBmo6RNOD0vWev0bVzMYN7XWNth0NERCpSe/1CPz8/+Pn5QSaTQQgBAwMuhajLztx/lqh93TmxCRFRVVLu7Kqvz3dwdZ0QAmfupQEAfN2rx2IiREQ1hVqJuk2bNsXuk0qlcHNzw/jx4+Hr61ue2EhD7sTkIDatAMaGEnRwM9N2OEREVAZqjfru0qULrl+/Dnt7e/Tv3x8DBgxAnTp1cP36dXh6eiInJwd9+/blwDIdceb58+kujcy58AYRURWjVos6NDQUGzduxKRJk5S2b9y4EQEBATh8+DA2bNiAL774AkOHDtVUrKSm0+z2JiKqsiRCCFHWSrVr18bjx49hYaE8MCktLQ2urq5ITExEcnIynJycqtWiHWVZ6FtXFMgEbN+/gdRsGS4vcEdbF3Z9ExFpW1nyiVr9oPr6+vj7778Lbb9w4YJicFlycjLq16+vzuFJg64/yUJqtgxWJvrwdjLVdjhERFRGanV9T58+HSNGjMCMGTPQpk0bCCFw7do1/PDDD5gzZw4AYP369Zz3Wwe8eC2rexNz6HOSEyKiKketRL148WI4OTnhu+++w8qVKwEATZo0wapVqzBx4kQAwJw5c+Dg4KDZaKnM+HyaiKhqU+sZ9cvkcjkAQE+v+o8mrmrPqHPy5aj1nxDk5AvcWuiBZvVMtB0SERGVMZ+Ue8KTmpCgq6qLjzKRky9Q19IAHg6cNpSIqCpilq3GXnR793S3VFrljIiIqg4m6mrsxUIcnN+biKjqYqKuptKyZbgcngkA6NmUiZqIqKpioq6mzj/MgEwOuNkawcVWqu1wiIhITUzU1dTpu3wti4ioOmCirqZOP5/opCefTxMRVWlM1NXQ07R8/BuZDfD5NBFRlcdEXQ2dffCsNe1Z3wR1LA21HQ4REZUDE3U1xNeyiIiqDybqauhFoma3NxFR1cdEXc08TsxFaHwu9PUAnyZM1EREVR0TdTVz5nlruq2zGSxN9LUdDhERlVOVSNT5+fm4efMm7t69i7y8vGLLhYWF4eLFi0hOTi5XmaqMz6eJiKoXnU7UQggsXrwYDRo0wOjRozF48GC4uLjg0KFDSuWys7MxdOhQeHp6YvLkyahXrx6+/fbbMpep6oQQLy3EwURNRFQdlHuZy4okhIAQAvfv34e1tTUAYMmSJfD390doaCgcHBwAAIsWLUJwcDBCQ0Nhb2+PI0eOYMiQIejYsSM6dOigcpmq7l5sDmLTCmBsKEGnhubaDoeIiDRAp1vUenp6WLRokSJJA8CUKVOQnZ2N4OBgxbaffvoJkyZNgr29PQBg8ODBaNGiBbZt21amMlXdi27vzg3NYWyo099aIiJSkU63qIty+fJlAECjRo0AAFFRUXj69Clat26tVK5169aKZK5KmaLk5uYiNzdX8TktLU2j16JpB0NSAAC9X+P83kRE1UWVanY9ffoUM2fOhL+/P5o0aQIAikFhtWvXVipra2ur2KdKmaIsW7YMVlZWii9HR0eNX5OmxKfnI/D5/N5vtq6l7XCIiEhDqkyiTk5ORr9+/dCgQQNs2rRJsd3IyAh4PljsZVlZWYp9qpQpyoIFC5Camqr4ioiI0Og1adKB4BTIBeDtZIqGdlzWkoiouqgSXd8pKSno06cPjI2NceLECZiZmSn2NWjQAHp6eoiKilKqExUVBWdnZ5XLFEUqlUIqrRpJ77frz3oGhrM1TURUreh8izo1NRV9+vSBgYEBTpw4AQsL5deOTE1N0blzZxw+fFixLSMjA3/++Sd69+6tcpmqLCGjQNHtzURNRFS96HSLOj8/H/369cPjx4+xadMmXL9+XbGvadOmitezvvzyS/j6+uLDDz9Ex44d8cMPP6Bu3bqYPHmyorwqZaqqA8HJkMmBVo4m7PYmIqpmJEIIoe0gipOamoqhQ4cWue/DDz/EwIEDFZ8vXryIH3/8EXFxcfD09MRHH32EOnXqKNVRpUxJ0tLSYGVlhdTUVFha6s7I6j6rH+DU3XR8NaweFvR30HY4RERUirLkE51O1LpGFxN1QkYB6n54AzI58N8lzdCojrG2QyIiolKUJZ/o/DNqKtnL3d5M0kRE1Q8TdRX32zWO9iYiqs6YqKuwhIwCnOFobyKiao2Jugo7GJICmRzwYrc3EVG1xURdhSm6vb3ZmiYiqq6YqKuoxIwCxdrT7PYmIqq+mKirqBfd3i0bmKCxPbu9iYiqKybqKupFt/dbbE0TEVVrTNRVELu9iYhqDibqKuhgSAoK2O1NRFQjMFFXQZzkhIio5mCirmKSMtntTURUkzBRVzEvur1bNDBBE3Z7ExFVezq9HjX9jxAC4Yl52P5PIsBJToiIagwmah2UlSvD7ehs3IjMxo2ILIREZuPfyCyk5cgVZd70tgZXKCUi0g6JRFJp52Ki1jHHt+yA1frJkIpctALQqphymYOB65UcGxERAbV6vwW3Zb9W2vn4jFrHpO37AVKRq+0wiIhIR7BFrUNiI6Lh+vQqAMBq7SU4N3HVdkhERPQKPUNppZ6PiVqH/PPrXjhBjsdWr+H19u20HQ4REekAdn3rkOygQwCA3NaDtR0KERHpCCZqHZGckAzXqAsAAM/XR2g7HCIi0hFM1Driwt59MBL5iDZzgUc7L22HQ0REOoKJWkeknj0AAEhrORB6evy2EBHRM8wIOiArIwtO4YEAgCZD3tJ2OEREpEOYqHVA0L5DMJVnI0Fqj9Y9u2g7HCIi0iFM1Drg6an9AID4ZgPY7U1EREqYFbQsPzcf9R+eAgA4939T2+EQEZGOYaLWsr8DTsCyIBUphrXQYUBvbYdDREQ6holayyKO/w4AiGncG4ZSQ22HQ0REOoaJWotkBTLUuXscAFC3D7u9iYioMCZqLbr651+wyY1Hpr4ZOvsN0nY4RESkg5iotehhwLNu7yeuPWBiZqLtcIiISAcxUWuJXC6H1b9HAQC1e76u7XCIiEhHMVFryc1/rsEh6wlyJUbo8iYTNRERFa3GJeqnT5/i1q1byMrK0moct/f/CgAIb9AFlrWttBoLERHprhqTqPPy8jB69Gg4OTlh8ODBqFOnDjZt2qS1eEyCAwAApt2GaS0GIiLSfTUmUS9duhSBgYH473//i7CwMGzZsgVTp07FtWvXKj2WBzfuwCntPgqgj05vcREOIiIqXo1J1Js3b8akSZPg6OgIABgxYgRee+01bNmypdJjCf79Wbd3WN12qFPfvtLPT0REVUeNSNQxMTGIiYlB27Ztlba3b98e169fr/R49C4fBgAYdBhS6ecmIqKqxUDbAVSGxMREAICNjY3SdltbW8W+ouTm5iI3N1fxOS0trdyxRD56ArfEEABAuxEjyn08IiKq3mpEi9rQ8Nkc2i8nXQDIzs5W7CvKsmXLYGVlpfh60W1eHkZSKR72+wx3m/vDqbFruY9HRETVm0QIIbQdREXLzMyEhYUFduzYgdGjRyu2v/HGG8jMzMSJEyeKrFdUi9rR0RGpqamwtLSslNiJiKj6SUtLg5WVlUr5pEa0qM3MzNChQwccPXpUsS07OxunT5+Gr69vsfWkUiksLS2VvoiIiCpTjXhGDQBLlixBv3790LRpU3Ts2BHfffcdrK2tMXXqVG2HRkREVKwa0aIGAF9fX5w8eRI3btzA4sWL4ejoiKCgILaSiYhIp9WIZ9SaUpZnCkRERMXhM2oiIqJqgomaiIhIhzFRExER6bAaM+pbE148ztfEDGVERFRzvcgjqgwTY6Iug/T0dADQyAxlRERE6enpsLKyKrEMR32XgVwuR3R0NCwsLCCRSFSu92JGs4iICI4WLwXvlWp4n1TD+6Qa3ifVaepeCSGQnp6OevXqQU+v5KfQbFGXgZ6eHho0aKB2fc5upjreK9XwPqmG90k1vE+q08S9Kq0l/QIHkxEREekwJmoiIiIdxkRdCaRSKRYuXAipVKrtUHQe75VqeJ9Uw/ukGt4n1WnjXnEwGRERkQ5ji5qIiEiHMVETERHpMCZqIiIiHcZErabQ0FBcu3YNWVlZGq2jznF12dOnT3HlyhU8ffpUY3WEEHj06BFu376NnJwcDUarPVlZWbh27RpCQ0MrpM7169dx4cKFckapfXK5HLdu3cKNGzcgk8k0WictLQ3Xrl1DcnKyBiPWnvDwcFy9ehUZGRkaqyOXyxEaGorr168jISFBg9FqT1JSEq5cuYKYmBiV62RkZODChQuIiIjQ6HGLJahMkpKShI+Pj7CwsBCNGzcWlpaWYs+ePeWuo85xdd2cOXOEVCoVHh4eQiqVijlz5pS7ztatW4Wrq6twdXUVr732mrC2thbr1q2rwKuoeHv27BGWlpaicePGwsLCQvj4+IikpCSN1Tl69KjQ19cX+vr6FXQFleP27duiUaNGom7duqJ+/frCyclJXLlypdx1CgoKxPvvvy9MTEyEl5eXcHJyEp988kkFX03FycjIEP379xdmZmaiadOmwtTUVGzevLncdS5duiQaNWokHBwchLe3tzAxMRGjRo0Subm5FXxFFWfx4sVKv2/Gjx8vCgoKii0fFRUlZsyYIRwcHISRkZFYsmSJRo5bGibqMnr77bdF8+bNRWpqqhBCiO+//14YGRmJsLCwctVR57i67KeffhKmpqYiJCRECCHE9evXhYmJidi+fXu56nz55ZciPDxc8fmXX34REolE/PPPPxV6PRXl0aNHwsjISPz4449CCCFSU1OFh4eHGDNmjEbqREVFiQYNGojZs2dX6UQtk8lEs2bNxJtvvinkcrkQQohx48YJZ2fnYhOFqnU++OADYW9vL+7cuaOo9+LeVkUzZ84Ubm5uIj4+XgghxI4dO4Senp64efNmueq0aNFCvP7664qE8+DBA2FiYlJl/1AOCAgQBgYG4ty5c0I8vx5ra2uxatWqYuucO3dOrFmzRqSkpAhnZ+ciE7U6xy0NE3UZpKenCyMjI6W/NAsKCoStrW2xf1mpUked4+q6bt26iZEjRypte/PNN4WPj49G6wghyv1DoE1ffPGFqFOnjpDJZIpt69evF1KpVGRkZJSrjkwmEz169BCrV68WmzZtqtKJ+u+//xYAFH/ECSHEw4cPBQBx/PhxtevEx8cr/dFT1eXn5wtLS0uxcuVKpe2urq7i/fffL1cdBwcHsXz5cqUyzs7O4osvvtDoNVSW119/XfTq1Utp27Rp00SzZs1Uql9coi7vcYvCZ9RlcOfOHeTl5aF169aKbfr6+vD29kZwcLDaddQ5rq4LDg5Wuh4AaNeuXYnXo06d//73v0hNTUWjRo00EHXlCw4ORqtWrZQm5W/Xrh1yc3Nx586dctX58ssvYWRkhFmzZlXwVVS84OBgGBgYoEWLFoptDRs2RO3atYv9/6FKnfPnzyMvLw+DBw9GZGQkQkJCFKvkVUWPHj1CWlpaoZ+jtm3bFnufVK2zdOlSrF27Ftu3b8epU6cwa9YsGBsbY9KkSRV0NRWruN83d+/eRW5urk4dl4m6DJKSkgAANjY2StttbGwU+9Spo85xdVlBQQHS09OLvJ60tLQiB/SoUycvLw/jx49Hq1atMGDAgAq4koqXlJRU5DXjpf8X6tQJCgrC2rVrsW3btjKt9KarkpKSULt27ULXUtrPXml1oqOjYWBggBUrVqB9+/YYO3Ys6tSpg88++6wCr6biVNTvKADo27cvPD09MW/ePHz00UfYsWMHpk2bBgcHhwq4kopX3M+RXC5HSkqKTh2Xq2eVgaGhIQAUGmmcnZ0NIyMjteuoc1xdpq+vDz09vSKvR09PD/r6+uWuU1BQgJEjRyI6Ohrnzp2DgUHV/K9saGhY5DUDKPH/VGl1Ro0ahXHjxiE0NBShoaF4+PAh8DyBu7m5oV69ehVyPRWlqGuGCj97pdUxNDREQUEB0tLS8OTJE+jr6+Ps2bPw9fVFmzZtMHTo0Aq4mopTUb+jZDIZevXqhdatWyMyMhIGBgYIDQ1Fu3btAABz5sypkOupSOr87GnruGxRl4GzszMAICoqSml7VFQUnJyc1K6jznF1mUQigaOjY5mupyx1CgoK4O/vj+vXryMwMBCOjo4VcBWVw9nZuchrBlDi/6nS6jg5OeHChQuYP38+5s+fj3379kEmk2H+/PlV8jUtZ2dnpKWlKXVL5+XlIT4+vsT7VFodFxcXAMCkSZMUfwx2794d7u7uOH/+fAVfleZV1O+ohw8f4t69e5g6darij+KGDRuiX79+OHz4cIVcS0Ur7ufIwsICtWrV0qnjMlGXQaNGjeDq6qr0HzMyMhLXrl1D7969Fdtu376N27dvq1xH1eNWJb1790ZAQABeTCUvhMCRI0eUricqKkopaahSRyaTYdSoUbhy5QrOnj2r+EVbVfXu3RtXr15Vetfy0KFDcHV1RcOGDQEAmZmZCAoKQmpqqsp1goKClL7mzZsHfX19BAUFYfjw4ZV+neXVo0cPGBgY4MiRI4ptJ0+eRF5eHnr16qXYdunSJTx58kTlOp07d4aZmZnSL9YXydzOzq6Srk5zbG1t4eXlpfS7JDk5GefOnVP6Obp//z7+/fdfleu8uBeRkZFK54uIiKiS9wnPf46OHz+OgoICxbZDhw4p/X96+vQpgoKCVH5nX9Xjlpnaw9BqqD179ggDAwOxbNkysX//ftGmTRvRtm1bpXfk+vbtK/r27VumOqqUqUpCQ0OFtbW1GDt2rDh8+LAYM2aMsLa2FqGhoYoy33zzjdJIZFXqjBkzRhgbG4sdO3aI8+fPK76q6mtsBQUFom3btqJt27Zi//79YtmyZcLAwED8+uuvijLBwcECgAgMDFS5zquq+qhvIYT48MMPhY2Njdi6davYsWOHcHBwEFOmTFEqY2Njo/QOtCp1Vq1aJerVqye2bdsmjh07JoYNGybs7e1FbGxspV2bJh07dkzo6+uLhQsXioMHD4quXbsKDw8PkZ2drSgzYsQI0b59+zLVGT16tKhbt67YtGmTOHHihJgxY4bQ09MTZ8+erfRr1ITY2FhRt25d8frrr4vDhw+L6dOnCxMTE6W3BLZt2yYAiOTkZCGEENnZ2YrfOXXr1hWTJ08W58+fV3qNTZXjlhVXz1LD8ePHsWXLFqSkpKB9+/b46KOPYGVlpdj/wQcfAABWrlypch1Vy1Ql9+7dw8qVK/Ho0SO4ubnhgw8+gLu7u2L/nj17sH79epw9e1blOoMGDSpyQIa/vz9mzJhRCVeleampqfj6669x6dIlWFtbY+LEiejfv79i/8OHDzF+/Hh8//33aNWqlUp1XhUQEICVK1cq3euqRi6XY/PmzTh06BDkcjkGDBiA6dOnK41PGDRoEIYOHYrJkyerXAcAfvvtN+zevRvZ2dnw9PTE3Llzq+wgKQAIDAzE+vXrkZiYiFatWmHevHmwtbVV7F+0aBFiYmKwYcMGlesUFBRg69at+PPPP5GSkgJXV1dMnToV3t7elX59mvL48WOsWLEC9+/fR4MGDTBnzhzFzxie/07+8ssvceLECZibmyMqKgojRowodJy2bdvi22+/Vfm4ZcVETUREpMP4jJqIiEiHMVETERHpMCZqIiIiHcZETUREpMOYqImIiHQYEzUREZEOY6ImIiLSYUzURDVAXFwc9uzZo5VzX79+HUFBQeU+zr///lvqhC137tzB6dOny30uIl3CRE1UzSQkJGDPnj2Qy+WKbTdv3oS/v79W4tm6davSLH3q2r17N5YuXVpimf3792PhwoXlPheRLmGiJqpm7t27B39/f+Tl5Wk7FCLSgKq5iC8RFSk9PV3R9fvbb7/B0NAQrq6uSmVu376NR48ewd3dHY0bN1Zsj4uLQ2BgIEaOHIkrV67g8ePH6NGjB2xsbIDn3coPHjxA/fr10apVq0LzZUdGRuLGjRuoVasWvL29YWxsrLQ/OzsbISEhSE1NRdu2bRXHfUEul+Py5cuIjo6Gm5sbvLy8Sr3e/Px8nDt3DkKIcs2lTKTLmKiJqpH09HTFOspHjhyBnp4efHx8FAl5yJAhiImJga2tLc6cOYMVK1Zgzpw5wEvd4zt37kRMTAwaNWoELy8vmJqaYtSoUbh06RLatGmD0NBQxfKRL9Yr/vbbb/H555+jc+fOyM3NRXx8PPbu3QsPDw/g+cIibdq0QYMGDZCQkIDw8HCcPn1akYyTkpIwYMAAREREoGXLlrh06RI6deqEffv2wcjIqMhrjY+PR8+ePZGamopmzZrh5s2bcHNzq5T7TFSpyrVOGBHpnPPnzwsASksUnjp1SgAQS5cuVWzbvHmzMDU1Ffn5+Upl5s2bp3S8uXPnCh8fH8Xx5HK5GDNmjBg6dKjis6Wlpdi/f7+iTlhYmLh+/boQQogZM2YIQ0NDxWchhBg8eLB46623FJ9nzpwpmjVrJlJTU4UQQkRERAg7OzuxatUqRZl58+YJX19fxefp06cLb29vkZGRIYQQ4v79+8LExER07txZA3eRSHfwGTVRDTJ9+nTFv7t3746srCxERkYqlZk1a5bi30II/PTTT/D09ERAQAB+++03/Pbbb3BwcEBgYCAAQCKRwNjYGLdu3YJMJgMAuLi4KHVFd+jQQemzj48P7t+/r/i8Z88ezJgxA5aWlgCABg0aYNy4cSWOVN+7dy/effddmJmZAQCaNGkCPz+/ct0fIl3ERE1Ug9SuXVvxb6lUCgDIyclRKvPyOszp6elISkrCrVu38Pvvv2Pfvn3Yv38/Hj9+jP79+yM/Px8A8PPPP2Pnzp2oU6cO/Pz8sHv3bry8gu7L531x7hfnzcrKQkJCQqFu64YNG+Lx48dFXkdWVhYSExPh4uKitP3V5/FE1QGfUROREolEovi3iYkJDAwMMGbMGLzzzjvF1unbty/u37+PR48e4cSJE5g9ezYePnyIzz//vNTzmZqawtTUFElJSUrbk5KSYGtrW2QdExMTmJiYIDk5WWn7q5+JqgO2qImqGXNzc6CIlrI6DA0N4evri40bNyq9lw0AUVFRAIC8vDzEx8cDANzc3PDuu+/irbfewsWLF1U+T+fOnbF//37FZyEE9u3bhy5duhRZXiKRoHPnzjh48KBiW15eHgICAsp8jUS6ji1qomqmUaNGsLKywieffIIuXbqUeyT0mjVr4OPjg65du8Lf3x8ymQyBgYGoXbs2tm7ditzcXHTs2BEDBgxAy5YtERcXh507d2L16tUqn+Prr79G586dMXbsWHTt2hUHDx5EREQEDh8+XGydpUuXolu3bjA2Nka7du2we/duZGVlletaiXQRW9RE1Yy5uTlOnz4NQ0NDHDlyBNevX0fdunUxYsQIpXKmpqYYMWKEYgBXUWXwfJDW7du38cYbb+DatWt4/PgxJk6ciK1btwIALCwsEBwcjCZNmuCff/5BQkICDh06hAkTJgAAWrduja5duxY65sCBAxWfvby8EBISgnr16uH8+fNo27Ytbty4gfr16yvKtGzZEj169FB8bt++Pf755x8YGxvj33//xfTp07Fx40b06tVLY/eSSBdIxMsjPoiIiEinsEVNRESkw5ioiYiIdBgTNRERkQ5joiYiItJhTNREREQ6jImaiIhIhzFRExER6TAmaiIiIh3GRE1ERKTDmKiJiIh0GBM1ERGRDmOiJiIi0mH/DxzAR0JoDJKFAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x320 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "target = ~is_control\n",
    "alphas = np.linspace(0.001, 0.1, 50)\n",
    "pv, qv = p.to_numpy(), q.to_numpy()\n",
    "plt.figure(figsize=(5, 3.2))\n",
    "plt.plot(\n",
    "    alphas, [int(np.nansum(pv[target] <= a)) for a in alphas], color=\"#06c\", label=\"empirical p (control-gene rate)\"\n",
    ")\n",
    "plt.plot(alphas, [int(np.nansum(qv[target] <= a)) for a in alphas], color=\"#c30\", label=\"target-decoy q (FDR)\")\n",
    "plt.xlabel(\"threshold\")\n",
    "plt.ylabel(\"genes called\")\n",
    "plt.legend(frameon=False, fontsize=8)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "calib-intro",
   "metadata": {},
   "source": [
    "## Checking the null\n",
    "\n",
    "If the p-value threshold is set honestly, holding some control genes out and setting the cutoff on the rest should call the held-out ones at about the nominal rate. Split them in half, set the cutoff on one half, and measure the rate on the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "calib",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:10.634411Z",
     "iopub.status.busy": "2026-10-21T00:50:10.634268Z",
     "iopub.status.idle": "2026-10-21T00:50:10.687880Z",
     "shell.execute_reply": "2026-10-21T00:50:10.686309Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nominal 0.01: realized 0.0214\n",
      "nominal 0.05: realized 0.0571\n"
     ]
    }
   ],
   "source": [
    "rng = np.random.default_rng(0)\n",
    "order = rng.permutation(np.where(is_control)[0])\n",
    "fit, test = order[: len(order) // 2], order[len(order) // 2 :]\n",
    "for a in (0.01, 0.05):\n",
    "    cutoff = np.quantile(score[fit], 1 - a)\n",
    "    print(f\"nominal {a}: realized {np.mean(score[test] >= cutoff):.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "calib-fig",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:10.691245Z",
     "iopub.status.busy": "2026-10-21T00:50:10.690864Z",
     "iopub.status.idle": "2026-10-21T00:50:10.774970Z",
     "shell.execute_reply": "2026-10-21T00:50:10.774080Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV4AAAFKCAYAAABVbDeXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVmRJREFUeJzt3Xdck9f+B/BPwgh7yVCmIshQRK6KIDhQ3NraXq3t1Vtt1eqtdfbaasfPUStW21tb6621tdZq6/aq1WoVNxVQlOHejMgSCGEHkuf8/lBSw0weQsL4vl+vvF7mJOfk+wB+eTjPeb5HwBhjIIQQojNCfQdACCHtDSVeQgjRMUq8hBCiY5R4CSFExyjxEkKIjlHiJYQQHaPESwghOkaJlxBCdMxQ3wG0VBzHITMzE5aWlhAIBPoOhxDSAjHGUFxcDGdnZwiF6p/HUuKtR2ZmJtzc3PQdBiGkFcjIyICrq6va76fEWw9LS0vg2RfUyspK3+EQQlqgoqIiuLm5KfOFuijx1qN6esHKyooSLyGkQZpOR9LFNUII0TFKvIQQomOUeAkhRMco8RJCiI5R4iWEEB2jxEsIITpGiZcQ0m6IJZU4c6cYYkklr9e1hdbxEkLahS0xeXhrRxo4BggFwOYpHpgebq/269okoM0u61ZUVARra2tIpVK6gYKQVk4sqYTH0mvgamQ7axMhhEIBOI5BWsGpvGYgBFJXB8DV1rjecfnmCTrjJYS0efdyZbWSLoBayfZ5Cg64nytrMPHyRYmXENLmudoY1WoTCoBTC7uho7URsqVVGPrlXZXkbCAEvBxFzRIPXVwjhLR5Oy8XqDw3ED6dwx3sYwnfjiYY7GOJzVM8YCD86/XvJns0y9ku6IyXENLWpeXLEHU8GwCw8TU3+HcyhZejqFZSnR5ujxHdrXA/V1bn69pEiZcQ0qa9u0+MiiqGwd0s8K9BDg1WEnO1NW7WhFuNphoIIW3WqVtF2H+1EAZC4OtX3VvMbjKUeAkhbVKVgmHu7gwAwNuDHBDgYqrvkJQo8RJC2qRvzuTiVlYF7C0MsWKcs77DUUGJlxDS5uQUVWH5b5kAgKiXXGBr3rIuZ1HiJYS0OUsOPEZRBYc+HmZ4s38HfYdTCyVeQkibEvewBD/F5gMAvnnNHUJhy7ig9jy9n3+XlZXh6NGjyMnJQUBAAAYNGtRoH4VCgWPHjuH27duYNGlSnduw3759G7GxsTA0NERoaCi8vLya6QgIIboillTiXq4M3vWss03Pl2HqT6kAgGmhHdCvi7keomycXhNvVlYWBg4cCHNzcwQFBWHlypUYPnw4duzYUW+fffv2YfHixXBzc8OFCxfQp08flcTLGMMLL7yAhw8fol+/figrK8Nbb72Fjz76CB9++KGOjowQom3PVw8TAHgpyAZ9O/+VWC+nluJ/iYWovuu3p2vLWcVQk14T7/vvvw9LS0vExsZCJBLh+vXrCAwMxIQJEzB+/Pg6+zg4OODs2bMwMDCo80yXMYZZs2Zh7Nixyra9e/filVdewcSJE9GtW7dmPSZCiPaJJZXKpAsADMCBxEIcSCyst8/i/WJM7G2rkxsiNKW3OV6FQoEDBw5g6tSpEImeFqLo0aMHwsPDsWfPnnr7DRo0CB4eHvW+LhQKVZJudR8AuH//vtbiJ4ToTn3VxUb1sMKbYR0wqnvtkozV1cVaIr2d8WZkZKC0tBQ+Pj4q7b6+vrh06ZJWP+vAgQMwNDREUFBQve+RyWSQyf76JhUVFWk1BkIIf8YGtS+QVRe6cbU1rrPebnNWF2sqvZ3xFhcXAwBsbGxU2m1sbJSvaUNycjIWL16MDz74AJ06dar3fVFRUbC2tlY+6prGIITox9oT2SrPa1YPc7U11ml1sabS2xmvmZkZUMeZpVQqhbm5dq5E3rp1C8OHD8ekSZOwfPnyBt+7dOlSLFq0SPm8qKiIki8hLcDx61IcTpbCUAicmO8NgUCg9+piTaW3xOvh4QGRSIQHDx6otD948ADe3t5NHv/OnTsYMmQIxowZg82bNzdaHEMkEinnmgkhLUOlnMP8PU/rLcwb4ogI34a319FVdbGm0ttUg6GhIcaMGYNffvkFHPd0+43U1FScO3cOL730kvJ9J06cwE8//aTR2Hfv3kVERARGjRqFH374AUIh3SdCSGu0/lQu7ubI4GRliGVjW1a9habQ63KytWvXon///hg+fDiCg4Oxa9cuRERE4NVXX1W+Z8+ePYiLi8O0adMAADdu3MCxY8cglUoBALt370ZCQgL69++P/v37o7y8HEOGDIFCoYCvry/+85//KMcaMWIEAgIC9HCkhBBNZRZW4pOjWQCAz152hZWpgb5D0hq9Jt6uXbvi+vXr2LlzJ3JycrB69WpMnDgRBgZ/fYFHjBgBPz8/5fPy8nJkZz+daH/33XcBANnZ2SgpKQGereOtTty5ubkqn1deXq6T4yKENN17+x+jRMYhpIs5/tnPTt/haBVt714P2t6dEP2JuV+CAevuQCAALi3xRZ/OLfPWX755giY/CSEtioJjmLsrHQAwPcy+xSbdptB7kRxCCKkmllTii5M5SMooh42ZAVaPbzsX1J5HiZcQ0iI8XwQHAEZ1t4KDpZG+w2oWvBIvYww3btyAWCwGALi5ucHf37/FbCRHCGldahbBAYA9VyRY+/fKVrEuV1MaJd7bt2/j66+/xu7du1FQUKDyWocOHfDqq69i7ty5teovEEJIQ+oqglNd5KYtJl61L67NnTsXffv2RWFhITZs2IC7d+9CKpVCKpXi7t27+Oqrr5Cfn48+ffpg3rx5zRs1IaRN8XYUoebfyy25yE1TqX3Ga2xsjEePHsHe3r7Wa1ZWVvD29sbkyZORl5eHqKgobcdJCGnDXG2N0dneGI/yKoFWUOSmqWgdbz1oHS8hunPtcTl6rrwJAwGwZ5Yngjubt4qkq9N1vFVVVbhw4YLy+eXLlzFjxgxERUVBLpfzGZIQ0o59f+EJAGB8kA1eDmqZu0ZoE69VDZ9++ikMDQ0xYMAAlJSUYPTo0ejVqxdOnTqFoqIimmoghKitvJLD9vinF+tnhNWeymyLeJ3xbtu2TVm05uTJk+jcuTNOnjyJ3377Db/++qu2YySEtGEHEiUoLFPA3c4Yw/zbx7Qer8Sbm5ur3DnizJkzGDNmDPCs6E1eXp52IySEtGnfxzzNGdPDOsBA2D7uBeCVeLt3746vv/4aKSkp2LVrF0aOHAkAuHnzJvz9/bUdIyGkjbqbU4Fzd0sgFABv9G8f0wzgm3g/++wzrFu3DoGBgRg+fDhCQkIAAOvXr8f8+fO1HSMhpI3a8uxsd2R3K7jZte0Las/jdXEtIiICeXl5kEqlsLP7q07m0qVLVWrnEkJIfSrlHH6KzQcAzAhvP2e7aEqRHAMDA5WkC4CmGQghajuSIkVusRxOVoYY29NGjR5th9qJd8qUKWoPumPHDr7xEELaieqLam/0t4eRQfu4qFZN7cRraEgVJAkh2pGWL8MfN4uAZ6sZ2hu1s6mmO/0SQkh9tl7MB2NAhI8lvBxN9B2OztHWP4QQnVJwDD/++XSaYWY7u6hWjXfivXbtGqZNm4Z+/fohODgY06ZNw/Xr17UbHSGkzfnjRhEyJFWwMzfAS0Ht66JaNV6J98iRIwgKCkJqaioiIiIwdOhQpKamolevXjhy5Ij2oySEtAliSSWijmUBAP7ZrwNMjNrnH928rph9+OGHWLt2LRYtWqTS/p///Acffvghxo4dq634CCFtRM091ezMDfQdkt7wqsdrbGyMJ0+ewNraWqVdKpXC0dERMplMmzHqBdXjJUR7xJJKeCy9prK9j4EQSF0d0KpLQOq0Hq+DgwOSk5NrtSclJcHBwYHPkISQNqyhPdXaI15TDdOnT8ekSZPw4YcfIjg4GAAQHx+PVatWYdasWdqOkRDSynk7iiAUoNYZb1vdU60xvBLvsmXLYGJigo8//hiFhYUAABsbGyxevBjvv/++RmNdvHgRGzduRE5ODgICArBkyRI4OTk12CchIQGbNm3C7du3sWHDBgQFBWllXEJI83C1NUYfDzNcSi0D2sGeao3hNdWQkZGBDz74AAUFBRCLxXj8+DEKCgrwwQcfwMBA/Qnz8+fPY/DgwXBzc8PcuXNx/fp1hIWFoaSkpN4+//d//4fZs2fD09MTf/75J6RSqVbGJYQ0n/wSOZLE5QCATf9wQ+rqAExvp2t4AQCMB6FQyKdbLeHh4WzixInK5yUlJczS0pJ98cUX9fYpKChgjDGWkZHBALAzZ85oZdyapFIpA8CkUqkGR0QIqcv66GyGtxJYr09uMI7j9B2O1vDNE7zOeJ2dnZGRkdGkhF9WVoaLFy/ihRdeULaZm5sjMjISJ0+erLefra1ts4xLCGkejDFlQZwZYfYQCNpXQZy68Eq8ixcvxpw5c5Cdnc37gzMyMsBxHFxcXFTaXVxckJ6ervNxZTIZioqKVB6EkKaLf1SKG5kVMDUSYHI/OzV66AdjDBkZGeCxwlZjvC6urV69Gjk5OejUqRNsbW1hbKw6Qa5OQq6qqgIAmJioFsgwNTVFZWUln7CaNG5UVBRWrFjB+3MJIXWrPtud2NsWNmYts8phZWUl4uPjkZ2djYiICNjbN+/8M6+vwueff97kD64uop6fn6/Snp+fX6vAui7GXbp0qcqdeEVFRXBzc+MdByEEKCpXYNdlCdCCd5mQSCS4ePEi5HI5BgwY0OxJF3wT7/Xr17FmzZo6X1uyZIlaYzg7O8PJyQlXrlxRucX40qVLCAsL4xNWk8YViUQQidrnmkJCmsuuywUoq+Tg29EE4V4W+g6nFolEgtOnT8Pa2hqDBw+Gubm5Tj6X92aXfF6r6Y033sAPP/yArKynRTMOHjyI69ev44033lAZ7/XXX9coPnXGJYQ0P+VFtfCWdVGteh7XxsYGPXv2REREhM6SLpqy51pdkpOTNbpleNmyZbhz5w68vLzg4eGB1NRUfPXVV+jXr5/yPffu3cPVq1eVz48dO4ZPP/1UOV87d+5cWFtb480338Sbb76p9riEkOaVlFGGhLQyGBkI8HpIy7moVlpairi4OAQEBMDR0RHe3t46j0GjxGtjY1PnvwGA4zgUFxdj4cKFao9nYmKCAwcOID09HTk5OejWrVutwjtLlixBcXGx8nnv3r3rnOZwd3fXaFxCSPP64dnZ7vheNnCwNNJ3OMCzC//x8fEwMjKCkZH+YtKoOtm+ffsAABMnTsTevXtVXjMyMkLnzp0RGBio/Sj1gKqTEcJfWSUH5/dSIC1X4MR8bwzz1+//IcYYbt68iZs3b6Jjx47o169frdVYfPDNExqd8U6YMAEAcObMGQwePFjzKAkh7cK+KxJIyxXo3MEYQ30t9R0O5HI50tLS0L17d/j5+el9vpnXHG910lUoFCrTANVqTkMQQtqX6mmG6WH2EAr1l+QkEglEIhHMzMwwfPjwFrNbOq9VDbdv38aAAQNgYmICW1vbWg9CSPt1O7sCF+6XQCgA3uivv63bHz16hNOnT+PmzZsA0GKSLvie8c6YMQO2trb4/fffKdES0o6IJZW4lyuDt6OozpKOYkkl/u9wJgBgTIA1XPRQ9lGhUCAxMRGPHj2Cp6cnevXqpfMYGsMr8SYmJiI9PR0dOujvtxkhRLee3zNNKAA2T/FQKe1Yc081DzvdJ13GGM6fP4+CggL06dMHXbp00XkM6uCVeN3c3FBRUaH9aAghLZJYUqmSVDkGzNyeht+vS2EuEqJUxuFAYqFKn2/PP8H7IzvqrNg5YwwCgQDe3t4wNzdv0X+N85rjnTt3LhYuXFhnEXJCSNtT155pDMCBxEJsjyuolXShwz3VGGO4ceMGEhMTAQCurq4tOumC7xnvJ598gpycHOzfvx/29rVvBWxKuUhCSMtT155pAgHwwciOsDEzQGGZAquPZeP53KyLPdWeryrWvXt35VlvS6e36mSEkNbD1dYYE/5miz1XnlYaq94z7fk53i72Isz6JQ0KTjd7qtWsKtaxY8dm+yxt45V4p0yZov1ICCEtFmMMKY+f7pm2MNIRiyKdaiXV6eH2GNHdCvdzZfCqZ9WDNqWmpkIkEum0qpi2NHlhW1ZWFhhjcHZ21k5EhJAW588HpbidXQEzYyGWj3WGlWndm9q62ho3a8JVKBSQSCSwt7dHz549AUCjDXZbCl4X1xQKBaKiomBrawtnZ2e4uLjA1tYWUVFRUCgU2o+SEKJX1XeivdrXtt6k29xKS0tx+vRpxMTEoKqqCgYGBq0y6YLvGe+KFSuwadMmLFu2DCEhIRAIBIiNjcXq1atRXl6OlStXaj9SQoheFJbJsSehAHi2WaU+ZGVlIT4+HsbGxhg8eLBeK4tpA6/Eu2XLFuzZs0elUE6/fv3Qq1cvTJkyhRIvIW3Ir5cKUF7F0N3ZBCGeup9LffToERISEtCpUycEBwdrpaqYvvFKvE+ePEFQUFCt9qCgIDx58kQbcRFCWojqaYaZetpFwsnJCT179kS3bt1axVIxdfCa4/Xz88O2bdtqtW/duhW+vr7aiIsQ0gJcSStFYkY5jA0FmNJPdyUCJBIJzp49i8rKSpiZmcHHx6fNJF3wPeNdtWoVXnrpJRw+fBjBwcEAgPj4eJw7dw4HDx7UdoyEED2p3jPt70E26GChm+peDx8+RGJiIqytrSGXy9vE1EJNvL6S48aNQ2JiIr744gucPHkSAoEA/v7+SExMREBAgPajJIToXKlMgV8vPb2oNnOA+nsp8qVQKHD16lWkpqYqq4q11lULjeH9KywgIAA//fSTdqMhhLQYexIkKK7g0NVBhEHezb81u0QigVgsRt++fdG5c+dm/zx94r2ONykpqVZ7UlISreMlpI14fmv25txFIj8/H4wx2NvbY/To0W0+6YJv4o2KisL+/ftrte/fvx+fffaZNuIihOjRjcxyxD4shYEQmBbaPBfVqquKnT59Gunp6QAAkah5i+q0FLwS76ZNm/Cvf/2rVvvs2bOxefNmbcRFCNGj6iVk43raoKO19m9WkMlkiImJwc2bN9G9e3e4u7tr/TNaMl5zvBKJpM6lHQKBALm5udqIixCiJxVVHH6Oyweerd3VtvLycpw+fbpVVhXTFl5nvP369cP69etrtX/55ZfK5WWEkNbpf4mFKChVwNXWCCO6W2l9fBMTE3h4eCAyMrJdJl3wPeP99NNPMWTIEFy4cAEDBw5U7nOUlJSE06dPaz9KQojOfHPm6V+tE4JsYaCli2rVG1A6OzvD2dkZPXr00Mq4rRWvM97Q0FBcunQJ3t7eOHr0KI4dO4Zu3brh0qVLCA0N1Xg8uVyOwsLaW4c0tQ/HcbQ9ESEaiDqWhYsPSwEAX5/JxZZnc71NUV1VLC0tDVVVVVqIsvXjlXjxbB3vtm3bcO3aNaSkpGDbtm0a3zzBGMP7778Pa2trdOrUCe7u7jh8+HCT+4jFYowdOxZmZmZwd3eHjY0N3n//fXAcx+tYCWkP4h6W4IODmcrnHANm/ZIGsaSS95hZWVk4efIkqqqqMHToUHh4eGgp2taNd+LVhi+//BKbN2/GuXPnUFJSggULFmDChAm4c+dOk/rMnDkTT548gVgshlQqxbFjx7BhwwZ8//33OjoyQlqX+EelGL3hXq32pmxYyXEcUlJSYG9vj8jISNjY2Ggh0rZBr4l3w4YNmDFjBvr06QMDAwMsWrQIrq6u+O6775rU586dOxg1ahTs7Z9ekQ0NDYW3t3eDCZ2Q9mpPQgEGf3EHkrLafxHy2bBSJpOhtLQUQqEQgwYNQlhYWJust9AUeku8T548QWpqKsLDw1XaBwwYgEuXLjWpz5w5c7B9+3acO3cOqamp2Lx5M8RiMaZOndpMR0NI68MYw6qjWZj0/SNUVDGMDbDGN6+6weBZVuCzYaVEIkF0dDQSEhKAZysY2lJVMW3htaphyZIlWLNmjcavPa+6bm/1WWk1BwcHxMbGNqnP22+/jYSEBERERMDCwgIymQwbNmxAYGBgvfHIZDLIZH/9SVVUVNToMRDSEokllbiXK4N3PRtOiiWVuJFZjs0X8nAg8ekF6oVDHbFugisMhAK82MuG14aVz1cV69u3r1aPqa3hlXg/++yzepNrQ689r/q3oFwuV2mXy+X1ViRSt8/LL7+M4uJiZGdnw9HREXFxcRg+fDgEAgFmzpxZ59hRUVFYsWJFo3ET0pJticnDWzvSwDFAKAA2T1Hdgv351wFAIAD++5o7Zg/6q/oYnw0rExMTcf/+/TZfVUxbtFpgMzk5GQ4O6pWPq96VOCcnR6U9Jyen3h2L1emTnZ2N48eP49ChQ3B0dAQAhISEYMKECfjxxx/rTbxLly7FokWLlM+Liorg5uam1rEQ0hKIJZUqSZVjwIztadhxKR+mRkKUV3E4e6dEpY8AwNie1k3+bAcHB9ja2raLAjfaoNEcr42NjfLKZPW/qx9WVlbKPdfUYW1tjcDAQJw8eVLZJpfLcerUKQwcOFDZVlRUhPz8fLX7mJubQyAQoLS0VOXzSkpKYG5e/35RIpEIVlZWKg9CWpODiYXKpPu8s3dKcOx6Ua2ki2fJme+qhaysLCQmJoIxBldXV0q6GtDojPeHH34AAEycOFH572pGRkbo3Llzg/OoNX300Uf4xz/+gZCQEISGhuLzzz8Hx3EqBXgWLVqEuLg4XL9+Xa0+lpaWeOGFF/DRRx/BwcEBnp6eOHHiBA4cOIAtW7ZocriEtBo7LxXg3X0ZtdqFAuDzCa6wMzdEQakc7+4Tgz2XnPmsWqiuKnbr1i106tQJHMfR1IKGNEq8EyZMAACcOXNGZYdhviZMmACZTIavvvoKK1euREBAAM6ePaucIsCzs9znL6ap02f79u1Ys2YNlixZgvz8fHh4eGDbtm2YPHlyk2MmpCVhjGHlkSwsP5IFAOjlaoprj8uhYH+tSnh+jtfKxACzfkmDguO3akEmkyE+Ph45OTno0aMHfH19adUCDwLGWB1/nJCioiJYW1tDKpXStANpkSqqOLy5LRU7L0sAAP8e5oQ1L7sgS1rV4KoEsaSS16oFALh58ybu3buHkJAQODk5ae1YWiu+eYJX4jU0bPhEueaqg9aIEi/RF3WWg116VIpPj2Xhano5DIXAt5M9MKMZSjji2Vl1cXExrKyswHEcZDIZTE1Nm+WzWhu+eYLXqoYjR46oPOc4Dvfu3cOaNWuwcOFCPkMSQngsBzM1FuC3t70w1K95Tg6qN6BMS0vDyJEjYWFhQUlXC7Q61XD27FmsXLmyTZSGpDNeomtiSSU8ll6rtTJBKHi63pYx1PlaWlSAxlMG6igpKUFsbCyKiorQu3dvWrVQB52e8danb9++uHLlijaHJKTduJcrq3M5GMcA1HN6VL0cTNuJNy8vDzExMTA2NsbQoUOpwI2WaTXx7t69G3Z2dtockpB2I0dau/yigQCIW+ILZxsjZBZWod+a2yrJmc9yMHVYWFjA1dUVPXv2pAI3zYBX4q2rerxEIkFOTg6VXiSEhyoFwye/ZwPP7iZjzxWp6dP56Y0/zjbG2DzFo0nLwRoik8mQnJyMnj17wsTEBH369NHKuKQ2Xol3xowZtdpsbW3Rv39/eHt7ayMuQtqVjWdycTOrAvYWhjizyBt5JYo6l3tND7fHiO5WvJeD1aegoACxsbGQy+Xo2rUrTExMtDIuqRut460HXVwjupJTVIVuH19HUQWHzVPcMXOAevVOtIExhkePHimrivXv3x9mZmY6+/zWTm8X17KyssAYq7ewDSGkYUv/9xhFFRx6u5vhzbDmWYtbn9LSUly9ehWenp4IDAykW391hFchdIVCgaioKNja2sLZ2RkuLi6wtbVFVFQUFAqF9qMkpI2Kf1SKrRefFoH65jU3re3q25iysjJwHAcLCwuMGDECf/vb3yjp6hCvM94VK1Zg06ZNWLZsGUJCQiAQCBAbG4vVq1ejvLwcK1eu1H6khLQxHMcwd1c6AGBqaAeEeFro5HOzsrIQHx8PHx8f+Pn5wdLSUiefS/7Ca47XxcUFv/zyS61COWfPnsWUKVMgFou1GaNe0BwvaW5bYvIwY3saLE2EuLuyBzpaGzXr59WsKhYcHExLxZpIp3O8T548QVBQUK32oKAg5fY8hJD6FZbJsfTgYwDA8rHOzZ50FQoF/vzzT6oq1kLwmuP18/PDtm3barVv3boVvr6+2oiLkDZt2W9ZeFIsh18nE8wd4qhGj6YRCoWwsrLCwIED4efnR0lXz3id8a5atQovvfQSDh8+jODgYABAfHw8zp07h4MHD2o7RkLalOhbRfjmTC4A4OtJbjAyaJ4kyBjDw4cPYWxsDDc3N/Tq1atZPodojtcZ77hx45CYmAhXV1ecPHkS0dHRcHNzQ2JiIsaOHav9KAlpI36IeYJh6+8pb/tNy699m7A2KBQKXL58GVevXkVBQUGzfAbhj26gqAddXCPaJpZUwn3ptVpb76Su1m51seqqYsXFxejduzc8PDy0NjZRxTdP8DrjrWn58uXaGIaQNi1FXIaapzkKjv9mk/VJSEhAVVUVhgwZQkm3hdJK4l2xYoU2hiGkTTt+vahWm7aqizHGUF5eDgAIDg5GZGQklXJswbSSeAkhDXvwRIbvYvKAZ8XLAe1VF5PJZLhw4QLOnTsHjuNgZmZG63NbOK3W4yWE1G3hngxUyhmG+Vliy+seePCkUivVxZ6vKhYSEgKhkM6lWgOtJN7ExERtDENIm3TsuhS/pUhhKAS+ftUdbnYiuNk1fXohNTUVV65cgY2NDUJDQ6mqWCuilcRL6wMJqZusisP83RkAgPlDneDbUXt1bkUiEbp06UJVxVohtRNvZGSk2oNGR0fzjYeQNmX9qVzcy5Who5Uh/m9MpyaPV1JSgocPHyIgIACdOnVCp05NH5PontoTQiEhIcqHp6cnTp06hZKSEvj5+cHPzw8lJSU4deoUPD09mzdiQlqJx5JKfPJ7FgDgs5ddYWXatLPSzMxMREdHQywWQybT7hI0oltqn/GuWrVK+e8JEyZg7dq1WLx4scp71q1bh0uXLmk3QkJaqcX7xSiVcQj1NMeUfvw3gX2+qpizszP69u1LqxZaOV53rjk4OOD+/fuwtrZWaZdKpfDy8tKoQlllZSVOnTqFnJwcBAQEoHfv3lrr8/jxY1y4cAFmZmaIjIzU6OID3blGmuL83WIM+uIuBAIg4QM//M2d/4Wv1NRUXL58maqKtUA6vXNNoVDUuZLh6tWr4DhO7XGePHmC3r17Y8GCBTh69CgiIyMxe/ZsrfT59NNP0a1bN+zcuRM7d+5E//79kZaWpnZshPCVmifDG9tSAQAzw+15J93q6QQPDw8MHTqUqoq1IbxWNcycORMTJ07E4sWL0bdvXzDGkJCQgHXr1uGtt95Se5ylS5cCz5ajmZmZ4cqVK+jbty/GjRuHMWPG8O6ze/duLF++HGfOnEF4eDjw7KyBylKQ5rYlJg8zt6eh+ifNr5PmqxiqN6BMSkrCwIEDYW9vDzs7/lMVpAViPCgUCvbll18yd3d3BoABYO7u7mz9+vVMoVCoPYalpSX74osvVNrDwsLYlClTmtTnb3/7G3v11Vf5HJqSVCplAJhUKm3SOKT9yCiQMeGsBIa3/noYzE5gGQUytceQy+UsPj6e7dmzh125coXJ5fJmjZk0Dd88weuMVygUYsGCBViwYAFKS0sBAObm5hqNkZGRgeLiYvj7+6u0+/v7IyEhgXefsrIyJCYm4l//+hdSUlKQmJgIZ2dnDBw4ECJR/YvWZTKZypXioqLa99UT0pDvzj9RlnusVl0ER5071MrKyhATE4OSkhIEBwdTgZs2rMn3F5qbm2ucdPFcYqtZyMPW1rbepKdOn4KCAjDGsG/fPkyaNAmnTp3C/Pnz4efnh/v379cbT1RUFKytrZUPNzc3jY+JtE+MMXxyNAurfs+u9ZomRXAMDQ1hYmJCVcXaAbXPeDUpcH7kyJFG32NqagoAKC4uVmkvKiqqd/WBOn2q35Obm4uUlBQYGRlBoVAgPDwc7777Lg4dOlTn2EuXLsWiRYtUxqTkSxojq+IwY3sadsQ/LTY+zM8Sp+8UQ8GpVwSHMYZbt26hc+fOMDMzw8CBA3UYPdEXtROvl5eXVj/Yw8MDxsbGePTokUr7o0eP6v0sdfp06NABdnZ2iIyMhJHR0w0EDQwMEBkZiR07dtQbj0gkanAqgpCanhRX4aVvH+DPB6UwEAL/fc0dbw10gFhSifu5skaL4MhkMsTHxyMnJwfm5uZ0ltuOqJ14169fr9UPNjIywsiRI7Fr1y7MnDkTAoEAmZmZOHPmDL777jvl+86dO4ecnBy88soravcZP348kpKSVD4vOTmZ7qojahFLKnEvVwbvehKnWFKJ6FtFWHY4E+mSKlibGmDfLE9E+j1dx+lqa9zonG5BQQEuXrwIhUKBgQMHwsnJqdmOh7Q8et365/bt2+jfvz/Cw8MREhKCn3/+GY6Ojjh9+jQMDZ/+TpgxYwbi4uJw/fp1tfs8fvwYoaGh6Nu3L8LCwhAfH4+jR4/i9OnTys05G0M3ULRPW2Ly8NaONHAMEAiAWQPsMcT3r+//6dtF+O58nnK5mL2FAc7/2wd+nUzV/oyqqiocPXoUlpaWVFWsleObJ3glXoVCgXXr1mHLli14+PAhFAoFAGD+/PlYsGABunTpovZYYrEYP//8s/IutNdff13ldsgdO3bg0aNH+Pjjj9Xug2dnFNu2bUNaWhrc3d3x2muvaVRQhBJv+1PXnmiNMRAAqVHq7Zkml8shEAhgYGCA/Px82NjYUFWxVk6niXf16tX46aef8NFHH2Hq1KnKGxN++eUXREdHY+vWrZoO2eJQ4m1fGGOYujUV2+Nr78gb6GoCGzNDSMrkSBFX1Hr9zKJuGOxj2eD4JSUluHjxIhwcHBAUFKTV2In+6DTxdu3aFfv27UNQUBAEAoEy8WZmZiIgIAD5+fmaDtniUOJtPyqqOLyxLRW7Lktqvfb8LsBiSSU8ll5TWaurzi7BmZmZuHTpEkQiEUJDQ2kvtDZEp7UaxGIxfH19AUDl3nGRSKS8oYKQ1iCnqAoRX9zFrssSGAqBqSF2MHj2v6LmcjBXW2NsnuJR7+s1McZw/fp1/Pnnn3BwcMDQoUMp6RKAb60GLy8vxMfHY/DgwSqJd+fOnejRo4c24yNE66pXLSg4hhnb05CWXwkbMwMcmN0VET6WWDXepd7lYNPD7TGiu5Vay8UEAgGqqqoQEBAAHx8fKnBDlHgl3vfeew9TpkzB8uXLAQBHjx7F8ePH8d133zW4VpYQfXt+1UK1rg4iHH3HCz7PtuVpbDlYY68XFBSgtLQUbm5u6NWrFyVcUguvxDt16lRwHIdVq1aB4ziMHTsW7u7u2Lx5M1555RXtR0mIFogllbWSLgD8b7anMuk2BWMMDx8+RFJSEjp06ABXV1dKuqROvDe7fOONN/DGG29AIpGA4zh06NBBu5ERomW3sypqJV0AyC9VNHlsuVyOq1evIi0tDV27dkVgYCAlXVKvJu8ybGtrq51ICGlGReUKfHo8q1a7JkVsGpKUlASxWExVxYhaeCXe4uJirF27Fn/++ScKCmqve6x5uy4h+pSWL8O4jQ9w7XE5jAyelmrkmHpFbBojl8thaGgIf39/eHl50aoFohZeiXfWrFmIjY3FP/7xDzrjJS1a/KNSvPjf+8gpkqOjlSEOz/FCJ2sjtVYlNKR6A8r09HQMGzYMZmZmdOsvURuvxPv7778jPj4ePj4+2o+IkEaoU8TmXq4Mt7PKsWifGBVVDIGupvhtjhfc7P5ak8uXTCZDXFwccnNzERAQoKwRQoi6eP3EWFpaomPHjtqPhpBGPL8cTCgANk/xwPRw+zpfrzY2wBq/zugCS5Om10WorirGcRwGDRoER0fHJo9J2h9etwzPnTsXrq6ueP/995snqhaAbhlueeq6ZRcAXGyMYGgggFzB8LiwSuU1AYBHq3vAo4N2ai3n5ubi+vXrCAkJoakFwjtPqH3G++qrryr/XVZWhm+++QZ79uyBl5dXrWUzu3btUjsAQtSVmFFW53Kwmsn2eQzAo7zKJiVeuVyOBw8eoFu3bnB0dERERAQtFSNNonbitbCwUPn39OnTmysmQmp5lCfDu3vFtdqFAuDQ213hZGWEnKIqvPjfB7WK2DRluVh1VbGSkhI4OTnBxsaGki5pMrUT7w8//NC8kRBSj9gHJXjx2wd4UiyHtakBiisUKsvBxvb8awnX5ikemPVLmtp7njXk+apiQ4cOhbW1tRaPirRnaifehIQE9OnTR+vvJaQhOy8V4I1tqZDJGYLcnq5MYIBWitg0JC8vD3/++SecnZ0RHBys3L+PEG1Q++Ja586dERQUhNmzZ2Po0KG1ltBUVlbixIkT+P7775GcnIzU1NTmilkn6OIaf+ou92ro9bs5FTiSIsWXp3IBAC8GWmPHm11goYWVCQ1RKBQwMDAAYwyPHz+Gi4sLTS2QejX7xbWbN2/iiy++wD//+U+UlpaiZ8+ecHJyAmMM2dnZSElJgZWVFd555x26uNaO1Vzu9ckLzpjY+6+bbPZekeDjw5lqvV7t38OcsOZlFxgImzcBFhQUIDY2FkFBQXB2doarq2uzfh5pvzReTiaTyXDq1Cn8+eefyMjIgEAggKurK8LDwzF06NBae5+1VnTGq7n6lns1hVAApKm5pxlf1VXFEhMTYWtrSxtQErU1+xlvNZFIhNGjR2P06NGadiVt3L1cWZ1J19xYqFxnW1rJafQ6x57O5zZX4n2+qpiXlxcCAwMhFPLamIUQtdG9jkRrvBxqJ0cDIXB7ZfcG9yxr7HVtVA+rD2MMRUVFVFWM6BT9aidaE/eoTOW5pnuWabqnWVNkZmaiuLgYRkZGGDp0KCVdolO8bhluD2iOVzOlMgX8lt1AhqQKC4c64oVAm3qXc4kllQ0u92rs9aao3oDy9u3b8PHxQc+ePbU6PmlfdDbHS0hd1hzPRoakCh4djPHpeBeYGtf/x1RT9zTjq2ZVMaquR/SFEi9psodPZFh3IgcA8J8Jrg0mXX1hjOHs2bOQyWRUVYzondqJNzIyUu1Bo6Oj+cZDWqGFezMgkzNE+lnipaCWtQMDYwyMMQiFQvTq1QuWlpa0VIzondqnJiEhIcqHp6cnTp06hZKSEvj5+cHPzw8lJSU4deoUPD09NQogMTERs2bNwvjx4/Hxxx9DIpFotc/BgwcxcuRIfPfddxrFRdRz/LoUh5OlMBQCX09ya1F3ecnlcly+fBmXL18GYwxOTk6UdEmLoHbiXbVqlfJRUFCAtWvXIi4uDhs2bMCGDRsQFxeHtWvXqpU4q8XFxSE0NBRGRkaYOHEizp49i7CwMJSVlWmlT2pqKubOnYvk5GTcunVL7biIeirlHObvyQAAzBviCL9OpvoOSamkpASnT5+GWCxGx44dW9QvBELAeLC3t2eFhYW12gsLC5m9vb3a40RERLDx48er9Dc3N2dff/11k/tUVVWxkJAQtnXrVhYYGMjmz5+vdlyMMSaVShkAJpVKNerXnnx2PIvhrQTm9O8kJi2T6zscpcePH7P//e9/7Pfff6/z55QQbeGbJ3hdBVEoFEhMTKzVfvXqVXBc7TuT6lJeXo7z589j/PjxyjZra2tERkbi+PHjTe7z0UcfwdXVFdOmTdPgyIg6xJJK7L1SgBVHMgEAn73sCivT5i1eo4m8vDw4OjoiMjKSSjmSFonXqoaZM2di4sSJWLx4Mfr27QvGGBISErBu3Tq89dZbao2RkZEBhUIBNzc3lXZXV1ecOXOmSX1OnjyJX3/9VaNt5mUyGWQymfJ5UVGR2n3bk5p7mnnaG+Of/ez0HRZkMhny8vLg4uKCgIAAAKDpBdJi8Uq8UVFR6NSpE7788kukp6cDANzd3fHRRx9h7ty5ao1RWVkJADA1VZ0XNDMzU77Gp09OTg6mTp2K7du3w85O/YQQFRWFFStWqP3+9kgsqay1kWRqfiUypVXNWsSmMfn5+YiNjQUAODk50a6/pMXj9RMqFAqxYMECLFiwAKWlpQAAc3NzjcawsXm67KigoEClPT8/H7a2trz77Nu3D6WlpVi3bh3WrVsHAHj48CEKCwtx+/Zt/P7773UWQVm6dCkWLVqkfF5UVFTrzLq9q6sITnMXsWkIYwwPHjxAUlKSsqoYJV3SGjT5p1TThFvN1dUV9vb2SExMxJgxY5TtiYmJ6N27N+8+48aNQ9euXVX63bt3Dz179sSsWbPq/fNTJBJBJGq+YixtgbejCIJnG0hWa+4iNg25e/cuUlJSqKoYaXV4X1xbs2YNvL29YWDw10WV+fPn49GjR2qP8/rrr2PLli3Iy8sDAPzxxx9ITEzE1KlTle9Zv349Zs+erXYfd3d3jBw5UuVhaWkJDw8PjBw5kub9msDUSAgz47++fs1ZxKYh1eVFPDw8EBISgqCgIEq6pFXh9dP62Wef4ccff8THH3+ssoohODgYK1euVHuclStXwtvbG97e3ggODsb48eOxevVqhIeHK99z/fp1xMTEaNSHNI+PD2eitJLBt6MJohd4I3V1AKaH2+s0hsePH+OPP/5ARUUFTExMaDqItEq8qpN17doV+/btQ1BQEAQCgfIMJDMzEwEBAcjPz9dovJs3byInJwf+/v5wcnJSee3GjRsoLCxEWFiY2n1qunjxIuzs7ODr66t2TFSdTFVSRhl6f3oLHAPOvtsNg7pZ6vTzn68q5uLigr59+9IGlETv+OYJXolXJBKhsLAQpqamEAqFyrPe/Px8uLi4oKKiQtMhWxxKvH9hjGHg53cRc78Er/a1xc4Zmt0W3lTPVxXr2bMnunXrRlNGpEXgmyd4TTV4eXkhPj4eqLFWcufOnejRowefIUkL9uulAsTcL4GZsRDr/q77DSBLS0tRXFyMQYMGwcfHh5IuafV4rWp47733MGXKFCxfvhwAcPToURw/fhzfffcdduzYoe0YiR4VVyiweP9jAMBHozvq7EIaYwxisRguLi6ws7PDqFGjVC7kEtKa8Uq8U6dOBcdxWLVqFTiOw9ixY+Hu7o7NmzfjlVde0X6URG9WHc1ClrQKXo4iLIpseC5dW+RyOa5cuYL09HT0798fLi4ulHRJm9LkrX8kEgk4jkOHDh20F1ULQHO8wJ3sCgSsvIkqBcORd7wwJqD56x4UFxcjNjYWJSUl6NOnD9zd3Zv9MwnhS6dzvJ07d8b06dMhk8lga2urknQ7d+7MZ0jSwjDGsGBPBqoUDGMCrHWSdEtLSxEdHQ2FQoGhQ4dS0iVtFq/Em5aWhtOnT2PQoEHIysqq9Rpp/bZezMfxG0UwMgDWv9K8F9Sq/+gyMzNDz549qaoYafN43+5z+fJlmJiYoE+fPrh06ZJ2oyJ6tencE0z/+ekvULkCOHe3pNk+q6KiAufPn0d6ejoEAgG6du1K63NJm8c78drb2yM6Ohrjx4/HoEGDsH37du1GRvRCLKnE27+mK58zALN+SYNYUnfFuKbIz89HdHQ0pFIpTExMtD4+IS1Vk4rkGBoaYuPGjejVqxdmzJiB5ORk7UVG9CLmfglqXm1VcNqtQPZ8VTE7OzuEhobWKvVJSFumlRp6M2fOhJ+fH/7+979rYziiRzviC2q1absCGcdxePToEbp27UpVxUi7xCvxZmRk1GoLDw9HcnIyrl+/ro24iB6culWEo9ekEAAQCJ7W2tVmBbLi4mIwxmBlZYWIiAiqnUvaLV4/+a6udV/l7tixIzp27NjUmIgeVCkY5u1++gv17cEOWDKyI+7nyuDlKNJK0n38+DEuXboEBwcHhIeHU9Il7ZraP/1jx44FABw5ckT57/ocOXKk6ZERndp4Jhc3syrQwdwAK19whp25oVYSLsdxuH79Ou7cuaOsKkZIe6d24vXy8qrz36T1yymqwrLfnu4YvHq8C+zMtXc2Gh8fD7FYTFXFCHlOk28Zbqva0y3Db25LxdaL+fibuxkuLfWFgbDpyZExBoFAgOzsbAiFQjg6OmolVkJaEr55giba2rn4R6XYevFp4fpvXnVrctKtXir25MkThISE0Jw/IXXQeI5XHTTH2zpwHMPcXU9vlng9xA6hXS2aNN7zVcW8vLyUZ72EEFW85nhJ27D1Yj4up5bB0kSIz15uWj2G4uJiXLx4EaWlpejXrx8VuCGkAWon3vXr1zdvJESnbmSW4d19T5ePLRvrjI7WTauPIBaLwXEchg4dSgVuCGkEzfG2Q1ti8jBze5ry1mALEb87xziOQ15eHhwdHeHr6wtvb29an0uIGnj/L5FKpTh37hzS09Mhl8tVXluwYIE2YiPNQCypxFs70lTqMczZmY4xAdYardutqKhAXFwc8vPzMWrUKJiZmVHSJURNvP6nXL58GWPGjIGxsTEeP36Mrl27Ii0tDXK5HL6+vpR4W7AtMXngaiwg1LQITn5+PmJjY8FxHAYOHAgzM7PmCZaQNorX35jvvvsu5s2bB7FYDAC4f/8+xGIxIiMjqVBOC8UYw5rj2Vh+JKvWa5oUwcnKysKZM2dgZmaGYcOGwcHBoRmiJaRt43UDhZWVFdLT02FjYwMDAwOUl5fD2NgYDx48wKBBg5QJuTVrSzdQVMo5zP4lXbled4iPBc7dK4GC+6sIzvRwe/XGqqzE/fv34evrS1XFSLun0xsoiouLYWNjAwBwdHRUrtu0t7dHXl4enyFJMykolePlTQ9w7m4JhAJg/StumDvEEWJJpdpFcIqLi3HlyhX07dsX5ubm8Pf311n8hLRFTb4aMmjQILz//vv417/+he3btyMwMJDXOBzHaXwG1ViftryAXyypxL1cGbzrSZxiSSXO3C7G//2WidT8SliaCLFrhidGP9u00tXWWK053eqqYiYmJlAoFM1yLIS0N7z+VoyKilL+e+3atcjNzcXo0aNx5coVbNq0SeOxnJycYGRkhICAAJw+fbrJff744w8MGzYM1tbWsLa2xpgxY3Dr1i2N4mrJtsTkwWPpNQz5z114LL2Gb889QXklp3x8e+4JPJZew+s/pSI1vxJ25gb48z1fZdJVB8dxSElJwcWLF+Hk5ITIyMhWP+VCSEuhtSI5fM4uN23ahMWLF+N///sfQkJCsG7dOqxbtw43btxAly5dePVRKBQYM2YMFi9ejH79+qGiogJvv/02Ll68iJs3b6qdPFrqHK9YUgmPpddqrUxoiIEASI0K0Gi5WHFxMaKjo+Hv709VxQipB988odfqZN26dcPo0aOVd8UxxuDh4YHXXnsNn332mdb6ZGRkwN3dHSdOnMCwYcPUiq2lJt5Tt4oQuf6exv3OLOqGwT6Wjb5PIpHAysoKBgYGkMlkEIm0t+UPIW0N3zzBa6pBoVBgzZo18Pb2hoGBgbJ9/vz5ePTokVpj5Ofn4969exg0aJCyTSAQYNCgQYiNjdVaHzxbAgUAtra2asXWkqU8Lq/VZiAAbq/wR/FXvXB7hT9qFhhTZ7kYYwz379/HqVOncO/e08ROSZeQ5sEr8X722Wf48ccf8fHHH4PjOGV7cHAwVq5cqdYYOTk5AFBrHaijo6PyNW30qaysxIIFCxAcHIzevXvXG49MJkNRUZHKo6WRlivw2R/ZwLM90fAsqX43xQM+HU1hYWIAn46m2DzFAwbC515vZM80uVyOS5cuITExEV5eXujWrZtOjoeQ9orXqoYtW7Zg3759CAoKwtSpU5XtERERmDdvnkZjPZ+4q583Np+obh+O4zBt2jSkp6cjJiamwXGjoqKwYsUKjWLXtZVHMpFTJIe3owjH5nohQ1JV53Kw6eH2GNHdSq3lYnK5HKdOnaKqYoToEK/EKxaL4evrCzz7U7+aSCRCaWmpWmN06tQJAJCbm6vSnpubW2/xbE36VCfds2fP4ty5c+jcuXOD8SxduhSLFi1SPi8qKoKbm5tax6ILt7LK8fXpp8f99SQ3dHU0QVdHk3rfr+5yMUNDQ3h4eMDZ2blFzWUT0pbxmmrw8vJCfHw8UCPx7ty5Ez169FBrDFtbW/j7++PMmTPKNo7jcObMGYSFhSnb5HI5KisrNerDcRzefPNNREdH48yZM/D29m40HpFIBCsrK5VHS8HY0x2A5RzwQqA1RvZoWtnF6qVi1fPxvr6+Lep4CWnzGA8//fQTc3FxYd9//z0TCoXsyJEj7J133mFGRkZs9+7dao+zbds2ZmJiwvbt28ceP37M5s+fzywtLVlGRobyPdOnT2fdu3dXuw/HcezNN99kjo6OLDk5mZWXlysfcrlc7dikUikDwKRSqdp9msv+KwUMbyUw0dtX2IPciiaNVV5ezs6cOcP27t3L7t69q7UYCWmP+OYJXomXMcZ+/PFH5unpyQAwAMzd3Z1t3bpV43E2btzIvL29mZWVFQsLC2NxcXEqr8+ePZv97W9/U7tPXl4eE4lEdT40ia+lJN5SmYJ5LE1heCuBfXRQ3KSx8vLy2G+//cYOHTrEcnNztRYjIe0V3zzBax1vUlISevXqBTxb98lxHDp06NAcJ+R601LW8S7/LRMrjmTBzdYIt1f2gJkxv8I0jDGcO3cOHMchNDQUpqamWo+VkPZGp0Vyevfurbxvvy2sjW2pUvNkyuVjX0x045V05XI5ysvLYWlpidDQUBgZGVFVMUL0jNf/QHd3d7VvlCD8iCWV+OePj1BRxRDhY4kJf7PReIzi4mKcOnUKFy9eBGMMIpGIki4hLQCv/4UfffQRZsyYgdu3b9daU0uarroITsyDp0vzBnez0LhWglgsRnR0NBhjCA0NpVoLhLQgvOZ4TUxMIJPJAAAGBga19tqqqKjQXoR6oq853rqK4BgIgdTV6he5uX37Nq5duwYXFxf07dsXRkZN20GYEFI3nc7x7tu3j083oob4R6VN3hOtQ4cO6NmzJ1UVI6SF4pV4x44dq/1ICK49Lse83Rm12tUpcpOfn48HDx6gb9++cHBwoL3QCGnBaD/uFuLYdSkmff8QxRUcHC0NkVciB8caL3JTXVUsOTkZdnZ2qKqqgrGx+nV3CSG6R4m3BdhwOhcL9mSAY8CgbhY4MLsryiq5RovcyOVyJCQkICMjA97e3ujZsyetWiCkFaDEqydiSSVuZ1VgR3w+tsUVAADe6N8Bmya7w9hQCDtzNDqnm5GRgczMTISEhLSogj6EkIZR4tWDLTF5eGtHmspFtDUvueC9EU5qXQwrKiqClZUVOnfuDCcnJ5iZmTVvwIQQraK/S3VMLKmslXSFAmByPzu16hCnpKTgjz/+QH5+PgQCASVdQlohOuPVsXu5slrLxTjW+HKxiooKxMXFIS8vD4GBgbCzs2v+YAkhzYISr455O4ogFKDWDRINLReTSqU4f/48AGDQoEG0VIyQVo6mGnTM1dZY4z3RzM3N0bFjR0RGRlLSJaQN0Ov27i1Zc98yLJZUNrhcTC6XIykpCT4+PrC0bHxbdkKI7un0lmHSdA3tiVZcXIyLFy+itLQULi4ulHgJaWMo8bYwYrEYly9fhqmpKSIjI2kvNELaIEq8LYhMJsPly5fRsWNH9OnTh6qKEdJGUeJtASoqKmBoaAiRSITIyEhYWGhef5cQ0nrQqgY9y8vLw8mTJ5GSkgIAsLS0pKRLSBtHZ7x68nxVsQ4dOsDPz0/fIRFCdIQSrx4wxhAfH09VxQhppyjx6oFAIICVlRVVFSOknaLEq0NisRgymQxdu3aFv7+/vsMhhOgJ/X2rAxzHITk5GbGxscjLywPdLEhI+0ZnvM2sZlUxb29vWrVASDun98TLcRzi4+ORk5ODHj16wMvLSyt9+IzbHJKTk1FcXExVxQghf2F6VFhYyEJCQpizszOLjIxkZmZm7P33329yHz7j1iSVShkAJpVKNT4ujuNYWVkZY4yxiooK5b8JIW0L3zyh18Q7Z84c1q1bN1ZYWMgYY+zChQsMAIuOjm5SHz7j1sT3C1pVVcViY2PZ4cOHWWVlpUZ9CSGtC988obeLa4wx/PLLL5g+fTqsra0BAOHh4QgODsaOHTt49+EzrrYUFxfj1KlTyMzMRK9evajWAiGkTnqb4xWLxSgsLESPHj1U2gMCApCUlMS7D59x8axAjUwmUz4vKirS6HiysrIQFxdHVcUIIY3S2xmvVCoFgFp7h3Xo0AGFhYW8+/AZFwCioqJgbW2tfGh6Y4NIJIKzszOGDh1KSZcQ0iC9JV6R6OkeYyUlJSrtJSUlMDEx4d2Hz7gAsHTpUkilUuUjIyNDo+Oxs7NDv379aHqBENIovU01uLu7w9DQEOnp6SrtaWlp8PT05N2Hz7h4lrCrkzYhhDQnvZ7xDh06FHv37lW25eXl4fTp0xgzZoyy7fLlyzh+/LjafdQdlxBC9EWvm10mJSUhPDwcL730EkJDQ/HDDz9AIBDg4sWLyrPPGTNmIC4uDtevX1e7jzrvaUxzb3ZJCGn9+OYJvdZq6NWrF65evYpOnTohPj4e//jHP3D+/HmV5BgcHIxRo0Zp1Eed9xBCiL7Q9u71oDNeQkhjWuUZLyGEtEeUeAkhRMco8RJCiI5R4iWEEB3Tez3elqr6mqOmNRsIIe1HdX7QdI0CJd56FBcXAwBtRkkIaVRxcbGyGqI6aDlZPTiOQ2ZmJiwtLdXaqqeoqAhubm7IyMholcvPWnv8aAPH0NrjRxs4Bk3jZ4yhuLgYzs7OEArVn7mlM956CIVCuLq6atzPysqqVf7AVWvt8aMNHENrjx9t4Bg0iV+TM91qdHGNEEJ0jBIvIYToGCVeLRGJRFi2bFmrrQfR2uNHGziG1h4/2sAx6Cp+urhGCCE6Rme8hBCiY5R4CSFExyjxEkKIjlHirUdxcTESEhKQlpam1T58xuXrzp07uHr1qsq29U3tI5PJkJKSgvT0dI1vk9RUYWEhEhISIBaLm6VPbGwsrl692sQo68cYw82bN5GUlAS5XK7VPiUlJbh69Sry8/O1GHFt+fn5uHz5MrKzs7XaJzMzEwkJCbX2RtQ2juNw7do1pKSkQKFQqN3v2rVriI+P1/q4SozUsmnTJmZmZsZ8fX2ZmZkZGzt2LCstLW1yHz7j8pGZmcl69+7NbG1tmaenJ+vQoQP7/fffm9RHKpWyd955h9nY2LCePXsyR0dH1rNnT5aUlKT1+BljbN26dczExIT5+fkxU1NTNmnSJCaTybTW56uvvmJCoZD5+Pg0S/z37t1j/v7+zMHBgbm7u7NOnTqxmJiYJvfhOI59+OGHzMzMjAUGBjIPDw+2cOHCZjmGDz74gIlEIubv789EIhF76623mEKhaFKfvLw8FhERwaysrJQ/b3369GFpaWlajz8lJYV17dqVderUibm4uDAPDw929erVBvv8+OOPrFevXszW1pZ16NBBa+PWRIm3hsTERCYQCNiePXsYY4zl5OQwd3f3Bn+41enDZ1y+Ro4cyfr378/Ky8sZY4wtW7aMWVlZsSdPnvDuc/fuXbZhwwbl6zKZjE2cOJF16dJF6/GfPXuWCQQC9scffzDGGEtLS2OOjo5sxYoVWumTmJjI3Nzc2LRp05ot8fbt25eNHj2ayeVyxhhjc+bMYR07dmzwF606fZYvX85sbW1ZYmIiY88S8caNG7Ue/759+5ixsTGLjY1ljDF28+ZNZmlp2eBnqdNn3rx5zMPDg0kkEsYYY6WlpSwwMJBNmjRJq/ErFArm6+vLJk2axDiOY4wxNmXKFObp6cmqqqrq7bdkyRJ25coV9uWXX9aZePmOWxMl3hrmzZvHfH19VdpWrVrFbG1t6/1tr04fPuPy8fjxYyYQCNjBgweVbSUlJczU1JR9++23WuvDGGPHjx9nANjjx4+1Fj9jjL3++ussJCREpe3f//438/DwaHKf4uJi5uPjw3777Tf2/vvvN0viTUlJYQBUzlYzMzOZUChke/fu5d2nqKiImZubs7Vr12o95ppGjx7Nxo4dq9I2bdo01rt37yb1ee2119jIkSNV3jN16lQ2ZMgQrcXOGGPnz59nANj169eVbXfu3GEA2MmTJxvtX1/ibeq41WiOt4bExET07t1bpS04OBgSiaTeeVl1+vAZl4+kpCQwxlQ+y9zcHH5+fkhMTNRaHwC4fPkyLCws4OTkpLX40cDXKi0tDRKJpEl95syZgyFDhmDs2LFajblmLABU4unUqRNcXV3r/Xqq0ycuLg6lpaUYN24csrKykJiYCKlU2mzHUNfXMyUlBRzH8e7z73//G0lJSVizZg1OnTqFr7/+GsePH8fHH3+s9fhFIhG6d++ubOvWrRusrKwa/JnW1biUeGsoKChAhw4dVNqqnxcUFPDuw2dcvvE/P/bzn9VQ/Jr2SU5ORlRUFD744AMYGBhoKfq/4mmO78GOHTtw6dIlfPHFF1qNt65YzMzMYGJiUiuehuJvrE9mZiYAYPPmzejduzemTZuGjh074t13322WY6jr61lVVaUsmcqnT/fu3TF58mSsWbMGixcvxvLlyzFu3Dj06dOn2eNHI98DXY5LibcGIyMjVFRUqLSVl5cDAIyNjXn34TMu3/gB1PlZDcWvSZ979+5h1KhR+Pvf/44lS5ZoLfbn49H29yA/Px//+te/8Pbbb+PKlSuIiYmBWCxGeXk5YmJi6j2T5hu/TCarteqjse9BY32qv0/p6elIS0tDcnIyzp49iw0bNuCXX37RWvxoxv8H77zzDn7//Xc8fPgQV69eRXp6Om7duoUpU6Y0e/xo5Hugy3Ep8dbg4eGBx48fq7RVP6+vKLo6ffiMyzf+58d+/rPc3d2b3Of+/fuIiIhAREQEtm7dqlatYk3V97UyMjJCx44defWpqKhAYGAgdu/ejSVLlmDJkiU4d+4cnjx5giVLluDWrVtajV+hUCAnJ0fZxnEcsrOzG/weNNanc+fOAIDp06crk3C/fv0QFBSECxcuaC1+NPD1tLe3h6mpKe8+R44cwWuvvQY7OzsAgIWFBaZOnYqjR4/yW5bVQPwSiQRlZWXKNplMhvz8/Hq/B7oclxJvDcOGDcPZs2dV/pw6dOgQ+vbtCxsbGwCAVCpFTEyM8ouvTh913qMNvXv3hp2dHQ4fPqxsu3btGh49eoRhw4Yp25KTk3H37l2N+jx8+BAREREYOHAgfv75Z61PMVQbNmwYTpw4obKW+NChQxg8eLAy4eTn5yMmJgaVlZVq9XFxcUFMTIzKY/LkyXB3d0dMTAz69++vtfgHDBgAkUik8vU8d+4cCgsLVb6eCQkJePTokdp9+vTpAzs7O5XkplAokJ2dDQcHB63Fj2dfz5rJ8NChQyrxZ2dnIyYmRnmWrk4fBweHWmusMzIyYGdnp9WfpyFDhkAoFOLIkSPKtt9//x1yuRxDhw5VtsXFxWm0lljdcRul9mW4dqK0tJT5+PiwiIgIdvDgQfbRRx8xAwMDduLECeV7Tp48yQCwa9euqd1Hnfdoy8aNG5mJiQn76quv2J49e5ifnx8bPny4ynt69+7NJk+erHafzMxM5u7uzgIDA9nZs2fZhQsXlI/i4mKtxl9QUMDc3d3Z6NGj2eHDh9miRYtUlikxxtjevXsZAJaRkaF2n5qaa1UDY4ytWLGCWVlZsc2bN7Nff/2Vubu7s3/84x8q7/Hw8GDz58/XqM93333HHBwc2ObNm9mxY8fYq6++ymxtbbW+DlYsFjN7e3v2yiuvsMOHD7MZM2YwCwsLduPGDeV7vv32WwZAucRQnT7ff/89MzQ0ZCtXrmQnTpxgn3/+OTMzM2OffPKJVuNnjLGFCxcye3t7tnXrVvbzzz8zJycn9vbbb6u8x9rami1btkz5/MaNG+zChQts3rx5zNraWvkz/vySPnXGbQwl3jrk5uayRYsWsSFDhrBJkyaxs2fPqrx++fJlFhYWxh4+fKh2H3Xfoy179uxhL774Ihs2bBj75JNPaq0fnT59Olu5cqXafaqPua7H7du3tR6/WCxmc+fOZUOGDGGTJ09mcXFxKq+fOXOGhYWFsdzcXLX71PTtt9/WSmzatHXrVjZ27Fg2fPhw9vnnn9e6mWPChAnsq6++0qgPY4wdPHiQ/f3vf2cjRoxgCxcuZOnp6c0S/4MHD9isWbPYkCFD2LRp01hKSkqtOMLCwlRibKwPe7YMcerUqSwyMpJNnjyZ7du3r1niVygU7Ntvv2WjRo1iI0eOZBs2bFCuka42cuRItmXLFuXz9957r86f8ef/r6szbmOoLCQhhOgYzfESQoiOUeIlhBAdo8RLCCE6RomXEEJ0jBIvIYToGCVeQgjRMUq8hBCiY5R4SZPExcU1uEVKSxtX3/bv31+rngFpfyjxkib55ptv8O2337aacTUVGxuLy5cva228qVOnanW81mbv3r3K8pbtmaG+AyCtW2hoKITCtvv7+6uvvoKFhQX69u2r71DahH/+85/Yt28fnJ2d9R2KXlHibeViYmJgZmYGb29vJCYmQiaTISQkBJaWlirvk8vliI2NRV5eHvz8/ODr61vnOF5eXkhMTERZWRkGDRoEMzMzFBYW4uLFizAzM0P//v1V6o727t1bpTSkOvEkJSXh9u3bAABbW1sEBgbWW+6xMcnJyUhNTYWfnx+6devG65jri7W6XqyJiQl27doFABgxYgRu3LgBMzMzeHp6Ii4uDgYGBsoKXBKJBLGxseA4DqGhoXUWzVZHeno6EhMT4ebmhsDAQBw5cgRBQUEqpQcrKipw8eJFlJeXIyAgoFZZwl27diEiIgJyuRwpKSmws7NDcHBwrVKejY1TU/W4paWlSE5Oho+PD/z9/XH06FEUFxdDKBTCxcUFQUFBMDMzU/Y7fPgwOI7DhQsXUFJSAjMzM7zwwgu8Ymj1ml6KgujTiy++yIKDg5m3tzcbNWoU8/f3Z66urirVqlJTU5mPjw/r0qULGzFiBLO0tGRvvvlmrXGCgoKYm5sbGz16NOvSpQtzd3dnW7duZR4eHmzMmDHM3d2dhYSEqGzqN3nyZDZ16lSN4tmxYwebNGkSmzRpEhsyZAgzNzdn//3vf1XiqTluTU+ePGEDBw5k9vb2bNSoUczHx0el0pe6x9xQrD/++CNzd3dnnp6eyngfPHig7Ofp6clGjRrFPvjgA8YYY/v372cWFhYsNDSUhYWFMXNzc7Zz506VzzQ3N2f/+9//Gvye/vjjj8zY2JgNHDiQ9evXj4WFhTErKyuVsS5cuMA6duzIgoOD2ZgxY5itrS1bunSpyjgA2OjRo5mXlxcbO3Yss7e3Zy+88ILKe9QZpyYAbNSoUaxLly7s5ZdfZgcOHGCMMTZ//nw2adIkNnHiRNazZ0/m5uam3JSTMcbefvttJhQK2YABA9ikSZOUFb34xNDaUeJt5V588UVmbW3NUlNTGWOMyeVy1qdPH7Zo0SLle8aPH88GDRrEKioqGHtW+k4kErHdu3erjGNjY6OsdFVSUsI6dOjA7O3tWWZmJmPPSi9aWFiw/fv3K/vVlXgbi6emc+fOMRMTE5adnV3vuDW99NJLLCgoiBUUFCjbnt+sU91jbizWSZMmsenTp9f6mpuamrK7d+8q2woKCpidnR1bs2aNsu3zzz9nVlZWKhXUGku8+fn5zNLSkm3YsEHZ9uGHHzIAysRbVFTE7O3t2fbt25XvSU1NZdbW1iobLgJgw4cPV1YPu3//PjMwMGDnz5/XaJyaALCwsDBlOcj6LFy4kA0ePFilTSQSsd9++035nG8MrV3bnZxrR0aPHq3cRcLAwADh4eG4c+cO8OxPuEOHDuHdd9+FSCQCAPj7++Oll15S/vn8/DjVu2GYm5sjICAAL7zwAjp16gQ8mxbw9fVVFlDnE081iUSCs2fPYs+ePcjMzIRQKERKSopaxyuVSnHo0CF8+OGHsLW1Vba/+OKLvI65sVjrMnLkSHh7eyufnzhxAuXl5ViwYIGybd68eWCM4dixY2odFwD88ccfAIDZs2cr2/7973+rvOe3335DaWkpRCIR9u7di7179yI+Ph6dO3fGmTNnVN47ffp05dRQ165d4ebmpjw+Tcap6a233qq1PxyeFcs/duwYdu/eDQsLC1y6dKnBcZoSQ2tGc7xtQPU2KtVEIpFyX6j09HQwxuDp6anynq5du9ZKCM8nsepx6mqra88pdeMBgJ9++gnz5s2Dn58fXFxcYGxsDI7jkJubq9bxisVicBxXa063mibH3Fis9an+ZVQtLS0NLi4uykSPZ/tzubu717uL9NWrV1V+ib344ovIyMiAi4sLDA3/+q9pY2OjsktJamoqDA0NsX//fpXxfH19lb9E1Dk+TcZp7Pg5jsPUqVNx6NAhBAcHw87OTrlFTmlpKczNzescpykxtGaUeNu46i1hau6AWlBQAHt7e53Hw3Ec5s6di//+97/45z//qWwzNTWttdFjfaqTUH5+fp2v6+KYa16gsre3r3OX2YY+Mzk5WXmGCwDDhw+HnZ0dCgsLVd6nUChUtoyysrICx3H49ddfm7SipCnj1Dz+6OhoHDhwAA8ePFBeKD1y5Aiio6Mb/L5q61ham/ZzpO2Ura0t/Pz8cODAAWWbTCbDb7/9hvDwcJ3HU1JSgpKSEvj4+CjbDh8+rNw7TR0uLi7w9/fHzz//rNL+5MkTQMvHbGFhodYZcP/+/VFYWIhz584p2y5evIjs7Ox693N74403sGvXLuXDzs4OoaGhyMnJUfkT/ciRIyr7mA0fPhzl5eXYsWOHyniVlZXIy8tT+9i0NQ6e7b9mY2MDJycnZdu+fftqva/m11ObMbQmdMbbDqxfvx5jx45FVVUVAgICsH37dohEIixcuFDnsVhZWSEiIgKzZs3CnDlzIBaLsXHjRpU/0dWxadMmjBo1CsXFxRg2bBgePnyIhIQEREdHA1o85j59+mDZsmXYvHkzrKysMGLEiDrf5+vrizlz5mDChAl47733IBQKsXbtWsycOROBgYFqf1737t3x+uuvY/z48Vi8eDEqKyuVX5/qs8xu3brhk08+wcyZM5GQkIDAwECkpqZi3759+OGHH9Q+q9fWOHi2CWRxcTFef/11RERE4PTp0zh69Git9/Xp0wcbN25ERUUFrKys8MILL2gthtaEznhbuQEDBqB3794qbYGBgYiIiFA+Hz58OOLj42FsbIw///wT48aNw+XLl2FhYdHgOIMHD0ZQUJBKW2RkJHr06KF8Hhoain79+mkUz6FDh/DKK6/gwoULKCsrw9mzZzFt2jTl9uV1jVvXcV+7dg1+fn6IjY2Fvb09Dh482ORjrhnrjBkzsGLFCly6dAkHDx6ERCKpsx+e3WyxceNG3L17F7du3cL69etr3X03YcIEuLq61ntcAPDDDz/g448/RnJyMkpKSnDy5EkAUFkL/cEHH+Ds2bMwNjbGhQsXYGpqiuPHjyMsLEz5nkmTJtVaHz1mzBiVuXF1xqmprnFdXV0RHx8PBwcHnDt3DoGBgTh58iQmTZqk3Bkaz+b3BwwYgD/++EM5zcInhtaO9lwjpIUpKChQuSgWHR2NESNGICMjo93f8dVW0FQDIS3Mt99+i7t372LAgAEQi8X46quvMGfOHEq6bQid8RLSwnAch127din/7I6IiMC4ceP0HRbRIkq8hBCiY3RxjRBCdIwSLyGE6BglXkII0TFKvIQQomOUeAkhRMco8RJCiI5R4iWEEB2jxEsIITpGiZcQQnTs/wFktF5RFBHZvQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 360x340 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nominal = np.linspace(0.005, 0.1, 30)\n",
    "realized = [np.mean(score[test] >= np.quantile(score[fit], 1 - a)) for a in nominal]\n",
    "plt.figure(figsize=(3.6, 3.4))\n",
    "plt.plot([0, 0.1], [0, 0.1], \"--\", color=\"#aaa\", lw=1)\n",
    "plt.plot(nominal, realized, \"o-\", ms=3, color=\"#06c\")\n",
    "plt.xlabel(\"nominal control-gene rate\")\n",
    "plt.ylabel(\"realized (held-out controls)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "calib-read",
   "metadata": {},
   "source": [
    "The held-out control genes come through above the rate the cutoff promises at 1% (2.1%) and close to it at 5% (5.7%). The control set is small, a few hundred genes, so the held-out rate is coarse and does not track the line exactly. The split shows the control set is internally consistent; whether unexpressed genes stand in for a silent expressed gene is a separate assumption the biological check below speaks to."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "corum-intro",
   "metadata": {},
   "source": [
    "## The hits in biological terms\n",
    "\n",
    "A hit list that is right should be enriched for genes that work together. CORUM records which genes share a protein complex, and {func}`~mantispy.ds.corum` returns it. Members of one complex tend to give the same phenotype when knocked out, so they should be over-represented among the hits and rare among the control genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "corum",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:10.778146Z",
     "iopub.status.busy": "2026-10-21T00:50:10.778011Z",
     "iopub.status.idle": "2026-10-21T00:50:10.949940Z",
     "shell.execute_reply": "2026-10-21T00:50:10.948449Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p<=1% hits (331) in a complex: 59%   odds ratio 4.7\n",
      "q<=1% hits (102) in a complex: 76%   odds ratio 10.4\n",
      "other tested genes in a complex: 23%\n"
     ]
    }
   ],
   "source": [
    "complexes = mt.ds.corum()\n",
    "in_complex = set(complexes[\"target\"].astype(str))\n",
    "membership = np.isin(genes, list(in_complex))\n",
    "hit = target & (np.nan_to_num(pv, nan=1.0) <= 0.01)\n",
    "q_hit = target & (np.nan_to_num(qv, nan=1.0) <= 0.01)\n",
    "\n",
    "\n",
    "def enrichment(mask):\n",
    "    \"\"\"The in-complex fraction and Fisher odds ratio for a hit mask, against the other tested genes.\"\"\"\n",
    "    table = [\n",
    "        [int((mask & membership).sum()), int((mask & ~membership).sum())],\n",
    "        [int((target & ~mask & membership).sum()), int((target & ~mask & ~membership).sum())],\n",
    "    ]\n",
    "    return np.mean(membership[mask]), fisher_exact(table)[0]\n",
    "\n",
    "\n",
    "for name, mask in [(\"p<=1%\", hit), (\"q<=1%\", q_hit)]:\n",
    "    frac, odds = enrichment(mask)\n",
    "    print(f\"{name} hits ({mask.sum():>3}) in a complex: {frac:.0%}   odds ratio {odds:.1f}\")\n",
    "print(f\"other tested genes in a complex: {np.mean(membership[target & ~hit]):.0%}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "corum-fig",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-21T00:50:10.953501Z",
     "iopub.status.busy": "2026-10-21T00:50:10.953248Z",
     "iopub.status.idle": "2026-10-21T00:50:11.018825Z",
     "shell.execute_reply": "2026-10-21T00:50:11.017921Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAE2CAYAAAB87RlzAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAANkZJREFUeJzt3Xl4TOf/PvB7YomQZGInRFQ0qZDQWpNYa0lpiViLVi2h9hJRBtVQxFJLS4vai1JC41MtagmtLaSlBLGWhFRs2SNDzPv3x0/Pt2m2mTEymbhf15XrynnOeZ55H2nvnDxnU4mIgIiILJaVuQsgIqLnwyAnIrJwDHIiIgvHICcisnAMciIiC8cgJyKycAxyIiILxyAnIrJwxc1dgDnodDrExcXBzs4OKpXK3OUQEWUjIkhJSYGjoyOsrPI+5n4pgzwuLg5OTk7mLoOIKF+xsbGoXr16ntuYPchTUlIQFhaG+Ph4eHh4wNfXV69+UVFROHToEEqXLg0/Pz+UL19e78+0s7MDnv0D2dvbG107EdGLkpycDCcnJyWv8mLWIL916xaaN2+OSpUqoWHDhli4cCF8fHywdevWPKc8xo0bh9WrV6Nnz54oXbo0Fi1ahB9++AG1a9fW63P/Gdve3p5BTkSFmj7Tv2YN8okTJ6JixYo4evQoSpQogbFjx6Ju3brYvn07evTokWOf1atXY9myZYiIiED9+vUBAPfu3cPjx48LuHoiosLBbFetPH36FGFhYejfvz9KlCgBAHBzc0OLFi0QGhqaa7/Fixfj3XffVUIcACpWrIhq1aoVSN1ERIWN2Y7IY2JikJ6eDldX1yztrq6uiIiIyLFPamoqoqKiMH78eBw6dAi///47HB0d0alTJ6jV6lw/S6vVQqvVKsvJyckm3BMiIvMy2xF5amoqAGQLYAcHB2XdfyUmJgIAVq5ciUmTJuH27dtYunQpXn31VZw9ezbXzwoJCYFarVa+eMUKERUlZgvyMmXKADkcHSclJSnrcusDAMePH8fChQtx9OhR1KtXDx9//HGun6XRaJCUlKR8xcbGmmw/iIjMzWxBXqNGDdjY2ODq1atZ2q9cuQI3N7cc+5QtWxaVKlWCj49PljO5zZs3R3R0dK6fZW1trVyhwitViKioMVuQFy9eHJ07d8aGDRuQmZkJALh27Rp+/fVXdOvWTdnup59+wjfffKMs9+zZEydOnMgy1vHjx1GnTp0CrJ6IqPBQmfOdnTdu3ICPjw9q1aqFxo0bIzQ0FPXr18fOnTuVW1IDAgJw4sQJREVFAQDu37+PFi1aKEfmEREROHfuHMLDw1G3bl29Pjc5ORlqtRpJSUk8OieiQsmQnDJrkANAQkICtm3bptzZ2aVLlyzPFQgLC0NsbCxGjx6ttD169AihoaG4efMmatSoAT8/vzyvWvkvBjkRFXYWFeTmYGyQu+yf9ELrKqqutZtj7hKILI4hOcXH2BIRWTgGORGRhWOQExFZOAY5EZGFY5ATEVk4BjkRkYVjkBMRWTgGORGRhTMqyC9cuJDrul9++eV56iEiIgMZFeSNGjXCV199laVNq9Vi3LhxePvtt01VGxER6cGoIF+5ciUmT56Md955B3fv3sW5c+fQqFEjhIWF4fDhw6avkoiIcmVUkPfr1w9nzpxBQkIC6tWrh8aNG6N+/fr4888/4e3tbfoqiYgoV0af7CxdujTs7OyQlpaGzMxMvP7667CzszNtdURElC+jgvynn36Ch4cHHj16hOjoaGzbtg0hISFo3749bt++bfoqiYgoV0YFub+/P8aOHYvw8HA4OTnB398fZ8+ehZWVFTw9PU1fJRER5aq4MZ2OHj2Kxo0bZ2lzdHTE3r17sWjRIlPVRkREejDqiPy/If4PlUqFwMDA562JiIgMYFSQiwhWrVoFLy8vVK5cWWmfOnUq4uLiTFkfERHlw6ggX7p0KYKDg9GzZ0/cvXtXaa9ZsyZmzpxpyvqIiCgfRgf5tm3bsk2j+Pr6Yvv27aaqjYiI9GBUkN+4cQMNGjQAns2L/8POzg6JiYmmq46IiPJlVJDXqFEDZ86cAf4T5GFhYXBzczNddURElC+jLj8cM2YMPvjgA8ybNw8AEBERgT179mDevHn48ssvTV0jERHlwaggHz16NDIyMjBo0CDodDo0a9YMarUa06dPx+DBg01fJRER5cqoIAeACRMmIDAwENevX4dOp0OtWrVQokQJ01ZHRET5MjrIAaBYsWJ49dVXTVcNEREZTO8gHzBggN6Drlu3zth6iIjIQHoHeWZm5outhIiIjKJ3kG/cuPHFVkJEREYx+sUSePbMlbi4OMTFxUFEjBpj0aJFcHZ2RqlSpdC4cWMcOXIkz+2HDRsGlUqV5atevXpG7gERkeUzKsifPHmC4OBglC1bFtWqVUO1atVQrlw5zJgxw6ApmNWrV2PKlCn46quvcPv2bbz55pt46623EBMTk2e/7t27Q0SUr6ioKGN2g4ioSDAqyAMDA7Fy5UrMmTMHkZGRiIyMREhICJYvX47x48frPc78+fMREBCAd955B+XLl8ecOXPg4OCAZcuWGVMWEdFLyajLDzds2IDdu3fDy8tLaWvYsCE8PT3RqVMnfPHFF/mO8fDhQ1y6dAmzZ89W2lQqFVq3bo1jx47l2XfPnj0oXbo01Go1WrRogXnz5qFmzZrG7AoRkcUz6ojc2toaderUydZep04dlCpVSq8x7ty5AwCoWLFilvZKlSop63JSu3ZtbN68GXfu3MHBgweRlJSEFi1aICkpKdc+Wq0WycnJWb6IiIoKo4K8ffv2WLx4cbb2L774Au3bt3+ugkQky4O4/isoKAidO3eGvb096tSpg++//x7379/H999/n2ufkJAQqNVq5cvJyem5aiQiKkyMmloREUyfPh3btm1Do0aNICKIjIzExYsX0a9fP4waNUrZdunSpTmOUaVKFQDAvXv3srTfu3cvy1uH8uPg4IDq1avjypUruW6j0WiyPDs9OTmZYU5ERYZRQZ6WlgY/Pz8AUKY0XF1d4erqitTUVKSmpuY7Rrly5eDm5obw8HB069YNePYL4tChQ3j//ff1riUxMRG3bt2Co6NjrttYW1vD2tpa7zGJiCyJUUEeFhZmkg8PCgrCRx99hLfeegteXl6YO3cuEhMTMWzYMGWbgIAAnDhxAlFRUdBqtXj33Xeh0WhQt25dxMTEYOzYsbCzs0O/fv1MUhMRkaV5rodmPa+AgAAkJydj+PDhiI+Ph4eHB/bs2QNnZ+cct7e2tkZAQACCgoJw+vRplC1bFi1atEBERAQqVapU4PUTERUGKjHylszIyEgcPXoUCQkJ2dYFBweborYXJjk5GWq1GklJSbC3t9e7n8v+SS+0rqLqWrs55i6ByOIYklNGHZHPmzcPkyZNgru7OxwcHLKtL+xBTkRUlBgV5AsXLsSuXbvQqVMn01dEREQGMeo68rS0NLRs2dL01RARkcGMCvK33nrLZFeuEBHR8zFqauWLL75A/fr1ERoaChcXl2x3Yn7++eemqo+IiPJhVJB/9tlnSEhIwJ07d5CRkWH6qoiISG9GBfmmTZvwyy+/4M033zR9RUREZBCj5shLlSqFJk2amL4aIiIymFFB3rJlS2zZssX01RARkcGMmlopVqwYhg4dim3btqF27drZTnbm9sRDIiIyPaOCXKvVokuXLgCA27dvm7omIiIygFmffkhERM/PqDlyIiIqPIx+jO3FixexYMECXLx4ESICd3d3jB8/Psd3eRIR0Ytj1BH5zz//DE9PT0RHR8PLyws+Pj6Ijo6Gp6cndu/ebfoqiYgoV0YdkU+ZMgWzZ8/GhAkTsrTPnz8fkydPRseOHU1VHxER5cOoI/Lz589j6NCh2dqHDBmCCxcumKIuIiLSk1FBXqFCBZw7dy5b+9mzZ1GhQgVT1EVERHoyampl4MCB6NWrF6ZNm6bcqh8REYEZM2Zg8ODBpq6RiIjyYFSQz5gxA9bW1vj444+RkpICALCzs0NQUBCmTJli6hqJiCgPRt+iP23aNEydOhWxsbFQqVSoXr06rKx4WToRUUEz+jpyALCysoKzs7PpqiEiIoMZdQh94MABBAQEZGsPCAjAwYMHTVEXERHpyaggnzBhAkaNGpWtfeTIkZg4caIp6iIiIj0ZFeQXLlxArVq1srW7uLggKirKFHUREZGejArymjVrYt++fdna9+7dCycnJ1PURUREejLqZOeoUaMwZMgQ3Lx5Ey1btoSI4Ndff8Vnn32GmTNnmr5KIiLKldFBnpKSghkzZiApKQkAoFarMXHixBznzomI6MUx+vJDjUaDCRMm4Nq1a1CpVKhVqxaKF3+uqxmJiMgIz5W8xYsXh5ub23MVkJCQgNDQUMTHx8PDwwNdunTJ9g7Q3Pz555/Ytm0bvL290alTp+eqg4jIUpn1VsybN2/Cw8MD3377LR48eIDRo0fDz88POp0u374pKSno1asXlixZgl9++aVA6iUiKozMGuQff/wxnJycEB4ejkWLFiE8PBy7d+/G1q1b8+374YcfokePHnjllVcKpFYiosLKbEGemZmJH3/8Ee+//74yt+7i4oKWLVtix44defZds2YNrly5guDg4AKqloio8DLb2cmYmBg8evQItWvXztL+6quv4vjx47n2i46OxqRJk3DkyBGUKFFCr8/SarXQarXKcnJy8nNUTkRUuBgU5M2bN9druyNHjuS7TVpaGgDA3t4+S7tarVbW/VdGRgZ69+6NmTNnwtXVVa9aACAkJATTp0/Xe3siIktiUJAfPXoUzs7OaN++PYoVK/ZcH2xrawsAynXo/0hMTFTW/dfGjRtx8+ZNxMTEYOrUqQCAO3fu4Pjx45g6dSpmzJiR46N0NRoNAgMDleXk5GTegUpERYZBQT59+nSsXbsWu3fvxoABAzBw4EC4uLgY9cE1atRA6dKlcfnyZfj6+irtly9fxmuvvZZjnwYNGiAoKChLm0qlQrFixVCqVKlcL1u0traGtbW1UXUSERV2KhERQzqICPbt24dVq1bhf//7H7y9vTF48GB069YNNjY2Bn143759ceXKFRw7dgwlSpTApUuXULduXWzZsgU9evQAAISFhSE2NhajR4/OcYwGDRqgdevWWLx4sd6fm5ycDLVajaSkpGxTO3lx2T9J723p/1xrN8fcJRBZHENyyuCrVlQqFTp06ICtW7fi9u3b6Ny5Mz777DM4OjoaXOjcuXNx7949+Pj4YMSIEWjTpg26du2K7t27K9vs2rULK1asMHhsIqKXxXNdtXL16lVER0fj77//NmqKxcnJCefOncMPP/yA+Ph4rFmzBr6+vlmmSPz9/ZUXPOdk1KhRfEsREb3UDA7ye/fuYcOGDVizZg3i4uLQr18/HD58GA0aNDCqADs7O/Tv3z/X9W+//Xae/XN6UxER0cvEoCDv3r07fvrpJzRv3hxTpkxBt27deBKRiMjMDDrZqVKp4OTkhBo1auS5nT7XkZsTT3YWLJ7sJDKcITll0BH5yJEjn7c2IiIyMYOCfOnSpS+uEiIiMopBQX7//v0c221tbVGqVClT1URERAYwKMgrVqyYY7tKpYKzszOCgoIwYsQIvV8MQUREz8+gIP/hhx9ybE9JScHZs2eh0WhgY2ODQYMGmao+IiLKh0FB3rVr1zzXN2vWDLNmzWKQExEVIJO+WKJ9+/a4evWqKYckIqJ8mDTIU1NTUaZMGVMOSURE+TBZkOt0OsybNw9vvPGGqYYkIiI9GDRH/u677+bYnpKSgvPnz+P+/fuF/q5OIqKixiTXkdvZ2eHdd9/FiBEj8r19n4iITMugIN+/f/+Lq4SIiIxi0pOdAPDXX3+ZekgiIsqDyYL8wYMHGDduXK7v2yQiohfDoCBPSkpC//79Ub58edja2iIgIABarRbr169H7dq1ERoaimXLlr24aomIKBuD5sgnT56M/fv3Y/jw4QCA1atXIzY2Fr/99humTp2KwMBAPjyLiKiAGRTkP/74I8LCwpR3aL799tvw9vbGjh074O/v/6JqJCKiPBg0tRIXF4eGDRsqy40bNwb0eK8mERG9OAYF+dOnT1GsWDFluXjx/39AX7JkSdNXRkREejFoagUAGjVqlG9bZGTk81VFRER6e+53djZr1syU9RARkYH4zk4iIgtn8NTKP+7cuYMrV65ARODq6ooqVaqYtjIiItKLwXd2XrlyBW3btkXVqlXRsmVLtGrVClWrVkXbtm1x5cqVF1MlERHlyqAj8ps3b8Lb2xtVq1bFypUr4ebmBgC4dOkSvvzyS3h7eyMyMhLOzs4vql4iIvoPlYiIvhsPGDAA6enp2LJlC6yssh7M63Q69O7dG7a2tli7du2LqNVkkpOToVarkZSUBHt7e737ueyf9ELrKqqutZtj7hKILI4hOWXQEfmePXvw22+/ZQtxALCyssKsWbPQunVrwyt+9osgp3HzIiJQqVRGfR4RUVFhUHI+fPgQNWvWzHV9zZo1c335RG5CQkJQuXJllChRAh4eHjh48GCe20dHR6Nv376oUKECbGxs8MYbb+CHH34w6DOJiIoSg4K8cuXKuHDhQq7rL1y4gKpVq+o93vLlyzF79mxs2rQJSUlJ6NatG9555508n2m+cuVK9OrVC9euXUNCQgL69u2LHj164PTp04bsChFRkWFQkPv5+WHs2LFIS0vLti4tLQ3jxo1Dly5d9B5v4cKFGDx4MNq1awdbW1sEBwejQoUKWL58ea59FixYgK5du0KtVsPGxgaBgYGwsrLCH3/8YciuEBEVGQbNkU+bNg1NmjSBi4sLAgIClJdIREdHY+XKlbCxscH333+v11gPHjzAlStX0KpVK6VNpVKhVatWOH78eJ59dTodHj9+jOTkZKxYsQIODg7w9fU1ZFeIiIoMg4K8UqVKOHHiBIKCgrBw4UI8evQIAGBjY4Pu3btj/vz5qFSpkl5jxcfHAwAqVqyY7TNOnjyZZ98DBw6gc+fO0Gq1cHBwwObNm1G9evVct9dqtdBqtcpycnKyXjUSEVkCg+/srFKlCjZu3IjHjx8jNjYWAODk5GT0ExB1Ol225fyuRGnfvj0yMjKQmpqKFStWwM/PDwcPHoSPj0+O24eEhGD69OlG1UdEVNgZfGdnUlIS8OzRtS4uLnBxcVFC/J91+vjnpOjdu3eztN+9e1fv2/1tbW0xfvx41K9fH+vWrct1O41Gg6SkJOXrn19ARERFgUFBvnfvXvTp0yfX9e+++y5++eUXvcYqW7Ys3N3dER4errTpdDqEh4dnObLOzMzE48eP8xzr0aNHWZ6T/l/W1tawt7fP8kVEVFQYFOSLFy9GYGBgrusDAwOxaNEivcebOHEi1qxZg+3btyMuLg6BgYFITU1V3gkKAMOGDcMbb7wBPLsypnPnzjhy5AgePnyI69evIzAwENHR0RgwYIAhu0JEVGQYNEd+5swZeHp65rre09MTf/75p97j9e/fH6mpqdBoNIiPj4eHhwf27duX5cRliRIlYG1tDQAoU6YMxo4dixkzZuD06dMoU6YMXn/9dRw5cgRNmzY1ZFeIiIoMg561Ym1tjaSkJJQqVSrH9RkZGXBwcEBGRoYpazQ5PmulYPFZK0SGMySnDJpaqVGjRp6vcTt16hSffEhEVMAMCvKuXbti0qRJSE1NzbYuJSUFkyZNQteuXU1ZHxER5cOgOXKNRoPGjRujdu3aGDJkCNzc3CAiuHz5MlauXAk7OztMmsTpByKigmRQkJcrVw7Hjh1DUFAQFixYkOXOzh49emDBggUoW7bsi6qViIhyYPCdnZUrV8aGDRuwevVqxMbGQqVSoXr16kbf2UlERM/H6Jcv/3NnJxERmZfBt+gTEVHhwiAnIrJwDHIiIgvHICcisnBGn+zUarX4/fffERMTg8zMzCzr3nvvPVPURkREejAqyM+fP4/OnTsjPj4e6enpyvMA8OzyRAY5EVHBMWpqZdy4cfD390dKSgoAIDExEZcuXULTpk0xbtw4U9dIRER5MCrIT506BY1GAysrK6hUKjx58gSurq5Yu3Ytli1bZvoqiYgoV0YFeWJiIipUqAAAqFChAuLi4oBn7+68c+eOaSskIqI8PfdVK97e3pg+fTrOnDmDKVOm4LXXXjNNZUREpBejTnaOHz9e+X7u3Lno2rUr1q5diypVqmDbtm2mrI+IiPJhVJB//vnnyvdubm64ePGi8hYLlUplyvqIiCgfRl9H/l9qtdpUQxERkQF4ZycRkYVjkBMRWTgGORGRhTPZHDlRQfjmd3NXYJmGNjR3BfQiGRXkWq0Wy5cvx9GjR/Hw4cNs6/fv32+K2oiISA9GBfmoUaOwc+dO+Pv7o1atWqavioiI9GZUkO/YsQMHDx5E/fr1TV8REREZxKiTnSVLluSROBFRIWFUkHfs2BGbN282fTVERGQwo6ZWHj9+jGHDhmH79u2oXbt2ttvyly5daqr6iIgoH0YFeWpqKrp06QIAuH379nMVcO/ePWzZsgXx8fHw8PBAjx49UKxYsTz7HD58GMeOHUPx4sXRvHlzeHl5PVcNRESWzKggDwsLM8mHX7t2DT4+Pqhbty6aNGkCjUaDNWvW4Oeff84xzHU6HZo1a4YyZcrAy8sL6enp8PX1xYABA/Dll1+apCYiIktj1huCJk6ciNq1a2Pfvn2wsrLChx9+CFdXV2zZsgX9+vXLtr1KpcLXX3+NRo0aKW2+vr7o1KkThg4dinr16hXwHhARmZ/eQT5gwAAAwLp165Tvc7Nu3bp8x8vMzMRPP/2EhQsXwsrq/59zrVmzJlq1aoWwsLBcg/zfIQ4Ar7/+OgDg1q1bDHIieinpHeSZmZk5fm+smzdvIiMjAy4uLlnaXVxccPz4cb3H+fbbb1GqVKlsAf9vWq0WWq1WWU5OTjayaiKiwkfvIN+4cWOO3xsrPT0dAGBvb5+lXa1WIy0tTa8xfv31V3zyySeYP3++8g7RnISEhGD69OnPWTERUeFktqcf2traAs9e5PxviYmJsLOzy7d/REQEOnfujPHjx2PMmDF5bqvRaJCUlKR8xcbGPmf1RESFh9mCvEaNGihTpgwuXbqUpT06Ohp16tTJs+/Jkyfh6+uLYcOGYfbs2fl+lrW1Nezt7bN8EREVFWYL8mLFiqFbt25Yv369Mn8dFRWFI0eOoGfPnsp227Ztw4IFC5TlU6dOoUOHDhg2bBjmzp1rltqJiAoTs15+OHfuXLRo0QJNmzbFG2+8gV27dqFPnz7w9/dXttm7dy9OnDiB8ePHIy0tDb6+vrC2tkZmZiaCgoKU7Xr16oUmTZqYaU+IiMzHrEFetWpVnD17Frt27UJ8fDz69++P1q1bZ9mmV69eaN68OQDAysoKkydPznEsGxubAqmZiKiwMfuLJUqXLo1evXrlur5Dhw7K9zY2NlmOwomIiC+WICKyeHyxBBGRheOLJYiILBxfLEFEZOH4YgkiIgtn9hdLEBHR8zHriyWIiOj5me0WfSIiMg2j7+y8ePEiFixYgIsXL0JE4O7ujvHjx+f7wCsiIjIto47If/75Z3h6eiI6OhpeXl7w8fFBdHQ0PD09sXv3btNXSUREuTLqiHzKlCmYPXs2JkyYkKV9/vz5mDx5Mjp27Giq+oiIKB9GHZGfP38eQ4cOzdY+ZMgQXLhwwRR1ERGRnowK8goVKuDcuXPZ2s+ePZvnK9eIiMj0jJpaGThwIHr16oVp06YpzwCPiIjAjBkzMHjwYFPXSEREeTAqyGfMmAFra2t8/PHHSElJAQDY2dkhKCgIU6ZMMXWNRESUB6OCvFixYpg2bRqmTp2K2NhYqFQqVK9eHVZWvCydiKigPdcbgqysrODs7Gy6aoiIyGB6B/mAAQMAAOvWrVO+z826deuevzIiItKL3kGemZmZ4/dERGReegf5xo0ble/Hjh2LRo0a5bhdZGSkaSojIiK9GHV2snHjxkatIyIi0zPpZSbJycmwtbU15ZBERJQPg65aCQgIyPF7ANDpdDh79iyPyImICphBQZ6amprj9wBQokQJdOzYEcOHDzdddURElC+DgnzLli3As2et8L2cRESFg1Fz5AxxIqLCw6ggP3DgQLY5cjybNz948KAp6iIiIj0ZFeQTJkzAqFGjsrWPHDkSEydONEVdRESkJ6OC/MKFC6hVq1a2dhcXF0RFRRlVSEZGhkHbP336FImJibzLlIheekYFec2aNbFv375s7Xv37oWTk5NBY3366adwcHCAra0tXn31VezZsyfP7ePi4hAcHAxnZ2eULVsWR44cMbh+IqKixKggHzVqFIYMGYKFCxciMjISp06dwoIFCzBkyBCMGTNG73GWLFmCxYsXY9euXUhPT8fAgQPRtWtXXL16Ndc+69evh4ggNDTUmNKJiIocox5jO2rUKKSkpGDGjBlISkoCAKjVakycODHHufPcLF68GAEBAWjevDkAYPLkyfjmm2+wfPlyfP755zn20Wg0AIBbt24ZUzoRUZFj9PPINRoNJkyYgGvXrkGlUqFWrVooXlz/4e7fv4/r16+jRYsWWdpbtmyJiIgIY8siInrpPNeLJYoXLw43Nzej+t69exd4dnPRv1WqVMnkQa7VaqHVapXl5ORkk45PRGRORge5VqvF77//jpiYmGxXjrz33nt6j6PT6bIsZ2ZmQqVSGVtWjkJCQjB9+nSTjklEVFgYFeTnz59H586dER8fj/T0dKjVamWuvHLlynoFuaOjIwAgPj4+S/vdu3eVdaai0WgQGBioLCcnJxt8dQ0RUWFl1FUr48aNg7+/P1JSUgAAiYmJuHTpEpo2bYpx48bpNYaDgwPq1auHAwcOKG06nQ4HDx5UTn4CwKNHj5TPMZa1tTXs7e2zfBERFRVGBfmpU6eg0WhgZWUFlUqFJ0+ewNXVFWvXrsWyZcv0Hmfy5MlYu3YtNm3ahOvXr2PkyJHQarVZnqA4evRoeHl5KcuPHz9GYmKiMs+dmpqKxMREg28oIiIqKowK8sTEROUkZYUKFRAXFwcAcHJywp07d/Qep0+fPli2bBnmzp0LLy8vXLlyBQcPHkTVqlWVbUqXLp3lCDo0NBQ1a9aEt7c31Go13nvvPdSsWROLFy82ZleIiCyeSkTE4E4qFf7p1rVrV5QrVw5jxozB2rVrcfjwYZw5c+ZF1GoyycnJyry+IdMsLvsnvdC6iqpr7eaYbKxvfjfZUC+VoQ3NXQEZypCcMupk5/jx45Xv586di65du2Lt2rWoUqUKtm3bZsyQRERkJKOC/N93Xbq5ueHixYvKbw1TXzpIRER5M2qOPKf3cqrVaoY4EZEZGBXkly9ffu5LAomIyDSMCnI/Pz+sXbvW9NUQEZHBjJoj1+l0+Oijj7B161a4u7ujZMmSWdbznZ5ERAXHqCBPTU2Fn58f8K+HXxERkXnoHeSHDh1C69atAQBhYWEvsiYiKsTmzp1r7hIszot+l7Hec+Rt2rRRvq9SpcqLqoeIiAykd5CXLVsWN27cAHJ4YiEREZmP3lMrvXv3Rt26dVG9enUAwGuvvZbrttHR0aapjoiI8qV3kH/99dfo3r07rl69iuHDhxv0bk4iInpx9A5ylUqFdu3aoV27djhx4gSDnIiokDDqhqB169aZvhIiIjKKUUFORESFB4OciMjCMciJiCwcg5yIyMIxyImILByDnIjIwjHIiYgsHIOciMjCMciJiCwcg5yIyMIxyImILByDnIjIwjHIiYgsHIOciMjCMciJiCyc3i+WeFH+/vtvbNy4EfHx8fDw8EDfvn1RokQJk/chIiqqzHpEfvnyZXh4eCA8PBx2dnaYNWsWOnTogKdPn5q0DxFRUWbWI/KJEyeiXr16+Omnn6BSqTB48GC4uLhg06ZN6N+/v8n6EBEVZWY7In/y5Al2796Nvn37QqVSAQCqV6+O1q1bY+fOnSbrQ0RU1JntiDwmJgZarRavvPJKlvZatWrh6NGjJusDAFqtFlqtVllOSkoCACQnJxtUsy5Nq8dW9F+G/jvn5VGqyYZ6qZjwR4CMjAzTDfaSMOb/gX/6iEi+25otyB89egQAsLOzy9Jub2+P9PR0k/UBgJCQEEyfPj1bu5OTk1G1k2HUWGzuEl56Y81dwEsuODjY6L4pKSlQq9V5bmO2ILe1tQUAJCYmZmlPSEiAvb29yfoAgEajQWBgoLKs0+nw8OFDlC9fXpmisWTJyclwcnJCbGxsnv8O9OLwZ2B+Re1nICJISUmBo6NjvtuaLchr1KgBW1tbXLx4EW+99ZbSfvHiRbi7u5usDwBYW1vD2to6S5uDg4NJ9qMwsbe3LxL/AVsy/gzMryj9DPI7Ev+H2U52WllZoWfPnli3bp0yZXLmzBkcO3YMvXv3Vrb77rvvMHv2bIP6EBG9TFSiz0z6C3L37l20bt0aANCgQQPs2bMH/v7+WL16tbJNQEAATpw4gaioKL37vGySk5OhVquRlJRUZI5ELA1/Bub3Mv8MzHodeaVKlXD69Gns3bsX8fHxGDNmDJo1a5Zlm379+qFt27YG9XnZWFtb49NPP802fUQFhz8D83uZfwZmPSInIqLnx4dmERFZOAY5EZGFY5ATEVk4BjkRkYVjkBMRWTgGeSH2+PFjc5dARBaAQV6IpKam4vz58wCAdevWYfjw4VnWa7VafPfdd+jTpw8++OADM1X5cktNTcXq1avh7++PCRMmmLucl9KDBw/w5ZdfomPHjli8mA9kA4O84N26dQu//fZbjuuioqKwYsUKAICPjw/27t2LnTt3YvTo0Xj8+DHat2+P9evXo0OHDpg1a1YBV05arRaNGzfG7t270atXLwa5GcTHx6NBgwY4c+YMhg4dioEDB5q7pMJBqEAdOnRIGjdunOv648ePS0JCggQHBwsA8fb2lvnz58udO3fE3t5e0tLSJDMzU3Q6XYHW/bIJDQ2V5s2bi7OzswwbNkwyMjIkOjpaXFxcJDMzkz+DArBixQpp1KiR1KpVS6ZOnSpPnz6VPXv2SOvWrUVE5MmTJ+YusdBgkBcwrVYr9vb28vDhQxERiYuLk+TkZGV9ly5d5LvvvpO7d+/K+++/L0uXLlXWjRw5UkqXLi0lSpQQGxsb+fDDD82yD0XNzZs3JSMjQ1n+4YcfpGnTpnLgwAG5cuWKeHt7y+effy46nU66d+8upUqVkpIlS0qZMmUkODjYrLUXBZmZmXLjxg3JzMxU2hYvXiwdO3aU48ePy7lz5+TVV1+V0NBQSUtLk5YtW4q1tbVYW1uLvb29rFy50qz1FwYM8hfg6NGjkpSUlOv6jh07ynfffSdvv/22lC9fXuzs7OT48eMiIrJkyRIZOHCgiIhs3rxZ/Pz8chzj77//FisrK3n8+PEL2ouibceOHXL8+HFp2rSpVKhQQSpXrixXrlwREZHu3bvL0KFDpW3btmJrayvNmjWTVatWZRvj999/l8qVK5uh+qJh1apVcvz4cXF3d5fy5cuLq6ur3L9/X0REGjZsKIGBgeLt7S22trbSpk0bCQsLyzbG9u3bpWnTpmaovnBhkD+np0+fyqNHj7K0tWzZMs8/uxcuXChly5aVtWvXiojI2rVrpV27diIicvHiRXFychIRkXv37omDg4MS1k+fPpWTJ0/KzJkzpU6dOjwi11Nqaqr89ttvcv78eaXto48+kooVK8q+fftEROTTTz+VgIAAERHp06ePtG3bVsLCwrL9Qn78+LH89ttvMnXqVHF2dpbPPvusgPfGMiUkJMihQ4fk2rVrSpufn5/UrFlTTp8+LSIiH3zwgcycOVNERLy8vKRXr16yd+9eSU9PzzJWenq6/PLLLxIYGChVqlSRdevWFfDeFD4M8udw8uRJcXR0FHt7e1m0aJGIiERFRcmwYcPy7Hfu3DkpVaqUspyRkSH29vbKFEv16tXl0qVLIiLi6+srvXr1En9/f/n222/Fw8NDxo0bJ0eOHHmh+1ZUbN26VapUqSItW7YUNzc30Wg0IiKya9cucXNzU7a7efOmcnS9detWcXV1lRs3bohOp5Nz587JwoUL5f79+xIRESGNGjUSjUajBBDl7auvvpJKlSrJm2++KTVr1lSmC5cuXSq+vr7KdocPH5YmTZqIiMi8efPEx8dH7t+/L5mZmXLy5EmZM2eOZGZmSmhoqLRo0UJmzJghly9fNtt+FSYMcj3Ex8fLpk2bZMmSJVnaBw0aJCtWrJDU1FSpVauWaLVamTVrluzevTvfMatUqSJXr15Vljt06CD/+9//RERkwIAByn/scXFxMn/+fNm8ebPyZydld/36dfnxxx9ly5YtMnToUJFnR4Hu7u5y7do1CQ8Pl6CgIClRooTEx8dLSkqK2NraSlpamjKGq6urnD17VnQ6nUycOFHKlCkj1tbWUqtWLRkyZIjExsaacQ/Ny8fHR86ePZvnNufPn5eff/5ZlixZItOnTxcRkUuXLknTpk0lJiZGdu/eLYMGDVL+yoyOjlb++pRnf+2o1WpJSEgQrVYrAwcOFBsbGylZsqTUqVNHxowZk+eU5cuMQZ6PS5cuiYODg7Rv314cHByynJAJDg6WmTNnyvnz56VKlSqSmpoqbdu2lYyMDNm3b58MHTpUvv/++xzHfe+992T58uXK8vz582XMmDEiz65s2bVrVwHsXdGwYsUKKVeunPTu3Vvc3d2VI+vffvtNypUrJ3Z2dtKoUSOZNGmSHDhwQJmqat68uezZs0cZZ+TIkbJw4UJlOSMjI9u02cvqww8/zPJv819Tp06VKlWqSJ8+fcTZ2Vk5st68ebNUqlRJypQpI82bN5fg4GA5cuSIPH36VEREnJycJDo6Whnn7bfflh07dijL6enpPA+kBwa5Hv4Jb0dHRzlx4oTSfv/+fXnzzTfllVdekeXLl0tcXJz07t1brl69Kk5OTvL1119L/fr1lROZ/7Zu3Trp0aOHsnz79m2JjIwsoD2yPCkpKRIXF5fjuooVK8off/whIiIPHz6UChUqyNmzZyU2NlZsbGzk5s2bOfYLDg6WoKAgZfnSpUtZQoX+z7Zt26Rjx445rktKSpLSpUsrP5/r169LqVKlJCEhQU6cOCEVK1aUxMTEHPv++69PEZHTp09LTEzMC9qLoosvljBA//794erqiqlTp+a4fsWKFbC1tYWLiwsWLVqE77//HnPnzoW1tTXGjh2bZdv79+/jjz/+QIcOHQqoess0a9YsNG7cGO+//z40Gk22f8fMzEyULl0ajx49QrFixQAAH3zwARo0aIBx48ahb9++uH37NoYOHYp79+4hPDwcjo6OWLZsGa5fv474+Hh4eXmZae8sR0JCApydnXH//n2ULFkyy7obN26gdevWuHHjhtLWqlUrjB07Fv7+/mjZsiXKli2LPn36ICYmBgcPHoSXlxc+/fRTnD17FlZWVqhXr54Z9qoIMfdvksLqwoULWa5yEBFZv369tGrVKtc+Xbp0kYcPH8qTJ0+kXr160qVLF1Gr1S/9SbHLly9Lv3798t3uxo0b8uOPP2Y5geXt7S3u7u5y69atXPvVr19fDhw4oCx37txZOnXqJPJs3nXOnDni7+8vI0aMkE2bNsndu3efe5+Ksj179oi/v7+0b98+y9Rg48aN5dChQ9m2f/LkiVSqVEm5fFOn04mnp6eMGDFC5NlfU5988on4+/vLRx99JKGhoZKQkFCAe1T0McifuXPnjmzcuFE++OADqVatmjg7O8s333yTZZvbt29LyZIlJTU1NVv/lJSULH96PnjwQHbs2CE3btwokPoLs4yMDLGzs8vzRFVgYKBUq1ZN/Pz85JVXXpGNGzeKPLssML9fAtu2bZNq1apJSEiIDBw4UHr16iVVq1Y1+X5Yug4dOsiff/6ZpU2r1cqRI0dkwIABotPpZPfu3eLo6CgbN25Urvj58ccfRUREo9HIlClTchz7888/Fzc3N5k/f7507dpVevXqJV5eXgWyXyTyUj5rZdOmTcpzMv7++294enqiTp06CA0NRZMmTRAeHo6IiAjY29vj3zNPjo6OqF27Ng4fPozz588jISFBWXfp0iWMGDFCWS5Xrhz8/f3h7OxcwHtXOLz//vu4cOEC8OyluF5eXggPD0dcXBxmzJiBkJAQpKenAwD27duH69ev49ixY+jZsyfc3d0xbdo0AEC7du0QHR2d52f16NEDK1euxM2bN9GwYUN88sknePr0aQHspWWpXbs2tm/fjrCwMOzYsQMrV65E165dMXr0aLRs2RJPnz7Fpk2boNFo0K9fP/Ts2RMhISFYunQpAKB9+/bYv39/jmOPHz8e06ZNw9WrV9GtWzcMGjSIP4OCZO7fJAUpMTFR+vTpI82bN5e//vpL5NmfgSdOnJDk5GTZvXu3jB8/Xjw9PaV8+fLSo0cP5SRNYmKirF+/XurUqSMlS5aUV155hddy52HYsGEyf/58iYiIkP3798u8efPk3XffFQ8PD9FoNPLmm2/KoEGDRJ5dsVOuXDlxcHAQPz8/WbJkiVy8eFHk2Z/tarVaHjx4kOfn6XQ6+eOPP2TmzJni6Ogoc+bMKZD9tBRPnjwRd3d3KV68uLRv3142btwoUVFRAiDLUfrgwYPliy++UJbPnDkjNWvWFHl29G5nZ5frtMiTJ0/k6NGjotFopFy5crJ58+YC2DOSoja1cv/+ffnjjz+yXCL4j19//VVeffVVmT17tnLp0z+Cg4PF3t5e2rRpI7Nnz5aTJ09m2+b69evSo0cPWb58eZbrvym7zMxM8fLyEpVKJfXq1ZMvvvhCTp8+LQCUK3gSExNFrVZLamqqhIWFibu7e44/NxGRd955R7Zu3ZrnZyYlJUm9evVk9OjRcurUqReyX5bir7/+kpiYGHn//felSZMmsmHDBhERCQ8PFzs7uyzPlalataryS1OePbqgbt26ys1pn3zyibz99tvK+kGDBilXCP3XhQsX5PXXX5eJEydmGZNePIsP8idPnsiECRNk8+bNUrVqVXF2dhZ3d3f5+++/lfVTpkwRDw+PXP8DvH37drbbgEl/ERERMmHCBJkyZYpyTiA0NFTs7e2Va4B1Op1UrFhRuWNVRKRt27by008/SWZmpnh4eEj37t1l1apVMnHiRGnYsKF8+eWXIs8uScvtEkLKKiUlRezs7KRdu3ayZs0aCQsLk3LlyimX9DVr1kwOHjyobP/ee+9leejU06dPpX///lKuXDlxdnYWT09PuX79uln2hfRn8XPkxYsXx86dO7FgwQJcunQJN27cQNu2bfHJJ58AAFavXo1Zs2bB2dkZq1atwvTp07Fs2TJs374dp0+fBp7NfdvY2Jh5TyzT5s2b0bVrVzg4OODhw4fw9vbG3bt30b17d7i5ueH48eMAAJVKhbZt22Z5Fvs/c67FihXD4cOH4eHhgV9//RU2NjZYsGABhg0bBgBo0KABatSoYbZ9LKw2btyIZs2aoU6dOpg3bx4AwNbWFg0aNICPjw8GDhwIPz8/9O7dG9u2bQOe/Zv/8ssvyhjt2rXDvn37lGUrKyusX78ekZGR2L9/P/7880+88sorZtg7MkRxcxdgCu3atYOVlRXs7OyAZydemjZtCgAYOHAgWrdujfj4eNy5cwfx8fG4desWIiMj8fjxY2zYsMHM1Vu2BQsWYNWqVejUqRMAICUlBStXrsSUKVOU0GjZsiXwLET27duHwYMHK8vBwcEAgLJly+LTTz81454UblFRUdi+fTsqVaqE4cOHY+fOnZg2bRpWrVqFMmXK4IMPPkD58uUxePBgtGvXDomJiUpfX19fLF++HIGBgejRowc6duyIsmXL4vLly5gxY0aO/w8wvC2Muf8kMIUdO3ZIixYtlOU7d+6Ira2tWWt6WXh4eEhERISyvHbtWunevbuIiERGRoqzs7OsX79eNBqNxMTEiLu7uxmrtUzDhg2T6tWry8cff6zcWTxgwIAsJyV37typPM712LFjUr9+fWVdUlKSqNVq0Wq1Is/uhxg3bpxs2LAh27kgskxF4oi8TZs26NmzJ6ZMmYLWrVtjwYIF8Pf3N3dZL4VWrVphw4YNaNKkCQAgMjISVatWBQA0bNgQo0ePxp49e/D666+jatWqiIqKMnPFluXAgQPYt28fLly4oPzFCQBqtRpJSUnKsqurK27evAkAaNKkCW7evIl79+6hYsWKsLe3x4cffogHDx6gatWq6N+/P/r372+W/aEXo8jcot+sWTO4u7vj3r17aNSoEcaPHw9bW1tzl1Xk3bp1Cz4+PnB1dcXjx49x69YtHD58GNWrVzd3aUXC119/jZMnT2LdunVZ2o8cOYI+ffrgyJEjqFatGoKCgvDXX39h586dAIDJkyejb9++vPX9JVEkjsjxbJ4cANasWWPuUl4q1atXx9mzZ7F3717Y29ujbdu2KFGihLnLKjLq1q2LJUuWQKfTwcrq/65NaN68OYYMGYK6detCRODt7Y1vv/1WWT979mwzVUzmUGSOyA8dOoRJkybhxIkT5i6FyKTatGkDa2tr+Pr6Ijo6GhcvXsSvv/4KAEhLS4NKpULp0qXNXSaZUZEJ8idPnuDy5cuoW7euuUshMqn09HSsWLEC165dg4uLC9q2bQtPT09zl0WFSJEJciKil5XF3xBERPSyY5ATEVk4BjkRkYVjkBMRWTgGORGRhWOQExFZOAY5EZGFY5ATEVk4BjkRkYVjkBMRWTgGORGRhWOQExFZuP8Hw7K1zP6AJ4AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 380x320 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fractions = [np.mean(membership[hit]), np.mean(membership[target & ~hit]), np.mean(membership[is_control])]\n",
    "plt.figure(figsize=(3.8, 3.2))\n",
    "plt.bar([\"p<=1% hits\", \"other genes\", \"control genes\"], fractions, color=[\"#2b7\", \"#9cf\", \"#888\"])\n",
    "plt.ylabel(\"fraction in a CORUM complex\")\n",
    "plt.xticks(rotation=15, ha=\"right\", fontsize=8)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "corum-read",
   "metadata": {},
   "source": [
    "Three fifths of the p-value hits belong to a known complex, against a quarter of the other tested genes, an odds ratio of 4.7. The honest q-value list is cleaner still: three in four of those genes sit in a complex. The top hits include proteasome, mRNA export and RNA polymerase genes, machines a cell cannot lose without changing shape."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "another-line",
   "metadata": {},
   "source": [
    "## Another cell line\n",
    "\n",
    "Only one line in the whole analysis names the cell line: the DepMap model id passed to {func}`~mantispy.io.unexpressed_genes`. For a screen run in a different line, change the model id, for example `\"ACH-001086\"` for HeLa. Nothing else changes, and nothing is hard-coded to U2OS."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "summary",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- Score each gene by how many features separate its wells from the negative controls, a Mann-Whitney test per feature.\n",
    "- Read the genes the screened cell line does not express from DepMap with {func}`~mantispy.io.unexpressed_genes`, keyed by the line's DepMap model id.\n",
    "- {func}`~mantispy.tl.empirical_fdr` turns the score and the null into a per-gene p-value, the control-gene rate PERISCOPE thresholds, and a target-decoy q-value, the honest false discovery rate.\n",
    "- On JUMP CRISPR the hit list is enriched for CORUM complex members, with an odds ratio of 4.7 at the p-value cutoff and 10.4 at the q-value cutoff."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.13.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}