notebooks/melanoma_sccoda.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8ea6915c",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Differential abundance between responder and non-responder lesions\n",
    "***\n",
    "Jan T. Schleicher, 2022\n",
    "\n",
    "In this notebook, we use scCODA and Wilcoxon rank-sum tests as benchmark methods to detect differentially abundant cell types between melanoma responder and non-responder lesions. This analysis is based on the data from [Sade-Feldman et al. 2018](https://doi.org/10.1016/j.cell.2018.10.038)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "98703599",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5be811c6",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import anndata as ad\n",
    "import pickle as pkl\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "\n",
    "from sccoda.util import comp_ana as mod\n",
    "from sccoda.util import cell_composition_data as dat\n",
    "from sccoda.util import data_visualization as viz\n",
    "\n",
    "%config Completer.use_jedi = False\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "68d77a36",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "out_dir = \"../output/differential_abundance\"\n",
    "os.makedirs(out_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff9e66e4",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Loading data\n",
    "\n",
    "We load the data in the form of a scanpy object which we created from the Seurat object. From this, we create the covariates data frame and compositional AnnData object required for scCODA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dafbc39e",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ubuntu/miniconda3/envs/sccoda/lib/python3.8/site-packages/sccoda/util/cell_composition_data.py:256: FutureWarning: X.dtype being converted to np.float32 from float64. In the next version of anndata (0.9) conversion will not be automatic. Pass dtype explicitly to avoid this warning. Pass `AnnData(X, dtype=X.dtype, ...)` to get the future behavour.\n",
      "  return ad.AnnData(X=count_data.values,\n"
     ]
    }
   ],
   "source": [
    "# read data\n",
    "adata = ad.read_h5ad(\"../data/melanoma_seurat.h5ad\")\n",
    "\n",
    "# create covariates dataframe\n",
    "cov_df = pd.DataFrame({\"Sample\": adata.obs.sample_name, \"Cond\": adata.obs.response_status})\\\n",
    "    .drop_duplicates(subset=\"Sample\")\\\n",
    "    .set_index(\"Sample\")\n",
    "cov_df.Cond = pd.Categorical(cov_df.Cond, categories=[\"R\", \"NR\"], ordered=True)\n",
    "\n",
    "# create compositional data from scanpy data\n",
    "comp_data = dat.from_scanpy(adata,\n",
    "                            cell_type_identifier=\"cluster\",\n",
    "                            sample_identifier=\"sample_name\",\n",
    "                            covariate_df=cov_df)\n",
    "comp_data.write(os.path.join(out_dir, \"melanoma_sccoda.h5ad\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "170e452e",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## scCODA\n",
    "\n",
    "We set up the scCODA model for compositional analysis from the compositional data, using automatic reference cell type detection. Then, we run MCMC and look at the results, setting the estimated false discovery rate to 0.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7c2fd9a1",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Automatic reference selection! Reference cell type set to 6\n",
      "Zero counts encountered in data! Added a pseudocount of 0.5.\n"
     ]
    }
   ],
   "source": [
    "# model setup and inference\n",
    "model = mod.CompositionalAnalysis(comp_data, formula=\"Cond\", reference_cell_type=\"automatic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e2799027",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████| 20000/20000 [01:59<00:00, 168.02it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MCMC sampling finished. (151.622 sec)\n",
      "Acceptance rate: 58.0%\n"
     ]
    }
   ],
   "source": [
    "# run MCMC\n",
    "sim_results = model.sample_hmc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c9429658",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compositional Analysis summary:\n",
      "\n",
      "Data: 42 samples, 11 cell types\n",
      "Reference index: 6\n",
      "Formula: Cond\n",
      "\n",
      "Intercepts:\n",
      "           Final Parameter  Expected Sample\n",
      "Cell Type                                  \n",
      "0                   -0.191        40.943199\n",
      "1                   -1.210        14.778682\n",
      "2                   -0.698        24.660079\n",
      "3                   -1.067        17.050606\n",
      "4                   -0.237        39.102473\n",
      "5                   -0.655        25.743591\n",
      "6                    0.036        51.376749\n",
      "7                   -0.158        42.316866\n",
      "8                   -0.567        28.111696\n",
      "9                   -0.009        49.116043\n",
      "10                  -0.732        23.835730\n",
      "\n",
      "\n",
      "Effects:\n",
      "                      Final Parameter  Expected Sample  log2-fold change\n",
      "Covariate  Cell Type                                                    \n",
      "Cond[T.NR] 0                -0.772828        17.498119         -1.226424\n",
      "           1                 0.000000        13.679821         -0.111468\n",
      "           2                 0.000000        22.826494         -0.111468\n",
      "           3                 0.000000        15.782818         -0.111468\n",
      "           4                 0.000000        36.195032         -0.111468\n",
      "           5                 0.739857        49.937856          0.955921\n",
      "           6                 0.000000        47.556661         -0.111468\n",
      "           7                 0.000000        39.170420         -0.111468\n",
      "           8                 0.000000        26.021467         -0.111468\n",
      "           9                 0.000000        45.464048         -0.111468\n",
      "           10                0.665019        42.902980          0.847952\n"
     ]
    }
   ],
   "source": [
    "sim_results.save(os.path.join(out_dir, \"melanoma_sccoda_results.pkl\"))\n",
    "sim_results.set_fdr(est_fdr=0.1)\n",
    "sim_results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1396a79c",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We extract the credible effects from the simulation results and create a data frame containing the cell type proportions, required for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "def8873c",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "credible_effects = sim_results.credible_effects()\n",
    "sccoda_relevant_indices = np.where(credible_effects)[0]\n",
    "cell_types = credible_effects.index.get_level_values(1).to_list()\n",
    "\n",
    "prop_df = pd.DataFrame(comp_data.X/np.sum(comp_data.X, axis=1, keepdims=True),\n",
    "                       columns=comp_data.var.index, index=comp_data.obs.index)\\\n",
    "    .merge(comp_data.obs.Cond, left_index=True, right_index=True)\\\n",
    "    .melt(id_vars=\"Cond\", var_name=\"cell type\", value_name=\"proportion\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5fdad2",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Wilcoxon rank-sum tests\n",
    "\n",
    "As another basic method, we use Wilcoxon rank-sum tests to identify differentially abundant cell types between the two conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ed4de79f",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>response_status</th>\n",
       "      <th>sample_name</th>\n",
       "      <th>cluster</th>\n",
       "      <th>frequency</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NR</td>\n",
       "      <td>Post_P10</td>\n",
       "      <td>G1</td>\n",
       "      <td>0.002525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NR</td>\n",
       "      <td>Post_P11</td>\n",
       "      <td>G1</td>\n",
       "      <td>0.003195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NR</td>\n",
       "      <td>Post_P12</td>\n",
       "      <td>G1</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NR</td>\n",
       "      <td>Post_P13</td>\n",
       "      <td>G1</td>\n",
       "      <td>0.253898</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NR</td>\n",
       "      <td>Post_P14</td>\n",
       "      <td>G1</td>\n",
       "      <td>0.057751</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>457</th>\n",
       "      <td>R</td>\n",
       "      <td>Pre_P26</td>\n",
       "      <td>G9</td>\n",
       "      <td>0.006006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>458</th>\n",
       "      <td>R</td>\n",
       "      <td>Pre_P28</td>\n",
       "      <td>G9</td>\n",
       "      <td>0.026706</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>459</th>\n",
       "      <td>R</td>\n",
       "      <td>Pre_P33</td>\n",
       "      <td>G9</td>\n",
       "      <td>0.002212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>460</th>\n",
       "      <td>R</td>\n",
       "      <td>Pre_P7</td>\n",
       "      <td>G9</td>\n",
       "      <td>0.073370</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>461</th>\n",
       "      <td>R</td>\n",
       "      <td>Pre_P8</td>\n",
       "      <td>G9</td>\n",
       "      <td>0.014981</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>462 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    response_status sample_name cluster  frequency\n",
       "0                NR    Post_P10      G1   0.002525\n",
       "1                NR    Post_P11      G1   0.003195\n",
       "2                NR    Post_P12      G1   0.000000\n",
       "3                NR    Post_P13      G1   0.253898\n",
       "4                NR    Post_P14      G1   0.057751\n",
       "..              ...         ...     ...        ...\n",
       "457               R     Pre_P26      G9   0.006006\n",
       "458               R     Pre_P28      G9   0.026706\n",
       "459               R     Pre_P33      G9   0.002212\n",
       "460               R      Pre_P7      G9   0.073370\n",
       "461               R      Pre_P8      G9   0.014981\n",
       "\n",
       "[462 rows x 4 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# perform testing for cluster frequencies\n",
    "meta_data = pd.read_csv(\"../data/melanoma_umap.txt\", sep=\"\\t\", index_col=0)\n",
    "cell_type_frequencies = meta_data.groupby([\"response_status\", \"sample_name\"], observed=True)[\"cluster\"]\\\n",
    "    .value_counts(normalize=True).rename(\"frequency\").unstack(fill_value=0).reset_index()\\\n",
    "    .melt(id_vars=[\"response_status\", \"sample_name\"], value_name=\"frequency\")\n",
    "cell_type_frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b789f2b4",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cluster_index</th>\n",
       "      <th>cluster</th>\n",
       "      <th>log2fc</th>\n",
       "      <th>p</th>\n",
       "      <th>padj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>G1</td>\n",
       "      <td>-2.158012</td>\n",
       "      <td>0.002407</td>\n",
       "      <td>0.013238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>G2</td>\n",
       "      <td>0.578780</td>\n",
       "      <td>0.780688</td>\n",
       "      <td>0.780688</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>G3</td>\n",
       "      <td>2.065498</td>\n",
       "      <td>0.010422</td>\n",
       "      <td>0.028662</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>G4</td>\n",
       "      <td>1.278635</td>\n",
       "      <td>0.089436</td>\n",
       "      <td>0.157844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>G5</td>\n",
       "      <td>-1.032665</td>\n",
       "      <td>0.100446</td>\n",
       "      <td>0.157844</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>G6</td>\n",
       "      <td>1.962359</td>\n",
       "      <td>0.005844</td>\n",
       "      <td>0.021429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>6</td>\n",
       "      <td>G7</td>\n",
       "      <td>0.252103</td>\n",
       "      <td>0.290052</td>\n",
       "      <td>0.398822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>7</td>\n",
       "      <td>G8</td>\n",
       "      <td>0.480612</td>\n",
       "      <td>0.372954</td>\n",
       "      <td>0.455833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>8</td>\n",
       "      <td>G9</td>\n",
       "      <td>0.283632</td>\n",
       "      <td>0.635986</td>\n",
       "      <td>0.699584</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>9</td>\n",
       "      <td>G10</td>\n",
       "      <td>-1.860139</td>\n",
       "      <td>0.042104</td>\n",
       "      <td>0.092629</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>10</td>\n",
       "      <td>G11</td>\n",
       "      <td>1.830245</td>\n",
       "      <td>0.000238</td>\n",
       "      <td>0.002614</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    cluster_index cluster    log2fc         p      padj\n",
       "0               0      G1 -2.158012  0.002407  0.013238\n",
       "1               1      G2  0.578780  0.780688  0.780688\n",
       "2               2      G3  2.065498  0.010422  0.028662\n",
       "3               3      G4  1.278635  0.089436  0.157844\n",
       "4               4      G5 -1.032665  0.100446  0.157844\n",
       "5               5      G6  1.962359  0.005844  0.021429\n",
       "6               6      G7  0.252103  0.290052  0.398822\n",
       "7               7      G8  0.480612  0.372954  0.455833\n",
       "8               8      G9  0.283632  0.635986  0.699584\n",
       "9               9     G10 -1.860139  0.042104  0.092629\n",
       "10             10     G11  1.830245  0.000238  0.002614"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clusters = [\"G\" + str(i) for i in range(1, 12)]\n",
    "differential_abundance_dict = {\"cluster_index\": [], \"cluster\": [], \"log2fc\": [], \"p\": []}\n",
    "\n",
    "for i, c in enumerate(clusters):\n",
    "    df = cell_type_frequencies[cell_type_frequencies.cluster == c]\n",
    "    # Wilcoxon test\n",
    "    _, p = stats.ranksums(*df.groupby('response_status')['frequency'].apply(lambda x: x.values))\n",
    "    differential_abundance_dict[\"cluster_index\"].append(i)\n",
    "    differential_abundance_dict[\"cluster\"].append(c)\n",
    "    differential_abundance_dict[\"p\"].append(p)\n",
    "    # log2fc\n",
    "    differential_abundance_dict[\"log2fc\"].append(np.log2(df.loc[df.response_status == \"NR\", \"frequency\"].mean() /\n",
    "                                                         df.loc[df.response_status == \"R\", \"frequency\"].mean()))\n",
    "    \n",
    "differential_abundance = pd.DataFrame(differential_abundance_dict)\n",
    "differential_abundance[\"padj\"] = multipletests(differential_abundance.p, method=\"fdr_bh\")[1]\n",
    "\n",
    "significant_indices = differential_abundance.loc[differential_abundance.padj < 0.1, \"cluster_index\"]\n",
    "\n",
    "differential_abundance.to_csv(os.path.join(out_dir, \"melanoma_differential_abundance_wilcoxon.csv\"))\n",
    "differential_abundance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b08686c1",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Plotting of results\n",
    "\n",
    "Finally, we create grouped boxplots to show the results of scCODA and the Wilcoxon rank-sum tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2075fd0c",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib.lines import Line2D\n",
    "sns.set_theme(style=\"ticks\", font_scale=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6d261bd6",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkcAAAGACAYAAACnRsZmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXy0lEQVR4nO3dd1zV1f8H8Ndl47gqCYYDEwxcoOJES1RygubIkYqYJioopZIpOTItTM2BmJp75kopxNyr/DoSNM0ZYGriQBFQNtzP7w9/3O4FLtw94PV8PHgI57Pe9wr3vu857885IkEQBBARERERAMDM0AEQERERGRMmR0REREQymBwRERERyWByRERERCSDyRERERGRDCZHRERERDKYHBERERHJYHJEREREJIPJkYwRI0ZgxIgRhg6DiIiIDMjC0AEYk0ePHhk6BCIiIjIw9hwRERERyWByRERERCSDyRERERGRDCZHRERERDKYHBERERHJYHJEREREJIPJEREREZEMJkdEREREMpgcEREREclgckREREQkg8kRERERkQyurUZal5GToXCbuZk5bCxtSt23snVlncRFRESkDCZHpHVVJlZRuK23e2/EhMRIf3aY4oDM3Ey5fYS1gs5iIyIiKguH1YiIiIhksOeItO5V5CuF28zNzOV+frrkqa7DISIiUonBk6N79+5h/fr1+PPPP/H333/D2dkZBw4cUOrY/fv3Y82aNXj48CHq16+P4OBg9OrVS8cRU1lUqRlifRERERkbgw+r/f333zh9+jTq168PFxcXpY87dOgQpk+fjm7dumHt2rVo3749Jk+ejN9//12H0RIREVF5JxIEwaDVrxKJBGZmr3O06dOn46+//lKq56hXr15wdXXF8uXLpW1jxozBy5cvsXv3brVi8fHxAQAcP35creOJiIjI9Bm856gwMVLFgwcPkJiYCD8/P7l2Pz8/XL16FSkpKdoKj4iIiCoYgydH6khMTAQAODs7y7W7uLhAEATpdiIiIiJVGbwgWx1paWkAALFYLNderVo1ue0lKRw6K8mjR4/g6OiohQiJiIjIVJlkz1EhkUgk93Nh+VTRdiIiIiJlmWTPkWwPUc2aNaXt6enpAIr3KMkqrdi6tF4lIiIiqhhMsueosNaoaG1RQkICRCJRsVokIiIiImWZZHJUr149ODs74+DBg3LtBw4cgIeHB+zs7AwUGREREZk6gw+rZWVl4fTp0wCAhw8f4tWrVzh06BAAoG3btrCzs0NYWBiioqJw48YN6XEhISGYPHkynJyc0KFDBxw/fhxnz57FunXrDPI4iIio4sjKyoKtra2hwyAdMXhy9Pz5c3zyySdybYU/b9myBe3atYNEIkFBQYHcPr169UJ2djZWr16N9evXo379+li6dCneeecdvcVORETl34oVKxAZGYl9+/ZhzZo1OHfuHKytrfHbb79hx44d2L17N+7evQtra2t4eXnhs88+Q7169aTH37hxA8uXL8fVq1fx8uVLVK9eHU2aNMFXX32FN998EwDg5uaG4cOHw9XVFRs3bsTDhw/h5OSE4OBg+Pr6ysVz584dLF26FJcuXUJ2djacnZ0xatQo9O/fX7rPhQsXMHLkSHz33Xf4+++/8dNPPyEzMxMeHh6YPXu2XPmJMvEJgqDUYy0vDD5DtjHhDNlERFRUYXJUp04d9O7dGx06dEBmZiZOnz6N/fv3w9/fHx06dEBaWhpWrlyJ9PR0/Pzzz6hZsyYyMzPRpUsX1K1bFx9//DFq1qyJ5ORkafJSuGyWm5sbHB0dUblyZQQFBcHW1hY7duzAb7/9huXLl6Nnz54AXtfafvDBB3jjjTcQHByMGjVq4JdffsGBAwcQGhqKsWPHAvgvOapTpw48PT3Rp08fvHr1CosXL4aVlRUOHjwIc3NzpeObNWtWmY+1XBFIqmvXrkLXrl0NHQYRERmRiIgIwdXVVVi+fLm07fLly4Krq6uwYcMGuX0fPXokeHh4CAsXLhQEQRCuXbsmuLq6CkePHi31Gq6uroKHh4eQnJwsbcvPzxd69uwpdOvWTdo2efJkoVmzZkJSUpLc8R9//LHQvHlzIT09XRAEQTh//rzg6uoqjB07Vm6/gwcPCq6ursLly5eVjk/Zx1qemGRBNhERkb51795d+v3JkychEonQt29f5OfnS79q1qyJRo0a4eLFiwCA+vXro1q1ali8eDF+/PFHxMfHKzy/l5eXXA+Mubk5evfujXv37uHx48cAgPPnz8PLy6vYhMX9+/dHVlYWLl++LNfetWtXuZ/d3NwAAElJSUrHp+xjLU8MXnNERERkChwcHKTfP3/+HIIgoEOHDiXuW1iHU7VqVWzduhWrV6/G0qVLkZaWBnt7ewwePBgTJkyApaWl9JiShqYK21JTU/Hmm28iNTUV9vb2CmNLTU2Va69evbrcz1ZWVgCA7OxspeNT9rGWJ0yOiIiIVFSjRg2IRCJs375dmnDIkm1zc3PD0qVLIQgCbt++jX379mHlypWwsbFBYGCgdL9nz54VO09hW2GSU716dSQnJxfb7+nTp9K4VFVWfKo81vKCw2pkFGJjYxEQEIDY2FhDh0JEVKbOnTtDEAQ8efIE7u7uxb4Kh69kiUQiNGrUCGFhYRCLxbh+/brc9nPnzsklSAUFBTh48CCcnJykd415eXnh/PnzePLkidyxP//8M2xtbdGiRQu1H5Oi+NR5rKaOPUdkFGJiYpCXl4eYmBi0atXK0OEQEZWqVatWGDJkCMLCwvDXX3+hTZs2sLW1RXJyMmJjY+Hq6ophw4bh5MmT2LFjB9577z3Uq1cPgiDgyJEjSE9PR8eOHeXOWaNGDQQEBMjdrZaYmIilS5dK9wkODsbJkycxcuRIBAcHo1q1aoiOjsapU6fw2WefoWrVqio9DmXiU/axlidMjsgo+Pr6IiEhodh8HkRExuqrr75C8+bNsWvXLvz444+QSCRwcHCAp6cnPDw8ALwueBaLxVi3bh2ePn0KS0tLNGjQAAsWLJCblwh4XTzdsGFDLFu2DI8ePUK9evWwePFi9O7dW7qPs7Mzdu7ciSVLluCrr75CdnY2XFxcEB4ejgEDBqj8GJSNT5nHWp5wniMZnOeIiIgMoXASyNmzZxs6FAJrjoiIiIjkMDkiIiIiksGaIyIiIgO7ffu2oUMgGew5IiIiIpLB5IiIiIhIBpMjIiIiDUVHR2PTpk0anWP79u3Yt2+fdgIijTA5IiIi0tCBAwewZcsWjc7x448/Yv/+/VqKiDTB5IiIiIhIBpMjIiKiMqSkpGDWrFnw9vZGs2bN0L59ewwdOhT/+9//4O/vj1OnTuHhw4dwc3OTfhWKjIzEoEGD0LZtW3h6eqJ///7Ys2cPZOdg7tq1K/7++29cvHhRenzXrl0BAPv27YObmxv+/fdfuZguXLgANzc3XLhwQdp248YNjBs3Dl5eXmjWrBneeecdBAYG4vHjxzp+hsoX3spPRERUhs8++ww3btzA5MmT8dZbbyE9PR03btxAamoq5syZg1mzZuHBgweIjIwsduzDhw8xZMgQ1K5dGwBw5coVzJ8/H0+ePMHEiRMBvE6gQkJCULVqVcyZMweA6qvdZ2Zm4qOPPkLdunUxe/Zs1KxZE8nJybhw4QIyMjI0fAYqFiZHREREZYiLi8OgQYMwePBgadt7770n/V4sFsPKygotWrQodmx4eLj0e4lEgrZt20IQBGzZsgXBwcEQiURo0qQJbGxsUKVKlRLPoYzExESkpqbi66+/lotNdm02Ug6TIyIiojJ4eHhg//79qF69Ojp06ICmTZvC0tJSqWPPnTuHNWvW4Nq1a3j16pXctufPn6NmzZpaibF+/fqoVq0aFi9ejOTkZLRp0wYNGzbUyrkrGtYcERERlWHp0qXo168f9u7diyFDhqBt27aYNm0akpOTSz3u6tWrGDNmDABg3rx5+PHHH7F3716MHz8eAJCdna21GKtWrYqtW7eicePGWLp0KXx9ffHOO+8gIiICeXl5WrtORcCeIyIiojLY2dnhiy++wBdffIGkpCScOHEC3333HZ4/f47169crPC4mJgYWFhZYs2YNrK2tpe3Hjh1T+tqFx+Xm5sq1v3jxoti+bm5uWLp0KQRBwO3bt7Fv3z6sXLkSNjY2CAwMVPqaFR17joiIiFRQu3ZtjBgxAh06dMCNGzcAvC6eLqkXSCQSwdzcHGZm/73dZmdn45dffim2r6Jz1KlTB0Dx9ddOnDihMEaRSIRGjRohLCwMYrEY169fV+7BEQD2HBEREZXq5cuXGDlyJPz8/ODs7IzKlSvj2rVr+O2339CtWzcAgKurK44cOYIdO3agWbNmEIlEcHd3h7e3NzZu3IipU6diyJAhSE1Nxfr160u8E83V1RUxMTE4ePAg6tatC2tra7i5ucHd3R0NGjTAwoULUVBQALFYjGPHjiE2Nlbu+JMnT2LHjh147733UK9ePQiCgCNHjiA9PR0dO3bUy3NVXogE2YkWKjgfHx8AwPHjxw0cCRERGYvc3Fx8/fXXuHz5Mh4+fIj8/Hw4OjrC19cXH3/8MWxtbZGWloZZs2bh3LlzePnypXRYCwB++uknrF27Fg8fPkStWrUwePBg6TDd8ePHUbduXQCvb/mfPXs2Ll++jIyMDNSpU0faO/TPP/9g3rx5uHLlCqysrODr64t3330XgYGB2LJlC9q1a4fExESsXLkSV69exdOnT2FpaYkGDRpg2LBh6N+/v8GeP1PE5EgGkyMiIiJizRERERGRDCZHRERERDKYHBERERHJYHJEREREJIPJEREREZEMJkdEREREMpgcEREREclgckREREQkg8kRERERkQwmR0REROXYihUr4ObmJv1q164dPvzwQ5w+fdrQoRktLjxLRESkAolEAolEYpBrm5mZwcxM9X4NGxsbbN68GQDw9OlT/PDDDxg/fjy2b98OT09PbYdp8pgcERERKUkikSBo4kSkp6Ya5Pri6tXxfWSkygmSmZkZWrRoIf25RYsW6NSpE6KiopgclYDJERERkZIkEgnSU1Px5ofdADORni8u4PGPRyGRSNTqPZLl4OAAOzs7JCUlaSm48oXJERERkarMRBBpmKCoSoD2hvIyMjKQlpYGJycnrZ2zPGFyREREVAHk5+cDAJKTk7F48WJUqVIFI0eONHBUxonJERERUTmXmZmJpk2bSn82NzfHqlWr8NZbbxkuKCPGW/mJiIjKORsbG+zduxd79uzBokWLYG9vj2nTpuHp06eGDs0oGUVydPfuXYwZMwYtWrSAl5cX5s+fj+zs7DKPy8zMxOLFi/Hee++hefPm6N69O1asWIHc3Fw9RE1ERGQazMzM4O7uDg8PD/Tt2xcrV65Eeno6Vq5caejQjJLBh9XS09MREBCA2rVrIyIiAikpKQgPD0dqaioWL15c6rFffvkljh07hsmTJ+Ptt9/G1atXERERgbS0NMycOVNPj4CIiMi0NGvWDL6+vti3bx8mTpwIe3t7Q4dkVAzec7Rz506kp6fj+++/R6dOndCvXz/MnDkT0dHRSEhIUHhcfn4+Dh06hFGjRsHf3x/t27dHYGAgPvjgAxw8eFCPj4CIiMj0BAUFIT8/Xzo5JP3H4MnRmTNn4OXlBTs7O2lbjx49YGVlVerU5oIgoKCgAFWrVpVrF4vFEARBZ/ESERFBIkCQSPT6BYl239ucnZ3h6+uLH3/8ES9fvtTquU2dwYfVEhISMHDgQLk2KysrODk5ldpzZGlpiQEDBmDr1q3w9PREw4YNce3aNezevRsjRozQddhERFQBmZmZQVy9Oh7/eNQg1xdXr67yBJCTJk3CpEmTStxWVvlKRWXw5Cg9PR1isbhYu1gsRlpaWqnHfvnll5gzZw4GDx4sbfP398fEiRMVHuPj46Nw26NHj+Do6KhE1EREVBGZmZnh+8hIk1tbjVRj8ORIEUEQIBKVPjX74sWLcerUKcybNw8NGjTA9evXERERAbFYjJCQED1FSkREFQkTlPLP4MmRWCxGenp6sfaXL1/CxcVF4XF37tzBhg0b8P3330t7g9q0aQORSISFCxdi+PDheOONN4odd/z4cYXnLK1XiYiIiCoGg6e+Li4uxWqLcnNzcf/+/VKTo/j4eABA48aN5dobN26M/Px8PHz4UPvBEhERUbln8OSoU6dOOH/+PF68eCFtO3r0KHJzc+Ht7a3wuDp16gAArl+/Ltf+119/AQDq1q2rg2iJiIiovDP4sNrQoUOxbds2BAUFISgoCM+fP8eCBQvQp08fuZ6jsLAwREVF4caNGwBeT2Dl4eGBOXPm4NmzZ2jQoAGuXbuG77//Hr1795abGoCIiIhIWQZPjsRiMTZv3oz58+dj0qRJsLGxgZ+fH0JDQ+X2k0gkKCgokP5sbm6O1atXY/ny5Vi7di2ePXsGR0dHjBgxAuPHj9f3wyAiIqJyQiRwxkSpwoLs0oq2iYiIqHwzeM0RERERkTFhckREREQkw+A1R0RERKZEIpGY1AzZK1asQGRkJFq3bo3t27cX27ZhwwZcvnwZAODm5ibdZmlpCUdHR3h7e2PixImoXr26xvGbCiZHRERESpJIJAgJDkJKWvHJi/XBrpoYESu/V2uG7kuXLuHcuXPw8vIqdT9/f3/4+fkhNzcXFy9exOrVq/HPP/9g3bp16oZtcpgcERERKUkikSAlLR0LO9SGeRlLXGlbgSBg2v+SIJFIVE6OKlWqhLfffhsrV64sMzlydHREixYtAABt27bF06dPsWvXLjx9+hQODg7qhm9SWHNERESkInORCOZmev7SMBkLDg7GH3/8gQsXLqh0XKNGjQC8Xpy9omByREREVAF4e3vD3d0dkZGRKh336NEjmJmZSVemqAiYHBEREVUQwcHBuHjxIi5evKhwH4lEgvz8fGRmZuL06dPYuXMnhg8fjpo1a+oxUsNickRERFRBdOnSBU2bNi2192jx4sVo2rQpWrZsicDAQLz99tuYPn26HqM0PCZHREREFUhQUBAuXLiAS5culbh95MiR2Lt3L7Zu3YpBgwYhNjYWy5cv13OUhsW71YiIiCqQ9957D40bN0ZkZCRatWpVbPubb74Jd3d3AK/vVktJScHGjRsxbNgwODo66jtcg2DPEVE5lZGTofArOy9bqX2JqHwKDg7GuXPnEBsbW+a+oaGhKCgowPr16/UQmXFgzxFROVVlYhWF23q790ZMSIz0Z4cpDsjMzSy2n7CW61ITlUfvvfce3NzccO7cOVSqVKnUfZ2dndG7d2/s3bsXwcHBqFGjhp6iNBz2HBEREamoQBBQINHzl6C9DysikQjBwcFK7x8cHIzc3Fxs27ZNazEYM5EgaPHZNnE+Pj4AgOPHjxs4EiLNlTYsZm5mDhtLmzL3rWxdWetxEZkyU14+hJTHYTWickqVxIZJEJFyzMzMELHye5NaeJZUx+SIiIhIBUxQyj8mR0TlkCAIyMrK0vg8tra2EOl5cU0iIkNjckRkRGJjYxEREYGQkJAS5x9RhiAI6Nevn8IJ3lTRpk0b7N+/nwkSEVUo7BckMiIxMTHIy8tDTExM2TuXgskMEZH62HNEZER8fX2RkJAAX19ftc8hEomwf/9+DqsREamJt/LL4K38RERExGE1IiIiIhlMjoiIiIhkMDkiIiIiksHkiIiIiEgGkyMiIiIiGUyOiIiIiGQwOSIiIiKSweSISEdiY2MREBCA2NhYQ4dCREQqYHJEpCPaWgqEiIj0i8uHEOmINpYC0RdBEJCTkyP9HvhvfTZra2suIUJEFQqXD5HB5UOoIhIEAXPnzsWdO3dK3O7q6oo5c+YwQSKiCkPtnqNXr17hzJkzSEpKQnZ2ttw2kUiE4OBgjYMjIiIi0je1eo7+/PNPBAYGIi0treSTikS4efOmxsHpG3uOqKIqHFbLycnBhAkTAACrVq2CtbU1h9WIqMJRq+coPDwctWrVwrp16+Dm5gYrKyttx2WSMnIyFG4zNzOHjaVNqftWtq6sk7iIyiISiWBjYyPXZm1tXayNiKgiUCs5un37Nr777ju4u7trOx6TVmViFYXberv3RkzIf3ctOUxxQGZuptw+wlqWfxERERmaWrfy29nZaTsOIiIiIqOgVs+Rv78/du7ciS5durAWQcaryFcKt5mbmcv9/HTJU12HQ0RERGpQKzmSSCRITExEv3790LlzZ1SvXl1uu0gkwqhRo7QQnmlRpWaI9UVERETGSa3kaOHChdLvb9++XWx7RU2OiIiIyPSplRzxVvfiBEFAVlaWRuewtbXlMCUREZGBqZUc1alTR6tB3L17F/Pnz0dsbCxsbW3h6+uL0NBQpW4jTk1NxbJly3Ds2DGkpaWhdu3a+OijjzB06FCtxlgaQRDQr18/XLp0SaPztGnTBvv372eCREREZEAara127949nD9/HqmpqahRowbatWuH+vXrq3SO9PR0BAQEoHbt2oiIiEBKSgrCw8ORmpqKxYsXl3psRkYG/P39YW1tjbCwMLzxxhu4d+8e8vLyNHlYamFCQ0REVD6olRwJgoB58+Zh586dkEgk0nYzMzMMGzYMM2fOVPpcO3fuRHp6OqKioqRTBJibmyM0NBQTJkyAi4uLwmPXrFmD7Oxs7NmzR9rL1K5dO3UekkZEIhH279/PYTUiIqJyQK3kaNOmTdixYwc+/PBD9O/fH7Vq1cKTJ08QFRWFHTt2oG7dukoXZJ85cwZeXl5ycyf16NEDYWFhOH36dKnJ0U8//YSRI0caxSy+IpEIlSpVMnQYREREpCG1kqM9e/ZgxIgRcj1EtWrVgoeHB8zMzLB7926lk6OEhAQMHDhQrs3KygpOTk5ISEhQeNyDBw/w7NkziMVijBs3DmfPnkXlypXRu3dvfP7550aRMBERGQsub0SkPLWSowcPHqBLly4lbuvSpQt27dql9LnS09MhFouLtYvFYoUL2wLAs2fPALyeVqBnz55Yu3Yt4uPjsWTJEuTl5WH+/PklHle4uGxJHj16BEdHR6VjJyIyFVzeiEh5aiVHVatWRVJSUonbkpKSUKWK4j9CZQmCUGr9TWGtk4uLC8LDwwEAXl5eyM/Px8KFC/HJJ5/A3t5e4ziIiIioYlErOerQoQOWLVuGxo0bo1mzZtL2mzdvIiIiAu+8847S5xKLxUhPTy/W/vLly1LrjQpn5W7fvr1ce/v27SGRSJCQkFBiclTaHE2l9SoREZkyLm9EpDy1kqOpU6diyJAhGDRoEBo2bAh7e3skJycjPj4eDg4OmDp1qtLncnFxKVZblJubi/v37xerRZJVr149WFpaFmsXhNddv2Zmaq2pS0RULnF5IyLlqZVBODo6IioqCh9//DFsbW3x77//wtbWFoGBgdi/fz/efPNNpc/VqVMnnD9/Hi9evJC2HT16FLm5ufD29lZ4nJWVFTp27Ihz587JtZ87dw4WFhZo2LCh6g+MiIiIKjy1J4G0s7NTqYdIkaFDh2Lbtm0ICgpCUFAQnj9/jgULFqBPnz5yw2phYWGIiorCjRs3pG3BwcEYNmwYpk2bhr59+yI+Ph4rVqzA8OHD5aYGICIiIlKWRjNka4NYLMbmzZsxf/58TJo0CTY2NvDz80NoaKjcfhKJBAUFBXJtHh4eWLNmDb777juMHz8e1atXx4gRI/DJJ5/o8yEQERFROSISCot0yjBjxgwEBQWhXr16mDFjRuknFYnwzTffaCVAfSosyObCulRRZWdnY/To0QCADRs2cL4wIqqQlO45unDhAgICAqTfl6YiLYEhCAJycnLkfgb+ew6sra0r1PNBRERk6pROjk6cOFHi9xWZIAiYO3cu7ty5o3AfV1dXzJkzhwkSERGRiVDrbrWkpCTk5eWVuC0/P1/hBJFERERExk6tgmwfHx/s2rULHh4exbbdunULgwYNws2bNzUOztiJRCLMmTNHOqyWk5ODCRMmAABWrVoFa2trDqsRmajY2FhEREQgJCQErVq1MnQ4RKRHavUclVbDLZFIKlQyIBKJYGNjAxsbG1hbW0vbra2tYWNjU6GeC6LyJCYmBnl5eYiJiSl7ZyIqV7Q6jXRubi7OnDmDGjVqaPO0RER65+vrC0tLS/j6+ho6FCLSM6WH1SIjI7Fy5UoAr3tLhgwZonDfQYMGaR4ZEZEBtWrVCps3bzZ0GERkAEonR+7u7hg2bBgEQcCOHTvQo0cP1KxZU24fKysruLq6ws/PT+uBEhEREemD0smRt7e3dK2zrKws6YSQREREROWJyjVH2dnZ+P333/H333/rIh4iIiIig1L5Vn4bGxvk5OTA1tZWF/EQEZGWcAZ/IvWoNc9R+/btce7cOXh5eWk7HiIi0gLO4E+kPrVu5R8/fjxiYmIQGRmJO3fu4MWLF0hNTZX7IiKi/8TGxiIgIACxsbGGDoWIyqBWz9GAAQMAyN/eX1RFmCGbiEhZspNK6mPGbc7gT6Q+tZKj4OBg/kEREanA19cXCQkJep1UsnAG/6IKZ/AnopKplRxNmjRJ23EQEZVrnFSSyHRovHxITk4Onj59KndHBBEREZGpUqvnCADi4uKwePFi/Pnnn5BIJDAzM0PLli0xdepUtGzZUpsxEhEREemNWsnRlStXEBAQALFYjMGDB8PBwQFPnjzB0aNHERAQgK1bt6J58+bajpWIiIhI59RKjiIiIuDm5oYtW7agUqVK0vZp06Zh5MiRiIiIwPr167UWJBERmb6MnAyF28zNzGFjaVPqvpWtK+skLqKi1O45+uabb+QSIwCoVKkSxowZgy+++EIrwRERUflRZWIVhdt6u/dGTEiM9GeHKQ7IzM2U20dYK+gsNiJZahVkSyQSWFlZlbjN2toaEolEo6CIiIiIDEWtnqNGjRrhxx9/RNeuXYtt27VrFxo1aqRxYEREVL68inylcJu5mbncz0+XPNV1OEQKqZUcjR07FsHBwejXrx/69u0Le3t7JCcn48CBA7h586bCWbOJiKjiUqVmiPVFZEhqJUc+Pj5YtGgRFi1ahIULF0rba9WqhUWLFpXYo0RERERkCtSe56hPnz7w8/NDYmIiUlNTUb16dTg7O3NZESIqF2JjYxEREYGQkBC9rIVW3gmCgKysLI3OYWtry/cY0gu1kyPg9bo9Li4u2oqFiMho6Huh2PJMEAT069cPly5d0ug8bdq0wf79+5kgkc6pnRz9+++/+OGHH3DhwgVpz1G7du0wduxY1KtXT5sxEpkMQRCkS+kIwuvbjgtfyLkCumkxxEKx5Rl/98mUqJUc3bx5EyNHjkRWVhZatmyJpk2bIjk5Gfv27cPBgwexdetWNG7cWNuxEhk1QRAwd+5c3Llzp8Ttrq6umDNnDt8kTAQXitUekUiE/fv3c1iNTIZaydE333wDOzs7bNy4EbVr15a2P3z4EKNHj8Y333yDrVu3ai3I8o61DURU3olEomITBxMZK7WSo6tXr+Lrr7+WS4wAoE6dOpg4cSJmzpypleAqCtY2lA8ikQhz5sxBTk4OcnJyMGHCBADAqlWrYG1tzWE1Uhk/OJEhvHjxAuvWrcOxY8eQlJQES0tLuLm5YejQofD19YWFxX+pw7Vr17BmzRpcunQJr169Qu3atdGjRw8EBgaiatWq0v3+/fdf+Pj4SH+2traGnZ0dmjRpgr59+6JHjx4lvj4WFBTg3XffxfPnz3Hs2DG9le2olRxVrVpV7kHLEovFqFJF8RTxVBxrG8oPkUgEGxsbuTZra+tibUTK4Acn0rcHDx5g5MiRyM3NxahRo+Du7o68vDz88ccfmD9/PiQSCfr37w8AOHz4MKZOnYpmzZph9uzZqFmzJq5du4YffvgBJ06cwNatW2FnZyd3/ilTpqBdu3bIy8tDUlISjh8/jk8++QRdu3bFihUr5BIvADh79iyeP38OADhw4ID0Q6euqZUc+fn5Yc+ePfD29i62bffu3XyTVxFrG4ioJI0aNcKtW7e46gDpTWhoKHJycvDTTz/B0dFR2v7uu+9ixIgRePLkCQDg6dOnmDFjBjw8PLB582ZYWloCANq2bYuOHTvigw8+wLx587B06VK589evXx8tWrSQ/vz+++9j165dmD17NtauXVss+YmOjka1atXg5OSE6OhovSVHaq2t1qRJE1y/fh0ffPABNm3ahAMHDmDTpk344IMPcPPmTbi7u+PIkSPSLyIiUt2tW7fk/iUqy99//42xY8eiXbt2aN68OXr06IG1a9dKt1++fBmjR4+Gp6cnWrZsiUGDBuHs2bMAgEuXLuHKlSsYP368XGJUyMHBAe7u7gCAPXv2ICMjA5MnT5YmRoUaNWqEPn364NChQ3j06FGZMQ8ZMgTu7u7Yvn27XHtWVhaOHTuGHj16oH///khISMDNmzdVfk7UoVbP0bRp0wAAjx49wl9//VXidtnbmPX1YIiIyhMOuZOqJkyYgDfeeANff/01qlSpgvv37+Px48cAXtewBQQEoEWLFpg/fz7EYjH++usvJCUlAQAuXrwIACWOChV18eJFVKtWDW3atClxu4+PD/bt24c//vgDffv2LfN8HTt2xOrVq/Hw4UPUqVMHAHDixAlkZmbC19cXrq6u+Oabb/DLL7/o5W54tZKjLVu2aDsOIiIqgkPuhqWNWb3Vpc60BSkpKXjw4AHCwsKky3i1b99eun3RokWoX78+Nm/eDHPz1wv9vvPOO9LthUNmJfUaFfXkyZNiN2XJKkxwChOzshRe89mzZ9Jjo6OjUatWLbRt2xZmZmbo2LEjYmJi8Nlnn8HMTK2BL6WplRy1bdtW23EQEREZDW3N6q0udWYDr1GjBurUqYMlS5YgLS0NXl5eePPNNwG8HqL6888/MWXKFGliVFThiI82KRt/0UlzX7x4gd9//x0jRoyQJkJ9+vTB6dOn8ccff6Bdu3Zaj1WWRqnXq1ev8Pvvv+PAgQM4e/YsXr16pa24iIiIDMrUpt4QiURYt24dnJ2d8dVXX8Hb2xsDBgzAH3/8gfT0dEgkEjg4OCg8vjCRUqZOqFatWtLhuJI8fPhQup8yCnuYatasCQA4dOgQ8vLy4O3tjfT0dKSnp6Nt27awtrZGdHS0UufUhNrLh6xfvx6RkZHIzs6WZny2trYICQnBRx99pLUAiYiI9E1bs3qrS93ZwJ2dnREREYG8vDxcvnwZS5Yswfjx43Hq1CmYmZnh6dOnCo8t7I05c+YM/P39S71O27Ztcf78eVy6dAmtW7cutv3EiRMwMzNTWJNU1O+//45atWpJh+oKE6BRo0YV2/fw4cOYPXs2rKyslDq3OtRKjqKiorBo0SJ06tQJ/fv3h4ODA54+fYqoqCgsXLgQNWrUQL9+/bQcKhERkf6Y8qzelpaWaNu2LQIDAzFhwgQ8e/YMLVq0wM8//4zRo0eXOLTWqlUrtGjRAmvWrEH37t2L9fo8e/YMjx49gru7OwYNGoT169dj2bJl2LRpk9z8RHfu3MEvv/yCnj17KlW/tGvXLvz111+YOnUqACApKQlxcXEYOnQoevfuLbfvnTt3MH/+fJw5cwbvvfeeOk+NUtRKjjZt2gQ/Pz8sXrxYrr1Xr14IDQ3F5s2bmRwREZFWcKZw5dy6dQvffvstevfujXr16uHVq1dYs2YN6tSpAycnJ0ydOhWjRo3CqFGjMGzYMFSrVg3Xr19HjRo18MEHHwAAFi9ejJEjR2LgwIEYNWoUmjVrhvz8fMTGxmLHjh2YPn063N3d4eDggPDwcEydOhX+/v7w9/dHzZo1cf36daxevRpOTk6YNWtWsRjv3buHK1euID8/H0lJSTh27BgOHz6Mbt26YcyYMQBeT/YoCAI+/vjjYjNit27dGmvWrMEvv/xifMlRYmIipkyZUuK2vn37YuLEiRoFRURkaLt370ZUVBT69euHwYMHGzqcCo0zhSvH3t4eNWvWxJo1a/DkyRNUrVoVrVu3xqJFi2Bubo7WrVtjy5YtWLZsGWbMmAEzMzO8/fbb+PTTT6XnqFevHn766SesW7cOe/bswfLly2FhYYFGjRrh008/hZ+fn3TfHj16oHbt2lizZg3mzp2LjIwMODo6YtCgQRg3blyJK2ksWbIEAGBlZQU7Ozs0bdoUy5cvl1s+JDo6Gq1atSpxqRBzc3P06dMH27dvx6tXr3S2IodayZGNjQ3S0tJK3JaWlsalEojI5B0+fFj6L5Mjw+J8T8p54403sGjRolL38fT0LHM6Hjs7O0ybNk06p2Fp3N3dERkZWeZ+devWxe3bt8vcD0CZBdeff/45Pv/8c6XOpS617lZr1aoVIiMjpXMiFEpOTsbKlStLLM4qzd27dzFmzBi0aNECXl5emD9/PrKzs1U6x9GjR+Hm5iaX1RIRqatHjx5y/5LhFM73xF4j0he1eo4mT56MoUOHonv37vDy8oK9vT2Sk5Nx/vx5WFhYKJVFFkpPT0dAQABq166NiIgIpKSkIDw8HKmpqcVqmhTJzs5GeHi49BZAIiJNDR48mD1GRBWUWsmRq6srfvrpJ6xYsQIXLlxAamoqqlevDh8fH0ycOBENGjRQ+lw7d+5Eeno6oqKipKv3mpubIzQ0FBMmTICLi0uZ51izZg1q166NunXrlricCREREZGyVE6OcnJysHLlSnTv3l1aWKWJM2fOwMvLS5oYAa+7scPCwnD69Okyk6P79+9j48aN2LlzJzZt2qRxPERERFSxqVxzZG1tjU2bNmltYqyEhIRiCZCVlRWcnJyQkJBQ5vFff/013n//fTRq1Egr8RAREVHFptawmouLC/7991+lZ74sTXp6OsRicbF2sVis8I64QidOnMDly5dx6NAhpa/n4+OjcNujR4+UmrCKiIiIyi+17lYLCgrCqlWrcP/+fW3HIyUIQqlTp+fk5OCbb77BpEmT5IbkiIjKi9jYWAQEBCA2NtbQoRBVKGr1HP3000/IyspC79694erqWmwhO5FIhFWrVil1LrFYjPT09GLtL1++LLXeaPPmzTAzM4Ovr6/0+Ly8PEgkEqSnp8PGxqbEdVeOHz+u8Jyl9SoREekbJz8kMgy1kqM7d+7A0tISDg4OSE1NRWpqqtx2VRbLc3FxKVZblJubi/v372PgwIEKj0tMTMS9e/fg5eVVbFubNm3w5Zdf4sMPP1Q6DiIiY8PJD4kMQ63k6MSJE1oLoFOnTli1ahVevHiBGjVqAHg9oWNubi68vb0VHjd27Fj0799fru2HH37A3bt3ER4ejrfeektrMRIREVHFoVbNEQAUFBTgwIEDmDVrFiZNmoTZs2fjwIEDyM/PV+k8Q4cORdWqVREUFITffvsNUVFRmDdvHvr06SM3rBYWFoYmTZpIf3ZxcUG7du3kvuzt7VGpUiW0a9eu2GrCRESmRnZYjags+/fvh5ubG/755x+59u3bt8PNzQ1Lly6Va8/IyECTJk2wcuVK+Pv7Y9y4cdJt+/btg5ubG1JSUvQRutFRq+coJSUFH3/8MW7cuAELCwtUr14dqamp2L17NzZs2IB169YpXSQtFouxefNmzJ8/H5MmTYKNjQ38/PwQGhoqt59EIkFBQYE64RIRmSQOq5EqCuvS4uLi5EZPLl++DFtbW8TFxcnt/+eff6KgoACenp7o0aMHzMzU7i8pd9RKjsLDw3H37l0sXrwYvXr1grm5OQoKCnDw4EHMmTMH4eHhZS5+J6tBgwZYv359qfssWLAACxYsKHMfIqLyonBNMSJlODk5wd7eHnFxcRgwYIC0PS4uDv3798e+ffuQl5cHS0tLabuFhQWaN2+OSpUqGSpso6RWmnjy5El8+umn8PPzg7m5OYDXS3706dMHISEhOHnypFaDJCIiorK1bNlSrofoyZMnePjwIfz9/VFQUICbN29Kt8XFxaFx48aoVKlSsWG1kuTm5mLp0qXw8fFBs2bN0KlTJ8yYMUNun2PHjqFfv35wd3dHx44dMXfuXGRkZAAAkpKS0KpVK3z77bdyx4wbNw4+Pj5y+4WEhKB169Zo3rw5Ro4ciWvXrskd07VrV3z11VfYtm0bunTpglatWiEoKEhrw4BqJUeCIODtt98ucZurqysEQdAoKCIiIlJdq1atkJiYKL2LPC4uDg4ODnB2dkbTpk2liZNEIsGVK1fg6emp9LknTZqETZs2YeDAgfjhhx8wbdo0aUIDvJ4qZ+LEiXjrrbcQGRmJoKAg/PzzzwgODgYA1K5dGzNnzsSmTZtw8eJFAMCuXbtw5swZfPvtt6hcuTJevXqFESNG4K+//sLs2bPx3XffITc3FyNHjix2Z/uJEydw8uRJzJ49G1988QUuXryIefPmafL0Sak1rNahQwf873//Q4cOHYptO3v2LNq2batxYERERIaWkZOhcJu5mTlsLG2U2tdMZAZbK1u19lWFp6cnBEFAXFwcunbtisuXL6Nly5YAXvcqxcbGYtSoUbhz5w4yMjKUnj/r7NmzOHXqFL777jv4+flJ22W/j4yMhLu7O5YtWyZtq1atGqZOnYoLFy6gXbt26N+/P44fP47p06cjMjISCxYswJgxY9C6dWsArwvBk5KSEB0dLe2E8fLyQteuXbF27Vq58hlBELBq1SrpnIb37t3D+vXrIZFINK6fUis5CgoKwsSJE1FQUIA+ffqgZs2aePbsGaKjo3HkyBFERkbKzX1UvXp1jYIkIiIyhCoTqyjc1tu9N2JC/ruT0GGKAzJzM0vc19vVG6c+OyX9+a3pb+HZq2cl7tu6fmv8MfMPteJt3LixtPi6a9euiIuLkyYwLVu2RHR0NABIZ11Xtufo3LlzsLW1VXhzQEZGBm7evIlp06bJtffs2ROff/45Ll26hHbt2gEAvvrqK/Tp0wdDhgyBs7MzQkJCpPtfunQJb7/9ttzoVOXKldGlSxdcunRJ7txt2rSRm+y5YcOGyMvLw/Pnz2Fvb6/U41JEreSocH6hjRs3YtOmTdL2wuE02UIwAHJjnEREpBuCICAnJ6fEbbLtivaxtrZWaRJfMj6WlpZwd3dHXFwcsrKycOvWLcyePRvA6+To2bNnuH//PuLi4qQF3MpITU2Fvb29wt+Ply9fQhAE1KxZU6698I522bVS7ezs0KFDB/zyyy8YPHiwXIKTnp5e7BwAULNmzWLrrRZdl7Ww0FzR77cq1EqOgoOD+QdERGRkcnJyMHr06DL3mzBhQontGzZsgI2NTYnbKqpXka8UbjM3M5f7+emSpwr3NRPJD/P8s+AfpfdVlaenJzZu3IjY2FiYm5ujcePGAAAHBwfUqVMHcXFxuHz5skolMNWrV0dycrLCdU+rVq0KkUiE58+fy7Xn5+cjNTUV1apVk7b9/vvviI6ORpMmTbBixQr06NFDmhBVq1YNiYmJxc7/7NkzuXPomlrJ0aRJk7QdBxERkdGpbF3Z4PuqytPTE6tXr8bWrVvh7u4u7VEBXvceHTp0CA8fPlRpvb4OHTpg7dq1+PXXX9G7d+9i2ytXrozGjRvj119/xUcffSRtP3LkCPLz86U1RWlpaQgLC4Ovry/mzp2Lvn37YubMmVi9ejWA1wXlhw8fRnx8PBo2bAgAyMzMxMmTJ9GlSxe1ng91qJUcERGRcXP4oAtEFvI9G4WlD7Kf/IX8Ajzdy+lXyhNPT0+YmZnh9OnTGDt2bLFthXd0qXKnWocOHeDt7Y2wsDDcv38fzZs3R2pqKo4cOSKdeXvixIkIDg7GlClT0K9fPzx48ABLliyBl5eXtN5o7ty5EAQBs2fPRpUqVbBgwQIEBARgz549GDRoEAYMGIBNmzZh3Lhx+PTTT1GpUiWsW7cOOTk5xR6LLjE5IiIqh0QW5jCzLPslXqKHWEi/qlatioYNG+LOnTvSO9UKtWzZEoIgoHr16nB2dlbpvCtWrEBkZCR27dqFyMhIvPHGG+jYsaN0u4+PD1asWIGVK1ciKCgIYrEYffv2la54cfDgQcTExGDt2rXSIbK2bdti1KhRCA8Ph5eXF+rWrYtt27ZhwYIF+PLLL5Gfnw8PDw9s2bJFbkkxXRMJnJRIysfHB8DruRrUkZ2dLR3v59h9xWaqvwumGje9Jvv/V2voe8olR3n5eLLzGAD+nxMV4kIqRERERDKYHBERERHJYHJEREREJIMF2UQGJjtxX9G7iTgpHxGR/jE5IoMpLSkAKkZiIAgC5s6dizt37pS43dXVFXPmzCn3zwMRkTFhckQAXq+zExERgZCQEJUmBlNXWUkBwMSAiIgMg8kRAQBiYmKQl5eHmJgYvSRH9JpIJMKcOXOQk5ODnJwc6bIOq1atgrW1dYXoPSMiMjZMjggA4Ovri4SEBIUrLmtbWUkBUDGG1YDXz0XRuWWsra053wwRkYEwOSIAr9ez2bx5s16vyaSAdImF7kSkLiZHRFTusNCdiDTBeY6IiIiIZLDniIjKHRa6E5EmmBwRVUCy9TiyZNtK2l7IFJILU61pY60UkeExOSKqgHJycqSrtytS2NtSEq7erhuslSIyDkyOVKToEzeg3KdufvIjMl7stSEigMmRypT5xA0o/tTNT9xkbOa1c4SV+X9v+iUt5QIAuQUCZl14pNfY9EnXvTbKzELPWiki48C71YgqOCtzEazNzaRfNhbmsLEwl2uzNjeTS6BIdbKz0JemsFaqcDJU4L9aKSZGRPrBniMNOHzQBSILc7m2kj51C/kFeLr3pF5jI/2oCIXNFYWue230PQs9EamPyZEGRBbmMLMs+ymU6CEWMgwWNpcvurzDzRCz0BORejisRkRGJzY2FgEBAYiNjTV0KERUAbHniEhLig6zKips5jBr2WTrcxQVLxMR6QqTIyIt4TCr9rA+h4gMickRERkd1ucQkSGx5sgIsL6CiIjIeDA5MgLKzn9CREREusfkyAj4+vrC0tKS9RVERERGgDVHRoD1FURERMaDyRERGQUu+kpExoLJEREZnK4XfSUiUgVrjoiIiIhksOeIiAxO14u+EhGpgskRERkFXS76SkSkCiZHRKR1LK4mIlNmFMnR3bt3MX/+fMTGxsLW1ha+vr4IDQ0t9RPjq1evsHHjRpw5cwZ3796FhYUFmjZtiilTpqBp06Z6jJ6IZLG4mohMncELstPT0xEQEICMjAxERETg888/R3R0NGbOnFnqcUlJSdi1axc6dOiApUuXIjw8HBKJBEOHDsX169f1FD0RERGVNwbvOdq5cyfS09MRFRUFOzs7AIC5uTlCQ0MxYcIEuLi4lHhc3bp1cfToUdja2krbOnToAB8fH2zbtg3h4eF6iZ+I5LG4mohMncF7js6cOQMvLy9pYgQAPXr0gJWVFU6fPq3wuEqVKsklRsDrWgYXFxc8ffpUZ/ESUdkKi6utra2lbYXF1UyMiMjYGTw5SkhIKNY7ZGVlBScnJyQkJKh0rszMTNy8eRPOzs7aDJGIiIgqEIMPq6Wnp0MsFhdrF4vFSEtLU+lcy5YtQ1ZWFkaMGKFwHx8fH4XbHj16BEdHR5WuSVRIyC+ARMn9iIjIeBk8OVJEEASVut+jo6OxefNmzJ49G/Xr19dhZET/KbxNHQCe7j2p8vHZ2dnS7wtvfS/6fVGs2SEi0i2DJ0disRjp6enF2l++fKmwGLuos2fPYsaMGRgzZgyGDx9e6r7Hjx9XuK20XiUiXQgKCiqxvbCIuSQbNmzgxIhERDpk8OTIxcWlWG1Rbm4u7t+/j4EDB5Z5/NWrVzFx4kT07NkTn332ma7CJCqRbA+OwwddILIwL/MYSXYukqPO6DKsCkd20smilOmRY28cEckyeHLUqVMnrFq1Ci9evECNGjUAAEePHkVubi68vb1LPTYhIQFjx46Fp6cnwsPD+eJGBiWyMIeZZdl/UrI1R/PaOcLK/L/f26KzSRfKLRAw68IjLUVa/uTk5GD06NFl7qeoR469cYrFxsYiIiICISEhaNWqlV6uWTTZ5SzrpG8GT46GDh2Kbdu2ISgoCEFBQXj+/DkWLFiAPn36yA2rhYWFISoqCjdu3AAAPH/+HGPGjIGlpSU+/vhjuYkfrays0KRJE70/FiJVWZmLYG2uzE2jypR6E2lfTEwM8vLyEBMTo5fkqKwZ1gHOsk66Z/DkSCwWY/PmzZg/fz4mTZoEGxsb+Pn5ITQ0VG4/iUSCgoL/PnHHx8fj0aPXn6RHjRolt2+dOnVw4sQJncdORManpOHNknrkJHn5SP7pFIDiw20civuPr68vEhIS4Ovra+hQiPTG4MkRADRo0ADr168vdZ8FCxZgwYIF0p/btWuH27dv6zo0IjIx6gxvllYAz6E4/ZKdYR0AZ1kngzD4JJBERKR9Qn4BJHn5ZX6VNe+W7LCavhTOsM5Z1slQjKLniIjIkIoWxgMlD8UZe2G8pvNuCYJQrBi6W7duiI+PR/fu3ZGdnc1eG6oQmBwRlWOGuNPIFLEw/rXSiqEjIiIAsBiaKgYmR0TlmL7vNCLDUmfeLSG/QNrLxISH6DUmR0TlGO80qriULUyX7QtjMTTRa0yOSG8UzWLMNcV0p1WrVti8ebNOr8HZqcuXwmLoogqLoYkqAiZHpDfKzGLMNcVMD2enJqLyhskRUQWXW6BckbGy+xERmTomR2QQXFPMsGRv+Z514bFGx8tSdnZq2SJgzk5NRMaGyRFpnTK3j/PW6fJJnSJgzk5NRMaGyRFpHW8fN36yvS7z2r0JKyUS1dwCibSXib02RFSeMTkirePt46bFytxMyV483VFmmJVDrKaPk5KSqWByRFqnj9vHqXxRbpi1/A2xVrRpENirTKaCyZEGhPwCpV6uy1rYkYgqpoo2DQJ7lclUMDlSkTYWdiTSJdneiKLDU6bW00Dli6Je5YrWg0bGj8kRUTlS2sKhwH+LhpI8Y5jrqWjdFVBxaq8qWg8aGT8mRyriwo5E5YOu5npSF6e3IDIeTI40oM6cLrrGu0EqNtmFQ01p0VBlemQ4Q3fFUJF70Mh4MDkqZ3g3iPrKS2JZ0sKhxr5oqKo9N9roteFcT8aJPWhkDAw7uQlpna+vLywtLXk3iBpkE0uqWArneirrS5kEiohMH3uOyhnOMaQ+3mZsOMr03LDXhoj0hckR0f9jYmk4xjBLNxFRISZHZBDGcOs0ERFRSZgckd4Y263TREREJWE/NhEREZEM9hwRAP3cxs5bp8s/rjdoemR7ZIsuz8GlO6iiYnJkBIxhfh19z49kjAW4xvD/YIo0XW+wvCq6XlhJ69wZBZm6PkXLc5S2jUt3UHnE5MgIGMPEjbyN3Tj+Hyqq3AIBspP6KZoR2VSUtcYd8Hqdu88//1yPUZkG3qxBxoDJkREwhsSEt7Ebx/+DKdJ0vUEAWl0GwlBDe7I9RYIgQCIpPQqJRKJwqMpQDLV0B2/WIGPD5MgIMDExHNk3tCZNmmD16tUQiUTIzs4GYBr1FEWHbwqVVi+iqzdlQ603qOnQnjbeXJVdWb5QfHw8goKCNL6uNnHpDqLXmByRVhT91AzIf9I0mvoKGcoOfcyZM8eoEyRl3pRLqyUxtMJFcQuVtGBuUcb4+0Tq480aZGyYHJHGWF9BmihtUVxlF8xVZ2hPkp2L5KgzAIA8CSAqUsOiSd1TSTGUdD7ZGOg1Y7xZgyoeJkdGgHdJGYZIJMKcOXOQk5OjsLfCFIbVZBV9Uy7pDblovY82lVTvoygGXVF2aE82Bm3Wz6gbAxEZDyZHRsDU75JSNsnQZ/GpsgmnSCQq1jOhbG+FMVLmTVmX1SK8lZ+IygMmRwZQtIC2W7duiI+PR/fu3ZGdnW1yvRWA8SUZpp5wkh7IDN2w7uk1XdxGb4r1iERMjvSstPqciIgIAKZRBGzseFu+/lhbW2PDhg0lbjPmJKPoG7SmdU+mSpe30atTj1h0zivZa5jqvFdkepgcUbnE6RH0p6Rew5KU9yRDU5z88DVt138RqYPJkZ7J1ucAJX+yNoVhNWOaW4c0o8zs1P/tR9pkDJMf6vI2emOsRyRSBpMjA1D0SduQn6xVvWPOVOfWUSepk2UKiauq+Eld+0x1AV5d3EavTD2iqQ7NUvnF5EgDxnDbsrbs3LkTeXl52LVrV7kuYNY0qTP2RTaVeVM2hd9HU6TpLN1fta0FayXmZ1K218ZQr0+6/ADCoVnSFyZHGiiPty2r002vzNw6ACe80xVN3pS///57uTcbZT6lA/ykrgva7rUx1OtTef8AQhUDkyM9K60OR9GnKV0M5RT9dDdgwACsWrUKAwcOVHk6AVOd8K7oIpul1dqU16GninyXlrZpOkt3eRuuJTJlTI5UpM7YuGx7aZ+YFG3T9icpY5hOwBiKgMvLIpuqvinLzpBtTG/IhQl7SR8YTK3Wy1AfGIq+PgmCgNzcXACvn8tPP/0UALBs2TJYW1vDysoKubm5OqsPtO/XSaleZSG/gL3KZFSMIjm6e/cu5s+fj9jYWNja2sLX1xehoaFKJQT79+/HmjVr8PDhQ9SvXx/BwcHo1auXzmLlbcvaUR57YorWeJT2RlBI27dvG3qGbHUpStgL37RVSdbLUy2gqmRfn0r7EFSYJOlijiHZYV51Eh5t3IFHpCmDJ0fp6ekICAhA7dq1ERERgZSUFISHhyM1NRWLFy8u9dhDhw5h+vTpCAwMRMeOHXHs2DFMnjwZVatWxTvvvKOnR6CaWa1ryQ3lAIpffOZdelLiOYoOicl+OiyJlZXV63P+/z5WVlb4/PPPy/xEKXuN0j65G+rOHFWeBysrK4hEIrn9tZ2YqFPjYajbt8uz8lgLqC/G8qFF9m+7PPQkkukxeHK0c+dOpKenIyoqCnZ2dgAAc3NzhIaGYsKECXBxcVF47PLly9GzZ09MnToVANC+fXvcvXsXERERRpUcyb6ZKUp4lD0eUK7gUV2FSVJRRYf2NL0zR7YQWN0iYE2fByYmuqPqEFnR+b+KfmDQ15uhIYf2tD0zdNHntKTz6WKOIdlY1R1WU9TjpU5PIpE6DJ4cnTlzBl5eXtLECAB69OiBsLAwnD59WmFy9ODBAyQmJmLKlCly7X5+fpgxYwZSUlLkzqkrynzCoeIUDTuqMhxpDImKohq00hI+2W1FexJLq70qTKx18YagzaRA3SEyZYesS6KNeXIMPbSni14bZZ5TXc4xxDoiMlUGT44SEhIwcOBAuTYrKys4OTkhISFB4XGJiYkAAGdnZ7l2FxcXCIKAxMREnSdHpY3py76gyo7pq6PYpywDJAVFr6npG7RIJCr1DRkwzq7zst7QFc3xApT8eNTpSSyNkF+AAkEAFA0FmpsV26bNpMBQSkoCVP2/0ObflSkN7anyOwzo/u9S2R4vY/59JNNn8OQoPT0dYrG4WLtYLEZaWprC4wq3FT22WrVqctuL8vHxUXjOR48ewdHRscyYSf7TpiAICA8PR3x8vML9GzZsiBkzZkhf0KysrEp9QwbKflM2thdHZZMMXTKGN2VjGCJTJ+EzdNyrVq2ClZVVqX9LRf+OtN0zrewHvtL+LjV9bSh8nnlDCxmSwZMjRQRBUOqFSNGnP328cSr7CQeAxi8WRa+rb6X1mAiCADOz0m+JNzMzg42NjfQ8xjAk9v3332PJkiUa/V+oQ5U3j6LX11YMsnSRFJjqm5u2hvY0+Tsv7W+p6N+RMdL0tYHIGIgEA79LeXl5YeDAgQgNDZVr9/X1RYsWLfD111+XeNzp06cRGBiIgwcPytUlXb16FYMGDcL27dvRunVrlWIp7FU6fvy4io9CNSXdZQUo/2akjbvVZH9WdEzRN2RlYyopQVU0jKHqMYqOL/y5rLvVZMnOQ6VuDKXFpOz/qzrHKDq+8Gdln4fyPDyh6fOqzeurEoOpxq3M+bTxN0akDwbvOXJxcSlWW5Sbm4v79+8Xq0WSVVhrlJiYKJccJSQkQCQSFatFMiaafqou6XhbW9syjyu6jzLHaBKTLo4p63h1HpM2ezhM+Xkobwzde6Xu9U01bn2dj0gftLv8sho6deqE8+fP48WLF9K2o0ePIjc3F97e3gqPq1evHpydnXHw4EG59gMHDsDDw0Mvd6oRERFR+WPw5Gjo0KGoWrUqgoKC8NtvvyEqKgrz5s1Dnz595HqEwsLC0KRJE7ljQ0JC8Ouvv2Lp0qW4cOECvvnmG5w9exYhISH6fhhERERUThh8WE0sFmPz5s2YP38+Jk2aBBsbG/j5+RWrQZJIJCgokJ9huVevXsjOzsbq1auxfv161K9fH0uXLjWqCSCJiIjItBi8INuY6Ksgm4iIiIyXwYfViIiIiIwJkyMiIiIiGUyOiIiIiGQwOSIiIiKSweSIiIiISAaTIyIiIiIZvJVfhru7OwoKCuDo6GjoUIiISIscHR2xbds2Q4dBJoI9RzKsra1hYaGbeTEfPXqER48e6eTcjIExmGIMhr4+Y2AMRIqw50hPjGGCScbAGIwpBkNfnzEwBiJF2HNEREREJIPJEREREZEMJkdEREREMpgcEREREclgckREREQkg8kRERERkQzeyk9EREQkgz1HRERERDKYHBERERHJYHJEREREJIPJEREREZEM3ayyWoGdPn0a27Ztw7Vr1/Dy5UtUq1YNHh4eGDp0KLy9vSESiXDt2jXs2LEDV65cwd27d+Ht7Y01a9boNYadO3fi6NGjuH37NjIzM9GgQQOMGTMGvXv31lsMP//8M7Zv345//vkHWVlZqF27Nvr27YuxY8fCyspKLzHIevz4MXr16oXMzEycO3cOdnZ2eolhxYoViIyMLHbsl19+iQ8//FAvMQBAQUEBNm7ciD179uDhw4eoXr06unbtiq+++kovMXTt2hUPHz4s8fhdu3ahRYsWOr0+AOzZswdbt27F/fv3IRaL0aFDB0yePBm1atVS+9qqxvDTTz9h7dq1+Pfff+Ho6IiRI0fC399fJ9dT5bVo//79WLNmDR4+fIj69esjODgYvXr10lsMBw8exK+//oorV67g6dOnmDZtGsaMGaPy80KkDCZHWrRkyRKsWbMG3bp1w+zZs2Fvb49nz57h2LFjmDBhAlatWoXOnTsjLi4Oly5dgoeHB3JycgwSw6pVq9CxY0cMHToUlStXxokTJzB58mS8ePECw4cP10sMaWlp6NSpE8aPH49KlSrh6tWriIyMxOPHjzFv3jy9xCBrwYIFqFSpEjIzMzW6tjox2NjYYPPmzXLH16tXT68xfPHFF/j9998RFBSEhg0b4tmzZ7h27ZreYoiMjERubq7csYsXL0ZCQgKaNWum8+vv3bsXM2fOxEcffYQZM2bg8ePHWLp0KW7duoV9+/bBzEz9jnZlYzh48CDCwsLg7++Pzp07IzY2FuHh4RCJRBgxYoTWr6fsa9GhQ4cwffp0BAYGomPHjjh27BgmT56MqlWr4p133tFbDA8ePECXLl2wa9cupZ8LIrUIpBUnT54UXF1dhYiIiBK3//nnn8K1a9cEQRCEgoICafuIESOEwMBAvcfw/PnzYtuDg4OFLl266C2GkixZskTw8PAQ8vPz9RrD//73P6Ft27bC+vXrBVdX1xKfH13FEBERIbRo0UKj62kaw2+//SY0adJE+Pvvvw0WQ1EZGRlCixYthC+//FIv1w8ICBCGDx8utz0mJkZwdXUV4uPj9RJDjx49hIkTJ8ptnzt3rtC2bVshNzdX69dT9rWoZ8+eQkhIiFzb6NGjhUGDBuktBtn9XF1dhXXr1pW4H5E2sOdISzZu3Ah7e3tMmDChxO0eHh7S7zX5BKqtGEoaMmrcuDFOnTqltxhKUr16deTn50MikcDc3FwvMeTl5WHevHmYNGkSKlWqpNY1NY1BF1SJYe/evWjfvj0aNmxosBiKOn78ODIzM9GnTx+9XD8/Px9VqlSR2y4WiwEAggbTwSkbQ1ZWFv755x989NFHctvfffddbN++HVeuXEGbNm20dj1AudeiBw8eIDExEVOmTJFr9/Pzw4wZM5CSklLs9UQXr4e6et0kKgl/27QgPz8fcXFxaN++PSwsDJNvaiOG2NhYuLi46D2G/Px8ZGVl4dKlS9i8eTM+/PBDWFpa6i2GLVu2wNzcXCv1PerGkJ2djfbt26NJkybo3bs3du/erdcY/vzzT7z11luYP38+WrVqBQ8PD4wbNw4PHjzQWwxFHThwAHXq1EHLli31cv3Bgwfj999/x8GDB/Hq1SskJCRgxYoVGiWNqsSQm5sLQRCK/e4X1t8lJCRo9XrKSkxMBAA4OzvLtbu4uEAQBOl2XcZApG/8zdWC1NRU5ObmwtHRUa5dEAQUFBRIfzYzM9PZpx9NYzh27BjOnj2LhQsX6jWG/Px8NG3aVLqtf//+CAsL01sMT548wcqVK7Fy5Uq1e6o0jcHJyQmhoaFo0qQJcnJyEB0djVmzZuHly5dqF5yqGkNycjL27dsHV1dXLFmyBJmZmViyZAkCAwMRHR2t1pucJr+TL168wNmzZzF69OhihfO6un7fvn2RnZ2NadOmIS8vDwDQqlWrEovldRFDtWrVUKNGDVy9ehUDBgyQbrty5QoAIC0tTavXU/a1qPC6hb1ohapVq1ZiXMbwekikKf5makFhl3vRF/HDhw+jadOm0q/58+cbZQzx8fGYMWMGevbsiffff1+vMVhYWGDv3r3Yvn07ZsyYgZMnT2LGjBl6i2HhwoXo2LEjvLy81L6mpjG8//77GDNmDLy8vNC5c2d899136NGjB1atWiV9k9Z1DIXHrFq1Ct7e3ujVqxeWLVuGxMREHDlyRG8xFPr111+Rl5cHPz8/ta6tzvWPHDmC8PBwjB8/Hlu3bsXixYvx/PlzBAcHy72p6zKG4cOHY9++fYiOjkZaWhpOnTqFLVu2lHgObVxPFUXPqehaxvB6SKQp9hxpQY0aNWBlZYXHjx/LtXt5eWHv3r0AoHDs3dAxPH78GB9//DFcXV016jXSJAZ3d3cAQOvWreHo6IiQkBCMGDFC2q6rGC5fvozDhw9j9+7dSE9PB/C67gMAMjIyYGtrC1tbW53GoEivXr1w+PBh3L9/X62hTlVjEIvFqFu3rlztSNOmTVG1alXEx8erfH11YpB14MABuLm5wdXVVa1rq3p9QRAwZ84cDB48GBMnTpTu+/bbb+P999/H8ePH0b17d53GAACBgYG4f/8+PvvsMwiCgEqVKiE0NBRfffUV7O3ttX49Zcj2ENWsWVPaXvg3U7RHyRheD4k0xeRICywsLODp6Ylz586hoKBAOjxTrVo16Ru8Nubt0XYML168wOjRo1GlShV8//33sLa21nsMRRUOsd2/f1+t5EiVGO7evYu8vDz079+/2Hnee+899O7dG0uXLtVpDIpoUgCsTgwuLi4Kb6NWd+hD3echKSkJcXFxxQqAdXn9lJQUpKSkoHHjxnLncHV1hbm5Oe7fv6/zGADA2toaixYtwowZM5CcnIx69epJk9PmzZtr/XrKKKw1SkxMlEvUExISIBKJitUiGcPrIZGmOKymJR999BGePn2K1atXm0QMGRkZGDt2LDIzM7Fu3Trpp0N9xlCS2NhYAJrN8aNsDO+++y62bNki9zV27FgAwMqVKxEcHKzzGBT59ddfIRaL4eTkpJcYOnfujNu3byMlJUXaVjhxX6NGjfQSQ6EDBw4AgEZDaqpe387ODra2tvjrr7/k2m/cuIGCggLUqVNH5zEUjcfNzQ2VKlXC9u3b0bp162JJiDavV5p69erB2dkZBw8elGs/cOAAPDw8Srzz1RheD4k0wZ4jLencuTMCAwMRERGBW7duoVevXnBwcMDLly9x6dIlJCcno3LlygBef0q9ePGi9PuMjAwcOnQIAODt7a3WUI6qMYSEhODmzZuYN28eHj9+LNcF3qRJE7U/2akSw/Dhw9GtWzc4OzvDzMwMV65cwYYNG/Duu+9qdKu7sjHY29sXG6oonKHZ09NToxmyVXkeBgwYgP79+6NBgwbIzs5GdHQ0jhw5grCwMLXv2lM1hiFDhmDbtm0YN24cxo8fj+zsbCxZsgRNmzZF165d9RJDoQMHDsDT0xO1a9dW+7qqXl8kEuHDDz/Eli1bUKVKFbRr1w5PnjzBihUrUK9ePXh7e+s8BuD1jNL3799Hw4YNkZaWhujoaFy4cAE//vijTq6n7GtRSEgIJk+eDCcnJ3To0AHHjx/H2bNnsW7dOr3FEB8fLzfEe+fOHRw6dAi2trYa/f8QlUQkaNp/T3JOnTqF7du34+rVq3j16hWqVauGpk2bom/fvvD19YWZmRkuXLiAkSNHlnj88ePHUbduXZ3H4ObmpvB4fcUQHh6O3377DUlJSbCwsEDdunXRr18/DBs2TCvd7srEUNS+ffswY8YMrS0fokwMn376Ka5evYpnz55BJBLB1dUV/v7+6Nu3r8bXVzYG4PVQ5tdff42LFy/CzMwMnTp1wowZM+Dg4KC3GOLj4+Hr64s5c+Zg2LBhGl9Xlevn5uZi06ZNiIqKki6f0rp1a0yePFnjvwdlYyi8Y/TevXuwsLBA27ZtMXXqVLXqzrT9WrR//36sXr1aunzIxIkTy1w+RJsxKFpmp06dOjhx4oSyTwuRUpgcEREREclgzRERERGRDCZHRERERDKYHBERERHJYHJEREREJIPJEREREZEMJkdEREREMpgcEREREclgckREREQkg8kRkYm5cOEC3NzccOHCBZ2cPz4+HitWrMC///6rk/MTERk7JkdEJCc+Ph6RkZHSdeaIiCoaJkdEpBdZWVmGDoGISClMjoiMUEJCAqZMmYIOHTqgWbNm6Ny5M6ZNm4bc3NwS9/f394e/v3+x9unTp6Nr165ybTt27EDfvn3RsmVLtGzZEj179sSSJUsAvF5495NPPgEAjBw5Em5ubnBzc8O+ffukx//vf/9DQEAAPD090bx5cwwdOhTnzp2Tu8aKFSvg5uaG69evIyQkBG3atEG3bt00ek6IiPTFwtABEJG8W7du4cMPP0SNGjUQEhKC+vXrIzk5GSdOnFCYHCkrJiYGc+fOhb+/Pz7//HOYmZnh3r17iI+PBwB07twZU6ZMwZIlSzB79mw0bdoUAODk5AQA+Pnnn/H555/Dx8cH3377LSwsLLBr1y6MGTMG69evh5eXl9z1Jk2ahN69e2Po0KHIzMzUKHYiIn1hckRkZMLDw2FhYYG9e/fCzs5O2t63b1+Nzx0XFwexWIyZM2dK22QTGjs7O9SvXx8A0LBhQ7Ro0UK6LSsrC9988w06d+6MlStXStu9vb3Rv39/LFmyBHv27JG7Xr9+/RASEqJx3ERE+sRhNSIjkpWVhT/++AO9evWSS4y0xd3dHenp6ZgyZQqOHTuGlJQUpY+9fPkyUlNT0b9/f+Tn50u/JBIJ3n33XVy7dq1Y71D37t21/RCIiHSOPUdERiQ9PR0FBQWoVauWTs7fr18/FBQUYM+ePQgJCYFEIoG7uzs+/fRTdOzYsdRjnz17BgCl9gSlpaWhUqVK0p8dHBy0EzgRkR4xOSIyItWqVYO5uTmePHmi0nFWVlZ49epVsfYXL14Uaxs4cCAGDhyIzMxM/PHHH1ixYgXGjRuHw4cPo06dOgqvUaNGDQDArFmz0Lx58xL3eeONN1SKm4jIGHFYjciI2NjYoE2bNjh06JBKQ1516tTBP//8I1ew/eLFC1y+fFnhMZUqVYK3tzfGjx+PvLw8aVG2lZUVACA7O1tuf09PT4jFYsTHx8Pd3b3Er8JjiYhMGXuOiIzMjBkz8OGHH2Lw4MEIDAyEk5MTnj9/jhMnTmDu3LklHvP+++9j165dCA0NxeDBg5Gamop169ahSpUqcvvNnDkTNjY28PT0hL29PZKTk/HDDz+gatWqcHd3BwC8/fbbAIDdu3ejcuXKsLa2Rt26dVGjRg3MnDkT06dPR1paGnr06IE33ngDKSkpuHXrFlJSUhTGR0RkSpgcERmZRo0aYe/evYiIiMB3332HjIwM2Nvbo3379gp7Zlq1aoVvv/0WP/zwA4KCglCvXj0EBwfjzJkzuHjxonS/1q1bY9++ffj111+RlpaGGjVqSI8tLACvV68ewsLCsGXLFowcORIFBQUIDw/HgAED8P7776N27dpYt24d5syZg4yMDNjZ2aFx48bo37+/Xp4fIiJdEwmCIBg6CCIiIiJjwZojIiIiIhlMjoiIiIhkMDkiIiIiksHkiIiIiEgGkyMiIiIiGUyOiIiIiGQwOSIiIiKSweSIiIiISAaTIyIiIiIZTI6IiIiIZDA5IiIiIpLB5IiIiIhIxv8Bx+rEEEYdjGAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# grouped boxplots\n",
    "ax = viz.boxplots(comp_data, feature_name=\"Cond\", plot_facets=False,\n",
    "                  y_scale=\"relative\", add_dots=False, figsize=(6,4),\n",
    "                  cmap=\"Set2\")\n",
    "ax.set_xticks(range(0, 11), labels=clusters, rotation=0)\n",
    "\n",
    "\n",
    "# add annotations for sccoda credible effects\n",
    "for i in sccoda_relevant_indices:\n",
    "    x_min = i - .25\n",
    "    x_max = i + .25\n",
    "    y = np.max(prop_df.loc[prop_df[\"cell type\"] == cell_types[i], \"proportion\"] + .05)\n",
    "    ax.plot([x_min, x_max], [y, y], linewidth=1.5, c=\"k\")\n",
    "\n",
    "# add annotations for Wilcoxon rank sum test significant differences\n",
    "for i in significant_indices:\n",
    "    x_min = i - .25\n",
    "    x_max = i + .25\n",
    "    y = np.max(prop_df.loc[prop_df[\"cell type\"] == cell_types[i], \"proportion\"] + .07)\n",
    "    ax.plot([x_min, x_max], [y, y], linewidth=1.5, c=\"darkgreen\", linestyle=\"--\")\n",
    "\n",
    "ax.set_xlabel(\"cluster\")\n",
    "ax.set_ylabel(\"proportion\")\n",
    "legend_1 = ax.legend(loc=\"center left\", title=\"response\\nstatus\", frameon=False, bbox_to_anchor=(1,.8))\n",
    "custom_lines = [Line2D([0], [0], color=\"k\", lw=1.5),\n",
    "                Line2D([0], [0], color=\"darkgreen\", ls=\"--\", lw=1.5)]\n",
    "ax.legend(custom_lines, [\"scCODA\", \"Wilcoxon\"], loc=\"center left\", frameon=False, bbox_to_anchor=(1,.4))\n",
    "ax.add_artist(legend_1)\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(out_dir, \"melanoma_cluster_frequencies.svg\"))\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sccoda",
   "language": "python",
   "name": "sccoda"
  },
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}