analysis/7_fragsys_FUNC_vs_UNK.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EXPLORING THE DIFFERENCES BETWEEN KNOWN FUNCTIONAL AND UNKNOWN LIGAND BINDING SITES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we compare the residues found in binding sites overlapping with functional annotations, or known functional (KF) sites, with those found in unknown functional (UF) sites."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IMPORTING NECESSARY PACKAGES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fragsys_analysis import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## READING INPUT DATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "main_dir = \"./../\"\n",
    "results_dir = os.path.join(main_dir, \"results\")\n",
    "figs_dir = os.path.join(main_dir, \"figs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_bs_ress = pd.read_pickle(os.path.join(results_dir, \"all_bs_ress.pkl\"))\n",
    "all_bss = pd.read_pickle(os.path.join(results_dir, \"all_bss.pkl\"))\n",
    "all_bs_ress[\"RSA\"] = all_bs_ress[\"RSA\"].apply(lambda x: x if x <= 100 else 100)\n",
    "func_df = load_pickle(os.path.join(results_dir, \"prot_func_dict_comb.pkl\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>bs_id</th>\n",
       "      <th>auto</th>\n",
       "      <th>manual</th>\n",
       "      <th>comb</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>H0Y4R8_0_BS0</td>\n",
       "      <td>UF</td>\n",
       "      <td>UF</td>\n",
       "      <td>UF</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O15178_0_BS0</td>\n",
       "      <td>UF</td>\n",
       "      <td>UF</td>\n",
       "      <td>UF</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>O15178_0_BS1</td>\n",
       "      <td>UF</td>\n",
       "      <td>UF</td>\n",
       "      <td>UF</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          bs_id auto manual comb\n",
       "0  H0Y4R8_0_BS0   UF     UF   UF\n",
       "1  O15178_0_BS0   UF     UF   UF\n",
       "2  O15178_0_BS1   UF     UF   UF"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func_df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#func_df[func_df.bs_id.str.contains(\"Q5T0W9\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "UF    249\n",
       "KF     44\n",
       "Name: comb, dtype: int64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func_df.comb.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "func_dict = {}\n",
    "for _, row in func_df.iterrows():\n",
    "    prot, group, bs = row.bs_id.split(\"_\")\n",
    "    func_dict[(prot, group, bs)] = row.comb # combination of automated feature extraction and manual one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "func = []\n",
    "unk = []\n",
    "for k, v in func_dict.items():\n",
    "    prot = k[0]\n",
    "    group = k[1]\n",
    "    bs_id = k[2]\n",
    "    sel_rows = all_bs_ress[(all_bs_ress.protein == prot)&(all_bs_ress.group == group)&(all_bs_ress[bs_id] == 1)]\n",
    "    if v == \"KF\":\n",
    "        func.append(sel_rows)\n",
    "    elif v == \"UF\":\n",
    "        unk.append(sel_rows)\n",
    "func_df = pd.concat(func)\n",
    "unk_df = pd.concat(unk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5767\n",
      "8456\n"
     ]
    }
   ],
   "source": [
    "print(len(func_df))\n",
    "print(len(unk_df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "func_df = func_df.drop_duplicates([\"protein\", \"UniProt_ResNum\", \"UniProt_ResName\"]).reset_index(drop = True)\n",
    "unk_df = unk_df.drop_duplicates([\"protein\", \"UniProt_ResNum\", \"UniProt_ResName\"]).reset_index(drop = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "691\n",
      "2050\n"
     ]
    }
   ],
   "source": [
    "print(len(func_df))\n",
    "print(len(unk_df))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RSA, MES, AND SHENKIN DO NOT FOLLOW A NORMAL DISTRIBUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NormaltestResult(statistic=27.074347971485764, pvalue=1.3209307109831062e-06)\n",
      "NormaltestResult(statistic=64.17466730550996, pvalue=1.1605077781952751e-14)\n"
     ]
    }
   ],
   "source": [
    "print(scipy.stats.normaltest(func_df.log_oddsratio.tolist()))\n",
    "print(scipy.stats.normaltest(unk_df.log_oddsratio.tolist()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NormaltestResult(statistic=78.14182283023712, pvalue=1.0757683664834876e-17)\n",
      "NormaltestResult(statistic=215.75923969317964, pvalue=1.4075907065839778e-47)\n"
     ]
    }
   ],
   "source": [
    "print(scipy.stats.normaltest(func_df.shenkin.tolist()))\n",
    "print(scipy.stats.normaltest(unk_df.shenkin.tolist()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NormaltestResult(statistic=99.0439945451912, pvalue=3.1107839652197576e-22)\n",
      "NormaltestResult(statistic=161.07691023320209, pvalue=1.0534001160039963e-35)\n"
     ]
    }
   ],
   "source": [
    "print(scipy.stats.normaltest(func_df.RSA.tolist()))\n",
    "print(scipy.stats.normaltest(unk_df.RSA.tolist()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FUNCTIONAL VS UNKNOWN: RSA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "691 2050\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAHACAYAAAD5vIKYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABAuklEQVR4nO3de1xUZeI/8M/MMBcQQQXloshFS1G8QrmoeGlbvOWqX9vITczvZvsjXRXZLl4z3ZTWXJc1b6tS5pqX70aZ+WVLLCVbWU0C80KYhoIKX8QSRIQZZp7fH8jsTswMAwdmBubzfr3mVZzzzDnPeSg/nvM853lkQggBIiIiaja5oytARETU1jFMiYiIJGKYEhERScQwJSIikohhSkREJBHDlIiISCKGKRERkUQMUyIiIoncHF0BZ2QwGHDz5k107NgRMpnM0dUhIiIHEULg7t27CAwMhFxu+f6TYWrGzZs3ERQU5OhqEBGRkygqKkKPHj0s7nd4mG7ZsgVvvvkmiouL0b9/f6SkpCAmJsZs2eLiYvz+979HdnY2vvvuOyxYsAApKSkmZXbs2IHdu3fj/PnzAIDIyEisXbsWjz76qM116tixI4C6xvPy8mrehRERUZtXUVGBoKAgYy5Y4tAwPXDgABITE7FlyxaMGDECf/3rXzFhwgRcvHgRPXv2bFC+pqYGXbt2xbJly/DnP//Z7DGPHz+OGTNmYPjw4dBoNFi3bh1iY2Nx4cIFdO/e3aZ61T/a9fLyYpgSEVGjXX4yR050P2zYMAwdOhRbt241bgsPD8fUqVORnJxs9btjxozB4MGDG9yZ/pRer0fnzp2xadMmzJo1y6Z6VVRUwNvbG+Xl5QxTIiIXZmseOGw0r1arRXZ2NmJjY022x8bG4uTJky12nqqqKuh0OnTp0sVimZqaGlRUVJh8iIiIbOWwMC0rK4Ner4efn5/Jdj8/P5SUlLTYeRYvXozu3bvj8ccft1gmOTkZ3t7exg8HHxERUVM4fADST59DCyFa7HWUdevWYd++fTh+/Dg0Go3FckuWLEFSUpLx5/oOZyIioO51ubKyMty5cwcGg8HR1aEWIpfL0alTJ/j6+lp97cUWDgtTX19fKBSKBnehpaWlDe5Wm2P9+vVYu3Ytjh49ioEDB1otq1aroVarJZ+TiNqfGzduYNWqVfj6668dXRVqJZGRkXj11VdtHqRqjsPCVKVSITIyEhkZGZg2bZpxe0ZGBqZMmSLp2G+++SZef/11fPrpp4iKipJaVSJyUVqtFjNmzECnTp2wZs0a9OjRAwqFwtHVohai1+tx/fp1bN68GTNmzMDRo0ehUqmadSyHPuZNSkpCfHw8oqKiEB0dje3bt6OwsBAJCQkA6h6/3rhxA7t37zZ+Jzc3FwBQWVmJW7duITc3FyqVCv369QNQ92h3xYoV2Lt3L0JCQox3vp6envD09LTvBRJRm3bt2jVUVVVh48aNGDx4sKOrQ62gf//+8PPzw5w5c1BYWIjevXs36zgODdO4uDjcvn0bq1evRnFxMSIiIpCeno7g4GAAdZM0FBYWmnxnyJAhxn/Pzs7G3r17ERwcjKtXrwKomwRCq9XiySefNPneypUr8dprr7Xq9RBR+6LX6wHA6pgLavvqf7+1tbXNPobDByDNnTsXc+fONbtv165dDbY19lpsfagSERHZC1eNISIiksjhd6ZERG1R0R0dblfZ7zUZHw85gjop7Xa+1hASEoLExEQkJiba7ZxXr15FaGgocnJyWrXfm2FKRNRERXd0GJpShOpa+83GqnGT4evEIJsDdfbs2Xj33XcbbP/uu++aPcjGVrt27UJiYiLu3Lljsv2rr75Chw4dWvXcjsIwJSJqottVBlTXCvTpqoSHsvXXPK7SCeTfqrsTDupk+/fGjx+Pd955x2Rb165dW7ZyTeDIc7c29pkSETWTh1IGT7W81T/NDWy1Wg1/f3+Tz3PPPYepU6ealEtMTMSYMWOMP48ZMwYLFizAyy+/jC5dusDf37/B2xB37tzBb3/7W/j5+UGj0SAiIgKHDx/G8ePH8d///d8oLy+HTCaDTCYzfjckJMRkcZLCwkJMmTIFnp6e8PLywlNPPYX/+7//M+5/7bXXMHjwYPztb39DSEgIvL298fTTT+Pu3bvGMp988glGjhyJTp06wcfHB0888QSuXLnSrPaSgnemrexymRaVWuuPgjxVMvT2bd6LwkREreHdd99FUlISTp06haysLMyePRsjRozAL37xCxgMBkyYMAF3797Fnj170KtXL1y8eBEKhQLDhw9HSkoKXn31VeTn5wOA2Xf8hRCYOnUqOnTogMzMTNTW1mLu3LmIi4vD8ePHjeWuXLmCgwcP4vDhw/jxxx/x1FNP4Y033sCaNWsAAPfu3UNSUhIGDBiAe/fu4dVXX8W0adOQm5sreYrApmCYtqLLZVoMSSmyqWxOYhADlYha1OHDh02CbMKECTb3WQ4cOBArV64EADz00EPYtGkTPvvsM/ziF7/A0aNHcfr0aeTl5eHhhx8GAISFhRm/6+3tDZlMBn9/f4vHP3r0KL755hsUFBQY50L/29/+hv79++Orr77CI488AqBuXuRdu3YZF+eOj4/HZ599ZgzT6dOnmxw3NTUV3bp1w8WLFxEREWHTtbYEhmkrqr8jtdavUt8X0tjdKxFRU40dO9ZkvegOHTpgyZIlNn33p3OaBwQEoLS0FEDdTHQ9evQwBmlz5OXlISgoyGRRkX79+qFTp07Iy8szhmlISIgxSH9aD6DuznXFihX417/+hbKyMuNCBIWFhQzT9qa+X8U8rkBBRK2jQ4cODUbuyuXyBpPf6HS6Bt9VKk1HDctkMmNQubu7S66bpRXCfrrdWj0AYPLkyQgKCsKOHTsQGBgIg8GAiIgIaLVayXVsCg5AIiJyIV27dkVxcbHJtvo5z201cOBAXL9+HZcuXTK7X6VSGaditKRfv34oLCxEUdG/u8IuXryI8vJyhIeH21SP27dvIy8vD8uXL8fPf/5zhIeH48cff7T9QloQ70yJiJqpSidgj6dLdedpGY899hjefPNN7N69G9HR0dizZw/Onz9vMu95Y0aPHo1Ro0Zh+vTp2LBhA3r37o1vv/0WMpkM48ePR0hICCorK/HZZ59h0KBB8PDwgIeHh8kxHn/8cQwcOBDPPPMMUlJSjAOQRo8ebfNqX507d4aPjw+2b9+OgIAAFBYWYvHixU1qj5bCMCUiaiIfDzk0bjLk32r4eLS1aNxk8PGQ/jBx3LhxWLFiBV5++WVUV1fjN7/5DWbNmoVz58416ThpaWl48cUXMWPGDNy7dw+9e/fGG2+8AQAYPnw4EhISjIuZmFtoRCaT4eDBg5g/fz5GjRoFuVyO8ePH46233rK5DnK5HPv378eCBQsQERGBPn36YOPGjSav+diLTDQ2c7wLqqiogLe3N8rLy+Hl5dXs4+TerEHMlusYEqiy2GdaWWNAzk0tTsztgcGBXKCcyJl8++23mDlzJvbs2YO+ffua7ON0gu2Htd+zrXnAO1MiomYI6qRs0mxE1L5xABIREZFEDFMiIiKJGKZEREQSMUyJiIgkYpgSERFJxDAlIiKSiGFKREQkEd8zJSJqjuoioLbMfudz8wU0QY2XI4dgmBIRNVV1EZA9GDBU2++ccg0QmWv3QB0zZgwGDx6MlJQUu563rWGYEhE1VW1ZXZC69wEUHo2Xl0pfBdzPf3AnbFuYWgrBgwcPYtq0aQ2WYSNpGKZERM2l8AAUHRsvR+0eByAREbmo1157DYMHD8bf/vY3hISEwNvbG08//TTu3r1r8TuffPIJvL29sXv3bgDA7NmzMXXqVKxfvx4BAQHw8fHBvHnzTBYc//HHHzFr1ix07twZHh4emDBhAr777jsAdYuBd+3aFWlpacbygwcPRrdu3Yw/Z2VlQalUorKyEkDdijM7d+7EtGnT4OHhgYceegiHDh1q0bZpKoYpEZELu3LlCg4ePIjDhw/j8OHDyMzMNC6l9lP79+/HU089hd27d2PWrFnG7ceOHcOVK1dw7NgxvPvuu9i1axd27dpl3D979mycOXMGhw4dQlZWFoQQmDhxInQ6HWQyGUaNGoXjx48DqAveixcvQqfT4eLFiwCA48ePIzIyEp6ensZjrlq1Ck899RS++eYbTJw4Ec888wx++OGHlm8gGzFMiYhcmMFgwK5duxAREYGYmBjEx8fjs88+a1Buy5YtSEhIwEcffYQpU6aY7OvcuTM2bdqEvn374oknnsCkSZOMx/juu+9w6NAh7Ny5EzExMRg0aBDee+893LhxAwcPHgRQ179bH6ZffPEFBg0ahMcee8y47fjx4w3WKJ09ezZmzJiB3r17Y+3atbh37x5Onz7dom3TFAxTIiIXFhISgo4d/93vGxAQgNLSUpMyaWlpSExMxJEjRzB27NgGx+jfvz8UCoXZY+Tl5cHNzQ3Dhg0z7vfx8UGfPn2Ql5cHoC5ML1y4gLKyMmRmZmLMmDEYM2YMMjMzUVtbi5MnT2L06NEm5xw4cKDx3zt06ICOHTs2qLc9MUyJiNohLy8vlJeXN9h+584dk0WulUrTBcdlMhkMBtNFzwcPHoyuXbvinXfeMTsK2NoxLI0aFkJAJpMBACIiIuDj44PMzExjmI4ePRqZmZn46quvcP/+fYwcOdLmczoCw5SIqB3q27cvzpw502D7V199hT59+jTpWL169cKxY8fw0UcfYf78+U36br9+/VBbW4tTp04Zt92+fRuXLl1CeHg4ABj7TT/66COcP38eMTExGDBgAHQ6HbZt24ahQ4ea3D07I4YpEVFz6asA/V07fKqaXLW5c+fiypUrmDdvHs6ePYtLly5h8+bNSE1NxUsvvdTk4z388MM4duyY8ZGvrR566CFMmTIFzz//PL788kucPXsWM2fORPfu3U36XseMGYO9e/di4MCB8PLyMgbse++916C/1BnxPVMioqZy862bkeh+vv3OKdfUnddGISEhOHHiBJYtW4bY2FhUV1fj4Ycfxq5du/CrX/2qWVXo06cPPv/8c4wZMwYKhQJ/+tOfbPreO++8g4ULF+KJJ56AVqvFqFGjkJ6ebvKoduzYsdDr9SbBOXr0aBw8eLBBf6kzkglOg9FARUUFvL29UV5ebtK30FS5N2sQs+U6hgSq4Kk2/xCgssaAnJtanJjbA4MD1c0+FxG1vG+//RYzZ87Enj170LdvX9OdnJu33bD2e7Y1D3hnSkTUHJog2Dq1H7V/7DMlIiKSiGFKREQkEcOUiIhIIoYpEZEFcnndH5H/OWk7tT/1v9/633dzMEyJiCzw9/cHAOTk5Di4JtSa6n+/AQEBzT4GR/MSEVng5eWFadOm4a233gIADBkypME0dtR26XQ65OTk4K233sK0adMkzbLEMCUismLJkiUAgI0bNzq4JtRapk2bZvw9NxfDlIjICrlcjmXLlmHBggUoLi526GTq1LLkcjkCAgJaZN5fhikRkQ06duzo9JOtk+NwABIREZFEDFMiIiKJGKZEREQSMUyJiIgkYpgSERFJxDAlIiKSiGFKREQkEcOUiIhIIoeH6ZYtWxAaGgqNRoPIyEicOHHCYtni4mL8+te/Rp8+fSCXy5GYmGi2XFpaGvr16we1Wo1+/frhww8/bKXaExEROThMDxw4gMTERCxbtgw5OTmIiYnBhAkTUFhYaLZ8TU0NunbtimXLlmHQoEFmy2RlZSEuLg7x8fE4e/Ys4uPj8dRTT+HUqVOteSlEROTCZEII4aiTDxs2DEOHDsXWrVuN28LDwzF16lQkJydb/e6YMWMwePBgpKSkmGyPi4tDRUUF/vGPfxi3jR8/Hp07d8a+fftsqldFRQW8vb1RXl4OLy8v2y/oJ3Jv1iBmy3UMCVTBU23+7y2VNQbk3NTixNweGByobva5iIio5dmaBw67M9VqtcjOzkZsbKzJ9tjYWJw8ebLZx83KympwzHHjxlk9Zk1NDSoqKkw+REREtnJYmJaVlUGv18PPz89ku5+fH0pKSpp93JKSkiYfMzk5Gd7e3sZPUFBQs89PRESux+EDkGQymcnPQogG21r7mEuWLEF5ebnxU1RUJOn8RETkWhy2BJuvry8UCkWDO8bS0tIGd5ZN4e/v3+RjqtVqqNXsryQiouZx2J2pSqVCZGQkMjIyTLZnZGRg+PDhzT5udHR0g2MeOXJE0jGJiIisceji4ElJSYiPj0dUVBSio6Oxfft2FBYWIiEhAUDd49cbN25g9+7dxu/k5uYCACorK3Hr1i3k5uZCpVKhX79+AICFCxdi1KhR+OMf/4gpU6bgo48+wtGjR/Hll1/a/fqIiMg1ODRM4+LicPv2baxevRrFxcWIiIhAeno6goODAdRN0vDTd06HDBli/Pfs7Gzs3bsXwcHBuHr1KgBg+PDh2L9/P5YvX44VK1agV69eOHDgAIYNG2a36yIiItfi0PdMnRXfMyUiIqANvGdKRETUXjBMiYiIJGKYEhERScQwJSIikohhSkREJBHDlIiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCRimBIREUnEMCUiIpKIYUpERCQRw5SIiEgihikREZFEDl0cnP4tv1Rrdb+nSobevio71YaIiJqCYepgigfPBua8X9po2ZzEIAYqEZETYpg6mLtSjqgeKugNlstU6QTyb+lQqRX2qxgREdmMYeoE3JWNdV1bSVoiInI4DkAiIiKSiGFKREQkEcOUiIhIIoYpERGRRAxTIiIiiRimREREEjFMiYiIJGKYEhERScQwJSIikohhSkREJBHDlIiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCRimBIREUnEMCUiIpKIYUpERCQRw5SIiEgihikREZFEDFMiIiKJGKZEREQSMUyJiIgkYpgSERFJxDAlIiKSiGFKREQkEcOUiIhIIoYpERGRRAxTIiIiiRimREREEjFMiYiIJGKYEhERScQwJSIikohhSkREJJHDw3TLli0IDQ2FRqNBZGQkTpw4YbV8ZmYmIiMjodFoEBYWhm3btjUok5KSgj59+sDd3R1BQUFYtGgRqqurW+sSiIjIxTk0TA8cOIDExEQsW7YMOTk5iImJwYQJE1BYWGi2fEFBASZOnIiYmBjk5ORg6dKlWLBgAdLS0oxl3nvvPSxevBgrV65EXl4eUlNTceDAASxZssRel0VERC7GzZEn37BhA5577jnMmTMHQN0d5aeffoqtW7ciOTm5Qflt27ahZ8+eSElJAQCEh4fjzJkzWL9+PaZPnw4AyMrKwogRI/DrX/8aABASEoIZM2bg9OnTFutRU1ODmpoa488VFRUtdYlEROQCHHZnqtVqkZ2djdjYWJPtsbGxOHnypNnvZGVlNSg/btw4nDlzBjqdDgAwcuRIZGdnG8Pz+++/R3p6OiZNmmSxLsnJyfD29jZ+goKCpFwaERG5GIeFaVlZGfR6Pfz8/Ey2+/n5oaSkxOx3SkpKzJavra1FWVkZAODpp5/GH/7wB4wcORJKpRK9evXC2LFjsXjxYot1WbJkCcrLy42foqIiiVdHRESuxKGPeQFAJpOZ/CyEaLCtsfL/uf348eNYs2YNtmzZgmHDhuHy5ctYuHAhAgICsGLFCrPHVKvVUKvVUi6DiIhcmMPC1NfXFwqFosFdaGlpaYO7z3r+/v5my7u5ucHHxwcAsGLFCsTHxxv7YQcMGIB79+7ht7/9LZYtWwa53OEDmImIqJ1xWLKoVCpERkYiIyPDZHtGRgaGDx9u9jvR0dENyh85cgRRUVFQKpUAgKqqqgaBqVAoIIQw3sUSERG1JIfepiUlJWHnzp14++23kZeXh0WLFqGwsBAJCQkA6voyZ82aZSyfkJCAa9euISkpCXl5eXj77beRmpqKF1980Vhm8uTJ2Lp1K/bv34+CggJkZGRgxYoV+OUvfwmFQmH3ayQiovbPoX2mcXFxuH37NlavXo3i4mJEREQgPT0dwcHBAIDi4mKTd05DQ0ORnp6ORYsWYfPmzQgMDMTGjRuNr8UAwPLlyyGTybB8+XLcuHEDXbt2xeTJk7FmzRq7Xx8REbkGmeCzzwYqKirg7e2N8vJyeHl5Nfs4uTdrMHvXPxHlr4W7yvJDgPuGDripC7G4v7LGgJybWpyY2wODAzlQiojIXmzNA4eP5m3PVNoryI2ZbFPZOQVHrAYqERE5L4ZpK1KISgBAQXUYdMLDbBm1/D6C1VfgLr9nz6oREVELYpjaQbXBHVp0ML/TYPtx8ku1Vvd7qmTo7atqQs2IiKglMEzbAMWD7tY575c2WjYnMYiBSkRkZwzTNsBdKUdUDxX0Vu5iq3QC+bd0qNRyPBkRkb0xTNsId2VjrwQ34XkxERG1KM6tR0REJBHDlIiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCRimBIREUnEMCUiIpKIYUpERCQRw5SIiEgihikREZFEDFMiIiKJGKZEREQSMUyJiIgkYpgSERFJxDAlIiKSiGFKREQkEcOUiIhIIoYpERGRRG6OrgC1rPxSrdX9nioZevuq7FQbIiLXwDBtJxQPnjHMeb+00bI5iUEMVCKiFsQwbSfclXJE9VBBb7BcpkonkH9Lh0qtsF/FiIhcAMO0HXFXNtYFbiVpiYio2TgAiYiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCRimBIREUnEMCUiIpKIYUpERCQRw5SIiEgihikREZFEDFMiIiKJGKZEREQSMUyJiIgkYpgSERFJxDAlIiKSiGFKREQkEcOUiIhIIoYpERGRRAxTIiIiiRimREREEjFMiYiIJGKYEhERSSQ5TIUQEEK0RF2IiIjapGaHaWpqKiIiIqDRaKDRaBAREYGdO3c2+ThbtmxBaGgoNBoNIiMjceLECavlMzMzERkZCY1Gg7CwMGzbtq1BmTt37mDevHkICAiARqNBeHg40tPTm1w3IiIiW7g150srVqzAn//8Z8yfPx/R0dEAgKysLCxatAhXr17F66+/btNxDhw4gMTERGzZsgUjRozAX//6V0yYMAEXL15Ez549G5QvKCjAxIkT8fzzz2PPnj345z//iblz56Jr166YPn06AECr1eIXv/gFunXrhvfffx89evRAUVEROnbs2JxLJSIiapRMNOMZra+vL9566y3MmDHDZPu+ffswf/58lJWV2XScYcOGYejQodi6datxW3h4OKZOnYrk5OQG5V955RUcOnQIeXl5xm0JCQk4e/YssrKyAADbtm3Dm2++iW+//RZKpbKplwYAqKiogLe3N8rLy+Hl5dWsYwBAfsFp9Lk+FnlV/aGFp9ky7rJ7eNj9POZf+xBXavo3+1y2qKwxIOemFifm9sDgQHWrnouIqD2wNQ+a9ZhXr9cjKiqqwfbIyEjU1tbadAytVovs7GzExsaabI+NjcXJkyfNficrK6tB+XHjxuHMmTPQ6XQAgEOHDiE6Ohrz5s2Dn58fIiIisHbtWuj1eot1qampQUVFhcmHiIjIVs0K05kzZ5rcTdbbvn07nnnmGZuOUVZWBr1eDz8/P5Ptfn5+KCkpMfudkpISs+Vra2uNd8Pff/893n//fej1eqSnp2P58uX405/+hDVr1lisS3JyMry9vY2foKAgm66BiIgIaGafKVA3AOnIkSP42c9+BgD417/+haKiIsyaNQtJSUnGchs2bLB6HJlMZvKzEKLBtsbK/+d2g8GAbt26Yfv27VAoFIiMjMTNmzfx5ptv4tVXXzV7zCVLlpjUuaKigoFKREQ2a1aYnj9/HkOHDgUAXLlyBQDQtWtXdO3aFefPnzeWsxaKvr6+UCgUDe5CS0tLG9x91vP39zdb3s3NDT4+PgCAgIAAKJVKKBQKY5nw8HCUlJRAq9VCpVI1OK5arYZazT5EIiJqnmaF6bFjxySfWKVSITIyEhkZGZg2bZpxe0ZGBqZMmWL2O9HR0fj4449Nth05cgRRUVHGwUYjRozA3r17YTAYIJfXPcW+dOkSAgICzAYpERGRVA6dASkpKQk7d+7E22+/jby8PCxatAiFhYVISEgAUPf4ddasWcbyCQkJuHbtGpKSkpCXl4e3334bqampePHFF41lXnjhBdy+fRsLFy7EpUuX8L//+79Yu3Yt5s2bZ/frIyIi19DsPtOWEBcXh9u3b2P16tUoLi5GREQE0tPTERwcDAAoLi5GYWGhsXxoaCjS09OxaNEibN68GYGBgdi4caPxHVMACAoKwpEjR7Bo0SIMHDgQ3bt3x8KFC/HKK6/Y/fqIiMg1NOs90/aO75kSERHQyu+ZEhER0b8xTImIiCRimBIREUnEMCUiIpKIYUpERCQRw5SIiEgihikREZFEDFMiIiKJHDoDEjlGfqnW6n5PlQy9fTmPMRGRrRimLkTx4DnEnPdLGy2bkxjEQCUishHD1IW4K+WI6qGC3mC5TJVOIP+WDpVazjJJRGQrhqmLcVc21k1uJWmJiMgsDkAiIiKSiGFKREQkEcOUiIhIIoYpERGRRAxTIiIiiRimREREEjFMiYiIJGKYEhERScQwJSIikogzIDmJINUVq/vvGzrgpi7EPpUhIqImYZg6mB4KAMDLAS82WnZOwREGKhGRE2KYOphWaJB3fxAU0Fsso5bfR7D6Ctzl9+xYMyIishXD1AlohcZ6Ac49T0Tk1DgAiYiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCTiaF5qNZfLtKjUCqtlPFUy9PZV2alGREStg2FKreJymRZDUopsKpuTGMRAJaI2jWFKraL+jrRPVyU8lDKzZap0Avm3dI3evRIROTuGaRvSFufv9VDK4Km21DXP2SiIqH1gmLYBnL+XiMi5MUzbAGecv7exwUX5pVq71IOIyBkwTNsIZ5q/tymDixR8+YqIXADDlJrMlsFFQF2QuiuZpkTU/jFM2xl7DVLq5XENgztq4a6yHJbOOCCKiKg1MEzbCXsOUlJpryA3ZrJNZTkgiohcAcO0nbDnICWFqAQAFFSHQSc8WvVcRERtAcO0HbH3IKVqgzu06GCXcxEROTOODiEiIpKIYUpERCQRw5SIiEgihikREZFEDFMiIiKJGKZEREQSMUyJiIgk4num1EAvj2twr7kNVCrN7tdoL9m5RkREzo1hSiZ6qK5ib8xk4DrqPmYEP/inQSgAy/PcExG5DIYpmfB4MP3ffeXDcFebn93onlbgQqkebm4aKBimRESO7zPdsmULQkNDodFoEBkZiRMnTlgtn5mZicjISGg0GoSFhWHbtm0Wy+7fvx8ymQxTp05t4Vq3fwaZB6DoaPZjkHui2tDI1IVERC7EoWF64MABJCYmYtmyZcjJyUFMTAwmTJiAwsJCs+ULCgowceJExMTEICcnB0uXLsWCBQuQlpbWoOy1a9fw4osvIiYmprUvg4iIXJxDw3TDhg147rnnMGfOHISHhyMlJQVBQUHYunWr2fLbtm1Dz549kZKSgvDwcMyZMwe/+c1vsH79epNyer0ezzzzDFatWoWwsDB7XAoREbkwh/WZarVaZGdnY/HixSbbY2NjcfLkSbPfycrKQmxsrMm2cePGITU1FTqdDkpl3ejT1atXo2vXrnjuuecafWwMADU1NaipqTH+XFFR0dTLaVOsLSDeTf69HWtCRNQ+OCxMy8rKoNfr4efnZ7Ldz88PJSUlZr9TUlJitnxtbS3KysoQEBCAf/7zn0hNTUVubq7NdUlOTsaqVauafA1tTVMWEBcPyhIRUeMcPppXJjMdDiqEaLCtsfL12+/evYuZM2dix44d8PX1tbkOS5YsQVJSkvHniooKBAUF2fz9tsKWBcT1QuCuVo5wf3c71oyIqG1zWJj6+vpCoVA0uAstLS1tcPdZz9/f32x5Nzc3+Pj44MKFC7h69SomT55s3G8w1K1S7ebmhvz8fPTq1avBcdVqNdRqtdRLahMaW0BcbxCoNghU6Syv7m1tX1t3uUyLSq2wWsZTJUNvX5WdakREbYHDwlSlUiEyMhIZGRmYNm2acXtGRgamTJli9jvR0dH4+OOPTbYdOXIEUVFRUCqV6Nu3L86dO2eyf/ny5bh79y7+8pe/tMu7zZZWf+Off0tnc9n24nKZFkNSimwqm5MYxEAlIiOHPuZNSkpCfHw8oqKiEB0dje3bt6OwsBAJCQkA6h6/3rhxA7t37wYAJCQkYNOmTUhKSsLzzz+PrKwspKamYt++fQAAjUaDiIgIk3N06tQJABpsJ/PkMhk6qgFh/eYMMlld2fak/o60T1clPJTmr61KJ5B/S9fo3SsRuRaHhmlcXBxu376N1atXo7i4GBEREUhPT0dwcN2EdcXFxSbvnIaGhiI9PR2LFi3C5s2bERgYiI0bN2L69OmOuoR2SS6T2W2awMbmAQZQN1mEe2/7VAiAh1IGT7Wlt8ba7yNuImo+hw9Amjt3LubOnWt2365duxpsGz16NL7++mubj2/uGOQcbJkH2CjqrF0DlYioKRwepuS6bJkHGPoq4H4+oL9rx5oRETUNw5QczjgPMBFRG+Xwie6JiIjaOoYpERGRRAxTIiIiiRimREREEnEAErUNVfnW99v5XVQiov/EMCUn9+A/0fznGi/Kd1GJyEEYpuTcFO6A5yMAai2X4buoRORgDFNyfgrnWg7OGadAJCLHYpgSNQGnQCQicximRE3AKRCJyByGKVEzcApEIvpPfM+UiIhIIt6ZUqsKUl2xuK+b/Hs71oSIqPUwTKlV6KEAALwc8GKjZcWDskREbRXDlFqFVmiQd38QFNBbLKMXAne1coT7O9erL0RETcUwpVajFRqr+/UGgWqDsFNtbNPL4xoe1mjhrjI/nICPponIHIYpOb0qnQF6g+X9coOAhZdUmkSlvYLcmMk2lbXp0TTnEyZyGQxTcmpVOgPOXK+xWsZTocVQb6Dox1oEeTb/XApRCQAoqA6DTniYLWPbo2nOJ0zkahim5NTq70g9lDLIZebLaGR1O6pqW+aRcbXBHVoL97o2PZrmfMJELodhSm2CXAYoLKSphYx1LCebT5iIWhfDlKgZqnRWOnEBKOSAh5JzohC5CoYpURM8eKKM/Fu6RstG9VAzUIlcBMOUXMf9y1b7KDXaS40eQi6ToaMaEFa6TQ0CqNIJqyOQiah9YZiSa7h/GTgzyGqR4Af/NAiF1Y5YuUxmvaPWyd6dJaLWxzAl11B/R+reB1CYf+3lnlbgQqkebm4aKJxyVBMROSuGKbmUe3p3GIT5116q9AZUG3SQ8KoqEbkohim5hKIfaxEEIP+WFpV665NAyJzprrSRfl4AnEmJyAkwTMkl1E/ooHGTwdozXJnsQZ+oM7Chn9eIMykRORTDlFyKTAYonCUsG2NDPy9nUiJyDgxTImen8Kh7lEtETothSg5nbTahxmYaIiJyBgxTcpimzCbUVp7MEpFrYpiSw9gymxDgZIOCWpq1NU8bWw+1KWwZFQxwZDBRMzFMyaEanU2oCTTaS0Cl0vI+p9KENU+l/m/alFHBANDvfwB1oOX9DFyiBhim1OYZhAIAEFz6/4BS82VsnSrQbmxZ8xQA4CZ9OTdbRgUDgEELVF0ALj7V+DH5Kg6RCYYptXk1QoPTdwagfzcFOqjMJ6VTThVo7zVPGxsVrAAXNSdqJoYptQvVBg0McnXdQqJmGOQGVBtq2u9UgS3V98pFzYmahWFK1KbZse+ViCzi/11EbZk9+16JyCKGKbUbzjb5Q2PnVMgBD6X5x9JNOo9BDb1B3fi5FDYcS2ewuqi53CBgfs0dItfGMKU2z9kmf2hKfaJ6qCUFapXOgDPXra+CY+u5bDmWp0KLod4PVuFptx3QRE3HMKU2z9kmf7ClPgYBVOmE1btAW9R/30Mpg9zCpdl6LluOpXnQfvWr8LQ6LkFHbQTDlNqFlpz8oSU0Wh9Dy4aRXAYorKVpCx3Lrk3MJeioDWGYEpFz4hJ01IYwTImoyaxN3QigZR+9cgk6agMYpkRkM1umbjTio1dyIQxTIrKZLVM38tEruSKGKRE1SWNTNxK5IoYpkQPZa2IHImpdDFMiB7DnxA4O09gE+y05SMme5yIyg2FK5AD2nNjB/pow+b7kQUr2PBeRZQ7/6+6WLVsQGhoKjUaDyMhInDhxwmr5zMxMREZGQqPRICwsDNu2bTPZv2PHDsTExKBz587o3LkzHn/8cZw+fbo1L4GoWeQyGRRyyx9LczA4vfrJ9z2HWP6496krK3WQkj3PRWSFQ8P0wIEDSExMxLJly5CTk4OYmBhMmDABhYWFZssXFBRg4sSJiImJQU5ODpYuXYoFCxYgLS3NWOb48eOYMWMGjh07hqysLPTs2ROxsbG4ceOGvS6LiBTudY9WLX4sTMLg7OcissChj3k3bNiA5557DnPmzAEApKSk4NNPP8XWrVuRnJzcoPy2bdvQs2dPpKSkAADCw8Nx5swZrF+/HtOnTwcAvPfeeybf2bFjB95//3189tlnmDVrltl61NTUoKbm3xN8V1RUtMTlEbUIZ1sNh4gacliYarVaZGdnY/HixSbbY2NjcfLkSbPfycrKQmxsrMm2cePGITU1FTqdDkplwxlZqqqqoNPp0KVLF4t1SU5OxqpVq5pxFUStx9lWwyEiyxz2mLesrAx6vR5+fn4m2/38/FBSUmL2OyUlJWbL19bWoqyszOx3Fi9ejO7du+Pxxx+3WJclS5agvLzc+CkqKmri1RC1vLpBSjJ4qqx/OqpldlkNh4gsc/hoXtlP/hAQQjTY1lh5c9sBYN26ddi3bx+OHz8OjUZj8ZhqtRpqtfXFlYkcwdlWw2lx1l5paex1FyIn4rAw9fX1hUKhaHAXWlpa2uDus56/v7/Z8m5ubvDx8THZvn79eqxduxZHjx7FwIEDW7byRCRRE15pcfzf+Yka5bD/SlUqFSIjI5GRkYFp06YZt2dkZGDKlClmvxMdHY2PP/7YZNuRI0cQFRVl0l/65ptv4vXXX8enn36KqKio1rkAImq++ldaUNtIQbe6skROzqF/5UtKSkJ8fDyioqIQHR2N7du3o7CwEAkJCQDq+jJv3LiB3bt3AwASEhKwadMmJCUl4fnnn0dWVhZSU1Oxb98+4zHXrVuHFStWYO/evQgJCTHeyXp6esLT09P+F0lE5jEkqR1xaJjGxcXh9u3bWL16NYqLixEREYH09HQEBwcDAIqLi03eOQ0NDUV6ejoWLVqEzZs3IzAwEBs3bjS+FgPUTQKh1Wrx5JNPmpxr5cqVeO211+xyXUTUBt2/3PjEDpyWkCyQCWFtQjPXVFFRAW9vb5SXl8PLy6vZx8kvOI0+18cir6o/tOBdMTmG3iBQqRUYEqhGR7XlAfx3awzIuVkDT1XdDExSjuVU9HeByhxgyJd1MyKZc/8ycGaQbcfjtIQuxdY8YM8+EVH9Hal7H8szJtm4TuvlMi0qtdbvUTxVMvT2VTWnpuSkGKZERPUUHnWPcpvpcpkWQ1Jse089JzGIgdqOMEyJiJrCyvuv+rs69PKohFuH3vBQmn9UXqUTyL+la/TuldoWhikRkU0afze2D4DcGODX+Z/gjizMQinOp9weMUyJiGxhw7ux92vuwV13CR7ye7jDG0+XwjAlIrJVI+/GGmRMUFfVRsa2ExEROS+GKRERkUQMUyIiIonYZ0rkIqp01keRNra/pVXpDNA3ckqFHPBQ8u/85PwYpkTtXP1Sv/m3dE0q35qqdAacuV5jU9moHuo2F6jBmu/hbjBf5/syAyo8VAB62LdS1KoYpkTtnFwmQ0c1YMss3DLZgwXJW1n9HamHUgYL0wDDIOomOGjs7tWZCCgAACuCXrJe8CHgovYrAP2knZCT8zsNhimRC5DLZIAd7jibSi6DxUn1YWh7r5kIuTtO3xmAjioDFBb+UqKUVSFU8z0UolLayTg5v1NhmBJRkzXWv9qSfZ32PFdLqDZo4GawsvJOSy3U1YKT87dJTnZXzjAlIps1pf9Val+nPc/VpkmcnL9NcsK7coYpEdnMlv7XlurrtOe5qI1xwrtyhikRNUmj/a8t2Ndpz3NRG+REd+Uu+lyEiIio5fDOlIhahbWBQ/aeIKLupJbXIbW6z1EaG2DjjHV2YQxTImpRTRk4ZI8JImxZh7RBWUdrygAbZ6mzi+NvgYhalK2TRNhrgghb1iGt49boEmt2Y8sAGwBOVWcXxzAlohbndJNEtNXAcaIBNmQdw5SIyAE02ktApdL8TvaHtjkMUyJq89rSLEkGUTd/b3Dp/wNKGyvNP6LbCv6miKjNaouzJGmhwek7A+Am01stVysUiPBQw0PRAidt7E6Xk+FLxjAlojarLc6SJJfJoFS6N1rn6hapcxNGMnMyfEkYpkTUprXFWZLsVmdbRjK358nw7YhhSkTkpFqkL9jWkcwtMejJhR8XM0yJiJyMffuCmzKphQ1c9HExw5SIyMnYtS/Y5kktGuHij4sZpkRETsiufcFtdVILJ+L4ceJERERtHO9MiYhsVKUzWH2s6pDVcFzYPa2AQW6+zeUGgQ52rAvDlIjIBlU6A85cr7GprH1Ww3FdRT/WIghA/i0tKvXmfyeeCi2Gej8o69n6dWKYEhHZoP6O1EMpg9xKWNptNRwXVlVb11+scZMBCvNtrXnwO6gv29oYpkRETSCXAQpraUp2I5MBCgt/cbH3b4hhSkQE9oc2prH2MfZRuug8wAxTInJ57A+1zpb20cj1eLQTXHYeYIYpEbk89odaZ0v7GIQ7Tt8ZgP7dFOigslCoHU/swDAlInqA/aHWWW0fg0CVQQODXF03abCLcb0rJiIiamG8MyUil2BtABEHFznX4KvLZVpUai2/0lL6gw597FgfWzBMiahda8oKLC7YHep0g68ul2nx5PYv4am4Z7FMH88CxA58UB8nWa6WYUpE7ZotK7AAHFzkLIOvtJWXkRsz2aaywokizHlqQkTUShpdgYWcZvCVQlQCAAqqw6ATHhbL6aGAVmjsVa1GMUyJiMjpVBvcobXrVPXSMEyJiNowqYODXH3wVUthmBIRtUFNGVjVlONR8zBMiYjaIFsHVtnCVQdftSSGKRFRG8WBVc6DYUpERC3GWh+sLSvLaLSXWr5SdsAwJSIiyWzpw7VlZZngB/80CEWbuut2+Ny8W7ZsQWhoKDQaDSIjI3HixAmr5TMzMxEZGQmNRoOwsDBs27atQZm0tDT069cParUa/fr1w4cfftha1SciItT34crgqbL8kSvqVpa5pxoMeA4x+7mnGozTdwagxoneIbWFQ8P0wIEDSExMxLJly5CTk4OYmBhMmDABhYWFZssXFBRg4sSJiImJQU5ODpYuXYoFCxYgLS3NWCYrKwtxcXGIj4/H2bNnER8fj6eeegqnTp2y12UREbkkuUwGhdzyRy4Dqg0aGOSedYuEm/kY5J6oNrStIAUAmRAtMRaseYYNG4ahQ4di69atxm3h4eGYOnUqkpOTG5R/5ZVXcOjQIeTl5Rm3JSQk4OzZs8jKygIAxMXFoaKiAv/4xz+MZcaPH4/OnTtj3759NtWroqIC3t7eKC8vh5eXV3MvD/kFp9Hn+ljkVfWHFp7NPg4RUXugNwhUagWGBKrRUW3+Xu5ujQE5N2vgqZJJmpFJhUqEe1xAfo9j6BP6aLOPY2seOKzPVKvVIjs7G4sXLzbZHhsbi5MnT5r9TlZWFmJjY022jRs3DqmpqdDpdFAqlcjKysKiRYsalElJSbFYl5qaGtTU/Hui5/LycgB1jShF5d1KVNwT0N6vhFboJR2LiKitMwgAeoGi/5NDozQflNU6AdQYoKuVQS+lz1R2HxVC1P05LOHP8vrvNnbf6bAwLSsrg16vh5+fn8l2Pz8/lJSUmP1OSUmJ2fK1tbUoKytDQECAxTKWjgkAycnJWLVqVYPtQUFBtl5OI/IaL0JERK3g5y1ylLt378Lb29vifoeP5pX95EVhIUSDbY2V/+n2ph5zyZIlSEpKMv5sMBjwww8/wMfHx+r3GlNRUYGgoCAUFRVJelzc3rBdLGPbmMd2sYxtY15LtYsQAnfv3kVgYKDVcg4LU19fXygUigZ3jKWlpQ3uLOv5+/ubLe/m5gYfHx+rZSwdEwDUajXUarXJtk6dOtl6KY3y8vLif+RmsF0sY9uYx3axjG1jXku0i7U70noOG82rUqkQGRmJjIwMk+0ZGRkYPny42e9ER0c3KH/kyBFERUVBqVRaLWPpmERERFI59DFvUlIS4uPjERUVhejoaGzfvh2FhYVISEgAUPf49caNG9i9ezeAupG7mzZtQlJSEp5//nlkZWUhNTXVZJTuwoULMWrUKPzxj3/ElClT8NFHH+Ho0aP48ssvHXKNRETkAoSDbd68WQQHBwuVSiWGDh0qMjMzjfueffZZMXr0aJPyx48fF0OGDBEqlUqEhISIrVu3Njjm3//+d9GnTx+hVCpF3759RVpaWmtfhlnV1dVi5cqVorq62iHnd1ZsF8vYNuaxXSxj25hn73Zx6HumRERE7YHDpxMkIiJq6ximREREEjFMiYiIJGKYEhERScQwbSVNXVquvUlOTsYjjzyCjh07olu3bpg6dSry800XBBZC4LXXXkNgYCDc3d0xZswYXLhwwUE1dpzk5GTIZDIkJiYat7lq29y4cQMzZ86Ej48PPDw8MHjwYGRnZxv3u2q71NbWYvny5QgNDYW7uzvCwsKwevVqGAz/XojbFdrmiy++wOTJkxEYGAiZTIaDBw+a7LelDWpqajB//nz4+vqiQ4cO+OUvf4nr169Lr5xdxgy7mP379wulUil27NghLl68KBYuXCg6dOggrl275uiq2c24cePEO++8I86fPy9yc3PFpEmTRM+ePUVlZaWxzBtvvCE6duwo0tLSxLlz50RcXJwICAgQFRUVDqy5fZ0+fVqEhISIgQMHioULFxq3u2Lb/PDDDyI4OFjMnj1bnDp1ShQUFIijR4+Ky5cvG8u4YrsIIcTrr78ufHx8xOHDh0VBQYH4+9//Ljw9PUVKSoqxjCu0TXp6uli2bJlIS0sTAMSHH35ost+WNkhISBDdu3cXGRkZ4uuvvxZjx44VgwYNErW1tZLqxjBtBY8++qhISEgw2da3b1+xePFiB9XI8UpLSwUA43vEBoNB+Pv7izfeeMNYprq6Wnh7e4tt27Y5qpp2dffuXfHQQw+JjIwMMXr0aGOYumrbvPLKK2LkyJEW97tquwghxKRJk8RvfvMbk23/9V//JWbOnCmEcM22+WmY2tIGd+7cEUqlUuzfv99Y5saNG0Iul4tPPvlEUn34mLeF1S8t99Ol4qwtLecK6pe169KlC4C6hd5LSkpM2kmtVmP06NEu007z5s3DpEmT8Pjjj5tsd9W2OXToEKKiovCrX/0K3bp1w5AhQ7Bjxw7jfldtFwAYOXIkPvvsM1y6dAkAcPbsWXz55ZeYOHEiANdum3q2tEF2djZ0Op1JmcDAQEREREhuJ4evGtPeNGdpufZOCIGkpCSMHDkSERERAGBsC3PtdO3aNbvX0d7279+Pr7/+Gl999VWDfa7aNt9//z22bt2KpKQkLF26FKdPn8aCBQugVqsxa9Ysl20XAHjllVdQXl6Ovn37QqFQQK/XY82aNZgxYwYA1/1v5j/Z0gYlJSVQqVTo3LlzgzJS/3xmmLaSpi4D15797ne/wzfffGN2fmRXbKeioiIsXLgQR44cgUajsVjO1drGYDAgKioKa9euBQAMGTIEFy5cwNatWzFr1ixjOVdrFwA4cOAA9uzZg71796J///7Izc1FYmIiAgMD8eyzzxrLuWLb/FRz2qAl2omPeVtYc5aWa8/mz5+PQ4cO4dixY+jRo4dxu7+/PwC4ZDtlZ2ejtLQUkZGRcHNzg5ubGzIzM7Fx40a4ubkZr9/V2iYgIAD9+vUz2RYeHo7CwkIArv3fzEsvvYTFixfj6aefxoABAxAfH49FixYhOTkZgGu3TT1b2sDf3x9arRY//vijxTLNxTBtYc1ZWq49EkLgd7/7HT744AN8/vnnCA0NNdkfGhoKf39/k3bSarXIzMxs9+3085//HOfOnUNubq7xExUVhWeeeQa5ubkICwtzybYZMWJEg9enLl26hODgYACu/d9MVVUV5HLTP64VCoXx1RhXbpt6trRBZGQklEqlSZni4mKcP39eejtJGr5EZtW/GpOamiouXrwoEhMTRYcOHcTVq1cdXTW7eeGFF4S3t7c4fvy4KC4uNn6qqqqMZd544w3h7e0tPvjgA3Hu3DkxY8aMdjeU31b/OZpXCNdsm9OnTws3NzexZs0a8d1334n33ntPeHh4iD179hjLuGK7CFG3glb37t2Nr8Z88MEHwtfXV7z88svGMq7QNnfv3hU5OTkiJydHABAbNmwQOTk5xtcObWmDhIQE0aNHD3H06FHx9ddfi8cee4yvxjgza0vLuQIAZj/vvPOOsYzBYBArV64U/v7+Qq1Wi1GjRolz5845rtIO9NMwddW2+fjjj0VERIRQq9Wib9++Yvv27Sb7XbVdKioqxMKFC0XPnj2FRqMRYWFhYtmyZaKmpsZYxhXa5tixY2b/XHn22WeFELa1wf3798Xvfvc70aVLF+Hu7i6eeOIJUVhYKLluXIKNiIhIIvaZEhERScQwJSIikohhSkREJBHDlIiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCRimBIREUnEMCVqB2bPng2ZTAaZTAY3Nzf07NkTL7zwQoPVMXJycvDEE0+gW7du0Gg0CAkJQVxcHMrKyhocc+3atVAoFHjjjTcaPX9NTQ3i4+Ph5eWFPn364PPPPzfZv27dOsyfP1/aRRI5MYYpUTsxfvx4FBcX4+rVq9i5cyc+/vhjzJ0717i/tLQUjz/+OHx9ffHpp58iLy8Pb7/9NgICAlBVVdXgeO+88w5efvllvP32242ee/v27cjOzkZWVhaef/55zJgxA/UzlRYUFGDnzp1Ys2ZNy10skZPh4uBE7YRarTau6dijRw/ExcVh165dxv0nT55ERUUFdu7cCTe3uv/1Q0ND8dhjjzU4VmZmJu7fv4/Vq1dj9+7d+OKLLzBq1CiL587Ly8Mvf/lL9O/fH2FhYXjppZdQVlaGrl274oUXXsAf//hHeHl5tewFEzkR3pkStUPff/89PvnkEyiVSuM2f39/1NbW4sMPP0Rj61ukpqZixowZUCqVmDFjBlJTU62WHzRoEL788kvcv38fn376KQICAuDr64s9e/ZAo9Fg2rRpLXJdRE5L8rozRORwzz77rFAoFKJDhw5Co9EYl6basGGDSbmlS5cKNzc30aVLFzF+/Hixbt06UVJSYlKmvLxceHh4iNzcXCGEEDk5OcLDw0OUl5dbPL9WqxVz584VISEhIioqSpw4cULcvn1bhIWFiWvXrolly5aJXr16idjYWHH9+vWWbwAiB+MSbETtwOzZs3Hjxg1s3boVVVVV2LlzJy5duoTDhw8bH+nWu337Nj7//HP861//wsGDB/HDDz/giy++wIABAwAA27Ztw+bNm3Hu3DnjdwYMGID58+fjt7/9bZPqNGTIEISGhmLp0qU4deoU1q1bh/PnzyMtLa1lLpzISTBMidqB2bNn486dOzh48KBx29ixYzFy5Ej84Q9/sPg9rVaLIUOGICoqCu+++y4A4NFHH8WZM2cgl/+7F8hgMOCRRx7BqVOnbKrP559/jsWLFyMrKwsvvfQS3NzcsG7dOly4cAGjRo3C7du3m3ehRE6KA5CI2qmVK1diwoQJeOGFFxAYGGi2jEqlQq9evXDv3j0AwLlz53DmzBkcP34cXbp0MZa7c+cORo0ahfPnzyMiIsLqeaurqzFv3jzs3bsXCoUCer3e2Eer0+mg1+tb6AqJnAcHIBG1U2PGjEH//v2xdu1aAMDhw4cxc+ZMHD58GJcuXUJ+fj7Wr1+P9PR0TJkyBUDdwKNHH30Uo0aNQkREhPEzcuRIREdHNzoQCQBWr16NSZMmYciQIQCAESNG4IMPPsA333yDTZs2YcSIEa130UQOwjAlaseSkpKwY8cOFBUVoV+/fvDw8MDvf/97DB48GD/72c/wP//zP9i5cyfi4+Oh1WqxZ88eTJ8+3eyxpk+fjj179kCr1Vo83/nz5/H3v/8dq1atMm578sknMWnSJMTExOCbb77BX/7ylxa/TiJHY58pERGRRLwzJSIikohhSkREJBHDlIiISCKGKRERkUQMUyIiIokYpkRERBIxTImIiCRimBIREUnEMCUiIpKIYUpERCQRw5SIiEii/w+cXZyB9lsdswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (5, 5), dpi = 100)\n",
    "RSA_bns = np.linspace(0, 100, 40)\n",
    "a = func_df.query('RSA == RSA').RSA\n",
    "b = unk_df.query('RSA == RSA').RSA\n",
    "print(len(a), len(b))\n",
    "sns.histplot(a, label = \"Functional\", fill = True, color = \"#0C7BDC\", stat = \"proportion\", element = \"step\", bins = RSA_bns)\n",
    "sns.histplot(b, label = \"Unknown\", fill = True, color = \"#FFC20A\", stat = \"proportion\", element = \"step\", bins = RSA_bns)\n",
    "plt.ylabel(\"p\")\n",
    "plt.xlabel(\"RSA %\")\n",
    "plt.legend(edgecolor = \"k\")\n",
    "#plt.savefig(os.path.join(figs_dir, \"FUNC_vs_UNK_rsa.png\"))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/1211748.1.all.q/ipykernel_24921/3228072165.py:1: UserWarning: p-value floored: true value smaller than 0.001\n",
      "  scipy.stats.anderson_ksamp([a, b])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Anderson_ksampResult(statistic=42.72297518369387, critical_values=array([0.325, 1.226, 1.961, 2.718, 3.752, 4.592, 6.546]), pvalue=0.001)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.anderson_ksamp([a, b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RSA average (ALL residues) =  25.4\n",
      "RSA average (BS residues) =  34.1\n",
      "Kolmogorov-Smirnov two sample test result:\n",
      "KS statistic = 0.17\tp = 0.0\n",
      "\n",
      "Mann-Whitney U rank test result:\n",
      "U = 561780.5\tp = 0.0\n",
      "\n",
      "T-test result:\n",
      "T = -7.85\tp = 0.0\n"
     ]
    }
   ],
   "source": [
    "print(\"RSA average (ALL residues) = \", round(statistics.mean(a),1))\n",
    "print(\"RSA average (BS residues) = \", round(statistics.mean(b),1))\n",
    "\n",
    "s, p = scipy.stats.ks_2samp(a, b, alternative=\"two-sided\")\n",
    "print(\"Kolmogorov-Smirnov two sample test result:\\nKS statistic = {}\\tp = {}\\n\".format(round(s,2), round(p,2)))\n",
    "\n",
    "s, p = scipy.stats.mannwhitneyu(a, b, alternative=\"two-sided\")\n",
    "print(\"Mann-Whitney U rank test result:\\nU = {}\\tp = {}\\n\".format(round(s,2), round(p,2)))\n",
    "\n",
    "s, p = scipy.stats.ttest_ind(a, b)\n",
    "print(\"T-test result:\\nT = {}\\tp = {}\".format(round(s,2), round(p,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FUNCTIONAL VS UNKNOWN: DIVERGENCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "468 1197\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAHACAYAAAD5vIKYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/RUlEQVR4nO3de1xVVeL///fhcBdFQ+WiqGipGHkJ5tMHE7UZBy9d1I+VmWnOaPMjLUW+TWnmx7Qpq/Hjh7G8jGaZXdSZsaz8MCWaMpaUqViWZmkqZjCGjaKh3M7+/UGc6XjOQWADBziv5+NxHsXaa++z9vbydu299loWwzAMAQCAWvPxdAMAAGjqCFMAAEwiTAEAMIkwBQDAJMIUAACTCFMAAEwiTAEAMIkwBQDAJF9PN6Axstls+u6779SyZUtZLBZPNwcA4CGGYej8+fOKioqSj4/7/idh6sJ3332n6OhoTzcDANBInDx5Uh07dnS7nTB1oWXLlpIqLl6rVq083BoAgKcUFhYqOjrangvuEKYuVN7abdWqFWEKALjiIz8GIAEAYBJhCgCASYQpAAAm8cwUAK7AZrOpoKBAZ8+elc1m83RzUEd8fHzUunVrtW3btsrXXqqDMAWAKpw6dUrz58/Xvn37PN0U1JP4+Hj993//tzp06FDrYxCmAOBGSUmJxo0bp9atW+vJJ59Ux44dZbVaPd0s1JHy8nJ9++23Wrp0qcaNG6etW7fK39+/VsciTAHAjRMnTqioqEhLlixR3759Pd0c1INrr71W4eHhmjJlinJzc3X11VfX6jgMQAIAN8rLyyVJgYGBHm4J6lPlr29ZWVmtj0GYAgBgEmEKAIBJPDMFgFo4ebZUZ4oa7jWZsGAfRbf2a7Dvqw9dunRRamqqUlNTG+w7jx8/rpiYGOXk5NTrc2/CFABq6OTZUl2fflKXyowG+85AX4v2pUZXO1AnTZqkl19+2an866+/rvUgm+pas2aNUlNTdfbsWYfyTz75RC1atKjX7/YUwhQAauhMkU2Xygz1aOenYL/6X/O4qNTQ4e8resLRrau/37Bhw/TSSy85lLVr165uG1cDnvzu+sYzUwCopWA/i0ICfOr9U9vADggIUEREhMNn8uTJGjVqlEO91NRUDR482P7z4MGDNX36dD388MO66qqrFBERoccff9xhn7Nnz+p3v/udwsPDFRgYqLi4OG3evFk7duzQb37zG507d04Wi0UWi8W+b5cuXZSenm4/Rm5urkaOHKmQkBC1atVKd955p/75z3/atz/++OPq27evXnnlFXXp0kWhoaG66667dP78eXudd999VwMGDFDr1q0VFhamW265RUePHq3V9TKDMK1nRwpKtP+7YpefIwUlnm4eALj08ssvq0WLFvr444/17LPPasGCBcrMzJRUMb3i8OHDtWvXLr366qs6ePCgnn76aVmtVvXv31/p6elq1aqV8vLylJeXp4ceesjp+IZhaNSoUfrhhx+UlZWlzMxMHT16VGPHjnWod/ToUW3atEmbN2/W5s2blZWVpaefftq+/ccff1RaWpo++eQTbdu2TT4+Pho9enSDT/vIbd56dKSgRP3ST1ZZJyc1Wle3rd2MGwBQlc2bNyskJMT+8/Dhw6v9zLJ3796aN2+eJOmaa67R888/r23btunXv/61tm7dqt27d+vQoUPq3r27JKlr1672fUNDQ2WxWBQREeH2+Fu3btVnn32mY8eOKTo6WpL0yiuv6Nprr9Unn3yiX/ziF5IqgnvNmjX2xbknTJigbdu26cknn5QkjRkzxuG4q1evVvv27XXw4EHFxcVV61zrAj3TenShpGJwQo92fuoX5e/w6dHOz6EOANS1m266Sfv377d/lixZUu19e/fu7fBzZGSkTp8+LUnav3+/OnbsaA/S2jh06JCio6PtQSpJvXr1UuvWrXXo0CF7WZcuXexBenk7pIqe6913362uXbuqVatWiomJkVRxC7kh0TNtAJXPVRyx8gSA+tWiRQunkbs+Pj4yDMd/xJeWljrt6+fnOGrYYrHYb50GBQWZbpthGLJYnJ8FX15eVTsk6dZbb1V0dLRWrVqlqKgo2Ww2xcXFqaSkYR+j0TMFAC/Srl075eXlOZTt37+/Rsfo3bu3vv32W3311Vcut/v7+9unYnSnV69eys3N1cmT/34UdvDgQZ07d06xsbHVaseZM2d06NAhPfbYY/rVr36l2NhY/etf/6r+idQheqYAUEtFpYYa4i5TxffUjV/+8pf64x//qLVr1yoxMVGvvvqqPv/8c/Xr16/axxg0aJAGDhyoMWPGaPHixbr66qv15ZdfymKxaNiwYerSpYsuXLigbdu2qU+fPgoODlZwcLDDMYYMGaLevXtr/PjxSk9PV1lZmaZOnapBgwYpISGhWu1o06aNwsLCtHLlSkVGRio3N1ezZs2q0fWoK4QpANRQWLCPAn0tOvy98+3R+hLoa1FYsPmbiUOHDtXcuXP18MMP69KlS/rtb3+riRMn6sCBAzU6zsaNG/XQQw9p3Lhx+vHHH3X11VfbR9n2799fKSkpGjt2rM6cOaN58+Y5vVpjsVi0adMmPfjggxo4cKB8fHw0bNgwPffcc9Vug4+Pj9avX6/p06crLi5OPXr00JIlSxxe82koFuPym+dQYWGhQkNDde7cObVq1arWx9n/XbGSln2rflH+Ts9MLxTblPNdiXZO7ai+UQFmmwygHnz55Ze655579Oqrr6pnz54O25hOsPmo6te5unlAzxQAaiG6tV+NZiNC88YAJAAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAk3jMFgNq4dFIqK2i47/NtKwVGX7leHRs8eLD69u3rsKg3nBGmAFBTl05Ke/tKtksN950+gVL8/moHqrsQ3LRpk0aPHu20cgzMIUwBoKbKCiqCNKiHZA2+cn2zyouki4d/6gk3fO8UV8YzUwCoLWuwZG3ZAJ/6CezHH39cffv21SuvvKIuXbooNDRUd911l86fP+92n3fffVehoaFau3atJGnSpEkaNWqUFi1apMjISIWFhWnatGkOa6T+61//0sSJE9WmTRsFBwdr+PDh+vrrryVVrF/arl07bdy40V6/b9++at++vf3n7Oxs+fn56cKFC5IqJsl/4YUXNHr0aAUHB+uaa67R22+/XafXpqYIUwDwYkePHtWmTZu0efNmbd68WVlZWfbVXy63fv163XnnnVq7dq0mTpxoL9++fbuOHj2q7du36+WXX9aaNWu0Zs0a+/ZJkyZpz549evvtt5WdnS3DMDRixAiVlpbKYrFo4MCB2rFjh6SK4D148KBKS0t18OBBSdKOHTsUHx+vkJAQ+zHnz5+vO++8U5999plGjBih8ePH64cffqj7C1RNhCkAeDGbzaY1a9YoLi5OSUlJmjBhgrZt2+ZUb9myZUpJSdFbb72lkSNHOmxr06aNnn/+efXs2VO33HKLbr75Zvsxvv76a7399tt64YUXlJSUpD59+ui1117TqVOntGnTJkkVz3crw/Qf//iH+vTpo1/+8pf2sh07djgtqzZp0iSNGzdOV199tZ566in9+OOP2r17d51em5ogTAHAi3Xp0kUtW7a0/xwZGanTp0871Nm4caNSU1O1ZcsW3XTTTU7HuPbaa2W1Wl0e49ChQ/L19dUNN9xg3x4WFqYePXro0KFDkirC9IsvvlBBQYGysrI0ePBgDR48WFlZWSorK9OuXbs0aNAgh+/s3bu3/f9btGihli1bOrW7IRGmANAMtWrVSufOnXMqP3v2rMO6nH5+jmukWiwW2WyO67T27dtX7dq100svveRyFHBVx3A3atgwDFksFklSXFycwsLClJWVZQ/TQYMGKSsrS5988okuXryoAQMGVPs7PYEwBYBmqGfPntqzZ49T+SeffKIePXrU6FjdunXT9u3b9dZbb+nBBx+s0b69evVSWVmZPv74Y3vZmTNn9NVXXyk2NlaS7M9N33rrLX3++edKSkrSddddp9LSUq1YsULXX3+9Q++5MSJMAaAZmjp1qo4ePapp06bp008/1VdffaWlS5dq9erV+v3vf1/j43Xv3l3bt2+33/KtrmuuuUYjR47Ufffdpw8++ECffvqp7rnnHnXo0MHh2evgwYP1+uuvq3fv3mrVqpU9YF977TWn56WNEe+ZAkBtlRc12u/p0qWLdu7cqTlz5ig5OVmXLl1S9+7dtWbNGt1xxx21akaPHj30/vvva/DgwbJarfqf//mfau330ksvacaMGbrllltUUlKigQMHKiMjw+FW7U033aTy8nKH4Bw0aJA2bdrk9Ly0MbIYTIPhpLCwUKGhoTp37pzDs4Wa2v9dsZKWfat+Uf4KCXC8CXCh2Kac70q0c2pH9Y0KMNtkAPXgyy+/1D333KNXX31VPXv2/PeGJjADEqrP7a+zqp8H9EwBoKYCoyuCzQvm5kX1EKYAUBuB0WJqP1RiABIAACZ5PEyXLVummJgYBQYGKj4+Xjt37nRbNy8vT3fffbd69OghHx8flyPKVq1apaSkJLVp00Zt2rTRkCFDPDorBgCg+fNomG7YsEGpqamaM2eOcnJylJSUpOHDhys3N9dl/eLiYrVr105z5sxRnz59XNbZsWOHxo0bp+3btys7O1udOnVScnKyTp06VZ+nAgDwYh4N08WLF2vy5MmaMmWKYmNjlZ6erujoaC1fvtxl/S5duuhPf/qTJk6cqNDQUJd1XnvtNU2dOlV9+/ZVz549tWrVKtlsNpdzTQJAVXx8Kv6K/PkKKGh+Kn99K3+9a8NjYVpSUqK9e/cqOTnZoTw5OVm7du2qs+8pKipSaWmprrrqqjo7JgDvEBERIUnKycnxcEtQnyp/fSMjI2t9DI+N5i0oKFB5ebnCw8MdysPDw5Wfn19n3zNr1ix16NBBQ4YMcVunuLhYxcXF9p8LCwvr7PsBNF2tWrXS6NGj9dxzz0mS+vXr5zQnLJqu0tJS5eTk6LnnntPo0aNNTVno8VdjKic6rvTzyY/NevbZZ7Vu3Trt2LFDgYGBbustXLhQ8+fPr5PvBNC8zJ49W5K0ZMkSD7cE9WX06NH2X+fa8liYtm3bVlar1akXevr0aafeam0sWrRITz31lLZu3eqwVI8rs2fPVlpamv3nwsJCRUfz/hiAiudoc+bM0fTp05WXl+fRlUlQt3x8fBQZGVknk+h7LEz9/f0VHx+vzMxMjR492l6emZnptPBsTf3xj3/UH/7wB7333ntKSEi4Yv2AgAAFBDClHwD3WrZs2ehXLoHnePQ2b1pamiZMmKCEhAQlJiZq5cqVys3NVUpKiqSKHuOpU6e0du1a+z779++XJF24cEHff/+99u/fL39/f/Xq1UtSxa3duXPn6vXXX1eXLl3sPd+QkBCFhIQ07AkCALyCR8N07NixOnPmjBYsWKC8vDzFxcUpIyNDnTt3llQxScPl75z269fP/v979+7V66+/rs6dO+v48eOSKiaBKCkp0e233+6w37x58/T444/X6/kAALyTxwcgTZ06VVOnTnW5bc2aNU5lV1rkpjJUAQBoKB6fThAAgKaOMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADDJ42G6bNkyxcTEKDAwUPHx8dq5c6fbunl5ebr77rvVo0cP+fj4KDU11WW9jRs3qlevXgoICFCvXr305ptv1lPrAQDwcJhu2LBBqampmjNnjnJycpSUlKThw4crNzfXZf3i4mK1a9dOc+bMUZ8+fVzWyc7O1tixYzVhwgR9+umnmjBhgu688059/PHH9XkqAAAvZjEMw/DUl99www26/vrrtXz5cntZbGysRo0apYULF1a57+DBg9W3b1+lp6c7lI8dO1aFhYX6+9//bi8bNmyY2rRpo3Xr1lWrXYWFhQoNDdW5c+fUqlWr6p/QZfZ/V6ykZd+qX5S/QgIc/91yodimnO9KtHNqR/WNCqj1dwAA6k9188BjPdOSkhLt3btXycnJDuXJycnatWtXrY+bnZ3tdMyhQ4dWeczi4mIVFhY6fAAAqC6PhWlBQYHKy8sVHh7uUB4eHq78/PxaHzc/P7/Gx1y4cKFCQ0Ptn+jo6Fp/PwDA+3h8AJLFYnH42TAMp7L6Pubs2bN17tw5++fkyZOmvh8A4F18PfXFbdu2ldVqdeoxnj592qlnWRMRERE1PmZAQIACAnhuCQCoHY/1TP39/RUfH6/MzEyH8szMTPXv37/Wx01MTHQ65pYtW0wdEwCAqnisZypJaWlpmjBhghISEpSYmKiVK1cqNzdXKSkpkipuv546dUpr166177N//35J0oULF/T9999r//798vf3V69evSRJM2bM0MCBA/XMM89o5MiReuutt7R161Z98MEHDX5+AADv4NEwHTt2rM6cOaMFCxYoLy9PcXFxysjIUOfOnSVVTNJw+Tun/fr1s///3r179frrr6tz5846fvy4JKl///5av369HnvsMc2dO1fdunXThg0bdMMNNzTYeQEAvItH3zNtrHjPFAAgNYH3TAEAaC4IUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkj4fpsmXLFBMTo8DAQMXHx2vnzp1V1s/KylJ8fLwCAwPVtWtXrVixwqlOenq6evTooaCgIEVHR2vmzJm6dOlSfZ0CAMDLeTRMN2zYoNTUVM2ZM0c5OTlKSkrS8OHDlZub67L+sWPHNGLECCUlJSknJ0ePPvqopk+fro0bN9rrvPbaa5o1a5bmzZunQ4cOafXq1dqwYYNmz57dUKcFAPAyvp788sWLF2vy5MmaMmWKpIoe5Xvvvafly5dr4cKFTvVXrFihTp06KT09XZIUGxurPXv2aNGiRRozZowkKTs7WzfeeKPuvvtuSVKXLl00btw47d69u2FOCgDgdTzWMy0pKdHevXuVnJzsUJ6cnKxdu3a53Cc7O9up/tChQ7Vnzx6VlpZKkgYMGKC9e/faw/Obb75RRkaGbr75ZrdtKS4uVmFhocMHAIDq8ljPtKCgQOXl5QoPD3coDw8PV35+vst98vPzXdYvKytTQUGBIiMjddddd+n777/XgAEDZBiGysrKdP/992vWrFlu27Jw4ULNnz/f/EkBALySxwcgWSwWh58Nw3Aqu1L9n5fv2LFDTz75pJYtW6Z9+/bpjTfe0ObNm/XEE0+4Pebs2bN17tw5++fkyZO1PR0AgBfyWM+0bdu2slqtTr3Q06dPO/U+K0VERLis7+vrq7CwMEnS3LlzNWHCBPtz2Ouuu04//vijfve732nOnDny8XH+90NAQIACAgLq4rTqzJGCEl0oMVxuC/G36Oq2/g3cIgCAOx4LU39/f8XHxyszM1OjR4+2l2dmZmrkyJEu90lMTNQ777zjULZlyxYlJCTIz89PklRUVOQUmFarVYZh2Huxjd2RghL1S6+6d5yTGk2gAkAj4dHRvGlpaZowYYISEhKUmJiolStXKjc3VykpKZIqbr+eOnVKa9eulSSlpKTo+eefV1pamu677z5lZ2dr9erVWrdunf2Yt956qxYvXqx+/frphhtu0JEjRzR37lzddtttslqtHjnPmqrskfZo56dgP8fb2kWlhg5/X+q21woAaHgeDdOxY8fqzJkzWrBggfLy8hQXF6eMjAx17txZkpSXl+fwzmlMTIwyMjI0c+ZMLV26VFFRUVqyZIn9tRhJeuyxx2SxWPTYY4/p1KlTateunW699VY9+eSTDX5+ZgX7WRQScPltaZtH2gIAcM9iNJV7nw2osLBQoaGhOnfunFq1alXr4+z/rlhJy75Vvyh/p1C8UGxTzncl2jm1o/pGBdTJfgCAulXdPPD4aF4AAJo6whQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwyXSYGoYhwzDqoi0AADRJtQ7T1atXKy4uToGBgQoMDFRcXJxeeOGFumwbAABNgm9tdpo7d67+93//Vw8++KASExMlSdnZ2Zo5c6aOHz+uP/zhD3XaSAAAGrNaheny5cu1atUqjRs3zl522223qXfv3nrwwQcJUwCAV6lVmJaXlyshIcGpPD4+XmVlZaYb5U0Ony6pVhkAoPGqVZjec889Wr58uRYvXuxQvnLlSo0fP75OGtbcWX96Wj3lb6evWAcA0LjVKkyligFIW7Zs0X/+539Kkj766COdPHlSEydOVFpamr3e5YGLCkF+Pkro6K9ym+vtVp+KOgCAxq9WYfr555/r+uuvlyQdPXpUktSuXTu1a9dOn3/+ub2exWKpgyY2X4QlADQPtQrT7du313U7AABosugaAQBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACZ5PEyXLVummJgYBQYGKj4+Xjt37qyyflZWluLj4xUYGKiuXbtqxYoVTnXOnj2radOmKTIyUoGBgYqNjVVGRkZ9nQIAwMt5NEw3bNig1NRUzZkzRzk5OUpKStLw4cOVm5vrsv6xY8c0YsQIJSUlKScnR48++qimT5+ujRs32uuUlJTo17/+tY4fP66//e1vOnz4sFatWqUOHTo01GkBALyMrye/fPHixZo8ebKmTJkiSUpPT9d7772n5cuXa+HChU71V6xYoU6dOik9PV2SFBsbqz179mjRokUaM2aMJOnFF1/UDz/8oF27dsnPz0+S1Llz54Y5IQCAV/JYz7SkpER79+5VcnKyQ3lycrJ27drlcp/s7Gyn+kOHDtWePXtUWloqSXr77beVmJioadOmKTw8XHFxcXrqqadUXl5ePycCAPB6HuuZFhQUqLy8XOHh4Q7l4eHhys/Pd7lPfn6+y/plZWUqKChQZGSkvvnmG73//vsaP368MjIy9PXXX2vatGkqKyvTf//3f7s8bnFxsYqLi+0/FxYWmjw7AIA38fgAJIvF4vCzYRhOZVeq//Nym82m9u3ba+XKlYqPj9ddd92lOXPmaPny5W6PuXDhQoWGhto/0dHRtT0dAIAX8liYtm3bVlar1akXevr0aafeZ6WIiAiX9X19fRUWFiZJioyMVPfu3WW1Wu11YmNjlZ+fr5KSEpfHnT17ts6dO2f/nDx50sypAQC8jMfC1N/fX/Hx8crMzHQoz8zMVP/+/V3uk5iY6FR/y5YtSkhIsA82uvHGG3XkyBHZbDZ7na+++kqRkZHy9/d3edyAgAC1atXK4QMAQHV59DZvWlqaXnjhBb344os6dOiQZs6cqdzcXKWkpEiq6DFOnDjRXj8lJUUnTpxQWlqaDh06pBdffFGrV6/WQw89ZK9z//3368yZM5oxY4a++uor/d///Z+eeuopTZs2rcHPDwDgHTz6aszYsWN15swZLViwQHl5eYqLi1NGRob9VZa8vDyHd05jYmKUkZGhmTNnaunSpYqKitKSJUvsr8VIUnR0tLZs2aKZM2eqd+/e6tChg2bMmKFHHnmkwc8PAOAdLEblCB7YFRYWKjQ0VOfOnTN1y3f/d8VKWvat+kX5KySgbm4CXCi2Kee7Eu2c2lF9owLq5JiVjhSU6EKJ698OIf4WXd3W9W1yAGiuqpsHHu2ZovE4UlCifulVD7zKSY0mUAHABcIUkmTvkfZo56dgP8fXj4pKDR3+vtRtrxUAvB1hCgfBfhYXt6RtLusCACp4fNIGAACaOsIUAACTCFMAAEwiTAEAMIkwBQDAJMIUAACTCFMAAEwiTAEAMIkwBQDAJMIUAACTCFMAAEwiTAEAMIkwBQDAJMIUAACTWIKtiTp8usTtthB/C4t4A0ADIkybGOtP9xKm/O10lfVyUqMJVABoIIRpExPk56OEjv4qd7Ned1GpocPfl+pCidGwDQMAL0aYNkFBflU96naTsgCAesMAJAAATCJMAQAwidu8qDZ3I4gZPQzA2xGmuKLqjCBm9DAAb0aY4oqqGkHM6GEAIExRTe5HEDN6GAAYgAQAgEmEKQAAJnGbF/XuSEGJ22eqjAQG0BwQpqhXRwpK1C/9ZJV1GAkMoKkjTFGvKnukPdr5KdjP4rCNkcAAmgvCFA0i2M+ikIDLH9EzEhhA88AAJAAATKJnCtS3i0ek8vOut1lbSkFXN2x7ANQ5wrSZYh7dRuLiEWlPn6rrJHxKoAJNHGHazDCPrqruCUoN2xusbEdQD8kafNm2Iuni4arbCqBJIEybGa+fR7c6PUGp4XuD1uCKEAfQLBGmzZBXz6NbVU9QojcIoF4Qpmie6AkCaECEKeqEuwFP7sqrU4fBUgCaCsIUplRnwNPP69V03yYzWMrdoKeiww3fFgANjjCFKVUNeKpk9XH9HLfZDJaq1qAn/qgBzRl/wmGa+wFPZva1qVvwCQUVn5Eu+DlvbkyTHVxp0JN8JWuQ+/2r6r02pvME4BZhikapo/9xvZ50q/StKj6uNLbJDmo86OmnP36HJ1ddrbGdJwAnHp+bd9myZYqJiVFgYKDi4+O1c+fOKutnZWUpPj5egYGB6tq1q1asWOG27vr162WxWDRq1Kg6bjXqW7DPj5Kki37dpZB+jp+gHhWVmvrrLdYgKeQXzufX3M4T8AIeDdMNGzYoNTVVc+bMUU5OjpKSkjR8+HDl5ua6rH/s2DGNGDFCSUlJysnJ0aOPPqrp06dr48aNTnVPnDihhx56SElJSfV9GqhHNstPvT2Hj6tbqU2UNcjF+TXD8wSaOY+G6eLFizV58mRNmTJFsbGxSk9PV3R0tJYvX+6y/ooVK9SpUyelp6crNjZWU6ZM0W9/+1stWrTIoV55ebnGjx+v+fPnq2vXrg1xKgAAL+axMC0pKdHevXuVnJzsUJ6cnKxdu3a53Cc7O9up/tChQ7Vnzx6VlpbayxYsWKB27dpp8uQrPIv6SXFxsQoLCx0+AABUl8fCtKCgQOXl5QoPD3coDw8PV35+vst98vPzXdYvKytTQUGBJOnDDz/U6tWrtWrVqmq3ZeHChQoNDbV/oqOja3g2AABv5vEBSBaLxeFnwzCcyq5Uv7L8/Pnzuueee7Rq1Sq1bdu22m2YPXu2zp07Z/+cPHmyBmcAAPB2Hns1pm3btrJarU690NOnTzv1PitFRES4rO/r66uwsDB98cUXOn78uG699Vb7dputYkYAX19fHT58WN26dXM6bkBAgAICAsyeEgDAS3ksTP39/RUfH6/MzEyNHj3aXp6ZmamRI0e63CcxMVHvvPOOQ9mWLVuUkJAgPz8/9ezZUwcOHHDY/thjj+n8+fP605/+xO1bmHKkoMTljExBxaXqUZ9f7G5SByZ0ABoNj07akJaWpgkTJighIUGJiYlauXKlcnNzlZKSIqni9uupU6e0du1aSVJKSoqef/55paWl6b777lN2drZWr16tdevWSZICAwMVFxfn8B2tW7eWJKdyoCaOFJSoX7rr2/99Wv5TH/SXLpYZCrLW5bdWY1IHJnQAGgWPhunYsWN15swZLViwQHl5eYqLi1NGRoY6d+4sScrLy3N45zQmJkYZGRmaOXOmli5dqqioKC1ZskRjxozx1CnAS1T2SHu081Own+Nz+2jfij9GtrpeLrZyUgeVOW9jXVagUfH4dIJTp07V1KlTXW5bs2aNU9mgQYO0b9++ah/f1TGA2gr2sygkwHHcXlAVA+ZMq2pOXwCNhsdH8wIA0NR5vGfa3HULPqHugSUK8nf+d8tFWwt9V9qlwdvjaiWWoOJSdQu+IOmaBm0PADQHhGk98i85qv1Jt1ZZZ8qxLQ0WqFWtxNJD0v4k6e7D7+qsmIIRAGqCMK1HVuOCJOnYpa4qNRwnLQ/wuajOAUcV9NPqKA3h5yuxBAW0cNh2sfhHBZV+pWCfH3W2hutxR/kdr/I8PNEDB4CGRJg2gEu2IJXIMbxU1yM/a8C+EotDWQ0T9CdRfsf1QkzyFes1ZA8cABoaYQpTKnukJ4q7qdjmPPLUEz1wAGhohCnqRLEtSBeNFs4b6rMH7mpmIHezBQFAPSJM0QRVY2YgfmsDaED8jYOmp6qZgSRJvkx2AKBBEaZomghLAI0IYQo0ZXW9oszFI1XP98tKNYBLhCnQJNXDijIXj0h7+ly5HivVAE4IU6Apqo8VZSrrB/WQrMEutrNSDeAOYQo0VfX13NjqPKkHgKqxagwAACbRM22kqprvlrlu4VF1PegJaAYI00aoOvPdMtctGl49DHoCmgnCtBGqar5b5rqFx9THoCegmSBMGzGX8916cLUZgMkyANcIUzjoHPiNgmzO49Ka3XNab37u525iBhYJAGqNMIUkyZBVkjQ3+vdu6zSP57RN67lfUalN5W7uRlh9pGC/Gg7Ir9bEDPy1ANQUf2ogSTJ8grT77HVq6W+T1WJx2NacntMW2QL0+dnr5Gspd9oWbL2kniHf6GTBvxQd7YHGXaao1KY93xZXWSehY0DNAvVKEzOwSABQK4Qp7C7ZAuVrs8jq4ximlc9po/2POu3jqswVd/Ua+vZxua3iPIP9LHI6zZ/+EVFUZrjct1vwCXUPLFGQv2N4tff5pt7aKsl1Ww2pqNRw22u9IiZmAOoUYYorKv/pFvDDkQ9dsU5t9vXE7WMfi5z+0WBxU1eS/EuOan/SrVUe03BzDcxy1VbZXAc+AM8gTHFFJUagDl3sI6ucb41KFYFZYgTWeN/Gevs4sOQr6YKfQ1mL4i8kSccudVWp4Xh7tNwwdL7ER7ER3B4FvBVhimpxF5am9m1kr/nYjIqeZefT/5902nFb55/+W1QeojKLY2iW2wxdoqcIeDXC1AtdKrPJVuyYZEWlnks2V89TK59DumqrVMuRrFdQbARq99nrdG17q1r4O95W/bHE0Beny+XrGyhrVfeDAXglwtSL+PyUPSf+VaYL5a5HiVoaMCiq8zz16A+GLtlct7XGI1mr4ZItUDafgIq0/hmbj02XbMUKqdNvA9BcEKZeJMC3IiCC/Cxy1b2yWCSfBkzTKz2LLTN85Osb6BRgpkeyAkAdI0y9kNXi4vUXD7nSs1irq44nzyerx9WMRsxyBNQLwhRodqoxyxN/9IE6xZ+oJqqqyRKa3Ty6qJmqVneRVK+zHHnznMfwaoRpE1OdQTtSc5lHt2ruRiDXx0hfM+p8ft3qaPApAZvWnMdAXSNMm5grDdpprBMh1KXKMVKHvy91W6c+RvrWRr3Mr9sYsdYpvBxh2gRVOWjHC0a4+lgsahkgGS7GITW2kb71Or9uY8ME+fBihCmaJB+LxfVkuo10pC/z6wLNWxO/twQAgOcRpgAAmMRt3mbKzNqjAICaIUybGTNrjwIAaocwbWbMrD0KAKgdwrQZIiwBoGExAAkAAJMIUwAATPJ4mC5btkwxMTEKDAxUfHy8du7cWWX9rKwsxcfHKzAwUF27dtWKFSsctq9atUpJSUlq06aN2rRpoyFDhmj37t31eQpohIpKbTpf7PxxN58vAJjh0TDdsGGDUlNTNWfOHOXk5CgpKUnDhw9Xbm6uy/rHjh3TiBEjlJSUpJycHD366KOaPn26Nm7caK+zY8cOjRs3Ttu3b1d2drY6deqk5ORknTp1qqFOCx7083l7c74rdvpUzufbgGugA/ACHh2AtHjxYk2ePFlTpkyRJKWnp+u9997T8uXLtXDhQqf6K1asUKdOnZSeni5Jio2N1Z49e7Ro0SKNGTNGkvTaa6857LNq1Sr97W9/07Zt2zRx4sT6PSF4XFXz9layWH6ajhAA6ojHwrSkpER79+7VrFmzHMqTk5O1a9cul/tkZ2crOTnZoWzo0KFavXq1SktL5efn57RPUVGRSktLddVVV7ltS3FxsYqL/72yR2FhYU1OxRQmV6h7bufthWex1imaMY+FaUFBgcrLyxUeHu5QHh4ervz8fJf75Ofnu6xfVlamgoICRUZGOu0za9YsdejQQUOGDHHbloULF2r+/Pm1OIvaY3IFeA/WOkXz5/H3TC2X3W4zDMOp7Er1XZVL0rPPPqt169Zpx44dCgx0/+7l7NmzlZaWZv+5sLBQ0dHR1Wp/bTG5Aiq5GxRVncFSTWKBdNY6hRfwWJi2bdtWVqvVqRd6+vRpp95npYiICJf1fX19FRYW5lC+aNEiPfXUU9q6dat69+5dZVsCAgIUEBBQi7Mwh7D0btVZ5Pzn9Wq6b6NadJy1TtHMeexPmr+/v+Lj45WZmelQnpmZqf79+7vcJzEx0an+li1blJCQ4PC89I9//KOeeOIJvfvuu0pISKj7xgN1oGKwlEUh/u4/LQMsLgdLVbVvsF9F/Waz6DjQBHj0Nm9aWpomTJighIQEJSYmauXKlcrNzVVKSoqkituvp06d0tq1ayVJKSkpev7555WWlqb77rtP2dnZWr16tdatW2c/5rPPPqu5c+fq9ddfV5cuXew92ZCQEIWEhDT8SQJVMDNYqqktkF4lBiehifNomI4dO1ZnzpzRggULlJeXp7i4OGVkZKhz586SpLy8PId3TmNiYpSRkaGZM2dq6dKlioqK0pIlS+yvxUgVk0CUlJTo9ttvd/iuefPm6fHHH2+Q8wJQXQxOQvPg8QFIU6dO1dSpU11uW7NmjVPZoEGDtG/fPrfHO378eB21DEC9Y3ASmgmPhynQ2LgaIest0xAWldqqfNZaL6OEGZyEZoAwBX5SnRGyzXnipKJSm/Z8W3zFeo1qlDDQSBCmwE+uNBVhc5+GsLJHGuxnkY+L07QZUlGpwShhwAXCFPgZpiKUfCyS1V2aAnCJMAXQdF084n6AEq/VoAERpgCapotHpD19qq7DazVoIIQpUEe8eRSwR1T2SIN6SNbgy7bxWg0aFmEKmOTto4A9zhpccUsX8CDCFDDJ20cBAyBMgTrBKGDAu/HmNQAAJhGmAACYxG1eADXiboRyVfP2VjXnb73M9ws0MMIUQLVUZ9Syq3l7qzPnL/P9oqkjTAFUS1Wjlquat7eqOX+Z7xfNBWEKoNrcjlquxry9Luf8Zb5fNBPcVwEAwCTCFAAAk7jNC3gZdyNrm+U8wkWHXZfXdkWZqlapMXNcNHmEKeBFqjOytnnMfPjTX22HJ7uvUtMVZaqzSk1tjotmgTAFvEhVI2ulRjqPsLvepbtySbIGSSG/kFTmvK22K8pUtUqNmeOiWSBMAS/kcmRto1ON3uXP613OGlSnrfn3cVmlBs4IUwCNU1W9Szvf+gvNulbV81aetTZ5hCmAxqupBOWVVOd5K89amzTCFGimXI3ObZYjdpuCqp638qy1WSBMgWamOnPoNrYxRl6D563NFmEKNDNVzaErNdIRu0ATR5gCzZDbOXQB1AumEwQAwCR6pgDwc+5eYalqkgh4PcIUACpVa8pA/tqEM35XAEClK00Z2JQmiUCDIkwB4HK8woIaYgASAAAmEaYAAJhEmAIAYBLPTAHUmdrOB+yujtVHCvZz/2/+olKbfY3W6vKxGWqhiukWLxY6LpQeVFyqHjU7nItGuXiFhtdqmj3CFIBptZ0PuDr7JXQMcBmoRaU27fm22MUeVQuxluj6UGnyX/6pT89/67CtT8t/6oP+0sUyQ0HWmh65Ouuv8lduc8WvLADTajsfcFX72QypqNRw2/OsLA/2s6gm65wH/VR5aPRJ9bzk57Atwje34rtrs7jOFddf5bWa5owwBVAnajsfsNv9bG6S2Wl/yVqTNLVU/LU3N/r3bqsYqnG3tAJh6bUIUwBepcQI1KGLfWRVudO2csPQ+RIfxUYQiqgZwhSA1ykxAl2Wl9sMXapmj7jOuRukZG0pBV3dsG1BjRGmAHCZ6oxAduVKo49dq8bApYRPCdRGjjAFgJ9UZ3TxlbgbfexWVQOXyouki4ddr2IjuV/hRqJH28A8PmnDsmXLFBMTo8DAQMXHx2vnzp1V1s/KylJ8fLwCAwPVtWtXrVixwqnOxo0b1atXLwUEBKhXr156880366v5AJqRitHFFoX41/wT7FeRxDV971VSRaBaW7r4uJps/yeVK9zkDHD92dOnog4ahEd7phs2bFBqaqqWLVumG2+8UX/+8581fPhwHTx4UJ06dXKqf+zYMY0YMUL33XefXn31VX344YeaOnWq2rVrpzFjxkiSsrOzNXbsWD3xxBMaPXq03nzzTd1555364IMPdMMNNzT0KQJoYmo7Krm6o49rpaqJIFytcFPZoz2/h55rA7EYhrs3w+rfDTfcoOuvv17Lly+3l8XGxmrUqFFauHChU/1HHnlEb7/9tg4dOmQvS0lJ0aeffqrs7GxJ0tixY1VYWKi///3v9jrDhg1TmzZttG7dumq1q7CwUKGhoTp37pxatWpV29PT4WO71ePbm3So6FqVKKTWxwG8UbnN0IUSQ/2iAtQywPkm2vlim3K+K1aIv6Vmr8bUkyu1t3YHvShd+KTqOiG/cH4lpzr7STyLrYbq5oHHeqYlJSXau3evZs2a5VCenJysXbt2udwnOztbycnJDmVDhw7V6tWrVVpaKj8/P2VnZ2vmzJlOddLT0+u0/QBQ72o7EcSV9rvSs1jUmMfCtKCgQOXl5QoPD3coDw8PV35+vst98vPzXdYvKytTQUGBIiMj3dZxd0xJKi4uVnHxv6clO3funKSKf5GYceH8BRX+aKjk4gWVGM7vtAFwz2ZIKjd08p8+CvRz7nleKjWkYptKyywq93zH9IrtrR/Fkn6s8V6+KlZrw9DBj15QoaVb3TerkQhu2Vm9e/7a1DEqc+BKN3E9PprXctkUY4ZhOJVdqf7l5TU95sKFCzV//nyn8ujoaPcNr5FDV64CAA1umacb0GScP39eoaGhbrd7LEzbtm0rq9Xq1GM8ffq0U8+yUkREhMv6vr6+CgsLq7KOu2NK0uzZs5WWlmb/2Waz6YcfflBYWFiVIXwlhYWFio6O1smTJ009e21uuC7ucW1c47q4x7Vxra6ui2EYOn/+vKKioqqs57Ew9ff3V3x8vDIzMzV69Gh7eWZmpkaOHOlyn8TERL3zzjsOZVu2bFFCQoL8/PzsdTIzMx2em27ZskX9+/d325aAgAAFBAQ4lLVu3bqmp+RWq1at+E3uAtfFPa6Na1wX97g2rtXFdamqR1rJo7d509LSNGHCBCUkJCgxMVErV65Ubm6uUlJSJFX0GE+dOqW1a9dKqhi5+/zzzystLU333XefsrOztXr1aodRujNmzNDAgQP1zDPPaOTIkXrrrbe0detWffDBBx45RwBA8+fRMB07dqzOnDmjBQsWKC8vT3FxccrIyFDnzp0lSXl5ecrNzbXXj4mJUUZGhmbOnKmlS5cqKipKS5Yssb9jKkn9+/fX+vXr9dhjj2nu3Lnq1q2bNmzYwDumAIB64/EBSFOnTtXUqVNdbluzZo1T2aBBg7Rv374qj3n77bfr9ttvr4vmmRIQEKB58+Y53UL2dlwX97g2rnFd3OPauNbQ18WjkzYAANAceHxuXgAAmjrCFAAAkwhTAABMIkwBADCJMK0nNV2ntblZuHChfvGLX6hly5Zq3769Ro0apcOHHZeRMgxDjz/+uKKiohQUFKTBgwfriy++8FCLPWfhwoWyWCxKTU21l3nrtTl16pTuuecehYWFKTg4WH379tXevXvt2731upSVlemxxx5TTEyMgoKC1LVrVy1YsEA2278XT/WGa/OPf/xDt956q6KiomSxWLRp0yaH7dW5BsXFxXrwwQfVtm1btWjRQrfddpu+/fZb840zUOfWr19v+Pn5GatWrTIOHjxozJgxw2jRooVx4sQJTzetwQwdOtR46aWXjM8//9zYv3+/cfPNNxudOnUyLly4YK/z9NNPGy1btjQ2btxoHDhwwBg7dqwRGRlpFBYWerDlDWv37t1Gly5djN69exszZsywl3vjtfnhhx+Mzp07G5MmTTI+/vhj49ixY8bWrVuNI0eO2Ot443UxDMP4wx/+YISFhRmbN282jh07Zvz1r381QkJCjPT0dHsdb7g2GRkZxpw5c4yNGzcakow333zTYXt1rkFKSorRoUMHIzMz09i3b59x0003GX369DHKyspMtY0wrQf/8R//YaSkpDiU9ezZ05g1a5aHWuR5p0+fNiQZWVlZhmEYhs1mMyIiIoynn37aXufSpUtGaGiosWLFCk81s0GdP3/euOaaa4zMzExj0KBB9jD11mvzyCOPGAMGDHC73Vuvi2EYxs0332z89re/dSj7r//6L+Oee+4xDMM7r83lYVqda3D27FnDz8/PWL9+vb3OqVOnDB8fH+Pdd9811R5u89axynVaL193tap1Wr1B5bJ2V111lSTp2LFjys/Pd7hOAQEBGjRokNdcp2nTpunmm2/WkCFDHMq99dq8/fbbSkhI0B133KH27durX79+WrVqlX27t14XSRowYIC2bdumr776SpL06aef6oMPPtCIESMkefe1qVSda7B3716VlpY61ImKilJcXJzp6+TxGZCam9qs09rcGYahtLQ0DRgwQHFxcZJkvxaurtOJEycavI0Nbf369dq3b58++eQTp23eem2++eYbLV++XGlpaXr00Ue1e/duTZ8+XQEBAZo4caLXXhdJeuSRR3Tu3Dn17NlTVqtV5eXlevLJJzVu3DhJ3vt75ueqcw3y8/Pl7++vNm3aONUx+/czYVpParqmanP2wAMP6LPPPnO52IA3XqeTJ09qxowZ2rJliwIDA93W87ZrY7PZlJCQoKeeekqS1K9fP33xxRdavny5Jk6caK/nbddFkjZs2KBXX31Vr7/+uq699lrt379fqampioqK0r333muv543X5nK1uQZ1cZ24zVvHarNOa3P24IMP6u2339b27dvVsWNHe3lERIQkeeV12rt3r06fPq34+Hj5+vrK19dXWVlZWrJkiXx9fe3n723XJjIyUr169XIoi42NtS924c2/Z37/+99r1qxZuuuuu3TddddpwoQJmjlzphYuXCjJu69Npepcg4iICJWUlOhf//qX2zq1RZjWsZ+v0/pzmZmZVa6p2twYhqEHHnhAb7zxht5//33FxMQ4bI+JiVFERITDdSopKVFWVlazv06/+tWvdODAAe3fv9/+SUhI0Pjx47V//3517drVK6/NjTfe6PT61FdffWVfRcqbf88UFRXJx8fxr2ur1Wp/Ncabr02l6lyD+Ph4+fn5OdTJy8vT559/bv46mRq+BJcqX41ZvXq1cfDgQSM1NdVo0aKFcfz4cU83rcHcf//9RmhoqLFjxw4jLy/P/ikqKrLXefrpp43Q0FDjjTfeMA4cOGCMGzeu2Q3lr66fj+Y1DO+8Nrt37zZ8fX2NJ5980vj666+N1157zQgODjZeffVVex1vvC6GYRj33nuv0aFDB/urMW+88YbRtm1b4+GHH7bX8YZrc/78eSMnJ8fIyckxJBmLFy82cnJy7K8dVucapKSkGB07djS2bt1q7Nu3z/jlL3/JqzGN2dKlS43OnTsb/v7+xvXXX29/JcRbSHL5eemll+x1bDabMW/ePCMiIsIICAgwBg4caBw4cMBzjfagy8PUW6/NO++8Y8TFxRkBAQFGz549jZUrVzps99brUlhYaMyYMcPo1KmTERgYaHTt2tWYM2eOUVxcbK/jDddm+/btLv9euffeew3DqN41uHjxovHAAw8YV111lREUFGTccsstRm5urum2sQQbAAAm8cwUAACTCFMAAEwiTAEAMIkwBQDAJMIUAACTCFMAAEwiTAEAMIkwBbzcpEmTNGrUqGbzPYAnEKZAEzVp0iRZLBY9/fTTDuWbNm1yWAHjz3/+s/r06aMWLVqodevW6tevn5555pmGbi7QrBGmQBMWGBioZ555xmkVjEqrV69WWlqapk+frk8//VQffvihHn74YV24cKGBWwo0b4Qp0IQNGTJEERER9qW4LvfOO+/ozjvv1OTJk3X11Vfr2muv1bhx4/TEE0841V20aJEiIyMVFhamadOmqbS01L6tpKREDz/8sDp06KAWLVrohhtu0I4dO+zb16xZo9atW+u9995TbGysQkJCNGzYMOXl5blt+969e9W+fXs9+eSTtb8AQCNBmAJNmNVq1VNPPaXnnntO3377rdP2iIgIffTRRzpx4kSVx9m+fbuOHj2q7du36+WXX9aaNWu0Zs0a+/bf/OY3+vDDD7V+/Xp99tlnuuOOOzRs2DB9/fXX9jpFRUVatGiRXnnlFf3jH/9Qbm6uHnroIZfft2PHDv3qV7/S/PnzNWfOnNqdPNCIEKZAEzd69Gj17dtX8+bNc9o2b948tW7dWl26dFGPHj00adIk/eUvf7Gvg1mpTZs2ev7559WzZ0/dcsstuvnmm7Vt2zZJ0tGjR7Vu3Tr99a9/VVJSkrp166aHHnpIAwYM0EsvvWQ/RmlpqVasWKGEhARdf/31euCBB+zH+Lm33npLt912m5YvX67777+/jq8G4BmEKdAMPPPMM3r55Zd18OBBh/LIyEhlZ2frwIEDmj59ukpLS3Xvvfdq2LBhDoF67bXXymq1Oux3+vRpSdK+fftkGIa6d++ukJAQ+ycrK0tHjx617xMcHKxu3bq5PEaljz/+WGPGjNHLL7+scePG1ek1ADzJ19MNAGDewIEDNXToUD366KOaNGmS0/a4uDjFxcVp2rRp+uCDD5SUlKSsrCzddNNNkiQ/Pz+H+haLxR62NptNVqtVe/fudQhcSQoJCbH/v6tjXL7CY7du3RQWFqYXX3xRN998s/z9/Wt9zkBjQpgCzcTTTz+tvn37qnv37lXW69WrlyTpxx9/rNZx+/Xrp/Lycp0+fVpJSUmm2ti2bVu98cYbGjx4sMaOHau//OUvTiEMNEXc5gWaieuuu07jx4/Xc889Zy+7//779cQTT+jDDz/UiRMn9NFHH2nixIlq166dEhMTq3Xc7t27a/z48Zo4caLeeOMNHTt2TJ988omeeeYZZWRk1Lid7du31/vvv68vv/xS48aNU1lZWY2PATQ2hCnQjDzxxBMOt1aHDBmijz76SHfccYe6d++uMWPGKDAwUNu2bVNYWFi1j/vSSy9p4sSJ+n//7/+pR48euu222/Txxx8rOjq6Vu2MiIjQ+++/rwMHDmj8+PEqLy+v1XGAxsJiXP5QAwAA1Ag9UwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJMIUwAATCJMAQAwiTAFAMAkwhQAAJP+f/zR80mTi+uhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (5, 5), dpi = 100)\n",
    "shenk_bns = np.linspace(0, 100, 50)\n",
    "a = func_df.query('occ >= 25 and occ_pct >= 0.50').abs_norm_shenkin\n",
    "b = unk_df.query('occ >= 25 and occ_pct >= 0.50').abs_norm_shenkin\n",
    "print(len(a), len(b))\n",
    "sns.histplot(a, label = \"Functional\", fill = True, color = \"#0C7BDC\", stat = \"proportion\", element = \"step\", bins = shenk_bns)\n",
    "sns.histplot(b, label = \"Unknown\", fill = True, color = \"#FFC20A\", stat = \"proportion\", element = \"step\", bins = shenk_bns)\n",
    "plt.ylabel(\"p\")\n",
    "plt.xlabel(\"NShenk\")\n",
    "plt.legend(edgecolor = \"k\")\n",
    "#plt.savefig(os.path.join(figs_dir, \"FUNC_vs_UNK_shenk.png\"))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/1211748.1.all.q/ipykernel_24921/3228072165.py:1: UserWarning: p-value floored: true value smaller than 0.001\n",
      "  scipy.stats.anderson_ksamp([a, b])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Anderson_ksampResult(statistic=104.75819493187241, critical_values=array([0.325, 1.226, 1.961, 2.718, 3.752, 4.592, 6.546]), pvalue=0.001)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.anderson_ksamp([a, b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shenkin average (FUNCTIONAL residues) =  20.0\n",
      "Shenkin average (UNKNOWN residues) =  32.0\n",
      "Kolmogorov-Smirnov two sample test result:\n",
      "KS statistic = 0.32\tp = 0.0\n",
      "\n",
      "Mann-Whitney U rank test result:\n",
      "U = 169182.0\tp = 0.0\n",
      "\n",
      "T-test result:\n",
      "T = -12.91\tp = 0.0\n"
     ]
    }
   ],
   "source": [
    "print(\"Shenkin average (FUNCTIONAL residues) = \", round(statistics.mean(a),0))\n",
    "print(\"Shenkin average (UNKNOWN residues) = \", round(statistics.mean(b),0))\n",
    "\n",
    "s, p = scipy.stats.ks_2samp(a, b, alternative=\"two-sided\")\n",
    "print(\"Kolmogorov-Smirnov two sample test result:\\nKS statistic = {}\\tp = {}\\n\".format(round(s,2), round(p,2)))\n",
    "\n",
    "s, p = scipy.stats.mannwhitneyu(a, b, alternative=\"two-sided\")\n",
    "print(\"Mann-Whitney U rank test result:\\nU = {}\\tp = {}\\n\".format(round(s,2), round(p,2)))\n",
    "\n",
    "s, p = scipy.stats.ttest_ind(a,b)\n",
    "print(\"T-test result:\\nT = {}\\tp = {}\".format(round(s,2), round(p,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FUNCTIONAL VS UNKNOWN: MES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "372 857\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAHACAYAAAD5vIKYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8x0lEQVR4nO3dfVhUdf7/8dcwAgOC1IoCFipqokbeBLsbFmp3eLO15nZjVpSbusvSZsr2bTVzK/ua3bguX/NuLctaK72utbt1+bVSq3QjbWrYVhqWqZjBKpaCogww5/cHMtvEzHBzmBluno/rmqs45z3nfD4enRdnzud8jsUwDEMAAKDFggLdAAAA2jvCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCpS6Ab0BY5HA598803ioyMlMViCXRzAAABYhiGKioq1KtXLwUFeT7/JEzd+OabbxQfHx/oZgAA2ohDhw7p/PPP97ieMHUjMjJSUt0fXrdu3QLcGgBAoJSXlys+Pt6ZC54Qpm7Uf7XbrVs3whQA0OglPwYgAQBgEmEKAIBJhCkAACZxzRQAGuFwOFRWVqbjx4/L4XAEujloJUFBQTrnnHMUHR3t9baXpiBMAcCLw4cP6+GHH9ZHH30U6KbAR5KTk/WHP/xB5513Xou3QZgCgAd2u11TpkzROeeco4ULF+r888+X1WoNdLPQSmpra/X1119r+fLlmjJlit566y2FhIS0aFuEKQB4cPDgQVVWVmrp0qUaPnx4oJsDH7jwwgsVExOj6dOnq7i4WAMGDGjRdhiABAAe1NbWSpJsNluAWwJfqj++NTU1Ld4GYQoAgEmEKQAAJnHNFABa4NDxah2r9N9tMt3DgxR/TrDf9ucLffv21axZszRr1iy/7fPAgQNKSEhQYWGhT697E6YA0EyHjlfr4pxDOlNj+G2fti4WfTQrvsmBOnXqVD3//PMNln/xxRctHmTTVGvXrtWsWbN0/Phxl+Xbt29X165dfbrvQCFMAaCZjlU6dKbGUGKPYIUH+/6Zx5XVhoqO1p0Jx5/T9PeNGzdOzz33nMuyHj16tG7jmiGQ+/Y1rpkCQAuFB1sUERrk81dLAzs0NFSxsbEur2nTpum6665zqZs1a5bGjBnj/HnMmDGaOXOm7rvvPv3oRz9SbGysHnroIZf3HD9+XL/61a8UExMjm82mpKQkbdq0SVu3btUvf/lLnThxQhaLRRaLxfnevn37Kicnx7mN4uJiTZw4UREREerWrZtuuukm/ec//3Guf+ihhzR8+HD95S9/Ud++fRUVFaWbb75ZFRUVzpo333xTl112mc455xx1795d11xzjfbt29eiPy8zCFMA/3X6S+lkoffX6S8D3Ur4wfPPP6+uXbvqX//6l5544gktWLBAeXl5kuqmVxw/fry2bdumdevWaffu3XrsscdktVo1cuRI5eTkqFu3biopKVFJSYnuvffeBts3DEPXXXedvv32W+Xn5ysvL0/79u3T5MmTXer27dun1157TZs2bdKmTZuUn5+vxx57zLn+1KlTys7O1vbt2/X2228rKChIkyZN8vu0j3zNC6DO6S+lHcOaVpvysRTm2+tuMG/Tpk2KiIhw/jx+/PgmX7McOnSoHnzwQUnSBRdcoGXLluntt9/W1Vdfrbfeeksffvih9uzZo4EDB0qS+vXr53xvVFSULBaLYmNjPW7/rbfe0r///W/t379f8fHxkqS//OUvuvDCC7V9+3b9+Mc/llQX3GvXrnU+nDsjI0Nvv/22Fi5cKEm6/vrrXba7Zs0a9ezZU7t371ZSUlKT+toaCFMAdWrPfnUWlihZwz3UVEqni/5bizbt8ssv18qVK50/d+3aVXPnzm3Se4cOHeryc1xcnI4cOSJJ2rVrl84//3xnkLbEnj17FB8f7wxSSRoyZIjOOecc7dmzxxmmffv2dQbpD9sh1Z25zp8/Xx988IHKysqcZ6TFxcWEKYAAsoZL1sjG69Dmde3atcHI3aCgIBmG6yjk6urqBu8NDnYdNWyxWJxBFRYWZrpthmHIYml4LfiHy721Q5KuvfZaxcfH6+mnn1avXr3kcDiUlJQku91uuo3NwTVTAOhEevTooZKSEpdlu3btatY2hg4dqq+//lp79+51uz4kJMQ5FaMnQ4YMUXFxsQ4dOuRctnv3bp04cUKDBw9uUjuOHTumPXv26IEHHtCVV16pwYMH67vvvmt6R1oRZ6YA0EKV1YYk3w90qdtP67jiiiv05JNP6oUXXlBqaqrWrVunTz/9VCNGjGjyNkaPHq1Ro0bp+uuv15IlSzRgwAB9/vnnslgsGjdunPr27auTJ0/q7bff1rBhwxQeHq7wcNdLB1dddZWGDh2qW2+9VTk5OaqpqVFWVpZGjx6tlJSUJrXj3HPPVffu3bV69WrFxcWpuLhYc+bMadafR2shTIE27Msyu07avX+QRoRYNCC6ZY+NQst0Dw+SrYtFRUcbfj3qK7YuFnUPN/9l4tixYzV//nzdd999OnPmjO68807dfvvt+uSTT5q1nY0bN+ree+/VlClTdOrUKQ0YMMA5ynbkyJHKzMzU5MmTdezYMT344IMNbq2xWCx67bXXdPfdd2vUqFEKCgrSuHHj9NRTTzW5DUFBQVq/fr1mzpyppKQkJSYmaunSpS63+fiLxfjhl+dQeXm5oqKidOLECXXr1i3QzUEn9WWZXSNyDjVeKKlwVrz5QD1ZKBVeJkWM8HzNtLairm7Ee3V1Hdznn3+u2267TevWrdOgQYNc1jGdYMfh7Tg3NQ84MwXaqPozUm+z7NTPjNPY2StaX/w5wc2ajQgdG2EKtHH1s+y4598b0wG4F/DRvCtWrFBCQoJsNpuSk5P17rvveqwtKSnRLbfcosTERAUFBXl88sDGjRs1ZMgQhYaGasiQIXr11Vd91HoAAAIcphs2bNCsWbM0b948FRYWKi0tTePHj1dxcbHb+qqqKvXo0UPz5s3TsGHuZ2opKCjQ5MmTlZGRoY8//lgZGRm66aab9K9//cuXXQEAdGIBDdMlS5Zo2rRpmj59ugYPHqycnBzFx8e7zNjxfX379tX//d//6fbbb1dUVJTbmpycHF199dWaO3euBg0apLlz5+rKK690mVwZAIDWFLAwtdvt2rlzp9LT012Wp6ena9u2bS3ebkFBQYNtjh071tQ2AQDwJmADkMrKylRbW6uYmBiX5TExMSotLW3xdktLS5u9zaqqKlVVVTl/Li8vb/H+AQCdT8AHIP1wbkZP8zX6cpuLFi1SVFSU8/X9iZcBAGhMwM5Mo6OjZbVaG5wxHjlypMGZZXPExsY2e5tz585Vdna28+fy8nICFYB3Zw5JNWX+21+XaMnm/8+lMWPGaPjw4Yw7aUTAwjQkJETJycnKy8vTpEmTnMvz8vI0ceLEFm83NTVVeXl5mj17tnPZ5s2bNXLkSI/vCQ0NVWhoaIv3CaCTOXNI2jlccpzx3z6DbFLyriYHqqcQfO211zRp0qQGT46BOQGdtCE7O1sZGRlKSUlRamqqVq9ereLiYmVmZkqqO2M8fPiwXnjhBed76p9ucPLkSR09elS7du1SSEiIhgwZIkm65557NGrUKD3++OOaOHGiXn/9db311lt67733/N4/AB1UTVldkHp79mtrqn+ObE2ZJL41a4sCes108uTJysnJ0YIFCzR8+HC98847ys3NVZ8+fSTVTdLww3tOR4wYoREjRmjnzp166aWXNGLECE2YMMG5fuTIkVq/fr2ee+45DR06VGvXrtWGDRv005/+1K99A9AJ1D/71ecv3wT2Qw89pOHDh+svf/mL+vbtq6ioKN18882qqPD88Pc333xTUVFRzpOcqVOn6rrrrtPixYsVFxen7t2766677nJ5Rup3332n22+/Xeeee67Cw8M1fvx4ffHFF5LqxrT06NFDGzdudNYPHz5cPXv2dP5cUFCg4OBgnTx5UlLduJhnnnlGkyZNUnh4uC644AK98cYbrfpn01wBH4CUlZWlAwcOqKqqSjt37tSoUaOc69auXautW7e61BuG0eB14MABl5obbrhBn3/+uex2u/bs2aNf/OIXfugJALQ/+/bt02uvvaZNmzZp06ZNys/Pdz795YfWr1+vm266SS+88IJuv/125/ItW7Zo37592rJli55//nmtXbtWa9euda6fOnWqduzYoTfeeEMFBQUyDEMTJkxQdXW1LBaLRo0a5fys/+6777R7925VV1dr9+7dkqStW7cqOTlZERERzm0+/PDDuummm/Tvf/9bEyZM0K233qpvv/229f+AmijgYQoACByHw6G1a9cqKSlJaWlpysjI0Ntvv92gbsWKFcrMzNTrr7/eYFzLueeeq2XLlmnQoEG65ppr9LOf/cy5jS+++EJvvPGGnnnmGaWlpWnYsGF68cUXdfjwYb322muS6q7v1ofpO++8o2HDhumKK65wLtu6dWuDx6pNnTpVU6ZM0YABA/Too4/q1KlT+vDDD1v1z6Y5CFMA6MT69u2ryMj/PnIvLi5OR44ccanZuHGjZs2apc2bN+vyyy9vsI0LL7xQVqvV7Tb27NmjLl26uFxq6969uxITE7Vnzx5JdWH62WefqaysTPn5+RozZozGjBmj/Px81dTUaNu2bRo9erTLPocOHer8/65duyoyMrJBu/2JMAWADqhbt246ceJEg+XHjx93eS5ncLDrM1ItFoscDtenEQ0fPlw9evTQc88953YUsLdteBo1/P37/5OSktS9e3fl5+c7w3T06NHKz8/X9u3bdfr0aV122WVN3mcgEKYA0AENGjRIO3bsaLB8+/btSkxMbNa2+vfvry1btuj111/X3Xff3az3DhkyRDU1NS4PGzl27Jj27t2rwYMHS5Lzuunrr7+uTz/9VGlpabroootUXV2tVatW6eKLL3Y5e26LCFMA6ICysrK0b98+3XXXXfr444+1d+9eLV++XGvWrNH//M//NHt7AwcO1JYtW5xf+TbVBRdcoIkTJ2rGjBl677339PHHH+u2227Teeed53LtdcyYMXrppZc0dOhQdevWzRmwL774YoPrpW0RDwcHgJaqrWyz++nbt6/effddzZs3T+np6Tpz5owGDhyotWvX6sYbb2xRMxITE/XPf/5TY8aMkdVq1R//+Mcmve+5557TPffco2uuuUZ2u12jRo1Sbm6uy1e1l19+uWpra12Cc/To0XrttdcaXC9tiywG02A0UF5erqioKJ04ccLl2gLgT7u+qVLaiq81oleIIkLdf4l0ssqhwm/sejfrfA3vZXIWr5OFUuFlUsSIunsb3amtqKsb8V5dXQf3+eef67bbbtO6des0aNCg/65oBzMgoek8Hmc1PQ84MwWA5rLF1wVbJ5ibF01DmAJAS9jixdR+qMcAJAAATCJMAQAwiTAFAMAkwhQAPAgKqvuI/P4TUNDx1B/f+uPdEoQpAHgQGxsrSSosLAxwS+BL9cc3Li6uxdtgNC8AeNCtWzdNmjRJTz31lKS65yn/cE5YtF/V1dUqLCzUU089pUmTJpmaspAwBQAv5s6dK0launRpgFsCX5k0aZLzOLcUYQq0Yf3DD2qgza6wEPdXZE5bHCoPD5F0vn8b1okEBQVp3rx5mjlzpkpKSgL6ZBK0rqCgIMXFxbXKJPqEKdBGhdj3aVfatY0XXiDttm+XNMTnberMIiMj2/yTSxA4hCnQRlmNk5Kk/Wf6qdoId1sTbKlUgu0rZy2AwCBMgTbujCNMdnV1u66W51QAbQK3xgAAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJjEpA0Amq+yyPt6a6QUNsA/bQHaAMIUQDOc/cgomtZ4acrHBCo6DcIUQNNZw6SIH0uq8VxTWymdLpJqK/zWLCDQCFMAzWMNC3QLgDaHAUgAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmEaYAAJhEmAIAYBJhCgCASYQpAAAmBTxMV6xYoYSEBNlsNiUnJ+vdd9/1Wp+fn6/k5GTZbDb169dPq1atalCTk5OjxMREhYWFKT4+XrNnz9aZM2d81QUAQCfXJZA737Bhg2bNmqUVK1bo0ksv1Z///GeNHz9eu3fvVu/evRvU79+/XxMmTNCMGTO0bt06vf/++8rKylKPHj10/fXXS5JefPFFzZkzR88++6xGjhypvXv3aurUqZKkP/3pT/7sHjqoL8vsOmk3vNZEhFg0IDrETy1qXFPaHFZVrUQ/tQfoaAIapkuWLNG0adM0ffp0SXVnlP/4xz+0cuVKLVq0qEH9qlWr1Lt3b+Xk5EiSBg8erB07dmjx4sXOMC0oKNCll16qW265RZLUt29fTZkyRR9++KF/OoUO7csyu0bkHGpSbeGs+DYRqE1t87DI/+i9kdLpGkNhVj80DOhAAhamdrtdO3fu1Jw5c1yWp6ena9u2bW7fU1BQoPT0dJdlY8eO1Zo1a1RdXa3g4GBddtllWrdunT788EP95Cc/0VdffaXc3FzdcccdHttSVVWlqqoq58/l5eUmeoaOrP7sLrFHsMKDLW5rKqsNFR2tbvRM0F+a0mZJiu9S93HgcPilWUCHErAwLSsrU21trWJiYlyWx8TEqLS01O17SktL3dbX1NSorKxMcXFxuvnmm3X06FFddtllMgxDNTU1+s1vftMgtL9v0aJFevjhh813Cp1GeLBFEaGehhy0zTTy3mYpzOI5aAF4F/ABSJYf/AM2DKPBssbqv79869atWrhwoVasWKGPPvpIr7zyijZt2qRHHnnE4zbnzp2rEydOOF+HDjXtazwAAKQAnplGR0fLarU2OAs9cuRIg7PPerGxsW7ru3Tpou7du0uS5s+fr4yMDOd12IsuukinTp3Sr371K82bN09BQQ1/fwgNDVVoaGhrdAsA0AkF7Mw0JCREycnJysvLc1mel5enkSNHun1Pampqg/rNmzcrJSVFwcHBkqTKysoGgWm1WmUYhvMsFgCA1hTQr3mzs7P1zDPP6Nlnn9WePXs0e/ZsFRcXKzMzU1Ld16+33367sz4zM1MHDx5Udna29uzZo2effVZr1qzRvffe66y59tprtXLlSq1fv1779+9XXl6e5s+fr5///OeyWhmiCABofQG9NWby5Mk6duyYFixYoJKSEiUlJSk3N1d9+vSRJJWUlKi4uNhZn5CQoNzcXM2ePVvLly9Xr169tHTpUudtMZL0wAMPyGKx6IEHHtDhw4fVo0cPXXvttVq4cKHf+wcA6BwCGqaSlJWVpaysLLfr1q5d22DZ6NGj9dFHH3ncXpcuXfTggw/qwQcfbK0mAgDgVcBH8wIA0N4RpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJXQLdAAC+1z/8oAba7AoL8fz7c8+gr/zYIkmnv5RqK7zXWCOlsAH+aQ9gAmEKdHAh9n3alXZtk+sNWX3YmrNOfyntGNa02pSPCVS0eYQp0MFZjZOSpP1n+qnaCPdYV2sYqrAHaXBsmO8bVX9GGpYoWT20qbZSOl3U+Nkr0AYQpkAnccYRJru6elxf6zB0xmH4sUWqC1JrpH/3CfgAA5AAADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABM6hLoBgAdVdERu9f1ESEWDYgO8VNrWk9ltUO1Ds/rgxyGuvqvOUCbQJgCrcx69vue6X890mht4az4dhWoldUO7fi6ymtNhNWui6OkQ9/VKD7CTw0DAowwBVpZWHCQUs4P8Xr2VlltqOhotU7aDf81rBXU9yk82KIgi/sam6VuRWVN++obYAZhCvhAWHBjwxG8JG07EGSRrB7S1EPGAh1awAcgrVixQgkJCbLZbEpOTta7777rtT4/P1/Jycmy2Wzq16+fVq1a1aDm+PHjuuuuuxQXFyebzabBgwcrNzfXV10AAHRyAQ3TDRs2aNasWZo3b54KCwuVlpam8ePHq7i42G39/v37NWHCBKWlpamwsFD333+/Zs6cqY0bNzpr7Ha7rr76ah04cEB//etfVVRUpKefflrnnXeev7oFAOhkAvo175IlSzRt2jRNnz5dkpSTk6N//OMfWrlypRYtWtSgftWqVerdu7dycnIkSYMHD9aOHTu0ePFiXX/99ZKkZ599Vt9++622bdum4OBgSVKfPn380yEAQKcUsDNTu92unTt3Kj093WV5enq6tm3b5vY9BQUFDerHjh2rHTt2qLq6WpL0xhtvKDU1VXfddZdiYmKUlJSkRx99VLW1tR7bUlVVpfLycpcXAABNFbAwLSsrU21trWJiYlyWx8TEqLS01O17SktL3dbX1NSorKxMkvTVV1/pr3/9q2pra5Wbm6sHHnhAf/zjH7Vw4UKPbVm0aJGioqKcr/j4eJO9AwB0JgEfzWuxuI79MwyjwbLG6r+/3OFwqGfPnlq9erWsVquSk5P1zTff6Mknn9Qf/vAHt9ucO3eusrOznT+Xl5cTqIBJNvte6WSw+5WVRU3fUGO11kgpbEDTtwf4QMDCNDo6WlartcFZ6JEjRxqcfdaLjY11W9+lSxd1795dkhQXF6fg4GBZrVZnzeDBg1VaWiq73a6QkIY3yIeGhio0NNRsl4Bm8zZL0pFvq5Xox7a0FodR92+vz5FfS43OW+HtI+jsuqJpje805WMCFQEVsDANCQlRcnKy8vLyNGnSJOfyvLw8TZw40e17UlNT9be//c1l2ebNm5WSkuIcbHTppZfqpZdeksPhUFBQ3bfYe/fuVVxcnNsgBQKhKbMkDYv8VukjJYtFUjua/6DKsOnD4xfpwp5WdQ3xdtdpF8ka5nm1NUyK+LGkGs81tZXS6SKptqKlzQVaRUC/5s3OzlZGRoZSUlKUmpqq1atXq7i4WJmZmZLqvn49fPiwXnjhBUlSZmamli1bpuzsbM2YMUMFBQVas2aNXn75Zec2f/Ob3+ipp57SPffco7vvvltffPGFHn30Uc2cOTMgfQTcacosSQNtdb8gBrXDaRDOOGxyBIX+97eGlvIWtkAbEtAwnTx5so4dO6YFCxaopKRESUlJys3Ndd7KUlJS4nLPaUJCgnJzczV79mwtX75cvXr10tKlS523xUhSfHy8Nm/erNmzZ2vo0KE677zzdM899+j3v/+93/sHeNPYLElhIQGfUwVAEwV8AFJWVpaysrLcrlu7dm2DZaNHj9ZHH33kdZupqan64IMPWqN5AAA0il99AQAwiTAFAMAkwhQAAJMIUwAATAr4ACQAbUtltef7dbytAzozwhSApLOTQ0gqOlrd5FoAdQhTAJKkIItFkaGS0chsSxZLXS2A/yJMATgFWSxqhxMuAQHHACQAAEwiTAEAMMl0mBqG4XymKAAAnVGLw3TNmjVKSkqSzWaTzWZTUlKSnnnmmdZsGwAA7UKLBiDNnz9ff/rTn3T33XcrNTVVklRQUKDZs2frwIED+t///d9WbSQAAG1Zi8J05cqVevrppzVlyhTnsp///OcaOnSo7r77bsIUANCptChMa2trlZKS0mB5cnKyampqTDcK+L4vy+w6afd+XT4ixKIB0SF+ahEAuGpRmN52221auXKllixZ4rJ89erVuvXWW1ulYYBUF6Qjcg41qbZwVjyBCiAgWjxpw5o1a7R582ZdcsklkqQPPvhAhw4d0u23367s7Gxn3Q8DF2iO+jPSxB7BCg92P5tAZbWhoqPVjZ69AoCvtChMP/30U1188cWSpH379kmSevTooR49eujTTz911lmYcgytJDzYoohQT4PPmXwdQGC1KEy3bNnS2u0AAKDdYgYkAABMIkwBADCJMAUAwCTCFAAAk3ieKfA9jU0QUXTE7sfWtG+V1d5HWVuDpPBgfp9Hx0CYAmc1Z4IIKxngUf0dcUVHqxutTTk/lEBFh0CYAmc1ZYIIqS5IwwgAj4IsFkWGSt6ezOgw6ibbqOUWYXQQhCnwA94niEBTBFkskrc5WxzMVoWOhTAF0P5VFnlfb42Uwgb4py3olAhTAO3Y2Y+wommNl6Z8TKDCZwhTAO2XNUyK+LEkL49+rK2UThdJtRV+axY6H8IUQPtmDQt0CwAmbQAAwCzCFAAAk/iaF50GsxsB8BXCFJ0CsxsB8CXCFJ0CsxsB8CXCFJ0KsxsB8AU+VQAAMIkwBQDAJMIUAACTCFMAAEwiTAEAMInRvEAHcPDbap0OrXK77si31Ur0c3uAzoYwBdoxy9lbZhe89a0+rvjabc2wyG+VPvJsLc/kBnyCMAXasSDVpemgHsEKigxxWzPQFuxSC6D1EaZABxAWEqQIw/0QiLAQhkYAvsa/MgAATOLMFPieXsEHFBZ0ymvNaUdXfVPd1z8NQuupLPK+3hophQ3wT1vQ4RCmwFm9gg/omYT0JtVO37+ZQG03zn7MFU1rvDTlYwIVLUKYAmfVn5EerOqvKkeY25rQoNPqE7qv0bNXtCHWMCnix5JqPNfUVkqni6TaCr81Cx0LYQr8QJUjTKeNru5XOvzbFrQSq/tfjoDWQpgCPsC1V6BzIUyBVsa1V6DzIUyBVsa1V6DzIUwBH+HaK9B5MGkDAAAmEaYAAJhEmAIAYBJhCgCASQEP0xUrVighIUE2m03Jycl69913vdbn5+crOTlZNptN/fr106pVqzzWrl+/XhaLRdddd10rtxoAgP8KaJhu2LBBs2bN0rx581RYWKi0tDSNHz9excXFbuv379+vCRMmKC0tTYWFhbr//vs1c+ZMbdy4sUHtwYMHde+99yotLc3X3QAAdHIBDdMlS5Zo2rRpmj59ugYPHqycnBzFx8dr5cqVbutXrVql3r17KycnR4MHD9b06dN15513avHixS51tbW1uvXWW/Xwww+rX79+/ugKAKATC1iY2u127dy5U+nprjPFpKena9u2bW7fU1BQ0KB+7Nix2rFjh6qrq53LFixYoB49emjatCY8JUJSVVWVysvLXV4AADRVwMK0rKxMtbW1iomJcVkeExOj0tJSt+8pLS11W19TU6OysjJJ0vvvv681a9bo6aefbnJbFi1apKioKOcrPj6+mb0BAHRmAR+AZLFYXH42DKPBssbq65dXVFTotttu09NPP63o6Ogmt2Hu3Lk6ceKE83Xo0KFm9AAA0NkFbDrB6OhoWa3WBmehR44caXD2WS82NtZtfZcuXdS9e3d99tlnOnDggK699lrneoejbt62Ll26qKioSP3792+w3dDQUIWGhprtEgCgkwrYmWlISIiSk5OVl5fnsjwvL08jR450+57U1NQG9Zs3b1ZKSoqCg4M1aNAgffLJJ9q1a5fz9fOf/1yXX365du3axde3AACfCOhE99nZ2crIyFBKSopSU1O1evVqFRcXKzMzU1Ld16+HDx/WCy+8IEnKzMzUsmXLlJ2drRkzZqigoEBr1qzRyy+/LEmy2WxKSkpy2cc555wjSQ2WAwDQWgIappMnT9axY8e0YMEClZSUKCkpSbm5uerTp48kqaSkxOWe04SEBOXm5mr27Nlavny5evXqpaVLl+r6668PVBcAAAj8I9iysrKUlZXldt3atWsbLBs9erQ++uijJm/f3TYAAGhNAR/NCwBAexfwM1N0XF+W2XXSbnitiQixaEB0iJ9ahLamstr7U9KtQVJ4ML/zo+0jTOETX5bZNSKnaffrFs6KJ1A7mfrbxYuOVnsvlJRyfiiBijaPMIVP1J+RJvYIVniw+0k4KqsNFR2tbvTsFR1PkMWiyFDJ8HLoHUbd35Fa7yevQJtAmMKnwoMtigj1dFbBp2RnFmSxSJ4nO6tLU6Cd4LsTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMYjQv0AHEh+xr0ToArYMwBdqxWlklSffF3dvkWgCtjzAF2jG7YdOe08NkVa3XulpZZTdsfmoV0PkQpkA7R0gCgUeYAgHEtU6gYyBMgQDgWifQsRCmQABwrRPoWAhTIEAISaDjYNIGAABMIkwBADCJMAUAwCTCFAAAkxiABLQA94cC+D7CFGgG7g8F4A5hCjQD94cCcIcwBZqJkATwQwxAAgDAJMIUAACTCFMAAEzimimANq2y2uF1vTVICg/2fl5QWe1QrZfNBDkMdW1J44CzCFMAbZLFUvffoqPVjdamnB/qMVArqx3a8XWV1/dHWO26OEo69F2N4iOa3VSAMIXv9A8/qIE2u8JC3H/InbY4VB4eIul8/zYM7UKQxaLIUMkwPNc4DKmy2vB61lm/LjzYoiCL+xrb2eSurPGyM8ALwhQ+EWLfp11p1zZeeIG0275d0hCftwntT5DFInkIQEl1adrkbUlWD2nqbRdAUxCm8AmrcVKStP9MP1Ub4W5rgi2VSrB95awFgPaKMIVPnXGEye5haEett+/vAKAd4dYYAABMIkwBADCJMAUAwCTCFAAAkxiABABn2ex7pZPBnguskVLYAP81CO0GYQqg03MYdQ9y73Pk19KRRopTPiZQ0QBhCqDTqzJs+vD4Rbqwp1VdQzxM4VBbKZ0ukmor/Ns4tAuEKQBIOuOwyREUWjdzPtBM/K0BAMAkwhQAAJMIUwAATCJMAQAwiQFIANAclUXe13MvaqdEmKLDKDpib9E6oGnOflwWTWu8lHtROx3CFO1e/Z0M0//a2N323PUAE6xhUsSPJdV4ruFe1E6LMEW7FxYcpJTzQ1Tr8F5nDaqrBVrMGhboFqCNIkzRIRCSAAKJTyAAAEwiTAEAMIkwBQDAJMIUAACTCFMAAEwiTAEAMIkwBQDAJMIUAACTAh6mK1asUEJCgmw2m5KTk/Xuu+96rc/Pz1dycrJsNpv69eunVatWuax/+umnlZaWpnPPPVfnnnuurrrqKn344Ye+7AIAoJMLaJhu2LBBs2bN0rx581RYWKi0tDSNHz9excXFbuv379+vCRMmKC0tTYWFhbr//vs1c+ZMbdy40VmzdetWTZkyRVu2bFFBQYF69+6t9PR0HT582F/dAgB0MgEN0yVLlmjatGmaPn26Bg8erJycHMXHx2vlypVu61etWqXevXsrJydHgwcP1vTp03XnnXdq8eLFzpoXX3xRWVlZGj58uAYNGqSnn35aDodDb7/9tr+6BQDoZAIWpna7XTt37lR6errL8vT0dG3bts3tewoKChrUjx07Vjt27FB1dbXb91RWVqq6ulo/+tGPPLalqqpK5eXlLi8AAJoqYGFaVlam2tpaxcTEuCyPiYlRaWmp2/eUlpa6ra+pqVFZWZnb98yZM0fnnXeerrrqKo9tWbRokaKiopyv+Pj4ZvYGANCZBXwAksVicfnZMIwGyxqrd7dckp544gm9/PLLeuWVV2Sz2Txuc+7cuTpx4oTzdejQoeZ0AQDQyQXsEWzR0dGyWq0NzkKPHDnS4OyzXmxsrNv6Ll26qHv37i7LFy9erEcffVRvvfWWhg4d6rUtoaGhCg0NbUEvAAAI4JlpSEiIkpOTlZeX57I8Ly9PI0eOdPue1NTUBvWbN29WSkqKgoODncuefPJJPfLII3rzzTeVkpLS+o0H0KZUVjtUUeX+VVndyFPjgVYQ0IeDZ2dnKyMjQykpKUpNTdXq1atVXFyszMxMSXVfvx4+fFgvvPCCJCkzM1PLli1Tdna2ZsyYoYKCAq1Zs0Yvv/yyc5tPPPGE5s+fr5deekl9+/Z1nslGREQoIiLC/50E4DP1V3eKjrofgOiuFvCFgIbp5MmTdezYMS1YsEAlJSVKSkpSbm6u+vTpI0kqKSlxuec0ISFBubm5mj17tpYvX65evXpp6dKluv766501K1askN1u1w033OCyrwcffFAPPfSQX/oFwD+CLBZFhkpnh054ZLHU1QK+EtAwlaSsrCxlZWW5Xbd27doGy0aPHq2PPvrI4/YOHDjQSi0D0B4EWSwSOYkAC/hoXgAA2ruAn5kCraFX8AGFBZ3yWnPa0VXfVPf1T4MAdCqEKdq9XsEH9ExCeuOFkqbv30ygAmh1hCnavfoz0oNV/VXlCHNbExp0Wn1C9zV69goALUGYosOocoTptNHV/UpuNQTgQwxAAgDAJM5MEXAHv63W6dAqt+uKjtj93Bp0Zo3NlmQNksKDOQdBQ4QpAqb+HvoFb32rjyu+9lpr5fMLPtScmZRSzg8lUNEAYYqACTp7p/2gHsEKigzxWGcNksL48IIPNWUmJYchVVYbquX6O9wgTBFwYSFBijAISwRWozMpORqZsxCdGp9gAACYRJgCAGASYQoAgElcMwWA1lZZ5H29NVIKG+CftsAvCFMAaDVnP1KLpjVemvIxgdqBEKYA0FqsYVLEjyXVeK6prZROF0m1FX5rFnyPMAWA1mR1/7AFdGwMQAIAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJSRvQIl+W2XXS7vn5jke+rVaiH9sDAIFEmKLZviyza0TOIa81wyK/VfpIyWKRZPKZyr2CDygs6JTH9fEh+8ztAAgEJsPvUAhTNFv9GWlij2CFB1vc1gy0BUuSguR+fVP1Cj6gZxLSm1RbK6upfQH+wWT4HRFhihYLD7YoItT9ZfewkNa5HF9/Rnqwqr+qHJ7nPK2VVXbD1ir7BHyKyfA7JMIU7UKVI0ynja6BbgbQOpgMv8NhNC8AACYRpgAAmESYAgBgEmEKAIBJDEBqJxqbJEGSIkIsGhAd4qcWtR5v94m29j2k/twX4ElltUO1Ds/rgxyGGG7XvhCm7UBTJkmoVzgrvt0Eav19offF3dvk2vawL8CbymqHdnxd5bUmwmrXxVHSoe9qFB/hp4bBFMK0HWjKJAmV1YaKjlY3evbaltgNm/acHiarar3WtcY9pP7cF+BN/RlpeLBFQR7mNLFZ6lZU1rSff8+dHWHajnibJEHy8p1RG+bP4CIk0ZYEWSSrhzQ1N28YAoEBSAAAmESYAgBgEmEKAIBJhCkAACYxAAkA2iibfa90MthzAc88bTMIUwBoYxxG3b3OfY78WjrSSDHPPG0TCFMAaGOqDJs+PH6RLuxpVdcQDzfK8MzTNoUwBYA26IzDJkdQqGRlaEt7QJiiodNfev1tN6yqWv3DT0q6wH9tAoA2jDCFq9NfSjuGeS1JlLQrTbql6E0dVz//tAsA2jDCFK7qz0jDEiVruNuS01WnFFa9V+FBp3ScqUMBgDCFB9bwumH3bjgsJCgAfB9h2k70Dz+ogTa7wkLcD0Y4bXGoPDxE0vn+bRgAgDBtD0Ls+7Qr7drGCy+Qdtu3Sxri8zYBAP6LMG0HrMZJSdL+M/1Ubbi/jhlsqVSC7StnLQDAfwjTduSMI0x2dXW7rtZo2nXML8vsXh8gHlZVrcQWtQ7oHCqrzT07uDnv91Yb5DA8fBogEAjTTuTLMrtG5BzyWjMs8j96b6R0usZQmNVPDQPaAcvZiYiKjla36vZauq8Iq10XR0ml//lMsd52xPy9fkGYdiL1Z6SJPYIVHuz+X3J8l7q/Eg5zv3wDHU6QxaLIUKmJXwJ5ZbHUbc/MvoLPfnzHfvNr6ZtGdsj8vT5HmHZC4cEWRYS6HxUc5u3XZaCTC7JYJD/9E2lsX9UO5u9tSwhTAGinmL+37eAIAABgUsDPTFesWKEnn3xSJSUluvDCC5WTk6O0tDSP9fn5+crOztZnn32mXr166b777lNmZqZLzcaNGzV//nzt27dP/fv318KFCzVp0iRfd6VN8PYw4bCqag2L/I8G2oI9Tv7QM+irJu+rj+0rhTncbyc+ZF+TtwMgwBp5uIWk1hvI5M99+VFAw3TDhg2aNWuWVqxYoUsvvVR//vOfNX78eO3evVu9e/duUL9//35NmDBBM2bM0Lp16/T+++8rKytLPXr00PXXXy9JKigo0OTJk/XII49o0qRJevXVV3XTTTfpvffe009/+lN/d9FvmvIw4URJ741s2vYMeR7KW79ufvz/NLqdWi/bAdAGNOHhFk5mBzL5c19+ZjGM1hib1jI//elPdfHFF2vlypXOZYMHD9Z1112nRYsWNaj//e9/rzfeeEN79uxxLsvMzNTHH3+sgoICSdLkyZNVXl6u//f//p+zZty4cTr33HP18ssvN6ld5eXlioqK0okTJ9StW7eWdq/VFO3/UIlfX649lRfKrgi3NbUOQzU1p70ORjhlN1R01K6wYIusHgYa1RqGKuxBGhwbpUgPg5QqqhzaU3pCkSEOj9uR6oLUbtga6R2Alqh1GDppNzSiV6jHf6uqrZBOFkoj3pMiRrivOVkoFV7m9eEWzoFM3rbTFP7cVytpah4E7MzUbrdr586dmjNnjsvy9PR0bdu2ze17CgoKlJ6e7rJs7NixWrNmjaqrqxUcHKyCggLNnj27QU1OTk6rtr8tamwwgiPIoZO1VZLVImuQhzB1GDrjaPz3qzMOm7o4PG8HQDvj5eEW7XpffhKwMC0rK1Ntba1iYmJclsfExKi0tNTte0pLS93W19TUqKysTHFxcR5rPG1TkqqqqlRVVeX8+cSJE5LqfiMx69+f56my4qCpbThO7VVciKHqM+WyO2rd1xiSag0d+k+QbB7uIT1TbUhVDlXXWFTrIQNbazsAfKsp/1a7qFLnGIZ2f/CMyi393dZ0M/ZpSI2h45XfqUZVbmuasp2m8Oe+JCk8so+GDrra1Dbqc6CxL3EDPgDJ8oOvCQ3DaLCssfofLm/uNhctWqSHH364wfL4+HjPDQ+IokA3AEC7tKIJNbtbaTtN4c99tY6KigpFRUV5XB+wMI2OjpbVam1wxnjkyJEGZ5b1YmNj3dZ36dJF3bt391rjaZuSNHfuXGVnZzt/djgc+vbbb9W9e3evIdyY8vJyxcfH69ChQ23i2qsvdZa+dpZ+SvS1I+os/ZRar6+GYaiiokK9evXyWhewMA0JCVFycrLy8vJcblvJy8vTxIkT3b4nNTVVf/vb31yWbd68WSkpKQoODnbW5OXluVw33bx5s0aO9DyMNTQ0VKGhoS7LzjnnnOZ2yaNu3bp1+L+49TpLXztLPyX62hF1ln5KrdNXb2ek9QL6NW92drYyMjKUkpKi1NRUrV69WsXFxc77RufOnavDhw/rhRdekFQ3cnfZsmXKzs7WjBkzVFBQoDVr1riM0r3nnns0atQoPf7445o4caJef/11vfXWW3rvvfcC0kcAQMcX0DCdPHmyjh07pgULFqikpERJSUnKzc1Vnz59JEklJSUqLi521ickJCg3N1ezZ8/W8uXL1atXLy1dutR5j6kkjRw5UuvXr9cDDzyg+fPnq3///tqwYUOHvscUABBYAR+AlJWVpaysLLfr1q5d22DZ6NGj9dFHH3nd5g033KAbbrihNZpnSmhoqB588MEGXyF3RJ2lr52lnxJ97Yg6Sz8l//c1oJM2AADQETDRPQAAJhGmAACYRJgCAGASYQoAgEmEaStbuHChRo4cqfDw8CZP/DB16lRZLBaX1yWXXOLbhprUkn4ahqGHHnpIvXr1UlhYmMaMGaPPPvvMtw1tBd99950yMjIUFRWlqKgoZWRk6Pjx417f016O6YoVK5SQkCCbzabk5GS9++67Xuvz8/OVnJwsm82mfv36adWqVX5qqXnN6evWrVsbHD+LxaLPP//cjy1uvnfeeUfXXnutevXqJYvFotdee63R97TXY9rcvvr6mBKmrcxut+vGG2/Ub37zm2a9b9y4cSopKXG+cnNzfdTC1tGSfj7xxBNasmSJli1bpu3btys2NlZXX321KioaeVBwgN1yyy3atWuX3nzzTb355pvatWuXMjIyGn1fWz+m9c8TnjdvngoLC5WWlqbx48e73Nv9ffXPE05LS1NhYaHuv/9+zZw5Uxs3bvRzy5uvuX2tV1RU5HIML7jgAj+1uGVOnTqlYcOGadmyZU2qb8/HtLl9reezY2rAJ5577jkjKiqqSbV33HGHMXHiRJ+2x1ea2k+Hw2HExsYajz32mHPZmTNnjKioKGPVqlU+bKE5u3fvNiQZH3zwgXNZQUGBIcn4/PPPPb6vPRzTn/zkJ0ZmZqbLskGDBhlz5sxxW3/fffcZgwYNcln261//2rjkkkt81sbW0ty+btmyxZBkfPfdd35onW9IMl599VWvNe35mH5fU/rq62PKmWkbsXXrVvXs2VMDBw7UjBkzdOTIkUA3qVXt379fpaWlLs+jDQ0N1ejRoz0+v7YtKCgoUFRUlMsMWpdccomioqIabXdbPqb1zxP+4fOBW/I84R07dqi6utpnbTWrJX2tN2LECMXFxenKK6/Uli1bfNnMgGivx9QMXx1TwrQNGD9+vF588UX985//1B//+Edt375dV1xxhcszVtu7+if5NPdZs4FWWlqqnj17Nljes2dPr+1u68fUF88Tbqta0te4uDitXr1aGzdu1CuvvKLExERdeeWVeuedd/zRZL9pr8e0JXx9TAM+nWB78NBDD7l93un3bd++XSkpKS3a/uTJk53/n5SUpJSUFPXp00d///vf9Ytf/KJF22wJX/dTav6zZn2lqX2VGrZZarzdbeWYNsYXzxNuq5rT18TERCUmJjp/Tk1N1aFDh7R48WKNGjXKp+30t/Z8TJvD18eUMG2C3/72t7r55pu91vTt27fV9hcXF6c+ffroiy++aLVtNoUv+xkbGyup7jfhuLg45/LGnjXrK03t67///W/95z//abDu6NGjzWp3oI6pJ756nnBb1JK+unPJJZdo3bp1rd28gGqvx7S1tOYxJUybIDo6WtHR0X7b37Fjx3To0CGX0PEHX/YzISFBsbGxysvL04gRIyTVXcvKz8/X448/7pN9etPUvqampurEiRP68MMP9ZOf/ESS9K9//UsnTpzw+ozcHwrUMfXEV88Tbota0ld3CgsL28zxay3t9Zi2llY9pj4Z1tSJHTx40CgsLDQefvhhIyIiwigsLDQKCwuNiooKZ01iYqLxyiuvGIZhGBUVFcbvfvc7Y9u2bcb+/fuNLVu2GKmpqcZ5551nlJeXB6objWpuPw3DMB577DEjKirKeOWVV4xPPvnEmDJlihEXF9em+2kYhjFu3Dhj6NChRkFBgVFQUGBcdNFFxjXXXONS0x6P6fr1643g4GBjzZo1xu7du41Zs2YZXbt2NQ4cOGAYhmHMmTPHyMjIcNZ/9dVXRnh4uDF79mxj9+7dxpo1a4zg4GDjr3/9a6C60GTN7euf/vQn49VXXzX27t1rfPrpp8acOXMMScbGjRsD1YUmqaiocP5blGQsWbLEKCwsNA4ePGgYRsc6ps3tq6+PKWHayu644w5DUoPXli1bnDWSjOeee84wDMOorKw00tPTjR49ehjBwcFG7969jTvuuMMoLi4OTAeaqLn9NIy622MefPBBIzY21ggNDTVGjRplfPLJJ/5vfDMdO3bMuPXWW43IyEgjMjLSuPXWWxsMr2+vx3T58uVGnz59jJCQEOPiiy828vPznevuuOMOY/To0S71W7duNUaMGGGEhIQYffv2NVauXOnnFrdcc/r6+OOPG/379zdsNptx7rnnGpdddpnx97//PQCtbp762z9++LrjjjsMw+hYx7S5ffX1MeURbAAAmMStMQAAmESYAgBgEmEKAIBJhCkAACYRpgAAmESYAgBgEmEKAIBJhCkAACYRpkAnMXXqVFksFmVmZjZYl5WVJYvFoqlTp7rU/vA1btw453sKCwt1zTXXqGfPnrLZbOrbt68mT57c4R7dBTQFYQp0IvHx8Vq/fr1Onz7tXHbmzBm9/PLL6t27t0vtuHHjVFJS4vJ6+eWXJdU9WeSqq65SdHS0/vGPf2jPnj169tlnFRcXp8rKSr/2CWgLeGoM0IlcfPHF+uqrr/TKK6/o1ltvlSS98sorio+PV79+/VxqQ0NDnY/O+6Ft27apvLxczzzzjLp0qfsYSUhI0BVXXOHbDgBtFGemQCfzy1/+Us8995zz52effVZ33nlns7YRGxurmpoavfrqq2J6b4AwBTqdjIwMvffeezpw4IAOHjyo999/X7fddluDuk2bNikiIsLl9cgjj0iqe6jy/fffr1tuuUXR0dEaP368nnzySbcPUgc6A77mBTqZ6Oho/exnP9Pzzz8vwzD0s5/9zO2D0i+//HKtXLnSZdmPfvQj5/8vXLhQ2dnZ+uc//6kPPvhAq1at0qOPPqp33nlHF110kc/7AbQlhCnQCd1555367W9/K0lavny525quXbtqwIABXrfTvXt33Xjjjbrxxhu1aNEijRgxQosXL9bzzz/f6m0G2jLCFOiExo0bJ7vdLkkaO3Zsq2wzJCRE/fv316lTp1ple0B7QpgCnZDVatWePXuc/+9OVVWVSktLXZZ16dJF0dHR2rRpk9avX6+bb75ZAwcOlGEY+tvf/qbc3FyXwU1AZ0GYAp1Ut27dvK5/8803FRcX57IsMTFRn3/+uYYMGaLw8HD97ne/06FDhxQaGqoLLrhAzzzzjDIyMnzZbKBNshiMawcAwBRujQEAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADCJMAUAwCTCFAAAkwhTAABMIkwBADDp/wNntB0/zQGyXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 500x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (5, 5), dpi = 100)\n",
    "mes_bns = np.linspace(-1.5, 1.5, 40)\n",
    "a = func_df.query('occ >= 25 and occ_pct >= 0.50 & human_occ > 10 & human_occ_pct > 0.50 & variants >= 2').log_oddsratio\n",
    "b = unk_df.query('occ >= 25 and occ_pct >= 0.50 & human_occ > 10 & human_occ_pct > 0.50 & variants >= 2').log_oddsratio\n",
    "print(len(a), len(b))\n",
    "sns.histplot(a, label = \"Functional\", fill = True, color = \"#0C7BDC\", stat = \"proportion\", element = \"step\", bins = mes_bns)\n",
    "sns.histplot(b, label = \"Unknown\", fill = True, color = \"#FFC20A\", stat = \"proportion\", element = \"step\", bins = mes_bns)\n",
    "plt.ylabel(\"p\")\n",
    "plt.xlabel(\"MES\")\n",
    "plt.legend(edgecolor = \"k\")\n",
    "#plt.savefig(os.path.join(figs_dir, \"FUNC_vs_UNK_MES.png\"))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/1211748.1.all.q/ipykernel_24921/3228072165.py:1: UserWarning: p-value floored: true value smaller than 0.001\n",
      "  scipy.stats.anderson_ksamp([a, b])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Anderson_ksampResult(statistic=45.72428785837081, critical_values=array([0.325, 1.226, 1.961, 2.718, 3.752, 4.592, 6.546]), pvalue=0.001)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.stats.anderson_ksamp([a, b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MES average (FUNCIONAL residues) =  -0.2\n",
      "MES average (UNKNOWN residues) =  -0.0\n",
      "Kolmogorov-Smirnov two sample test result:\n",
      "KS statistic = 0.26\tp = 0.0\n",
      "\n",
      "Mann-Whitney U rank test result:\n",
      "U = 111507.0\tp = 0.0\n",
      "\n",
      "T-test result:\n",
      "T = -8.35\tp = 0.0\n"
     ]
    }
   ],
   "source": [
    "print(\"MES average (FUNCIONAL residues) = \", round(statistics.mean(a),1))\n",
    "print(\"MES average (UNKNOWN residues) = \", round(statistics.mean(b),1))\n",
    "\n",
    "s, p = scipy.stats.ks_2samp(a, b, alternative=\"two-sided\")\n",
    "print(\"Kolmogorov-Smirnov two sample test result:\\nKS statistic = {}\\tp = {}\\n\".format(round(s,2), round(p,2)))\n",
    "\n",
    "s, p = scipy.stats.mannwhitneyu(a, b, alternative=\"two-sided\")\n",
    "print(\"Mann-Whitney U rank test result:\\nU = {}\\tp = {}\\n\".format(round(s,2), round(p,2)))\n",
    "\n",
    "s, p = scipy.stats.ttest_ind(a, b)\n",
    "print(\"T-test result:\\nT = {}\\tp = {}\".format(round(s,2), round(p,2)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:myenv2]",
   "language": "python",
   "name": "conda-env-myenv2-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}