radiogenomic_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CBTN learning lab for multi-modal imaging\n",
    "## Radiogenomic Machine Learning\n",
    "#### October 10, 2024\n",
    "\n",
    "This notebook illustrates how to run an unsupervised machine learning analysis on Flywheel using processed tabular data.\n",
    "\n",
    "- pLGG subjects with multi-parametric images selected from CBTN.\n",
    "- Molecular subtype information downloaded from: https://github.com/d3b-center/OpenPedCan-analysis/tree/dev/analyses/molecular-subtyping-LGAT/results\n",
    "\n",
    "Based on: [Haldar et al. 2023 Neoplasia](https://pubmed.ncbi.nlm.nih.gov/36566592/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import flywheel\n",
    "\n",
    "random_state = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load CSV feature sets\n",
    "first we need to load our group-level feature files that are attached to the \"CBTN_proc_data\" project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 338,
   "metadata": {},
   "outputs": [],
   "source": [
    "# access the Flywheel instance\n",
    "your_api_key = 'enter-your-api-key'\n",
    "fw = flywheel.Client(your_api_key)\n",
    "\n",
    "# get the project container\n",
    "proj_container = fw.lookup('cbtn-summit/CBTN_proc_data')\n",
    "\n",
    "# loop throug project attachments and download feature CSVs\n",
    "for file_attachment in proj_container.files:\n",
    "    if file_attachment.name == 'molec_features.csv':\n",
    "        proj_container.download_file(file_attachment.name, file_attachment.name)\n",
    "        molec_df = pd.read_csv('molec_features.csv')\n",
    "        \n",
    "    elif file_attachment.name == 'radiomic_features.csv':\n",
    "        proj_container.download_file(file_attachment.name, file_attachment.name)\n",
    "        rad_df = pd.read_csv('radiomic_features.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter to only subjects with molec & radiomic data\n",
    "X = molec_df.merge(rad_df, how='left', on='SubjectID')\n",
    "y = X[['OpenPedCan Molecular Subtype']]\n",
    "X = X.drop(columns='OpenPedCan Molecular Subtype')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare radiomic features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.impute import SimpleImputer\n",
    "\n",
    "# pre-process the radiomic features\n",
    "X_radiomic = X.copy()\n",
    "X_radiomic = X_radiomic.set_index('SubjectID')\n",
    "\n",
    "# apply imputing\n",
    "X_radiomic = X_radiomic.replace([np.inf, -np.inf], np.nan)\n",
    "imputer = SimpleImputer()\n",
    "imputer.fit(X_radiomic)\n",
    "X_radiomic.iloc[:, :] = imputer.transform(X_radiomic)\n",
    "\n",
    "# apply standardization\n",
    "scaling = StandardScaler()\n",
    "scaling.fit(X_radiomic)\n",
    "X_radiomic.iloc[:, :] = scaling.transform(X_radiomic)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N subjects: 192\n",
      "N radiomic features: 929\n"
     ]
    }
   ],
   "source": [
    "print(f'N subjects: {X_radiomic.shape[0]}')\n",
    "print(f'N radiomic features: {X_radiomic.shape[1]}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run PCA to reduce dimensionality for radiomic features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "# Apply PCA\n",
    "pca = PCA(n_components=0.9, random_state=random_state) # Choose the number of components you want to keep\n",
    "pca_data = pca.fit_transform(X_radiomic)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N principal components: 51\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print(f'N principal components: {pca_data.shape[1]}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run unsupervised KMeans clustering on PCs & visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/sklearn/cluster/_kmeans.py:870: FutureWarning:\n",
      "\n",
      "The default value of `n_init` will change from 10 to 'auto' in 1.4. Set the value of `n_init` explicitly to suppress the warning\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAABKK0lEQVR4nO3dd3hUZfbA8e+ZnkbvIB1BQUWK2LGAFXvvvevq6v7Wtru6ay9rWXtF3bVir1gAQVGkI00EQem9pE4/vz/uJCSZSSVhEjif58lD5t6Ze9+bhHvu284rqooxxhhTmivdBTDGGNPwWHAwxhiTxIKDMcaYJBYcjDHGJLHgYIwxJokFB2OMMUksOJidmoh8KyKXNoByHCQiC9Jchrkickgl+xvEz8psHxYcTBIR+V1EhpV6faaIbBKRoSLSVURURGaU+0wrEQmLyO/bvcCVEBGfiNwpIgtFpCBxbS+LSNc6PMeFIvL9thxDVb9T1d51VaZalqGvqn4LkPiZ/a+2xxKRQ0QkLiL5IpInIgtE5KJS+6v1exGRV0QkKiLta1sWUzsWHEylROQC4CngWFUdX2pXpoj0K/X6bGDJdi1c9bwLHI9TvqbAXsA04PB0Fqo0EfGkuwz1ZKWqZgNNgJuBF0Rk98S+Kn8vIpIFnAJsAc7djuU2AKpqX/ZV5gv4HRgGXAGsBwaV2tcVUOBvwEOltk8Fbgd+L7WtA/AesA4ncPyp1L59gB+BzcAq4EnAV2q/AlcCCxPveQqQxL6ewHicm8Z64O0KrmMYUATsUsm1fgtcmvj+TuB/Ka7Vk3h9IbAYyEtczznAbkAQiAH5wObEe/3Aw8BSYA3wLJCR2HcIsBznhrka+G/xtnK/g78APyeu820gUGr/XxM/t5XApYly9kxxfYcCs0u9/hqYUur1d8CJ5X7vRwFhIJK4plmlflZ3ARMTP4OvgFYV/FzLXE9i2zrg1Or8XhLvPx9YBlwPzEn3/4ud7ctqDqYiVwH/Ag5X1akp9v8POFNE3ImnwWzgp+KdIuICPgFmAR1xnghvEJEjE2+JAX8GWgH7JfZfXe4cI4DBwJ7A6UDxZ+/CuTE1BzoBT1RwDcOAyaq6rJrXXKHEU+x/gKNVNQfYH5ipqvNxgtiPqpqtqs0SH7kf2BXojxPMOgL/KHXIdkALoAtweQWnPR3nRt0N52dwYaIsRwE3Jq6vJ86NuCKTgF6JZj9v4jgdRCRHRDKAQTgBooSqjgbuxQm62aq6V6ndZwMXAW0AH04Aq5SIuETkJKAZMJvq/14uAN4E3gL6iMjAqs5l6o4FB1OR4Tg3ltkV7F8OLMD5j34+ztNvaYOB1qr6L1UNq+pi4AXgTABVnaaqk1Q1qqq/A88BQ8sd435V3ayqS4FxODdacJ5ouwAdVDWoqhW197fEebquK3Ggn4hkqOoqVZ2b6k0iIjg3/D+r6kZVzcO52Z5Z7lh3qGpIVYsqON9/VHWlqm7ECbT9E9tPB0aq6lxVLcSp8aSUOPYU4GBgIE6wnggcAOwLLFTVDdW49mIjVfXXxHHfKVWmVDqIyGac2t0dwHmquoBq/F5EpDNOrecNVV0DjMH5OzPbiQUHU5GrcJ58X0zc7FJ5Dedp9iySg0MXEjeH4i/gNqAtgIjsKiKfishqEcnFuXm2KneM1aW+L8SpnYDTpCLA5MQIm4srKN8GoE46MlW1ADgDp5awSkQ+E5E+Fby9NZAJTCt17aMT24utU9VgFaet6Po74DS3FKvqCXw8Tu3i4MT33+IE4qGJ1zVRUZlSWamqzVS1har2V9W3Etur83s5D5ivqjMTr18Hzk7Ufsx2YMHBVGQNTlPPQcDTFbznPeBYYHHi6b60ZcCSxM2h+CtHVY9J7H8G+AXopapNcAJHRUGoDFVdraqXqWoHnH6Rp0WkZ4q3fgPsIyKdqnNcoADnpl6sXbnzfqmqw3FubL/g1ITAae8vbT1Om3rfUtfeVJ3OWSr4TE2swmlOK7ZLFe8vHxzGU3VwqM90zdX5vZwPdE88PKwGHsF5eDimks+YOmTBwVRIVVfiBIijROTRFPsLgMNwOkTLmwzkicjNIpKR6JvoJyKDE/tzgFwgP/EEflV1yyUip5W6sWzCuZHFU5TvG5wO2A9EZKCIeBJt7VdWUNuYCRwsIp1FpClwa6lzthWRExJ9DyGcjtric64BOomIL3HeOE7geFRE2iQ+37FUf8u2ege4SER2E5FM4O9VvP8HoDfOIIDJieawLsAQYEIFn1kDdE30HdWpqn4vIrIf0CNR3v6Jr37AG1jT0nZjwcFUKlEjOAw4VUTuS7F/qqr+lmJ7DKdDuT/OyJ71wIs4wxbB6cg8G2fUyws4o3GqazDwk4jkAx8D1yf6NFI5Ffg8cfwtwBycTthvUpT568T7fsYZVvlpqd0unE7glcBGnKfu4oA2FpgLrBaR9YltNwOLgEmJZrNvcG7Q20xVv8DpHB9XfI7ErlAF7y8ApgNzVTWc2Pwj8Ieqrq3gNKMS/24Qkel1Ue5yKvu9XAB8pKqzE7XE1aq6GngcGCEiLeqhPKac4qGBxphGSkR2w7m5+lU1mu7ymB2D1RyMaYRE5CQR8YtIc+AB4BMLDKYuWXAwpnG6AlgL/IYzZ6TafTbGVIc1KxljjEliNQdjjDFJdoiEX61atdKuXbumuxjGGNOoTJs2bb2qtk61b4cIDl27dmXq1FTpf4wxxlRERP6oaJ81KxljjEliwcEYY0wSCw7GGGOSpDU4iEgzEXlXRH4Rkfkisp+ItBCRrxPLB36dmORjjDFmO0p3zeFxYLSq9sFZJnA+cAswRlV74eRwvyWN5TPGmLSLq/LzmtXMWLWSaDwpx2S9SNtopUTWy4NJrG6VSAgWFpET2Lqy1as4uedv3v4lNMaY9Pt5zWou//RDCsJhQPC6XTx+1LEc1LlrvZ43nTWHbjhryo4UkRki8mIiHXJbVS1eJWo1icVhyhORy0VkqohMXbdu3XYqsjHGbD+FkQjnfTCKtQUFFEQiFETCbA4GufLTj1hbkF+v505ncPAAA4BnVHVvnIVWyjQhqZPbI2V+D1V9XlUHqeqg1q1TzuEwxphG7evFi4inSHEUV+XDX+bX67nTGRyWA8tVtXhR+ndxgsUaEWkPkPi3onzzxhizQ9tUVJSyjyEUi7GhsLBez5224JBYvGOZiBQvgHI4MA9n8ZYLEtsuAD5KQ/GMMSbthnTahVRLuGd6vRzQuUu9njvd6TOuA15PLK+4GLgIJ2C9IyKXAH8Ap6exfMYYkza7tWrNUT168eVvCymKOst1ZHg89G/bngN35OCgqjNxlgYs7/DtXBRjjGmQHj7iaA75tRtvz51NJB7nlD67c/JufXGlqFHUpXTXHIwxxlTCJcLxvXfj+N67bd/zbtezGWOMaRQsOBhjjEliwcEYY0wSCw7GGGOSWHAwxhiTxIKDMcaYJBYcjDHGJLHgYIwxJokFB2OMMUksOBhjjEliwcEYY0wSCw7GGGOSWHAwxhiTxIKDMcaYJBYcjDHGJLHgYIwxJokFB2OMMUksOBhjjEliwcEYY0wSCw7GGNMIBaMRZq5exZLNm+rl+J56Oaoxxph6M2rubP45YRwuEaLxOD2at+CF406kXXZOnZ3Dag7GGNOIzFi1kjvHj6UwEiE/HCYYjfLL+nVc9NH7qGqdnceCgzHGNCIjZ04nGI2W2RZTZemWLfyyYX2dnceCgzHGNCJrCvJJVT/wuIQNhYV1dh4LDsYY04gc1rU7AXdyd3E4FmPPtm3r7DwWHIwxphE5a4+9aJWVid/tLtmW4fFw/ZD9aeIP1Nl5bLSSMcY0Ik38fj496zxemTmDrxcvokVGJhf3H8DQrt3q9DxSl73b6TJo0CCdOnVquothjDGNiohMU9VBqfZZs5IxxpgkFhyMMcYkseBgjDEmiQUHY4wxSdIeHETELSIzROTTxOtuIvKTiCwSkbdFxJfuMhpjzM4m7cEBuB6YX+r1A8CjqtoT2ARckpZSGWPMTiytwUFEOgHHAi8mXgtwGPBu4i2vAiempXDGGLMTS3fN4THgr0A88bolsFlVi7NKLQc6pvqgiFwuIlNFZOq6devqvaDGGLMzSVtwEJERwFpVnVabz6vq86o6SFUHtW7duo5LZ4wxO7d0ps84ADheRI4BAkAT4HGgmYh4ErWHTsCKNJbRGGN2SmmrOajqraraSVW7AmcCY1X1HGAccGribRcAH6WpiMYYs9NKd59DKjcDN4rIIpw+iJfSXB5jjNnpNIisrKr6LfBt4vvFwD7pLI8xxuzsGmLNwRhjTJpZcDDGGJPEgoMxxpgkFhyMMcYkseBgjDEmiQUHY4wxSSw4GGOMSWLBwRhjTBILDsYYY5JYcDDGGJPEgoMxxpgkFhyMMcYkseBgjDEmiQUHY4wxSSw4GGOMSWLBwRhjTJIGsdiPMTuK75f+wX3fj2fRxg20zcrmT0P249Td+6W7WMbUmAUHY+rIpOXLuPzTDwlGowAsz8vljm/HkB8Oc2H/AWkunTE1Y81KxtSRh374riQwFCuKRnnspx+IxeNpKlXjVxAOEyr3czX1z2oOxtSR3zZuTLk9GI2yJRSkRUbmdi5R4zZrzWpuHfMVCzdswCUwvEdP7j1sOE38gXQXbadgNQdj6kjnpk1Tbve63HZDq6GVebmc8/47/LJ+HTGNE4nH+fq3RVz40fvpLtpOw4KDMXXkxv0OJOApWxnP8Hi4ctBgPC77r1YTr8+eRSRWtikuEo+zYP165q5dk6ZS7VzsL9aYOnJI1278e/hRdMxpggBN/QH+NGQ/rh40JN1Fa3R+3bCBSDyWtN3tEpbmbklDiXY+1udgTB06uldvju7Vm3AshtflQkTSXaRGaWD7Dkxc9kdSB38kFqd3y1ZpKtXOxWoOxtQDn9ttgWEbnNF3DzK9Xlxs/RkG3B6GdulK9+Yt0liynYcFB2NMg9M8I4OPzjyXY3rtSo7PR+vMLK4cNJgnjh6R7qLtNKxZyRjTIHXMacJ/LBikjQWHOqTxjWjeIxD8EsQDGSch2dchkpHuohljTI1YcKgjqiF0w6kQWw1EQYGC/6LhqdDibWt/NsY0KtbnUFeCoyG+ESg9uiIE0V8hMjVdpTLGmFqx4FBHNDwbtDDFjhhEftn+BTLGmG1gwaGueLpCqr4F8YC703YvjjHGbAsLDnVEMk4AfEDpvgU3uJqD/+A0lcoYY2onbcFBRHYRkXEiMk9E5orI9YntLUTkaxFZmPi3ebrKWBPiykFavgnePQE34AHffkiLNxFxp7t4xhhTI+kcrRQFblLV6SKSA0wTka+BC4Exqnq/iNwC3ALcnMZyVpt4eiItR6HxAhCXDWE1xjRaaas5qOoqVZ2e+D4PmA90BE4AXk287VXgxLQUcBuIK6tMYNDIz8Q3nE189R7E1x5EvOBVVDWNJTTGmMo1iD4HEekK7A38BLRV1VWJXauBthV85nIRmSoiU9etW1en5VENE899kPiagcRX70Z8w1nEw7Nrd6zIr+iG8xLDWUMQXwN5j6B5D9dpmY0xpi6lPTiISDbwHnCDquaW3qfO43XKR2xVfV5VB6nqoNatW9fq3BpZiAbHorEVZbdvvgkK/weaB8QgMg02nkJ8/YlodFHNzpH/JBAqt7UICl9D4/m1KrcxxtS3tM6QFhEvTmB4XVWLl3haIyLtVXWViLQH1tb1eTWeh266HCJznaGmGkEDRyBNH3BmOIe+JfmGDkTnoRvOhNbfIK5m1TtZdB6QYv1g8UJsBbh61/5CjDGmnqRztJIALwHzVfWRUrs+Bi5IfH8B8FFdn1tz/w6R2UAQNB8IQfBrtOBliC0B8VXy4TBa9GH1T+buQdnhrcXHiYC7fc0Kbowx20k6m5UOAM4DDhORmYmvY4D7geEishAYlnhdZ1SdQADhcnuCTlOSuyto+X3l3hddWO3zSfbVgL/c1gBknIC4mlT7OMYYsz2lrVlJVb8n5SM1AIfX34nDVNCNAVqIeHZBfQdC+FsgeZlCyATPntU+nfj2guZPo7n/gthSkABknI3k3FiLwhtjzPax02VlFVcO6u4Csd/K7XGB/0A0XuiMKEoZt9zgykYyjqvZOf0HIq2/QjUMeBBJ+zgAY4yp1E55l5Km9yTyIHkTW/wgTZDsv6AFTyeajaLlPuWCwLFIy/cRV2aV51ANo9Hf0PjGrecVnwUGY0yjsNPVHADENwBafoIW/heiv4F3byTrHMTVItHZnGKkEi6kyR2IK6fK48cLR0He/UAMNIr6hyJNH0Bc2XV8JcYYUz92yuAAIJ7OSJPbU+xJMezU+QQV9lWUoqGJkHsXENy6MfQtuuUv0OwJCH2DhiaAqzWScSri6VyL0jduGl0K4YkgWeA/zIKmMQ3QThscKhQYAYWvA5FSGwU8fao1ukgLXqBMYADnWKHx6IbTIPZ7Yt0HD1rwCjR7DAkcVlelb/DiuQ9B4WuAgLiAf0Cz5xD/kHQXzRhTSpUN4CLSRER6pNhe/SE7jYhk/ymxNkNWYksmSFOk2UPVO0B0ZQU7Ys6qcCULAkWBILrl/1CNVPCZHYuGJjnDhQnhzDEpBC1EN1+V6Kw3xjQUldYcROR04DFgbWI284WqOiWx+xVgQL2WLg3ElQ0tP4LQODQyF3F3gsDRiCur6g8DeHpA+PcKdpbv5AZQiMwB397Je2LrIL4ePN0QCVT3EhosLXoXKEq9Mzxpp1n34tvfl/DKrOlsKiriiO49OX+vvcnxl58LY0x6VdWsdBswMJHKYh/gvyJyq6p+QMVzFBo9EQ8EhiOB4TX/cOBwCI+pwQfiIM6NQeP5zg009B1EF0N8bWK2dhzNvgFX1kU1L09DUlkNaSepPT01ZRJPT/mJoqjzoPDrhg28O38un5x1Htm+SmbmG7OdVRUc3MUZUlV1sogcCnwqIrtQnd7ZnYhzY/8YIjNxWuvKd2z7EtvK1R6kOXh2Q+Ob0fUnQXwDZfosim+a+Y+h7s5IoP7mB9Y3yTgWDY0Hyq21rVHw7fh9DpuDRTw5eRKh2NbJlaFYlDUF+bwzdzYX7z0wjaUzpqyq+hzySvc3JALFIThrLvStx3I1KhpdjK47FPLug+CHJAeGAHi6QeAUnDkVmU6fhrRAmj+HiDgd2fF1JHdmF5+kCM1/rl6vo975h4H/AOf6AefZxA9N7t4pRizNWr0anzt5VcBgNMqYJeUnZRqTXlXVHK6iXPORquaJyFHA6fVWqkZGt9wGmkvqypQ48yhaPItIBhq9DMJTEmtLH4gUJ/lLme+pnOivdVzy7UvEBc2ehPAkNDTWmXiYccJOM5y3RWYmsRSLPAmQ6fVyzeefMGXFclpmZnLloH04ftc+OPkpzY5MVfl14wYE6NWiZYP5nVcVHApwFtspv4jBPsCkeilRA6TxPLToIwj/BK5mkHk2Lu9uzj4tgsgsKm5lU4gtKVkZTjydIdXNUKrz5BxEY+sQd+3Wr2gIRAT8+yH+/erl+BrPc9KfuDs2uGVa+7VuQ/vsHH7fvKlMkPC53fywbCmhWIy4KuuLCrltzNcsz93CNYP3TWOJTX2btWY1V3/2MVtCTotBM3+Ap489nj3btktzyapuVnoMyE2xPTexb4en0T/QdYdB3l0Q+hKK3oYNJxDfeEliCGo10mGIt+q3ZF0AVHEzk4waZYTdmahGiG/5G7p2P3TDaeiafYnnP9mglmMVEV458RR6tWhJhsdDts9HptfLnm3aEU4EhmJF0QhPT/mJwsjO0VG/M8oNhTjv/VGsys+jMBKhMBJhZX4e534wirxQqiwN21dVNYe2qpq0Pqaqzk4s7bnD09y/gW5J3hGeiOY/i2RfA7hJncEVnPTcp5U9ZjwXzbsPgp87nbGuNuDaBbx7QmQGKTuuwXmvu9O2XdAOSvMegqKPgfDWlOv5L6DSGsk6I61lK61jThM+P+cCftu4gS2hELu3bs2IN/6bsrnJ7XLx++ZN7N66TRpKaurb5wsXpPy9x+PK54t+5Yy+e6ShVFtV9djbrJJ9DavOXg9UY07/QEpxKHor0aSU3MnocINvAJJ1caljxtGN5zg3Mi0CIhBfAdFJzrHcHSH7JpLXgPCBb/BO0z5fE6pRKHyL5M78Iih8Ph1FqlKPFi0Z0L4DAY+XXZo2TfmecCxG26wdv6N+Z7WusIBgNLlmGIxGWV9YkIYSlVVVcJgqIpeV3ygilwLT6qdIDYlQ6Y9Ii5yV5KSC4ODZHWk+cmunM0D4R4gtp2x6jmJBiK0CcSHNXwR3F5zKnQ8CRyHNnqj1laSbxvMTWWoLq35zjQ8eJPUEQxJDgxu2qwYNIcNTthLvd7s5vFt3WmZWnQHYNE6DO3Qiw5vc5BzwehjcIf0tBFU1K90AfCAi57A1GAzCGbR/Uj2Wq0EQcaH+YRAanfoNnr5oeBZoqijvhozTk0YeaHRh4mZWkSAERyMtL4JWX4HmgQTKBphaUI06TVYadWozsn1m5KrG0Lx7oHBUYr3uGJp1IZL957oblSFZ4GoN8VXJ+7zprZpXxz4dO/HgsKO4c/xY8sNhFOWYXr25+9Bh6S6aqUdDOnZiQLsOTF21gmBiUmSGx8PA9h0Y3KFjmksHUp0Ou8Tkt36Jl3NVdWy9lqqGBg0apFOnTq2XY2t8E7rhVIgtK7fH7cxe1jCp+xsEWryFK5EWQ+P5iXbxt6k482uC7xBcLbY2h6iqM1qq4DknnYa3H2RdgfiGVOsGq+GZ6KYrcNrjY0AEvEOQnD8jvr2c5rPgp2jhe4AimadA4DikohpRDcTzHoXCVxJNaAmSAdk34co6f5uPX0yDY9HNN7C1ackF+JGWbyDexjElJ67K2oJ8mvgDZKZ4ojQ7nkgsxttzZ/PO3NmICKfv3o/T++6BN8V8mPogItNUdVDKfZUFB3ES+lwJ9ARmAy+pagX19/Spz+AAiX6C4BdQ9F7iJueFyHQqn5cgkHUNrpw/OTf3DadCdD4VNn+UfCwDafYEUirPUDz/eSh4quwNFsDVDmn6IOKveLijxgvRdQc6zV9J/JD9Z4hMg/D3pY6fAf4DkGZPbdPTvaqiawekrlm52uJq812tj53yfOFpaP7TEFsCnn5I9rWId9c6PYcxO5LKgkNVzUqv4jSOfwccDeyG09S0UxFxIRnHQsaxAMTXHkKVE9bwbB1nH5maWJa0ssDgAryQdVmZwKAaTB0YAOKrnRpBqw8RT7eUR9WCZyto9gIIQf6/cTrUSzd1FUFootMM5RvgZEyNLgJXC8Rdk/HX0VJZaMuXfVMNjlOWxlZD6HunBuI/pCQpovgGIi1eqvVxjTFbVRUcdlfVPQBE5CVgcv0XqTGoolkIcJYVPdr5NroItLLPuMF/KNL0bsTVouyu2AoqHzcQRgteQ5rekbRHo4uh4EWqToOVakx1CMKTiEcXJla1w1nVzrsX0vwJxNW8imOCiBd1d4bYH8k7vbtX+flU4vnPQP7TibUgXIBCs2fqbVKdMdtbMBrBLa7t1rRUkaqCQ8mQGlWNNpRp3WmXcTwUjCT1iKNEDaDJPxGPM+JA41tIfQMu5kVy/pocGMDpaK00Y2kMYkucyV7hH5ysrhpGMkagwUlUHcjiiTKX7zfxo/FNkP8cZdJsR6ajm65BWr5RxXEd0uQf6KZr2FozEcCP5NxWrc+XpuFZkP8MECoT73Tz1dDmxx0irbnZec1du4ZbxnzFL+vX4RLhyB69uPuwYTTxp+fvuqrgsJeIFM+QFiAj8VoAVdWql0bbAUnWlWjR+07ncBI3ZF8Nhf8lnv+QMxw1MocKn94lG2n6b8TTNfVuVxM04zgo+pTUAcYPvsFo3oNQ9EZJ85OGJwIBqg4OFUzeE3HW105afyEKkdlodBni2aWKY4P4D4IWr6L5TyXW694t0RdQ85qDFr1H6uY8cZqZAja6x6TX8twtbCgqonfLlgQ81R9UsCY/nzPfe4eCiPP3HVPly98Wsix3C++ffnZa8i1VGhxUNb31mgZKXNmopz+Ev0m1F/KfpKRWUdU4e/WDf2jl52vyTxS/M+muzM3e7Qzj9B0EG8+hTPDQwkQZKpu9nYoPJAtp/iS65Z8VFMibuK6qgwOA+PZGWrxYgzJUQEOkDnaa2GdMemwsKuTKzz5m9po1eN0uYnHl1gMP5tw9+1fr86/PnkkkXvb/aSQeZ+GGDcxZt5Y92rSth1JXrhqJgUxKvr4kz2IG58m2JvlwNqHhyucTivhwNb0TWo0F776J8wbAfxjS6kMkOo/Uay9FqFlgAHBD6/FOp2+sgjTSGgNPxaOAYrEYY9/8ntuOvZc7TnqQHz+ZWic5jiRwVKl036XLEwX//tt8fGNq68rPPmbm6lWEYlHyw2GKohHu+348E5el6G9LYeHGDYRjyf9XXSL8sbn2gze2RVXNSqYCknEmWjAyMc+h+MbnxbkZV6fDulgcQhPAP7jSd2lsBWw4OVEjCAEeCE+E2AqUai5hWg2qQaZ//g6fP/cG/kAHbnhoOV6/srVWG4CcGxBX6pm7qsodJz3ErHFzCBY4T/PTv/mZ4ecP5U9PJU22L/kM4R/Q4KeA20nj7Uvx8/APBd/BEBqP09zlBryQc1u1OsiNKRaKRvlpxXJUlSGdOtWoCai8ZVu2MHvNGqLxsv/vi6JRXpw+jQN26VLlMfq3bc/4338nGCs7ojGqcfq0Sk8WZgsOtSTultDyHTT3X04qb/GDf0RisZ+qhrmWE/2lyrdo7kOJBIDFf4BRZ/TQlttA2lPhIkE1FI0qs8e8zPefOX+QcyZnceU/V7DfkXnOegyZZ1W6XOmMsXOY9e3cksAAECwI8eUr33LidcfQuU/ZmZ+qiubeDsHPEv0lghZ9gmaeg6vJX8u811kP4vFEIPkaXFlIxkmIp2edXLvZOXy39Heu+fyTkteqyuNHjeCwbt1rdbwNRYV43S5CKSrpawpSzS9Kdka/PXhhxlTCRTHiiYdNv9vD/p0607NFy1qVa1tZs9I2EE93XC1ewdVuPq62M3E1u9sZyUT50QU+kFYVH8jlJF7T+Ga06GPn5hjfmglWY2sgNIaUNZLY7xD9obJSVrA9dXeS1wsHH7ex5PWaZX4euLYLsyZmA17EXXk/w5TRMwjmpwhUqkz/5mdGvzyWy/vfxDldr+KJa18kb81EKPqMrfM4FCdh3n+dobjlr0YE8R+Aq+mduHL+zwKDqZFNRUVc+elH5IfDJV8FkQjXfvEJ6wpql+yud8tWxOLJzaY+t5tDunSt1jGaBTL46IxzOapXL7K8PlplZHLZgIE8fezxtSpTXbCaQx2TJneiuKDoI0CcGkXOzUjGSeiagSStn4wfyTyDeOFHkHsbTrOU84emWdeBuw3k3kXFQ2GrastPNUy1eLuSKuBEwmUDSrDQxeSxOfQ/cCOUmqCXSk6LbDw+D9Fw2eqx2+Pmx0+mMu+HBSW1is9e+IZ2bd7m5MuDKUKYOs1tnto9zRmTyueLUq+mqKp8unABF/UfUONjZni9/PWAA3lw4ncUJXIk+VxumvoDXLJ3ysnHKXVs0oQnjz6uxuevLxYc6piID2l6N5qTWDrU1QoRDxqeTkWdw+pqA7kXkDSDuuA/OL+ibclYUlGHdASnViGUDjDBQpj4RRPa7hJizTKnw93rU3KaCWSeX2XK8GHnHswbd7+XVGIFfh4/r0zQiEVi5G2KEo+5cLvLl9PtzIA2pg7lh0NE4ikeiGKxbVpg54K9BtC9eQtenD6VdQUFDO3ajUv3HtSos+pWK/FeQ1ffuZW2lWoU3Xyjs5Jc+Sd9yQbfvhBKNSx2exDnJqwxVMOEEq07IrB6qY+7LuvKmuV+Xpp5Oe12PaLMJ1XDzpBWVwtE/AQLQ/w8fh6/TF7Iu//+BHGJc7kuoXOfjvw69TfisbL/Mdt0DPPyxF/x+soHhwDSZrx1NJs6NW/dWk4d9WZJFtRiGR4Pb5x8Onu1a5+mkqXHtuRWMtsoXjgK8h5I5DdKFYgTzSfp5O4NsWUI6wmUetDZpVeIf3+wiDlTmtC2V+l8T4oWvAAFzySyvArLVw7jmkNWI65EX4bAWbeehNvj5tV/vM2vUxcRjyVf//rVAcZ8ciRHnTLGSemNADGk6WMWGEyd2711G07YdTc++fUXChML7WR6vRzRvWeDWLe5IbGaQz3S0LfopuuoKnVGcuK77S0Dp+kqeX6GaiLxYJspiCsHgHjB25B3L6VnTwcLXbz9ZGveeGzrfzBfho+cFtlsWLGx/GFL+DN9PDX5fjr3aeZkhsUFvgNLkumZhmvWmtWMXvgrbpeL43r3oXfLSgZdNCCqyrjfl/DuvDkoyim79eXwbj3SMgs53WqdsjudROQo4HGcO+eLqnp/Re9tqMEhvv6ERJruVBJj9DNOg+AHFaTUbiCkNV9++i/ef/wzCjYX8vRXU2jSLC/pbQW5Lk7u04/iEVL+DB/RaIxYJHW/R4v2zbjpxavZ5+i967P0ph7cM+Fb3pgzi2A0iksEr9vN9UP244qB+6S7aDuFwkiEGatXku31sWfbdrUObI2uWUmcVWaeAoYDy4EpIvKxqs5Lb8lqKLqk4n2BY5Csy8HdES16uwYH9VKzGdjbbuI3/Xj6+pEEC50aUCAjdSDLyIrj9WnJaKdoNIamGOIH0L57W1759T+4XNs2mjoSjvD1q+P55vUJ+AM+jr1iOAecuM9O+RS4vcxZu4Y35swqGZkTUyUWjfLYpB8Y0asPHZvslCnXtptRc2dz5/ixeFwu4qo0C2Qw8oST63w+REOd57APsEhVF6tqGHgLOCHNZao5qST2ZpwNnu7o5j9T9XDU0soHBh/4Dnc6tuuBKuwxYAK+wNYx4Evmp84SuWaFl0h465+Ux+Om7wF98HjLzqkIZPo57S/HbXNgiMVi3HzE3Tz951eYPWE+U7+axQPnP8ET19ZBHidToS9/W0gomnoE3ZglFaRcMXVi7to13DF+LEXRKHmJORor83I574NRxFKMwtoWDTU4dARKr8u5PLGthIhcLiJTRWTqunXrtmvhqs07sIIdbsTTFS14zZldnXTDb4rTD1AZFwTORtr+jDR/JNExvG2iETflWxlFwB+Ic/TZWxMIPn9nB4KFQum/xUjYw8v3bk0TEMjyc9g5B3H7mzfQc0B3fAEvmU0y8Pq9DDvvYI69fPg2l/enT6ezaPpiQoXlZmOPHMfyhSnWkzZ1wud2pwzsxc1L29OUlcu5+ZsvufHLzxm7ZHGd5PBqyF6fPSspB5MC+eEIU1auqNNzNdTgUCVVfV5VB6nqoNat05N7pCqS8xeSk/N5nfWZ3S2h6B1Sd0QXQvY14GpLxTOc4xB8Cw2Odlaca/IvnJnZtfuVxmIZfPxK6p+jP0PZtf/WyXtzJmdz65m9WLa4M7hagncA3rYvcvjF/+SgU4awz9F7065bG74cOY6zOl3BwmmLUSBUGGKvobtz6QPnbnOtAWDKlzMoSjEbW1wufv527jYf36R2bK/eeFL8/lRhePftN2P93z9+z4Ufvse78+bw4YL5/OmLT7n+y8926ACxsaiIeIrrE2BzsG4HtTTU4LCCsvmgOyW2NSri7Y20eBk8vXDmE2RC5nlI03sS76hocpuzLKmrzXeQdXUlZ4jDlluIR5dCZCZIjvPl2QOy/gTuPYEMkIrbgKMRWPWHj0VzW9Bj99TpA0JBYcm8sjWZRXObIM1H4mrzI66Wb+Hy78/+xw/m/0Zew4Kpv/HHvOXEY3FUlVg0RiQYIRaNM2v8XO49+7FKrqlqy39dyat3vM1vs/7A7UnxBOt2kdMyZ5vOYSrWvXkLbj7gIPxuNxkeD5keL363mweGHUmr7TTpa+mWzbw4fSpF0WhJo2xhNMLYxYuZvGL5dilDOgzr3oOMFEkCI/EYgzp0TPGJ2muQHdLAFKCXiHTDCQpnAment0i1I77BSKvPUI0C7rIdpYERUPASSYn63O3A1cH5Pvh5FSdww4ZTEvMoEsEmuhDiayG+xVl7QXNTfvTjkS0ZeX974jGIxVzstb+fPfbLo3QR43GIhoXP/lu2s2vEFUfQtW/ZPEuTv5jB41c9z5Z1qc8HEAlFmTl2DutXbqRVhxQr31Xhsxe+5ukbXiEWiRGLxVJ213g8LoYcs3OOgAoHw0x4dxKLZixml94dOfSsA8nMqfuZ5hfsNYAje/Ri3O9LcItweLce23U28IQ/fk856KAoGuGbJb8xpFP11hppbI7btQ+vzprBb5s2lkzky/B4uWrQ4DoPzA0yOCSWJL0W+BJnzOfLqtqo2wkkRee0ZF2GhsZAdDlOzqUAiBtp9u+tf/hafiW2cjSCM4+idC0kCPFEFVNTVzW//7wpL97dnlDR1jbi6ROyefL2DvzpvpXOR+Mwf1omj/+1E5vWlX1aKX/DGfPGdzx04ZPEolV3inm8HhZNX0yTljn4/NVPlbx53Raevn4k4WDyaC1/pg+X20VmTiZ3f3ILvoCv2sfdUWxet4XrhtzGlvW5FOUHCWT5eflvb/L4xHvo1KvuZ/62y87hrH57VvoeVWXKyhXMWL2Stlk5HNmjJxne2qfHLpbp9eJKERzcLhfZvh33d+/3eHj3tLN4d/5cPlu4gCZ+P+ft2b9aacFrqsHOc6iJhjrPoTpUIxD8Gg1PBXcnJPPEMmtJx3PvgsK3qHj4qouarR/huHJYL5bMS37S8HjjvD5tHs1axZgyLou7L+tGsLBsJ6PH5+Efo25iv+MGJa5BObnlReRvrn5WS3+mD3G5OPPmEzn7tpNZOH0xL/z1v/w6bTHN2jTl7FtP5ogLDynzdPjJs1/x9A0jk5L6icC+xw3inNtPodfA7nXSn9EYPXTRU4x547sy80rEJfQ7oA+PjP/Xdi9POBbjko/fZ8bqVYSiUQIeD16Xm7dOPYNdt3HCXG4oxP4vPVcyy7mY3+3hi3POp2szm11fHZXNc9g5/xc1ICJeJOMYXE3/gSv74jKBAUCyrwFXa5JHL/mdbb6hVJR+uyLRCDz91ULenTuHs29Yjcu19QHB5fEza8a5/Dq3H8sWZhOPJT+dRcNRFk7fOmQxf3MBBVvKZ5utXKgwTDA/yFv3fcCrd7zNTUPvYOa4uRTmFrFy0WqeuO4l3rzv/ZL3r126jhdu/m9SYAAnjXfbLq3pPbjnThsYACZ+ODlpwqHGlbk/LCAcrOEaI3XgtVkzmLZqJYWRCDFVCiIRtoSCZdZSqK0mfj/PjjiBLK+XbJ+PLK+PgNvDPYcNs8BQRxpks5LZSlwtoNVnaNGHzrBXVxNwtUbcHSBwBMTXo+snUTqVRVWK+7Nymsc449q1NGke49k7nM4sr8/LwWffwvIFK3n1oVsJV5Cp8vW738ft8XD2bSc7TUxlk7uW0aZLa1SVDcs3EC83KS5YGGLUvz8hUq6pKFQY4rU7R7Fo1h8cctp+THh3Uup1IgBvwMuw88quwx0JR3jv0c/4/IVviIQjDD1tf8752ynkNK+f+SANgcudOjCKiJMEcTsbNW9OUoI7BVbk5bI8dwudmjTdpuMf2LkLUy67iu+X/kEkHueAXbrQxJ9q6V5TGxYcGgFxZSFZ50DWOWh0CVr0ORr63lkkByBwDIQnQHwTNU3vHchUjjlvA6891I7CfDeHnXUgbrebzrt1olnrpqwuWJvyc/FYnFfveJtNqzdz7ROXsOfQvswaNyfpfV6fh5fmPorL7WJE1jkpjxUuSv1UG4vG+G7Uj0z5YgaRYDhpDkaxU248jt6DepTZdudJDzHr27mEEsf+6KnRTPp0Gs/PerjW/RGqyvxJv7LytzX06N+Vbv0qT1++vR1+zkF89vw3REJbA63b42LgEXvi9W17O39NpRpyWZ19NRHweBm2HYfP7kx23jp4IxQveBldfzwUPA6hLyDyg/MVfM9JnV1L0YjQumMYj9fNcVcfCThPm3e8/xe8voqfHzSujB45lkUzl3Dn+38hq1lyH0b3/l2Y9+MCRj38MTktavfUHswPVtjR7fa6Oee2k8tsWzh9MbPGzysJDOA0hW1ctYnxo36sVRlyN+Zx1cC/cvORd/Ofq1/gun1v5bZj7yUc2r6pTCpz8T1n0W2PXcjIDuD1e8jIyaD1Lq248YWr0lKek/rsjj/FpLi2Wdnsso21BlP/LDikkcbWo0UfocEv0SpGJWl0GeQ9SkUZXlXjRMK1myXt8SjrVviIRmJM/XJmyfae/bvRtE3leXIioSg/fTad7KZZdOvXOan54tepi/nbiPt47Y63U05Y2xYen4f23dpwTterOaX1xTx25fPkbsxjwZTfSNXGVZQfZO7EqtfrTuXRy57jj7nLCOYHKcoPEioMM2vcHN64+91tvIq6k5GdwZM/3c+/PrqZS+8/l9vfuJ5Xfv0PLdunpw3+4r0HsHvrNmQmRidleDxk+3w8cfQIy33VCFizUprEC16BvIcBDyUTC5o9i/iHpP5AaAyV5WAKFQm1+f8WLBQ+/19LCvOdJ7yXb3uD/Y4fTPtubbjrjEdYv7zidNvgNLWsXbqezeu2sGDKoqREexpXIiGnqSue6Ffw+r1lmj5qyuv34HK7EYE1f6wrOf6XI8cyY+xsrnj4fNwpnlh9AS/te7St8fnCoQiTPp1KtFxnbzgY4fMXx3DhXWfV7kLKWbtsPU9fP5Ipo2fg9ro5/OyDuOzB82o0T0FE6H9oP/of2q9OyrQtAh4vo047iwl//M70VStom53Dcbv2pok/dW4u07BYcEgDjcyDvEdwJr+FS+75uvlKaPODkw4jiYuKUmnE4xAOQfX/zwkaV3I3u3nv2da881Sbkj3RSIzxb09EEX76bFrV1xJXvn5tPPF4vNr5A8UluNyupFXhqsuf6eeyB87l2RtfLQkMxWXftHozhXlFZLfIIlgYKnMOj9fDERccWuPzxaKxClMypJpzURuFeUVcu88tbFmf55Q5GOHLV8axcMYSnvjx3kb7pO0S4ZCu3Tika7d0F8XUkDUrpYEWvU/SrOhiFa0KF0idqC4ShisO3ZXT99iDOy/qSrBQKMoXQkVSYQduMJjBzaf34PR+/Xj7ybaolr3xzP7+F964512i1WymioQijH5pLJEUw0xTCReFtyn/TWFuEcH8ENFI8vmK8oMsmrGER8f/i9327YXH58Hr97JL7w48OOYOmrepeVt3RlaAbnsmTzJyuV0MObbmC9KnMuZ/EyjKC5YJZpFQlD/mLmfej7/WyTmMqQmrOaSDFpFy4ppqhTOaxd0ebfI3yL0biKEaJRIWRt7flqULnZrG9AlNOHtAXw44eguZ2TGGntyU3fdeQunkfsFCF6893JTZkyqeaj/r27kp5xPUVGW1g4rWeaiODj3a0q57Gzw+T5maAzi1ik69OtCmc2se++5ucjfkEQlHt7nd/aYXr+KmoXcQjUQJByP4M31k5mRw6f3nbtNxi/06bXHJehmlqSq/z1lK3/17AzDn+/mM+vcnrFu2gYFH7MnJN4yoVcAzpioWHNJAAkegwc9Ay08ci4HvwAo/58o8A/UPheBXrF78B387dTrLFpb9FRbkuvnq7RYEMv0cfN7/QeBNCH5GJCLEYyE+erkl7z3bmq0TE1JPcqsLXfp2Ytn8FUlt9alkNskgGo4Rj8edGkElsSNYEGTIMQPIaZ5NqDBcEoBEBF/Ay+Cj+zP65bEU5hUx8Ii96LJbp22+lp79uzFyweN8/sI3/DFvOX2G9OKoiw4lq2ndLGfafc/O+DN9hArL1ihdLqFTbyfP1pevjOOJa18sec/vc5cy+uVxPD/rYZq3bVYn5TCmmKXPSAPVOLr5egh/lwgQLsAHOTfiyrqwWsf46tVvefK6l1KnrHYLJ19/LFc+fIFzvvgm7jn9ZqZ8tb6k47lUaRBxRizF4i7cnm3rLC6t824dWf37ugrnMZTWon1znpn2ABM/mEw0EiOnRTYPX/w0sWhyYHG5XTw/62Gymmby4IVP8fMEZ4HAXQd059grhvPkdS8BEIvGEZdwzKWHc/VjFzXodvv8zQVc0Ota8jYVlNSqPD4Pnft05NkZDxGNRDm1zSUU5pYd1ebxujn+mqO46pEL01Bq09g1umVCd3QiLmj2OIS/Q4NfgmQiGScj3t2rfYxue3ROmm0Mzs3ilBtHcOl9W5s7xNWcBTO8KQKDw+2NE4+78Hg9ddbB6va6ycjJqFZg8Aa8DD//YFq0a85xVznzLCZ/MaPCfgkRYcPKTXTZfRce/PofJR3PHq+b09pdSrCgbPPM6JfHMuTYgQw6Yq9tv7B6kt0si//8eC+PX/UCs76di8vt4uBT9+XaJy5BRFj2y8qUP49oJMaU0TO56pE0FNrs0Cw4pImIC/xDEf/Qqt+cQq8B3ek9uAe//LSw5IbucgnZzbI465aTkt6//wmD+fCJL1L0AQjRsDMSKlbB0o/F/Jl+MrL9bF5bcUruYrFIjJWLVuPL8CUHCAF/ho9YNIbX56VL31049++nlXnLHgfv5oyASkEEeu69dfRLINNJmTBl9IyU7w8WhPjylXE1Dg5/zF/OK39/i3k//krLDs05+7aTOfCkCoYa14GOPdvz4Nf/IB6POykvStV0mrTMrnCAgPU5mPpgo5UasXs+u40RVxxBdrMs/Bk+9j9hH56acn/KdvCzbj0Jf0ZFaSOqbm4JZAfod2AfbnzhqkpnTZdWkFsIqmXmX3i8brrs1ombX72OS+49h7s+voX//HBPyQ2+2Myxcyo8TyA7kPIpurKhsfEUzVOVWbZgBdcNuZWJH05m46pNLJy2mPvPe4IPn6xifY064HK5kprAWnVsSd/9e6dcj/vUm46r9zKZnY/VHBqxQKafqx69kKsevbDK9zZr3ZSL7z2L5/6SOrNpef4MH5fcfw65G/LJ35RP611a0axVE77/4KdqD1mNR+N03qMLHp+HX6f+hriEISMGcubNJ+H2uCjILeSnz6ezfsVGDjplSJmcR1+9+m3SSKRihVsKeeLal/jbW38us33PQ/qmTLMRyApw+DkHV6vMxV678x1ChaEyo6pChSFG3v4Wx14+PC25iv7+zo3ccdKD/DptMR6vm1gkxrn/OLUkdboxdcmCw07kyAsP5fW73mPLhrySm55IYmBQuQfxrKaZnHDNUWxZn8f/HXYna5euR4FIMIy4pNpDUdt3b8s/Rt1EqCjE0l9Wcs+Zj3LToXc4TU2JAVMZ2QFevv0Nnph0Ly3aOUNOK60FxJWJH05GVcs8YWdkBbj51Wt54PwniMfiRCJRApl+hhw7gH2PG1j9HxQw94cFKft04qqs+WN9vSyeU5UmLXN4dMJdrFq8ho2rN9Ntj871ssqbMWDBYaeSkZ3Bf368l0cue4afJ8xHRNj78H6sWryWDSs3EiwI4fF58Hjc3PK/P+FyuXjooqdYvnBV0joB1eH2ujnxuqNLXt88/F/kbczf+obEvbcoP0g4GOHZG1/ltjduIByK0Gefnkz7alaZ5HmlFa9P/fOEefw8fh7NWjdh6Bn7c9Ap+9J7n56Me/N78jcXMuSYvel7QJ8aj1Rq26U165YlJzOMRWI0ryLfVH1r370t7bvXPA2IMTVhQ1l3UuFQBJdLEiOUwox/50emf/Mzbbu25uhLDqdtl9YU5RdxcsuLUs9TkK1LhaoqwcIQGiv7t3TBv87g3L+dCsC4tyby6BXPUpRXcfI9f4aPG1+8kseufB6gzByG0lxuFwOH7wkIs7+bR6gwjC/Dh8sl3PvF7fQ7oE8tfypbTRk9g3+e+m9CpSam+TJ8HHzaftz8yrXbfHxjGgIbymqSlF672RfwMfz8oQw/v+zIqcomr+W0yOb2N27A4/Wwac1mHr746TJP+R6vmznfzS95vWn15mr1dTxy2bNlJ4IlHvjdXjexcIxAlp+M7Az2HNqX/931bsnNu/jff536b95a8VzJinCxaIzJX8xg7dL19B7cg96De1arFjH4qL259omLee4vrxEJRYjHlUPPPIA/PXVplZ+tTCwWS9nhbExDY8HBlBGLxnjrgQ/56KnRFOUW4fV7k4KE2+PmgBMGM3C4MzT0uv1uS2r+iUZi/DxhPpvWbKZ522b0PaA3LrebihYj8vg8tO7cipWLVpfdkeiTOPFPxxAuCtG1b2eGnrE/txx5V5mn+mLBgiC/zfydXgO6s+aPddxw0N8p2FJALBLH5Rb67t+buz65pVodykdddBjDzxvKhpUbyWmRTUZ27dv3l/6ygseueI45E3/B7XEz9LT9uO7JS+pshrUxdc2Cgynj4Uuf4bt3fyx5ei9+wvX6nTxGgSw/2c2yuPies0s+k7s+9bwHj9dN3qYCmrdtRu/BPRkwbA+mfzO77E1dnFFX7bq1oVXHlixfsDLpOOISeg/qwQEn7rN1WyVP3sX77jn7MTau3FimY3n2978w6t+fcPatJ1f08TLcHjdtOreu1nsrsmV9LtfvfzsFWwpQddKTTBj1I8sWrOTJn+6zWoRpkCw4mBLrlm9g/Ns/lEmfoap4fG52HdSDtl3b0He/3gw77+Ayo2QGHdmfT5/7Oql/IBqN0bFnu5LXd7z7Fz5/4Rs+f3EM0UiU3oN60ql3B3r078rA4Xvy+QtjmPPd/KQEdJFwlN3327XMtqMuOozfZv6eNBva7XGzZuk6mrZuwsJpi5NGHIWLwox+aWy1g0NdGP3yWMKhSJksuZFwlKXzl/PL5EXsNqTXdiuLMdVlwcGU+GPecnyB5NxK0XAMEeHW//4p5edO/78T+OSZL5N3xJX5Py0s6SB2e9wcd9WRJSkyyht+/sF8+MTnrFqytmRWdSDLz8k3HJuUWG74+UP54eMpzBgz2ymvCLFIjEg4woMXPAlQYfqNukosWF2/zfojZRoREWH5gpUWHEyDZDOkTYkOPdqmnODm9rjo1q9zhZ9b+ssKAlnJKw1FIzG+eHFMtc/vz/DzxKR7ueCfZ9B7cE8GHrEXt7/5Zy5KsdKa2+Pmnx/8lfu//DtHXzKspGkmVBimMLeIwtyilCOdvH4PB5++X7XLVBd6D+qOPzN5dno8rnTtt8t2LYsx1WXBwZTo0KMdex3SF1+gbGet1+/l5D+PqPBz4aJw0trR4Dy5FxXUbN3ojOwMTv/L8Tz5033cP/pv7Dui4slrIk4HczQaTTkPw5/hw5/pw5dIG5KRHaBd1zac9/dTa1SmbXXkRYcRyAzgKvUz8gW89B7cg14Dum/XshhTXdasZMr4x6ibePr6kYx5fQLRSIwuu3fi+mcur3RG8F6H9E2ZFC6Q5Wfoafsnbc/blM9HT37BlNEzadmxBafcMKJkMZvaKMoPpmxCcrlcXP34ReRtymfV4jX03b8PB50yZLunvshulsWTk+/jmT+/wtQvZ+L1eznigkO4+N6zq/6wMWlik+BMSrFojGgkij/DX/Wbgc9f/Ianrh9JNBwlHosTyPazx4G7cdcnt+B2b00Wl7shjysH/B9b1uUSDkYQAV+Gn2ufuJijLjqsVmWd8O6PPHTR0wTL1VJ8AS9vLn+OJi1yanVcY3Z0NgnO1Jjb48btSb3+QyrHXDqM3fbdldEvjyV/UwEHnjyEfY7Zu0xgAHjv0U/ZvHZLSVI9VWcC29PXj+Swsw4sk3yvug44aR9GvzyW2d//QjA/iMsleP1eLn/ovAYVGGKxGD99Np2JH0wmq2kmR118GN1TrE1tTENgwcHUmeJ1mmeNn8uvU39j/YqNHHv5sJLZygCTPp2WMtuquIQlc5bRe1CPGp/X7XZz1ye38NOn0/nu/UlkN83iyIsPpWf/blV/eDuJxWL8/fgHmP3dPIL5IVxuF5+/8A1XPnohIy4fnu7iGZPEgoOpE0X5RVw96GY2rNpUMlT0ub+8xi+TF/J/L19T8r6mrVMnrYtGYjRpkV3r87vdbvY/YTD7nzC41seoTz98OIXZE+aXzMuIx+KEisI8c8NIhp62HznNa3/txtQHG61k6sRXr37L5nW5ZeYQhApDfPvWRFYtXlOy7ZQ/jyCQVbYfw+1x0X3PLjt0ptHx7/yQ1CcCTtqQmePmpqFExlTOgoOpEzPHzU2Z68jtdbNgyqKS10OOGcC5fz8VX8BHVtNM/Jl+uu/ZhX9+8H/bs7jbXUZ2oMI0GYEUcyCMSTdrVjJ1on33Nni87uRMrgqtOrUss+mMv57IiCuGs2jG7zRr25Quu3XajiWtuVBRiE+f+5pxb03En+Hj+KuO5ODT9qtRTqSjLjmccW//kBRAXW4X/Q/rV9dFNmabpaXmICIPicgvIvKziHwgIs1K7btVRBaJyAIRSZ1nwTQ4x115JG5v2WcNl9tFq04tUs5hyGqaxV6H9G3wgSEaiXLj0DsYefubLJi8iJ/Hz+PhS57miWtfrNFx+u7fm3P+dgregJeM7ACZORlkNc3knk9vTcuSo8ZUJS3zHETkCGCsqkZF5AEAVb1ZRHYH3gT2AToA3wC7qmqly5DZPIeGYcbY2Tx04VPkbswjHovTe59e3P7mDbTq0CLdRau1cW9N5JHLnklK8OcLeHlh9iN06NGugk+mtmHVJmaMmU1mTgaDjtyrVkN3jakrDW6eg6p+VerlJKA4n8EJwFuqGgKWiMginEDx43YuoqmFvQ/bg9f/eIbVv68lkOlPSpbXGE39cmZSYACnVjT7u/k1Dg4t2zdn2LkH11XxjKk3DaFD+mLgi8T3HYFlpfYtT2xLIiKXi8hUEZm6bt26ei6iqS4RoX23tjtEYABo0b45Hm/yM5TL5aJZBcNyjdkR1FtwEJFvRGROiq8TSr3ndpylwV6v6fFV9XlVHaSqg1q33rbFWIypyDGXHo7bU/a/iZPyw8fAI/ZKU6mMqX/11qykqsMq2y8iFwIjgMN1a8fHCqB0DuNOiW3GpEX77m35+zs3cv/5TxCPxonH4zRv24y7PrklZY3CmB1FujqkjwIeAYaq6rpS2/sCb7C1Q3oM0Ms6pE26RSNRFs1Ygj/DR9d+nW1pT7NDaHAd0sCTgB/4OvGfbJKqXqmqc0XkHWAeTnPTNVUFBmO2B4/XQ599bMU2s/NI12ilnpXsuwe4ZzsWx5hGLx6PM+2rWSyZvZQOPdux74iB1uxlton99RjTyBVsKeDGoXewavEawsEIvgwvOc2zefyHexr1HBOTXg1hKKsxZhu8dNsbLP1lBUX5QWLRGEV5Qdav2Mgjlz2b7qKZRsyCgzGN3Li3JpbJhgtOSvDpX/9MJBxJU6lMY2fBwZhGTuOpRxyqKjvAKsAmTSw4GNPIHXjyEDzessuxulzCngfvhs9vSf1M7VhwMKaRu/zB82jTuRUZ2QEAAtkBmrRqwo0vXpXmkpnGzEYrGdPINWmZw4tzH+WHD6fw26zf6dirPUNP359Apr/qDxtTgbTMkK5rNkPaGGNqrrIZ0tasZIwxJokFB2OMMUksOBhjjEliwcEYY0wSCw7GGGOSWHAwxhiTxIKDMcaYJBYcjDHGJLHgYIwxJokFB2OMMUksOBhjjEliwcEYY0wSCw7GGGOSWHAwxhiTxNZzMKaRWzRzCWNe/45oOMrBp+5HvwP7ICLpLpZp5Cw4GNOIvf3gh/z3n6OIhCKowuiXx3L4uQdzwzOXp7toppGzZiVjGqm1S9fx2p3vECoKE48rqkqwIMSY/01g3o8L0l0808hZcDCmkfrp8xkpm49ChWG+/2ByGkpkdiQWHIxppLx+L+JKDg4utwt/hi8NJTI7EgsOxjRS+58wCI0nrwHv9ro59KwD01AisyOx4GBMI9WkRQ63vn49/gwfGdkBAll+fAEvVzx0Hp37dEx38UwjZ6OVjGnEDjhxH95a8TyTPp1GNBJjn6P706Jd83QXy+wALDgY08hlN8ti2LkHp7sYZgdjzUrGGGOSWHAwxhiTxIKDMcaYJBYcjDHGJLHgYIwxJomoJk+iaWxEZB3wR7rLUUutgPXpLkQd2VGuZUe5DrBraYga0nV0UdXWqXbsEMGhMRORqao6KN3lqAs7yrXsKNcBdi0NUWO5DmtWMsYYk8SCgzHGmCQWHNLv+XQXoA7tKNeyo1wH2LU0RI3iOqzPwRhjTBKrORhjjEliwcEYY0wSCw5pIiKnichcEYmLyKBy+24VkUUiskBEjkxXGatLRI5KlHWRiNyS7vLUhIi8LCJrRWROqW0tRORrEVmY+LfB58AWkV1EZJyIzEv8XV2f2N4YryUgIpNFZFbiWv6Z2N5NRH5K/J29LSKNYrk7EXGLyAwR+TTxulFchwWH9JkDnAxMKL1RRHYHzgT6AkcBT4uIe/sXr3oSZXsKOBrYHTgrcQ2NxSs4P+fSbgHGqGovYEzidUMXBW5S1d2BfYFrEr+HxngtIeAwVd0L6A8cJSL7Ag8Aj6pqT2ATcEn6ilgj1wPzS71uFNdhwSFNVHW+qi5IsesE4C1VDanqEmARsM/2LV2N7AMsUtXFqhoG3sK5hkZBVScAG8ttPgF4NfH9q8CJ27NMtaGqq1R1euL7PJybUUca57WoquYnXnoTXwocBryb2N4orkVEOgHHAi8mXguN5DosODQ8HYFlpV4vT2xrqBpbeaujraquSny/GmibzsLUlIh0BfYGfqKRXkuiKWYmsBb4GvgN2Kyq0cRbGsvf2WPAX4F44nVLGsl1WHCoRyLyjYjMSfHVaJ6sd3bqjPVuNOO9RSQbeA+4QVVzS+9rTNeiqjFV7Q90wqmd9klviWpOREYAa1V1WrrLUhu2TGg9UtVhtfjYCmCXUq87JbY1VI2tvNWxRkTaq+oqEWmP8/Ta4ImIFycwvK6q7yc2N8prKaaqm0VkHLAf0ExEPImn7sbwd3YAcLyIHAMEgCbA4zSS67CaQ8PzMXCmiPhFpBvQC5ic5jJVZgrQKzECw4fTmf5xmsu0rT4GLkh8fwHwURrLUi2JtuyXgPmq+kipXY3xWlqLSLPE9xnAcJw+lHHAqYm3NfhrUdVbVbWTqnbF+X8xVlXPobFch6raVxq+gJNw2htDwBrgy1L7bsdpY10AHJ3uslbjWo4Bfk2U+fZ0l6eGZX8TWAVEEr+PS3DahccAC4FvgBbpLmc1ruNAnCajn4GZia9jGum17AnMSFzLHOAfie3dcR6UFgGjAH+6y1qDazoE+LQxXYelzzDGGJPEmpWMMcYkseBgjDEmiQUHY4wxSSw4GGOMSWLBwRhjTBILDsZsAxGJicjMxMz3USKSmdjeTkTeEpHfRGSaiHwuIrsm9o0Wkc3FWTqNaYgsOBizbYpUtb+q9gPCwJWJCWkfAN+qag9VHQjcyta8Rg8B56WnuMZUjwUHY+rOd0BP4FAgoqrPFu9Q1Vmq+l3i+zFAXnqKaEz1WHAwpg6IiAdnTYvZQD+gUSZbM6aYBQdjtk1GIrX0VGApTn4jYxo9y8pqzLYpUie1dAkRmcvWxGrGNEpWczCm7o0F/CJyefEGEdlTRA5KY5mMqRELDsbUMXWyWZ4EDEsMZZ0L3IezEhsi8h1ONs7DRWS5iByZvtIak5plZTXGGJPEag7GGGOSWHAwxhiTxIKDMcaYJBYcjDHGJLHgYIwxJokFB2OMMUksOBhjjEny/y3kIlX9Dub2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Using K-Means as an example\n",
    "kmeans = KMeans(n_clusters=3, random_state=random_state) # Choose the number of clusters you want\n",
    "\n",
    "kmeans.fit(pca_data)\n",
    "labels = kmeans.labels_ \n",
    "\n",
    "# Add labels to the dataframe\n",
    "X = X_radiomic\n",
    "X['cluster'] = labels\n",
    "\n",
    "# Visualize the clusters\n",
    "plt.scatter(pca_data[:, 0], pca_data[:, 1], c=labels)\n",
    "plt.title('KMeans Clustering with PCA')\n",
    "plt.xlabel('PC1')\n",
    "plt.ylabel('PC2')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot clustering results as a function of molecular subtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='Molecular Subtype'>"
      ]
     },
     "execution_count": 344,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAFKCAYAAAAE+AT0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAx5UlEQVR4nO3deZxcVZn/8c+XJJCwxLCEiAZMUAYIW4DIgKCyCIowEvZtNDAgMiIQcQZRZkaU3XFYFEZ+YTFhVAKCQFgMWwAVWZKQkAQiAwSQRJbAsAuEJM/vj3MrqVSql3TfW1238n2/Xv1K1a2q+5zudD916pxzn6OIwMzMymeVnm6AmZl1jRO4mVlJOYGbmZWUE7iZWUk5gZuZlVTvRgZbb731YsiQIY0MaWZWelOnTn01IgbWHm9oAh8yZAhTpkxpZEgzs9KT9Hy94x5CMTMrKSdwM7OScgI3Myupho6Bm1lr+/DDD5k7dy7vv/9+TzellPr27cvgwYPp06dPp57vBG5muZk7dy5rrbUWQ4YMQVJPN6dUIoLXXnuNuXPnMnTo0E69xkMoZpab999/n3XXXdfJuwskse66667QpxcncDPLlZN3163oz84J3MyspDwGbmaFGXLabbme77nz9unwOS+99BKjR49m8uTJDBgwgEGDBnHRRRdxwAEHMGvWrBWOOXbsWPbaay8+9rGPdaXJAFx22WVceuml9OrVizXXXJMxY8YwbNiwLp+vonkT+BkfqXPszca3w8xKIyLYf//9GTVqFOPHjwfgscce4+WXX+7yOceOHcuWW265Qgl84cKF9O69NL0eccQRHH/88QBMmDCBU045hYkTJ3a5TRUeQjGzlnHvvffSp0+fJckSYJtttmHDDTdccn/s2LF861vfWnJ/33335b777mPRokUcddRRbLnllmy11VZceOGFXH/99UyZMoUjjzyS4cOH89577zF16lQ+//nPs/322/PFL36RF198EYBdd92V0aNHM2LECC6++OJl2tW/f/8lt999993c5gmatwduZraCZs2axfbbb9+l106fPp158+YtGWZ54403GDBgAJdccgk/+clPGDFiBB9++CEnnngiN998MwMHDuTaa6/l9NNP56qrrgJgwYIFbdZ7uvTSS7ngggtYsGABkyZN6to3WMM9cDMzYOONN2bOnDmceOKJTJw4cZlec8WTTz7JrFmz2HPPPRk+fDhnnXUWc+fOXfL4oYce2ub5TzjhBJ555hnOP/98zjrrrFza7ARuZi1jiy22YOrUqe0+p3fv3ixevHjJ/cq667XXXpvHHnuMXXfdlcsuu4xjjz12uddGBFtssQXTp09n+vTpzJw5kzvvvHPJ42ussQYAp59+OsOHD2f48OHLneOwww7jpptu6sJ3tzwncDNrGbvvvjsffPABY8aMWXJsxowZvPDCC0vuDxkyhOnTp7N48WJeeOEFHnnkEQBeffVVFi9ezIEHHshZZ53Fo48+CsBaa63F22+/DcCmm27K/PnzefDBB4FUOuDxxx9frh1nn332kiQP8NRTTy157LbbbmOTTTbJ5fv1GLiZFaYzy/7yJIkbb7yR0aNHc/7559O3b1+GDBnCRRddtOQ5O++8M0OHDmXYsGFsvvnmbLfddgDMmzePo48+eknv/NxzzwXgqKOO4vjjj6dfv348+OCDXH/99Zx00km8+eabLFy4kNGjR7PFFlu0265LLrmEu+++mz59+rD22mszbty4fL7fiMjlRJ0xYsSI6PSGDl5GaFY6s2fPZvPNN+/pZpRavZ+hpKkRMaL2uR5CMTMrqU4lcEkDJF0v6c+SZkvaSdI6ku6S9FT279pFN9bMzJbqbA/8YmBiRGwGbAPMBk4D7omITYB7svtmZtYgHSZwSR8BPgdcCRARCyLiDWA/oDISPw4YWUwTzcysns70wIcC84FfSJom6QpJawCDIuLF7DkvAYOKaqSZmS2vMwm8N7Ad8POI2BZ4l5rhkkhLWeouZ5F0nKQpkqbMnz+/u+01M7NMZ9aBzwXmRsTD2f3rSQn8ZUkbRMSLkjYAXqn34ogYA4yBtIwwhzabWVnUWw7crfN1vJS4GcvJ/v73v2f06NHMmDGD8ePHc9BBB3X5XNU67IFHxEvAC5I2zQ7tATwBTABGZcdGATfn0iIzsy6qlJPdddddeeaZZ5g6dSrnnntut8vJ/vWvf12h1yxcuHCZ+xtttBFjx47liCOO6HI76unsKpQTgV9JmgEMB84BzgP2lPQU8IXsvplZj2nWcrJDhgxh6623ZpVV8r30plOX0kfEdGC5q4BIvXEzs6bQzOVki+BaKGZmLFtOdp999mGvvfZa7jnV5WQBFi1axAYbbLDk8fbKyRbBCdzMWsYWW2zB9ddf3+5zOione8cdd3DZZZdx3XXXLelZV1TKyVaqEdaqLid7221pP9BKRcIiuBaKmbWMZi0nWxT3wM2sOA2uINqs5WQnT57M/vvvz+uvv84tt9zCD37wg7qJf4W/X5eTNbO8uJxs97mcrJnZSsAJ3MyspJzAzcxKygnczKyknMDNzErKCdzMrKS8DtzMCrPVuK1yPd/MUTM7fE4zlpO94IILuOKKK+jduzcDBw7kqquu4hOf+ESXz1fhHriZtYxmLSe77bbbMmXKFGbMmMFBBx3Eqaee2uX2VHMCN7OW0azlZHfbbTdWX311AHbccUfmzp2by/frIRQzaxllKCd75ZVXsvfee3epjbWcwM3MaEw52V/+8pdMmTKF+++/P5c2O4GbWcto5nKyd999N2effTb3338/q622Wpe+v1oeAzezltGs5WSnTZvGN77xDSZMmMD666+f2/frHriZFaYzy/7y1KzlZP/1X/+Vd955h4MPPhhImxxPmDCh+9+vy8maWV5cTrb7XE7WzGwl4ARuZlZSTuBmlqtGDsu2mhX92TmBm1lu+vbty2uvveYk3gURwWuvvUbfvn07/ZpOrUKR9BzwNrAIWBgRIyStA1wLDAGeAw6JiNdXsM1LDDnttmXuP9f578HMmsTgwYOZO3cu8+fP7+mmlFLfvn0ZPHhwp5+/IssId4uIV6vunwbcExHnSTotu//dFTifmbWYPn36MHTo0J5uxkqjO0Mo+wHjstvjgJHdbo2ZmXVaZxN4AHdKmirpuOzYoIh4Mbv9EjCo3gslHSdpiqQp/lhlZpafzg6h7BIR8yStD9wl6c/VD0ZESKo7axERY4AxkC7k6VZrzcxsiU71wCNiXvbvK8CNwA7Ay5I2AMj+faWoRpqZ2fI6TOCS1pC0VuU2sBcwC5gAjMqeNgq4uahGmpnZ8jozhDIIuFFS5fm/joiJkiYD10k6BngeOKS4ZpqZWa0OE3hEzAG2qXP8NWCPIhplZmYd85WYZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJdXpBC6pl6Rpkm7N7g+V9LCkpyVdK2nV4pppZma1VqQHfjIwu+r++cCFEfEp4HXgmDwbZmZm7etUApc0GNgHuCK7L2B34PrsKeOAkQW0z8zM2tDZHvhFwKnA4uz+usAbEbEwuz8X+Hi9F0o6TtIUSVPmz5/fnbaamVmVDhO4pH2BVyJialcCRMSYiBgRESMGDhzYlVOYmVkdvTvxnJ2Br0j6MtAX6A9cDAyQ1DvrhQ8G5hXXTDMzq9VhDzwivhcRgyNiCHAYMCkijgTuBQ7KnjYKuLmwVpqZ2XK6sw78u8Apkp4mjYlfmU+TzMysMzozhLJERNwH3JfdngPskH+TzMysM3wlpplZSTmBm5mVlBO4mVlJOYGbmZWUE7iZWUk5gZuZlZQTuJlZSTmBm5mVlBO4mVlJOYGbmZWUE7iZWUk5gZuZlZQTuJlZSTmBm5mVlBO4mVlJOYGbmZWUE7iZWUk5gZuZlZQTuJlZSTmBm5mVlBO4mVlJOYGbmZWUE7iZWUk5gZuZlVSHCVxSX0mPSHpM0uOSfpgdHyrpYUlPS7pW0qrFN9fMzCo60wP/ANg9IrYBhgNfkrQjcD5wYUR8CngdOKawVpqZ2XI6TOCRvJPd7ZN9BbA7cH12fBwwsogGmplZfZ0aA5fUS9J04BXgLuAZ4I2IWJg9ZS7w8TZee5ykKZKmzJ8/P4cmm5kZdDKBR8SiiBgODAZ2ADbrbICIGBMRIyJixMCBA7vWSjMzW84KrUKJiDeAe4GdgAGSemcPDQbm5ds0MzNrT2dWoQyUNCC73Q/YE5hNSuQHZU8bBdxcUBvNzKyO3h0/hQ2AcZJ6kRL+dRFxq6QngPGSzgKmAVcW2E4zM6vRYQKPiBnAtnWOzyGNh5uZWQ/wlZhmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl5QRuZlZSTuBmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl1WECl7ShpHslPSHpcUknZ8fXkXSXpKeyf9cuvrlmZlbRmR74QuA7ETEM2BE4QdIw4DTgnojYBLgnu29mZg3SYQKPiBcj4tHs9tvAbODjwH7AuOxp44CRBbXRzMzqWKExcElDgG2Bh4FBEfFi9tBLwKA2XnOcpCmSpsyfP787bTUzsyqdTuCS1gRuAEZHxFvVj0VEAFHvdRExJiJGRMSIgQMHdquxZma2VKcSuKQ+pOT9q4j4bXb4ZUkbZI9vALxSTBPNzKyezqxCEXAlMDsiLqh6aAIwKrs9Crg5/+aZmVlbenfiOTsDXwVmSpqeHfs+cB5wnaRjgOeBQwppoZmZ1dVhAo+IPwJq4+E98m2OmZl1lq/ENDMrKSdwM7OScgI3MyspJ3Azs5JyAjczKykncDOzknICNzMrKSdwM7OScgI3MyspJ3Azs5JyAjczKykncDOzknICNzMrKSdwM7OS6kw9cCuhIafdtsz9587bp4daYmZFcQ/czKyknMDNzErKCdzMrKScwM3MSsoJ3MyspJzAzcxKygnczKyknMDNzEqqwwt5JF0F7Au8EhFbZsfWAa4FhgDPAYdExOvFNTPZatxWyx2bOWpm0WHNzJpSZ3rgY4Ev1Rw7DbgnIjYB7snum5lZA3WYwCPi98D/1RzeDxiX3R4HjMy3WWZm1pGu1kIZFBEvZrdfAga19URJxwHHAWy00UZdDJc/1woxs7Lr9iRmRAQQ7Tw+JiJGRMSIgQMHdjecmZlluprAX5a0AUD27yv5NcnMzDqjqwl8AjAquz0KuDmf5piZWWd1mMAlXQM8CGwqaa6kY4DzgD0lPQV8IbtvZmYN1OEkZkQc3sZDe+TcFjMzWwG+EtPMrKScwM3MSsoJ3MyspJzAzcxKygnczKyknMDNzErKCdzMrKScwM3MSsoJ3MyspJzAzcxKygnczKykurqhg5lZl9RupgLeUKWr3AM3MyspJ3Azs5LyEMrK4oyP1Dn2ZuPbUUL+yG/Nyj1wM7OScgI3MyspD6FUFDzEUPsx3B/Bu67ukEbfI5Z/ooeIuq7278E/y6bkHriZWUk5gZuZlZSHUFZiW43bapn7M0fN7NJ5VspVGi20qqf+kFQPNKQgrTx86R64mVlJOYGbmZWUE7iZWUl1awxc0peAi4FewBURcV4urbLya6Ex4kZbKecUatTOz0DX52h6XIF/C13ugUvqBVwK7A0MAw6XNCyXVpmZWYe6M4SyA/B0RMyJiAXAeGC/fJplZmYdUUR07YXSQcCXIuLY7P5Xgb+PiG/VPO844Ljs7qbAk10Itx7wapca2jWtHK+VvzfHc7xWjfeJiBhYe7DwdeARMQYY051zSJoSESNyatJKHa+VvzfHc7yVLV53hlDmARtW3R+cHTMzswboTgKfDGwiaaikVYHDgAn5NMvMzDrS5SGUiFgo6VvAHaRlhFdFxOO5tWxZ3RqCcbwei+V4jud4Bcbr8iSmmZn1LF+JaWZWUk7gZmYl5QRuZlZSK30Cl7RhO4/t28i2tLJspVLRMVYvOkYjSapTRGPJY7mvXZa0TTuP/XPe8az7mi6BSzqg6vbaDQh5l6QhddrxT6RCXbmStF17XwXE273q9tCaxw5Y/hXdivUfbRz/CHBnnrFqzv8ZSU8Af87ubyPpv4uKVyf+30m6vIBT313vb0DSXsCNBcS7UdL2deL9EPh6AfEq5x8naUDV/bUlXVVUvDbasGcB5+wtSdntDSUdJGnbXINERFN9AY/Wu11gvC8D/wtsUnXse8BMYHAB8e6t+noLmFR1f1Ijf555/3xJSfrsmmODgOnAfxT4f/gw6aKyaVXHZhUQZ+vse5wFnAVsANwAzAW+XUC8r2c/u4FVx44AngW2LiDe9sAcYKfsvoDLst/N/gX+/03rzLEiv4C/FPB/93/AX7Lb/0uqF/Uk8N284jTjlmpq43YhIuJ2SR8Av5M0EjiWVKjrcxHxegHxdqvcljQtInZv7/k5aO/nmffP9yvA9ZIuiIhTJG0C/A74SURclnOsZUTEC1lnp2JRAWEuB34OPAh8iZRcxwFHRsT7eQeLiMslvQ9MynrdhwLHA7tFxHMFxJua/Q3cKOkElva6vxQRH+Qdr8oqktau/L1JWocCynxIautCQwHr5hxuNPBJYC1gNqmWyavZMN9k4Pw8gjRjAu+XfcxYBeib3V7ylxkRj+YdMCLukXQ0cB/wJ2D3Iv4g64VucIzaeLnGj4j3Je0PXCvpGuAzwOiIKOLjfrUXJH0GCEl9gJNJfzR5Wy0ixma3n5R0ckScWkCcJSLif7IkPo3Um9slIgopvpQlzrnAKOAm4G7gW8AaktaIiP8rIi7wX8CDkn5D+ls/CDi7gDifBf4ReKfmuEidtjwtyN6QXpf0dOX/LCL+JmlBXkGaMYG/CFyQ3X6p6jakhJNrj1XS29l5BawG7AG8ko1dRUT0zzNeD9g463mo6jbZ/aFtv2zFSTolu/kwcCrwB2Bo5XhEXNDWa7vpeNJ8xcdJ9XjuBE4oIE5th+KD6vt5dy4kzWTp7+bqpF7ipKrfza3zjAdMZemb+tvA3wOPZPED2DjneABExNWSprD0b/uAiHiigFAPAX+LiPtrH5DUlSqp7anuiK5a9XsiILcto0t1JaakPhHxYU+3ozsk/YylfySHkcbFloiIk3KO9/n2Hq/3y9yNWD/oINYP84rVEyTdR9ufWiLv4TBJn2jv8Yh4Pu94eZ+zg3j9I+KtrOe/nAJ7/IWTdG97j1cPpXYrTrMn8Ky3sTtp8mbfiBiU8/n7knpwnwJmkGq6LMwzRk28Ue09HhHjCoy9Thaj0D8MSesV9TG/jXgDSeO1Q6j6VBkR/9SoNhQlG5P+FDAzIu4oONajEZH7Sqh24t0aEftKepZl3xgrnzAK6fG30ZYHImLnHM/3sYj4a17nazNOsyZwSTuSkvZIYB3SR+IJeU8sSroW+JD0cX9v4PmIODnPGD1J0kbAj0lvgm+S/jj6k1a/nJbnZFi2bv4XpJ/nYuCQiPhTXudvJ+6fSP9/U6mavIyIG3KOc05EfD+7vWdE3JXn+evE+29gC9K8zB7ALRFxZoHxpkVEvsvcSkLSCxHR5jUhXTjf7aS8dR8wEfhjER3Dpkvgks4BDiZN2FxDWu86JSJyHa+tijczIrbKbvcGHim6F5L1wk8m7VAEacLtpxFxdQGxHgQuAq6PiEXZsV6kn/HoiNgxx1gzSEn7z5L+HvhxRLQ7hJNT3OkRMbwBcZb0UBvRW5U0C9gmIhZlqxf+EBHLrdPOMd4r1AzpVct7eK8q7j0RsUdHx4ok6S8RsVHO5+wL7ErqGO5MymkTgYkR8Zc8YjTjJOaxpDWTPyf1OD6QVOS7zJIx9UglcgsMtSR5jwZOAR4l9Yi3A/5TUkTE/+Qccr2IuLb6QJbIx0vKuze3MCL+nMV4WNJaOZ+/LbdK+nJE3N6geI2yoPKmm61eKHpZ7XukTzENkSW41YH1sguWKt9ff9KEdN7x2rpwTUC/vONlK9kmZl+VC+n2Bi6R9NGI6PbKl2bsgfcC9gQOJ31svBf4ArBhIR9BpEXAu5W7pP/Iv7F0HC7XVSiSHgIOqx26ULoadHyePeLsvONJFxSMA17IDm9IWiq2XkQckmOsuSy7auiU6vtFrULJVhKtAXxAekMu6v+u8v0J+DbLfq+5f3+S/gY8XblLWlf8NEu/v1xXofTAGPjJpM7Mx0irhyoJ/C3g8oi4JOd4v2jv8Yg4Os94bbThs6TcNjrSZvDdO1+zJfBqklYD9iV9w58F7omII3q2Vd0j6YmIGLaij3Uj3qrAMcB+LO3VzCPtnnRlnhdodLAKJSLiR3nF6gmNXmXTA6tQHqrXgZC0CnB4RPwqz3hV5z8xIn5WxLmbQbaE8AjSsOWzwA15vTk1dQKvln0cH13kJE5VrDWA/Um/tPvkfO6pbY1jtvdY2Un6dERMzvmcm2Xj7XV7jXmvy24WknYh/W7mutZdqWbNN0lv9BOAu0gX8nwHeCwi9sszXk3sLYFhVK2RLmhO6PPA6xExQ9IhwOeAZ4D/zrkz83ekjufhpF3orwX+JSLafVNe4TjNlsCzIZRDSL9EEyNiVra64ftAv6JmybOe6j6kd8ovkmpc/DYibsk5TvXH4mUeAjaOiDVyjteb1AMfybI98JtJPfDC1tVLGsbSX+I3IufdvyVdHhFfb2PNbRHrsrcAPhkRE7L7FwKVioGXFPmGUacX99u8e62SbgZeJ5UK2ANYn/R7eXJETM8zVk3cH5Am+4YBt5PGif8YEQflHOdSUj2bvqSaJGuSxqd3BlaJiCNzjLWYtDLqmIh4Ojs2J++lkc2YwMeSxmgfIV0J9ldgBGnJ200FxNuLlGD2Io23Xwv8LCKG5B0ri9foj8XXAG+QxsDnZocHk8bA14mIQ3OON4SlSftD4BPAiDyXK/YUSbcA51aWRipVQPx30kTcgRExMud4DenFVcWrXpHVi3RV9EZRcFkJpStOtyEVsNpG0iDglxGRa4XAyhBlNnk6D1g/W+EjYEble88p1kjShXo7k94kxgNX5L2arhlXoYwgVVpbnP2gXyL1el4rKN5E0jvlLhHxLICk3MvIVuSdoDth+4j4u5pjc4GHJP1vnoGyJYv9Sb+sB0bEU5KeLSp5t7OqAICI+G3OITeIZde1vxXZWnNJ38g5FqTyuH8gXcBW6cV9u4A4FdUrshZJmlt08s68l/29L5TUH3iF1InL2/uwpGbP81UrfEJSrp9Es87mTdlw7H6kydr1Jf0cuDEicimv3IwJfEFELIYlP+g5BSZvSEv4DiPVXp5DSj69igqmpbVXlnuIYmqv/J+kg0kTJ4uzNqxC+iied7XFl0nDNIOAgcBTFFuw6x/aeSyAvBP4Mssiayb81s85FsABpN/NeyVVenFFLiXcRtJb2W2R6nm8RXG/mxVTlOqBX05axvgOaRgnb+sr1eVR1W2y+wMLiEdEvAv8Gvh1tlTyYOC75FQfvxmHUBq6dKom9mdIH1cPBB4jvVOOKSpeI2RDGueTrsR8nfRzHMDSKzGfzTneR0iJ53BgkyzWFyPikTzj9IRsrP20iHi45viOwHkRsWtBcSu9uMNJ/49Xk2Mvrplkv6/9I2JGAeduuVo9zZjAGzpG3EYbViGtPT8sWqCeRoWkdQEK/kRTHW8QaUL6MNJYahEfi8l6b19j+VooeRcG24E0Dj2WdBEWpE0QRgGHNuJNqqoXd2g08ErFoknamuX//3L9BCVpx4h4KM9zthNrK9Inio+TauJ/N5bWO38kcriIB5owgdcjaT3gtSigsW0tQaso+1I0pVoor2TDUQKOIg0bPUG6WKKwwl017fhEUW++SrVQHiLtorS4cjwKKAwmaX3S0rotskOPA5dGxMsFxKpbpa8iSlytr5rS9mlbk36Wlf+/yLvzJOlR0mYKp0bEm3meu06sP5J2bXqIdHX50cBXIuIZ5VhzpukSeOXjKOnqwTOB/wHWI9XV/VpETMw53mLSFlmV6nnVY4y5L0VrNKV6GjtEuhT7fNKQ1E1ktZfz/CNR2zuekMX6Sl6xauI25ApCSRtFTjUsOhlvMWnCufImW/u72bBqfUUq4gK2NuKsApxEWut+ZuRftqI61mMRsU3V/d2AMcBXSWvOc/l9bcYEPoW05vsjpG9474h4SNJmwDV5rwOXNJq0A8ibpEmiGyOidseO0qr+45A0Ffh01WTmMr9kOcSaT7pc/xrSpg7LTLhFjrXHa+J+mzTxdSvpcvpKvFx7qFq2mNUNEXFgnuevE+8iYDfgAdLP9I9FfArtaZKuBP4ritnEoV68YaRJ0lVYumFGrpO0kh4jbcv4ZtWxrUnXl6wTEbls4dZ0u9IDvSPizoj4DfBSZcwqsiJJeYuIiyJiF+BE0tKleyRdJ2l4EfF6wAtaujP9c2TLsyrj4Tn7KOnNd0vSDjl7Aq9GxP1FJe/MAuA/SX+UU7OvKQXEqX5DKrz3GxGjgeHAb0g9t2mSfqxUFKmVXE3aUu1JSTMkzVSqbJk7SceQLmI7nTRZ2j8i1ipghc35wObVB7KJ2T3IcXVUMy4jXFx1+72axwrrfUTEnOxKtH6kP5a/I21aW3bHAldLOoP0KWO6pOmk1SGntP2yFZetq50ITFSqY3M4cJ+kH0bOhYlqfAf4VBS/iUR7+4sWEzD1uO+VNI00GXwmaXnm5Y2I3yBXkv7mlpnDyFs2V/Ic8NmIeKmoOAAR8es2jld2qc9FMw6hVKoD1lYGBOgbEX1yjrcx6Q9jP9LH//HAbRFR++ZRapI2J70p9SaNq06uDKXkHGc1UkmCw0mrCiaQdjmal3esqph3AiMj4m9Fxcji1PvdhILWSVctHzyUtE75t8B1jRyHbwRJD0bETg2I84WIuLvoOFmsX9D+9nvH5BKn2RJ4o2UTRTNIH6veouaHHsVtxNsQSvUffh0RDzQg1tWk4ZPbSaVxZxUdM4t7I2lVyL0sOwZeyAYEjSLpXVJvezx1LorKe5ldT1HaeWgAcAvL/v/lvYzwp+09nufvi6R68yMbksoQ94qIwbnEKUMCV7Z3XkHnPoN2Pg6XcXF/NaWay4cBGwDXkSaCpxUUazFLa6vX2+OwkCv51MY+o0UsI2wkpbpA7fXiWuIaBdWv013EMsIFpBVn15FqLNVOshfy+5J9yv8+qfLhhaQict2uBQ7lSeC5rZvsRhu+FxHn9mQbuiO7QOqw7KsfaVXDNRGRaz2URlMqunR35LTLdzORtBPwUCuuPOkJ2cT9waQhqYWki7Kuj4g3Coq3GfBvwLakSfZf5n3dRTOuQqmnkB7jCjq4pxvQHRHxfEScn70RHk4qLzu7Z1vVfdnE6eLsEv5W8zVgqqTxko6S9NGeblARJA2WdKOkV7KvGyTlMsRQLSJei4jLsjf7o0nDNk9I+mresST9hjSU+CCpVO4EoL+kdTq6QGuF4vjNvXOa4VNAdyjVBd+b1APfg7Rb9jURcXNPtisP2eqhbUkbEFSGcEo/Bl6R9eT2JtWp/whprH8i8ED2BlZqku4iFXyqXFjzj8CRkXM52ap425E6MXuSlpzmvgZd0nMsHf6q/FsZssntIqymS+CS1gROJRWUGkxa4/sMcFlEjO3BdjV0v8C8SKrsL/plUo318cDNkaqktYRWHQOvR1I/0sU9ewM7Rc6bZPQESdMjYnhHx3KI8yPSCqnZpL+DiXkPaXSyHevkdZFZMybwm4EbgbtJhZDWIP2w/w2YFxHf76F2lbIHLmkSqXdzQ2TFdNp43trtPd7sssS2UUQ82dNtyYukAyorMcr+/9MeSfcAvyDNy0DqcBwdORfryibZn2Xp8s/qnnFEjpVOJV0REcfWOT6Y9MaxZS5xmjCB19YQmBwRn87qGDwREZsVELMXcFJEXNjOc74fEefkHbtZlPUTBoCkfwB+AqwaEUOzq2h/FAXVXmmUmkv3S/v/05Fsgv1nQGUt+AOkv8dc17urgZVOJY0j7SvwtVhaumIYqdzDj/IaTWjGScx3lTZtRdJXSEWtyH4IhRSzz8YRD+/gOS2bvDNFbhRQtDOAHUhbxxFp/8ZWKPSkNm63lGyC/SsRMTD7GlnQxUpjSLXq+2Uxl/nKOdZRpJ7+tZJ6Ke01cAfpjWlsXkGa8VL644ErJG1CKi/5TwCSBgKXFhj3AUmXkJYWVU+Elbqc7Aporo9iK+bDiHhTWibHFXZJdgP1U9rMeBWgb3Z7yTfZKr+bkn5MKr36Hmlydmvg2xHxy5xDjQK+BJyhtN/ow1m8u/OeE8qWfh6XXTx0H2lv2IMj53rkTTeE0lPUoJ3Nm1WZP6Jn1ezuAU4jTX6fBPSJiON7tGHd1MbvZEXL/G5WJiwl7Q/sS6rR8/vIsVJmnZirkDZN35u0Kus94M6I+HFO5/8ZSysdHkHaAGTJst28Vkg1Yw98OZImFf3L2ooXgqygMn9EP5FUXe4D0kTYHaSiT6XW3u+kpFxrAvWwSh7aB/hNnU9TucuGZB/Mvv5DadOYL+YYYkobt3PVdD3wOmUkRSrC9CRAnjPFNXEHAecAH4uIvbMJh50i4soi4jVKZ1cy5Lm0yYqhlNV2J/Xo9o2IQT3cpFxIOo90Ydl7pLmMAcCtEfH3OcfpS7oK83VS3ZVTgc+SlimfGcVXs8xdMybwCaSiUpUxMQF/AHaB4vbElPQ70lKm0yNim+zCl2kRsVUR8RplZVjJIOkWlh/Df5PU8/l/EfF+41uVH6Vdqo4gJbl1gBOACa20rDC7OvHNiFikVIVxrci55Kuk64APSUuT1ybVRbmFlFuGR871lrLrE04GNs0OzQZ+GhFX5xWj6VahZEu/biDNGG8TEc+RJqmKmCmutl5EXEc2+ZUt8C/9VW6sHCsZ5pB25Lk8+3oLeJv0ya20dbMlnSPpKeBsUsXMbYH5ETGulZI3pN2TsuQ9JiLezTt5Z4ZFxJGkHbg2jYgTImJiRPwb2UYnecmS92hSrfqPkTY3PhU4Oc9L95tyDDwiblSq8Xym0g4aqzYg7LtZsZuAJb2eQjc+bZCVYSXDZyLi01X3b6m6fuDxHmtV9x0L/C/wc+CWiPhAUnN9ZM5fkVeWLoDUOZP015rH8u6s/TOwf9YBrZikVGZ2PEvLBnRLUyZwgGxZzymStmHpAv8inUIqOPNJSQ+QCugf1IC4RXsRqNQ0f6nqNqQ3q1ZYybCmqjYclrQRsGb2WC5lO3vIBqR6HYcDF2WrUvpJ6t0Tl4A3yCsFnntwtqxPVbfJ7n8851j9a5I3ABHxnKTcyio3bQKv8gzQS9KAKKjsI6SeqKTPk8arBDwZER8WFa9RVpKVDN8B/ijpGdL/3VDgm9lYamnrocTyW9TtSyoFPE/SPRFxRI82sAAR8aUCT/+vVbdrV4bkvVKkvR29ctvtqxknMf87Ir6Z3d6FVMfjGeBTwDci4vYCY3+GtA3Ykje2PCccmkELr2RYDaiUWXiy7BOX7ZG0FjA6Ikq/VLIt2Vj4cTmf89CIuDbPc7YT62/A0/UeAjaOiDVyidOECbx61cS9wHey3vHGpP0ACxkjk/Q/wCdJGxlXxsMirwX3PW1lWMkAxfzh95SsRs8hpI/3EyNilqR9Sbu79IsSFlerprbrYgt4LHLadqwq3q2kztk3I2JOnueuE2sTYBBpn91qGwIvRUS95L7Cmn0IpX9lki3SrvFFrpoZQZqlbq53tG6SdA5pM4q/kC5y+SEwJVqw1Gqm9OVVq1xJ+oN/BPhpNvE2AjgtIm7qyYblZD7wPMuujqpcvbh+3sEiYl9JI4HbJP2aNDm8uOrxPK+DuBD4Xu3KuWz8+0LgH/II0owJfLPsYh4BQyoXn2TJu8jVKLOAj5Im/VrJyraSochJsEYbAWwdEYuzi1BeAj4ZEa/1cLvyMgfYo17hKkm1PddcRMRNkp4Ffg8cw7KbLuRZAG1QRMysE3+mpCF5BWnGBL55zf13sn/XAf4j72BVF4GsRdpe6RGW3Rm71CVJWclWMhQ8CdZoC7JLvomI9yXNaaHkDXAR6YKaepUHc6lJUi2bJ/k30uqyIyPi1rxjVBnQzmP98grSdGPgjZatPGlTRNzfqLYUrWolw+GkS4hbciUDtMZYeM1EmEhzNE9TwAYEKwNJT5IuEjwzInJbCdJGrGuASRFxec3xY4E9I+LQXOKUKYFL+l1E7F3Quc+PiO92dKxVtMJKhkZPgjWaGrgBQU+QdE5kO2xJ2jMi7io43rDIee/LdmINIu0stoC07yakIbFVSRf45HKladMlcKUNR+s+RCpws0FBcZerEyJpRtl7Oa28kkHSItqeBPt4RDTiCt6GyqrmvdYKk+2NrtOT1VlqUxHDpZJ2Ayrbpz0eEZPyPH8zjoFPBu6nft2OAXkHk/TPwDdJV2BWV0JcC/hT3vF6QCuvZGj4JFgjZUs/zyPtSnUm6fLr9YBVJH0tIib2ZPtKaCfSsr5rSJs5FF4bKCLuBdqr694tzZjAZ5Mu2Hmq9oGC/ih/DfwOOJf0x/K57PgfI2JaAfEarZVXMlxEAyfBesAlpE9KHwEmAXtHxEOSNiMlobIn8PUlnUK2bDC7vUREXFD/ZV32UZZO6B8B3AZcExGlrZfTdNUISfsbttWuE/MOFhFvZjULHgJ+SerhDATGSco9Xg9YZiUD0DIrGSLi0oh4rI3Hftbo9hSgd0TcGRG/IV388RBARPy5h9uVl8tJn3TXrLpd/ZWriFiUVR8cBexImhC+T9K38o7VKE03Bt4eSYMi4uWCzj2DtIHDu9n9NYAHW2AMvGVXMjR6EqzR2hsjbsSYcSvKVmLtQ+qFDyEVsLsqIub1ZLu6qukTuKQBpH0OjwA2j4iPFRRnJvDpSg2NbLhhcpR/Q4eWXcnQ6EmwRssmad8lvdn2I+1yTna/b0SUuhiZpPau64i8V0hJupo0oXg7MD4iZuV5/p7QjGPgSOoH7EdK2tuSPk6NJF09VZRfAA9LujG7P5I0AVhq9RJ0K61kaGUR0aun21CwejvBr0G6QnJd8t/X9B+zmCcDJ2npvpuVT6O5lXltlKbrgWc1Cj4L3EkqfD4JeDoihjYg9nZkW7cBf2iFScz2VjIApV7JIGkuqb65gG+zbK3zIibBrCDZdQknk5L3dcB/RUQrlUUoRDP2wIeRNh2dDcyOtM1SQ95lssJZrbBDTbVWXslQmfiqvW0lkV2MdQpwJKl2+3atViGzSE3XAwfIksvhpB2kXyVtsrBlUROYrUzS9IgYnt2eHRGbVz02rcwX8li5SfpP4ADS/reXRsQ7HbzEajRlAq8maXvSWPjBwNyI+EwPN6lUWnklQ6MnwSxfkhaTCsctZGlVQCjxmHSjNX0Cr8h2ktklIv7Q020pk1ZeySDpO3UOL5kEi4g16zxu1jKaPoG36hZgli9PgtnKqBmvxATS6gmlXaOfB24mLSHcrP1X2cpG0jqSzgJmkCblt4uI7zp528qg6RK4pHMkPQWcTfqj3BaYHxHjPDtt1bJJsMnA28BWEXGGf0dsZdJ0QyiSXiFtAXYRS7cAmxMReW53ZC3Ak2C2smvGdeAr1RZg1nUR0XSfIM0aqel64NVWpi3AzMxWVFMn8GrZKoP9I+Lqnm6LmVkzaLohlNqi7mZmVl/TJXDar2dRjo8LZmYN0IwJ/KqIqLt1WrYZr5mZ0YTrwIG7JA2pPSjpaODixjfHzKw5NWMCPwW4U9ImlQOSvpcd/3yPtcrMrMk03RBKRNwu6QPgd5JGAscCOwCf81V2ZmZLNe0yQkmfBW4E/gQcUtmr0szMkqZL4JLeJq02EbAa8CGwCF8ebWa2jKZL4GZm1jnNOIlpZmad4ARuZlZSTuCWO0kh6ZdV93tLmi/p1g5ed5SkS3JuyxmS/qUbr19d0q8kzZQ0S9IfJbW7VZuk+ySNqHN8V0ne09Vy03TLCK0lvAtsKalfRLxHKg88r4fb1Cl1yhafDLwcEVtlj29Kmljvil2Bd0grq8y6zT1wK8rtwD7Z7cOBayoPZNug3SRphqSHJG1d+2JJAyXdIGly9rVzdnxNSb/IesQzJB2YHX+n6rUHSRpb55xfz871WHbu1bPjYyVdJulh4Mc1L9uAqjefiHgy22RkiKRZVef+F0lnVL3uq5KmZ732HbKri48Hvp0d/6ykZyX1yV7fv3I/68FfXP367DlrSLpK0iOSpknar6P/BGttTuBWlPHAYZL6AlsDD1c99kNgWkRsDXwfqFci+GLgwoj4NHAgcEV2/N+BNyNiq+z1k1agTb+NiE9HxDbAbNIGyBWDgc9ERG01zKuA70p6UNJZ1VcId2D1iBgOfJNU3+c54LLsexoeEX8A7mPpm9xhWfs+rPf67NjpwKSI2AHYDfhPSWt0sj3WgjyEYoWIiBlZr/NwUm+82i6kpExETJK0rqTa9f1fAIZJqtzvn409f4GU7CpxVuTq3C2zDZAHAGsCd1Q99puIWFTn+5guaWNgryz2ZEk7Ae91EOua7PW/z3rXA+o85wrgVOAm4Gjg6x28fi/gK1Vj+n2BjUhvRrYScgK3Ik0AfkIa+113BV+7CrBj7RW4VQm9VvUFDX3beM5YYGREPCbpqKxdFe+2eeKId4DfAr/N9uH8MnAty36CrY1Ze4HFchdcRMQD2VDMrkCviJjVzvMrF7cdGBFPttVWW7l4CMWKdBXww4iYWXP8D8CRkFZmAK9GxFs1z7kTOLFyR9Lw7OZdwAlVx9fObr4saXNJqwD7t9GetYAXs3HnIzvzDUjauRJD0qrAMOB54GVg/ezTQ2Xrv2qHZq/ZhTTk8ybwNsvXu78a+DXwi068/g7gRGXvYpK27cz3YK3LCdwKExFzI+KndR46A9he0gzgPGBUneecBIzIJiqfIE0AApwFrJ1N7j1GGgsGOA24lbTC48U2mvTvpLH4B4A/d/Lb+CRwv6SZwDRgCnBDNlb9I+AR0ptK7fnelzSNNO5dGWu/Bdi/MomZHfsVsDZVk7ztvP5MoA8wQ9Lj2X1biflSerMeJOkgYL+I+GrVsfuAf4mIKT3WMCsFj4Gb9RBJPwP2Jo2pm60w98DNzErKY+BmZiXlBG5mVlJO4GZmJeUEbmZWUk7gZmYl9f8BTHlxqBosnFwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# for each subtype, calculate the number of subjects per cluster\n",
    "# create a new DF with these results\n",
    "if 'SubjectID' not in X.columns.tolist():\n",
    "    X = X.reset_index(drop=False)\n",
    "\n",
    "k_clusters = X[['SubjectID','cluster']]\n",
    "df = k_clusters.merge(molec_df)\n",
    "\n",
    "unique_clusters = df['cluster'].unique().tolist()\n",
    "unique_subtypes = df['OpenPedCan Molecular Subtype'].unique().tolist()\n",
    "results=[]\n",
    "for subtype_ in unique_subtypes:\n",
    "    subtype_info = [subtype_]\n",
    "    cluster_labels = []\n",
    "    for cluster_ in unique_clusters:\n",
    "        cluster_count = len(df[(df['OpenPedCan Molecular Subtype']==subtype_) & (df['cluster']==cluster_)])\n",
    "        subtype_info = subtype_info + [cluster_count]\n",
    "        cluster_labels.append(f'Cluster-{cluster_+1}')\n",
    "    results.append(subtype_info)\n",
    "\n",
    "cols = ['Molecular Subtype'] + cluster_labels\n",
    "plot_df = pd.DataFrame(results, columns=cols)\n",
    "\n",
    "# plot a bar graph\n",
    "plot_df.plot.bar(x='Molecular Subtype', y=cluster_labels) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run statistics to assess effect of subtype on clustering results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "F-statistic: 2.2246894808469766\n",
      "p-value: 0.014968998373574493\n",
      "There is a significant effect of molecular subtype on cluster assignments.\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import f_oneway\n",
    "\n",
    "# Group data by condition\n",
    "grouped_data = [df[df['OpenPedCan Molecular Subtype'] == cond]['cluster'] for cond in df['OpenPedCan Molecular Subtype'].unique()]\n",
    "\n",
    "# Perform ANOVA\n",
    "f_statistic, p_value = f_oneway(*grouped_data)\n",
    "\n",
    "print('F-statistic:', f_statistic)\n",
    "print('p-value:', p_value)\n",
    "\n",
    "if p_value < 0.05:\n",
    "    print('There is a significant effect of molecular subtype on cluster assignments.')\n",
    "else:\n",
    "    print('There is no significant effect of molecular subtype on cluster assignments.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save fitted PCA and KMeans models as new project attachments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 346,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['kmeans_model.joblib']"
      ]
     },
     "execution_count": 346,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import joblib\n",
    "\n",
    "# save trained models to disk\n",
    "joblib.dump(pca, 'pca_model.joblib')\n",
    "joblib.dump(kmeans, 'kmeans_model.joblib')\n",
    "\n",
    "# save as a project attachment\n",
    "proj_container.upload_file('pca_model.joblib')\n",
    "proj_container.upload_file('kmeans_model.joblib')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}