joint_rpca.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "94488c0a-820d-4cf0-884b-5b5937f0e1bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gemelli import rpca\n",
    "import biom \n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import qiime2 as q2\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from qiime2.plugins import emperor\n",
    "\n",
    "# plotting\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cf4cbbce-93b8-4eef-84b5-850bfe8249e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "352 x 1177 <class 'biom.table.Table'> with 336820 nonzero entries (81% dense)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wol2_bt = biom.load_table('/projects/thdmi/feature-tables/wol2_filtered_table_zebra.biom')\n",
    "ko_bt = biom.load_table('/projects/thdmi/feature-tables/kegg_filtered_table_zebra.biom') \n",
    "pathways_bt = biom.load_table('/projects/thdmi/feature-tables/pathways_filtered_table_zebra.biom')\n",
    "\n",
    "metadata = q2.Metadata.load('/projects/thdmi/metadata/consolidated_metadata_subset.tsv').to_dataframe()\n",
    "wol2_bt.filter(metadata.index.intersection(wol2_bt.ids()))\n",
    "ko_bt.filter(metadata.index.intersection(ko_bt.ids()))\n",
    "pathways_bt.filter(metadata.index.intersection(pathways_bt.ids()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bb9cc247-c67b-40b7-899a-ad7af8e32773",
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata_subset = metadata.loc[metadata.index.intersection(wol2_bt.ids())]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "08698d2a-a11b-4699-8f05-a2c79ff9198d",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_, test_ = train_test_split(metadata_subset, test_size = 0.5,\n",
    "                                 stratify=metadata_subset[['thdmi_cohort', 'age_cat']],\n",
    "                                 random_state=6)\n",
    "metadata_subset['train_test'] = 'test'\n",
    "metadata_subset.loc[train_.index, 'train_test'] = 'train'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9dfb0a5d-713f-4600-999a-a8cbd485b6b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/adilmore/gemelli/gemelli/preprocessing.py:418: RuntimeWarning: divide by zero encountered in log\n",
      "  mat = np.log(matrix_closure(mat))\n",
      "/home/adilmore/gemelli/gemelli/preprocessing.py:418: RuntimeWarning: divide by zero encountered in log\n",
      "  mat = np.log(matrix_closure(mat))\n",
      "/home/adilmore/gemelli/gemelli/preprocessing.py:418: RuntimeWarning: divide by zero encountered in log\n",
      "  mat = np.log(matrix_closure(mat))\n",
      "/home/adilmore/gemelli/gemelli/rpca.py:740: FutureWarning: Passing a set as an indexer is deprecated and will raise in a future version. Use a list instead.\n",
      "  sample_metadata = sample_metadata.loc[shared_all_samples, :]\n",
      "/home/adilmore/gemelli/gemelli/rpca.py:900: RuntimeWarning: Removing 0 features(s) in table(s) but not the ordination.\n",
      "  warnings.warn('Removing %i features(s) in table(s)'\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAGdCAYAAADjWSL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8J0lEQVR4nO3deXxU5d3///dkmQkJyRAIySQSwiZhSaABlMWF1UAUUEHRYnNDa7G9q1gKfBW0WrS3Ru2ttLdUay11xeKvheCCRgKyiCyyRXYIEiRIQliSyQJkPb8/EkYGwhKYMDnJ6/l4nIeZc66ZfC6OMm/PXGc+FsMwDAEAAJiMj7cLAAAAuBKEGAAAYEqEGAAAYEqEGAAAYEqEGAAAYEqEGAAAYEqEGAAAYEqEGAAAYEp+3i6gvlRVVenw4cMKDg6WxWLxdjkAAOAyGIahoqIiRUVFycfn4tdaGm2IOXz4sKKjo71dBgAAuALZ2dlq06bNRcc02hATHBwsqfoPISQkxMvVAACAy1FYWKjo6GjX+/jFNNoQc+YjpJCQEEIMAAAmczlLQVjYCwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQAwAATIkQU0cnyyrUbsZitZuxWCfLKrxdDgAATRYhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmBIhBgAAmFKdQ8yqVas0atQoRUVFyWKxaNGiRW7HLRZLrduf/vQn15hBgwadd/z+++93e538/HwlJyfLbrfLbrcrOTlZBQUFVzRJAADQ+NQ5xJSUlKhnz56aM2dOrcdzcnLctn/+85+yWCwaO3as27hJkya5jXvjjTfcjo8fP14ZGRlKS0tTWlqaMjIylJycXNdyAQBAI+VX1yckJSUpKSnpgscdDofb448++kiDBw9Whw4d3PYHBgaeN/aMXbt2KS0tTevWrVPfvn0lSW+++ab69++vPXv2KDY2tq5lAwCARqZe18QcOXJEixcv1oMPPnjesXnz5iksLEzdu3fX9OnTVVRU5Dq2du1a2e12V4CRpH79+slut2vNmjW1/q7S0lIVFha6bQAAoPGq85WYunjnnXcUHBysMWPGuO1/4IEH1L59ezkcDm3fvl0zZ87Ut99+q/T0dElSbm6uwsPDz3u98PBw5ebm1vq7UlJS9Mwzz3h+EgAAoEGq1xDzz3/+Uw888IACAgLc9k+aNMn1c1xcnK6//nr16dNHmzdvVq9evSRVLxA+l2EYte6XpJkzZ2rq1Kmux4WFhYqOjvbENAAAQANUbyHmq6++0p49e/Thhx9ecmyvXr3k7++vzMxM9erVSw6HQ0eOHDlv3NGjRxUREVHra9hsNtlstquuGwAAmEO9rYmZO3euevfurZ49e15y7I4dO1ReXq7IyEhJUv/+/eV0OvXNN9+4xqxfv15Op1MDBgyor5IBAICJ1PlKTHFxsfbt2+d6nJWVpYyMDLVs2VJt27aVVP1Rzr///W+9/PLL5z3/u+++07x583T77bcrLCxMO3fu1LRp05SQkKCbbrpJktS1a1eNGDFCkyZNct16/dBDD2nkyJHcmQQAACRdwZWYjRs3KiEhQQkJCZKkqVOnKiEhQU8//bRrzPz582UYhn7605+e93yr1aply5Zp+PDhio2N1aOPPqrExEQtXbpUvr6+rnHz5s1TfHy8EhMTlZiYqB49eui99967kjkCAIBGyGIYhuHtIupDYWGh7Ha7nE6nQkJCPPa6J8sq1O3pLyRJO58drkBrva6NBgCgSanL+ze9kwAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCnVOcSsWrVKo0aNUlRUlCwWixYtWuR2fOLEibJYLG5bv3793MaUlpZq8uTJCgsLU1BQkEaPHq1Dhw65jcnPz1dycrLsdrvsdruSk5NVUFBQ5wkCAIDGqc4hpqSkRD179tScOXMuOGbEiBHKyclxbZ999pnb8SlTpig1NVXz58/X6tWrVVxcrJEjR6qystI1Zvz48crIyFBaWprS0tKUkZGh5OTkupYLAAAaKb+6PiEpKUlJSUkXHWOz2eRwOGo95nQ6NXfuXL333nsaNmyYJOn9999XdHS0li5dquHDh2vXrl1KS0vTunXr1LdvX0nSm2++qf79+2vPnj2KjY2ta9kAAKCRqZc1MStWrFB4eLg6d+6sSZMmKS8vz3Vs06ZNKi8vV2JiomtfVFSU4uLitGbNGknS2rVrZbfbXQFGkvr16ye73e4aAwAAmrY6X4m5lKSkJN17772KiYlRVlaWnnrqKQ0ZMkSbNm2SzWZTbm6urFarQkND3Z4XERGh3NxcSVJubq7Cw8PPe+3w8HDXmHOVlpaqtLTU9biwsNCDswIAAA2Nx0PMfffd5/o5Li5Offr0UUxMjBYvXqwxY8Zc8HmGYchisbgen/3zhcacLSUlRc8888xVVA4AAMyk3m+xjoyMVExMjDIzMyVJDodDZWVlys/PdxuXl5eniIgI15gjR46c91pHjx51jTnXzJkz5XQ6XVt2draHZwIAABqSeg8xx48fV3Z2tiIjIyVJvXv3lr+/v9LT011jcnJytH37dg0YMECS1L9/fzmdTn3zzTeuMevXr5fT6XSNOZfNZlNISIjbBgAAGq86f5xUXFysffv2uR5nZWUpIyNDLVu2VMuWLTVr1iyNHTtWkZGROnDggJ544gmFhYXp7rvvliTZ7XY9+OCDmjZtmlq1aqWWLVtq+vTpio+Pd92t1LVrV40YMUKTJk3SG2+8IUl66KGHNHLkSO5MAgAAkq4gxGzcuFGDBw92PZ46daokacKECXr99de1bds2vfvuuyooKFBkZKQGDx6sDz/8UMHBwa7nzJ49W35+fho3bpxOnTqloUOH6u2335avr69rzLx58/Too4+67mIaPXr0Rb+bBgAANC0WwzAMbxdRHwoLC2W32+V0Oj360dLJsgp1e/oLSdLOZ4cr0OrxtdEAADRZdXn/pncSAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUIMAAAwJUJMHa357ri3SwAAACLE1JnV98c/snX7CTQAAHgLIaaO+rQLdf38h4936FRZpRerAQCg6SLEXIXsE6f052V7vV0GAABNEiHmKv3jqyxt/8Hp7TIAAGhyCDFXYUScQ5VVhh77z1aVV1Z5uxwAAJoUQsxVeOL2LrI389fOnELNXZ3l7XIAAGhSCDFXIay5Tb+/o6skaXb6XmUdK/FyRQAANB2EmKt0T+82urlTmEorqjRz4VYZhuHtkgAAaBIIMVfJYrHo+bvjFeDvo3X7T+jDDdneLgkAgCaBEOMBbVsFatptsZKk5z7bpbzC016uCACAxo8Q4yE/v6mderSxq+h0hZ7+aIe3ywEAoNEjxHiIn6+PXhjTQ34+FqXtyFXa9hxvlwQAQKNGiPGgblEh+tXADpKkpz/aIeepci9XBABA40WI8bDJQ65Xh7Ag5RWV6oXPd3m7HAAAGi1CjIcF+PsqZUy8JOlf32Rr7Xd0ugYAoD7UOcSsWrVKo0aNUlRUlCwWixYtWuQ6Vl5erscff1zx8fEKCgpSVFSU/uu//kuHDx92e41BgwbJYrG4bffff7/bmPz8fCUnJ8tut8tutys5OVkFBQVXNMlrrW+HVhrft60kaebCrTpdTqdrAAA8rc4hpqSkRD179tScOXPOO3by5Elt3rxZTz31lDZv3qyFCxdq7969Gj169HljJ02apJycHNf2xhtvuB0fP368MjIylJaWprS0NGVkZCg5Obmu5XrNjKQuigix6cDxk/rz0kxvlwMAQKPjV9cnJCUlKSkpqdZjdrtd6enpbvteffVV3XjjjTp48KDatm3r2h8YGCiHw1Hr6+zatUtpaWlat26d+vbtK0l688031b9/f+3Zs0exsbF1LfuaCwnw1x/vjNND723Sm1/t18gekYq7zu7tsgAAaDTqfU2M0+mUxWJRixYt3PbPmzdPYWFh6t69u6ZPn66ioiLXsbVr18put7sCjCT169dPdrtda9asqfX3lJaWqrCw0G3ztsTuDt0eX93pesbCraqg0zUAAB5TryHm9OnTmjFjhsaPH6+QkBDX/gceeED/+te/tGLFCj311FNasGCBxowZ4zqem5ur8PDw814vPDxcubm5tf6ulJQU1/oZu92u6Ohoz0/oCswa3V0hAX7a/gOdrgEA8KR6CzHl5eW6//77VVVVpddee83t2KRJkzRs2DDFxcXp/vvv13/+8x8tXbpUmzdvdo2xWCznvaZhGLXul6SZM2fK6XS6tuzshtHDKDw4QL+/o5sk6ZX0vTpAp2sAADyiXkJMeXm5xo0bp6ysLKWnp7tdhalNr1695O/vr8zM6gWwDodDR44cOW/c0aNHFRERUetr2Gw2hYSEuG0Nxb192mhAx1YqrajSE6nb6HQNAIAHeDzEnAkwmZmZWrp0qVq1anXJ5+zYsUPl5eWKjIyUJPXv319Op1PffPONa8z69evldDo1YMAAT5dc7ywWi1LGVHe6XvPdcf174yFvlwQAgOnVOcQUFxcrIyNDGRkZkqSsrCxlZGTo4MGDqqio0D333KONGzdq3rx5qqysVG5urnJzc1VWViZJ+u677/Tss89q48aNOnDggD777DPde++9SkhI0E033SRJ6tq1q0aMGKFJkyZp3bp1WrdunSZNmqSRI0ea4s6k2sS0CtLU2zpLkv5n8U46XQMAcJXqHGI2btyohIQEJSQkSJKmTp2qhIQEPf300zp06JA+/vhjHTp0SD/5yU8UGRnp2s7cVWS1WrVs2TINHz5csbGxevTRR5WYmKilS5fK19fX9XvmzZun+Ph4JSYmKjExUT169NB7773noWl7xy9uaq/46+wqPF2hWZ/Q6RoAgKthMRrpAo3CwkLZ7XY5nU6Pro85WVahbk9/IUna+exwBVrr9lU7Ow47NXrO16qsMvRGcm8N7177d+UAANAU1eX9m95J11j3KLseurW60/VTi7bT6RoAgCtEiPGC3w69Xu1dna53e7scAABMiRDjBe6drg9q3X46XQMAUFeEGC/p16GVfnrjmU7X2+h0DQBAHRFivGhGUheFB9uUdaxE/7eMTtcAANQFIcaL7M389ce74iRJb6zarx2HnV6uCAAA8yDEeNnw7g4lxdV0ul6wjU7XAABcJkJMA/BMTafrbT849c+v6XQNAMDlIMQ0AOEhAXryjq6Sqjtdf3+cTtcAAFwKIaaBGNcnWv07tNLpcjpdAwBwOQgxDcSZTtc2Px99ve+4/r2JTtcAAFwMIaYBaRcWpN/VdLp+bvEu5RXR6RoAgAshxDQwv7y5vbpHhch5qlzPfLzT2+UAANBgEWIaGD9fH704tod8fSxavC1HS3bkerskAAAaJEJMAxR3nV2TbqnpdP3RdhWeptM1AADnIsQ0UFOGXa92rQJ1pLBUL9LpGgCA8xBiGqgAf189X9Ppet76g1pPp2sAANwQYhqwAR3DdP8N0ZLodA0AwLkIMQ3czNu7qnWwTfuPlejVL+l0DQDAGYSYBs7ezF9/vLO7JOmNlfu183ChlysCAKBhIMSYwIi4SI3o7lBFlaEZC7fS6RoAABFiTOOZO7srOMBPWw859faaA94uBwAAryPEmERESICevL260/X/Ltmjg8dPerkiAAC8ixBjIvfdEK1+HVrS6RoAABFiTKW603UP2fx8tHrfMS3Y/IO3SwIAwGsIMSbTPixIU4ZVd7r+46c7dbSo1MsVAQDgHYQYE/rlLe3VLbKm0/UnO7xdDgAAXkGIMSF/Xx+9dE91p+tPt+Zo6c4j3i4JAIBrjhBjUnHX2fXLm9tLkn6/aLuK6HQNAGhiCDEmNmVYZ8W0ClRu4Wm9mEanawBA00KIMbFmVl+l3F3d6fr9dQe14cAJL1cEAMC1Q4gxuQGdwnRfn+pO148v2EqnawBAk0GIaQSeONPp+miJ/rp8n7fLAQDgmiDENAL2QH89O7q60/XrK77Trhw6XQMAGj9CTCMxIs6hxG4R1Z2uF2xVZRUtCQAAjRshppGwWCz6411xCrb56dtDTr31dZa3SwIAoF4RYhqRiJAAzazpdP3ykr3KPkGnawBA41XnELNq1SqNGjVKUVFRslgsWrRokdtxwzA0a9YsRUVFqVmzZho0aJB27HD/avzS0lJNnjxZYWFhCgoK0ujRo3Xo0CG3Mfn5+UpOTpbdbpfdbldycrIKCgrqPMGm5v4botW3fUudKq+k0zUAoFGrc4gpKSlRz549NWfOnFqPv/TSS3rllVc0Z84cbdiwQQ6HQ7fddpuKiopcY6ZMmaLU1FTNnz9fq1evVnFxsUaOHKnKyh9vDx4/frwyMjKUlpamtLQ0ZWRkKDk5+Qqm2LT4+FiUMiZeVj8ffZV5TAvpdA0AaKyMqyDJSE1NdT2uqqoyHA6H8cILL7j2nT592rDb7cbf/vY3wzAMo6CgwPD39zfmz5/vGvPDDz8YPj4+RlpammEYhrFz505DkrFu3TrXmLVr1xqSjN27d19WbU6n05BkOJ3Oq5nieUpKy42Yxz81Yh7/1CgpLffoa3vSX5dnGjGPf2r0fOYL42jRaW+XAwDAZanL+7dH18RkZWUpNzdXiYmJrn02m00DBw7UmjVrJEmbNm1SeXm525ioqCjFxcW5xqxdu1Z2u119+/Z1jenXr5/sdrtrzLlKS0tVWFjotjVlk27poK6RISo4Wa5nPtnp7XIAAPA4j4aY3NxcSVJERITb/oiICNex3NxcWa1WhYaGXnRMeHj4ea8fHh7uGnOulJQU1/oZu92u6Ojoq56Pmfn7+ujFsfHysUiffHtYy3bR6RoA0LjUy91JFovF7bFhGOftO9e5Y2obf7HXmTlzppxOp2vLzs6+gsoblx5tWuiXt3SQRKdrAEDj49EQ43A4JOm8qyV5eXmuqzMOh0NlZWXKz8+/6JgjR86/cnD06NHzrvKcYbPZFBIS4rZB+t2wzmrbMlA5ztP60xd7vF0OAAAe49EQ0759ezkcDqWnp7v2lZWVaeXKlRowYIAkqXfv3vL393cbk5OTo+3bt7vG9O/fX06nU998841rzPr16+V0Ol1jcHmaWX2VMqa60/V7677XRjpdAwAaiTqHmOLiYmVkZCgjI0NS9WLejIwMHTx4UBaLRVOmTNHzzz+v1NRUbd++XRMnTlRgYKDGjx8vSbLb7XrwwQc1bdo0LVu2TFu2bNHPfvYzxcfHa9iwYZKkrl27asSIEZo0aZLWrVundevWadKkSRo5cqRiY2M9N/sm4qZOYbq3dxsZRnWn69IKOl0DAMzPr65P2LhxowYPHux6PHXqVEnShAkT9Pbbb+uxxx7TqVOn9Jvf/Eb5+fnq27evlixZouDgYNdzZs+eLT8/P40bN06nTp3S0KFD9fbbb8vX19c1Zt68eXr00UdddzGNHj36gt9Ng0t78o6uWr7nqL47WqK/frlPUxMJgwAAc7MYRuP8StfCwkLZ7XY5nU6Pro85WVahbk9/IUna+exwBVrrnAO9ZvHWHD38wWb5+Vj06aM3q4uDdUMAgIalLu/f9E5qQm6Pd+g2V6frbXS6BgCYGiGmCbFYLPrjndWdrjOyC/TOmgPeLgkAgCtGiGliHPYAzbi9iyTpT1/sodM1AMC0CDFN0E9vaKsbazpdP7loO52uAQCmRIhpgs7udL1q71EtyqDTNQDAfAgxTVTH1s3126HXS5Ke/WSnjheXerkiAADqhhDThD10awd1cQQr/2S5nv2UTtcAAHMhxDRh/r4+eumeHvKxSB9lHNaXu+l0DQAwD0JME9ejTQs9eHN7SdLvU7eruLTCyxUBAHB5CDHQ727rrOiWzXTYeVp/Stvt7XIAALgshBgo0OqnlLt7SJLeXfe9Nn1Pp2sAQMNHiIEk6ebrw3SPq9P1NjpdAwAaPEIMXH5/R1eFNbdqX16xXlv+nbfLAQDgoggxcGkRaNWs0d0lSa+t2Ke9R4q8XBEAABdGiIGbO+IjNaxruMorDT32n610ugYANFiEGLixWCz6411xal7T6frdtQe8XRIAALUixOA8kfZmejzpx07Xh/LpdA0AaHgIMajVAze21Q3tQnWyrFJPptLpGgDQ8BBiUKvqTtc9ZPX10cq9R/VRxmFvlwQAgBtCDC6oU3hzPTq0kyTpmU920OkaANCgEGJwUQ/d2tHV6fqPdLoGADQghBhclNXPRy+Mre50vSjjsJbvyfN2SQAASCLE4DL8JLqFfn4Tna4BAA0LIQaXZVpiZ7UJbaYfCk7pf7/Y4+1yAAAgxODyBFr99Pzd8ZKkd9Ye0Kbv871cEQCgqSPE4LLd2rm1xvaq7nQ9Y8FWlVVUebskAEATRohBnfz+jq5qFWRVZl6xXluxz9vlAACaMEIM6iQ0yKo/1HS6/uvyfcqk0zUAwEsIMaizUT0iNbRLdafrxxfQ6RoA4B2EmDoKtPrpwAt36MALdyjQ6uftcrzi7E7Xmw8W6P1133u7JABAE0SIwRWJatFMj4+IlSS9lLZbPxSc8nJFAICmhhCDK/ZA3xj1iQlVSVmlfp+6jU7XAIBrihCDK+bjY9ELY+Nl9fXR8j1H9fG3dLoGAFw7hBhclU7hwXpkyJlO1zt1oqTMyxUBAJoKQgyu2q8HdlRsRLBOlJTpf+h0DQC4RggxuGrVna7jZbFIC7f8oBV0ugYAXAOEGHhEQttQ/XxAdafrJ1O3q4RO1wCAeubxENOuXTtZLJbztocffliSNHHixPOO9evXz+01SktLNXnyZIWFhSkoKEijR4/WoUOHPF0qPGxaYmdd16Km0/USOl0DAOqXx0PMhg0blJOT49rS09MlSffee69rzIgRI9zGfPbZZ26vMWXKFKWmpmr+/PlavXq1iouLNXLkSFVWVnq6XHhQkM1Pz4+p7nT99poD2nKQTtcAgPrj8RDTunVrORwO1/bpp5+qY8eOGjhwoGuMzWZzG9OyZUvXMafTqblz5+rll1/WsGHDlJCQoPfff1/btm3T0qVLPV0uPGxg59Yak3BdTafrbXS6BgDUm3pdE1NWVqb3339fv/jFL2SxWFz7V6xYofDwcHXu3FmTJk1SXt6PC0E3bdqk8vJyJSYmuvZFRUUpLi5Oa9asueDvKi0tVWFhodsG73hqZDe1CrJqz5Ei/W3ld94uBwDQSNVriFm0aJEKCgo0ceJE176kpCTNmzdPX375pV5++WVt2LBBQ4YMUWlpqSQpNzdXVqtVoaGhbq8VERGh3NzcC/6ulJQU2e121xYdHV0vc8KlhQZZ9fSobpKkOV/u0748Ol0DADyvXkPM3LlzlZSUpKioKNe+++67T3fccYfi4uI0atQoff7559q7d68WL1580dcyDMPtas65Zs6cKafT6dqys7M9Ng/U3eieURrSJVxllVV6fME2VdHpGgDgYfUWYr7//nstXbpUv/zlLy86LjIyUjExMcrMzJQkORwOlZWVKT/ffVFoXl6eIiIiLvg6NptNISEhbhu850yn6yCrrzZ9n6/319PpGgDgWfUWYt566y2Fh4frjjvuuOi448ePKzs7W5GRkZKk3r17y9/f33VXkyTl5ORo+/btGjBgQH2Vi3pwXYtmemxEF0nSi5/v1mE6XQMAPKheQkxVVZXeeustTZgwQX5+fq79xcXFmj59utauXasDBw5oxYoVGjVqlMLCwnT33XdLkux2ux588EFNmzZNy5Yt05YtW/Szn/1M8fHxGjZsWH2Ui3qU3C9Gvdq2qO50vWg7na4BAB5TLyFm6dKlOnjwoH7xi1+47ff19dW2bdt05513qnPnzpowYYI6d+6stWvXKjg42DVu9uzZuuuuuzRu3DjddNNNCgwM1CeffCJfX9/6KBf1yMfHohfH9pDV10df7s7TJ1tzvF0SAKCRsBiN9H+NCwsLZbfb5XQ6WR/TAPxlaaZmL92rVkFWLZ06UKFBVm+XBABogOry/k3vJFwT/z2oozpHNNfxkjL9cTGdrgEAV48Qg2uiutN1j+pO15t/0Kq9R71dEgDA5AgxuGZ6tQ3VhP7tJElPpG6j0zUA4KoQYnBN/b/hsbquRTMdyj+lV9L3erscAICJEWJwTQXZ/PTc3XGSpLe+zlJGdoF3CwIAmBYhBtfcoNhw3Z1wnaoMacaCrXS6BgBcEUIMvOKpkd3UMsiq3blFeoNO1wCAK0CIgVe0DLLqDzWdrl/9cp/25RV7uSIAgNkQYuA1o3tGaVBsa5VVVmnmwq10ugYA1AkhBl5jsVj03N3xCrT6asOBfM375qC3SwIAmAghBl51XYtmemx4rKTqTtc5TjpdAwAuDyEGXpfcv516tW2h4tIK/T6VTtcAgMtDiIHX+fpY9MLYHvL3tWjZ7jx9SqdrAMBlIMSgQegcEazfDOokSZr18Q7ll5R5uSIAQENHiEGD8ZvBHXV9eHWn6/9ZvMvb5QAAGjhCDBoMm5+vq9P1gs2H9FUmna4BABdGiEGD0jvGvdP1yTI6XQMAakeIQYMzfXisouwByj5xSq8sodM1AKB2hBg0OM1tfnpuTLwk6Z9fZ+lbOl0DAGpBiEGDNDg2XHf+JEpVhvT4gq0qr6TTNQDAHSEGDdbTI7spNNBfu3OL9PdV+71dDgCggSHEoMFq1dymp2s6Xf9lWaa+O0qnawDAjwgxaNDu+sl1Gti5tcoqqjRzwbar7nR9sqxC7WYsVrsZi7nzCQBMjhCDBq2603WcAq2++ubACX1Ap2sAQA1CDBq8NqGB+n81na5f+Hy3cp2nvVwRAKAhIMTAFP6rfzv9JLqm0/UiOl0DAAgxMAlfH4teuqe60/XSXUf02bZcb5cEAPAyQgxMo3NEsP67ptP1Hz7eroKTdLoGgKaMEANTeXhwR3UKb65jxWV6jk7XANCkEWJgKjY/X704Nl4Wi/TvTYe0OvOYt0sCAHgJIQam0zumpZL7xUiSZqZu1amySi9XBADwBkIMTOmxEV0UWdPpevZSOl0DQFNEiIEpNbf56bm74yRJ//hqv7YeKvBuQQCAa44QA9Ma0iVCo3ue6XS9jU7XANDEEGJgak+P6qYWgf7alVNIp2sAaGIIMTC1sOY2PT3yx07X++l0DQBNBiEGpnd3wnW65fowlVVUacbCq+90DQAwB4+HmFmzZslisbhtDofDddwwDM2aNUtRUVFq1qyZBg0apB07dri9RmlpqSZPnqywsDAFBQVp9OjROnTokKdLRSNhsVj0/N3xaubvq2+yTmj+hmxvlwQAuAbq5UpM9+7dlZOT49q2bdvmOvbSSy/plVde0Zw5c7RhwwY5HA7ddtttKioqco2ZMmWKUlNTNX/+fK1evVrFxcUaOXKkKiv5PhDULrploKbXdLpO+WwXna4BoAmolxDj5+cnh8Ph2lq3bi2p+irMn//8Zz355JMaM2aM4uLi9M477+jkyZP64IMPJElOp1Nz587Vyy+/rGHDhikhIUHvv/++tm3bpqVLl9ZHuWgkJg5op57RLVRUWqGnPqLTNQA0dvUSYjIzMxUVFaX27dvr/vvv1/791XeNZGVlKTc3V4mJia6xNptNAwcO1Jo1ayRJmzZtUnl5uduYqKgoxcXFucbUprS0VIWFhW4bmhZfH4teHBsvPx+L0nce0efb6XQNAI2Zx0NM37599e677+qLL77Qm2++qdzcXA0YMEDHjx9Xbm71m0pERITbcyIiIlzHcnNzZbVaFRoaesExtUlJSZHdbndt0dHRHp4ZzKCLI0S/GdRRkvT0RzvkPFnu5YoAAPXF4yEmKSlJY8eOVXx8vIYNG6bFixdLkt555x3XGIvF4vYcwzDO23euS42ZOXOmnE6na8vOZnFnU/XwkE7q2DpIx4pL9dxnO71dDgCgntT7LdZBQUGKj49XZmam6y6lc6+o5OXlua7OOBwOlZWVKT8//4JjamOz2RQSEuK2oWmq7nTdQ5L0/208pK/30ekaABqjeg8xpaWl2rVrlyIjI9W+fXs5HA6lp6e7jpeVlWnlypUaMGCAJKl3797y9/d3G5OTk6Pt27e7xgCX0qfdWZ2uF26j0zUANEIeDzHTp0/XypUrlZWVpfXr1+uee+5RYWGhJkyYIIvFoilTpuj5559Xamqqtm/frokTJyowMFDjx4+XJNntdj344IOaNm2ali1bpi1btuhnP/uZ6+Mp4HI9NiJWkfYAHTxxUn+m0zUANDp+nn7BQ4cO6ac//amOHTum1q1bq1+/flq3bp1iYqr/r/ixxx7TqVOn9Jvf/Eb5+fnq27evlixZouDgYNdrzJ49W35+fho3bpxOnTqloUOH6u2335avr6+ny0UjFhzgr/+5K04PvrNRb361X6N6RqlD6yBvlwUA8BCL0Ui/TKOwsFB2u11Op5P1MU3cIx9s1qdbc9QtMkT/eqivej5T/VHlzmeHK9Dq8RwPALgKdXn/pncSGr1Zo7urRaC/duYU6p01B7xdDgDAQwgxaPTCmtv0+zuqO13/dfl3Xq4GAOAphBg0CWN7VXe6Lq2o8nYpAAAPIcSgSTi70zUAoHEgxKDJiG4ZqEeHdnI93vx9vqqqGuW6dgBoErg7CU1K0elyxc9a4nrcKsiqgbGtNaRLuG65vrXszfy9WB0AoC7v39xfiibF1+fH/lvNbX46XlKmhZt/0MLNP8jXx6I+MaEa0iVcQ7qEq1N480v29AIAeA8hBk3W6scHa2dOoZbvztOXu/P03dESrc86ofVZJ5Ty+W61CW2mwbHVgaZ/x1YKYD0NADQofJyEJuVkWYW6Pf2FpPO/7O7g8ZNavqc60Kzdf1xlZ93JFODvowEdwzS45irNdS2aXfPaAaAp4OMk4Aq0bRWoCQPaacKAdjpZVqE1+47ryz15Wr47TznO0/qy5orNU5I6RzSvDjSx4eodEyo/X9bIA8C1RogBahFo9dOwbhEa1i1ChmFoz5Eifbm7OtBs+j5fe48Ua++RYr2xcr9CAvx0a+fqxcEDO7dWq+Y2b5cPAE0CIQa4BIvFoi6OEHVxhOg3gzqp4GSZVu49quW787Ry71HlnyzXp1tz9OnWHFksUs82LVyLg7tHhbA4GADqCSEGqKMWgVbd+ZPrdOdPrlNllaGM7ALX4uCdOYXKyC5QRnaBXknfq/BgmwbHhmtwl3DdfH2Ymtv4Tw4APIW/UYGr4OtjUe+YUPWOCdX04bHKdZ52LQ7+et8x5RWV6sON2fpwY7b8fS26sX1L1x1PHVo393b5AGBqhBjAgxz2AP30xrb66Y1tVVpRqW+yTrjW0hw4flJf7zuur/cd1/8s3qV2rQI1uEu4BseGq2+HlrL5cQs3ANQFt1ijSbnYLdb1bf/R4upAsydP32SdUHnlj//pBVp9dVOnMA2pCTUOe8A1qwsAGhJusQYaoA6tm6tD6+b65S0dVFxaodWZx7S8JtTkFZUqfecRpe88IknqFhlSHWi6tNZPokPdvmkYAFCNEAN4QXObn0bEOTQizqGqKkM7cwpd30Pz7aEC7cwp1M6cQs1Zvk+hgf4a2Lm1Btfcwt0i0Ort8gGgQSDEAF7m42NR3HV2xV1n16NDr9fx4lKt3HtUX+7O06qaW7gXZRzWoozD8rFIvWNCXWtpujiCm9wt3N78SBBAw8J//UAD06q5TWN6tdGYXm1UUVmlTd/nu745eO+RYm04kK8NB/L1UtoeRdkDNKjmm4MHdGrFGzqAJoW/8YAGzM/XR307tFLfDq00M6mrDuWf1PI91V+09/W+YzrsPK0P1h/UB+sPyurno/4dWrkWB7dtFejt8gGgXhFiABNpExqo5H4xSu4Xo9PllVr73XHXWpofCk5p5d6jWrn3qP6gHerYOqhmcXC4bmjXUv70dwLQyBBiAJMK8PetXhvTJVzPGoYy84pd3xy88ft8fXe0RN8dzdKbX2Up2Oanm68Pc62laR1MfycA5keIARoBi8WizhHB6hwRrF8N7CjnqXJ9lVm9OHjlnqM6XlKmz7fn6vPtuZKkHm3srm8Ojr/OLh9u4cY1xOJseAr/5gCNkL2Zv0b2iNLIHlGqqjL07aEC11qabT84tfVQ9faXZZkKa27VwM7VgeaWzmEKCfD3dvkAcFkIMUAj5+NjUULbUCW0DdXU2zorr/C0Vuypvkqzet8xHSsu04LNh7Rg8yH5+VjUp12oqwt3x9bNm9wt3ADMgxADNDHhIQEad0O0xt0QrbKKKm04cKJ6Lc2ePO0/WqJ1+09o3f4Tev6z3Ypu2czVhbt/h1YK8Ke/E4CGgxADNGFWPx/d1ClMN3UK0+9HdtP3x0tcdzut339C2SdO6d213+vdtd8rwN9HN3UMcy0mvq5FM2+XD6CJI8QAcIlpFaSf39ReP7+pvUpKK/T1vmOutTS5hae1bHeelu3OkyR1cQRrUM3i4F5tW8iPW7iBJqOhLM4mxACoVZDNT4ndHUrs7pBhGNqVU6Tle6qv0mw5mK/duUXanVukv638TvZm/rq1c2sN6dJaAzuHq2UQ/Z0A1D9CDIBLslgs6hYVom5RIXp4cCfll5Rp5d6jWr4nTyv3HlXByXJ98u1hffLtYVksUkJ0C9damu5RISwOBlAvCDEA6iw0yKq7Eq7TXQnXqaKyShnZBa61NLtzi7T5YIE2HyzQy+l7FRFicwWamzuFKcjGXzsAPIO/TQBcFT9fH/Vp11J92rXUYyO66HDBKdct3F/vO6YjhaWavyFb8zdky+rro74dWrrW0rQPC/J2+QBMjBADwKOiWjTT+L5tNb5vW50ur9T6rBOudggHT5zUV5nH9FXmMf3x051qHxbk+ubgG9qHyubHLdwALh8hBkC9CfD31cDOrTWwc2v9YVQ37T9W4go032SdUNaxEmUdy9I/v85SkNVXN3UKczWtjAgJ8Hb5ABo4QgyAa8Jisahj6+bq2Lq5fnlLBxWdLtfqzGP6cneelu85qmPFpVqy84iW7DwiSeoeFeIKND3btJAv/Z0AnIMQA8ArggP8lRQfqaT4SFVVGdpxuLB6cfCePG09VKAdhwu143ChXv1yn1oGWTWwc2sN7hKuG9qFert0AA2Ex0NMSkqKFi5cqN27d6tZs2YaMGCAXnzxRcXGxrrGTJw4Ue+8847b8/r27at169a5HpeWlmr69On617/+pVOnTmno0KF67bXX1KZNG0+XjCYk0OqnAy/c4e0ycA4fH4vi29gV38au3w67XseKS7Wi5kv2Vu09qhMlZUrd8oNSt/zgdkXmH1/tV3Obv2z+PrL5+crq5yOba/Ot2X/W47P2W3196N4NmJzHQ8zKlSv18MMP64YbblBFRYWefPJJJSYmaufOnQoK+vFOhBEjRuitt95yPbZa3b8ca8qUKfrkk080f/58tWrVStOmTdPIkSO1adMm+fqy+A9ozMKa23RP7za6p3cblVdWadP3+a61NJl5xa5xr6RnXtXv8fe1nBVufGTz95XV1+es8FNzrCb0uAej6mPWs557Zr9bmPL/8fWt54Qrq68P36EDXAWPh5i0tDS3x2+99ZbCw8O1adMm3Xrrra79NptNDoej1tdwOp2aO3eu3nvvPQ0bNkyS9P777ys6OlpLly7V8OHDPV02gAbK39dH/Tq0Ur8OrTTz9q7ae6RQibO/kiTdlRClyiqprKJSpRVVKi2vUumZnyuqVFZx1uOaY1XGj69dXmmovLJCxaVempzkHm78fOp4VenHkGQ995i/j6y+54cuVyCreZ6fj4UgBdOq9zUxTqdTktSyZUu3/StWrFB4eLhatGihgQMH6rnnnlN4eLgkadOmTSovL1diYqJrfFRUlOLi4rRmzZpaQ0xpaalKS3/8m6iwsLA+pgPAy9qEBrp+fv7u+Dr1bDEMQxVVRk24qQk45T/+XNt+tyBUUaXS8rN+riUklZ17rLxKZZXuzzvbmX1FqvDYn1Fd+Fjk9hHbmRB13tWjC11V8nN/jrXW4OXrFsCqjB//DKrOTpVAHdVriDEMQ1OnTtXNN9+suLg41/6kpCTde++9iomJUVZWlp566ikNGTJEmzZtks1mU25urqxWq0JD3RfwRUREKDc3t9bflZKSomeeeaY+pwPA5CwWi/x9LfL39VGQzTs1GIZRHWpcAemcYFReWRN63ENSWS2hy+01LhjIfnz9M/vKKs8KEYZ0qrxSp8orvfLnETdrSc1HddUhKaAm+ATUBKYAf99zfj47IFU/Dqh5fPY/A/xr2+cerrgCZX71GmIeeeQRbd26VatXr3bbf99997l+jouLU58+fRQTE6PFixdrzJgxF3w9wzAu+C/dzJkzNXXqVNfjwsJCRUdHX+UMAMCzLJYz63C8t7avqso4KyhV1umq0vlXp86+2nT+6/0Y1NyPVZ51Baassvr5RaXX9mrUuWHp3KB0duipddxZwevcIGXzcw9fZ475+/LxnSfVW4iZPHmyPv74Y61ateqSdxRFRkYqJiZGmZnVi/QcDofKysqUn5/vdjUmLy9PAwYMqPU1bDabbDYv/a8VAJiIj49FAT7Vb7KS/zX//SfLKtTt6S8kSV8/PlgWi0Wna64Wnf3P0+Vnh6fqx27Ha46drrnadPqs55e6vc6PP5/96VVtH+/VtzMf350dlM6EHVutV59queLkFqgucMWpJmSdGevn63NN53mteDzEGIahyZMnKzU1VStWrFD79u0v+Zzjx48rOztbkZGRkqTevXvL399f6enpGjdunCQpJydH27dv10svveTpkgEAXhIaZK3TuqarcWZN1NkBye2ftQapMz+fFaRc4ensgHV++Dr7n2e4f3xXfk3mLUl+PpbzApLtrKBU+xUp9322s64+NZRrSR7/N+fhhx/WBx98oI8++kjBwcGuNSx2u13NmjVTcXGxZs2apbFjxyoyMlIHDhzQE088obCwMN19992usQ8++KCmTZumVq1aqWXLlpo+fbri4+NddysBAFAXZ6+JCr6GXS0MwzhvfdKFwk6tV6TOCktnwtMln1vuvvapospQRVmlSsq8s/apvng8xLz++uuSpEGDBrntf+uttzRx4kT5+vpq27Ztevfdd1VQUKDIyEgNHjxYH374oYKDg13jZ8+eLT8/P40bN871ZXdvv/023xEDADAVi8XiurKhZtfu47uqqjPh6SIBqZZ9pWd9PFdaXvvVp5NlFfr2kPOazeVC6uXjpItp1qyZvvjii0u+TkBAgF599VW9+uqrnioNAIAmw8fHomZWXzWzev5//s9e1+RNjXOlDwAAaPQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJQIMQAAwJSuTf9zAPCQQKufDrxwh7fLANAAcCUGAACYEldiAADXFFfT4ClciQEAAKbElRgAAFAnDeVqGldiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKRFiAACAKfl5u4D6YhiGJKmwsNDLlQAAgMt15n37zPv4xTTaEFNUVCRJio6O9nIlAACgroqKimS32y86xmJcTtQxoaqqKh0+fFjBwcGyWCwefe3CwkJFR0crOztbISEhHn3thoD5mV9jn2Njn5/U+OfI/MyvvuZoGIaKiooUFRUlH5+Lr3pptFdifHx81KZNm3r9HSEhIY32X06J+TUGjX2OjX1+UuOfI/Mzv/qY46WuwJzBwl4AAGBKhBgAAGBKhJgrYLPZ9Ic//EE2m83bpdQL5md+jX2OjX1+UuOfI/Mzv4Ywx0a7sBcAADRuXIkBAACmRIgBAACmRIgBAACmRIgBAACmRIi5gNdee03t27dXQECAevfura+++uqi41euXKnevXsrICBAHTp00N/+9rdrVOmVqcv8VqxYIYvFct62e/fua1jx5Vu1apVGjRqlqKgoWSwWLVq06JLPMdP5q+v8zHb+UlJSdMMNNyg4OFjh4eG66667tGfPnks+z0zn8ErmaKbz+Prrr6tHjx6uL0Hr37+/Pv/884s+x0znr67zM9O5q01KSoosFoumTJly0XHeOIeEmFp8+OGHmjJlip588klt2bJFt9xyi5KSknTw4MFax2dlZen222/XLbfcoi1btuiJJ57Qo48+qgULFlzjyi9PXed3xp49e5STk+Parr/++mtUcd2UlJSoZ8+emjNnzmWNN9v5q+v8zjDL+Vu5cqUefvhhrVu3Tunp6aqoqFBiYqJKSkou+ByzncMrmeMZZjiPbdq00QsvvKCNGzdq48aNGjJkiO68807t2LGj1vFmO391nd8ZZjh359qwYYP+/ve/q0ePHhcd57VzaOA8N954o/HrX//abV+XLl2MGTNm1Dr+scceM7p06eK271e/+pXRr1+/eqvxatR1fsuXLzckGfn5+degOs+SZKSmpl50jNnO39kuZ35mPn+GYRh5eXmGJGPlypUXHGPmc2gYlzdHs5/H0NBQ4x//+Eetx8x+/gzj4vMz67krKioyrr/+eiM9Pd0YOHCg8dvf/vaCY711DrkSc46ysjJt2rRJiYmJbvsTExO1Zs2aWp+zdu3a88YPHz5cGzduVHl5eb3VeiWuZH5nJCQkKDIyUkOHDtXy5cvrs8xrykzn72qY9fw5nU5JUsuWLS84xuzn8HLmeIbZzmNlZaXmz5+vkpIS9e/fv9YxZj5/lzO/M8x27h5++GHdcccdGjZs2CXHeuscEmLOcezYMVVWVioiIsJtf0REhHJzc2t9Tm5ubq3jKyoqdOzYsXqr9UpcyfwiIyP197//XQsWLNDChQsVGxuroUOHatWqVdei5HpnpvN3Jcx8/gzD0NSpU3XzzTcrLi7uguPMfA4vd45mO4/btm1T8+bNZbPZ9Otf/1qpqanq1q1brWPNeP7qMj+znTtJmj9/vjZv3qyUlJTLGu+tc9hou1hfLYvF4vbYMIzz9l1qfG37G4q6zC82NlaxsbGux/3791d2drb+93//V7feemu91nmtmO381YWZz98jjzyirVu3avXq1Zcca9ZzeLlzNNt5jI2NVUZGhgoKCrRgwQJNmDBBK1euvOAbvdnOX13mZ7Zzl52drd/+9rdasmSJAgICLvt53jiHXIk5R1hYmHx9fc+7KpGXl3deyjzD4XDUOt7Pz0+tWrWqt1qvxJXMrzb9+vVTZmamp8vzCjOdP08xw/mbPHmyPv74Yy1fvlxt2rS56FiznsO6zLE2Dfk8Wq1WderUSX369FFKSop69uypv/zlL7WONeP5q8v8atOQz92mTZuUl5en3r17y8/PT35+flq5cqX+7//+T35+fqqsrDzvOd46h4SYc1itVvXu3Vvp6elu+9PT0zVgwIBan9O/f//zxi9ZskR9+vSRv79/vdV6Ja5kfrXZsmWLIiMjPV2eV5jp/HlKQz5/hmHokUce0cKFC/Xll1+qffv2l3yO2c7hlcyxNg35PJ7LMAyVlpbWesxs5682F5tfbRryuRs6dKi2bdumjIwM19anTx898MADysjIkK+v73nP8do5rNdlwyY1f/58w9/f35g7d66xc+dOY8qUKUZQUJBx4MABwzAMY8aMGUZycrJr/P79+43AwEDjd7/7nbFz505j7ty5hr+/v/Gf//zHW1O4qLrOb/bs2UZqaqqxd+9eY/v27caMGTMMScaCBQu8NYWLKioqMrZs2WJs2bLFkGS88sorxpYtW4zvv//eMAzzn7+6zs9s5++///u/DbvdbqxYscLIyclxbSdPnnSNMfs5vJI5muk8zpw501i1apWRlZVlbN261XjiiScMHx8fY8mSJYZhmP/81XV+Zjp3F3Lu3UkN5RwSYi7gr3/9qxETE2NYrVajV69ebrc+TpgwwRg4cKDb+BUrVhgJCQmG1Wo12rVrZ7z++uvXuOK6qcv8XnzxRaNjx45GQECAERoaatx8883G4sWLvVD15TlzO+O524QJEwzDMP/5q+v8zHb+apubJOOtt95yjTH7ObySOZrpPP7iF79w/f3SunVrY+jQoa43eMMw//mr6/zMdO4u5NwQ01DOocUwalbeAAAAmAhrYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCkRYgAAgCn9/xZszKThz5IHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ordination, dm, cv = rpca.joint_rpca(tables=[wol2_bt, ko_bt, pathways_bt], n_test_samples = 200, sample_metadata=metadata_subset, train_test_column='train_test', n_components=3)\n",
    "corr_table = rpca.feature_correlation_table(ordination)\n",
    "plt.errorbar(x=cv.index, y=cv['mean_CV'], yerr=cv['std_CV']/2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8aeca454-258b-4923-8de4-9f6a781a668e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/projects/thdmi/basics/results/joint_rpca_dm.qza'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q2.Artifact.import_data('DistanceMatrix', dm).save('/projects/thdmi/basics/results/joint_rpca_dm.qza')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2831c5d3-1d40-4037-862b-f392b3c1bbd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiime2.plugins import diversity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "db5d9f4d-c8dc-44c9-a434-a2352805973d",
   "metadata": {},
   "outputs": [],
   "source": [
    "jrpca_dm = q2.Artifact.import_data('DistanceMatrix', dm)\n",
    "jrpca_pcoa = diversity.methods.pcoa(jrpca_dm, 3).pcoa\n",
    "jrpca_emp = emperor.visualizers.plot(jrpca_pcoa, q2.Metadata(metadata_subset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0193c9a6-927e-472e-be05-3baa4e54d0f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/projects/thdmi/basics/results/joint_rpca_plot.qzv'"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jrpca_emp.visualization.save('/projects/thdmi/basics/results/joint_rpca_plot.qzv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "b1bf8838-1297-4117-87dd-1951110a2e42",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_21126/4084200567.py:15: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  food_data.replace({np.nan: 0}, inplace=True)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "627 x 1177 <class 'biom.table.Table'> with 219444 nonzero entries (29% dense)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wol2_bt = biom.load_table('/projects/thdmi/feature-tables/wol2_filtered_table_zebra.biom')\n",
    "ko_bt = biom.load_table('/projects/thdmi/feature-tables/kegg_filtered_table_zebra.biom') \n",
    "pathways_bt = biom.load_table('/projects/thdmi/feature-tables/pathways_filtered_table_zebra.biom')\n",
    "\n",
    "foods_q2 = q2.Metadata.load('/projects/thdmi/metadata/consolidated_metadata_subset.tsv')\n",
    "foods = foods_q2.to_dataframe()\n",
    "#foods = foods.loc[foods['types_of_plants'].notna()]\n",
    "foods.rename_axis(index='SampleID', inplace=True)\n",
    "\n",
    "common_ids = foods.index.intersection(wol2_bt.ids())\n",
    "wol2_bt.filter(common_ids)\n",
    "\n",
    "collapsed_foods = foods.columns[foods.columns.str.startswith('Collapsed')]\n",
    "all_foods = foods.columns[foods.columns.str.startswith('GramsPerDay')]\n",
    "ultraprocessed = foods.columns[foods.columns.str.startswith('Processed')]\n",
    "\n",
    "food_data = foods[collapsed_foods.union(all_foods).union(ultraprocessed)]\n",
    "food_data.replace({np.nan: 0}, inplace=True)\n",
    "food_data_q2 = q2.Artifact.import_data('FeatureTable[Frequency]', food_data)\n",
    "food_data_bt = food_data_q2.view(biom.Table)\n",
    "food_data_bt.filter(common_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "7badac12-9393-466e-bd07-b0e75098ed3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_, test_ = train_test_split(foods, test_size = 0.5,\n",
    "                                 stratify=foods[['thdmi_cohort', 'age_cat']],\n",
    "                                 random_state=6)\n",
    "foods['train_test'] = 'test'\n",
    "foods.loc[train_.index, 'train_test'] = 'train'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "87da5421-cea8-411f-a850-4c94df78fb4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "wol2_table_q2 = q2.Artifact.import_data('FeatureTable[Frequency]', wol2_bt)\n",
    "food_table_q2 = q2.Artifact.import_data('FeatureTable[Frequency]', food_data_bt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "79642fad-efe9-41a3-865a-355a32bd4b9d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'food_table.qza'"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wol2_table_q2.save('wol2_table.qza')\n",
    "food_table_q2.save('food_table.qza')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "349fdaae-a1e6-4591-8b50-7f16677ecd11",
   "metadata": {},
   "outputs": [],
   "source": [
    "foods['train_test'].replace({'train': 'Train', 'test': 'Test'}, inplace=True)\n",
    "foods[['thdmi_cohort', 'age_cat', 'train_test', 'antibiotic_history', 'covid_level_of_wellbeing', 'types_of_plants', 'bmi_cat']].to_csv('mmvec_metadata.tsv', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "680cd5aa-1d86-4383-b370-92bf0c582e53",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/adilmore/gemelli/gemelli/preprocessing.py:418: RuntimeWarning: divide by zero encountered in log\n",
      "  mat = np.log(matrix_closure(mat))\n",
      "/home/adilmore/gemelli/gemelli/preprocessing.py:418: RuntimeWarning: divide by zero encountered in log\n",
      "  mat = np.log(matrix_closure(mat))\n",
      "/home/adilmore/gemelli/gemelli/rpca.py:740: FutureWarning: Passing a set as an indexer is deprecated and will raise in a future version. Use a list instead.\n",
      "  sample_metadata = sample_metadata.loc[shared_all_samples, :]\n",
      "/home/adilmore/gemelli/gemelli/rpca.py:900: RuntimeWarning: Removing 0 features(s) in table(s) but not the ordination.\n",
      "  warnings.warn('Removing %i features(s) in table(s)'\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAsklEQVR4nO3de3xU9Z3/8fdkJplcSIZ7hpAgQSMCCReDRVALyq0IoqKgxbq60l0sSo2CKNJadLeJUgXaUnV1WaHyoyBIqOuVoBhlqRUC4SqCCpJAhoCGSULCTC7n9wc4GgiXCYFzJnk9H4/zMHPOZ8Ln268y737PmXNshmEYAgAAsJAwsxsAAAA4GQEFAABYDgEFAABYDgEFAABYDgEFAABYDgEFAABYDgEFAABYDgEFAABYjsPsBhqitrZWBw4cUGxsrGw2m9ntAACAc2AYhsrKypSQkKCwsDOvkYRkQDlw4ICSkpLMbgMAADRAQUGBEhMTz1gTkgElNjZW0vEBxsXFmdwNAAA4F6WlpUpKSgp8jp9JSAaU70/rxMXFEVAAAAgx53J5BhfJAgAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGgAAAAyyGg/EiFv1qdH39bnR9/WxX+arPbAQCg2QoqoFRXV+s3v/mNkpOTFRUVpS5duujpp59WbW1toMYwDM2cOVMJCQmKiorSoEGDtH379jq/x+fzafLkyWrbtq1iYmI0evRoFRYWNs6IAABAyAsqoDz77LN66aWXNG/ePH3++eeaNWuW/vCHP+jPf/5zoGbWrFmaPXu25s2bp/Xr18vtdmvo0KEqKysL1GRkZCg7O1tLlizR2rVrVV5erlGjRqmmpqbxRgYAAEKWI5jif/zjH7r55ps1cuRISVLnzp31t7/9TRs2bJB0fPVk7ty5mjFjhsaMGSNJWrhwoeLj47V48WJNnDhRXq9X8+fP12uvvaYhQ4ZIkhYtWqSkpCStXr1aw4cPb8zxAQCAEBTUCsq1116rDz74QLt27ZIkbd68WWvXrtWNN94oSdqzZ488Ho+GDRsWeI/T6dTAgQO1bt06SVJeXp6qqqrq1CQkJCg1NTVQAwAAmregVlAee+wxeb1eXXHFFbLb7aqpqdHvf/97/fznP5ckeTweSVJ8fHyd98XHx+ubb74J1ERERKhVq1an1Hz//pP5fD75fL7A69LS0mDaBgAAISaoFZSlS5dq0aJFWrx4sTZu3KiFCxfqueee08KFC+vU2Wy2Oq8Nwzhl38nOVJOVlSWXyxXYkpKSgmkbAACEmKACyqOPPqrHH39cd955p9LS0nT33Xfr4YcfVlZWliTJ7XZL0ikrIcXFxYFVFbfbLb/fr5KSktPWnGz69Onyer2BraCgIJi2AQBAiAkqoFRUVCgsrO5b7HZ74GvGycnJcrvdysnJCRz3+/3Kzc3VgAEDJEnp6ekKDw+vU1NUVKRt27YFak7mdDoVFxdXZwMAAE1XUNeg3HTTTfr973+vTp06qUePHtq0aZNmz56t++67T9LxUzsZGRnKzMxUSkqKUlJSlJmZqejoaI0fP16S5HK5NGHCBE2ZMkVt2rRR69atNXXqVKWlpQW+1QMAAJq3oALKn//8Z/32t7/VpEmTVFxcrISEBE2cOFFPPvlkoGbatGmqrKzUpEmTVFJSon79+mnVqlWKjY0N1MyZM0cOh0Pjxo1TZWWlBg8erAULFshutzfeyAAAQMiyGYZhmN1EsEpLS+VyueT1ehv1dE+Fv1rdn3xfkrTj6eGKjggqvwEAgDMI5vObZ/EAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLCSqgdO7cWTab7ZTtgQcekCQZhqGZM2cqISFBUVFRGjRokLZv317nd/h8Pk2ePFlt27ZVTEyMRo8ercLCwsYbEQAACHlBBZT169erqKgosOXk5EiSxo4dK0maNWuWZs+erXnz5mn9+vVyu90aOnSoysrKAr8jIyND2dnZWrJkidauXavy8nKNGjVKNTU1jTgsAAAQyoIKKO3atZPb7Q5sb731li699FINHDhQhmFo7ty5mjFjhsaMGaPU1FQtXLhQFRUVWrx4sSTJ6/Vq/vz5ev755zVkyBD16dNHixYt0tatW7V69eoLMkAAABB6GnwNit/v16JFi3TffffJZrNpz5498ng8GjZsWKDG6XRq4MCBWrdunSQpLy9PVVVVdWoSEhKUmpoaqKmPz+dTaWlpnQ0AADRdDQ4oK1eu1JEjR3TvvfdKkjwejyQpPj6+Tl18fHzgmMfjUUREhFq1anXamvpkZWXJ5XIFtqSkpIa2DQAAQkCDA8r8+fM1YsQIJSQk1Nlvs9nqvDYM45R9JztbzfTp0+X1egNbQUFBQ9sGAAAhoEEB5ZtvvtHq1av1y1/+MrDP7XZL0ikrIcXFxYFVFbfbLb/fr5KSktPW1MfpdCouLq7OBgAAmq4GBZRXX31V7du318iRIwP7kpOT5Xa7A9/skY5fp5Kbm6sBAwZIktLT0xUeHl6npqioSNu2bQvUAAAAOIJ9Q21trV599VXdc889cjh+eLvNZlNGRoYyMzOVkpKilJQUZWZmKjo6WuPHj5ckuVwuTZgwQVOmTFGbNm3UunVrTZ06VWlpaRoyZEjjjQoAAIS0oAPK6tWrtW/fPt13332nHJs2bZoqKys1adIklZSUqF+/flq1apViY2MDNXPmzJHD4dC4ceNUWVmpwYMHa8GCBbLb7ec3EgAA0GTYDMMwzG4iWKWlpXK5XPJ6vY16PUqFv1rdn3xfkrTj6eGKjgg6vwEAgNMI5vObZ/EAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLIaAAAADLCTqg7N+/X7/4xS/Upk0bRUdHq3fv3srLywscNwxDM2fOVEJCgqKiojRo0CBt3769zu/w+XyaPHmy2rZtq5iYGI0ePVqFhYXnPxoAANAkBBVQSkpKdM011yg8PFzvvvuuduzYoeeff14tW7YM1MyaNUuzZ8/WvHnztH79erndbg0dOlRlZWWBmoyMDGVnZ2vJkiVau3atysvLNWrUKNXU1DTawAAAQOhyBFP87LPPKikpSa+++mpgX+fOnQM/G4ahuXPnasaMGRozZowkaeHChYqPj9fixYs1ceJEeb1ezZ8/X6+99pqGDBkiSVq0aJGSkpK0evVqDR8+vBGGBQAAQllQKyhvvvmm+vbtq7Fjx6p9+/bq06ePXnnllcDxPXv2yOPxaNiwYYF9TqdTAwcO1Lp16yRJeXl5qqqqqlOTkJCg1NTUQM3JfD6fSktL62wAAKDpCiqgfP3113rxxReVkpKi999/X/fff79+/etf669//askyePxSJLi4+PrvC8+Pj5wzOPxKCIiQq1atTptzcmysrLkcrkCW1JSUjBtAwCAEBNUQKmtrdWVV16pzMxM9enTRxMnTtS//du/6cUXX6xTZ7PZ6rw2DOOUfSc7U8306dPl9XoDW0FBQTBtAwCAEBNUQOnQoYO6d+9eZ1+3bt20b98+SZLb7ZakU1ZCiouLA6sqbrdbfr9fJSUlp605mdPpVFxcXJ0NAAA0XUEFlGuuuUZffPFFnX27du3SJZdcIklKTk6W2+1WTk5O4Ljf71dubq4GDBggSUpPT1d4eHidmqKiIm3bti1QAwAAmregvsXz8MMPa8CAAcrMzNS4ceP02Wef6eWXX9bLL78s6fipnYyMDGVmZiolJUUpKSnKzMxUdHS0xo8fL0lyuVyaMGGCpkyZojZt2qh169aaOnWq0tLSAt/qAQAAzVtQAeWqq65Sdna2pk+frqefflrJycmaO3eu7rrrrkDNtGnTVFlZqUmTJqmkpET9+vXTqlWrFBsbG6iZM2eOHA6Hxo0bp8rKSg0ePFgLFiyQ3W5vvJEBAICQZTMMwzC7iWCVlpbK5XLJ6/U26vUoFf5qdX/yfUnSjqeHKzoiqPwGAADOIJjPb57FAwAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeAAgAALIeA8iMffF4c+Pnbcp+JnQAA0LwRUH7k2svaBH5+8s3tMgzDxG4AAGi+CCg/4gy3B35es/OQlq4vMLEbAACaLwLKGTz91g598+1Rs9sAAKDZIaCcxlWdW6nCX6OHl+aruqbW7HYAAGhWCCinkTkmTbFOhzbuO6KXcr8yux0AAJoVAsppdGwZpadu7iFJmrt6t7YUHjG3IQAAmhECyhnc2qejbkxzq7rW0MNL81XprzG7JQAAmgUCyhnYbDb9/pY0tY916qtDR/XMu5+b3RIAAM0CAeUsWsVE6A9je0mSFv7jG+XuOmRyRwAANH0ElHMw8PJ2uqf/JZKkR5dtVslRv8kdAQDQtBFQztHjI7qpS7sYFZf5NGPlVu4yCwDABURAOUdREXbNvaO3HGE2vbPVo+xN+81uCQCAJouAEoSeiS310OAUSdLv/r5dhSUVJncEAEDTREAJ0q8GXao+nVqqzFetKa9vVm0tp3oAAGhsBJQgOexhmjOut6Ij7Prnnu/032u/NrslAACaHAJKA3RuG6PfjuouSXru/V36vKjU5I4AAGhaCCgNdOdVSRrSrb38NbV6eGm+jlVxl1kAABpLUAFl5syZstlsdTa32x04bhiGZs6cqYSEBEVFRWnQoEHavn17nd/h8/k0efJktW3bVjExMRo9erQKCwsbZzQXkc1mU9aYnmoTE6GdnjLNztlldksAADQZQa+g9OjRQ0VFRYFt69atgWOzZs3S7NmzNW/ePK1fv15ut1tDhw5VWVlZoCYjI0PZ2dlasmSJ1q5dq/Lyco0aNUo1NaG3AtEu1qlnbuspSXrlk6/1j6++NbkjAACahqADisPhkNvtDmzt2rWTdHz1ZO7cuZoxY4bGjBmj1NRULVy4UBUVFVq8eLEkyev1av78+Xr++ec1ZMgQ9enTR4sWLdLWrVu1evXqxh3ZRTK0e7zuvCpJhiFNeT1fpceqzG4JAICQF3RA2b17txISEpScnKw777xTX399/Fsse/bskcfj0bBhwwK1TqdTAwcO1Lp16yRJeXl5qqqqqlOTkJCg1NTUQE0o+u2o7urUOloHvMf0u79vP/sbAADAGQUVUPr166e//vWvev/99/XKK6/I4/FowIAB+vbbb+XxeCRJ8fHxdd4THx8fOObxeBQREaFWrVqdtqY+Pp9PpaWldTYriXE6NOeO3gqzSdmb9uutLQfMbgkAgJAWVEAZMWKEbrvtNqWlpWnIkCF6++23JUkLFy4M1NhstjrvMQzjlH0nO1tNVlaWXC5XYEtKSgqm7Ysi/ZJWeuD6yyRJM7K3yeM9ZnJHAACErvP6mnFMTIzS0tK0e/fuwLd5Tl4JKS4uDqyquN1u+f1+lZSUnLamPtOnT5fX6w1sBQUF59P2BfPrwSlK6+iSt7JKjy7nLrMAADTUeQUUn8+nzz//XB06dFBycrLcbrdycnICx/1+v3JzczVgwABJUnp6usLDw+vUFBUVadu2bYGa+jidTsXFxdXZrCjcHqY5d/SW0xGmT3Yf1l//sdfslgAACElBBZSpU6cqNzdXe/bs0T//+U/dfvvtKi0t1T333CObzaaMjAxlZmYqOztb27Zt07333qvo6GiNHz9ekuRyuTRhwgRNmTJFH3zwgTZt2qRf/OIXgVNGTcFl7VvoiRu7SZKy3t2pL4vLzvIOAABwMkcwxYWFhfr5z3+uw4cPq127drr66qv16aef6pJLLpEkTZs2TZWVlZo0aZJKSkrUr18/rVq1SrGxsYHfMWfOHDkcDo0bN06VlZUaPHiwFixYILvd3rgjM9G/9L9EH+ws1se7Diljab5W/OoaRTi4aS8AAOfKZhhGyF0oUVpaKpfLJa/X26ineyr81er+5PuSpB1PD1d0RFD5rY6Dpcc0fO7HOlJRpQevv0xTh3dtrDYBAAhJwXx+83/rL5D4uEhl3pomSXrhoy+V9813JncEAEDoIKBcQDemddCYPh1Va0gPL92scl+12S0BABASCCgX2Mybe6hjyyjt+65C//G/O8xuBwCAkEBAucDiIsP1/LhestmkpRsKtGr76e+YCwAAjiOgXARXd2mjf7+uiyRp+oqtOlTmM7kjAACsjYBykTwy7HJd4Y7Vt0f9evyNLQrBL08BAHDREFAuEqfDrrl39laEPUwf7CzWkvXWvF0/AABWQEC5iK5wx+nRE/dD+Y+3dmjv4aMmdwQAgDURUC6yCdcm6+ourVXhr9HDr+eruqbW7JYAALAcAspFFhZm0/PjeivW6dCmfUf0wkdfmd0SAACWQ0AxQceWUXr6lh6SpD9+sFubC46Y2xAAABZDQDHJLb07amTPDqqpNfTw6/mq9NeY3RIAAJZBQDGJzWbT729JVXycU18fOqqsdz83uyUAACyDgGKiltERem5sL0nSX//xjT76otjkjgAAsAYCismuS2mnewd0liQ9unyLSo76zW0IAAALIKBYwOMjrtCl7WJ0qMynJ7K3cpdZAECzR0CxgMhwu/54Zx85wmx6d5tHKzbuN7slAABMRUCxiNSOLj089HJJ0u/e3K6C7ypM7ggAAPMQUCxk4k+7KP2SVir3VWvK65tVU8upHgBA80RAsRCHPUxzxvVWTIRdn+39Tq988rXZLQEAYAoCisV0ahOtJ2/qLkl6ftUX2nGg1OSOAAC4+AgoFjSub5KGdo9XVY2hjKWbdKyKu8wCAJoXAooF2Ww2ZY1JU9sWEdp1sFzPvf+F2S0BAHBREVAsqm0Lp569rack6b/X7tG6Lw+b3BEAABcPAcXCBneL189/0kmSNGXZZnkrq0zuCACAi4OAYnG/GdlNndtEq8h7TL/7+zaz2wEA4KIgoFhcjNOh2Xf0lj3MppX5B/S/mw+Y3RIAABccASUEXNmplR64/jJJ0ozsrfJ4j5ncEQAAFxYBJURMvuEy9Ux0qfRYtaYu26xa7jILAGjCCCghItwepjl39FZkeJjWfnlYC/+x1+yWAAC4YAgoIeTSdi0048ZukqRn3t2p3QfLTO4IAIALg4ASYn5x9SUaeHk7+aprlbE0X/7qWrNbAgCg0RFQQozNZtMfbu+pltHh2n6gVHNX7zK7JQAAGt15BZSsrCzZbDZlZGQE9hmGoZkzZyohIUFRUVEaNGiQtm/fXud9Pp9PkydPVtu2bRUTE6PRo0ersLDwfFppVtrHRSrr1jRJ0ku5X2n93u9M7ggAgMbV4ICyfv16vfzyy+rZs2ed/bNmzdLs2bM1b948rV+/Xm63W0OHDlVZ2Q/XS2RkZCg7O1tLlizR2rVrVV5erlGjRqmmhofinasRaR1025WJqjWkR17PV9kx7jILAGg6GhRQysvLddddd+mVV15Rq1atAvsNw9DcuXM1Y8YMjRkzRqmpqVq4cKEqKiq0ePFiSZLX69X8+fP1/PPPa8iQIerTp48WLVqkrVu3avXq1Y0zqmbid6O7q2PLKBV8V6n/eGuH2e0AANBoGhRQHnjgAY0cOVJDhgyps3/Pnj3yeDwaNmxYYJ/T6dTAgQO1bt06SVJeXp6qqqrq1CQkJCg1NTVQczKfz6fS0tI6G6S4yHDNHtdLNpv0+oZCvb/dY3ZLAAA0iqADypIlS7Rx40ZlZWWdcszjOf4BGR8fX2d/fHx84JjH41FERESdlZeTa06WlZUll8sV2JKSkoJtu8nq16WN/v2nXSRJ01dsVXEZd5kFAIS+oAJKQUGBHnroIS1atEiRkZGnrbPZbHVeG4Zxyr6Tnalm+vTp8nq9ga2goCCYtpu8R4Zerm4d4vTdUb8eW75FhsFdZgEAoS2ogJKXl6fi4mKlp6fL4XDI4XAoNzdXf/rTn+RwOAIrJyevhBQXFweOud1u+f1+lZSUnLbmZE6nU3FxcXU2/MDpsGvuHb0V4QjTmi8OafFn+8xuCQCA8xJUQBk8eLC2bt2q/Pz8wNa3b1/dddddys/PV5cuXeR2u5WTkxN4j9/vV25urgYMGCBJSk9PV3h4eJ2aoqIibdu2LVCD4HV1x2ra8K6SpP9863N9fajc5I4AAGg4RzDFsbGxSk1NrbMvJiZGbdq0CezPyMhQZmamUlJSlJKSoszMTEVHR2v8+PGSJJfLpQkTJmjKlClq06aNWrduralTpyotLe2Ui24RnPuuSdaHO4u17qtv9fDrm/XG/f3lsHMvPgBA6AkqoJyLadOmqbKyUpMmTVJJSYn69eunVatWKTY2NlAzZ84cORwOjRs3TpWVlRo8eLAWLFggu93e2O00K2FhNj03tpeGz/1YmwuOaN6aL5Ux5HKz2wIAIGg2IwSvqCwtLZXL5ZLX623U61Eq/NXq/uT7kqQdTw9XdESj57eL4u/5+/XQknzZw2x641cD1DuppdktAQAQ1Oc36/9N0M29O+qmXgmqqTX08NJ8VfirzW4JAICgEFCaqP+8OVXuuEjtOXxUme98bnY7AAAEhYDSRLmiw/Xc2F6SpEWf7tOaL4pN7ggAgHNHQGnCrk1pq3+9prMkadryLfruqN/chgAAOEcElCbusZ9doZT2LXSozKfpK7jLLAAgNBBQmrjIcLvm3NFb4Xab3t9+UMvzCs1uCQCAsyKgNAOpHV2B+6E89b87VPBdhckdAQBwZgSUZuL+gZeq7yWtVO6r1iOv56umllM9AADrIqA0E/Ywm+bc0VsxEXat31uilz/+2uyWAAA4LQJKM5LUOlq/G91DkjQ75wtt2+81uSMAAOpHQGlmxqYnalj3eFXVHL/L7LGqGrNbAgDgFASUZsZmsylrTJratnBqd3G5Zr33hdktAQBwCgJKM9SmhVOzbk+TJP3P/+3R/3152OSOAACoi4DSTN1wRbzu6tdJkjTl9c3yVlSZ3BEAAD8goDRjM0Z2U3LbGHlKj+m3f99mdjsAAAQQUJqx6AiHZo/rJXuYTW9uPqC/5+83uyUAACQRUJq9Pp1a6cHrL5Mk/XblNh04UmlyRwAAEFAg6cEbLlOvpJYqPVatqcs2q5a7zAIATEZAgcLtYZozrpeiwu1a99W3enXdXrNbAgA0cwQUSJK6tGuhGSO7SZKefW+ndh0sM7kjAEBzRkBBwF39Oun6ru3kr65VxpJ8+atrzW4JANBMEVAQYLPZ9OztPdUqOlw7iko1Z/Uus1sCADRTBBTU0T42UlljekqSXsr9Sp/t+c7kjgAAzREBBaf4WapbY9MTZRjSI6/nq+wYd5kFAFxcBBTU68mbuiuxVZQKSyr11P/uMLsdAEAzQ0BBvWIjwzXnjt6y2aTleYV6b1uR2S0BAJoRAgpO66rOrXX/wEslSdNXbFVx6TGTOwIANBcEFJzRw0MuV/cOcSqpqNK0N7bIMLjLLADgwiOg4IwiHGGae2dvRTjC9NEXh7Ton/vMbgkA0AwQUHBWl8fH6rGfXSFJ+v3bO/T1oXKTOwIANHUEFJyTfx3QWddc1kbHqmr18NJ8VdVwl1kAwIVDQME5CQuz6bmxvRQX6dDmQq/mffil2S0BAJqwoALKiy++qJ49eyouLk5xcXHq37+/3n333cBxwzA0c+ZMJSQkKCoqSoMGDdL27dvr/A6fz6fJkyerbdu2iomJ0ejRo1VYWNg4o8EF1cEVpf+8NU2SNG/Nl9q0r8TkjgAATVVQASUxMVHPPPOMNmzYoA0bNuiGG27QzTffHAghs2bN0uzZszVv3jytX79ebrdbQ4cOVVnZD0/GzcjIUHZ2tpYsWaK1a9eqvLxco0aNUk1NTeOODBfE6F4JGt0rQTW1hh5emq8Kf7XZLQEAmiCbcZ7fG23durX+8Ic/6L777lNCQoIyMjL02GOPSTq+WhIfH69nn31WEydOlNfrVbt27fTaa6/pjjvukCQdOHBASUlJeueddzR8+PBz+jNLS0vlcrnk9XoVFxd3Pu3XUeGvVvcn35ck7Xh6uKIjHI32u5sSb0WVfvbHj1XkPabx/Top88SqCgAAZxLM53eDr0GpqanRkiVLdPToUfXv31979uyRx+PRsGHDAjVOp1MDBw7UunXrJEl5eXmqqqqqU5OQkKDU1NRATX18Pp9KS0vrbDCPKzpcz43tJUla/M99+nDnQZM7AgA0NUEHlK1bt6pFixZyOp26//77lZ2dre7du8vj8UiS4uPj69THx8cHjnk8HkVERKhVq1anralPVlaWXC5XYEtKSgq2bTSyay5rqwnXJkuSpi3fqm/LfSZ3BABoSoIOKF27dlV+fr4+/fRT/epXv9I999yjHTt+eJiczWarU28Yxin7Tna2munTp8vr9Qa2goKCYNvGBfDo8K66PL6FDpf7NH3FVu4yCwBoNEEHlIiICF122WXq27evsrKy1KtXL/3xj3+U2+2WpFNWQoqLiwOrKm63W36/XyUlJaetqY/T6Qx8c+j7DeaLDLdrzh29FW63adWOg1qWx7exAACN47zvg2IYhnw+n5KTk+V2u5WTkxM45vf7lZubqwEDBkiS0tPTFR4eXqemqKhI27ZtC9QgtPRIcOmRoV0lSU+9uV37vq0wuSMAQFMQ1NdUnnjiCY0YMUJJSUkqKyvTkiVL9NFHH+m9996TzWZTRkaGMjMzlZKSopSUFGVmZio6Olrjx4+XJLlcLk2YMEFTpkxRmzZt1Lp1a02dOlVpaWkaMmTIBRkgLrx//2kXrdlZrM/2fqdHXs/X0on9ZQ8782k9AADOJKiAcvDgQd19990qKiqSy+VSz5499d5772no0KGSpGnTpqmyslKTJk1SSUmJ+vXrp1WrVik2NjbwO+bMmSOHw6Fx48apsrJSgwcP1oIFC2S32xt3ZA0QHeHQ3mdGmt1GyLGH2fT8uF4a8cdPtOGbEr2U+5UeuP4ys9sCAISw874Pihku1H1QcH6W5xVq6rLNcoTZtPKBa5Ta0WV2SwAAC7ko90EBTnbblR31sx5uVdcayliar2NV3B0YANAwBBQ0GpvNpswxaWoX69SXxeV69r2dZrcEAAhRBBQ0qtYxEZp1e09J0qv/t1ef7D5kckcAgFBEQEGju75re9199SWSpKnLNutIhd/kjgAAoYaAggviiRu7qUvbGB0s9ek3K7dxl1kAQFAIKLggoiKO32XWHmbTW1uK9ObmA2a3BAAIIQQUXDC9klrq1zekSJJ+s3Kb9h+pNLkjAECoIKDggnrg+kvVO6mlyo5Va+rrm1Vby6keAMDZEVBwQTnsYZpzR29Fhdv1j6+/1f/83x6zWwIAhAACCi645LYx+s2obpKkWe99oS88ZSZ3BACwOgIKLorxP+mkG65oL39NrR5askm+au4yCwA4PQIKLgqbzaZnbktT65gI7fSUaXbOrgvy51T4q9X58bfV+fG3VeGvviB/BgDgwiOg4KJpHxuprDFpkqSXP/5an379rckdAQCsioCCi2p4D7fG9U2UYUhTXt+s0mNVZrcEALAgAgouuidv6qGk1lHaf6RST725w+x2AAAWREDBRdfC6dCccb0VZpPe2Fiod7cWmd0SAMBiCCgwRd/OrfWrQZdKkqZnb1Vx6TGTOwIAWAkBBaZ5aPDlSu0YpyMVVXp0+RYeKAgACCCgwDQRjjDNGddbTkeYcncd0qJPvzG7JQCARRBQYKqU+Fg9PuIKSdLv3/lcXx0qN7kjAIAVEFBgunv6d9Z1KW11rKpWDy/NV1VNrdktAQBMRkCB6cLCbPrD7b3kigrXlkKv/vzBbrNbAgCYjIACS3C7IvX7W1MlSfPWfKmN+0pM7ggAYCYCCixjVM8E3dI7QbWG9PDSfB318SwdAGiuCCiwlKduTlWCK1LffFuh/3z7c7PbAQCYhIACS3FFheu5cb1ks0l/+2yfVu84aHZLAAATEFBgOQMubatfXpssSXp8xRYdLveZ3BEA4GIjoMCSpgzrqivcsTpc7tfjb2zlLrMA0MwQUGBJkeF2zbmjtyLsYVr9+UG9vqHA7JYAABcRAQWW1a1DnKYMu1yS9NT/7tA33x41uSMAwMVCQIGl/fK6LvpJcmtV+Gv08NJ8VXOXWQBoFggosDR7mE2zx/VSC6dDG/cd0Uu5X5ndEgDgIggqoGRlZemqq65SbGys2rdvr1tuuUVffPFFnRrDMDRz5kwlJCQoKipKgwYN0vbt2+vU+Hw+TZ48WW3btlVMTIxGjx6twsLC8x8NmqTEVtF6anQPSdLc1bu1tdBrckcAgAstqICSm5urBx54QJ9++qlycnJUXV2tYcOG6ejRH64NmDVrlmbPnq158+Zp/fr1crvdGjp0qMrKygI1GRkZys7O1pIlS7R27VqVl5dr1KhRqqmpabyRoUkZc2VH3ZjmVnWtoYylm3Ssin9XAKApsxnn8f3NQ4cOqX379srNzdVPf/pTGYahhIQEZWRk6LHHHpN0fLUkPj5ezz77rCZOnCiv16t27drptdde0x133CFJOnDggJKSkvTOO+9o+PDhZ/1zS0tL5XK55PV6FRcX19D2EWJKjvo1fO7HKi7z6d4BnTXzxKrKj1X4q9X9yfclSTueHq7oCMfFbhMAcBrBfH6f1zUoXu/xpfbWrVtLkvbs2SOPx6Nhw4YFapxOpwYOHKh169ZJkvLy8lRVVVWnJiEhQampqYEaoD6tYiL0h7G9JEkL1u3Vx7sOmdwRAOBCaXBAMQxDjzzyiK699lqlph5/Cq3H45EkxcfH16mNj48PHPN4PIqIiFCrVq1OW3Myn8+n0tLSOhuap4GXt9O/9L9EkvTo8s06UuE3uSMAwIXQ4IDy4IMPasuWLfrb3/52yjGbzVbntWEYp+w72ZlqsrKy5HK5AltSUlJD20YTMH1EN3VpF6ODpT7NyN7GXWYBoAlqUECZPHmy3nzzTa1Zs0aJiYmB/W63W5JOWQkpLi4OrKq43W75/X6VlJSctuZk06dPl9frDWwFBdxVtDmLirBr7h295Qiz6e2tRVqZv9/slgAAjSyogGIYhh588EGtWLFCH374oZKTk+scT05OltvtVk5OTmCf3+9Xbm6uBgwYIElKT09XeHh4nZqioiJt27YtUHMyp9OpuLi4Ohuat56JLfXQ4BRJ0pMrt2v/kUqTOwIANKagAsoDDzygRYsWafHixYqNjZXH45HH41Fl5fEPB5vNpoyMDGVmZio7O1vbtm3Tvffeq+joaI0fP16S5HK5NGHCBE2ZMkUffPCBNm3apF/84hdKS0vTkCFDGn+EaLJ+NehS9enUUmW+ak15PV+1tZzqAYCmIqjvYL744ouSpEGDBtXZ/+qrr+ree++VJE2bNk2VlZWaNGmSSkpK1K9fP61atUqxsbGB+jlz5sjhcGjcuHGqrKzU4MGDtWDBAtnt9vMbDZoVhz1Mc8b11o1/+kSffv2d5q/do7uu7mR2WwCARnBe90ExC/dBwY/97bN9mr5iqyLsYVo68Wrd+sLxr6tzHxQAsJaLdh8UwAruvCpJQ7q1l7+mVo+9scXsdgAAjYCAgpBns9mUNaan2sREaNfBcrPbAQA0AgIKmoR2sU49c1vPOvt4Xg8AhC4CCpqMod3jdXt6x8Drn876SNNXbNXGfSXczA0AQgwXyaJJKTnqU5//WH3K/svat9Dt6Ym6tU9HxcdFmtAZAICLZNFsOcN/+Kr6/9zbV7f26ajI8DB9WVyuZ97dqf5ZH+hfX/1M72wtkq+aU0AAYFV8BxNN1tVd2uiGK+L19M099PaWIi3LK1TeNyVa88UhrfnikFpGh+uW3h11e3qieiTEnfV5UQCAi4eAgiYvNjJcd/6kk+78SSd9dahcy/MKtWJjoQ6W+rRg3V4tWLdXV7hjNbZvkm7pnaA2LZxmtwwAzR7XoKBJqfBXq/uT70s6843aamoNfbL7kJblFSpn+0H5a2olSY4wmwZ3a6/b05M0qGs7hds5CwoAjSWYz29WUNAs2cNsGtS1vQZ1ba8jFX69ufmAlucVakuhV+9vP6j3tx9U2xZO3donQWP7Juny+Niz/1IAQKMhoKDZaxkdoX/p31n/0r+zdnpKtXxDoVbm79fhcp9e+WSPXvlkj3olunR7eqJG9+ooV3S42S0DQJPHKR40Ked6iudsqmpqtWZnsZbnFerDncWqPvGk5AhHmIZ1j9fYvkm69rK2sodxYS0AnCtO8QDnKdwepmE93BrWw63D5T6t3LRfy/MKtdNTpre2FOmtLUVyx0VqzJXHvwXUpV0Ls1sGgCaFgAKcRdsWTv3yui6acG2ytu0v1fK8Aq3MPyBP6TG98NFXeuGjr9T3klYa2zdRI3smqIWT/6wA4HzxNylwjmw2m9ISXUpLdOmJkd20ekexlucVKHfXIW34pkQbvinRzDd3aESqW7f3TdTVyW0UxikgAGgQAgrQAE6HXSN7dtDInh10sPSYVmzcr2V5Bfr60FGt2LRfKzbtV1LrKN12ZaJuuzJRSa2jzW4ZAEIKF8miSWmsi2QbwjAMbdx3RMvzCvXW5gMq81UHjvXv0kZj+yZqRGoHRUXYz/BbAKDpCubzm4CCJsXMgPJjlf4avb/do2V5BVr31bf6/r+yFk6HRvXsoLF9E3Vlp1bcXh9As8K3eACTRUXYdUufjrqlT0cVllRoxcbj3wLa912Flqwv0JL1BerSNka3pR8/BeR28YRlAPgxVlDQpFhlBaU+tbWGPtv7nZZtKNQ7W4tUWXX8acphNum6lHYa2zdRQ7rFKzKcU0AAmiZWUAALCguz6eoubXR1lzZ66uYeemdrkZZvKNRne79T7q5Dyt11SK6ocI3ulaCxfROV1tHFKSAAzRYBBTBBC6dD4/omaVzfJO09fFTL8wr1xsZCFXmP6bVPv9Frn36jrvGxGts3UTf37qh2sTxhGUDzwikeNClWPsVzNjW1htZ9dVjLNhTqve0e+at/eMLyoK7tNbZvom64oj1PWAYQsjjFA4Qge5hN16W003Up7eStrNL/nnjCcn7BEa3+/KBWf35QbWIidEuf47fX79aBcA6g6SKgABbkigrXL66+RL+4+hLtPlim5XmFWrFpvw6V+TR/7R7NX7tHqR3jNDY9STf3TlDL6AizWwaARsUpHjQpoXyK52yqa2qVu+uQlucVavXnB1VVc+IJy/YwDeneXmPTk3RdSls5OAUEwKI4xQM0QQ57mAZ3i9fgbvH67qhff8/fr2UbCrWjqFTvbPXona0etY91asyVibo9PVGXtecJywBCFwEFCEGtYyL0r9ck61+vSdb2A14tzyvUyk37VVzm00u5X+ml3K/Up1NLjU1P0qheHRQXGW52ywAQFAIKEOJ6JLjUI8Gl6SO66cOdB7VsQ6E+2nVIm/Yd0aZ9R/T0W9v1sx5uje2bpP5deMIygNBAQAGaiAhHmH6W2kE/S+2g4rJjWrnp+Cmg3cXlWpl/QCvzD6hjyyjddmVH3Z6epE5teMIyAOviIlk0KU35ItmGMAxDmwu9WrahQG9uPqCyYz88YblfcmuN7ZukEaluxTib9/9OAC4OLpIFIEmy2WzqndRSvZNa6rejumvVjoNatqFAa788rH/u+U7/3POdfvf3bboxrYPG9k3SVZ3NfcIyARPA9/ivH2gmIsPtGt0rQaN7JejAkUqt2Fio5XmF2vtthZblFWpZXqE6t4nW7emJGnNlohJaRpndMoBmLOgbJnz88ce66aablJCQIJvNppUrV9Y5bhiGZs6cqYSEBEVFRWnQoEHavn17nRqfz6fJkyerbdu2iomJ0ejRo1VYWHheAwFw7hJaRunBG1K0ZuogLbu/v8b1TVRMhF17v63Qc6t26ZpnP9Td8/+pv+fv17ETT10GgIsp6IBy9OhR9erVS/Pmzav3+KxZszR79mzNmzdP69evl9vt1tChQ1VWVhaoycjIUHZ2tpYsWaK1a9eqvLxco0aNUk0NfxECF5PNZtNVnVtr1u29tP43Q/T82F66uktrGYb0ye7DemhJvq76/WrNyN6qTftKFIKXrAEIUUGf4hkxYoRGjBhR7zHDMDR37lzNmDFDY8aMkSQtXLhQ8fHxWrx4sSZOnCiv16v58+frtdde05AhQyRJixYtUlJSklavXq3hw4efx3AANFR0hEO3pSfqtvRE7fu2Qss3FuqNvELtP1Kp//fPffp//9ynlPYtdHt6om69sqPax0aa3TKAJqxR74m9Z88eeTweDRs2LLDP6XRq4MCBWrdunSQpLy9PVVVVdWoSEhKUmpoaqDmZz+dTaWlpnQ3AhdOpTbQeGXq5Ppl2vRb/sp9u7dNRkeFh2l1crqx3d6p/1oeasGC93t1aFHjqMgA0pka9SNbj8UiS4uPj6+yPj4/XN998E6iJiIhQq1atTqn5/v0ny8rK0lNPPdWYrQI4B2FhNg24rK0GXNZWT93cQ29vKdKyDQXauO+IPthZrA92FqtVdLhu7t1RY/smqkeCy+yWATQRF+RbPCd/TdEwjLN+dfFMNdOnT9cjjzwSeF1aWqqkpKTzbxTAOYuLDNfPf9JJP/9JJ311qPz4E5Y3FupgqU8L1u3VgnV71b1DnG5PT9QtfTqqdQxPWAbQcI0aUNxut6TjqyQdOnQI7C8uLg6sqrjdbvn9fpWUlNRZRSkuLtaAAQPq/b1Op1NOp7MxW0UTFR3h0N5nRprdRpN3absWeuxnV2jK0Mv1yZeHtXxDoXJ2HNSOolI9/dYOZb37uQZfEa+xfRM18PJ2PGEZQNAa9W+N5ORkud1u5eTkBPb5/X7l5uYGwkd6errCw8Pr1BQVFWnbtm2nDSgArMlhD9P1XdvrL3ddqc9mDNbTN/dQWkeXqmoMvbfdowkLN+jqrA+V+c7n2n2w7Oy/EABOCHoFpby8XF9++WXg9Z49e5Sfn6/WrVurU6dOysjIUGZmplJSUpSSkqLMzExFR0dr/PjxkiSXy6UJEyZoypQpatOmjVq3bq2pU6cqLS0t8K0eAKGnZXSE/qV/Z/1L/87a6SnVsg3Hn7B8uNynlz/+Wi9//LV6JbXU2PRE3dQrQa4onrAM4PSCfhbPRx99pOuvv/6U/ffcc48WLFggwzD01FNP6b/+679UUlKifv366S9/+YtSU1MDtceOHdOjjz6qxYsXq7KyUoMHD9YLL7xwzteV8CweIDRU1dRqzc5iLcsr1JqdxaquPf7XTYQjTMN7uDU2PVHXXNZW9hNPWOZW96GPOcSZBPP5zcMCAVwUh8t9Wrlpv5bnFWqn54fTPR1ckbrtyuP3X4mPc/LhFuIIKDgTHhYIwHLatnDql9d10YRrk7Vtf6mW5RXo7/kHVOQ9pnlrvtS8NV/qyk4tzW4TaPasEjIJKAAuKpvNprREl9ISXXrixm764PNiLcsr0Me7DmnjviOBuv5ZHyouyqEWznDFRjoUF+lQC6dDsZHhahHpUGzk8Z9jnT/83OLEz3GR4Ypx2vn2EBDCCCgATBMZbtfInh00smcHebzHtHT9Ps1ZvVuS5K2skreySlJlg39/dIRdsT8KNrHfB5sToadF5I/2/6jm+wAUFxkupyPsrPdxAtD4CCgALMHtitS//bRLIKC8+eA1qq41VHasWmXHqlR+rDrwc5nvh5/LAz//sM934vb7Ff4aVfhrdFC+BvflCLOdskJTJ+xE/rDKExtZ93hgxcfpCFwIDODcEFAAWNJl7Vs0+Ny3v7r2RHCpqhNcyo5V191/ItyU/6iu3Fet0hPBxzCk6lpDJRVVKqmoOq/xxETYT1qh+fHKzckB6NSQExvpYDUHzQoBBUCTE+EIU2tHxHndbr+21lBFVc3ZQ86x+lZzqk4EnerAwxSP+mt01F8jz3k86zTcbjspuNRzHc7prs/5frUnwqEwVnMQAggoAFCPsDCbWjiPB4EO5/EMRF91zY9OT1WrzPej1ZqTVnLqC0DlJ45LUlWNoe+O+vXdUf95je3HKzWnXJ9zItic7focp8N+Xj0AZ0NAAYALyOmwy9nCrjYtGv48sdpaQ+X+6tNfh3NO1+dUqarm+G2vyn3HQ1CRt+HjinCE1Xt6Kjrih+Dypw92KzLcLkeYTfawsBP/tMlhP/HPk/fXOV7P/rCwet7/o/1hNtntdfeH2U59gC1CAwEFACwuLMymuMhwxUWe3+MBjlXV1FmhKT92/DTU6U5P1Xd9TvmJ1Rx/da2+rfbr2zOs5ryU+/V59dtYTg1AYScFnCAC0BkD1In9dpvCbCftt1+gQFbneP3jCtWARkABgGYiMtyuyHC72sU2fDWnptYIBJUfn4YqPfFzyVG/ns/ZJUm6q18n2WzH31NdYxz/Z+33/6yt+7rmNPt/XF9zmv21RmB1qD7VJ+ob/l2u0GY/XZA6zYqWVS5RIqAAAM6ZPcwmV1T4iYc9Rp1yvMJfHQgoM0Z2u6h3Ia09UwCqNU4EnPoCUG09AerE/lMC1En76xyvZ//3f14jBLTaWtU9XnPS8dNktJoTNed35dLFR0ABADQJYWE2RQT+73/zu4i3ttZQjXGOgewMoemor1r3L9po9nAIKAAANAVhYTaFyabw88xmFf7qxmnoPPGgCgAAYDkEFAAAYDkEFAAAYDkEFAAAYDlcJAvAMqIjHNr7zEiz2wBgAaygAAAAy2EFBQDQaFgFQ2MhoAAAgACrhExO8QAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMshoAAAAMsxNaC88MILSk5OVmRkpNLT0/XJJ5+Y2Q4AALAI0wLK0qVLlZGRoRkzZmjTpk267rrrNGLECO3bt8+slgAAgEXYDMMwzPiD+/XrpyuvvFIvvvhiYF+3bt10yy23KCsr64zvLS0tlcvlktfrVVxc3IVuFQAANIJgPr9NWUHx+/3Ky8vTsGHD6uwfNmyY1q1bd0q9z+dTaWlpnQ0AADRdDjP+0MOHD6umpkbx8fF19sfHx8vj8ZxSn5WVpaeeeuqU/QQVAABCx/ef2+dy8saUgPI9m81W57VhGKfsk6Tp06frkUceCbzev3+/unfvrqSkpAveIwAAaFxlZWVyuVxnrDEloLRt21Z2u/2U1ZLi4uJTVlUkyel0yul0Bl63aNFCBQUFio2NrTfQnI/S0lIlJSWpoKCgSV7f0tTHJzX9MTK+0NfUx8j4Qt+FGqNhGCorK1NCQsJZa00JKBEREUpPT1dOTo5uvfXWwP6cnBzdfPPNZ31/WFiYEhMTL2SLiouLa7L/4klNf3xS0x8j4wt9TX2MjC/0XYgxnm3l5HumneJ55JFHdPfdd6tv377q37+/Xn75Ze3bt0/333+/WS0BAACLMC2g3HHHHfr222/19NNPq6ioSKmpqXrnnXd0ySWXmNUSAACwCFMvkp00aZImTZpkZguncDqd+t3vflfnmpempKmPT2r6Y2R8oa+pj5HxhT4rjNG0G7UBAACcDg8LBAAAlkNAAQAAlkNAAQAAlkNAAQAAltMsA8oLL7yg5ORkRUZGKj09XZ988skZ63Nzc5Wenq7IyEh16dJFL7300kXqtGGCGd9HH30km812yrZz586L2PG5+/jjj3XTTTcpISFBNptNK1euPOt7Qmn+gh1fqM1fVlaWrrrqKsXGxqp9+/a65ZZb9MUXX5z1faE0hw0ZYyjN44svvqiePXsGbuDVv39/vfvuu2d8TyjNX7DjC6W5q09WVpZsNpsyMjLOWGfGHDa7gLJ06VJlZGRoxowZ2rRpk6677jqNGDFC+/btq7d+z549uvHGG3Xddddp06ZNeuKJJ/TrX/9ab7zxxkXu/NwEO77vffHFFyoqKgpsKSkpF6nj4Bw9elS9evXSvHnzzqk+1OYv2PF9L1TmLzc3Vw888IA+/fRT5eTkqLq6WsOGDdPRo0dP+55Qm8OGjPF7oTCPiYmJeuaZZ7RhwwZt2LBBN9xwg26++WZt37693vpQm79gx/e9UJi7k61fv14vv/yyevbsecY60+bQaGZ+8pOfGPfff3+dfVdccYXx+OOP11s/bdo044orrqizb+LEicbVV199wXo8H8GOb82aNYYko6Sk5CJ017gkGdnZ2WesCbX5+7FzGV8oz59hGEZxcbEhycjNzT1tTSjPoWGc2xhDfR5btWpl/Pd//3e9x0J9/gzjzOML1bkrKyszUlJSjJycHGPgwIHGQw89dNpas+awWa2g+P1+5eXladiwYXX2Dxs2TOvWrav3Pf/4xz9OqR8+fLg2bNigqqqqC9ZrQzRkfN/r06ePOnTooMGDB2vNmjUXss2LKpTm73yE6vx5vV5JUuvWrU9bE+pzeC5j/F6ozWNNTY2WLFmio0ePqn///vXWhPL8ncv4vhdqc/fAAw9o5MiRGjJkyFlrzZrDZhVQDh8+rJqamlOemBwfH3/Kk5W/5/F46q2vrq7W4cOHL1ivDdGQ8XXo0EEvv/yy3njjDa1YsUJdu3bV4MGD9fHHH1+Mli+4UJq/hgjl+TMMQ4888oiuvfZapaamnrYulOfwXMcYavO4detWtWjRQk6nU/fff7+ys7PVvXv3emtDcf6CGV+ozZ0kLVmyRBs3blRWVtY51Zs1h6be6t4sNputzmvDME7Zd7b6+vZbRTDj69q1q7p27Rp43b9/fxUUFOi5557TT3/60wva58USavMXjFCevwcffFBbtmzR2rVrz1obqnN4rmMMtXns2rWr8vPzdeTIEb3xxhu65557lJube9oP8VCbv2DGF2pzV1BQoIceekirVq1SZGTkOb/PjDlsVisobdu2ld1uP2U1obi4+JR0+D23211vvcPhUJs2bS5Yrw3RkPHV5+qrr9bu3bsbuz1ThNL8NZZQmL/JkyfrzTff1Jo1a5SYmHjG2lCdw2DGWB8rz2NERIQuu+wy9e3bV1lZWerVq5f++Mc/1lsbivMXzPjqY+W5y8vLU3FxsdLT0+VwOORwOJSbm6s//elPcjgcqqmpOeU9Zs1hswooERERSk9PV05OTp39OTk5GjBgQL3v6d+//yn1q1atUt++fRUeHn7Bem2IhoyvPps2bVKHDh0auz1ThNL8NRYrz59hGHrwwQe1YsUKffjhh0pOTj7re0JtDhsyxvpYeR5PZhiGfD5fvcdCbf7qc6bx1cfKczd48GBt3bpV+fn5ga1v37666667lJ+fL7vdfsp7TJvDC3oJrgUtWbLECA8PN+bPn2/s2LHDyMjIMGJiYoy9e/cahmEYjz/+uHH33XcH6r/++msjOjraePjhh40dO3YY8+fPN8LDw43ly5ebNYQzCnZ8c+bMMbKzs41du3YZ27ZtMx5//HFDkvHGG2+YNYQzKisrMzZt2mRs2rTJkGTMnj3b2LRpk/HNN98YhhH68xfs+EJt/n71q18ZLpfL+Oijj4yioqLAVlFREagJ9TlsyBhDaR6nT59ufPzxx8aePXuMLVu2GE888YQRFhZmrFq1yjCM0J+/YMcXSnN3Oid/i8cqc9jsAophGMZf/vIX45JLLjEiIiKMK6+8ss7X/+655x5j4MCBdeo/+ugjo0+fPkZERITRuXNn48UXX7zIHQcnmPE9++yzxqWXXmpERkYarVq1Mq699lrj7bffNqHrc/P9V/pO3u655x7DMEJ//oIdX6jNX31jk2S8+uqrgZpQn8OGjDGU5vG+++4L/P3Srl07Y/DgwYEPb8MI/fkLdnyhNHenc3JAscoc2gzjxJUuAAAAFtGsrkEBAAChgYACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAsh4ACAAAs5/8DXU8QQ1ZeL9cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ordination, dm, cv = rpca.joint_rpca(tables=[wol2_bt, food_data_bt], n_test_samples = 200, sample_metadata=foods, train_test_column='train_test', n_components=3)\n",
    "corr_table = rpca.feature_correlation_table(ordination)\n",
    "plt.errorbar(x=cv.index, y=cv['mean_CV'], yerr=cv['std_CV']/2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "9bd01f75-4b7f-4f3f-8eb1-9f67cc24bffc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "615 x 1177 <class 'biom.table.Table'> with 219444 nonzero entries (30% dense)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "food_data_bt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "69428cc1-a918-44e1-9c3e-79c878b3743f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'emperor_biplot.qzv'"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "biplot = q2.Artifact.import_data(\"PCoAResults % Properties('biplot')\", ordination)\n",
    "\n",
    "emp_biplot = emperor.visualizers.biplot(biplot, foods_q2, number_of_features = 10).visualization\n",
    "emp_biplot.save('emperor_biplot.qzv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "deae866c-57be-4f0a-bb25-63d14b92c663",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_table.loc[corr_table.index.intersection(ultraprocessed)][corr_table.index.intersection(wol2_bt.ids(axis='observation'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "392775e4-7ac3-4781-bb2c-ffad2cf6be2a",
   "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>featureid</th>\n",
       "      <th>G000006605</th>\n",
       "      <th>G000006785</th>\n",
       "      <th>G000006845</th>\n",
       "      <th>G000006865</th>\n",
       "      <th>G000006925</th>\n",
       "      <th>G000007265</th>\n",
       "      <th>G000007325</th>\n",
       "      <th>G000007525</th>\n",
       "      <th>G000007745</th>\n",
       "      <th>G000007785</th>\n",
       "      <th>...</th>\n",
       "      <th>G902794445</th>\n",
       "      <th>G902798515</th>\n",
       "      <th>G902798535</th>\n",
       "      <th>G902798655</th>\n",
       "      <th>G902798845</th>\n",
       "      <th>G902799635</th>\n",
       "      <th>G902800075</th>\n",
       "      <th>G902801795</th>\n",
       "      <th>G902801825</th>\n",
       "      <th>G902801975</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Processed_Calories_Nova_processed_culinary_ingredients</th>\n",
       "      <td>0.999549</td>\n",
       "      <td>0.999796</td>\n",
       "      <td>0.999636</td>\n",
       "      <td>-0.991809</td>\n",
       "      <td>-0.991203</td>\n",
       "      <td>0.999559</td>\n",
       "      <td>0.999693</td>\n",
       "      <td>-0.999956</td>\n",
       "      <td>0.999745</td>\n",
       "      <td>-0.949400</td>\n",
       "      <td>...</td>\n",
       "      <td>0.999749</td>\n",
       "      <td>0.998335</td>\n",
       "      <td>0.999840</td>\n",
       "      <td>0.999658</td>\n",
       "      <td>0.999605</td>\n",
       "      <td>0.997464</td>\n",
       "      <td>0.999948</td>\n",
       "      <td>0.996862</td>\n",
       "      <td>0.996540</td>\n",
       "      <td>0.998806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Calories_Nova_processed_foods</th>\n",
       "      <td>-0.628812</td>\n",
       "      <td>-0.658344</td>\n",
       "      <td>-0.664031</td>\n",
       "      <td>0.718708</td>\n",
       "      <td>0.732671</td>\n",
       "      <td>-0.635961</td>\n",
       "      <td>-0.651258</td>\n",
       "      <td>0.652493</td>\n",
       "      <td>-0.662388</td>\n",
       "      <td>0.855622</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.655911</td>\n",
       "      <td>-0.656404</td>\n",
       "      <td>-0.651154</td>\n",
       "      <td>-0.657972</td>\n",
       "      <td>-0.641590</td>\n",
       "      <td>-0.686127</td>\n",
       "      <td>-0.657242</td>\n",
       "      <td>-0.666889</td>\n",
       "      <td>-0.675265</td>\n",
       "      <td>-0.670311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Calories_Nova_ratio_any_processed_unprocessed</th>\n",
       "      <td>0.999649</td>\n",
       "      <td>0.999606</td>\n",
       "      <td>0.999437</td>\n",
       "      <td>-0.987259</td>\n",
       "      <td>-0.987037</td>\n",
       "      <td>0.999939</td>\n",
       "      <td>0.999869</td>\n",
       "      <td>-0.999598</td>\n",
       "      <td>0.999446</td>\n",
       "      <td>-0.945420</td>\n",
       "      <td>...</td>\n",
       "      <td>0.999729</td>\n",
       "      <td>0.999458</td>\n",
       "      <td>0.999795</td>\n",
       "      <td>0.999698</td>\n",
       "      <td>0.999976</td>\n",
       "      <td>0.997899</td>\n",
       "      <td>0.999369</td>\n",
       "      <td>0.998393</td>\n",
       "      <td>0.997919</td>\n",
       "      <td>0.999167</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Calories_Nova_ratio_ultra_processed_unprocessed</th>\n",
       "      <td>0.998660</td>\n",
       "      <td>0.999770</td>\n",
       "      <td>0.999835</td>\n",
       "      <td>-0.989393</td>\n",
       "      <td>-0.989820</td>\n",
       "      <td>0.999287</td>\n",
       "      <td>0.999805</td>\n",
       "      <td>-0.999491</td>\n",
       "      <td>0.999759</td>\n",
       "      <td>-0.954368</td>\n",
       "      <td>...</td>\n",
       "      <td>0.999821</td>\n",
       "      <td>0.999795</td>\n",
       "      <td>0.999687</td>\n",
       "      <td>0.999887</td>\n",
       "      <td>0.999551</td>\n",
       "      <td>0.999361</td>\n",
       "      <td>0.999421</td>\n",
       "      <td>0.999263</td>\n",
       "      <td>0.999114</td>\n",
       "      <td>0.999935</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Calories_Nova_ultra_processed_foods</th>\n",
       "      <td>-0.999571</td>\n",
       "      <td>-0.998286</td>\n",
       "      <td>-0.997798</td>\n",
       "      <td>0.988823</td>\n",
       "      <td>0.987195</td>\n",
       "      <td>-0.999025</td>\n",
       "      <td>-0.998394</td>\n",
       "      <td>0.998875</td>\n",
       "      <td>-0.998057</td>\n",
       "      <td>0.936735</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.998288</td>\n",
       "      <td>-0.996197</td>\n",
       "      <td>-0.998667</td>\n",
       "      <td>-0.998044</td>\n",
       "      <td>-0.998756</td>\n",
       "      <td>-0.993949</td>\n",
       "      <td>-0.998709</td>\n",
       "      <td>-0.993857</td>\n",
       "      <td>-0.993122</td>\n",
       "      <td>-0.996286</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Grams_Subgroup_table_salt</th>\n",
       "      <td>0.999803</td>\n",
       "      <td>0.999793</td>\n",
       "      <td>0.999596</td>\n",
       "      <td>-0.990149</td>\n",
       "      <td>-0.989571</td>\n",
       "      <td>0.999852</td>\n",
       "      <td>0.999838</td>\n",
       "      <td>-0.999939</td>\n",
       "      <td>0.999686</td>\n",
       "      <td>-0.946938</td>\n",
       "      <td>...</td>\n",
       "      <td>0.999810</td>\n",
       "      <td>0.998721</td>\n",
       "      <td>0.999921</td>\n",
       "      <td>0.999724</td>\n",
       "      <td>0.999863</td>\n",
       "      <td>0.997471</td>\n",
       "      <td>0.999839</td>\n",
       "      <td>0.997299</td>\n",
       "      <td>0.996883</td>\n",
       "      <td>0.998890</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Grams_Subgroup_table_sugar</th>\n",
       "      <td>0.995554</td>\n",
       "      <td>0.998346</td>\n",
       "      <td>0.998689</td>\n",
       "      <td>-0.995988</td>\n",
       "      <td>-0.996637</td>\n",
       "      <td>0.996189</td>\n",
       "      <td>0.997641</td>\n",
       "      <td>-0.997986</td>\n",
       "      <td>0.998647</td>\n",
       "      <td>-0.967969</td>\n",
       "      <td>...</td>\n",
       "      <td>0.998094</td>\n",
       "      <td>0.996827</td>\n",
       "      <td>0.997765</td>\n",
       "      <td>0.998189</td>\n",
       "      <td>0.996722</td>\n",
       "      <td>0.998576</td>\n",
       "      <td>0.998386</td>\n",
       "      <td>0.996359</td>\n",
       "      <td>0.996780</td>\n",
       "      <td>0.998464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Grams_Subgroup_tofu__soy_or_meat_substitute</th>\n",
       "      <td>0.998563</td>\n",
       "      <td>0.999750</td>\n",
       "      <td>0.999831</td>\n",
       "      <td>-0.989504</td>\n",
       "      <td>-0.989975</td>\n",
       "      <td>0.999213</td>\n",
       "      <td>0.999770</td>\n",
       "      <td>-0.999453</td>\n",
       "      <td>0.999749</td>\n",
       "      <td>-0.954936</td>\n",
       "      <td>...</td>\n",
       "      <td>0.999797</td>\n",
       "      <td>0.999786</td>\n",
       "      <td>0.999649</td>\n",
       "      <td>0.999868</td>\n",
       "      <td>0.999492</td>\n",
       "      <td>0.999428</td>\n",
       "      <td>0.999393</td>\n",
       "      <td>0.999290</td>\n",
       "      <td>0.999164</td>\n",
       "      <td>0.999955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Grams_Subgroup_ultra_processed_legume_products</th>\n",
       "      <td>-0.977651</td>\n",
       "      <td>-0.983596</td>\n",
       "      <td>-0.984440</td>\n",
       "      <td>0.997927</td>\n",
       "      <td>0.998781</td>\n",
       "      <td>-0.978154</td>\n",
       "      <td>-0.981378</td>\n",
       "      <td>0.983175</td>\n",
       "      <td>-0.984510</td>\n",
       "      <td>0.982643</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.982713</td>\n",
       "      <td>-0.978651</td>\n",
       "      <td>-0.982049</td>\n",
       "      <td>-0.982830</td>\n",
       "      <td>-0.979172</td>\n",
       "      <td>-0.984667</td>\n",
       "      <td>-0.984531</td>\n",
       "      <td>-0.978103</td>\n",
       "      <td>-0.979625</td>\n",
       "      <td>-0.983354</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Processed_Grams_Subgroup_vegetables</th>\n",
       "      <td>-0.999801</td>\n",
       "      <td>-0.999593</td>\n",
       "      <td>-0.999383</td>\n",
       "      <td>0.987520</td>\n",
       "      <td>0.987127</td>\n",
       "      <td>-0.999990</td>\n",
       "      <td>-0.999844</td>\n",
       "      <td>0.999658</td>\n",
       "      <td>-0.999422</td>\n",
       "      <td>0.944376</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.999704</td>\n",
       "      <td>-0.999225</td>\n",
       "      <td>-0.999809</td>\n",
       "      <td>-0.999647</td>\n",
       "      <td>-0.999987</td>\n",
       "      <td>-0.997563</td>\n",
       "      <td>-0.999432</td>\n",
       "      <td>-0.998005</td>\n",
       "      <td>-0.997494</td>\n",
       "      <td>-0.998966</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>107 rows × 1783 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "featureid                                           G000006605  G000006785  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.999549    0.999796   \n",
       "Processed_Calories_Nova_processed_foods              -0.628812   -0.658344   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999649    0.999606   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.998660    0.999770   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.999571   -0.998286   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.999803    0.999793   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.995554    0.998346   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.998563    0.999750   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.977651   -0.983596   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.999801   -0.999593   \n",
       "\n",
       "featureid                                           G000006845  G000006865  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.999636   -0.991809   \n",
       "Processed_Calories_Nova_processed_foods              -0.664031    0.718708   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999437   -0.987259   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999835   -0.989393   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.997798    0.988823   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.999596   -0.990149   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.998689   -0.995988   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999831   -0.989504   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.984440    0.997927   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.999383    0.987520   \n",
       "\n",
       "featureid                                           G000006925  G000007265  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...   -0.991203    0.999559   \n",
       "Processed_Calories_Nova_processed_foods               0.732671   -0.635961   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...   -0.987037    0.999939   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...   -0.989820    0.999287   \n",
       "Processed_Calories_Nova_ultra_processed_foods         0.987195   -0.999025   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                  -0.989571    0.999852   \n",
       "Processed_Grams_Subgroup_table_sugar                 -0.996637    0.996189   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...   -0.989975    0.999213   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...    0.998781   -0.978154   \n",
       "Processed_Grams_Subgroup_vegetables                   0.987127   -0.999990   \n",
       "\n",
       "featureid                                           G000007325  G000007525  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.999693   -0.999956   \n",
       "Processed_Calories_Nova_processed_foods              -0.651258    0.652493   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999869   -0.999598   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999805   -0.999491   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.998394    0.998875   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.999838   -0.999939   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.997641   -0.997986   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999770   -0.999453   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.981378    0.983175   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.999844    0.999658   \n",
       "\n",
       "featureid                                           G000007745  G000007785  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.999745   -0.949400   \n",
       "Processed_Calories_Nova_processed_foods              -0.662388    0.855622   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999446   -0.945420   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999759   -0.954368   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.998057    0.936735   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.999686   -0.946938   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.998647   -0.967969   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999749   -0.954936   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.984510    0.982643   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.999422    0.944376   \n",
       "\n",
       "featureid                                           ...  G902794445  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...  ...    0.999749   \n",
       "Processed_Calories_Nova_processed_foods             ...   -0.655911   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...  ...    0.999729   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...  ...    0.999821   \n",
       "Processed_Calories_Nova_ultra_processed_foods       ...   -0.998288   \n",
       "...                                                 ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                 ...    0.999810   \n",
       "Processed_Grams_Subgroup_table_sugar                ...    0.998094   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...  ...    0.999797   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...  ...   -0.982713   \n",
       "Processed_Grams_Subgroup_vegetables                 ...   -0.999704   \n",
       "\n",
       "featureid                                           G902798515  G902798535  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.998335    0.999840   \n",
       "Processed_Calories_Nova_processed_foods              -0.656404   -0.651154   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999458    0.999795   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999795    0.999687   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.996197   -0.998667   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.998721    0.999921   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.996827    0.997765   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999786    0.999649   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.978651   -0.982049   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.999225   -0.999809   \n",
       "\n",
       "featureid                                           G902798655  G902798845  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.999658    0.999605   \n",
       "Processed_Calories_Nova_processed_foods              -0.657972   -0.641590   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999698    0.999976   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999887    0.999551   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.998044   -0.998756   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.999724    0.999863   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.998189    0.996722   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999868    0.999492   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.982830   -0.979172   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.999647   -0.999987   \n",
       "\n",
       "featureid                                           G902799635  G902800075  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.997464    0.999948   \n",
       "Processed_Calories_Nova_processed_foods              -0.686127   -0.657242   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.997899    0.999369   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999361    0.999421   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.993949   -0.998709   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.997471    0.999839   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.998576    0.998386   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999428    0.999393   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.984667   -0.984531   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.997563   -0.999432   \n",
       "\n",
       "featureid                                           G902801795  G902801825  \\\n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.996862    0.996540   \n",
       "Processed_Calories_Nova_processed_foods              -0.666889   -0.675265   \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.998393    0.997919   \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999263    0.999114   \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.993857   -0.993122   \n",
       "...                                                        ...         ...   \n",
       "Processed_Grams_Subgroup_table_salt                   0.997299    0.996883   \n",
       "Processed_Grams_Subgroup_table_sugar                  0.996359    0.996780   \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999290    0.999164   \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.978103   -0.979625   \n",
       "Processed_Grams_Subgroup_vegetables                  -0.998005   -0.997494   \n",
       "\n",
       "featureid                                           G902801975  \n",
       "Processed_Calories_Nova_processed_culinary_ingr...    0.998806  \n",
       "Processed_Calories_Nova_processed_foods              -0.670311  \n",
       "Processed_Calories_Nova_ratio_any_processed_unp...    0.999167  \n",
       "Processed_Calories_Nova_ratio_ultra_processed_u...    0.999935  \n",
       "Processed_Calories_Nova_ultra_processed_foods        -0.996286  \n",
       "...                                                        ...  \n",
       "Processed_Grams_Subgroup_table_salt                   0.998890  \n",
       "Processed_Grams_Subgroup_table_sugar                  0.998464  \n",
       "Processed_Grams_Subgroup_tofu__soy_or_meat_subs...    0.999955  \n",
       "Processed_Grams_Subgroup_ultra_processed_legume...   -0.983354  \n",
       "Processed_Grams_Subgroup_vegetables                  -0.998966  \n",
       "\n",
       "[107 rows x 1783 columns]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "corr_table.loc[corr_table.index.intersection(ultraprocessed)][corr_table.index.intersection(wol2_bt.ids(axis='observation'))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3d7d86b-17ae-479d-aa89-755f73d2de80",
   "metadata": {},
   "source": [
    "# metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8c21af8d-54c2-41a8-9473-b3bc2719191e",
   "metadata": {},
   "outputs": [],
   "source": [
    "wol2_table = q2.Artifact.load('wol2_table.qza').view(pd.DataFrame)\n",
    "food_table = q2.Artifact.load('food_table.qza').view(pd.DataFrame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e56acb34-6f43-4010-b9f3-17355641121f",
   "metadata": {},
   "outputs": [],
   "source": [
    "wol2_metadata = q2.Artifact.import_data(type='FeatureData[Taxonomy]', view='/projects/wol/qiyun/wol2/taxonomy/lineages.txt', view_type='HeaderlessTSVTaxonomyFormat') \n",
    "wol2_md = wol2_metadata.view(pd.DataFrame)\n",
    "wol2_md[['domain', 'phylum', 'class', 'order', 'family', 'genus', 'species']] = wol2_md['Taxon'].str.split(';', expand=True)\n",
    "wol2_md.drop(columns=['Taxon'], inplace=True) \n",
    "wol2_md.to_csv('wol2_taxonomy.tsv', sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "23d4c374-35ca-483e-9d5d-78cc2ef2a84f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fermented foods\n",
    "# non-fermented collapsed \n",
    "# GramsPerDay \n",
    "# Processed grams \n",
    "# Processed calories\n",
    "\n",
    "fermented = food_table.columns[food_table.columns.str.startswith('Collapsed_Fermented')]\n",
    "fermented_names = {i:'_'.join(i.split('_')[1:]) for i in fermented} \n",
    "collapsed = food_table.columns[food_table.columns.str.startswith('Collapsed_Foods')]\n",
    "collapsed_names = {i:'_'.join(i.split('_')[2:]) for i in collapsed}\n",
    "individual_foods = food_table.columns[food_table.columns.str.startswith('GramsPerDay')]\n",
    "individual_food_names = {i:'_'.join(i.split('_')[1:len(i.split('_'))-2]) for i in individual_foods}\n",
    "processed = food_table.columns[food_table.columns.str.startswith('Processed_')]\n",
    "processed_names = {i:i for i in processed}\n",
    "\n",
    "fermented_names.update(collapsed_names)\n",
    "fermented_names.update(individual_food_names)\n",
    "fermented_names.update(processed_names)\n",
    "food_md = pd.DataFrame.from_dict(fermented_names, orient='index', columns=['Name'])\n",
    "names = []\n",
    "for i in food_md.index: \n",
    "    if i in fermented: \n",
    "        names.append('Fermented') \n",
    "    elif i in collapsed: \n",
    "        names.append('Collapsed_Foods')\n",
    "    elif i in individual_foods: \n",
    "        names.append('Individual_Foods')\n",
    "    else: \n",
    "        names.append('Processed')\n",
    "food_md['Category'] = names\n",
    "food_md.rename_axis(index='Feature ID', inplace=True)\n",
    "food_md.to_csv('food_metadata.tsv', sep='\\t')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qiime2-2022.11",
   "language": "python",
   "name": "qiime2-2022.11"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}