{
"cells": [
{
"cell_type": "markdown",
"id": "116bcc5c",
"metadata": {},
"source": [
"# Summarizing complex landscapes\n",
"\n",
"In the previous section, we illustrate of how to infer a complete\n",
"combinatorial landscape from empirical data using either `VCregression` or `SeqDEFT`. \n",
"In this section, we illustrate how to obtain low-dimensional summary statistics of the\n",
"inferred landscapes that characterize the complexity and patterns of genetic interactions\n",
"across the different sequence positions. We do so by computing the variance explained by\n",
"genetic interactions of different orders, involving each possible site and pairs of sites.\n",
"\n",
"\n",
"## GB1 protein landscape\n",
"\n",
"In this section, we first show how we can use `gpmap-tools` to understand the complexity and patterns of genetic interactions across multiple sites in a previously inferred genotype-phenotype map (GB1) using `VCregression` and [previously published high throughput experimental data](https://doi.org/10.7554/eLife.16965)."
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "5e5d14fe",
"metadata": {},
"outputs": [],
"source": [
"# Import required libraries\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from gpmap.datasets import DataSet\n",
"from gpmap.summary import GPmapSummarizer"
]
},
{
"cell_type": "markdown",
"id": "fb22b7af",
"metadata": {},
"source": [
"We use the inferred complete combinatorial landscape as provided in `DataSet`s"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "8e0581bc",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" y \n",
" \n",
" \n",
" seq \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" AAAA \n",
" 0.296301 \n",
" \n",
" \n",
" AAAC \n",
" -2.713474 \n",
" \n",
" \n",
" AAAD \n",
" -2.912992 \n",
" \n",
" \n",
" AAAE \n",
" -4.548719 \n",
" \n",
" \n",
" AAAF \n",
" -3.276738 \n",
" \n",
" \n",
" ... \n",
" ... \n",
" \n",
" \n",
" YYYS \n",
" -4.662925 \n",
" \n",
" \n",
" YYYT \n",
" -3.223102 \n",
" \n",
" \n",
" YYYV \n",
" -3.001718 \n",
" \n",
" \n",
" YYYW \n",
" -4.723318 \n",
" \n",
" \n",
" YYYY \n",
" -4.876429 \n",
" \n",
" \n",
"
\n",
"
160000 rows × 1 columns
\n",
"
"
],
"text/plain": [
" y\n",
"seq \n",
"AAAA 0.296301\n",
"AAAC -2.713474\n",
"AAAD -2.912992\n",
"AAAE -4.548719\n",
"AAAF -3.276738\n",
"... ...\n",
"YYYS -4.662925\n",
"YYYT -3.223102\n",
"YYYV -3.001718\n",
"YYYW -4.723318\n",
"YYYY -4.876429\n",
"\n",
"[160000 rows x 1 columns]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gb1 = DataSet('gb1')\n",
"gb1.landscape"
]
},
{
"cell_type": "markdown",
"id": "0dc71915",
"metadata": {},
"source": [
"Firts, we define a `GPmapSummarizer` object with the configuration of sequence space and the associated functional values for every possible sequence `f`. \n",
"\n",
"> **Note**: Function values must be provided for sequences in alphabetical order."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "abd76f3b",
"metadata": {},
"outputs": [],
"source": [
"gb1s = GPmapSummarizer(n_alleles=20, seq_length=4, f=gb1.landscape['y'].values)"
]
},
{
"cell_type": "markdown",
"id": "ebe2678b",
"metadata": {},
"source": [
"### How to compute the variance components of interactions of order `k`\n",
"\n",
"Now, we only need to call the method `calc_V_k_variance_components` to obtain a `pd.DataFrame` with the total and percentage variance explained by interactions of every possible order"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "8519786b",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" k \n",
" variance \n",
" variance_perc \n",
" variance_perc_cum \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 1 \n",
" 180067.068327 \n",
" 52.811580 \n",
" 52.811580 \n",
" \n",
" \n",
" 1 \n",
" 2 \n",
" 128538.741722 \n",
" 37.698921 \n",
" 90.510501 \n",
" \n",
" \n",
" 2 \n",
" 3 \n",
" 31449.759748 \n",
" 9.223849 \n",
" 99.734350 \n",
" \n",
" \n",
" 3 \n",
" 4 \n",
" 905.762731 \n",
" 0.265650 \n",
" 100.000000 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" k variance variance_perc variance_perc_cum\n",
"0 1 180067.068327 52.811580 52.811580\n",
"1 2 128538.741722 37.698921 90.510501\n",
"2 3 31449.759748 9.223849 99.734350\n",
"3 4 905.762731 0.265650 100.000000"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v_k_vcs = gb1s.calc_V_k_variance_components()\n",
"v_k_vcs"
]
},
{
"cell_type": "markdown",
"id": "61d81198",
"metadata": {},
"source": [
"which can be easily represented graphically as follows"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "3d30ef73",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEhCAYAAADbKq0YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABClUlEQVR4nO3dd1RU1/YH8O+AzID0XhQUAUWkiA0RBA3EEkWIvjw1tthiDLYYY3lPfYkx6jP5JaZobIkxUWJMYhdRwxNEBFSKogiCoPReZ5A2c39/sLhxBJQpcCn7sxZrec9czuxB2d577jln8xiGYUAIIRxS4ToAQgihREQI4RwlIkII5ygREUI4R4mIEMI5SkSEEM5RIiKEcI4SESGEc5SICCGco0RECOEcp4no+vXr8PPzg4WFBXg8Hs6cOSP1OsMw2Lp1K8zNzaGhoQFfX1+kpqZKnVNaWoo5c+ZAR0cHenp6WLx4MYRCYQd+CkKIojhNRCKRCC4uLti7d2+Lr+/evRvffPMN9u/fj5iYGGhqamLixImoqalhz5kzZw4ePHiAq1ev4sKFC7h+/TrefffdjvoIhBBlYDoJAMzp06fZY4lEwpiZmTGff/4521ZeXs4IBALm119/ZRiGYZKSkhgAzO3bt9lzLl26xPB4PCYnJ6fDYieEKKYX14mwNRkZGcjPz4evry/bpqurCzc3N0RFRWHWrFmIioqCnp4eRowYwZ7j6+sLFRUVxMTE4M0332yx79raWtTW1rLHDQ0NePjwISwtLaGiQsNmpGNIJBKUlZVhyJAh6NWr0/4qdohO++nz8/MBAKamplLtpqam7Gv5+fkwMTGRer1Xr14wMDBgz2nJzp078cknnyg5YkLkk5CQABcXF67D4FSnTUTtadOmTVi7di17nJWVBUdHR9y6dQvm5uYcRkZ6kpKSEpw4cQJmZmZch8K5TpuImv5yCgoKpJJDQUEBhg4dyp5TWFgo9X0NDQ0oLS196V+uQCCAQCBgj3V1dQEA5ubm6Nu3r7I+AiEvJRAIoK6uTsMB6MTziKytrWFmZobQ0FC2rbKyEjExMXB3dwcAuLu7o7y8HLGxsew5//vf/yCRSODm5tbhMRNC5MPpFZFQKERaWhp7nJGRgYSEBBgYGMDKygpr1qzB9u3bYWdnB2tra2zZsgUWFhYICAgAAAwePBiTJk3C0qVLsX//ftTX12PFihWYNWsWLCwsOPpUhBCZcfnI7tq1awyAZl8LFixgGKbxEf6WLVsYU1NTRiAQMD4+PkxKSopUHyUlJczs2bMZLS0tRkdHh1m4cCFTVVUlUxxZWVkMACYrK0tZH42QVyosLGQ+/vhjprCwkOtQOMdjGNo8Pzs7G5aWlsjKyqIxIiKlpKQE8fHxqKiogK6uLlxdXWFoaKiUvouKirBv3z68//77MDY2btP3PH36FDdv3kRubi6EQiFmzpwJe3t79nWGYRAWFoa4uDjU1NTA0tISU6ZMkYr52bNnuHTpElJSUsDj8TB48GBMnjwZfD5fKZ9LHp12sJr0LO35Cy+v+Ph4nD9/Xqrt5s2bmDZtGvvApKPV1dXB1NQUQ4cOxcmTJ5u9HhkZiZiYGAQEBEBfXx/Xrl3DsWPHEBgYyM5VOnXqFKqqqjBv3jxIJBKcPXsW58+fx4wZMzr647AoERHOyfMLzzAMGIaBRCJp8etlr7Xl9crKSly9erXF9z537hysrKxgYGCg7B/FK9nZ2cHOzq7F1xiGQUxMDLy8vNirpICAAHzxxRdITk6Go6MjioqKkJaWhqVLl7LjqJMnT8bx48cxYcIEaGtrd9hneR4lIsKpkpISnDt3rsXXzp49iytXrgBAi4mES3FxcVKz/hVRV1cnNdNfVVVVrpnW5eXlEAqFGDBgANumrq6Ovn37snPlsrOzoa6uLvUwZ8CAAeDxeMjOzsbgwYMV+zByokREOBUWFvbS1589e9YxgciooqJCaX0dPnxY6tjb2xvjxo2TuZ+mXSc0NTWl2jU1NSESidhzXnxdRUUFGhoanO5aQYmIcCIrKwvXrl1DRkbGS89TU1ODtrY2VFRUmn3xeLwW2192Tlu/JzU19aWxNU2CVYYlS5bAyMiIPVZVVVVa310FJSLSoXJychAWFiY1f6w1PB4Po0aNUtotkCwGDhyIvXv3tnoLOGzYMKW9F5/Pl5rpLy8tLS0AjdvrPD/WIxKJ2DWbWlpa7NVRE4lEgmfPnrHfz4VOO7OadC95eXn49ddfcfjwYakkpKOj89LvU+YvvCwMDQ0xbdo08Hi8Zl/Tpk3jZKD6VfT09KClpYX09HS2rba2lp2eAgB9+/ZFTU0NcnNz2XMyMjLAMAynU1foioi0q4KCAoSFhSE5OVmqXVdXF15eXnBxcUFiYmKLA9Zc/8IPHToUVlZWiIuLY6cVDBs2jNOY6urqUFpayh6XlZUhPz8fGhoa7DY5ERERMDQ0hJ6eHq5duwZtbW32KZqxsTFsbW1x/vx5TJ06FWKxGMHBwXB0dOTsiRlAiYi0k6KiIoSHh+PBgwdS7dra2vDy8oKrqys7FtIZf+GbGBgYcHJr2Jrc3FwcPXqUPW56quji4oKAgAB4eHigvr4e58+fR01NDaysrDB37lypp3DTp09HcHAwfv75Z6kJjVyimdWgmdXKVFJSgvDwcCQmJkq1a2lpwdPTE8OHD+/xm4A1kWdmdXdF/yKIUpSVleH69eu4e/eu1ACvpqYmPDw8MGLECKipqXEYIenMKBERhZSXlyMiIgIJCQmQSCRsu4aGBjw8PDBy5EhO1zCRroESEZFLZWUlIiIiEBcXJ5WA1NXV4e7uDjc3N6U8kiY9AyUiIhOhUIgbN27gzp07EIvFbLtAIMDo0aMxevRoqKurcxgh6YooEZE2EYlEiIyMxO3bt9HQ0MC2q6mpwc3NDWPGjIGGhgaHEZKujBIReanq6mpERUUhJiYG9fX1bHuvXr0watQojBkzptnaJUJkRYmItKimpgZRUVGIjo5GXV0d266qqooRI0bA09OT0yUBpHuhRESk1NbWIiYmBlFRUVKlvVVVVTFs2DCMHTuW0xm4pHuiREQANC4duHXrFm7evCm19YaKigpcXV0xduxYpa44J+R5lIh6mBe3ZHV0dERGRgZu3LiB6upq9jwejwcXFxd4eXlBX1+fw4hJT0CJqAd5fkvWptnPkZGRUufweDw4OTnBy8uL8z2jSc9BiaiHKCkpwfnz51+6xaqjoyO8vLx6/Lon0vEoEfUQ8fHxL3196NCh8Pf376BoCJFGG6P1EBUVFa1eDfF4PKlJioR0NEpEPcSrnnjREzHCJUpEPcSrCgJytSUr6XokEgnS09Nx584dtgxSVVWV1MRXWdEYUQ9RWVnZrI3H4wHgfktW0nWUl5fj+PHjqKioQENDA2xsbCAQCHDjxg2IxWJMnTpVrn4pEfUQN27cYP88cOBA8Pn8TrUlK+kaQkJCYGFhgffeew+7d+9m2wcPHtysWq8sKBH1ANnZ2WyNLgMDA8ycORMqKnRXTmSXmZmJRYsWNau9pqen1+JVd1vRv8YeICIigv2zp6cnJSEiN4ZhWnz6WllZqdBGeG26Ivrmm2/a3OGqVavkDoYoX35+Ph49egSgsYaYs7MzxxGRrszGxgbR0dHw8/Nj2+rq6hAWFgY7Ozu5+21TIvrqq6+kjouKilBdXQ09PT0AjQNYvXv3homJCSWiTub5saExY8b0yHLGRHkmTJiAY8eOYe/evWhoaMCpU6dQUlKC3r17Y8aMGXL326ZE9HwN8KCgIOzbtw8//PADBg0aBABISUnB0qVLsWzZMrkDIcpXXFzM1hXT1NSkR/REYTo6Onjvvfdw//59FBQUoK6uDq6urnByclKoSovMgwVbtmzBt99+yyYhABg0aBC++uorbN68We5AWiIWi7FlyxZYW1tDQ0MDNjY2+PTTT6XuURmGwdatW2Fubg4NDQ34+voiNTVVqXF0Vc8vaB09ejSV8yFKoaKiAmdnZ7z++uuYMmUKhg0bpvC/LZkTUV5eXovLAcRiMQoKChQK5kX//e9/8f333+O7777Dw4cP8d///he7d+/Gt99+y56ze/dufPPNN9i/fz9iYmKgqamJiRMnSm3q1ROVl5fj3r17ABora4wcOZLjiEh3EBER0eK6xfj4eKlhAFnJnIh8fHywbNkyxMXFsW2xsbFYvny50kvz3rx5E/7+/pgyZQr69++Pf/zjH5gwYQJu3boFoPFqaM+ePdi8eTP8/f3h7OyMn3/+Gbm5uThz5oxSY+lqIiMj2TI/o0aNotI+RCliY2NhZGTUrN3Y2BixsbFy9ytzIvrxxx9hZmaGESNGQCAQQCAQYNSoUTA1NcXhw4flDqQlY8aMQWhoKPvU5+7du7hx4wZbpzsjIwP5+flSCVBXVxdubm6Iiopqtd/a2lpUVlayX1VVVUqNm2tVVVXs/1pNVTYIUQahUNjiXuWampoK/R7JPKHR2NgYwcHBePToEZKTkwEA9vb2GDhwoNxBtGbjxo2orKyEvb09VFVVIRaL8dlnn2HOnDkAGh9NA4CpqanU95mamrKvtWTnzp345JNPlB5vZxEVFcXWHBsxYgR69+7NcUSku9DV1UVWVlazXTszMzMV2stc7pnV/fv3B8MwsLGxQa9e7TNB++TJkzh+/DiCgoIwZMgQJCQkYM2aNbCwsMCCBQvk7nfTpk1Yu3Yte5yTkwMHBwdlhMy56upq3LlzB0Djhvfu7u4cR0S6k2HDhiEkJARisRjW1tYAGu9Mrl69qtC/NZkzSHV1NVauXImjR48CAB49eoQBAwZg5cqV6NOnDzZu3Ch3MC/66KOPsHHjRsyaNQsA4OTkhKdPn2Lnzp1YsGABzMzMAAAFBQUwNzdnv6+goOClq82bbimbKDI1vbN5vv6Yq6srVdwgSjVmzBhUV1cjODiYveru1asXPDw8MHbsWLn7lTkRbdq0CXfv3kVYWBgmTZrEtvv6+uLjjz9WaiKqrq5uthxBVVWVHYS1traGmZkZQkND2cRTWVmJmJgYLF++XGlxdBW1tbXsQL6Kigo8PDw4joh0NzweD6+//jq8vb1RVFQENTU1GBgYKHxXJPN3nzlzBr/99htGjx7NbiMBAEOGDMHjx48VCuZFfn5++Oyzz2BlZYUhQ4YgPj4eX375JRYtWgSg8YeyZs0abN++HXZ2drC2tsaWLVtgYWGBgIAApcbSFdy5c4edtuDs7MzOfCfdg0QiQVhYGBITEyEUCqGtrc1WWmn6XWQYBmFhYYiLi0NNTQ0sLS0xZcoUpRdC4PP56NOnj9L6kzkRFRUVwcTEpFm7SCSSSkzK8O2332LLli14//33UVhYCAsLCyxbtgxbt25lz1m/fj1EIhHeffddlJeXw9PTEyEhIVBXV1dqLJ1dfX291JNCuhrqfiIjI3Hnzh0EBATAxMQEubm5OHv2LNTV1dkno5GRkYiJiUFAQAD09fVx7do1HDt2DIGBgUoZy62rq8ONGzeQkZEBkUjUbAHs6tWr5epX5shGjBiBixcvYuXKlQD+3lzr8OHDSh8Y1dbWxp49e7Bnz55Wz+HxeNi2bRu2bdum1PfuauLj4yESiQA0Xp22NNeDdG1ZWVkYNGgQ+4RaT08P9+/fR05ODoDGq6GYmBh4eXnB3t4eABAQEIAvvvgCycnJcHR0VDiG8+fP48mTJ3B2dlbq+KPMiWjHjh2YPHkykpKS0NDQgK+//hpJSUm4efMmwsPDlRYYaTuxWCy1nMPT05PDaIis6urq2C1XgcZx0JauXiwtLREbG4uSkhIYGhoiPz8fmZmZmDBhAoDG2fRCoRADBgxgv0ddXR19+/ZFVlaWUhJRamoq3n77bVhZWSnc1/NkTkSenp5ISEjArl274OTkhCtXrmDYsGGIioqCk5OTUoMjbXPv3j32yZ+dnR37NJF0DS9OBPb29sa4ceOanefp6Yna2lp89913UFFRgUQiwWuvvcZu7SIUCgE0Ti58nqamJnu1rCgNDQ1oaGgopa/nyXXTaGNjg0OHDik7FiIHiUQitcZHkUeohBtLliyRupVubauWBw8eIDExETNmzICxsTHy8/Nx+fJlaGtrv7I4grKMHz8eYWFhCAgIUOoiarkSkUQiQVpaGgoLC9lH6U28vLyUEhhpm6SkJJSWlgJonGRqaWnJcUREVnw+v01rAa9evQoPDw/2FsvU1BQVFRW4ceMGhg4dyi69EIlEUuM3IpGo2eoDeUVFRaG0tBRffPEF9PT0mk2vkXcrIJkTUXR0NN5++208ffq02Yg5j8djJzmR9scwjNQ2sHQ11L3V19c3ezLN4/HY30M9PT1oaWkhPT2dvT2vra1FdnY2RowYoZQYnt/+R5lkTkTvvfce++TM3Nxc6Y/sSds9evQIhYWFAIA+ffqwU+5J9zRw4EBERERAV1cXJiYmyMvLQ3R0NHtbxuPx4ObmhoiICBgaGkJPTw/Xrl2DtrY2+xRNUS2NXSmDzIkoNTUVf/zxB2xtbdsjHtJGLV0N0X8K3dvkyZNx7do1BAcHs7dfw4cPh7e3N3uOh4cH6uvrcf78edTU1MDKygpz585tt/WgyiJzdG5ubkhLS6NExLGMjAx2/oipqWm77H5AOheBQIBJkyZJLa16EY/Hw/jx4zF+/Ph2iUEikSA6OhoPHjxARUVFs6GYDRs2yNWvzIlo5cqV+PDDD5Gfn9/iPrVUJaJjvFgiiK6GSEcIDw9HXFwc3N3dce3aNYwdOxbl5eVITk6WujKTlcyJqGmn/qb1XsDfA2Y0WN0xsrKy8OTJEwCAoaFht9nChHR+iYmJ8PPzw8CBAxEeHg5HR0cYGBjA1NQU2dnZcm/CJ3Mier6iB+HG81dDHh4eVDCRdBihUMhOBeDz+eyM8IEDB+LatWty9ytzIurXr5/cb0YUl5+fz1Yp0dXVpVth0qF0dHRQVVUFXV1d6Ovr4/HjxzA3N0dOTo5CNfPalIjOnTuHyZMnQ01NDefOnXvpudOmTZM7GPJqz18NUcFE0tHs7e2RkZGBvn37YtSoUTh9+jTi4+NRUVGB0aNHy91vmxJRQEAA8vPzYWJi8tJ9fmiMqH0VFxcjKSkJQOP6IVdXV44jIj3N84UqHB0doauri+zsbBgYGCg02bFNiej5ZRwvLukgHef5NWXu7u5UMJFwztLSUinLijr3LCfCerFgorKm7BPyKikpKbC1tYWqqipSUlJeeq68V0VyJSKRSITw8HBkZmairq5O6rVVq1bJFQh5ucjISHZNkZubGxVMJB3mxIkTWLduHTQ1NXHixIlWz+PxeFK7p8pC5kQUHx+PN954A9XV1RCJRDAwMEBxcTF69+4NExMTSkTt4PmCiXw+nwomkg71n//8p8U/K5PME1A++OAD+Pn5oaysDBoaGoiOjsbTp08xfPhwfPHFF+0RY4/3YsHE9tiYipBXEYvF+Pnnn1FSUqL0vmVORAkJCfjwww+hoqICVVVV1NbWwtLSErt378a//vUvpQfY01HBRNJZqKqqoqCgoF36ljkRqampsTN5TUxMkJmZCeDvUrREuZ4vmDhs2LAW644T0lGcnJzYYQJlknmMyNXVFbdv34adnR28vb2xdetWFBcX45dfflHK5tzkb1QwkXQ2EokEd+7cQXp6OszNzcHn86Venzhxolz9ylXFo6qqCgDw2WefYf78+Vi+fDns7Ozw448/yhUEadnt27elCibq6upyHBHp6YqKitjy7k1bFCuDXHXNmpiYmCAkJERpwZC/1dfXIzo6mj2mEkGkM1iwYEG79EvLtjupuLg4qYKJyi4ZTEhn0qYrIldX1zZvvBUXF6dQQKTxMenNmzfZY7oaIp1Jbm5uqzs0zpw5U64+27zolXScu3fvsgUTBw4cSAUTSadx//59nD59Gra2tnj8+DFsbGxQUlICoVCIwYMHy91vmxJRe82mJM1JJBKp8tFUIoh0JhEREZg4cSJGjRqFnTt3YtKkSdDT08OFCxcUmloi9xjRnTt38Msvv+CXX35BbGys3AEQaQ8ePGCfRlhbW6Nv374cR0TI38rKythCDaqqqqirqwOPx8Po0aMVGpaR+alZdnY2Zs+ejcjISOjp6QFoXBk+ZswYnDhxgn5xFMAwDJWPJp2auro6uz2strY2CgsLYWpqipqaGnbirTxkviJasmQJ6uvr8fDhQ5SWlqK0tBQPHz6ERCLBkiVL5A6ENG630FQwsW/fvujfvz+3ARHygn79+iE9PR0A4ODggJCQEJw7dw5//vmnQgU+Zb4iCg8Px82bN6X2HRk0aBC+/fZb+h9cAVQwkXQFb7zxBhoaGgAAXl5eUFVVRVZWFgYPHgwvLy+5+5U5EVlaWrZ4CSYWi2FhYSF3ID1deno6cnNzATQWTLSzs+M4IkKae37nBx6Pp7SpJTLfmn3++edYuXIluyIcaBy4Xr16dbtsA5KTk4O5c+fC0NAQGhoacHJyknpvhmGwdetWmJubQ0NDA76+vmyVi66EroZIV/Dzzz8jISGBHSdSFpmviN555x1UV1fDzc2Nrafd0NCAXr16YdGiRVKFFxVdi1JWVgYPDw+MHz8ely5dgrGxMVJTU6Gvr8+es3v3bnzzzTc4evQorK2tsWXLFkycOBFJSUlQV1dX6P07SmZmJp4+fQqgsWCiIvMxCGlPxsbGCA0NxcWLFzFw4EA4OTnBzs5O4WoyMieiPXv2KPSGsvjvf/8LS0tLHDlyhG17fkCMYRjs2bMHmzdvhr+/P4DGjG1qaoozZ85g1qxZHRarIp5/Uubp6UkFE0mnNXnyZEyaNAnp6elITEzEmTNnwOPx4ODgACcnJ7kfsMiciNpr0VtLzp07h4kTJ+Ktt95CeHg4+vTpg/fffx9Lly4F0Fh1Nj8/X6rEia6uLtzc3BAVFdVqIqqtrZW6tGzaTYALeXl5UgUTnZycOIuFkLbg8XiwsbGBjY0NGhoakJKSgoiICMTHx8u9Z7XM//X+9NNPLbY3NDRg06ZNcgXRmvT0dHz//fews7PD5cuXsXz5cqxatQpHjx4F0Fj1FABbAreJqakp+1pLdu7cCV1dXfaLy9rxz18NUcFE0pUIhULcuXMHkZGRKCgoUOhhlcxXRKtWrcLFixdx8OBBdqwmJSUFb7/9NkpKSrBz5065g3mRRCLBiBEjsGPHDgCNi2/v37+P/fv3K3RltmnTJqxdu5Y9zsnJ4SQZFRUVUcFE0qXU1tYiKSkJ9+/fx5MnT6Cvrw8nJyf84x//gIGBgdz9ylXFY+7cuXBycsKRI0fw6NEjrF+/HgEBAdi3b5/cgbTE3Ny8WYIYPHgw/vzzTwBgF4MWFBSwmzU1HQ8dOrTVfgUCgVQ5nqYFph3t+TVlVDCRtEVlZSX++usvpKWlob6+HgYGBvD392evRhiGQVhYGOLi4lBTUwNLS0tMmTJFadvIfPHFF1BXV8eQIUPg4+OjtCk7MiciGxsbREZGYs2aNZg0aRJUVVVx9OhRzJ49WykBPc/Dw6NZQbdHjx6hX79+ABoHrs3MzBAaGsomnsrKSsTExGD58uVKj0eZysrKqGAikcmzZ8/w448/wtraGnPmzEHv3r1RWloq9XQ4MjISMTExCAgIgL6+Pq5du4Zjx44hMDCQfcqtiFmzZmHAgAFKn14i1+OZixcv4sSJE3B3d4eenh5++OEHdjKeMn3wwQeIjo7Gjh07kJaWhqCgIBw8eBCBgYEAGgfN1qxZg+3bt+PcuXNITEzE/PnzYWFh0em3LqGCiURWkZGR0NXVhb+/P/r06QN9fX3Y2Niwt0QMwyAmJgZeXl6wt7eHqakpAgICUFVVheTkZKXEYGNj0y5z3GRORMuWLcNbb72FDRs2ICIiAvfu3QOfz4eTkxNOnjyp1OBGjhyJ06dP49dff4WjoyM+/fRT7NmzB3PmzGHPWb9+PVauXIl3330XI0eOhFAoREhISKeeQ1RVVYWEhAQAVDCRAHV1deyT3NraWnYJxYtSUlJgbm6O33//HZ9//jkOHDggtfNFeXk5hEIhBgwYwLapq6ujb9++nb7CjszXak2Xfi4uLgAax2mCg4Oxd+9eLFq0CP/85z+VGuDUqVMxderUVl/n8XjYtm0btm3bptT3bU83b96kgomEdfjwYaljb29vjBs3rtl5ZWVluHPnDtzd3eHp6Ync3FyEhIRAVVUVQ4cOhVAoBND44ON5mpqa7LbDnZXMiSg2NrbF24jAwECp+TykZdXV1ez/Yr169aKCiQRLliyBkZERe9zaFA6GYWBhYQEfHx8AjQ9zCgsLERsb+9KHM12BzLdmAoEAjx8/xubNmzF79mx224pLly61eklJ/hYdHc0uGnZ1daWCiQR8Pp99kisQCFodVNbW1oaxsbFUm5GRESoqKgCA/bf04tWPSCRqdpWkDMr8fZc5EYWHh8PJyQkxMTE4deoUezl49+5d2lL2FahgIlGEpaVls7rzJSUlbL07PT09aGlpsfsFAY3/5rKzs2FpaamUGBiGQXh4OL788kvs2LEDZWVlAID//e9/Cu3QKHMi2rhxI7Zv346rV69KVXl87bXXpOpwkeZu377NLi2hgolEVqNHj0Z2djYiIiJQWlqKxMRExMXFYeTIkQAax0vd3NwQERGBlJQUFBQU4PTp09DW1oa9vb1SYrh+/Tru3r0LX19fqVtIExMThUpRyzxGlJiYiKCgoGbtJiYmKC4uljuQ7q6+vh5RUVEAlLuPC+k5+vTpg5kzZyI0NBTh4eHQ19fHxIkT4ezszJ7j4eGB+vp6nD9/HjU1NbCyssLcuXOVMocIaLzzmTp1KgYMGICLFy+y7WZmZgr9/sscnZ6eHvLy8pptCxkfH48+ffrIHUh3FxcXh+rqagCNW2xSwUQij4EDB7Kb17eEx+Nh/PjxGD9+fLu8f1VVVYtLORiGaVbjTBYy35rNmjULGzZsQH5+Png8Hlv+Zt26dZg/f77cgXRnLxZMpC11SVdlbGyMzMzMZu1JSUlSy6xkJfMV0Y4dOxAYGAhLS0uIxWI4ODhALBbj7bffxubNm+UOpDt7sWDii7sFENJVeHl54cyZM6isrATDMHj48CGKi4tx7949hZZ58ZimdQYyysrKQmJiIoRCIVxdXbv0HstNTxWysrJaLIekyJR2FRUVrFixgr2cPXToEHJyctr0vXL+1ZAuoqioCPv27cP777/f7LF8Z/b06VNcv34d+fn5qKurg7m5Oby9vWFjYyN3n3KPYFlaWirtkWB3NmTIEDYJpaentzkJEdJZ9evXD/PmzVNqn7QnaTsxMDCAj48PpkyZwrZdv36dw4gIUdy5c+fw5MkTpfernGd6RMrQoUMxbdo0AGD3n2YYhq2MS0hXVV1djWPHjkFTUxNDhgyBs7Mzuy+YIigRKZmBgQGmTZvW4gb406ZNQ2ZmpsLVTQjhyqxZs/Ds2TMkJSUhMTER0dHRMDIygpOTE5ycnOT+z5ZuzZSste1emwa8aTtY0tVpaGhg+PDheOedd7BmzRq4uLjg3r17+Oabb+TuU65EFBERgblz58Ld3Z0dfP3ll1+kNoLvqV71PwLdnpHuQiwWIzc3Fzk5OSgvL1doAbfMt2Z//vkn5s2bhzlz5iA+Pp5dO1VRUYEdO3YgODhY7mC6g/LycoVeJ6Szy8jIQGJiIh4+fAiGYTB48GDMnj272WoLWciciLZv3479+/dj/vz5OHHiBNvu4eGB7du3yx1IdxEfHw8PDw8wDCM1/6hpTpAiCwMJ4dqXX36JZ8+ewdbWFn5+fhg4cKBS1rHJ3ENKSgq8vLyatevq6tL/9mgss33u3DlMmzYNEolEKhmdO3eOBqpJl+bt7Y0hQ4YofStmmRORmZkZ0tLSmpWWvXHjhtReuT1ZQkICMjMz4erqCj09PZSXlyM+Pp6SEOnyhg8f3i79ypyIli5ditWrV+PHH38Ej8dDbm4uoqKisG7dOmzZsqU9YuySSktLERoaynUYhCjst99+Q0BAAAQCAX777beXnjtz5ky53kPmRLRx40ZIJBL4+PiguroaXl5eEAgEWLduHVauXClXEISQzuv52zCBQNAu5YTkXvRaV1eHtLQ0CIVCODg4dOm9l9tz0asiaNFr99ZVF722B5nnEVVUVKC0tBR8Ph8ODg4YNWoUtLS0UFpaylnpZkJIxzh69ChqamqatdfW1uLo0aNy9yvXxmjPP7ZvcvLkScyaNUvuQAghnd+TJ09a3ImxoaGhxQ3T2krmRBQTE9PiNpTjxo1DTEyM3IEQQjqvgoICFBQUAGi8pWw6LigoQF5eHuLi4qCtrS13/zIPVrdWEre+vh7Pnj2TOxBCSOe1f/9+8Hg88Hi8Fm/B1NTUMHnyZLn7lzkRjRo1CgcPHsS3337bLND2mmNACOHW6tWrAQBff/01li5dit69e7OvqaqqQlNTs8UdJ9pKriUevr6+uHv3Llv6NjQ0FLdv38aVK1fkDoQQ0nk1LdZuryKqMiciDw8PREVF4fPPP8fJkyehoaEBZ2dn/PDDD11632pCSNsVFRWhoqKi2cD1oEGD5OpPrtVqQ4cOxfHjx+V6Q0JI11VWVobffvsNBQUF4PF47Fy3prl2W7dulatfuRKRRCJBWloaCgsLIZFIpF5raUEsIaR7CAkJgZ6eHubPn4+vv/4aS5YswbNnz3DlyhW8/vrrcvcrcyKKjo7G22+/jadPnzab+cvj8RSq9kgI6dyysrKwYMEC9O7dm32KZmVlBR8fH4SEhGDZsmVy9StzInrvvfcwYsQIXLx4Eebm5pwtfyCEdDyGYcDn8wEAvXv3RlVVFYyMjKCrq4vi4mK5+5X5eVtqaip27NiBwYMHQ09PD7q6ulJf7WnXrl3g8XhYs2YN21ZTU4PAwEAYGhpCS0sLM2bMYCdeEUKUy8TEhP396tOnD27evInMzExcv34d+vr6cvcrcyJyc3NDWlqa3G8or9u3b+PAgQNwdnaWav/ggw9w/vx5/P777wgPD0dubi6mT5/e4fER0hOMHTuWHZIZP348ysrKcOTIEaSmpnbshMaVK1fiww8/RH5+PpycnKCmpib1+ouJQhmEQiHmzJmDQ4cOSW1HW1FRgR9++AFBQUF47bXXAABHjhzB4MGDER0djdGjRys9FkJ6MltbW/bPBgYGWLFiBZ49ewZ1dXWFhmlkTkQzZswAACxatIhta3qM116D1YGBgZgyZQp8fX2lElFsbCzq6+vh6+vLttnb28PKygpRUVGtJqLa2lp2038AqKqqUnrMhPQUGhoaCvchcyLKyMhQ+E1lceLECcTFxeH27dvNXsvPzwefz29WosfU1BT5+fmt9rlz50588sknyg6VkA5148YNhIaGws3NDZMmTQLQuAr+8uXLePDgARoaGmBra4s33nhDof3CXrUr4/M6bIfGfv36yfVG8sjKysLq1atx9epVpW7WvWnTJqxdu5Y9zsnJgYODg9L6J6S95eTkIDY2FqamplLtISEhSE1NxVtvvQWBQIBLly7h5MmTUncwslL2RvktkbsOSFJSEjIzM1FXVyfV3lTzXRliY2NRWFiIYcOGsW1isRjXr1/Hd999h8uXL6Ourg7l5eVSV0UFBQUvrcctEAggEAjYY9rQjXQldXV1OHXqFPz8/HD9+nW2vaamBvHx8ZgxYwZbY8zf3x979+5FdnZ2i7uPtoW/v79S4n4ZmRNReno63nzzTSQmJrY4xVuZY0Q+Pj5ITEyUalu4cCHs7e2xYcMGWFpaQk1NDaGhoezYVUpKCjIzM+Hu7q60OAhpT3V1dVJjlqqqqi+tFRYcHAw7OzsMGDBAKhHl5eVBIpFIVdNpmuPT2jbInYXMiWj16tWwtrZGaGgorK2tcevWLZSUlODDDz/EF198odTgtLW14ejoKNWmqakJQ0NDtn3x4sVYu3YtDAwMoKOjg5UrV8Ld3Z2emJEu4/Dhw1LH3t7eGDduXIvn3r9/H3l5eVi6dGmz14RCIVRVVZvdSmlqakIoFCol1q+//vqlrzdtFyIrmRNRVFQU/ve//8HIyAgqKipQUVGBp6cndu7ciVWrVnV4JdOvvvoKKioqmDFjBmprazFx4kTs27evQ2PgCm3q3z0sWbIERkZG7LGqqmqL51VUVCAkJATz5s1TSnVVebi5uUkdSyQS5OfnIy0tDWPGjJG7X5k/jVgsZreENDIyQm5uLgYNGoR+/fohJSVF7kDaKiwsTOpYXV0de/fuxd69e9v9vQlpD3w+X2rMsjV5eXkQiUQ4cOAA28YwDJ4+fYpbt25h7ty5EIvFqKmpkboqEolESquy09qdxq1bt5CXlyd3vzInIkdHR9y9exfW1tZwc3PD7t27wefzcfDgQar0Skg7sra2xvLly6Xazp49CyMjI3h4eEBHRwcqKipIT09nnwIXFxejoqIClpaW7RqbnZ0dQkND5R7YljkRbd68GSKRCACwbds2TJ06FWPHjoWhoaFM8w0IIbIRCAQwMTGRalNTU4OGhgbb7urqiitXrkBDQ4N9fN+3b992H6hOSkpSaGKjzIlo4sSJ7J9tbW2RnJyM0tJS6Ovr00p8Qjg2adIkXL58GSdPnoRYLIaNjQ2mTJmitP6fvy1sIhQKIRKJFHofpYx4GRgYKKMbQoiM3nnnHanjXr16YcqUKUpNPs97cStYHo8HTU1N9O/fX2rAXVZtSkTTp0/HTz/9BB0dnVeubD916pTcwRBCOrfWphUoqk2JSFdXl73tau89hwghnZ9IJIJIJGo2lePFJSdt1aZEdOTIEQCNjwo/+eQTGBsbK2XFLSGka8nNzcWZM2dQXFzc4lbRHbJ5PsMwsLW1xYMHD6h0ECE90Llz52BoaIhp06YpbW4SIGMiUlFRgZ2dHUpKSigREdIDlZWV4Z///KfSH1DJvFXsrl278NFHH+H+/ftKDYQQ0vlZW1u/dK8vecn8+H7+/Pmorq6Gi4sL+Hx+s7Gi0tJSpQVHCOlcpk2bhjNnzqCwsBAmJibN1sV1WKXXPXv2yPVGhJCuLysrC5mZmUhNTW32WocNVgPAggUL5HojQkjXd+nSJTg7O8PLy4u7weoX1dTUNNuhUUdHR6GACCGd17NnzzB69GilJiFAjsFqkUiEFStWwMTEBJqamtDX15f6IoR0X4MHD8aTJ0+U3q/MV0Tr16/HtWvX8P3332PevHnYu3cvcnJycODAAezatUvpARJCOg8DAwOEhoYiMzOzxcHqFzdOayuZE9H58+fx888/Y9y4cVi4cCHGjh0LW1tb9OvXD8ePH8ecOXPkCoQQ0vnFx8eDz+fj6dOnePr0abPXOywRlZaWshug6ejosI/rPT09m23aRAjpXuTdk/pVZB4jGjBgAFtk0d7eHidPngTQeKX0YqFDQghpC5mviBYuXIi7d+/C29sbGzduhJ+fH7777jvU19fjyy+/bI8YCSGdxNmzZ1/6eodtFfvBBx+wf/b19UVycjJiY2Nha2sLZ2dnuYIghHQNNTU1UsdisRiFhYWoqalhizrKQ+ZElJWVJbURd79+/Tq0DDUhhDst1bZnGAYXLlxQaCGszGNE/fv3h7e3Nw4dOoSysjK535gQ0j3weDy4u7sjOjpa7j5kTkR37tzBqFGjsG3bNpibmyMgIAB//PGHVMlcQkjPUlZWBolEIvf3y3xr5urqCldXV+zevRthYWEICgrCu+++C4lEgunTp+PHH3+UOxhCSOd2+fJlqWOGYSAUCpGamgoXFxe5+5X5iqgJj8fD+PHjcejQIfz111+wtrbG0aNH5Q6EENL55efnS30VFhYCACZMmIBJkybJ3a/ci16zs7MRFBSEoKAg3L9/H+7u7lT2mZBurr1235A5ER04cABBQUGIjIyEvb095syZg7Nnz9KTM0J6gKaxIENDQ6n2kpISqKqqyj2pWeZbs+3bt8PNzQ2xsbG4f/8+Nm3aREmIkB7i7NmzyMrKataek5ODM2fOyN2vzFdEmZmZVFqakB4qLy8P06ZNa9bet29fBAcHy92vzFdElIQI6bl4PF6LU3Vqamqa1TmThdxPzQghPU+/fv1w48YNqTlDEokEN27cgJWVldz9KrRVLCGkZ/H19cWRI0fw3XffsWPDT58+RW1trUJP1Dr1FdHOnTsxcuRIaGtrw8TEBAEBAUhJSZE6p6amBoGBgTA0NISWlhZmzJiBgoICjiImpHszNjbG8uXLMWTIEIhEItTW1sLFxYXdPlpeCl0RFRcXIyYmBmKxGCNHjoS5ubki3TUTHh6OwMBAjBw5Eg0NDfjXv/6FCRMmICkpCZqamgAadwO4ePEifv/9d+jq6mLFihWYPn06IiMjlRoLIaSRtrY2fHx8lNqn3Inozz//xOLFizFw4EDU19cjJSUFe/fuxcKFC5UWXEhIiNTxTz/9BBMTE8TGxsLLywsVFRX44YcfEBQUhNdeew0AcOTIEQwePBjR0dEYPXq00mIhhGsRERFITk5GcXExevXqBUtLS/j6+sLIyIg9p6GhAZcvX8aDBw/Q0NAAW1tbvPHGG0qvuqFsbb41EwqFUseffPIJbt26hVu3biE+Ph6///47/v3vfys9wOdVVFQAALvdQGxsLOrr6+Hr68ueY29vDysrK0RFRbXaT21tLSorK9mvqqqqdo2bEGV4+vQpRo4cicWLF2PevHmQSCQ4duyYVEmvkJAQPHr0CG+99RbeeecdVFVVsbuodmZtTkTDhw+X2p2tV69e7DoTACgoKACfz1dudM+RSCRYs2YNPDw84OjoCKBx3Qufz282m9PU1PSl9bl37twJXV1d9svBwaHd4iZEWebOnYuhQ4fCxMQEZmZm8Pf3R0VFBfLy8gA0jpfGx8dj4sSJsLa2hoWFBfz9/ZGVlYXs7GyOo3+5Nieiy5cv4+DBg3jzzTeRm5uLr7/+GjNnzoSZmRmMjIywceNG7Nu3r90CDQwMxP3793HixAmF+9q0aRMqKirYr6SkJCVESIh86urqUFtby341NDS06fua5vNoaGgAaJxsKJFI2OIWAGBkZARdXd0WZ0N3Jm0eI+rfvz8uXryIX3/9Fd7e3li1ahXS0tKQlpYGsVgMe3t7qKurt0uQK1aswIULF3D9+nX07duXbTczM0NdXR3Ky8ulrooKCgpgZmbWan8CgQACgYA9rqysbJe4CWmLw4cPSx17e3tj3LhxL/0ehmEQEhICS0tL9mmVUCiEqqpqs99DTU3NZkMrylBdXY3s7GwwDAMLCwtoa2vL3ZfMg9WzZ8/G5MmTsW7dOowbNw4HDx7E0KFD5Q7gZRiGwcqVK3H69GmEhYU12xN3+PDhUFNTQ2hoKGbMmAEASElJQWZmJtzd3dslJkKUbcmSJVIDzi8WLWzJxYsXUVhYiEWLFrVnaK1KSkrCuXPnYGhoCIlEguLiYrzxxhtwdXWVqz+ZElFwcDAePnwIFxcXHD58GOHh4ZgzZw4mT56Mbdu2sZeIyhIYGIigoCCcPXsW2tra7LiPrq4uNDQ0oKuri8WLF2Pt2rUwMDCAjo4OVq5cCXd3d3piRroMPp8vdYX+KsHBwUhNTcU777wDHR0dtl1LSwtisRg1NTVSV0UikUjhp2Z1dXVSY8Dh4eFYunQpuwr/0aNHOH/+vNyJqM1jRB9++CEWLlyI27dvY9myZfj000/h7e2NuLg4qKurw9XVFZcuXZIriNZ8//33qKiowLhx42Bubs5+/fbbb+w5X331FaZOnYoZM2bAy8sLZmZmOHXqlFLjIKQzYBgGwcHBSE5Oxvz586Gvry/1urm5OVRUVJCens62FRcXo6KiQqrghTwOHjyI5ORk9lhFRQUikYg9FolEbbqSaw2PaeNKNUNDQ1y5cgXDhw9HaWkpRo8ejUePHrGvJyUlYdmyZYiIiJA7GK5kZ2fD0tISWVlZUmNQTbha6Puqv5rOGhdpm6KiIuzbtw/vv/8+jI2NX3n+xYsXkZiYiFmzZkndygkEAqipqQEALly4gLS0NPj7+0MgELAXB4sXL1Yo1vLycgQHB0NVVRVvvPEGSktL8eeff0IikUAikYDH4yEgIAB2dnZy9d/mWzNNTU1kZGRg+PDhyMrKajYg5uDg0CWTECFdxZ07dwCg2ZbM/v7+7DjtpEmTcPnyZZw8eRJisRg2NjaYMmWKwu+tp6eHt99+G4mJifjpp58watQorFy5EqWlpWAYBkZGRujVS/6FGm2+Ijp+/DiWLl0KPT09VFdX4+jRo3JXdexs6IpINnRFpByyXhF1FjU1Nbhy5QoKCwsxderUlz6hbqs2p7A5c+Zg0qRJSE9Ph52dHdW5J6SHSU1NRVFREczMzDBt2jQ8efIEp06dgq2tLcaPH8/eHspDptX3hoaGGDlyJCUhQnqYy5cv4+zZs8jNzcWFCxcQHh6O/v37Y9myZejVqxcOHDiA1NRUufun/YgIIa909+5dzJ07FxYWFnj27BkOHz4Mb29vqKqq4rXXXoOTkxMuXLgg92B1p96PiBDSOaipqaG8vBxA4+LzFwemjY2NFdp5g66ICCGv5OPjg9OnT+PSpUuor69HQECAUvunREQIeSVnZ2fY2tqirKwMhoaGSl9XSomIENImvXv3Ru/evdulbxojIoRwjhIRIYRzlIgIIZyjREQI4RwlIkII5ygREUI4R4mIEMI5mkdElI62JyGyoisiQgjnKBERQjhHiYgQwjlKRIQQzlEiIoRwjhIRIYRzlIgIIZyjREQI4RwlIkII5ygREUI4R4mIEMI5SkSEEM5RIiKEcI4SESGEc5SICCGco0RECOFct0lEe/fuRf/+/aGurg43NzfcunWL65AIaRe3bt3Cnj17sH37dhw+fBg5OTlch6SwbpGIfvvtN6xduxb/+c9/EBcXBxcXF0ycOBGFhYVch0aIUt2/fx9XrlyBt7c3li1bBlNTUxw7dgwikYjr0BTSLRLRl19+iaVLl2LhwoVwcHDA/v370bt3b/z4449ch0aIUkVHR2PYsGFwdXWFsbExpk6dCjU1NcTHx3MdmkK6/J7VdXV1iI2NxaZNm9g2FRUV+Pr6IioqqsXvqa2tRW1tLXtcUVEBAMjLy2vfYGWUnZ3NdQgt6qpxWVpadlAk0rKyslpsLykpQU1NDWpqaqT+PaqqqqJXr+a/mmKxGLm5ufD09GTbeDweBgwY0Gn/Ttqqyyei4uJiiMVimJqaSrWbmpoiOTm5xe/ZuXMnPvnkk2bto0aNapcY5cXVL86rUFyyaUtc6urq7J+9vb0xbty4ZudUV1eDYRhoampKtWtqaqK4uFjhOLnU5RORPDZt2oS1a9eyxw0NDXj48CEsLS2hoqK8u9Wqqio4ODggKSkJ2traSutXURSXbNorLolEgpKSEgwcOFDqCkhVVVVp79FVdPlEZGRkBFVVVRQUFEi1FxQUwMzMrMXvEQgEEAgEUm0eHh5Kj62yshIA0KdPH+jo6Ci9f3lRXLJpz7isrKzafG7v3r3B4/GaDUyLRCJoaWkpNa6O1uUHq/l8PoYPH47Q0FC2TSKRIDQ0FO7u7hxGRohyqaqqwsLCAunp6WwbwzBIT09H3759OYxMcV3+iggA1q5diwULFmDEiBEYNWoU9uzZA5FIhIULF3IdGiFKNXr0aJw5cwYWFhbo06cPoqOjUV9fj6FDh3IdmkK6RSKaOXMmioqKsHXrVuTn52Po0KEICQlpNoDd0QQCAf7zn/80uw3kGsUlm84Ul6OjI6qrqxEWFgahUAgzMzPMmTOny9+a8Riq00sI4ViXHyMihHR9lIgIIZyjREQI4RwlIkII5ygRtYPr16/Dz88PFhYW4PF4OHPmDNchAWhc2jJy5Ehoa2vDxMQEAQEBSElJ4TosfP/993B2doaOjg50dHTg7u6OS5cucR2WlF27doHH42HNmjVch9ItUSJqByKRCC4uLti7dy/XoUgJDw9HYGAgoqOjcfXqVdTX12PChAmcbyHRt29f7Nq1C7Gxsbhz5w5ee+01+Pv748GDB5zG1eT27ds4cOAAnJ2duQ6l+2JIuwLAnD59muswWlRYWMgAYMLDw7kOpRl9fX3m8OHDXIfBVFVVMXZ2dszVq1cZb29vZvXq1VyH1C3RFVEP1rT9iYGBAceR/E0sFuPEiRMQiUSdYolOYGAgpkyZAl9fX65D6da6xcxqIjuJRII1a9bAw8MDjo6OXIeDxMREuLu7o6amBlpaWjh9+jQcHBw4jenEiROIi4vD7du3OY2jJ6BE1EMFBgbi/v37uHHjBtehAAAGDRqEhIQEVFRU4I8//sCCBQsQHh7OWTLKysrC6tWrcfXqVam9gkj7oCUe7YzH4+H06dMICAjgOhTWihUrcPbsWVy/fh3W1tZch9MiX19f2NjY4MCBA5y8/5kzZ/Dmm29K7Q0kFovB4/GgoqKC2traHrlvUHuhK6IehGEYrFy5EqdPn0ZYWFinTUJA463j89undjQfHx8kJiZKtS1cuBD29vbYsGEDJSElo0TUDoRCIdLS0tjjjIwMJCQkwMDAQKaNsJQtMDAQQUFBOHv2LLS1tZGfnw8A0NXVhYaGBmdxbdq0CZMnT4aVlRWqqqoQFBSEsLAwXL58mbOYtLW1m42daWpqwtDQsFOMqXU7HD+165auXbvGAGj2tWDBAk7jaikmAMyRI0c4jWvRokVMv379GD6fzxgbGzM+Pj7MlStXOI2pJfT4vv3QGBEhhHM0j4gQwjlKRIQQzlEiIoRwjhIRIYRzlIgIIZyjREQI4RwlIkII5ygREUI4R4mIEMI5SkSkVePGjetSezS3d7zr1q3rVLsodCeUiJTgnXfekekfaGf7BW8tnlOnTuHTTz/t+IA6qYSEhC5fY76zokTUhdXV1bVr/wYGBtDW1m7X95BHe37ul/V99+5dSkTthBJROxg3bhxWrVqF9evXw8DAAGZmZvj4448BNF49hYeH4+uvvwaPxwOPx8OTJ08ANO7Bs3PnTlhbW0NDQwMuLi74448/pPpdsWIF1qxZAyMjI0ycOBEAEBISAk9PT+jp6cHQ0BBTp07F48ePpWKSSCTYvXs3bG1tIRAIYGVlhc8+++yl8Tx/pVRbW4tVq1bBxMQE6urq8PT0bLaF6ss+d2va2u+Ln1skEmH+/PnQ0tKCubk5/u///q9Z36/6eb7sZ/qi7OxsFBcXs4movLwcfn5+8PT0ZLdTIQrgevl/d7BgwQLG39+fPfb29mZ0dHSYjz/+mHn06BFz9OhRhsfjMVeuXGHKy8sZd3d3ZunSpUxeXh6Tl5fHNDQ0MAzDMNu3b2fs7e2ZkJAQ5vHjx8yRI0cYgUDAhIWFsf1qaWkxH330EZOcnMwkJyczDMMwf/zxB/Pnn38yqampTHx8POPn58c4OTkxYrGYjWn9+vWMvr4+89NPPzFpaWlMREQEc+jQoZfG8/y2F6tWrWIsLCyY4OBg5sGDB8yCBQsYfX19pqSkpE2fuzVt7ffFz718+XLGysqK+euvv5h79+4xU6dOZbS1taW26XjVz/NlP9MXnT9/ntHT02MYhmHu3bvH2NraMsuWLWPq6upa/Wyk7SgRKUFLicjT01PqnJEjRzIbNmxgX39xX5uamhqmd+/ezM2bN6XaFy9ezMyePZv9PldX11fGU1RUxABgEhMTGYZhmMrKSkYgEDCHDh1q8fzW9tlpahcKhYyamhpz/Phx9rW6ujrGwsKC2b17d5s/94tk6ff5z11VVcXw+Xzm5MmTbFtJSQmjoaHBfo62/Dxb6rs1n376KePt7c0cP36c0dfXZw4ePPjK7yFtRzs0tpMXi/GZm5ujsLCw1fPT0tJQXV2N119/Xaq9rq4Orq6u7PHw4cObfW9qaiq2bt2KmJgYFBcXQyKRAAAyMzPh6OiIhw8fora2Fj4+PnJ9lsePH6O+vh4eHh5sm5qaGkaNGoWHDx9KnSvL55al3+c/9+PHj1FXVwc3Nze2zcDAAIMGDWKP2/rzfLHv1iQkJODevXtYsWIFLl682ClKHXUnlIjaiZqamtQxj8djE0RLhEIhAODixYvo06eP1GsCgYD9s6amZrPv9fPzQ79+/XDo0CFYWFhAIpHA0dGRHXjtyG1gZf3cbdXS536Ztv4829p3QkICpk+fjqCgIJSXl8sUC3k1GqzmAJ/Ph1gslmpzcHCAQCBAZmYmbG1tpb4sLS1b7aukpAQpKSnYvHkzfHx8MHjwYJSVlUmdY2dnBw0NDYSGhrY5nufZ2NiAz+cjMjKSbauvr8ft27cVKvcjb782NjZQU1NDTEwM21ZWVoZHjx6xx/L+PFtSVVWF9PR0BAYG4rvvvsOsWbM6TTns7oKuiDjQv39/xMTE4MmTJ9DS0mIfk69btw4ffPABJBIJPD09UVFRgcjISOjo6GDBggUt9qWvrw9DQ0McPHgQ5ubmyMzMxMaNG6XOUVdXx4YNG7B+/Xrw+Xx4eHigqKgIDx48wOLFi1uMR0Xl7/+jNDU1sXz5cnz00UdsAYDdu3ejuroaixcvlvvnIG+/WlpaWLx4MT766CMYGhrCxMQE//73v6Vilvfn2ZK7d+9CVVUVDg4OcHV1xf379+Hn54dbt27ByMhI7s9P/kaJiAPr1q3DggUL4ODggGfPniEjIwP9+/fHp59+CmNjY+zcuRPp6enQ09PDsGHD8K9//avVvlRUVHDixAmsWrUKjo6OGDRoEL755huMGzdO6rwtW7agV69e2Lp1K3Jzc2Fubo733nvvpfE8b9euXZBIJJg3bx6qqqowYsQIXL58Gfr6+gr9LOTt9/PPP4dQKISfnx+0tbXx4YcfsiW0m8jz82xJQkIC7O3t2Vu6zz//HA8fPsT06dPx119/gc/ny/ahSTO0eT4hhHM0RkQI4RwlIkII5ygREUI4R4mIEMI5SkSEEM5RIiKEcI4SESGEc5SICCGco0RECOEcJSJCCOcoERFCOPf/mQhwpK2ne50AAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 1, figsize=(3, 3))\n",
"\n",
"# Plot percentage variance explained\n",
"axes.bar(x=v_k_vcs[\"k\"], height=v_k_vcs[\"variance_perc\"], color=\"black\")\n",
"axes.set(\n",
" xlabel=\"Interaction order $k$\",\n",
" ylabel=\"% variance explained\",\n",
" xticks=np.arange(1, 10),\n",
" ylim=(0, 100),\n",
")\n",
"\n",
"# Plot cumulative percentage variance explained\n",
"axes = axes.twinx()\n",
"axes.scatter(v_k_vcs[\"k\"], v_k_vcs[\"variance_perc_cum\"], color=\"grey\", s=25)\n",
"axes.plot(v_k_vcs[\"k\"], v_k_vcs[\"variance_perc_cum\"], color=\"grey\", lw=2)\n",
"axes.tick_params(axis=\"y\", colors=\"grey\")\n",
"axes.spines[\"right\"].set_color(\"grey\")\n",
"axes.set(ylim=(0, 105))\n",
"axes.set_ylabel(\"% cumulative variance\", color=\"grey\")\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "cff27023",
"metadata": {},
"source": [
"However, this only provides information about the complexity of the landscape, but not about which sites are involved in more or less complex genetic interactions and whether interactions are widespread or happening between specific sets of sites. \n",
"\n",
"### How to compute the variance components of interactions involving a subsets of sites $U$\n",
"\n",
"This can be done by calling the method `calc_V_U_variance_components` to obtain a `pd.DataFrame` with the total and percentage variance explained involving every possible combination of sites\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "4338153a",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" U \n",
" k \n",
" variance \n",
" variance_perc \n",
" variance_perc_cum \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" {0} \n",
" 1 \n",
" 28303.313131 \n",
" 8.301033 \n",
" 8.301033 \n",
" \n",
" \n",
" 1 \n",
" {1} \n",
" 1 \n",
" 16149.051370 \n",
" 4.736329 \n",
" 13.037362 \n",
" \n",
" \n",
" 2 \n",
" {2} \n",
" 1 \n",
" 77458.399775 \n",
" 22.717649 \n",
" 35.755012 \n",
" \n",
" \n",
" 3 \n",
" {3} \n",
" 1 \n",
" 58156.304052 \n",
" 17.056569 \n",
" 52.811580 \n",
" \n",
" \n",
" 4 \n",
" {0, 1} \n",
" 2 \n",
" 5138.287598 \n",
" 1.507000 \n",
" 54.318581 \n",
" \n",
" \n",
" 5 \n",
" {0, 2} \n",
" 2 \n",
" 25776.797274 \n",
" 7.560035 \n",
" 61.878616 \n",
" \n",
" \n",
" 6 \n",
" {0, 3} \n",
" 2 \n",
" 17682.891637 \n",
" 5.186187 \n",
" 67.064803 \n",
" \n",
" \n",
" 7 \n",
" {1, 2} \n",
" 2 \n",
" 8761.142103 \n",
" 2.569541 \n",
" 69.634344 \n",
" \n",
" \n",
" 8 \n",
" {1, 3} \n",
" 2 \n",
" 5422.459340 \n",
" 1.590344 \n",
" 71.224688 \n",
" \n",
" \n",
" 9 \n",
" {2, 3} \n",
" 2 \n",
" 65757.163770 \n",
" 19.285813 \n",
" 90.510501 \n",
" \n",
" \n",
" 10 \n",
" {0, 1, 2} \n",
" 3 \n",
" 3723.022582 \n",
" 1.091919 \n",
" 91.602420 \n",
" \n",
" \n",
" 11 \n",
" {0, 1, 3} \n",
" 3 \n",
" 2675.360421 \n",
" 0.784652 \n",
" 92.387072 \n",
" \n",
" \n",
" 12 \n",
" {0, 2, 3} \n",
" 3 \n",
" 18786.824440 \n",
" 5.509957 \n",
" 97.897030 \n",
" \n",
" \n",
" 13 \n",
" {1, 2, 3} \n",
" 3 \n",
" 6264.552305 \n",
" 1.837320 \n",
" 99.734350 \n",
" \n",
" \n",
" 14 \n",
" {0, 1, 2, 3} \n",
" 4 \n",
" 905.762731 \n",
" 0.265650 \n",
" 100.000000 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" U k variance variance_perc variance_perc_cum\n",
"0 {0} 1 28303.313131 8.301033 8.301033\n",
"1 {1} 1 16149.051370 4.736329 13.037362\n",
"2 {2} 1 77458.399775 22.717649 35.755012\n",
"3 {3} 1 58156.304052 17.056569 52.811580\n",
"4 {0, 1} 2 5138.287598 1.507000 54.318581\n",
"5 {0, 2} 2 25776.797274 7.560035 61.878616\n",
"6 {0, 3} 2 17682.891637 5.186187 67.064803\n",
"7 {1, 2} 2 8761.142103 2.569541 69.634344\n",
"8 {1, 3} 2 5422.459340 1.590344 71.224688\n",
"9 {2, 3} 2 65757.163770 19.285813 90.510501\n",
"10 {0, 1, 2} 3 3723.022582 1.091919 91.602420\n",
"11 {0, 1, 3} 3 2675.360421 0.784652 92.387072\n",
"12 {0, 2, 3} 3 18786.824440 5.509957 97.897030\n",
"13 {1, 2, 3} 3 6264.552305 1.837320 99.734350\n",
"14 {0, 1, 2, 3} 4 905.762731 0.265650 100.000000"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v_U_vcs = gb1s.calc_V_U_variance_components()\n",
"v_U_vcs\n"
]
},
{
"cell_type": "markdown",
"id": "68963729",
"metadata": {},
"source": [
"Whereas the number of possible combinations is relatively small in this case, it scales exponentially with sequence length, making it in practical in other cases to interpret the values in this `DataFrame` directly. Still, we can compute some useful aggregated statistics that we can represent easily as heatmaps.\n",
"\n",
"### How to compute the variance explained by interactions of order `k` involving each site\n",
"\n",
"`calc_sites_variance_perc` enables aggregation of these variances by whether they involve each possible site and the total number of sites involved in those genetic interactions. "
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "8c90e3d4",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" 1 \n",
" 2 \n",
" 3 \n",
" 4 \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 8.301033 \n",
" 14.253222 \n",
" 7.386529 \n",
" 0.26565 \n",
" \n",
" \n",
" 1 \n",
" 4.736329 \n",
" 5.666886 \n",
" 3.713892 \n",
" 0.26565 \n",
" \n",
" \n",
" 2 \n",
" 22.717649 \n",
" 29.415389 \n",
" 8.439197 \n",
" 0.26565 \n",
" \n",
" \n",
" 3 \n",
" 17.056569 \n",
" 26.062344 \n",
" 8.131930 \n",
" 0.26565 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 1 2 3 4\n",
"0 8.301033 14.253222 7.386529 0.26565\n",
"1 4.736329 5.666886 3.713892 0.26565\n",
"2 22.717649 29.415389 8.439197 0.26565\n",
"3 17.056569 26.062344 8.131930 0.26565"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sites = gb1s.calc_sites_variance_perc(v_U_vcs).T\n",
"sites"
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "5f0fba20",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAASUAAADoCAYAAACpZvKpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAArgElEQVR4nO3deVQUd7YH8G+zNDuNGBZRUBQCgwKCRILiksAoJM/EaDYlahBMwiAiiFHiFpUIMXFNiLghxidD3JOMRqMYUBkggOAyLgnICA6rsgkMW1PvDx4tLaDdRUNXwf2cU+dIdXn7lkevv/rVbxEwDMOAEEI4QkXZCRBCSEdUlAghnEJFiRDCKVSUCCGcQkWJEMIpVJQIIZxCRYkQwilqyk6gp1pbW1FUVAQ9PT0IBAJlp0NIlxiGwePHj2FmZgYVFWoLPAvvi1JRURHMzc2VnQYhMiksLMSwYcOUnQan8b4o6enpAQAKCgqgr6+v5GwI6VpNTQ0sLCwkf19J93hflNof2fT19akoEc6jLobno4dbQginUFEihHAKFSVCCKdQUSKEcAoVJUIIp1BRIoRwCu+HBBDS3zQ0NKCpqanbz4VCITQ1Nfswo75FRYkQDmloaICWltYzrzE1NUV+fn6/LUz0+EYIh7S3kAQCQbdHSUnJM1tSfEctJUI4qL0APY1hGPT3vT6oKBHCQaqqqt0WpdbWViVk9MTOnTtlvnbJkiVyx6eiRAgHqaiodFuUlG3btm1SP5eXl6O+vh4GBgYAgKqqKmhra8PY2JhVUaI+JUI4SEVFpdtD2fLz8yXHF198gbFjx+L27duoqKhARUUFbt++DWdnZ2zcuJFVfAHfN6OsqamBSCRCVVUVrRJAOKumpgYGBgaorq5+5t/T9r/POjo63baU6urqnhunr4waNQrHjh2Dk5OT1PmsrCy8/fbbyM/PlzsmPb4RwkFcfnzrqLi4GC0tLZ3Oi8VilJaWsoqp/LZgB1FRURAIBFi6dKmyUyFEqbj8+NaRh4cHPv74Y1y9elVyLisrCwEBAfD09GQVkzN3mJGRgd27d8PBwUHZqRCidM8ap8QlsbGxMDU1hYuLCzQ0NKChoYHx48fDxMQE+/btYxWTE49vtbW18PHxwd69exEREfHMaxsbG9HY2Cj5uaamprfTI6TPddcqUvZwgKcZGRnhzJkz+OOPP3Dnzh0AgK2tLV588UXWMTnRUgoMDMTrr78uU3MvMjISIpFIctCmAaQ/4svjW7sRI0bAxsYGr732Wo8KEsCyKJWVlfXoSztKSEjA1atXERkZKdP14eHhqK6ulhyFhYUKy4UQruBLUaqvr4efnx+0tbUxevRoFBQUAACCgoIQFRXFKiarO3z77bchFou7/KyrnvjuFBYWIjg4GIcPH5Z5cqGGhoZkkwDaLID0V3wpSuHh4bh27RqSkpKk/g17enrihx9+YBWT1R0aGBh0OVLz0aNHcvW4Z2VloaysDM7OzlBTU4OamhqSk5Oxc+dOqKmpdVv4COnv+FKUTp06hW+//Rbu7u5SnfCjR49GXl4eq5is7vD777/H+fPnERsbKzl3+/ZtjB8/Hjo6OjLH8fDwwI0bN5CTkyM5XFxc4OPjg5ycHKiqqrJJjxDe40tRKi8vh7GxcafzdXV1rN8Usnr7ZmBggOPHj2Pq1KkYM2YMKisr8d5778HPzw9fffWVzHH09PQwZswYqXM6OjoYPHhwp/OEDCTdvf7n2pAAFxcXnD59GkFBQQCe5Ldv3z64ubmxiilzUZo1axbGjh0rOezt7fHtt9/itddeQ0NDA7755hv4+vqySoIQIk1VVbXLVhHXitKmTZvg7e2NW7duoaWlBTt27MCtW7fwz3/+E8nJyaxiylyURo0ahcuXL+Pbb7/Fw4cPMWjQIDg6OoJhGMydOxfOzs5obm6Guro6q0TaJSUl9ej3E9IfcPFRrSvu7u7IyclBVFQU7O3t8euvv8LZ2Rmpqamwt7dnFZPVhNz//Oc/Uv1AOTk5uHfvHtTU1GBra4tr166xSoYNmpBL+EDeCblWVlZd9qmKxWLk5uZyZkJub2DVpzR06FAMHToUr7/+uuRcbW0tcnJy+rQgEdJf8aVPCWgbZZ6bm4uysrJOI84nT54sdzyFTTPR1dWFu7s73N3dFRWSkAFLVVWVF2+f09LSMHfuXNy/f7/TCgYCgYDVsB5OzH0jhEjjS1H65JNPJG/ghgwZopCWHBUlQjiIL49vf/75J44dOwYrKyuFxeR+9z4hA1B7S6mrg0tcXV2Rm5ur0Jhyt5Sam5vh5eWFmJgYWFtbKzQZQkgbLhagrgQFBWHZsmUoKSmBvb19pyFBbNZHk7soqaur4/r163J/ESFEdioqKl0WJa4thzt79mwAwMKFCyXnBAIBGIbp247uDz74APv372e9NAEh5Nn40qfEZmOA52FVlFpaWhAbG4sLFy5g3LhxnSbhbt26VSHJETJQdff4xrWW0vDhwxUek1VRunnzJpydnQEAf/zxh9RnXKvkhPARl4vSTz/9BG9vb6irq+Onn3565rVvvPGG3PFZFaXffvuNzW8jhMiIy49vM2fORElJCYyNjTFz5sxur6PBk4T0I1xuKXWcStIbGxmwLkqXL1/G7t27kZeXh2PHjmHo0KE4dOgQLC0tlTLVJDc3F7q6un3+vfIwMTFRdgoyu3//vrJTkAkfXpsDbXND5cHlotTbWBWl48ePY968efDx8UF2drZky6Pq6mps2rQJZ86cUWiShAw03S1dwsXlTOrq6pCcnIyCggI0NTVJfdbVstnPw6ooRUREICYmBvPnz0dCQoLk/MSJE5+7bxsh5PkEAgEvFnnLzs7Ga6+9hvr6etTV1cHQ0BAPHz6EtrY2jI2NWRUlVmX37t27XS5J0L6uESGkZ/gyzSQkJAQzZsxAZWUltLS0kJaWhvv372PcuHH4+uuvWcVkVZRMTU27nO9y5coVjBw5klUihJAn+FKUcnJysGzZMskI9MbGRpibm2Pz5s347LPPWMVkVZQWLVqE4OBgpKenQyAQoKioCIcPH0ZYWBgCAgJYJUIIeaJ9SEBXB5eoq6tLHjONjY0lm1GKRCLWG8Wy6lNauXIlWltb4eHhgfr6ekyePBkaGhoICwuT7GpACGGPLx3dTk5OyMjIgLW1NaZMmYK1a9fi4cOHOHToEOsdiVjdoUAgwKpVq1BRUYGbN28iLS0N5eXl2LhxI6skCCHS+LLv26ZNmzBkyBAAwBdffIFBgwYhICAA5eXl2LNnD6uYPRo8KRQKYWdn15MQhJAucHlEd0cuLi6SXxsbG+Ps2bM9jilzUQoNDZU5KE3IJaRn+PL41htkLkrZ2dlSP1+9ehUtLS2wsbEB0DYxV1VVFePGjVNshoQMQN29aeuNaR3ycnJykrnFdvXqVbnjy1yUOk7C3bp1K/T09HDw4EEMGjQIAFBZWQlfX19MmjRJ7iQIIdK4PHjyWZNwFYFVn9KWLVvw66+/SgoSAAwaNAgRERGYNm0ali1bprAECRmIFNGnFBkZiRMnTuDOnTvQ0tLChAkT8OWXX0qebgCgoaEBy5YtQ0JCAhobGzF9+nR89913z5ynuW7dOvluRk6sHlBrampQXl7e6Xx5eTkeP37c46QIGegU8fYtOTkZgYGBSEtLw/nz59Hc3Ixp06ahrq5Ock1ISAh+/vlnHD16FMnJySgqKsKsWbPkzjczMxOHDh3CoUOHkJWVJffv74hVS+mtt96Cr68vtmzZgvHjxwMA0tPTsXz5clY3RAiRpoiO7qffhMXFxcHY2BhZWVmYPHkyqqursX//fsTHx+PVV18FABw4cAB/+ctfkJaWhpdffvm53/HgwQPMmTMHKSkpMDAwAABUVVVhwoQJSEhIwLBhw2TOtx2rllJMTAy8vb0xd+5cDB8+HMOHD8fcuXPh5eWF7777jk1IQkgHzxvRXVNTI3W0r9TxLNXV1QAAQ0NDAEBWVhaam5vh6ekpucbW1hYWFhZITU2VKU9/f380Nzfj9u3bqKioQEVFBW7fvo3W1lb4+/vLe9sAWG6xNGPGDMTExOCrr75CXl4eAGDUqFGd1uomhLDzvJaSubm51Pl169bh888/7zZea2srli5diokTJ0pGWpeUlEAoFEpaOO1MTExQUlIiU57Jycn45z//KdVPZWNjg2+++Yb1S68ebbGko6PDal8nQsizPa8oFRYWQl9fX3JeQ0PjmfECAwNx8+ZNXLlyRaF5mpubo7m5udN5sVgMMzMzVjFZPb61b7GkCLt27YKDgwP09fWhr68PNzc3/PLLLwqJTQhfPa+ju/3fS/vxrKK0ePFi/OMf/8Bvv/0m1cdjamqKpqamTssNlZaWwtTUVKY8v/rqKwQFBSEzM1NyLjMzE8HBwayXLlH6FkvDhg1DVFQUrK2twTAMDh48iDfffBPZ2dkYPXo0m/QI4T1FDAlgGAZBQUE4efIkkpKSYGlpKfX5uHHjoK6ujsTERMmmknfv3kVBQQHc3Nxk+o4PP/wQ9fX1cHV1hZpaWzlpaWmBmpoaFi5cKLVJZUVFhUwxlb7F0owZM6R+/uKLL7Br1y6kpaVRUSIDliLevgUGBiI+Ph4//vgj9PT0JP1EIpEIWlpaEIlE8PPzQ2hoKAwNDaGvr4+goCC4ubnJ9OYNALZv3y5zPrLi1BZLYrEYR48eRV1dXbeVurGxUepNQ01NTa/kQogyKWJE965duwAAU6dOlTp/4MABfPjhhwCAbdu2QUVFBbNnz5YaPCmrBQsWyHytrDgxu+/GjRvQ1dWFhoYGPvnkE5w8ebLb1QciIyMhEokkx9NvIQjpDxSxyBvDMF0e7QUJADQ1NREdHY2KigrU1dXhxIkTMvcnAW1jn7rS0tKC8PBwmeN0xLooVVVVYcuWLfD394e/vz+2bt0qGQchLxsbG+Tk5CA9PR0BAQFYsGABbt261eW14eHhqK6ulhxsV7cjhMv4sp7SkiVL8M4776CyslJy7u7du3B1dcXf//53VjFZ3WFmZiZGjRqFbdu2SQZMbdu2DaNGjWI1K1goFMLKygrjxo1DZGQkHB0dsWPHji6v1dDQ6PTmgZD+hi9FKTs7Gw8ePIC9vT3Onz+P6OhoODs7w9bWFteuXWMVk1WfUkhICN544w3s3btXqsfd398fS5cuxaVLl1gl0661tVWmEaqE9FfdLV3CtY0DRo0ahZSUFCxduhReXl5QVVXFwYMHMWfOHNYxWRWlzMxMqYIEAGpqavj000+lVqKTRXh4OLy9vWFhYYHHjx8jPj4eSUlJOHfuHJvUCOkX+LLyJACcPn0aCQkJcHNzwx9//IH9+/djypQpfTt4Ul9fX7JrQUeFhYXQ09OTK1ZZWRnmz58PGxsbeHh4ICMjA+fOncNf//pXNqkR0i/w5fHt448/xjvvvIMVK1bg8uXLuH79OoRCIezt7XHkyBFWMVm1lN577z34+fnh66+/xoQJEwAAKSkpWL58udzNNkWNDCekP+HLcrgpKSlIT0+Ho6MjgLZR4mfOnEF0dDQWLlyId999V+6YrIrS119/DYFAgPnz56OlpQVA25y4gIAAREVFsQlJCOmAL49vWVlZXU5xCQwMlFp9QB6sipJQKMSOHTsQGRkptUqAtrY2qyQIIdL40lLS0NBAXl4eDhw4gLy8POzYsQPGxsb45ZdfYGFhwSpmj+5QW1sb9vb2sLe3p4JEiAK1b4P99MG1opScnAx7e3ukp6fjxIkTqK2tBQBcu3aN9bK53LpDQggA/nR0r1y5EhERETh//jyEQqHk/Kuvvoq0tDRWMbl1h4QQAIqZZtIXbty4gbfeeqvTeWNjYzx8+JBVTCpKhHAQX1pKBgYGKC4u7nQ+OzsbQ4cOZRWTW3dICAHAn6L0/vvvY8WKFSgpKYFAIEBraytSUlIQFhaG+fPns4rJ6u0bACQmJiIxMRFlZWWddu2MjY1lG5YQAv4MCdi0aRMCAwNhbm4OsVgMOzs7iMVizJ07F6tXr2YVk1VRWr9+PTZs2AAXFxcMGTKEc39QhPAdX+a+CYVC7N27F2vXrsWNGzdQW1sLJycnWFtbs47JqijFxMQgLi4O8+bNY/3FhJDu8WWcUjtzc3OFrW3Gqig1NTVJppcQQhSPb0VJkVjdob+/P+Lj4xWdCyHk//FlSEBvYNVSamhowJ49e3DhwgU4ODhAXV1d6nN5djMhhHSmiDW6+YpVUbp+/TrGjh0LoG1nk44Gwh8aIb1tID++cWo3k55ISUmBlpaWstN4Jm9vb2WnILP2/3S47um9zLjq6WEzz8OnonT58mXs3r0beXl5OHbsGIYOHYpDhw7B0tIS7u7ucsfj3h0SQngzePL48eOYPn06tLS0kJ2dLVnGurq6Gps2bWIVk/XgyaqqKuzfvx+3b98GANjZ2cHPzw8ikYhtSELI/+PLOKWIiAjExMRg/vz5SEhIkJyfOHEiIiIiWMXkxG4mhBBpfGkp3b17F5MnT+50XiQSoaqqilVMTu5mQshAx5dpJqampsjNzcWIESOkzl+5cgUjR45kFVPpu5kQQjrjy+PbokWLEBwcjNjYWAgEAhQVFSE1NRVhYWFYs2YNq5isilL7bia2trZS59nsZkII6Ywv45RWrlyJ1tZWeHh4oL6+HpMnT4aGhgbCwsIQFBTEKqbSdzMhhHTGl8c3gUCAVatWYfny5cjNzUVtbS3s7Oygq6vLOqZCdjNhGAZCoZB2MyFEQfjy+FZdXQ2xWAxDQ0PY2dlJzldUVEBNTQ36+vpyx2TVld++m0llZSVycnJw7do1yRu4rrZbIYTIhy9v395//32poQDtjhw5gvfff59VTJlbSqGhodi4cSN0dHQQGhr6zGtp7hshPcOXllJ6enqX/96nTp2KVatWsYopc1HKzs5Gc3Oz5Nfd4dozLyF8xJc+pcbGRsmGtB01Nzfjv//9L6uYMheljvPduDj3jZD+hC9z38aPH489e/bgm2++kTofExODcePGsYrJqqO7oKAA5ubmXVbtgoIC1jtjEkLa8OXxLSIiAp6enrh27Ro8PDwAtK3fn5GRgV9//ZVVTFZl19LSEuXl5Z3OP3r0iDeztgnhMr4s8jZx4kSkpqbC3NwcR44cwc8//wwrKytcv34dkyZNYhWTVUuJYZgu/3Bqa2uhqanJKhFCyBN8eXwD2pa5OXz4sMLiyVWU2t+6CQQCrFmzBtra2pLPxGIx0tPTebMODyFcpqKi0uWjGheLUmtrK3Jzc7vcbq2rybrPI1dRan/rxjAMbty4IbV3uFAohKOjI8LCwuRKIDIyEidOnMCdO3egpaWFCRMm4Msvv4SNjY1ccQjpT/jSUkpLS8PcuXNx//59MAwj9ZlAIIBYLJY7plxFqf2tm6+vL3bu3KmQeW7JyckIDAzESy+9hJaWFnz22WeYNm0abt26BR0dnR7HJ4SP+FKUPvnkE7i4uOD06dMK2wOSVZ+StbU1jh49ioULF0qdj42NRXl5OVasWCFzrLNnz0r9HBcXB2NjY2RlZbFq+hHSH/BlnNKff/6JY8eOwcrKSmExWZXdPXv2dFohAABGjx6NmJiYHiVUXV0NADA0NOzy88bGRtTU1EgdhPQ3ippmcunSJcyYMQNmZmYQCAQ4deqU1OcMw2Dt2rUYMmQItLS04OnpiT///FPm+K6ursjNzZUrp+dhVZRKSkowZMiQTueNjIxQXFzMOpnW1lYsXboUEydOxJgxY7q8JjIyEiKRSHIoaldOQrhEUUMC6urq4OjoiOjo6C4/37x5M3bu3ImYmBikp6dDR0cH06dPR0NDg0zxg4KCsGzZMsTFxSErKwvXr1+XOthg9fhmbm6OlJSUTmOSUlJSYGZmxioRAAgMDMTNmzdx5cqVbq8JDw+XmntXU1NDhYn0O4rqU/L29u52Fx2GYbB9+3asXr0ab775JgDg+++/h4mJCU6dOiXThNrZs2cDgFRXjkAgkAwb6vWO7naLFi3C0qVL0dzcjFdffRVA2yjOTz/9FMuWLWMTEosXL8Y//vEPXLp0CcOGDev2Og0NDVqJgPR7zytKT3dbsPl3kZ+fj5KSEnh6ekrOiUQiuLq6IjU1VaailJ+fL9d3yoJVUVq+fDkePXqEv/3tb2hqagIAaGpqYsWKFQgPD5crFsMwCAoKwsmTJ5GUlEQjwgnB8zu6n346WLduHT7//HO5vqOkpAQAYGJiInXexMRE8tnzDB8+XK7vlAWroiQQCPDll19izZo1uH37NrS0tGBtbc2qBRMYGIj4+Hj8+OOP0NPTk/xhiEQizm8uSUhveV5RKiwslFpATdlPD7du3UJBQYGkkdLujTfekDsW633fAEBXVxcvvfRST0Jg165dANrWX+nowIED+PDDD3sUmxC+et4a3fr6+qxWdezI1NQUAFBaWir14qq0tFTmmRn37t3DW2+9hRs3bkj6kjrm2Wd9Su0UUR2fHgVKCOmbwZOWlpYwNTVFYmKipAjV1NQgPT0dAQEBMsUIDg6GpaUlEhMTYWlpid9//x2PHj3CsmXL8PXXX7PKi1VR6o3qSAh5QlGDJ2tra6XGEeXn5yMnJweGhoawsLDA0qVLERERAWtra1haWmLNmjUwMzPDzJkzZYqfmpqKixcv4oUXXpAUUnd3d0RGRmLJkiXPXBCyO6zKbnt1LCsrg7a2Nv71r3/h0qVLcHFxQVJSEpuQhJAOFDVOKTMzE05OTnBycgLQNqneyckJa9euBQB8+umnCAoKwkcffYSXXnoJtbW1OHv2rMyrfYjFYsl0sxdeeAFFRUUA2jrA7969K1eu7Vi1lHqjOhJCnlDU49vUqVOf2UUiEAiwYcMGbNiwQe4cAWDMmDG4du0aLC0t4erqis2bN0MoFGLPnj2sd8hl1VLqjepICHmCL4u8rV69WrJcyYYNG5Cfn49JkybhzJkz2LlzJ6uYrFpKvVEdCSFP8GWVgOnTp0t+bWVlhTt37qCiogKDBg1iXUBZFaXVq1ejrq4OQFt1/J//+R9MmjQJgwcPxg8//MAqEULIE3wpSl3pbjK9rFgVpd6ojoSQJ7i8dMmsWbMQFxcHfX19zJo165nXnjhxQu74chel5uZmeHl5ISYmBtbW1pLzPa2OhJAnuFyURCKRJA+RSKTw+HIXJXV1ddZLEhBCZMPlx7cDBw4AaBv4vH79ehgZGSl0ShirO/zggw+wf/9+hSVBCJGmqEXeehPDMLCyssKDBw8UGpdVn1JLSwtiY2Nx4cIFjBs3rtNa2l3tLU4I6V9UVFRgbW2NR48eSXXl9BSronTz5k04OzsDAP744w+pz7jwzEsI3z1vQi5XREVFYfny5di1a1e3q8XKi1VRat/VhBDSO7jc0d3R/PnzUV9fD0dHRwiFwk59SxUVFXLHZL1KwOXLl7F7927cu3cPR48exdChQ3Ho0CFYWlrC3d2dbVhCCLjd0d3R9u3bFR6TVVE6fvw45s2bBx8fH1y9ehWNjY0A2nYi2bRpE86cOaPQJGXh4uICXV3dPv9eeXTcvJPrulvXmWvGjx+v7BRk0tjYiKioKGWnoXALFixQeExWZTciIgIxMTHYu3cv1NXVJecnTpyIq1evKiw5QgYqPrx9e1pDQ4NCtj9jdYd3797tcqNIkUiEqqoqVokQQp7gy4Tcuro6LF68GMbGxtDR0cGgQYOkDjZYFSVTU9MuN6C7cuUKTcglRAH4UpQ+/fRTXLx4Ebt27YKGhgb27duH9evXw8zMDN9//z2rmKy3WAoODkZsbCwEAgGKioqQmpqKsLAwrFmzhlUihJAn+PL27eeff8b333+PqVOnwtfXF5MmTYKVlRWGDx+Ow4cPw8fHR+6YrIrSypUr0draCg8PD9TX12Py5MnQ0NBAWFgYgoKC2IQkhHTAl7dvFRUVkqcjfX19yRAAd3d3mdf5fhqrOxQIBFi1ahUqKipw8+ZNpKWloby8HBs3bmSVBCFEGl8e30aOHCnZkNLW1hZHjhwB0NaCMjAwYBWTVVEqKCgAwzAQCoWws7PD+PHjJa/jCwoKWCVCCHmCL0XJ19cX165dA9D2BBUdHQ1NTU2EhIRg+fLlrGKyenyztLREcXExjI2Npc4/evQIlpaWtJsJIT3Elz6lkJAQya89PT1x584dZGVlwcrKCg4ODqxisipKDMN0+YdTW1sr8y4IhJDu8aVPqbCwUGoL8eHDh/d4K2+5ilJoaCiAtmq9Zs0aaGtrSz4Ti8VIT0+XeWdNQkj3+NJSGjFiBNzd3fHBBx/g7bffZj02qSO5ilL71kkMw+DGjRtS0yaEQiEcHR0RFhbW46QIIdwrQF3JzMxEfHw8NmzYgKCgIHh5eeGDDz7AjBkzoKGhwSqmXEWpfXUAX19f7Ny5U7LNEiFEsfjSUmrf6HLz5s1ISkpCfHw8PvroI7S2tmLWrFmIjY2VOyarPqUDBw4gMTERiYmJKCsrk+z71I5NIoSQJ/jSp9ROIBDglVdewSuvvIKAgAD4+fnh4MGDrGoBqzvcsGEDpk2bhsTERDx8+BCVlZVSByFkYHnw4AE2b96MsWPHSoYIRUdHs4rFqqW0a9cuxMXFYd68eay+lBDybHx5fNu9ezfi4+ORkpICW1tb+Pj44Mcff+zRGzhWRampqQkTJkxg/aWEkGfjS1GKiIjAnDlzsHPnTjg6OiokJqui5O/vj/j4eJp8S0gv4UtRKigoUHhOrIpSQ0MD9uzZgwsXLsDBwUFqoTdAvt1MLl26hK+++gpZWVkoLi7GyZMnMXPmTDZpEdJv8KUo9UY+rIrS9evXJYMkb968KfWZvEnW1dXB0dERCxcufO4WwIQMFHwpSr1B6buZeHt782Y9aEL6ChUlHmlsbJRsVACA9TrAhBBukqsoyfp4deLECVbJyCIyMhLr16/vtfiEcAEfW0oPHz5Eeno6xGIxXnrpJQwZMoRVHLmKkkgkYvUlihQeHi6ZGAy0tZQ6zlImpD/gW1E6fvw4/Pz88OKLL6K5uRl3795FdHQ0fH195Y4lV1E6cOCA3F+gaBoaGqwn+hFCFKO2tlZqn8X169fj999/x4svvggAOH36NBYtWsSqKHFzIg0hAxzXV54cN24cfvzxR8nPampqKCsrk/xcWlrKevNVpXd019bWSm3XlJ+fj5ycHBgaGsLCwkKJmRGiPFx/fDt37hwCAwMRFxeH6Oho7NixA++99x7EYjFaWlqgoqKCuLg4VrGVXpQyMzPxyiuvSH5u7y9asGAB65sihO+4XpRGjBiB06dP4+9//zumTJmCJUuWIDc3F7m5uRCLxbC1tWW9Cq3SH9+mTp0KhmE6HVSQyEDG9ce3dnPmzEFGRgauXbuGqVOnorW1FWPHju3RsthKL0qEkM4UWZSio6MxYsQIaGpqwtXVFb///rtCcjxz5gy2bNmCzMxM7Nu3D5s3b4aPjw+WL1+O//73v6zjUlEihIMUVZR++OEHhIaGYt26dbh69SocHR0xffp0qU5pNpYtWwZfX19kZGTg448/xsaNGzFlyhRcvXoVmpqacHJywi+//MIqNhUlQvqxrVu3Sl7N29nZISYmBtra2j1eHTYuLg5nzpxBQkICMjIycOjQIQBta/Vv3LgRJ06cwKZNm1jFpqJECAc9fvy42wNoGzTc8eg49apdU1MTsrKy4OnpKTmnoqICT09PpKam9ig/HR0dyc64hYWFnfqQ7OzscPnyZVaxqSgRwiFCoRCmpqYwNzeHSCTqdJibm0NXV7fT55GRkZ1iPXz4EGKxGCYmJlLnTUxMUFJS0qM8IyMjMX/+fJiZmWHKlCnYuHFjj+J1pPQhAYSQJzQ1NZGfn4+mpqZur+lqM9i+nuXg4+MDLy8v3Lt3D9bW1jAwMFBYbCpKhHCMpqamQnaafuGFF6CqqorS0lKp86WlpTA1Ne1x/MGDB2Pw4ME9jvM0enwjpJ8SCoUYN24cEhMTJedaW1uRmJgINzc3JWb2bNRSIqQfCw0NxYIFC+Di4oLx48dj+/btqKurYzVRtq9QUSKkH3vvvfdQXl6OtWvXoqSkBGPHjsXZs2c7dX5zCRUlQvq5xYsXY/HixcpOQ2bUp0QI4RQqSoQQTuH94xvDMADatmriuvbRuHzQ3Nys7BRk0tVIZi5qz7P97yvpnoDh+Z/SgwcPaI1uwhuFhYUYNmyYstPgNN4XpdbWVhQVFUFPT09ha820b0ZQWFgIfX19hcTsLXzJdaDnyTAMHj9+DDMzM6ioUK/Js/D+8U1FRaXX/ufR19fn9D+gjviS60DOkwu7AfEBlWxCCKdQUSKEcAoVpS5oaGhg3bp1vNhfji+5Up5EVrzv6CaE9C/UUiKEcAoVJUIIp1BRIoRwChUlQginDOiitGvXLjg4OEgGyrm5uUntVZWXl4e33noLRkZG0NfXx7vvvttpaVFliIqKgkAgwNKlSyXnGhoaEBgYiMGDB0NXVxezZ89Weq5d5blnzx5MnToV+vr6EAgEqKqqUkpun3/+eaf91GxtbTtdxzAMvL29IRAIcOrUqb5PdAAa0EVp2LBhiIqKQlZWFjIzM/Hqq6/izTffxL/+9S/U1dVh2rRpEAgEuHjxIlJSUtDU1IQZM2agtbVVaTlnZGRg9+7dcHBwkDofEhKCn3/+GUePHkVycjKKioowa9YsJWXZfZ719fXw8vLCZ599pqTMnhg9ejSKi4slx5UrVzpds337ds5tld3vMUTKoEGDmH379jHnzp1jVFRUmOrqaslnVVVVjEAgYM6fP6+U3B4/fsxYW1sz58+fZ6ZMmcIEBwdL8lJXV2eOHj0qufb27dsMACY1NZUzeXb022+/MQCYysrKPs+PYRhm3bp1jKOj4zOvyc7OZoYOHcoUFxczAJiTJ0/2SW4D3YBuKXUkFouRkJCAuro6uLm5obGxEQKBQGoQnaamJlRUVLr8H7UvBAYG4vXXX5faXBAAsrKy0NzcLHXe1tYWFhYWPd50kI3u8uSaP//8E2ZmZhg5ciR8fHxQUFAg+ay+vh5z585FdHS0Qnb+ILLj/YTcnrpx4wbc3NzQ0NAAXV1dnDx5EnZ2djAyMoKOjg5WrFiBTZs2gWEYrFy5EmKxGMXFxX2eZ0JCAq5evYqMjIxOn5WUlEAoFHbae0sRmw7K61l5comrqyvi4uJgY2OD4uJirF+/HpMmTcLNmzehp6eHkJAQTJgwAW+++aayUx1wBnxRsrGxQU5ODqqrq3Hs2DEsWLAAycnJsLOzw9GjRxEQEICdO3dCRUUFc+bMgbOzc58vPVFYWIjg4GCcP39eIfuB9Ra+5AkA3t7ekl87ODjA1dUVw4cPx5EjR2BkZISLFy8iOztbiRkOYMp+fuQaDw8P5qOPPpI6V15eLun7MDExYTZv3tynOZ08eZIBwKiqqkoOAIxAIGBUVVWZCxcudNk/Y2FhwWzdupUzeba0tEiuVXafUldcXFyYlStXMsHBwZKcO96HiooKM2XKFGWn2e8N+JbS01pbWzstsfrCCy8AAC5evIiysjK88cYbfZqTh4cHbty4IXXO19cXtra2WLFiBczNzaGuro7ExETMnj0bAHD37l0UFBT06aaDz8tTVVW1z3KRV21tLfLy8jBv3jy8++678Pf3l/rc3t4e27Ztw4wZM5SU4cAxoItSeHg4vL29YWFhgcePHyM+Ph5JSUk4d+4cAODAgQP4y1/+AiMjI6SmpiI4OBghISGwsbHp0zz19PQwZswYqXM6OjoYPHiw5Lyfnx9CQ0NhaGgIfX19BAUFwc3NDS+//DKn8iwpKUFJSQlyc3MBtPXp6enpwcLCAoaGhn2Wa1hYGGbMmIHhw4ejqKgI69atg6qqKubMmQMjI6MuO7ctLCxgaWnZZzkOVAO6KJWVlWH+/PkoLi6GSCSCg4MDzp07h7/+9a8A2lob4eHhqKiowIgRI7Bq1SqEhIQoOeuubdu2DSoqKpg9ezYaGxsxffp0fPfdd8pOq5OYmBisX79e8vPkyZMBtP0H8OGHH/ZZHg8ePMCcOXPw6NEjGBkZwd3dHWlpaTAyMuqzHEjXaOkSQgin0DglQginUFEihHAKFSVCCKdQUSKEcAoVJUIIp1BRIoRwChUlQginUFEihHAKFSUCALTcK+EMKkoDRHl5OQICAmBhYQENDQ2Ymppi+vTpSElJAQAUFxdLlvP497//DYFAgJycHCVmTAaqAT33bSCZPXs2mpqacPDgQYwcORKlpaVITEzEo0ePAIBWVyTcoey1U0jvq6ysZAAwSUlJ3V6DDmtQA5A6Oq4htHfvXsbW1pbR0NBgbGxsmOjo6F7Ongw01FIaAHR1daGrq4tTp07h5Zdfllp3vCu///47xo8fjwsXLmD06NEQCoUAgMOHD2Pt2rX49ttv4eTkhOzsbCxatAg6OjpYsGBBX9wKGQCoT2kAUFNTQ1xcHA4ePAgDAwNMnDgRn332Ga5fv97l9e3LdwwePBimpqaSdY7WrVuHLVu2YNasWbC0tMSsWbMQEhKC3bt399m9kP6PitIAMXv2bBQVFeGnn36Cl5cXkpKS4OzsjLi4OJl+f11dHfLy8uDn5ydpeenq6iIiIgJ5eXm9mzwZUGg9pQHM398f58+fx/379yEQCHDy5EnMnDkT//73v2FpaYns7GyMHTsWAFBaWgpTU1P87//+L1xdXaXiqKqq0oqMRGGoT2kAs7Oz63JsUnsfklgslpwzMTGBmZkZ7t27Bx8fn75KkQxAVJQGgEePHuGdd97BwoUL4eDgAD09PWRmZmLz5s1d7mtmbGwMLS0tnD17FsOGDYOmpiZEIhHWr1+PJUuWQCQSwcvLC42NjcjMzERlZSVCQ0OVcGekX1L26z/S+xoaGpiVK1cyzs7OjEgkYrS1tRkbGxtm9erVTH19PcMwTKdtqffu3cuYm5t32lbo8OHDzNixYxmhUMgMGjSImTx5MnPixIk+viPSn1GfEiGEU+jtGyGEU6goEUI4hYoSIYRTqCgRQjiFihIhhFOoKBFCOIWKEiGEU6goEUI4hYoSIYRTqCgRQjiFihIhhFP+D+ii8ziK8xJKAAAAAElFTkSuQmCC",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 1, figsize=(3, 3))\n",
"im = axes.imshow(\n",
" sites.T.iloc[::-1, :],\n",
" cmap=\"Greys\",\n",
" vmin=0,\n",
")\n",
"axes.set(\n",
" xticks=[0, 1, 2, 3],\n",
" xticklabels=[39, 40, 41, 54],\n",
" xlabel=\"Site\",\n",
" yticks=[0, 1, 2, 3],\n",
" yticklabels=[1, 2, 3, 4][::-1],\n",
" ylabel=\"Interaction order $k$\",\n",
" aspect='equal'\n",
")\n",
"plt.colorbar(im, shrink=0.6, label='% variance explained')\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "f6f038e5",
"metadata": {},
"source": [
"This clearly shows that some sites (41 and 54) have a much larger contribution to the phenotypic variance than others (40), and that they contribute mostly through their additive, pairwise and even threeway interactions.\n",
"\n",
"### How to compute the variance explained by interactions between every pair of sites\n",
"\n",
"Another interesting way to aggregate the variance explain by every subset of sites $U$ is by computing the variance explain by interactions of any order involving a pair of sites. \n"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "9e2ad948",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" site1 \n",
" site2 \n",
" variance \n",
" variance_perc \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0 \n",
" 1 \n",
" 12442.433333 \n",
" 7.733298 \n",
" \n",
" \n",
" 1 \n",
" 0 \n",
" 2 \n",
" 49192.407027 \n",
" 30.574370 \n",
" \n",
" \n",
" 2 \n",
" 0 \n",
" 3 \n",
" 40050.839230 \n",
" 24.892646 \n",
" \n",
" \n",
" 3 \n",
" 1 \n",
" 2 \n",
" 19654.479720 \n",
" 12.215774 \n",
" \n",
" \n",
" 4 \n",
" 1 \n",
" 3 \n",
" 15268.134797 \n",
" 9.489546 \n",
" \n",
" \n",
" 5 \n",
" 2 \n",
" 3 \n",
" 91714.303247 \n",
" 57.002842 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" site1 site2 variance variance_perc\n",
"0 0 1 12442.433333 7.733298\n",
"1 0 2 49192.407027 30.574370\n",
"2 0 3 40050.839230 24.892646\n",
"3 1 2 19654.479720 12.215774\n",
"4 1 3 15268.134797 9.489546\n",
"5 2 3 91714.303247 57.002842"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pairs = gb1s.calc_site_pairs_variance_perc(v_U_vcs)\n",
"pairs\n"
]
},
{
"cell_type": "markdown",
"id": "a7d1f007",
"metadata": {},
"source": [
"We can arrange these values on a matrix and represent them in a heatmap"
]
},
{
"cell_type": "code",
"execution_count": 54,
"id": "3e1cdd04",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" site2 \n",
" 1 \n",
" 2 \n",
" 3 \n",
" \n",
" \n",
" site1 \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 7.733298 \n",
" 30.574370 \n",
" 24.892646 \n",
" \n",
" \n",
" 1 \n",
" NaN \n",
" 12.215774 \n",
" 9.489546 \n",
" \n",
" \n",
" 2 \n",
" NaN \n",
" NaN \n",
" 57.002842 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
"site2 1 2 3\n",
"site1 \n",
"0 7.733298 30.574370 24.892646\n",
"1 NaN 12.215774 9.489546\n",
"2 NaN NaN 57.002842"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = pd.pivot_table(\n",
" pairs, index=\"site1\", columns=\"site2\", values=\"variance_perc\"\n",
")\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "4334860c",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAADiCAYAAADjymrMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAmHklEQVR4nO3deVxU1f8/8NcgDOvMsMgiCYhLjAsoYCKp4EfNJaNUHlbI56uhthghgpqiJrnCwxZXckET7aGhlKXmlkGiEiibimkURILBgMoiiywO9/eHP+bjOKDDZWDuHd/Px+M+Hsy545tzanx7zplzzxEwDMOAEEI4RE/bFSCEkCdRYiKEcA4lJkII51BiIoRwDiUmQgjnUGIihHAOJSZCCOfoa7sCna25uRnFxcUQiUQQCATarg55jjEMg+rqatjb20NPj/oET6Pziam4uBgODg7argYhCkVFRejZs6e2q8FpOp+YRCIRACAnJ0fxsy755ptvtF2FTqGL/5g8ePAA8+bN08nPoabpfGJqGb6JRCKIxWIt10bzjIyMtF2FTmFiYqLtKnQamlJ4NhroEkI4hxITIYRzKDERQjhH5+eYCOGr+vp6NDY2tnpPKBTq7PwiQImJEE6qr6+HsbFxm/ft7OxQUFCgs8mJhnKEcFBLT0kgELR6yWSyNntTuoB6TIRwWEsiehzDMND1jWcpMRHCYd26dWs1MTU3N2upRo9s2bJF7ffOnz+/3fEpMRHCYXp6eq0mJm3buHGj0us7d+6grq4O5ubmAIDKykqYmJjAxsaGVWKiOSZCOKytOSZtKygoUFzr1q3DkCFDcPPmTZSXl6O8vBw3b96Eh4cH1qxZwyo+JSZCOKxbt26tXlzyySefYOvWrXBxcVGUubi4YOPGjVixYgWrmDSUI4TDuDqUe1xJSQkePnyoUi6Xy1FaWsoqJvWYCOEwPT29Vi8uGTt2LN5//31kZWUpyjIzMzFv3jyMGzeOVUxutZAQooSrc0yP+/rrr2FnZ4ehQ4fC0NAQhoaGGDZsGGxtbbF7925WMWkoRwiHtdZD0vZSgSdZW1vj5MmT+PPPP/HHH38AAKRSKV588UXWMSkxEcJhXBy6taVXr15gGAZ9+vSBvn7HUgs/WkzIc4oPQ7m6ujrMmTMHJiYmGDhwIAoLCwEAISEhiI6OZhWTEhMhHMaH5QIRERG4evUqzp07p/RQ8bhx43Do0CFWMWkoRwiH8WEo9+OPP+LQoUMYPny4Um9u4MCByM/PZxWTEhMhHMaHxHTnzh3Y2NiolNfW1rIednK7xYQ85/gwxzR06FCcOHFC8bqlfrt374a3tzermFpNTNu3b4ebmxvEYjHEYjG8vb1x6tQpxf38/HxMnToV1tbWEIvFePPNN1mvJCWEj/gwx7R+/XosW7YM8+bNw8OHD7F582aMHz8ee/fuxbp161jF1Gpi6tmzJ6Kjo5GZmYmMjAyMGTMGb7zxBn7//XfU1tZi/PjxEAgESEpKQkpKChobG+Hn58e5dRyEdBY+rPweOXIkrly5gocPH8LV1RU///wzbGxskJqaCk9PT1YxtTrH5Ofnp/R63bp12L59O9LS0vDvv//in3/+QXZ2tuI8uH379sHCwgJJSUmsl7oTwjdcG7q1pk+fPoiNjdVYPM5MfsvlciQkJKC2thbe3t7Iz8+HQCCAoaGh4j1GRkbQ09PDxYsX20xMDQ0NaGhoULy+f/9+p9edkM7CxaFba5qbm5GXl4eysjKVEY2Pj0+742k9MeXk5MDb2xv19fUwMzPDDz/8gAEDBsDa2hqmpqZYsmQJ1q9fD4ZhsHTpUsjlcpSUlLQZLyoqCqtWrerCFhDSefiQmNLS0jBjxgzcunVLZecDgUAAuVze7phaH6y6uLjgypUruHTpEubNm4dZs2bhxo0bsLa2RkJCAo4fPw4zMzNIJBJUVlbCw8PjqWPsiIgIVFVVKa6ioqIubA0hmsWHye8PPvgAQ4cOxfXr11FeXo6KigrFVV5eziqm1ntMQqEQffv2BQB4enoiPT0dmzdvxs6dOzF+/Hjk5+fj7t270NfXh7m5Oezs7NC7d+8247U83UyILmhteQDX5pz++usvfPfdd4q/x5qg9R7Tk5qbm5XmiACge/fuMDc3R1JSEsrKyvD6669rqXaEdC0+9Ji8vLyQl5en0Zha7TFFRERg0qRJcHR0RHV1NQ4ePIhz587hzJkzAIC9e/eif//+sLa2RmpqKkJDQxEWFqa0hSchuoyLiehJISEhWLhwIWQyGVxdXWFgYKB0383Nrd0xtZqYysrKMHPmTJSUlEAikcDNzQ1nzpzBK6+8AgDIzc1FREQEysvL0atXLyxfvhxhYWHarDIhXYoPQzl/f38AwOzZsxVlAoEADMOwnvzWamLas2fPU+9HR0ez3jaBEF2gp6en0mPi2p7fBQUFGo+p9clvQkjbWhvKcS0xOTk5aTwmJSZCOIyrienYsWOYNGkSDAwMcOzYsae+l82XVZSYCOEwrs4xTZkyBTKZDDY2NpgyZUqb7+PlHBMh5Om42mN6/LGTznionhITIRzG1cTU2SgxEcJhXB3KPam2thbJyckoLCxEY2Oj0r358+e3Ox4lJkI4rLXlAlzbjyw7Oxuvvvoq6urqUFtbC0tLS9y9excmJiawsbFhlZg490gKIeR/+PBISlhYGPz8/FBRUQFjY2OkpaXh1q1b8PT0xOeff84qJiUmQjiMD4npypUrWLhwoaJ319DQAAcHB2zYsAHLli1jFZMSEyEcxofDCAwMDBRbEdnY2CgOvJRIJKy3HaI5JkI4rLU9vrm257e7uzvS09PRr18/+Pr6YuXKlbh79y6++eYbDBo0iFVMbrWQEKKED4cRrF+/Hj169ADwaN9+CwsLzJs3D3fu3MGuXbtYxaQeEyEcxoflAkOHDlX8bGNjg9OnT3c4JrdSLyFEiaZ7TNHR0RAIBFiwYIGirL6+HsHBwbCysoKZmRn8/f21fn4j9ZgI4bDWvoVju44pPT0dO3fuVNm4LSwsDCdOnEBCQgIkEgk++ugjTJs2DSkpKW3Gcnd3V7vnlpWV1e66UmIihMMEAoFKD4nNUK6mpgaBgYGIjY3F2rVrFeVVVVXYs2cPDh48iDFjxgD4386xaWlpGD58eKvxnvbgriZQYiKEw542x/TkmYlPO4gjODgYkydPxrhx45QSU2ZmJpqampTOaZRKpXB0dERqamqbiSkyMpJVe9RFiYkQDnvacgEHBwel8sjISHz66acqMeLj45GVlYX09HSVezKZDEKhEObm5krltra2kMlk7aprRkYGbt68CQAYMGAA6+PBAUpMhHDa0xJTUVERxGKxory13lJRURFCQ0Nx9uxZGBkZdUodb9++jYCAAKSkpCgSXGVlJV5++WXEx8ejZ8+e7Y5J38oRwmFPW/ktFouVrtYSU2ZmJsrKyuDh4QF9fX3o6+sjOTkZW7Zsgb6+PmxtbdHY2IjKykqlP1daWgo7Ozu16jh37lw0NTXh5s2bKC8vR3l5OW7evInm5mbMnTuXVbufmx6ThYWF0r8uumL69OnarkKn6Kx/3bWpurq63X+moyu/x44di5ycHKWyoKAgSKVSLFmyBA4ODjAwMEBiYqLitJPc3FwUFhbC29tbrd+RnJyM3377TelYNRcXF2zduhWjRo1Su66Pe24SEyF81NHEJBKJVB4LMTU1hZWVlaJ8zpw5CA8Ph6WlJcRiMUJCQuDt7d3mxPeTHBwc0NTUpFIul8thb2+vdl0fR0M5QjisKx5J2bhxI1577TX4+/vDx8cHdnZ2OHLkiNp//rPPPkNISAgyMjIUZRkZGQgNDWW97YmA0fF9Ou/fvw+JRIKqqiqdHMrdunVL21XoFLo6lOvXr59an8WWz+3cuXMhFAqV7jU2NmL37t2c+UxbWFigrq4ODx8+hL7+o0FYy8+mpqZK7y0vL1crJg3lCOEwPuwusGnTJo3HpMRECIdpauV3Z5o1a5bGY3Ir9RJClPBho7i4uLhWyx8+fIiIiAhWMSkxEcJhfNiPaf78+Zg+fToqKioUZbm5ufDy8sK3337LKia3WkgIUcKHxJSdnY3bt2/D1dUVZ8+eRUxMDDw8PCCVSnH16lVWMWmOiRAOa23bE64dRtCnTx+kpKRgwYIFmDhxIrp164Z9+/YhICCAdcx2pd4HDx7g4sWLuHHjhsq9+vp67N+/n3VFCCGq+DDHBAAnTpxAfHw8vL29YW5ujj179qC4uJh1PLUT059//on+/fvDx8cHrq6u8PX1RUlJieJ+VVUVgoKCWFeEEKKKD0O5999/H9OnT8eSJUtw4cIFXLt2DUKhEK6urjh8+DCrmGq3cMmSJRg0aBDKysqQm5sLkUiEESNGKI5qIYRoHh8SU0pKCi5duoSFCxdCIBDAzs4OJ0+exOrVqzF79mxWMdWeY/rtt9/wyy+/oHv37ujevTuOHz+ODz/8EKNGjcKvv/6qssKTENJxfDiMIDMzs9WdDYKDg5U2oGsPtVPvgwcPFMvNgUf/cbZv3w4/Pz/4+vrizz//ZFUBQkjb+NBjMjQ0RH5+PlasWIGAgACUlZUBAE6dOoWHDx+yiql2C6VSqdJDei22bduGN954A6+//jqrChBC2tZy7PbjF9cSU3JyMlxdXXHp0iUcOXIENTU1AICrV6+y3oJX7RZOnTq1zcVS27ZtQ0BAAHT8eWBCuhwfvpVbunQp1q5di7Nnzyo9cDxmzBikpaWxiql2YoqIiMDJkyfbvP/VV1+xPlaGENI6PgzlcnJyMHXqVJVyGxsb3L17l1VMbrWQEKKED4nJ3NxcaelQi+zsbLzwwgusYnKrhYQQJXxITG+//TaWLFkCmUwGgUCA5uZmpKSkYNGiRZg5cyarmNxqISFECR/mmNavXw+pVAoHBwfU1NRgwIAB8PHxwcsvv4wVK1awiknPyhHCYXx4Vk4oFCI2NhYrV65ETk4Oampq4O7ujn79+rGOSYmJEA7jww6WLRwcHFQO4WSLVQu/+eYbjBgxAvb29oo9pzdt2oSjR4+yrkh0dDQEAgEWLFigKKuvr0dwcDCsrKxgZmYGf39/lJaWsv4dhPANH4ZynaHdiWn79u0IDw/Hq6++isrKSsjlcgCPZubZ7v2bnp6OnTt3ws3NTak8LCwMx48fR0JCApKTk1FcXIxp06ax+h2E8BEfJr87Q7tbuHXrVsTGxmL58uVKY92hQ4eqHKynjpqaGgQGBiI2NhYWFhaK8qqqKuzZswdffvklxowZA09PT+zduxe//fYb60VbhPBNy57fj1/UY2pFQUEB3N3dVcoNDQ1RW1vb7goEBwdj8uTJKg/7ZWZmoqmpSalcKpXC0dERqampbcZraGjA/fv3lS5C+Ip6TGpydnbGlStXVMpPnz6N/v37tytWfHw8srKyEBUVpXJPJpNBKBTC3NxcqdzW1hYymazNmFFRUZBIJIpLU5NxhGgDXxLThQsX8N///hfe3t74999/ATyai7548SKreO1uYXh4OIKDg3Ho0CEwDIPLly9j3bp1iIiIwMcff6x2nKKiIoSGhuLAgQMaPdwwIiICVVVViquoqEhjsQnpanxITN9//z0mTJgAY2NjZGdno6GhAcCj6Zj169ezitnu5QJz586FsbExVqxYgbq6OsyYMQP29vbYvHkz3n77bbXjZGZmoqysDB4eHooyuVyO8+fPY9u2bThz5gwaGxtRWVmp1GsqLS2FnZ1dm3ENDQ1b3RuGED7iwzqmtWvXYseOHZg5cybi4+MV5SNGjMDatWtZxWS1jikwMBCBgYGoq6tDTU0NbGxs2h1j7NixKpPlQUFBkEqlWLJkCRwcHGBgYIDExET4+/sDeHQkTGFhIby9vdlUmxDe4cNGcbm5ufDx8VEpl0gkqKysZBWz3YlpzJgxOHLkCMzNzWFiYgITExMAj85anzJlCpKSktSKIxKJMGjQIKUyU1NTWFlZKcrnzJmD8PBwWFpaQiwWIyQkBN7e3hg+fHh7q00IL/FhgaWdnR3y8vLQq1cvpfKLFy+id+/erGK2OzGdO3cOjY2NKuX19fW4cOECq0q0ZePGjdDT04O/vz8aGhowYcIEfPXVVxr9HYRwGR+Gcu+++y5CQ0Px9ddfQyAQoLi4GKmpqVi0aBE++eQTVjHVTkzXrl1T/Hzjxg2lb8bkcjlOnz7NeouDFufOnVN6bWRkhJiYGMTExHQoLiF81bKO6ckyLlm6dCmam5sxduxY1NXVwcfHB4aGhli0aBFCQkJYxVQ7MQ0ZMkQx3h0zZozKfWNjY2zdupVVJQghrePDHJNAIMDy5cuxePFi5OXlKXYYMDMzYx1T7cRUUFAAhmHQu3dvXL58GdbW1op7QqEQNjY2nOtiEsJ3fBjKVVVVQS6Xw9LSEgMGDFCUl5eXQ19fH2KxuN0x1U5MTk5OAEDb5xLShfgw+f3222/Dz88PH374oVL54cOHcezYsaduyd0WtRLTsWPHMGnSJBgYGODYsWNPfS+dlkKI5vBhKHfp0iV8+eWXKuWjR4/G8uXLWcVUKzFNmTIFMpkMNjY2mDJlSpvvEwgEit0GCCEdx4ehXENDQ6vnxzU1NeHBgwesYqrVJ2xublYsomxubm7zoqREiGbx4ZGUYcOGYdeuXSrlO3bsgKenJ6uYtIMlIRzGhx7T2rVrMW7cOFy9ehVjx44FACQmJiI9PR0///wzq5hqp97U1FT89NNPSmX79++Hs7MzbGxs8N577yke3iOEaAYfdrAcMWIEUlNT4eDggMOHD+P48ePo27cvrl27hlGjRrGKqXaPafXq1Rg9ejRee+01AI8OuZszZw7eeecd9O/fH5999hns7e3x6aefsqoIIUQVH76VAx6tczxw4IDG4qmdmK5cuYI1a9YoXsfHx8PLywuxsbEAHm1EHhkZSYmJEA3S09NTGbpxMTE1NzcjLy8PZWVlKkuKWnvA91nUTkwVFRWwtbVVvE5OTsakSZMUr1966SXa+4gQDePDcoG0tDTMmDEDt27dAsMwSvfYflOvduq1tbVFQUEBAKCxsRFZWVlKT/lXV1fDwMCg3RUghLStZfL7yYtLPvjgAwwdOhTXr19HeXk5KioqFFd5eTmrmGonpldffRVLly7FhQsXEBERARMTE6WJrWvXrqFPnz6sKkEIaZ0mJr+joqLw0ksvQSQSKdYi5ubmKr2nI0el/fXXX1i/fj369+8Pc3Nzpa2tJRJJu+raQu3EtGbNGujr68PX1xexsbGIjY2FUChU3P/6668xfvx4VpUghLROE+uYkpOTERwcjLS0NJw9exZNTU0YP3680uEhHTkqzcvLC3l5ee2q07OoPcfUvXt3nD9/HlVVVTAzM1PpTiYkJHToaWJCiKqnzTE9eQJQW9tKnz59Wul1XFwcbGxskJmZCR8fH8VRaQcPHlTsHLJ37170798faWlpz9yYMSQkBAsXLoRMJoOrq6vKlM6T50Wqo90LLNvqmllaWrb7lxNCnu5pywWePAFI3W/Fq6qqAPzv7+yzjkp7VmJq2fp69uzZijKBQACGYVhPftPKb0I47GmJqaioSGlLEXUO4WhubsaCBQswYsQIxRbWbI9Ka9HypZgmUWIihKfEYnG79zoKDg7G9evXWZ/31pqWLZE0iRITIRymyZXfH330EX766SecP38ePXv2VJTb2dmxOirtSTdu3EBhYaHKmQBstkKixEQIh2liz2+GYRASEoIffvgB586dg7Ozs9J9T0/PDh2V9vfff2Pq1KnIyclRzC09Xk+aY3oOdUY3mgu4trpZWzTRYwoODsbBgwdx9OhRiEQixbyRRCKBsbExJBJJh45KCw0NhbOzMxITE+Hs7IzLly/j3r17WLhwIT7//PN21bUFJSZCOEwTj6Rs374dwKMdJR+3d+9evPPOOwA6dlRaamoqkpKS0L17d0UiHTlyJKKiojB//nxkZ2e3q74AJSZCOE0TienJ59da05Gj0uRyOUQiEYBH6x2Li4vh4uICJycnlRXm6qLERAiH8WHbk0GDBuHq1atwdnaGl5cXNmzYAKFQiF27dnXdSbyEkK7Dh90FVqxYoXi8ZfXq1XjttdcwatQoWFlZ4dChQ6xiUmIihMP40GOaMGGC4ue+ffvijz/+QHl5OSwsLFgnUUpMhHAYHxJTazr6iBolJkI4jKtDuWnTpiEuLg5isfiZuxAcOXKk3fEpMRHCYVxNTBKJRFEPtnsuPQ0lJkI4jKtDub179wJ4tBRh1apVsLa2hrGxscbia7+FhJA2cf3AS4Zh0LdvX9y+fVujcbnTQkII7+jp6aFfv364d++eZuNqNBohRKNaHuJ9/OLCHNPjoqOjsXjxYly/fl1jMWmOiRAO4+rk9+NmzpyJuro6DB48GEKhUGWuic1JKZSYCOEwrk5+P27Tpk0aj0mJiRDSIbNmzdJ4TEpMhHAYH3pMj6uvr1fZwbK92/8CNPlNCKdp4sDLzlZbW4uPPvoINjY2MDU1hYWFhdLFBiUmQjiMD4np448/RlJSErZv3w5DQ0Ps3r0bq1atgr29Pfbv388qJg3lCOEwPnwrd/z4cezfvx+jR49GUFAQRo0ahb59+8LJyQkHDhxAYGBgu2NSj4kQDuP6ym/g0XKAlg3hxGKxYnnAyJEjcf78eVYxudVCQogSPgzlevfurTj0UiqV4vDhwwAe9aSePERTXZSYCOEwPiSmoKAgXL16FQCwdOlSxMTEwMjICGFhYVi8eDGrmDTHRAiH8WGOKSwsTPHzuHHj8McffyAzMxN9+/aFm5sbq5ic6TFFR0dDIBBgwYIFirJdu3Zh9OjREIvFEAgEqKys1Fr9CNEGPswxFRUVKb12cnLCtGnTWCclgCOJKT09HTt37lRpSF1dHSZOnIhly5ZpqWaEaBcfhnK9evWCr68vYmNjUVFRoZGYWk9MNTU1CAwMRGxsrMpirAULFmDp0qVqnQbaoqGhAffv31e6COEzLiclAMjIyMCwYcOwevVq9OjRA1OmTMF3332HhoYG1jG1npiCg4MxefJkjBs3TiPxoqKiIJFIFJeDg4NG4hKiDXzoMbm7u+Ozzz5DYWEhTp06BWtra7z33nuwtbXF7NmzWcXUamKKj49HVlYWoqKiNBYzIiICVVVViuvJ8S8hfMKHxNRCIBDgP//5D2JjY/HLL7/A2dkZ+/btYxVLa9/KFRUVITQ0FGfPnoWRkZHG4hoaGsLQ0FBj8QjRJj58K9fi9u3bOHjwIA4ePIjr16/D29ub1ZHjgBYTU2ZmJsrKyuDh4aEok8vlOH/+PLZt24aGhgZ069ZNW9UjhBP4kJh27tyJgwcPIiUlBVKpFIGBgTh69CicnJxYx9RaYho7dixycnKUyoKCgiCVSrFkyRJKSoSAH4lp7dq1CAgIwJYtWzB48GCNxNRaYhKJRBg0aJBSmampKaysrBTlMpkMMpkMeXl5AICcnByIRCI4Ojp2+KRPQviAD4mpsLBQ43XS+rdyT7Njxw64u7vj3XffBQD4+PjA3d0dx44d03LNCOkafJj87oz6CBiGYTQelUPu378PiUSCqqoqVjvpEe3g2l8+TVLns9jyuS0pKVF57/3799GjRw+d/kxzusdECHk+0UO8hHAYH+aYOgMlJkI4jG+J6e7du7h06RLkcjleeukl9OjRg1UcSkyEcBifEtP333+POXPm4MUXX0RTUxNyc3MRExODoKCgdseiOSZCCCs1NTVKr1etWoXLly/j8uXLyM7ORkJCApYvX84qNiUmQjiMy8sFPD09cfToUcVrfX19lJWVKV6XlpZCKBSyik1DOUI4jMtDuTNnziA4OBhxcXGIiYnB5s2b8dZbb0Eul+Phw4fQ09NDXFwcq9iUmAjhMC4npl69euHEiRP49ttv4evri/nz5yMvLw95eXmQy+WQSqWsH9CnoRwhHMbloVyLgIAApKen4+rVqxg9ejSam5sxZMiQDu0aQomJEA7TZGKKiYlBr169YGRkBC8vL1y+fLnD9Tt58iS++OILZGRkYPfu3diwYQMCAwOxePFiPHjwgHVcSkyEcJimEtOhQ4cQHh6OyMhIZGVlYfDgwZgwYYLSZHV7LVy4EEFBQUhPT8f777+PNWvWwNfXF1lZWTAyMoK7uztOnTrFLjij46qqqhgATFVVlbarQtoBgM5e6nwWWz63lZWVTHNzs9JVWVnZ7s/0sGHDmODgYMVruVzO2NvbM1FRUaz+/zAMw1haWjIZGRkMwzDMvXv3mH79+ind//3335mRI0eyik09JkI4rLq6utULgMqhG21t/t/Y2IjMzEylffX19PQwbtw4pKamsq6bqamp4gTeoqIilTmlAQMG4MKFC6xiU2IihIOEQiHs7Ozg4OCgdLhGywEbZmZmKvfa2jv/7t27kMvlsLW1VSq3tbWFTCZjXceoqCjMnDkT9vb28PX1xZo1a1jHehItFyCEg4yMjFBQUIDGxsZW7zMMozLX1NV73QcGBmLixIn4+++/0a9fP5ibm2ssNiUmQjjKyMhIIwd1dO/eHd26dUNpaalSeWlpKezs7DoU28rKClZWVh2K0RoayhGi44RCITw9PZGYmKgoa25uRmJiIry9vbVYs7ZRj4mQ50B4eDhmzZqFoUOHYtiwYdi0aRNqa2tZPfnfFSgxEfIceOutt3Dnzh2sXLkSMpkMQ4YMwenTp1UmxLlC5/f8rqqqgrm5OYqKinR2f2RdJJFItF2FTlNZWanT7dMEne8xtaz5cHBw0HJNCHmkurqaEtMz6HyPqbm5GcXFxRCJRJ3+8OP9+/fh4OCgk70zXW1bV7aLYRhUV1fD3t4eenr0vdPT6HyPSU9PDz179uzS3ykWi3XqL+/jdLVtXdUu6imph9I2IYRzKDERQjiHEpMGGRoaIjIysssfDegKuto2XW0X3+n85DchhH+ox0QI4RxKTIQQzqHERAjhHEpMhBDOocTUQdHR0RAIBFiwYIGirL6+HsHBwbCysoKZmRn8/f1V9sLhutbatWvXLowePRpisRgCgQCVlZVaq197fPrppyqb+UulUpX3MQyDSZMmQSAQ4Mcff+z6ihIFSkwdkJ6ejp07d8LNzU2pPCwsDMePH0dCQgKSk5NRXFyMadOmaamW7ddWu+rq6jBx4kQsW7ZMSzVjb+DAgSgpKVFcFy9eVHnPpk2bOHdm23OL7QkJz7vq6mqmX79+zNmzZxlfX18mNDSUYRiGqaysZAwMDJiEhATFe2/evMkAYFJTU7VUW/W11a7H/frrrwwApqKiosvrx0ZkZCQzePDgp74nOzubeeGFF5iSkhIGAPPDDz90Sd1I66jHxFJwcDAmT56sdPIEAGRmZqKpqUmpXCqVwtHRsUMnUnSVttrFd3/99Rfs7e3Ru3dvBAYGorCwUHGvrq4OM2bMQExMTIe3miWaofMP8XaG+Ph4ZGVlIT09XeWeTCaDUChU2Zi9oydSdIWntYvPvLy8EBcXBxcXF5SUlGDVqlUYNWoUrl+/DpFIhLCwMLz88st44403tF1V8v9RYmqnoqIihIaG4uzZsxrZKJ4rdLVdADBp0iTFz25ubvDy8oKTkxMOHz4Ma2trJCUlITs7W4s1JE+ioVw7ZWZmoqysDB4eHtDX14e+vj6Sk5OxZcsW6Ovrw9bWFo2NjSrfWGniRIrO9Kx2yeVybVdRY8zNzfHiiy8iLy8PSUlJyM/Ph7m5uaLdAODv74/Ro0drt6LPMeoxtdPYsWORk5OjVBYUFASpVIolS5bAwcEBBgYGSExMhL+/PwAgNzcXhYWFnD2RAnh2u7p166almmleTU0N8vPz8X//93948803MXfuXKX7rq6u2LhxI/z8/LRUQ0KJqZ1EIhEGDRqkVGZqagorKytF+Zw5cxAeHg5LS0uIxWKEhITA29sbw4cP10aV1aJOu2QyGWQyGfLy8gAAOTk5EIlEcHR0hKWlZZfXWV2LFi2Cn58fnJycUFxcjMjISHTr1g0BAQGwtrZutSfr6OgIZ2dnLdSWAJSYOsXGjRuhp6cHf39/NDQ0YMKECfjqq6+0Xa0O27FjB1atWqV47ePjAwDYu3cv3nnnHS3V6tlu376NgIAA3Lt3D9bW1hg5ciTS0tJgbW2t7aqRNtC2J4QQzqHJb0II51BiIoRwDiUmQgjnUGIihHAOJSZCCOdQYiKEcA4lJkII51BiIoRwDiUmQlvJEs6hxKTj7ty5g3nz5sHR0RGGhoaws7PDhAkTkJKSonhPSUmJYmuQf/75BwKBAFeuXOnw7z5//jz8/Pxgb29PyY+0Cz0rp+P8/f3R2NiIffv2oXfv3igtLUViYiLu3buneE9nbcdSW1uLwYMHY/bs2bza85xwgLb39iWdp6KiggHAnDt37qnvw2N7XANQunx9fRXvi42NZaRSKWNoaMi4uLgwMTExatcFtI82aQfqMekwMzMzmJmZ4ccff8Tw4cNhaGj4zD9z+fJlDBs2DL/88gsGDhwIoVAIADhw4ABWrlyJbdu2wd3dHdnZ2Xj33XdhamqKWbNmdXZTyHOG5ph0mL6+PuLi4rBv3z6Ym5tjxIgRWLZsGa5du9bmn2nZCsTKygp2dnaKfZYiIyPxxRdfYNq0aXB2dsa0adMQFhaGnTt3dklbyPOFEpOO8/f3R3FxMY4dO4aJEyfi3Llz8PDwQFxcnNoxamtrkZ+fjzlz5ih6YWZmZli7di3y8/M7r/LkuUVDueeAkZERXnnlFbzyyiv45JNPMHfuXERGRqq9uVtNTQ0AIDY2Fl5eXkr3dGnLXcIdlJieQwMGDGjzq/uWOaXHDx+wtbWFvb09/v77bwQGBnZFFclzjhKTDrt37x6mT5+O2bNnw83NDSKRCBkZGdiwYUObZ6jZ2NjA2NgYp0+fRs+ePWFkZASJRIJVq1Zh/vz5kEgkmDhxIhoaGpCRkYGKigqEh4e3GqumpkaxPzgAFBQU4MqVK7C0tISjo2OntJnoCG1/LUg6T319PbN06VLGw8ODkUgkjImJCePi4sKsWLGCqaurU7wPT3yVHxsbyzg4ODB6enpKywUOHDjADBkyhBEKhYyFhQXj4+PDHDlypM3f33KU+JPXrFmzOqG1RJfQnt+EEM6hb+UIIZxDiYkQwjmUmAghnEOJiRDCOZSYCCGcQ4mJEMI5lJgIIZxDiYkQwjmUmAghnEOJiRDCOZSYCCGc8/8AYK8D/WA//uQAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 1, figsize=(3, 3))\n",
"im = axes.imshow(\n",
" m,\n",
" cmap=\"Greys\",\n",
" vmin=0,\n",
")\n",
"axes.set(\n",
" xticks=[0, 1, 2],\n",
" xlabel=\"Site 1\",\n",
" yticks=[0, 1, 2],\n",
" xticklabels=[40, 41, 54],\n",
" yticklabels=[39, 40, 41],\n",
" ylabel=\"Site 2\",\n",
" aspect=\"equal\",\n",
")\n",
"plt.colorbar(im, shrink=0.6, label=\"% variance explained\")\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "ca23a0e1",
"metadata": {},
"source": [
"This plot clearly shows that sites 41 and 54 interact strongly with each other, but also with position 39.\n",
"\n",
"## Shine-Dalgarno fitness landscape\n",
"\n",
"In the previous example, the total number of subsets of sites $U$ is still relatively small given the short sequence length. To better illustrate the usefulness of this approach, we will use a 9-nucleotide long complete combinatorial landscape inferred with `VCregression` from [previously published high throughput experimental data](https://genome.cshlp.org/content/30/5/711) in the *dmsC* gene context. The inferred landscape is also available as a `DataSet` for illustration purposes.\n"
]
},
{
"cell_type": "code",
"execution_count": 58,
"id": "998509f8",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" f \n",
" \n",
" \n",
" \n",
" \n",
" AAAAAAAAA \n",
" 0.563857 \n",
" \n",
" \n",
" AAAAAAAAC \n",
" 0.701511 \n",
" \n",
" \n",
" AAAAAAAAG \n",
" 0.619188 \n",
" \n",
" \n",
" AAAAAAAAU \n",
" 0.628772 \n",
" \n",
" \n",
" AAAAAAACA \n",
" 0.594683 \n",
" \n",
" \n",
" ... \n",
" ... \n",
" \n",
" \n",
" UUUUUUUGU \n",
" 0.615910 \n",
" \n",
" \n",
" UUUUUUUUA \n",
" 0.539058 \n",
" \n",
" \n",
" UUUUUUUUC \n",
" 0.548782 \n",
" \n",
" \n",
" UUUUUUUUG \n",
" 0.536610 \n",
" \n",
" \n",
" UUUUUUUUU \n",
" 0.564960 \n",
" \n",
" \n",
"
\n",
"
262144 rows × 1 columns
\n",
"
"
],
"text/plain": [
" f\n",
"AAAAAAAAA 0.563857\n",
"AAAAAAAAC 0.701511\n",
"AAAAAAAAG 0.619188\n",
"AAAAAAAAU 0.628772\n",
"AAAAAAACA 0.594683\n",
"... ...\n",
"UUUUUUUGU 0.615910\n",
"UUUUUUUUA 0.539058\n",
"UUUUUUUUC 0.548782\n",
"UUUUUUUUG 0.536610\n",
"UUUUUUUUU 0.564960\n",
"\n",
"[262144 rows x 1 columns]"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sd = DataSet('dmsc')\n",
"sd.landscape"
]
},
{
"cell_type": "markdown",
"id": "4043bcbf",
"metadata": {},
"source": [
"Lets start again by defining the `GPmapSummarizer` object and compute the variance explained by every possible subset of sites"
]
},
{
"cell_type": "code",
"execution_count": 59,
"id": "cc1ce857",
"metadata": {},
"outputs": [],
"source": [
"sds = GPmapSummarizer(n_alleles=4, seq_length=9, f=sd.landscape[\"f\"].values)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "d1fcbd68",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" U \n",
" k \n",
" variance \n",
" variance_perc \n",
" variance_perc_cum \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" {0} \n",
" 1 \n",
" 1574.656812 \n",
" 5.500762 \n",
" 5.500762 \n",
" \n",
" \n",
" 1 \n",
" {1} \n",
" 1 \n",
" 2687.622462 \n",
" 9.388695 \n",
" 14.889458 \n",
" \n",
" \n",
" 2 \n",
" {2} \n",
" 1 \n",
" 2312.377111 \n",
" 8.077847 \n",
" 22.967305 \n",
" \n",
" \n",
" 3 \n",
" {3} \n",
" 1 \n",
" 1806.662119 \n",
" 6.311229 \n",
" 29.278533 \n",
" \n",
" \n",
" 4 \n",
" {4} \n",
" 1 \n",
" 1051.784760 \n",
" 3.674209 \n",
" 32.952742 \n",
" \n",
" \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" ... \n",
" \n",
" \n",
" 506 \n",
" {0, 1, 2, 4, 5, 6, 7, 8} \n",
" 8 \n",
" 2.485503 \n",
" 0.008683 \n",
" 99.956748 \n",
" \n",
" \n",
" 507 \n",
" {0, 1, 3, 4, 5, 6, 7, 8} \n",
" 8 \n",
" 2.545558 \n",
" 0.008892 \n",
" 99.965641 \n",
" \n",
" \n",
" 508 \n",
" {0, 2, 3, 4, 5, 6, 7, 8} \n",
" 8 \n",
" 2.358090 \n",
" 0.008238 \n",
" 99.973878 \n",
" \n",
" \n",
" 509 \n",
" {1, 2, 3, 4, 5, 6, 7, 8} \n",
" 8 \n",
" 2.796804 \n",
" 0.009770 \n",
" 99.983648 \n",
" \n",
" \n",
" 510 \n",
" {0, 1, 2, 3, 4, 5, 6, 7, 8} \n",
" 9 \n",
" 4.680847 \n",
" 0.016352 \n",
" 100.000000 \n",
" \n",
" \n",
"
\n",
"
511 rows × 5 columns
\n",
"
"
],
"text/plain": [
" U k variance variance_perc \\\n",
"0 {0} 1 1574.656812 5.500762 \n",
"1 {1} 1 2687.622462 9.388695 \n",
"2 {2} 1 2312.377111 8.077847 \n",
"3 {3} 1 1806.662119 6.311229 \n",
"4 {4} 1 1051.784760 3.674209 \n",
".. ... .. ... ... \n",
"506 {0, 1, 2, 4, 5, 6, 7, 8} 8 2.485503 0.008683 \n",
"507 {0, 1, 3, 4, 5, 6, 7, 8} 8 2.545558 0.008892 \n",
"508 {0, 2, 3, 4, 5, 6, 7, 8} 8 2.358090 0.008238 \n",
"509 {1, 2, 3, 4, 5, 6, 7, 8} 8 2.796804 0.009770 \n",
"510 {0, 1, 2, 3, 4, 5, 6, 7, 8} 9 4.680847 0.016352 \n",
"\n",
" variance_perc_cum \n",
"0 5.500762 \n",
"1 14.889458 \n",
"2 22.967305 \n",
"3 29.278533 \n",
"4 32.952742 \n",
".. ... \n",
"506 99.956748 \n",
"507 99.965641 \n",
"508 99.973878 \n",
"509 99.983648 \n",
"510 100.000000 \n",
"\n",
"[511 rows x 5 columns]"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v_U_vcs = sds.calc_V_U_variance_components()\n",
"v_U_vcs"
]
},
{
"cell_type": "markdown",
"id": "bda04e69",
"metadata": {},
"source": [
"We can see that the number of components now is much larger than before, so we really need to aggregate these components in a meaningful way to obtain low-dimensional summary statistics reflecting the patterns of genetic interactions in the inferred genotype-phenotype map.\n",
"\n",
"### How to compute the variance explained by interactions of order `k` involving each site"
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "9eaf5911",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" \n",
" 0 \n",
" 1 \n",
" 2 \n",
" 3 \n",
" 4 \n",
" 5 \n",
" 6 \n",
" 7 \n",
" 8 \n",
" \n",
" \n",
" \n",
" \n",
" 1 \n",
" 5.500762 \n",
" 9.388695 \n",
" 8.077847 \n",
" 6.311229 \n",
" 3.674209 \n",
" 2.277808 \n",
" 1.837005 \n",
" 0.835968 \n",
" 0.399951 \n",
" \n",
" \n",
" 2 \n",
" 6.390370 \n",
" 10.088863 \n",
" 10.262287 \n",
" 8.163453 \n",
" 6.898441 \n",
" 4.512101 \n",
" 3.034212 \n",
" 1.908863 \n",
" 1.153861 \n",
" \n",
" \n",
" 3 \n",
" 5.799007 \n",
" 8.710946 \n",
" 10.348244 \n",
" 11.095941 \n",
" 10.542178 \n",
" 8.887247 \n",
" 5.887798 \n",
" 3.522035 \n",
" 2.133580 \n",
" \n",
" \n",
" 4 \n",
" 2.999159 \n",
" 4.047676 \n",
" 4.776184 \n",
" 5.714224 \n",
" 5.930684 \n",
" 5.423643 \n",
" 4.125461 \n",
" 2.856694 \n",
" 1.932568 \n",
" \n",
" \n",
" 5 \n",
" 1.208228 \n",
" 1.508030 \n",
" 1.640607 \n",
" 1.763638 \n",
" 1.809694 \n",
" 1.735487 \n",
" 1.502939 \n",
" 1.196867 \n",
" 0.925827 \n",
" \n",
" \n",
" 6 \n",
" 0.441740 \n",
" 0.488117 \n",
" 0.516701 \n",
" 0.522051 \n",
" 0.525214 \n",
" 0.525015 \n",
" 0.496194 \n",
" 0.446754 \n",
" 0.389411 \n",
" \n",
" \n",
" 7 \n",
" 0.187551 \n",
" 0.194456 \n",
" 0.197175 \n",
" 0.198214 \n",
" 0.199869 \n",
" 0.199935 \n",
" 0.199599 \n",
" 0.192142 \n",
" 0.179375 \n",
" \n",
" \n",
" 8 \n",
" 0.070391 \n",
" 0.071923 \n",
" 0.071268 \n",
" 0.071478 \n",
" 0.071555 \n",
" 0.071555 \n",
" 0.071691 \n",
" 0.071658 \n",
" 0.069766 \n",
" \n",
" \n",
" 9 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" 0.016352 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
" 0 1 2 3 4 5 6 \\\n",
"1 5.500762 9.388695 8.077847 6.311229 3.674209 2.277808 1.837005 \n",
"2 6.390370 10.088863 10.262287 8.163453 6.898441 4.512101 3.034212 \n",
"3 5.799007 8.710946 10.348244 11.095941 10.542178 8.887247 5.887798 \n",
"4 2.999159 4.047676 4.776184 5.714224 5.930684 5.423643 4.125461 \n",
"5 1.208228 1.508030 1.640607 1.763638 1.809694 1.735487 1.502939 \n",
"6 0.441740 0.488117 0.516701 0.522051 0.525214 0.525015 0.496194 \n",
"7 0.187551 0.194456 0.197175 0.198214 0.199869 0.199935 0.199599 \n",
"8 0.070391 0.071923 0.071268 0.071478 0.071555 0.071555 0.071691 \n",
"9 0.016352 0.016352 0.016352 0.016352 0.016352 0.016352 0.016352 \n",
"\n",
" 7 8 \n",
"1 0.835968 0.399951 \n",
"2 1.908863 1.153861 \n",
"3 3.522035 2.133580 \n",
"4 2.856694 1.932568 \n",
"5 1.196867 0.925827 \n",
"6 0.446754 0.389411 \n",
"7 0.192142 0.179375 \n",
"8 0.071658 0.069766 \n",
"9 0.016352 0.016352 "
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sites = sds.calc_sites_variance_perc(v_U_vcs)\n",
"sites"
]
},
{
"cell_type": "code",
"execution_count": 65,
"id": "917f5350",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAFACAYAAABTBmBPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6pUlEQVR4nO3dd1RU1/o38O+IDEUBS0BEARELF4INS1Bs0YsQY4leY8GoWK6X2FuUxC6KGkvUGFGjookaNZaYa4sSxRK7YLkiihIBEVERRiBShvP+wcv5OQKROdMY+X7WOms5Zw7P3gd1ntl7n723TBAEAURERAAqGboCRERUfjApEBGRiEmBiIhETApERCRiUiAiIhGTAhERiZgUiIhIxKRARESiyoaugKYKCgqQnJwMKysryGQyQ1eHiPRMEAS8fPkSDg4OqFSJ33M1ZfRJITk5GY6OjoauBhEZWGJiIurWrWvoahg9o08KVlZWAAr/QVhbWxu4NkSkbwqFAo6OjuJnAWnG6JNCUZeRtbU1kwJRBcbuY+1gBxwREYmYFIiISMSkQEREIiYFIiISMSkQEZGISYGIiERMCkREJDJ4Unj58iUmTpwIZ2dnWFhYoG3btrh8+bKhq0VEVCEZPCmMHDkSx48fxw8//ICbN2/C19cXXbt2xaNHjwxdNSKiCkcmCIJgqML/+usvWFlZ4ZdffkH37t3F815eXvD390dISMhbYygUCtjY2CAjI4MzmokqIH4GaJdBl7nIz8+HUqmEubm5ynkLCwucPXu2xJ/JyclBTk6O+FqhUOi0jkREFYlBu4+srKzg7e2NBQsWIDk5GUqlEj/++CPOnz+Px48fl/gzoaGhsLGxEQ+ukEpEpD0G7T4CgPv372P48OE4ffo0TExM0KJFCzRq1AhXr15FTExMsetLaik4Ojqy6UhUQbH7SLsMvkqqq6srIiMjkZWVBYVCgdq1a6N///6oX79+idebmZnBzMxMz7UkIqoYDP70UZEqVaqgdu3aePHiBY4dO4ZevXoZukpERBWOwVsKx44dgyAIaNy4MeLi4jBt2jS4ubkhMDDQ0FUjIqpwDN5SyMjIwJgxY+Dm5oYhQ4bAx8cHx44dg6mpqaGrRkRU4Rh8oFlTHGQiqti08Rnw6tUr5ObmlulauVxe7DH6d4nBu4+IiAzp1atXsLCwKPP19vb2iI+Pf2cTg8G7j4iIDKmsLYQiKSkpav+MMWFLgYgIgEwmg0wm+9trBEGAkfe4vxWTAhERypYUADApEBFVBCYmJmVqKRQUFOipRobBpEBEhLK3FN51TApERGBSKPLOJIWKMABERMVp6/89k0KhdyYpEBFpgkmhEJMCERHKPtD8rmNSICICWwpFmBSIiMCkUIRJgYgITApFmBSIiMCkUIRJgYgIQKVKlVCp0t+vEfquz2YGmBSIiACUraVQEVoSTApERGBSKGLw/RSUSiVmzZoFFxcXWFhYwNXVFQsWLKgQzwMTUflRlBTedqjj9OnT6NGjBxwcHCCTyXDgwAGV9wVBwOzZs1G7dm1YWFiga9euuHfvnhbvSn0GTwpLlizBunXr8O233yImJgZLlizB0qVLsWbNGkNXjYgqEF0khaysLDRt2hRr164t8f2lS5di9erVCAsLw8WLF1GlShV069YNr1690sYtSWLw7qM//vgDvXr1Qvfu3QEA9erVw86dO3Hp0iUD14yIKhITE5O3DjSrmxT8/f3h7+9f4nuCIOCbb77BzJkz0atXLwDAtm3bUKtWLRw4cAADBgxQqyxtMXhLoW3btoiIiMDdu3cBANevX8fZs2dL/UUSEemCOi0FhUKhcuTk5KhdXnx8PFJSUtC1a1fxnI2NDdq0aYPz589r7b7UZfCWwowZM6BQKODm5gYTExMolUosXLgQAQEBJV6fk5Oj8hegUCj0VVUieoepM9Ds6Oiocn7OnDmYO3euWuWlpKQAAGrVqqVyvlatWuJ7hmDwpLB7925s374dO3bsgIeHB6KjozFx4kQ4ODhg6NChxa4PDQ3FvHnzDFBTInqXqZMUEhMTYW1tLZ43MzPTad30yeDdR9OmTcOMGTMwYMAAeHp64rPPPsOkSZMQGhpa4vXBwcHIyMgQj8TERD3XmIjeRep0H1lbW6scUpKCvb09AODJkycq5588eSK+ZwgGTwrZ2dnFBndMTExKnTloZmZW7C+EiEhTJiYmZTq0xcXFBfb29oiIiBDPKRQKXLx4Ed7e3lorR10G7z7q0aMHFi5cCCcnJ3h4eCAqKgorVqzA8OHDDV01IqpAdDF5LTMzE3FxceLr+Ph4REdHo0aNGnBycsLEiRMREhKChg0bwsXFBbNmzYKDgwN69+4t5Ra0wuBJYc2aNZg1axY+//xzpKamwsHBAaNHj8bs2bMNXTUiqkB0kRSuXLmCzp07i68nT54MABg6dCjCw8PxxRdfICsrC//+97+Rnp4OHx8fHD16FObm5urfgJbIBCOfOqxQKGBjY4P09HR2JRFVQAqFAtWqVUNGRoakz4Ciz5D69eu/tXtIqVTiwYMHkssyBgZvKRARlQdlWSXVyL9DlwmTAhERmBSKMCkQEYFJoQiTAhERuHR2ESYFIiKwpVCESYGICGwpFGFSICKC8SSF1atXl/na8ePHqx2fSYGICND6Mha6snLlSpXXT58+RXZ2NqpVqwYASE9Ph6WlJezs7Cp2UhAEoUL09xGRKm39vzeWMYX4+Hjxzzt27MB3332HTZs2oXHjxgCA2NhYjBo1CqNHj5YU3+AL4hERlQe62I5T12bNmoU1a9aICQEAGjdujJUrV2LmzJmSYr4zLQUiIk3IZLK3thRKW73ZUB4/foz8/Pxi55VKZbElucuKLQUiIhhnS6FLly4YPXo0rl27Jp67evUqgoKCVLb5VAeTAhER9L+fgjZs3rwZ9vb2aNmyJczMzGBmZobWrVujVq1a+P777yXFZPcRERGM55HU19na2uLw4cO4e/cu7ty5AwBwc3NDo0aNJMdkUiAiQtmePnrb+4ZSr149CIIAV1dXVK6s2ce6pDtMTU3VqFAiovLGGMcUsrOzMWLECFhaWsLDwwMJCQkAgHHjxmHx4sWSYkpKCv/617+gVCpLfK+kkXAiovLOGJNCcHAwrl+/jlOnTqns1ta1a1fs2rVLUkxJSaFatWolzpR7/vy55BFvIiJDMsaB5gMHDuDbb7+Fj4+PSsLy8PDA/fv3JcWUlBS2bduG48ePY/PmzeK5mJgYtG7dGlWqVFErVr169UrMxmPGjJFSNSIiSYrGFN52lCdPnz6FnZ1dsfNZWVmSWzWSWwp79+7FtGnTcOnSJRw7dgze3t7o3bs3fv31V7ViXb58GY8fPxaP48ePAwD69esnpWpERJIYY/dRy5YtcejQIfF1Uf2+//57eHt7S4pZ5mHqPn36oFmzZuLh6emJb7/9Fh999BFevXqFNWvWIDAwUO0K2NraqrxevHgxXF1d0bFjR7VjERFJZYyPpC5atAj+/v64ffs28vPzsWrVKty+fRt//PEHIiMjJcUsc0vB1dUVZ86cwciRI1GvXj3UrFkTGzduhCAIGDRoEFq0aIG8vDxJlSiSm5uLH3/8EcOHDy/1l5+TkwOFQqFyEBFpqmiZi787yltS8PHxQXR0NPLz8+Hp6YnffvsNdnZ2OH/+PLy8vCTFLHNL4euvvxb//OjRI0RHRyM6Oho1a9bEyZMnsWnTJlSuXBlubm64fv26pMocOHAA6enpGDZsWKnXhIaGYt68eZLiExGVpiwDyeVt7SOg8Av7xo0btRZP0iyHOnXqoE6dOujevbt4LjMzE9HR0ZITAgBs2rQJ/v7+cHBwKPWa4OBgTJ48WXytUCjg6OgouUwiIsA4u4+AwkQVFxeH1NTUYkmrQ4cOasfT2ozmqlWrwsfHBz4+PpJ+/uHDhzhx4gT27dv3t9cVre9BRKRNxjij+cKFCxg0aBAePnxYbK8HmUxW6nyyv1NulrnYsmUL7OzsVFofRET6Yowthf/85z/iE0i1a9fWSv3KRVIoKCjAli1bMHToUI3X7SAiksIYk8K9e/fw888/o0GDBlqLWS7aQidOnEBCQgKGDx9u6KoQUQVljJPX2rRpg7i4OK3GVPsO8/Ly0KVLF9y7d09rlfD19YUgCBot90pEpAldTF5TKpWYNWsWXFxcYGFhAVdXVyxYsEBrez2PGzcOU6ZMQXh4OK5evYobN26oHFKo3VdjamoquTAiovJKF91HS5Yswbp167B161Z4eHjgypUrCAwMhI2NTYnrx6mrb9++AKDSyyKTySAIgn4HmgcPHoxNmzZJXpqViKi80UVS+OOPP9CrVy/xAZp69eph586duHTpkuR6vi4+Pl4rcV4nKSnk5+dj8+bNOHHiBLy8vIotgrdixQqtVI6ISF/USQpvrqRQ2qPybdu2xYYNG3D37l00atQI169fx9mzZ7X2Gens7KyVOK+TlBRu3bqFFi1aAADu3r2r8l55G50nIioLdWY0vzlhds6cOZg7d26x62fMmAGFQgE3NzeYmJhAqVRi4cKFCAgIkFzPgwcPwt/fH6ampjh48ODfXtuzZ0+140tKCidPnpTyY0RE5VpZv9QmJibC2tpafF3ahNrdu3dj+/bt2LFjBzw8PBAdHY2JEyfCwcEBQ4cOlVTH3r17IyUlBXZ2dujdu3ep1xn95DUiIkNSp/vI2tpaJSmUZtq0aZgxYwYGDBgAAPD09MTDhw8RGhoqOSm8vpSFLtZikvzQ7ZkzZzB48GB4e3vj0aNHAIAffvgBZ8+e1VrliIj0RRePpGZnZxeb22BiYlIuF9YrIqmlsHfvXnz22WcICAhAVFQUcnJyAAAZGRlYtGgRDh8+rNVKEhHpmi6ePurRowcWLlwIJycneHh4ICoqCitWrNDqRN2srCxERkYiISEBubm5Ku9JeexVUlIICQlBWFgYhgwZgp9++kk8365dO4SEhEgJqTGlUimp/0wqbU0+IdI1fT/8oe/ytPWtWxcL4q1ZswazZs3C559/jtTUVDg4OGD06NGYPXu2JlUVRUVF4aOPPkJ2djaysrJQo0YNPHv2DJaWlrCzs5OUFCR1H8XGxpa4JKuNjQ3S09OlhCQiMihddB9ZWVnhm2++wcOHD/HXX3/h/v37CAkJgVwu10qdJ02ahB49euDFixewsLDAhQsX8PDhQ3h5eWHZsmWSYkpKCvb29iWut3H27FnUr19fUkWIiAzJGPdojo6OxpQpU1CpUiWYmJggJycHjo6OWLp0Kb788ktJMSUlhVGjRmHChAm4ePEiZDIZkpOTsX37dkydOhVBQUGSKkJEZEjGmBRMTU3FLi07OzskJCQAKOy1SUxMlBRT0pjCjBkzUFBQgC5duiA7OxsdOnSAmZkZpk6dinHjxkmqCBGRIRnj0tnNmzfH5cuX0bBhQ3Ts2BGzZ8/Gs2fP8MMPP+D999+XFFNSS0Emk+Grr75CWloabt26hQsXLuDp06dYsGCBpEoQERmaMS6dvWjRItSuXRsAsHDhQlSvXh1BQUF4+vQpNmzYICmmRpPX5HI53N3dNQlBRFQuGGNLoWXLluKf7ezscPToUY1jljkpTJ48ucxBuSAeERkbY0wKulDmpBAVFaXy+tq1a8jPz0fjxo0BFC6MZ2JiAi8vL+3WkIhID4wlKTRv3rzM9bh27Zra8cucFF5fBG/FihWwsrLC1q1bUb16dQDAixcvEBgYiPbt26tdiUePHmH69Ok4cuQIsrOz0aBBA2zZskWlaUREpEvGkhT+bhE8bZA0prB8+XL89ttvYkIAgOrVqyMkJAS+vr6YMmVKmWO9ePEC7dq1Q+fOnXHkyBHY2tri3r17KrGJiHRNFzOadWHOnDk6jS8pKSgUCjx9+rTY+adPn+Lly5dqxVqyZAkcHR2xZcsW8ZyLi4uUahERaaQ8tASkuHLlCmJiYgAA7u7uGnXjS0p7n3zyCQIDA7Fv3z4kJSUhKSkJe/fuxYgRI9CnTx+1Yh08eBAtW7ZEv379YGdnh+bNm2Pjxo1SqkVEJJkxTl5LSkpC+/bt0bp1a0yYMAETJkxAq1at4OPjg6SkJEkxJSWFsLAw+Pv7Y9CgQXB2doazszMGDRoEPz8/fPfdd2rFevDgAdatW4eGDRvi2LFjCAoKwvjx47F169YSr8/JyYFCoVA5iIg0ZYxJYeTIkcjLy0NMTAzS0tKQlpaGmJgYFBQUYOTIkZJiqt19lJeXhx49eiAsLAxff/017t+/DwBwdXUttldzWRQUFKBly5ZYtGgRgMKR9Vu3biEsLKzETShCQ0Mxb948tcshIvo7xjLQ/LrIyEj88ccf4lOgANC4cWOsWbNG0kM/gISWgqmpKW7cuAEAqFKlCpo0aYImTZpISggAULt27WIT4P7xj3+Ia3i8KTg4GBkZGeIhdX0PIqLXGeOMZkdHR+Tl5RU7r1Qq4eDgICmmpDscPHgwNm3aJKnAN7Vr1w6xsbEq5+7evQtnZ+cSrzczMxO3wivrlnhERG9jjEnh66+/xrhx43DlyhXx3JUrVzBhwgTJS2dLevooPz8fmzdvxokTJ+Dl5VWslaDOjOZJkyahbdu2WLRoET799FNcunQJGzZskLxuBxGRFMbYfTRs2DBkZ2ejTZs2qFy58OM8Pz8flStXxvDhw1V2eEtLSytTTElJ4datW2jRogWAwm/1r1P3l9aqVSvs378fwcHBmD9/PlxcXPDNN98gICBAStWIiCQxxqTwzTffaD2mpKTw+uxmbfj444/x8ccfazUmEZE6jDEplPQwjqbKVwcZEZGBGOMjqeHh4SWez8/PR3BwsKSYkpNCeno6li9fjpEjR2LkyJFYsWIFMjIypIYjIjIoYxxoHj9+PPr164cXL16I52JjY9GmTRvs3LlTUkxJd3jlyhW4urpi5cqV4oSJlStXwtXVVdKqfEREhmaMLYWoqCgkJSXB09MTx48fx9q1a9GiRQu4ubnh+vXrkmJKGlOYNGkSevbsiY0bN6qMeI8cORITJ07E6dOnJVWGiMhQjHFMwdXVFefOncPEiRPh5+cHExMTbN26FQMHDpQcU3JLYfr06WJCAIDKlSvjiy++UHlelojIWBhjSwEADh06hJ9++gne3t6oVq0aNm3ahOTkZMnxJCUFa2vrEmccJyYmwsrKSnJliIgMxRiTwujRo9GvXz9Mnz4dZ86cwY0bNyCXy+Hp6Yndu3dLiimp+6h///4YMWIEli1bhrZt2wIAzp07h2nTpmnUbCEiMpRKlSrBxMTkrdeUJ+fOncPFixfRtGlTAIC9vT0OHz6MtWvXYvjw4fj000/VjikpKSxbtgwymQxDhgxBfn4+gMI1kYKCgrB48WIpIYmIDMoYxxSuXr0KMzOzYufHjBmDrl27SoopKSnI5XKsWrUKoaGhKqukWlpaSqqENrx69QpyuVxv5RUlQ30pKCjQa3mGKFMQBL2WB+j/P/nr43Dvaplv+7atbbm5uVqJY4xJwczMDPfv38eWLVtw//59rFq1CnZ2djhy5AicnJwkxdSoLWRpaQlPT094enoaNCEQEWlKV2MKjx49wuDBg1GzZk1YWFjA09NTaw/kREZGwtPTExcvXsS+ffuQmZkJALh+/brkbTvLVwcZEZGB6CIpFO1Bb2pqiiNHjuD27dtYvny51vagnzFjBkJCQnD8+HGVnpIPP/wQFy5ckBRT/21ZIqJyqCwzltUdaNb1HvQ3b97Ejh07ip23s7PDs2fPJMVkS4GICLppKeh6D/pq1arh8ePHxc5HRUWhTp06kmIyKRARQb2k8OY+8Tk5OSXGVHcPenUNGDAA06dPR0pKCmQyGQoKCnDu3DlMnToVQ4YMkRRTcvdRREQEIiIikJqaWuwplc2bN0sNS0RkEOp0Hzk6OqqcnzNnDubOnVvsenX3oFfXokWLMGbMGDg6OkKpVMLd3R1KpRKDBg3CzJkzJcWUlBTmzZuH+fPno2XLlqhdu3a5e0yLiEhd6jySmpiYqLIVcElzBYDS96Dfu3evhrUtJJfLsXHjRsyePRs3b95EZmYmmjdvjoYNG0qOKSkphIWFITw8HJ999pnkgomIyhN1Wgpl3R9e3T3opXJ0dCzWepFK0phCbm6uuLwFEdG7QBcDzZMmTcKFCxewaNEixMXFYceOHdiwYQPGjBmjo7vQnKSkMHLkyBIfg5Ji7ty5xX7pbm5uWolNRFRWukgKRXvQ79y5E++//z4WLFhQ7vegl9R99OrVK2zYsAEnTpxAkyZNYGpqqvL+ihUr1Irn4eGBEydO/F+lDLAUABGRLsZHjW0Pekmfvjdu3ECzZs0AALdu3VJ5T8ovtXLlyrC3t5dSFSIirTDGtY90QVJSOHnypFYrce/ePTg4OMDc3Bze3t4IDQ0tdTGnnJwclWeCFQqFVutCRBWTLmY068OZM2ewfv163L9/Hz///DPq1KmDH374AS4uLvDx8VE7nsHvsE2bNggPD8fRo0exbt06xMfHo3379nj58mWJ14eGhsLGxkY8tDXiTkQVW1FSeNtRnuzduxfdunWDhYUFoqKixC/MGRkZ4twIdUm+w/T0dCxfvhwjR47EyJEjsWLFCmRkZKgdx9/fH/369UOTJk3QrVs3HD58GOnp6aXuGhQcHIyMjAzxSExMlHoLREQiY9x5LSQkBGFhYdi4caPK2G67du1w7do1STEl79Hs6uqKlStXIi0tDWlpaVi5ciVcXV0lV6RItWrV0KhRI8TFxZX4vpmZmfiMcFmfFSYiehtjTAqxsbHo0KFDsfM2NjZIT0+XFFNSUpg0aRJ69uyJP//8E/v27cO+ffsQHx+Pjz/+GBMnTpRUkSKZmZm4f/8+ateurVEcIiJ1GGNSsLe3L/EL9NmzZ1G/fn1JMSW3FKZPn67y6GjlypXxxRdfqL15xNSpUxEZGYk///wTf/zxBz755BOYmJhwr2ci0isTE5MyHeXJqFGjMGHCBFy8eBEymQzJycnYvn07pk6diqCgIEkxJT19ZG1tjYSEhGKTzBITE2FlZaVWrKSkJAwcOBDPnz+Hra0tfHx8cOHCBdja2kqpGhGRJMb4SOqMGTNQUFCALl26IDs7Gx06dICZmRmmTp2KcePGSYopKSn0798fI0aMwLJly8TlLs6dO4dp06ap/Q3/p59+klIFIiKtMsakIJPJ8NVXX2HatGmIi4tDZmYm3N3dUbVqVckxJSWFZcuWQSaTYciQIcjPz4cgCJDL5QgKCsLixYslV4aIyFCMMSlkZGRAqVSiRo0aKquxpqWloXLlypIexJE0piCXy7Fq1Sq8ePEC0dHRuH79uvgEUmlLyBIRlWfGOE9hwIABJfa27N69GwMGDJAUs8wthcmTJ2PBggWoUqUKJk+e/LfXqrv2ERGRoclksrd+6Je3lsLFixdL/Lzt1KkTvvrqK0kxy5wUoqKikJeXJ/65NOXtl0ZEVBbG2H2Uk5OD/Pz8Yufz8vLw119/SYpZ5qTw+npH2l77iIjI0Ixx7aPWrVtjw4YNWLNmjcr5sLAweHl5SYopaaA5ISEBjo6OJWbNhISEUhezIyIqr4yxpRASEoKuXbvi+vXr6NKlCwAgIiICly9fxm+//SYppqS05+LigqdPnxY7//z5c7i4uEiqCBGRIRnjjOZ27drh/PnzcHR0xO7du/Hrr7+iQYMGuHHjBtq3by8ppqSWgiAIJf5yMjMzYW5uLqkiRESGZIzdRwDQrFkzbN++XWvx1EoKRU8dyWQyzJo1C5aWluJ7SqUSFy9eFDff0bekpCSNJmyoq6SWki49e/ZMr+UBwOPHj/VaXm5url7LA6D2DHxNNWjQQK/lAYCDg4Ney6tWrZpeyyttmX11GWtSKCgoQFxcHFJTU1FQUKDyXkmL5b2NWkmh6KkjQRBw8+ZNyOVy8T25XI6mTZti6tSpaleCiMjQjHFM4cKFCxg0aBAePnwIQRBU3pPJZFAqlWrHVCspFD11FBgYiNWrV+v9WxYRka4YY1L4z3/+g5YtW+LQoUOoXbu2VuonaUyhYcOG2LNnD4YPH65yfvPmzXj69CmmT5+uccWIiPTJGJPCvXv38PPPP2u1W1JSB9mGDRuKrZAKAB4eHggLC9O4UkRE+maMS2e3adOm1A3JpJLUUkhJSSlxExxbW1u9D04SEWmDMbYUxo0bhylTpiAlJQWenp4qW3ICQJMmTdSOKSkpODo64ty5c8XmJJw7d07vTzoQEWmDMSaFvn37AoBKV75MJhOnDeh8oLnIqFGjMHHiROTl5eHDDz8EUDiL7osvvsCUKVOkhCQiMihjTArx8fFajykpKUybNg3Pnz/H559/Lj5bbm5ujunTpyM4OFirFSQi0gdjXCXV2dlZ6zElJQWZTIYlS5Zg1qxZiImJgYWFBRo2bMi9FIjIaBnr5DUAuH37NhISEopNAO3Zs6fasSQlhSJVq1ZFq1atNAmhYvHixQgODsaECRPwzTffaC0uEdHbGGP30YMHD/DJJ5/g5s2b4lgC8H/1lDKmoFHau337No4ePYqDBw+qHFJcvnwZ69evlzRaTkSkKV3vvLZ48WLIZDJMnDhRa3WeMGECXFxckJqaCktLS/zvf//D6dOn0bJlS5w6dUpSTEktBW1np8zMTAQEBGDjxo0ICQmRUiUiIo3osqWgqy+958+fx++//4733ntPTFo+Pj4IDQ3F+PHj/3ZDtNJISnvazk5jxoxB9+7d0bVr17dem5OTA4VCoXIQEZVXr3/prV69ulZjK5VKcbmh9957D8nJyQAKB6BjY2MlxZSUFM6fP4/58+eXmp3U8dNPP+HatWsIDQ0t0/WhoaGwsbERD0dHRym3QESkQp39FN78YpqTk1NqXHW+9Krr/fffx/Xr1wEUzm5eunQpzp07h/nz56N+/fqSYkpKCtrKTomJiZgwYQK2b99e5n0YgoODkZGRIR6JiYnq3wAR0RvUGVNwdHRU+XJa2pdadb/0qmvmzJnictnz589HfHw82rdvj8OHD2P16tWSYkoaUyjKTi4uLmJ2ksvl2LBhg1rZ6erVq0hNTUWLFi3Ec0qlEqdPn8a3336LnJycYmuNmJmZ8dFXItI6dcYUEhMTYW1tLZ4v6TOp6Evv8ePHdbb5WLdu3cQ/N2jQAHfu3EFaWhqqV68uefxDUlKYOXMmsrKyABRmp48//hjt27dHzZo1sWvXrjLH6dKlC27evKlyLjAwEG5ubpg+fXq5W3yKiN5d6iQFa2trlaRQEilferWhRo0aGv28pKSgrexkZWWF999/X+VclSpVULNmzWLniYh0SdtPH+nqS2+fPn0QHh4Oa2tr9OnT52+v3bdvn9rx1U4KeXl58PPzQ1hYGBo2bCie1zQ7EREZkraTgq6+9NrY2Ij1sLGxkRynNGonBVNTU9y4cUPrFSkidcIFEZEmjGVG85YtWwAUbos8b9482NrawsLCQmvxJT19NHjwYGzatElrlSAiMjR1HkmV6tSpU1pbwkcQBDRo0ABJSUlaiVdE0phCfn4+Nm/ejBMnTsDLywtVqlRReX/FihVaqRwRkb4YS0uhSKVKldCwYUM8f/5cpStfU5KSwq1bt8QR9bt376q8V55+aURE6jC2z6/Fixdj2rRpWLdundYezpGUFE6ePKmVwomIygtjaykAwJAhQ5CdnY2mTZtCLpcXG1tIS0tTO6bkpbPPnDmD9evX48GDB9izZw/q1KmDH374AS4uLvDx8ZEalojIIIxxPwVdbDEgKSns3bsXn332GQICAnDt2jVx3Y+MjAwsWrQIhw8f1moliYiouKFDh2o9pqS0FxISgrCwMGzcuBGmpqbi+Xbt2uHatWtaqxwRkb7o4+kjXXr16pVWVpCW1FKIjY1Fhw4dip23sbFBenq6pIpo6uTJkzpbX6Qk+m4NGSLZPn/+XK/lFS2yqE/NmzfXa3m9evXSa3lA4Q6J73J5+fn5WoljjGMKWVlZmD59Onbv3l3i/1e97bxmb2+PuLi4YufPnj0reblWIiJDMsaWwhdffIHff/8d69atg5mZGb7//nvMmzcPDg4O2LZtm6SYkloKo0aNwoQJE7B582bIZDIkJyfj/PnzmDp1KmbNmiWpIkREhmSMA82//vortm3bhk6dOiEwMBDt27dHgwYN4OzsjO3btyMgIEDtmJKSwowZM1BQUIAuXbogOzsbHTp0gJmZGaZOnYpx48ZJCUlEZFDG2H2UlpYm9s5YW1uLj6D6+PggKChIUkxJaU8mk+Grr75CWloabt26hQsXLuDp06dYsGCBpEoQERmaMXYf1a9fH/Hx8QAANzc37N69G0BhC6JatWqSYkpKCgkJCRAEAXK5HO7u7mjdurU4uJSQkCCpIkREhmSMSSEwMFDcjnPGjBlYu3YtzM3NMWnSJEybNk1STEndRy4uLnj8+DHs7OxUzj9//hwuLi6SRryJiEg9kyZNEv/ctWtX3LlzB1evXkWDBg3QpEkTSTElJQVBEErMmJmZmXp9LJSISFuMcaA5MTERjo6O4mtnZ2c4OztrFFOtpDB58mQAhc2sWbNmwdLSUnxPqVTi4sWLaNasmUYVIiIyBGMcaK5Xrx58fHwwePBg/Otf/0L16tU1jqlWUoiKigJQ2FK4efMm5HK5+J5cLkfTpk0xdepUjStFRKRvxpgUrly5gh07dmD+/PkYN24c/Pz8MHjwYPTo0QNmZmaSYqqVFIpWRw0MDMTq1au1MgN13bp1WLduHf78808AgIeHB2bPng1/f3+NYxMRlZUxJoXmzZujefPmWLp0KU6dOoUdO3bg3//+NwoKCtCnTx9s3rxZ7ZiSxhS2bNmCiIgIREREIDU1FQUFBSrvq1ORunXrYvHixWjYsCEEQcDWrVvRq1cvREVFwcPDQ0r1iIjUZoxJoYhMJkPnzp3RuXNnBAUFYcSIEdi6daukpCBp1GT+/Pnw9fVFREQEnj17hhcvXqgc6ujRowc++ugjNGzYEI0aNcLChQtRtWpVXLhwQUrViIgkM6bHUV+XlJSEpUuXolmzZuIUgbVr10qKJamlsG7dOoSHh+Ozzz6TVGhplEol9uzZg6ysLHh7e5d4TU5OjrhUNwDJKwESEb3OGFsK69evx44dO3Du3Dm4ubkhICAAv/zyi0ZPIElKCrm5uWjbtq3kQt908+ZNeHt749WrV6hatSr2798Pd3f3Eq8NDQ3FvHnztFY2EZGxCgkJwcCBA7F69Wo0bdpUKzEldR+NHDkSO3bs0EoFAKBx48aIjo7GxYsXERQUhKFDh+L27dslXhscHIyMjAzxSExM1Fo9iKjiMsYZzQkJCVi6dKnWEgIgsaXw6tUrbNiwASdOnECTJk1UNtoBgBUrVqgVTy6Xo0GDBgAALy8vXL58GatWrcL69euLXWtmZib5USsiotIYY/eRLuojKSncuHFDnKR269Ytlfe0UcmCggKVcQMiIl0zxqSgC5KSQtF8BW0IDg6Gv78/nJyc8PLlS+zYsQOnTp3CsWPHtFYGEdHbMCkUkpQUtCk1NRVDhgzB48ePYWNjgyZNmuDYsWP45z//aeiqEVEFwqRQSK2k0KdPnzJdt2/fvjLH3LRpkzpVICKiEjx79gwXL16EUqlEq1atULt2bUlx1EoKNjY2kgohIirvdNFSCA0Nxb59+3Dnzh1YWFigbdu2WLJkCRo3bqxJVYvZu3cvRowYgUaNGiEvLw+xsbFYu3YtAgMD1Y6lVlLYsmWL2gUQERkDXSydHRkZiTFjxqBVq1bIz8/Hl19+CV9fX9y+fRtVqlSRXNfMzExxYzMAmDdvHi5duoRGjRoBAA4dOoRRo0ZJSgrla3FwIqJ3yNGjRzFs2DB4eHigadOmCA8PR0JCAq5evapRXC8vL/zyyy/i68qVKyM1NVV8/eTJE5VVrNVh8IFmIqLyQB8DzRkZGQCAGjVqaBTn2LFjGDNmDMLDw7F27VqsWrUK/fv3h1KpRH5+PipVqoTw8HBJsZkUiIigXlJ4c821skyqLSgowMSJE9GuXTu8//77GtW1Xr16OHToEHbu3ImOHTti/PjxiIuLQ1xcHJRKJdzc3CTvgsnuIyIiqLfMhaOjI2xsbMQjNDT0rfHHjBmDW7du4aefftJanQcOHIjLly/j+vXr6NSpEwoKCtCsWTONtkVmS4GICOq1FBITE2FtbS2ef1srYezYsfjvf/+L06dPo27duppXFsDhw4cRExODpk2b4vvvv0dkZCQCAgLg7++P+fPnw8LCQlJcthSIiKBeS8Ha2lrlKC0pCIKAsWPHYv/+/fj999/h4uKilbpOmTIFgYGBuHz5MkaPHo0FCxagY8eOuHbtGszNzdG8eXMcOXJEUux3pqWQkZGh1/WS/ve//+mtLAAwMTHRa3kA0K9fP72W5+TkpNfyAKBjx456LU9bHwrqeO+99/RaniaPWkpRnmcZjxkzBjt27MAvv/wCKysrpKSkACic8yX1mzwAhIeH47fffoOXlxfS0tLwwQcfYNasWZDL5ViwYAEGDhyI0aNHS9rWmC0FIiLoZunsdevWISMjA506dULt2rXFY9euXRrVtUqVKoiPjwdQ2JX15hiCu7s7zpw5Iyn2O9NSICLShC4eSRUEQZMqlSo0NBRDhgzB+PHjkZ2dja1bt2otNpMCERGMa0G8gIAA+Pn54cGDB2jYsCGqVaumtdhMCkRE/195+dAvi5o1a6JmzZpaj8ukQEQE42op6BIHmomISMSWAhER2FIowqRARAQmhSJMCkREYFIoYvAxhdDQULRq1QpWVlaws7ND7969ERsba+hqERFVSAZPCkU7E124cAHHjx9HXl4efH19kZWVZeiqEVEFoosZzcbI4N1HR48eVXkdHh4OOzs7XL16FR06dDBQrYiIKiaDJ4U3vW1nopycHJWF797c7IKISAqOKRQyePfR68qyM1FoaKjK5haOjo56riURvYvYfVSoXCWFsuxMFBwcjIyMDPFITEzUYw2JiN5t5ab7qKw7E5VlL1QiInWx+6iQwZOCIAgYN24c9u/fj1OnThlkExIiIipk8KSgq52JiIjUwZZCIYOPKehqZyIiInVwoLmQwVsKutqZiIiI1GfwlgIREZUfBm8pEBGVBxxTKMSWAhERidhSICICWwpFmBSIiMCkUITdR0REJGJLgYgIbCkUeWeSgp2dnV5nQNva2uqtLABwc3PTa3kA4OHhodfy/vnPf+q1PABwcHDQa3mGWLfL2tpar+VVrqzfj5Xc3Fy9lveue2eSAhGRJthSKMSkQEQEJoUiHGgmItKhtWvXol69ejA3N0ebNm1w6dIlQ1fpbzEpEBFBNwvi7dq1C5MnT8acOXNw7do1NG3aFN26dUNqaqqO7kJzTApERDqyYsUKjBo1CoGBgXB3d0dYWBgsLS2xefNmQ1etVBxTICIC8PLly7e2BF6+fAkAUCgUKudL2hEyNzcXV69eRXBwsHiuUqVK6Nq1K86fP6+lWmsfkwIRVWhyuRz29vZwdHQs0/VVq1Ytdu2cOXMwd+5clXPPnj2DUqlErVq1VM7XqlULd+7c0ajOusSkQEQVmrm5OeLj48s830EQhGItindp33gmBSKq8MzNzWFubq7VmO+99x5MTEzw5MkTlfNPnjyBvb29VsvSJoMPNJ8+fRo9evSAg4MDZDIZDhw4YOgqERFpTC6Xw8vLCxEREeK5goICREREwNvb24A1+3sGTwpZWVlo2rQp1q5da+iqEBFp1eTJk7Fx40Zs3boVMTExCAoKQlZWFgIDAw1dtVIZvPvI398f/v7+hq4GEZHW9e/fH0+fPsXs2bORkpKCZs2a4ejRo8UGn8sTgycFdeXk5CAnJ0d8/eajYURE5cnYsWMxduxYQ1ejzAzefaSu0NBQ2NjYiEdZHyMjIqK3M7qkEBwcjIyMDPFITEw0dJWIiN4ZRtd9VNLMQSIi0g6jaykQEZHuGLylkJmZibi4OPF1fHw8oqOjUaNGDTg5ORmwZkREFY/Bk8KVK1fQuXNn8fXkyZMBAEOHDkV4eLiBakVEVDEZPCl06tQJgiAYuhpERASOKRAR0WuYFIiISMSkQEREIiYFIiISMSkQEZGISYGIiERMCkREJDL4PAVNFc1xePXqlV7LVSqVei0vLy9Pr+UB+v+dZmVl6bU8AHj58qVeyyvrPsDGrHJl/X6sFP0dcr6TdsgEI/9NJiUlcflsIkJiYiLq1q1r6GoYPaNPCgUFBUhOToaVlRVkMlmZf06hUMDR0RGJiYmwtrbWYQ0rVpm8x3ejTGO6R0EQ8PLlSzg4OKBSJfaIa8rou48qVaqk0bcDa2trvf2jr0hl8h7fjTKN5R5tbGx0VJuKh2mViIhETApERCSqsEnBzMwMc+bM0esubhWhTN7ju1FmRbhHKpnRDzQTEZH2VNiWAhERFcekQEREIiYFIiISMSkQEZGowiSFffv2wdfXFzVr1oRMJkN0dHSxa0aPHg1XV1dYWFjA1tYWvXr1wp07d3RSXlpaGsaNG4fGjRvDwsICTk5OGD9+PDIyMiSVV5YyAWDDhg3o1KkTrK2tIZPJkJ6ertPyXr16hTFjxqBmzZqoWrUq+vbtiydPnkgu801PnjzBsGHD4ODgAEtLS/j5+eHevXtai1+SzMxMjB07FnXr1oWFhQXc3d0RFham0zJlMlmJx9dff62zMmNiYtCzZ0/Y2NigSpUqaNWqFRISEnRW3rBhw4rdn5+fn87Ko5JVmKSQlZUFHx8fLFmypNRrvLy8sGXLFsTExODYsWMQBAG+vr6SFr97W3nJyclITk7GsmXLcOvWLYSHh+Po0aMYMWKE2mWVtUwAyM7Ohp+fH7788kvJ5ahT3qRJk/Drr79iz549iIyMRHJyMvr06aNx2UDh8ga9e/fGgwcP8MsvvyAqKgrOzs7o2rWrThfXmzx5Mo4ePYoff/wRMTExmDhxIsaOHYuDBw/qrMzHjx+rHJs3b4ZMJkPfvn11Ut79+/fh4+MDNzc3nDp1Cjdu3MCsWbNgbm6uk/KK+Pn5qdznzp07dVoelUCoYOLj4wUAQlRU1FuvvX79ugBAiIuL00t5u3fvFuRyuZCXlye5vLKWefLkSQGA8OLFC43K+rvy0tPTBVNTU2HPnj3iuZiYGAGAcP78eY3LjY2NFQAIt27dEs8plUrB1tZW2Lhxo8bxS+Ph4SHMnz9f5VyLFi2Er776SmdlvqlXr17Chx9+qLP4/fv3FwYPHqyz+CUZOnSo0KtXL72WScVVmJaCurKysrBlyxa4uLjobRXWjIwMWFtb633pYV25evUq8vLy0LVrV/Gcm5sbnJyccP78eY3j5+TkAIDKt9dKlSrBzMwMZ8+e1Th+adq2bYuDBw/i0aNHEAQBJ0+exN27d+Hr66uzMl/35MkTHDp0SKNW5d8pKCjAoUOH0KhRI3Tr1g12dnZo06YNDhw4oJPyXnfq1CnY2dmhcePGCAoKwvPnz3VeJqliUnjDd999h6pVq6Jq1ao4cuQIjh8/DrlcrvNynz17hgULFuDf//63zsvSl5SUFMjlclSrVk3lfK1atZCSkqJx/KIEExwcjBcvXiA3NxdLlixBUlISHj9+rHH80qxZswbu7u6oW7cu5HI5/Pz8sHbtWnTo0EFnZb5u69atsLKy0lo33JtSU1ORmZmJxYsXw8/PD7/99hs++eQT9OnTB5GRkTopEyjsOtq2bRsiIiKwZMkSREZGwt/fX+97l1R072RS2L59u/jBXrVqVZw5c6bMPxsQEICoqChERkaiUaNG+PTTT9+62Ywm5QGFSwZ3794d7u7umDt3bpl+RtMy1aXv8spShwsXLmDfvn24e/cuatSoAUtLS5w8eRL+/v5aW0K5pPtes2YNLly4gIMHD+Lq1atYvnw5xowZgxMnTuiszNdt3rwZAQEBWuvff7O82NhYAECvXr0wadIkNGvWDDNmzMDHH3+stQH1ku5xwIAB6NmzJzw9PdG7d2/897//xeXLl3Hq1CmtlEll8270U7yhZ8+eaNOmjfi6Tp06Zf5ZGxsb2NjYoGHDhvjggw9QvXp17N+/HwMHDtRJeS9fvoSfnx+srKywf/9+mJqalunnNClTCinl2dvbIzc3F+np6SqthSdPnsDe3l4rdbCwsEB0dDQyMjKQm5sLW1tbtGnTBi1btlQ7flnL7NKlC/bv34/u3bsDAJo0aYLo6GgsW7ZMpatMm2UWOXPmDGJjY7Fr1y6NyymtPFtbW1SuXBnu7u4q1/3jH//QWrdcWf491a9fH++99x7i4uLQpUsXrZRLb/dOJgUrKytYWVlpHEcQBAiCIPZda7s8hUKBbt26wczMDAcPHlTrm5+27lGX5Xl5ecHU1BQRERHiUzKxsbFISEiAt7e3VutQtJ7+vXv3cOXKFSxYsEDt+GUpU6FQIC8vr1hLxMTEBAUFBTop83WbNm2Cl5cXmjZtqpWySiuvVatWYouhyN27d+Hs7KyzMt+UlJSE58+fo3bt2lopk8rmnUwKJUlLS0NCQgKSk5MBQPwHb29vD3t7ezx48AC7du2Cr68vbG1tkZSUhMWLF8PCwgIfffSR1stTKBTw9fVFdnY2fvzxRygUCigUCgCF39RMTEy0XiZQ2M+fkpKCuLg4AMDNmzdhZWUFJycn1KhRQ6vl2djYYMSIEZg8eTJq1KgBa2trjBs3Dt7e3vjggw/Uvr+S7NmzB7a2tnBycsLNmzcxYcIE9O7dW2eDvtbW1ujYsSOmTZsGCwsLODs7IzIyEtu2bcOKFSt0UmYRhUKBPXv2YPny5TotBwCmTZuG/v37o0OHDujcuTOOHj2KX3/9VWddOZmZmZg3bx769u0Le3t73L9/H1988QUaNGiAbt266aRMKoWBn37Smy1btggAih1z5swRBEEQHj16JPj7+wt2dnaCqampULduXWHQoEHCnTt3dFJe0SOhJR3x8fE6KVMQBGHOnDklXrNlyxadlPfXX38Jn3/+uVC9enXB0tJS+OSTT4THjx9Lur+SrFq1Sqhbt65gamoqODk5CTNnzhRycnK0Fr8kjx8/FoYNGyY4ODgI5ubmQuPGjYXly5cLBQUFOi13/fr1goWFhZCenq7Tcops2rRJaNCggWBubi40bdpUOHDggM7Kys7OFnx9fQVbW1vB1NRUcHZ2FkaNGiWkpKTorEwqGZfOJiIi0Tv59BEREUnDpEBERCImBSIiEjEpEBGRiEmBiIhETApERCRiUiAiIhGTAhklmUyml6WciSoaJgUql54+fYqgoCA4OTnBzMwM9vb26NatG86dOwegcCcyf39/AMCff/5Z6nagRKSeCrP2ERmXvn37Ijc3F1u3bkX9+vXx5MkTREREiJuuSFlllYjejstcULmTnp6O6tWr49SpU+jYsWOJ18hkMuzfvx+9e/eGTCZTea9jx47iwm3ff/89li9fjvj4eNSrVw/jx4/H559/rutbIDJabClQuVO08cqBAwfwwQcfwMzM7G+vv3TpElq3bo0TJ07Aw8ND3Clv+/btmD17Nr799ls0b94cUVFRGDVqFKpUqYKhQ4fq41aIjA7HFKjcqVy5MsLDw7F161ZUq1YN7dq1w5dffokbN26UeL2trS0AoGbNmrC3txeXAJ8zZw6WL1+OPn36wMXFBX369MGkSZOwfv16vd0LkbFhUqByqW/fvkhOTsbBgwfh5+eHU6dOoUWLFggPDy/Tz2dlZeH+/fsYMWKEyraPISEhuH//vm4rT2TEOKZARmPkyJE4fvw4Hj58qDKm8Oeff8LFxQVRUVFo1qwZgP/b8vPHH39U2fYRKNwlzcXFxQB3QFT+cUyBjIa7u3uJcxOKxhCUSqV4rlatWnBwcMCDBw8QEBCgryoSGT0mBSp3nj9/jn79+mH48OFo0qQJrKyscOXKFSxduhS9evUqdr2dnR0sLCxw9OhR1K1bF+bm5rCxscG8efMwfvx42NjYwM/PDzk5Obhy5QpevHiByZMnG+DOiIyAIbd9IyrJq1evhBkzZggtWrQQbGxsBEtLS6Fx48bCzJkzhezsbEEQBAGAsH//fvFnNm7cKDg6OgqVKlUSOnbsKJ7fvn270KxZM0EulwvVq1cXOnToIOzbt0/Pd0RkPDimQEREIj59REREIiYFIiISMSkQEZGISYGIiERMCkREJGJSICIiEZMCERGJmBSIiEjEpEBERCImBSIiEjEpEBGRiEmBiIhE/w+MVjl/fzJiyQAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"orders = np.arange(1, 10)\n",
"positions = np.arange(-13, -4)\n",
"ticks = np.arange(9)\n",
"\n",
"fig, axes = plt.subplots(1, 1, figsize=(4, 4))\n",
"im = axes.imshow(\n",
" sites.iloc[::-1, :],\n",
" cmap=\"Greys\",\n",
" vmin=0,\n",
")\n",
"axes.set(\n",
" xticks=ticks,\n",
" xticklabels=positions,\n",
" xlabel=\"Site\",\n",
" yticks=ticks,\n",
" yticklabels=orders[::-1],\n",
" ylabel=\"Interaction order $k$\",\n",
" aspect=\"equal\",\n",
")\n",
"plt.colorbar(im, shrink=0.6, label=\"% variance explained\")\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "f7c279bf",
"metadata": {},
"source": [
"This plot shows that whereas some positions have large additive and pairwise contributions, other sites have more context-dependent contributions to phenotypic variation, with the largest contributions arising from interactions of order 3 or higher.\n",
"\n",
"### How to compute the variance explained by pairwise and higher order interactions between every pair of sites\n",
"\n",
"In this example, we do an finer aggregation and separate the variance explained by pairwise and higher order interactions between pairs of sites on the upper and lower triangular parts of the heatmap representation.\n"
]
},
{
"cell_type": "code",
"execution_count": 88,
"id": "76874f0a",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" site2 \n",
" 0 \n",
" 1 \n",
" 2 \n",
" 3 \n",
" 4 \n",
" 5 \n",
" 6 \n",
" 7 \n",
" 8 \n",
" \n",
" \n",
" site1 \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0.000000 \n",
" 20.235400 \n",
" 19.297067 \n",
" 12.736994 \n",
" 8.528527 \n",
" 7.706634 \n",
" 5.887429 \n",
" 4.593244 \n",
" 3.814681 \n",
" \n",
" \n",
" 1 \n",
" 8.562998 \n",
" 0.000000 \n",
" 28.937652 \n",
" 22.225871 \n",
" 13.779223 \n",
" 9.638418 \n",
" 7.624355 \n",
" 5.565630 \n",
" 4.245642 \n",
" \n",
" \n",
" 2 \n",
" 9.734242 \n",
" 13.710318 \n",
" 0.000000 \n",
" 28.779777 \n",
" 20.837386 \n",
" 13.780875 \n",
" 7.749691 \n",
" 5.786126 \n",
" 4.398405 \n",
" \n",
" \n",
" 3 \n",
" 1.380292 \n",
" 9.409153 \n",
" 6.562947 \n",
" 0.000000 \n",
" 31.867450 \n",
" 23.175051 \n",
" 12.769771 \n",
" 6.189255 \n",
" 5.449330 \n",
" \n",
" \n",
" 4 \n",
" 1.371878 \n",
" 2.865711 \n",
" 4.951972 \n",
" 7.331872 \n",
" 0.000000 \n",
" 31.061032 \n",
" 19.615363 \n",
" 10.645316 \n",
" 6.161446 \n",
" \n",
" \n",
" 5 \n",
" 1.202204 \n",
" 1.320697 \n",
" 2.437803 \n",
" 3.856483 \n",
" 4.213326 \n",
" 0.000000 \n",
" 21.838515 \n",
" 13.053678 \n",
" 7.791366 \n",
" \n",
" \n",
" 6 \n",
" 1.367760 \n",
" 1.333584 \n",
" 0.959993 \n",
" 1.554637 \n",
" 3.281415 \n",
" 1.401816 \n",
" 0.000000 \n",
" 13.716746 \n",
" 7.937204 \n",
" \n",
" \n",
" 7 \n",
" 0.635882 \n",
" 0.887281 \n",
" 0.527592 \n",
" 0.344046 \n",
" 1.860773 \n",
" 1.435976 \n",
" 1.094349 \n",
" 0.000000 \n",
" 9.259406 \n",
" \n",
" \n",
" 8 \n",
" 0.129675 \n",
" 0.408224 \n",
" 0.274862 \n",
" 0.711385 \n",
" 0.446727 \n",
" 1.349362 \n",
" 0.584660 \n",
" 0.498108 \n",
" 0.000000 \n",
" \n",
" \n",
"
\n",
"
"
],
"text/plain": [
"site2 0 1 2 3 4 5 \\\n",
"site1 \n",
"0 0.000000 20.235400 19.297067 12.736994 8.528527 7.706634 \n",
"1 8.562998 0.000000 28.937652 22.225871 13.779223 9.638418 \n",
"2 9.734242 13.710318 0.000000 28.779777 20.837386 13.780875 \n",
"3 1.380292 9.409153 6.562947 0.000000 31.867450 23.175051 \n",
"4 1.371878 2.865711 4.951972 7.331872 0.000000 31.061032 \n",
"5 1.202204 1.320697 2.437803 3.856483 4.213326 0.000000 \n",
"6 1.367760 1.333584 0.959993 1.554637 3.281415 1.401816 \n",
"7 0.635882 0.887281 0.527592 0.344046 1.860773 1.435976 \n",
"8 0.129675 0.408224 0.274862 0.711385 0.446727 1.349362 \n",
"\n",
"site2 6 7 8 \n",
"site1 \n",
"0 5.887429 4.593244 3.814681 \n",
"1 7.624355 5.565630 4.245642 \n",
"2 7.749691 5.786126 4.398405 \n",
"3 12.769771 6.189255 5.449330 \n",
"4 19.615363 10.645316 6.161446 \n",
"5 21.838515 13.053678 7.791366 \n",
"6 0.000000 13.716746 7.937204 \n",
"7 1.094349 0.000000 9.259406 \n",
"8 0.584660 0.498108 0.000000 "
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pairs_pw = sds.calc_site_pairs_variance_perc(v_U_vcs.loc[v_U_vcs['k'] == 2, :])\n",
"pairs_ho = sds.calc_site_pairs_variance_perc(v_U_vcs, min_k=3)\n",
"\n",
"m_bottom = pairs_pw.pivot(index=\"site1\", columns=\"site2\", values=\"variance_perc\")\n",
"m_bottom = m_bottom.reindex(ticks).T.reindex(ticks).fillna(0)\n",
"m_top = pairs_ho.pivot(index='site1', columns='site2', values='variance_perc')\n",
"m_top = m_top.reindex(ticks).T.reindex(ticks).fillna(0).T\n",
"m = m_top + m_bottom\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 89,
"id": "63840aa5",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAE2CAYAAACKiF6uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9K0lEQVR4nO3deVxU1f8/8NeArDKMguCIAqIkLqkYGpJryofFEhe+pWYlStmCWlKWlOUeavvHSItUsjCXQsVULC0QUzQRXFJRED9ACKYICAQo3N8fPJjfDMMyM8wMM/h6Ph738eDeuXPe56LMe849554jEgRBABERkRyTtq4AEREZHiYHIiJSwuRARERKmByIiEgJkwMRESlhciAiIiVMDkREpITJgYiIlHRo6wq0hdraWuTn50MsFkMkErV1dYhISwRBwN27d+Hk5AQTE373bY0HMjnk5+fD2dm5ratBRDqSm5uLHj16tHU1jNoDmRzEYjGAuv9Atra2Oo+3ePFinceQZ2FhobdYI0aM0FssAHBxcdFbLDMzM73FAgBra2u9xdLn/xFTU1O9xbp79y4GDBgg+xsnzT2QyaH+VpKtra1ekoM+/xD1HU+fH2gAYGNjo7dY+k4OHTt21FssS0tLvcXSZ3Kox9vFrcebckREpITJgYiIlDA5EBGREiYHIiJSwuRARERKmByIiEiJQSWHuLg4+Pn5wd7eHiKRCOnp6UrnvPTSS+jduzesrKzg4OCASZMm4fLly/qvLBFRO2ZQyaG8vBwjR47E2rVrmzzHy8sLW7ZswaVLl3Do0CEIggA/Pz/U1NTosaZERO2bQT0E99xzzwEArl+/3uQ5c+fOlf3cs2dPrFq1CoMHD8b169fRu3dvXVeRiOiBYFAtB3WVl5djy5YtcHNz41xJRERaZJTJ4csvv4SNjQ1sbGxw8OBB/PrrrzA3N2/y/KqqKpSWlipsRETUtDZLDrGxsbIPeBsbGyQnJ6v83pkzZyItLQ1JSUno06cPnn76aVRWVjZ5fmRkJCQSiWxjK4OIqHlt1ucQFBQEb29v2X737t1Vfm/9h/xDDz2E4cOHo3Pnzti9ezdmzJjR6PkREREIDw+X7ZeWljJBEJGSyspKVFdXq3Suubm5Xicw1Lc2Sw5isVgr0+oKggBBEFBVVdXkORYWFnqfGZWIjEtlZSWsrKxUPl8qlSI7O7vdJgiDGq1UVFSEnJwc5OfnAwAyMjIA1P0jSKVSXLt2DTt27ICfnx8cHByQl5eHNWvWwMrKChMmTGjLqhORkVO1xVCvoKAA1dXV7TY5GFSHdHx8PIYMGYInnngCADB9+nQMGTIEGzduBFA3B31ycjImTJgAd3d3TJs2DWKxGMePH4ejo2NbVp2I2gmRSAQTE5NmtwdhvQiDajmEhIQgJCSkydednJxw4MAB/VWIiB44IpFIpQ9/QRD0UJu2Y1DJgYioranSMhAEAbW1tXqqUdtgciAikqNqcmjvmByIiOSoelupvWNyICKSw+RQx6BGKxERtbX65NDSpo4NGzZg0KBBsLW1ha2tLXx8fHDw4EHZ65WVlQgLC4O9vT1sbGwQHByMwsJCbV+aWpgciIjktDSMtX5TR48ePbBmzRqkpqbi9OnTGDduHCZNmoS//voLALBw4ULs27cPu3btQlJSEvLz8zF16lRdXJ7KeFuJiEiOLjqkJ06cqLC/evVqbNiwASkpKejRowc2bdqEbdu2Ydy4cQCALVu2oF+/fkhJScHw4cPVuwAteaCTw+XLl2FjY6PzOJ9++qnOY8hTZ56q1iooKNBbLAB48skn9RarV69eeoulb+p+820NfU5do43hpbruc6ipqcGuXbtQXl4OHx8fpKam4t69e/D19ZWd07dvX7i4uODEiRNMDkREhkCd5NBw+v/m5nE7f/48fHx8UFlZCRsbG+zevRv9+/dHeno6zM3N0alTJ4Xzu3btqvcvX/LY50BEJEedDmlnZ2eF5QAiIyObLNfDwwPp6ek4efIkXnnlFcyaNQsXL17U12WpjS0HIiI59XMrNaf+9lVubi5sbW1lx5u7hWZubg53d3cAgJeXF/788098/vnnmDZtGqqrq1FcXKzQeigsLIRUKm3FlbQOWw5ERHLUGa1UPzS1flOnf6W2thZVVVXw8vKCmZkZjhw5InstIyMDOTk58PHx0fr1qYotByIiOar0OajbYR0REYHAwEC4uLjg7t272LZtGxITE3Ho0CFIJBKEhoYiPDwcdnZ2sLW1xfz58+Hj49NmndEAkwMRkQJdJIebN2/i+eefx40bNyCRSDBo0CAcOnQI//nPfwDUjWg0MTFBcHAwqqqq4O/vjy+//FLja9AGJgciIjm6SA6bNm1q9nVLS0tERUUhKipKrXJ1icmBiEiOJk9At0dMDkREcpgc6hjUbyAuLg5+fn6wt7eHSCRCenq6wutFRUWYP38+PDw8YGVlBRcXFyxYsAAlJSVtU2Eiand0MfGeMTKo5FBeXo6RI0di7dq1jb6en5+P/Px8fPTRR7hw4QJiYmKQkJCA0NBQPdeUiNorJoc6BnVb6bnnngMAXL9+vdHXH374Yfz000+y/d69e2P16tV49tlncf/+fXToYFCXQ0RGSBcd0sbI6D9NS0pKYGtr22xiqKqqQlVVlWy/4XwoRET12OdQx6h/A7du3cLKlSsxd+7cZs+LjIxUmP/E2dlZTzUkImOji/UcjFGbXWFsbCxsbGxkW3JyslrvLy0txRNPPIH+/ftj2bJlzZ4bERGBkpIS2Zabm9uKmhNRe8Y+hzptdlspKCgI3t7esn111iC4e/cuAgICIBaLsXv3bpiZmTV7fnPT6BIRNfQgfPi3pM2Sg1gshlgsVvt9paWl8Pf3h4WFBeLj42FpaamD2hHRg0qV20bqrgRnjAyqQ7qoqAg5OTnIz88HUDczIQBIpVJIpVKUlpbCz88PFRUV+P7771FaWirrXHZwcICpqWmb1Z2I2geOVqpjUMkhPj4es2fPlu1Pnz4dALB06VIsW7YMZ86cwcmTJwFANi96vezsbPTs2VNvdSWi9snU1JRfNGFgySEkJAQhISFNvj527NgHojlHRG2HLYc6BpUciIjaGvsc6jA5EBHJYcuhDpMDEZEcY0kO//3vf1U+d8GCBWqXz+RARCTHWDqkP/30U4X9f/75BxUVFejUqRMAoLi4GNbW1nB0dNQoObT/Z8CJiNRgLNNnZGdny7bVq1fD09MTly5dQlFREYqKinDp0iU88sgjWLlypUblt/0VEhEZEGOcPuO9997D+vXr4eHhITvm4eGBTz/9FEuWLNGoTN5WIiKSIxKJWmwZ1NbW6qk2qrlx4wbu37+vdLympgaFhYUalflAJ4d///1XL/cWf/zxR53HkPf333/rLVaPHj30FgsAbGxs9BpPnwYPHqy3WPr8cNPnsE9txDKWDml548ePx0svvYRvvvkGjzzyCAAgNTUVr7zyCnx9fTUqk7eViIjk1HdIt7QZks2bN0MqlWLo0KGyiUYfffRRdO3aFd98841GZT7QLQciooaMseXg4OCAAwcO4MqVK7h8+TIAoG/fvujTp4/GZbLlQEQkRxejlSIjIzFs2DCIxWI4Ojpi8uTJsolF640dO1ap0/vll19WK07Pnj3h4eGBCRMmtCoxAEwOREQKdDFaKSkpCWFhYUhJScGvv/6Ke/fuwc/PD+Xl5Qrnvfjii7hx44ZsW7dunUrlV1RUIDQ0FNbW1hgwYABycnIAAPPnz8eaNWvUqms9JgciIjm6SA4JCQkICQnBgAEDMHjwYMTExCAnJwepqakK51lbW8uWKJBKpbC1tVWp/IiICJw9exaJiYkKa9z4+vpix44datW1HpMDEZEcdTqk69eUqd+qqqpUilFSUgIAsLOzUzgeGxuLLl264OGHH0ZERAQqKipUKm/Pnj344osvMHLkSIXENWDAAGRlZalURkPskCYikqNKn0L9687OzgrH69eeaU5tbS1ef/11jBgxAg8//LDs+DPPPANXV1c4OTnh3LlzePvtt5GRkYG4uLgW6/zPP//A0dFR6Xh5ebnGnedMDkREctQZrZSbm6tw60eVterDwsJw4cIFHDt2TOH43LlzZT8PHDgQ3bp1w/jx45GVlYXevXs3W+bQoUOxf/9+zJ8/X6F+33zzDXx8fFqsU2MMKjnExcVh48aNSE1NRVFREdLS0uDp6alwztdff41t27bhzJkzuHv3Lu7cuSObaIqIqLXUSQ62trYq9wsAwLx58/Dzzz/j6NGjLT5A6u3tDQDIzMxsMTl88MEHCAwMxMWLF3H//n18/vnnuHjxIo4fP46kpCSV6yfPoPocysvLMXLkSKxdu7bJcyoqKhAQEIB33nlHjzUjogdF/fQZzW3q3qoRBAHz5s3D7t278dtvv8HNza3F96SnpwMAunXr1uK5I0eORHp6Ou7fv4+BAwfil19+gaOjI06cOAEvLy+16lrPoFoOzz33HADg+vXrTZ7z+uuvAwASExN1XyEieuCo8gS0utOPhIWFYdu2bdi7dy/EYjEKCgoAABKJBFZWVsjKysK2bdswYcIE2Nvb49y5c1i4cCFGjx6NQYMGqRSjd+/eiI6OVqtezTGo5EBE1NZ08YT0hg0bANQ96CZvy5YtCAkJgbm5OQ4fPozPPvsM5eXlcHZ2RnBwsFozqtbW1iIzMxM3b95USl6jR49Wq77AA5IcqqqqFIaYlZaWtmFtiMiQqTNaSVUtTQjo7Oyscd8AAKSkpOCZZ57B//73P6VYIpEINTU1apfZZn0OsbGxsLGxkW3Jyck6ixUZGQmJRCLbGg4/IyKqZ4zrObz88ssYOnQoLly4gKKiIty5c0e2FRUVaVRmm7UcgoKCZL3xANC9e3edxYqIiEB4eLhsv7S0lAmCiBpljBPvXb16FT/++CPc3d21VmabJQexWAyxWKyXWPVT2BIRtcQYk4O3tzcyMzPbR3JoTFFREXJycpCfnw8AslkL6+cZAYCCggIUFBQgMzMTAHD+/HmIxWK4uLgoPYpORKQuXfQ56Nr8+fPxxhtvoKCgAAMHDoSZmZnC66qOeJJnUMkhPj4es2fPlu1Pnz4dgOIj6Rs3bsTy5ctl59T3wtf3+hMRtYYxthyCg4MBAHPmzJEdE4lEEARB4w5pg0oOISEhLX7AL1u2rMW5S4iINGWMySE7O1vrZRpUciAiamvGmBxcXV21XiaTAxGRHGPpc4iPj0dgYCDMzMwQHx/f7LlBQUFql8/kQEQkp35upZbOaWuTJ09GQUGBbNnRprSLPgciorZmLLeV5KfIUHeuJ1UwORARyTGW5KBrTA5ERHKMNTmUl5cjKSkJOTk5qK6uVnhtwYIFapfH5EBEJMcYk0NaWhomTJiAiooKlJeXw87ODrdu3YK1tTUcHR2ZHNQ1YMAAtVZx0lRLc8NrW05Ojt5i5eXl6S0WoN8/ylu3buktFgBYW1vrLZYmT8xqquG3WF26d+9eq8swltFK8hYuXIiJEydi48aNkEgkSElJgZmZGZ599lm89tprGpVpWFdIRNTGjHFW1vT0dLzxxhswMTGBqakpqqqq4OzsjHXr1mm8aiaTAxGRHGNMDmZmZrLWjKOjo+zugUQiQW5urkZlPtC3lYiIGjLGPochQ4bgzz//xEMPPYQxY8bg/fffx61bt/Ddd9/h4Ycf1qhMthyIiOQYY8vhgw8+QLdu3QAAq1evRufOnfHKK6/gn3/+wddff61RmWw5EBHJMcYO6aFDh8p+dnR0REJCQqvLZHIgIpJjjLeVdIHJgYhIjrEkhyFDhqhcjzNnzqhdPpMDEZEcY0kOzU22pw0GlRzi4uKwceNGpKamoqioCGlpafD09FQ4p7KyEm+88Qa2b9+Oqqoq+Pv748svv0TXrl3bptJE1K4YS3JYunSpTss3qF6V8vJyjBw5EmvXrm3ynIULF2Lfvn3YtWsXkpKSkJ+fj6lTp+qxlkTUntVP2d3cpm5yiIyMxLBhwyAWi2VTbGdkZCicU1lZibCwMNjb28PGxgbBwcEoLCxUK87p06fx3Xff4bvvvkNqaqpa723IoFoOzz33HADg+vXrjb5eUlKCTZs2Ydu2bRg3bhyAurWj+/Xrh5SUFAwfPlxfVSWidkoXLYekpCSEhYVh2LBhuH//Pt555x34+fnh4sWL6NixI4C6L7779+/Hrl27IJFIMG/ePEydOhV//PFHi+Xn5eVhxowZ+OOPP9CpUycAQHFxMR577DFs374dPXr0UKu+gIG1HFqSmpqKe/fuwdfXV3asb9++cHFxwYkTJ9qwZkTUXujiOYeEhASEhIRgwIABGDx4MGJiYpCTkyP7dl//xfeTTz7BuHHj4OXlhS1btuD48eNISUlpsfwXXngB9+7dw6VLl1BUVISioiJcunQJtbW1eOGFFzT6PRhUy6ElBQUFMDc3l2XGel27dkVBQUGT76uqqkJVVZVsv7S0VFdVJCIjp07LoeFniYWFBSwsLFqMUVJSAgCws7MD0PIX35buiiQlJeH48ePw8PCQHfPw8MD69esxatSoFuvTmDZrOcTGxsLGxka2JScn6yxWZGQkJBKJbHN2dtZZLCIybuq0HJydnRU+WyIjI1ssv7a2Fq+//jpGjBghm9pC0y++9ZydnRudkbampgZOTk4qXLWyNms5BAUFwdvbW7bfvXv3Ft8jlUpRXV2N4uJihV9iYWEhpFJpk++LiIhAeHi4bL+0tJQJgogaZWpq2uI0+/Wv5+bmKkz7r0qrISwsDBcuXMCxY8daV1E5H374IebPn4+oqCjZ09KnT5/Ga6+9ho8++kijMtssOYjFYojFYrXe4+XlBTMzMxw5cgTBwcEAgIyMDOTk5MDHx6fJ96na1CMiUue2kq2trVprwsybNw8///wzjh49qtBJrOkX33ohISGoqKiAt7c3OnSo+1i/f/8+OnTogDlz5mDOnDmyc4uKilSqq0H1ORQVFSEnJwf5+fkAIBvqJZVKIZVKIZFIEBoaivDwcNjZ2cHW1hbz58+Hj48PRyoRkVboYrSSIAiYP38+du/ejcTERLi5uSm8rukX33qfffaZWvVRhUElh/j4eMyePVu2P336dAB1D3ssW7YMAPDpp5/CxMQEwcHBCg/BERFpgy6SQ1hYGLZt24a9e/dCLBbL+hEkEgmsrKxa/cV31qxZatVHFQY1lDUkJASCICht9YkBACwtLREVFYWioiKUl5cjLi5OpWYXEZEqdDGUdcOGDSgpKcHYsWPRrVs32bZjxw7ZOZ9++imefPJJBAcHY/To0ZBKpYiLi1Op/JiYmEaP379/HxEREWrVtZ5BJQciorbW0tPRqkzp3VBjX3oFQUBISIjsnNZ88V2wYAGeeuop3LlzR3YsIyMD3t7e+OGHH9Sqaz0mByIiOca42E9aWhry8vIwcOBA/Prrr4iKisIjjzyCvn374uzZsxqVaVB9DkREbc1YJt6T17t3b/zxxx94/fXXERAQAFNTU3z77beYMWOGxmWq1XL4999/cezYMVy8eFHptcrKSmzdulXjihARGQJjbDkAwP79+7F9+3b4+PigU6dO2LRpk2zkpyZUTg5XrlxBv379MHr0aAwcOBBjxozBjRs3ZK+XlJQojDQiIjJGxpgcXnrpJTz11FN4++23kZycjHPnzsHc3BwDBw7Ezp07NSpT5eTw9ttv4+GHH8bNmzeRkZEBsViMESNGICcnR6PARESGyMTERPaUdFOboa0h/ccff+DkyZN44403IBKJIJVKceDAAaxYsULhATh1qNzncPz4cRw+fBhdunRBly5dsG/fPrz66qsYNWoUfv/9d9m0s0RExswY+xxSU1MbnQUiLCxMYTI/daic/v7991/ZY9lA3S9nw4YNmDhxIsaMGYMrV65oVAEiIkNijLeVLCwskJWVhSVLlmDGjBm4efMmAODgwYO4f/++RmWq3HLo27cvTp8+jX79+ikc/+KLLwDUTaRHjXvooYf0Gq+iokJvsa5du6a3WEDdeHF90fd8XC4uLnqL1blzZ73F0uckl5p+EMozxpZDUlISAgMDMWLECBw9ehSrV6+Go6Mjzp49i02bNuHHH39Uu0yVWw5Tpkxp8mGKL774AjNmzNDrHy4RkS4YY8th8eLFWLVqFX799VeYm5vLjo8bN06lxYIao3JyiIiIwIEDB5p8/csvv0Rtba1GlSAiMhS6eEJa186fP48pU6YoHXd0dMStW7c0KtOwrpCIqI0ZY8uhU6dOCo8W1EtLS1NprZzGMDkQEckxxuQwffp0vP322ygoKIBIJEJtbS3++OMPvPnmm3j++ec1KpPJgYhIjjHeVvrggw/Qt29fODs7o6ysDP3798fo0aPx2GOPYcmSJRqVybmViIjkGONoJXNzc0RHR+P999/H+fPnUVZWhiFDhrRqpCSTAxGRHFVaBobWcqjn7OystaHDGl3hd999hxEjRsDJyQn/+9//ANQtU7d3716tVIqIqK0YY5+DLqidHDZs2IDw8HBMmDABxcXFqKmpAVDXW66LdUwbKiwsREhICJycnGBtbY2AgABcvXpV53GJ6MHA5FBH7eSwfv16REdH491334Wpqans+NChQ3H+/HmtVq4hQRAwefJkXLt2DXv37kVaWhpcXV3h6+uL8vJyncYmogfHg54YAA36HLKzszFkyBCl4xYWFjr/gL569SpSUlJw4cIFDBgwAEBdS0YqleKHH37ACy+8oNP4RNT+GWOHtC6o3XJwc3NDenq60vGEhASleZe0raqqCkDdWqv1TExMYGFhgWPHjjX7vtLSUoWNiKgxxnpbKTk5Gc8++yx8fHzw999/A6jrH27us7E5aieH8PBwhIWFYceOHRAEAadOncLq1asRERGBt956S6NKqKpv375wcXFBREQE7ty5g+rqaqxduxZ5eXmNPh1YLzIyEhKJRLbpcyIwIjIuLa3lUL8Zkp9++gn+/v6wsrJCWlqa7It0SUkJPvjgA43KVDs5vPDCC1i7di2WLFmCiooKPPPMM9iwYQM+//xzTJ8+XaNKNCU2NhY2NjayLSUlBXFxcbhy5Qrs7OxgbW2N33//HYGBgc0OLYuIiEBJSYlsy83N1Wo9iaj9MMaWw6pVq7Bx40ZER0fDzMxMdnzEiBE4c+aMRmVqNJR15syZuHr1KsrKylBQUIC8vDyEhoZqVIHmBAUFIT09XbYNHToUXl5eSE9PR3FxMW7cuIGEhATcvn0bvXr1arIcCwsL2NraKmxERI3RRXI4evQoJk6cCCcnJ4hEIuzZs0fh9ZCQEKXyAwICVC4/IyMDo0ePVjoukUhQXFysVl3rqZ0cxo0bJwtmbW0NR0dHAEBpaSnGjRunUSWaIhaL4e7uLtusrKxkr0kkEjg4OODq1as4ffo0Jk2apNXYRPRg0kVyKC8vx+DBgxEVFdXkOQEBAbhx44Zsa2qJhMZIpVJkZmYqHT927FizX5ybo/ZopcTERFRXVysdr6ysRHJyskaVUMeuXbvg4OAAFxcXnD9/Hq+99homT54MPz8/nccmovZPF09IBwYGIjAwsNlzLCwsIJVK1Sq33osvvojXXnsNmzdvhkgkQn5+Pk6cOIE333wT7733nkZlqpwczp07J/v54sWLKCgokO3X1NQgISFB46lh1XHjxg2Eh4ejsLAQ3bp1w/PPP6/xxRMRNdRW02ckJibC0dERnTt3xrhx47Bq1SrY29ur9N7FixejtrYW48ePR0VFBUaPHg0LCwu8+eabmD9/vkb1UTk5eHp6yppTjd0+srKywvr16zWqhDoWLFiABQsW6DwOET2Y1HnOoeGweAsLC42Wlw0ICMDUqVPh5uaGrKwsvPPOOwgMDMSJEydUGhklEonw7rvvYtGiRcjMzJTNzGpjY6N2XeqpnByys7MhCAJ69eqFU6dOwcHBQfaaubk5HB0dDW54FxGRutRJDg2HxS9duhTLli1TO6b8SM+BAwdi0KBB6N27NxITEzF+/PgW319SUoKamhrY2dmhf//+suNFRUXo0KGDRoNwVE4Orq6uAMClQImoXVPntlJubq7CB68mrYbG9OrVC126dEFmZqZKyWH69OmYOHEiXn31VYXjO3fuRHx8fLNLPDdFpeQQHx+PwMBAmJmZIT4+vtlzg4KC1K4EEZGhUKfloKuh8Xl5ebh9+za6deum0vknT57EJ598onR87NixePfddzWqg0rJYfLkySgoKICjoyMmT57c5HkikUg2SysRkTESiUQtthzUHcpaVlamMNQ0Ozsb6enpsLOzg52dHZYvX47g4GBIpVJkZWXhrbfegru7O/z9/VUqv6qqCvfv31c6fu/ePfz7779q1bWeSl3utbW1sucZamtrm9yYGIjI2OlimdDTp09jyJAhsklLw8PDMWTIELz//vswNTXFuXPnEBQUhD59+iA0NBReXl5ITk5W+TbVo48+iq+//lrp+MaNG+Hl5aVWXetxJTgiIjm6mJV17NixEAShydcPHTqkVnkNrVq1Cr6+vjh79qysj+LIkSP4888/8csvv2hUpsrp78SJE/j5558Vjm3duhVubm5wdHTE3LlzZZM9EREZK2OcW2nEiBE4ceIEnJ2dsXPnTuzbtw/u7u44d+4cRo0apVGZKrccVqxYgbFjx+LJJ58EAJw/fx6hoaEICQlBv3798OGHH8LJyUmjYVxtRV//yObm5jqPIU+f69vq+9o0nSdGE/r+sqOL+cmacuvWLb3F6tKli95iVVZWtroMY13PwdPTE7GxsVorT+XkkJ6ejpUrV8r2t2/fDm9vb0RHRwOoG++r6RhfIiJDocqU3Ib4TFdtbS0yMzNx8+ZNpUcOGpuUryUqJ4c7d+6ga9eusv2kpCSFuUKGDRvGqbCJyOgZY8shJSUFzzzzDP73v/8p9W1oOopU5fsPXbt2RXZ2NgCguroaZ86cwfDhw2Wv3717V2EecSIiY2SMfQ4vv/wyhg4digsXLqCoqAh37tyRbUVFRRqVqXLLYcKECVi8eDHWrl2LPXv2wNraWqGj49y5c+jdu7dGlSAiMhTG2HK4evUqfvzxR7i7u2utTJVbDitXrkSHDh0wZswYREdHIzo6WqEzcvPmzZw2m4iMni6ec9A1b2/vRtdzaA2VWw5dunTB0aNHUVJSAhsbG6UOmV27drVqBkAiIkPQVlN2t8b8+fPxxhtvoKCgAAMHDlS6xT9o0CC1y1T7ITiJRNLocTs7O7WDExEZGmO8rRQcHAwAmDNnjuyYSCSCIAgad0jzCWkiIjnGmBzqBwtpk9Elh7KyMixevBh79uzB7du34ebmhgULFuDll19u66oRUTugi4n3dK1+SQVtMrrkEB4ejt9++w3ff/89evbsiV9++QWvvvoqnJycOF04EbWaMbYc6l28eBE5OTmorq5WOK7JZ6PRJYfjx49j1qxZGDt2LABg7ty5+Oqrr3Dq1CkmByJqNWPskL527RqmTJmC8+fPy/oagP+fxHT6EJyheOyxxxAfH4+///4bgiDg999/x5UrVziMloi0whiHsr722mtwc3PDzZs3YW1tjb/++gtHjx7F0KFDkZiYqFGZRtdyWL9+PebOnYsePXqgQ4cOMDExQXR0dLNzh1RVVSlMotZwUXAionrGeFvpxIkT+O2339ClSxdZ8ho5ciQiIyOxYMECpKWlqV2mYaW/BmJjY2FjYyPbkpOTsX79eqSkpCA+Ph6pqan4+OOPERYWhsOHDzdZTmRkJCQSiWxruCg4EZExq6mpgVgsBlD3TFp+fj6Auo7qjIwMjco06JZDUFAQvL29Zfvdu3fH+PHjsXv3bjzxxBMA6h7uSE9Px0cffQRfX99Gy4mIiEB4eLhsv7S0lAmCiBpljC2Hhx9+GGfPnoWbmxu8vb2xbt06mJub4+uvv0avXr00KtOgk4NYLJZlQ6DuQ/3evXtK9/tMTU2VpqiVZ2FhofJye0T0YDPGDuklS5agvLwcQN3aO08++SRGjRoFe3t77NixQ6MyDTo5NGRra4sxY8Zg0aJFsLKygqurK5KSkrB161Z88sknbV09ImoHjLHl4O/vL/vZ3d0dly9fRlFRETp37qxxXY0qOQB1iwxFRERg5syZKCoqgqurK1avXs2H4IhIK4wxOTSmtVMaGV1ykEql2LJlS1tXg4jaKWNJDlOnTkVMTAxsbW0xderUZs+Ni4tTu3yjSw5ERLpkLMlBIpHI6tHUhKitweRARCTHWJJD/R0UQRCwfPlyODg4wMrKSmvlG1aXOxFRG9PFMqFHjx7FxIkT4eTkBJFIhD179ii8LggC3n//fXTr1g1WVlbw9fXF1atXVSpbEAS4u7sjLy9PrTq1hMmBiEiOLpJDeXk5Bg8ejKioqEZfX7duHf773/9i48aNOHnyJDp27Ah/f39UVla2WLaJiQkeeugh3L59W606tViuVksjImoHtJkYACAwMBCrVq3ClClTlF4TBAGfffYZlixZgkmTJmHQoEHYunUr8vPzlVoYTVmzZg0WLVqECxcuqF23prDPgYhIjjp9Dg3nadPkgdvs7GwUFBQozPAgkUjg7e2NEydOYPr06S2W8fzzz6OiogKDBw+Gubm5Ut9DUVGRWnUCmByIiBSo84R0w2l4li5dimXLlqkVr6CgAADQtWtXheNdu3aVvdaSzz77TK2YqmByICLSUG5uLmxtbWX7bTVNz6xZs7ReJpNDO2RmZqa3WKampnqLBUBh6nVda7ialq5t2rRJb7F+/PFHvcWqX3jGWGKpc1vJ1tZWITloQiqVAgAKCwvRrVs32fHCwkJ4enqqXV5lZaXS/11N6sgOaSIiOboYrdQcNzc3SKVSHDlyRHastLQUJ0+ehI+Pj0pllJeXY968eXB0dETHjh3RuXNnhU0TTA5ERHJ0kRzKysqQnp6O9PR0AHWd0Onp6cjJyYFIJMLrr7+OVatWIT4+HufPn8fzzz8PJycnTJ48WaXy33rrLfz222/YsGEDLCws8M0332D58uVwcnLC1q1b1fwN1OFtJSIiObqYsvv06dN4/PHHZfv168vMmjULMTExeOutt1BeXo65c+eiuLgYI0eOREJCAiwtLVUqf9++fdi6dSvGjh2L2bNnY9SoUXB3d4erqytiY2Mxc+ZMteoLMDkQESnQxfQZY8eObbY/RCQSYcWKFVixYoVa5dYrKiqSLepja2srG7o6cuRIvPLKKxqVydtKRERy9N3noA29evVCdnY2AKBv377YuXMngLoWRadOnTQqk8mBiEiOMSaH2bNn4+zZswCAxYsXIyoqCpaWlli4cCEWLVqkUZm8rUREZOQWLlwo+9nX1xeXL19Gamoq3N3dMWjQII3KZHIgIpJjLFN2y8vNzVV4WtvV1RWurq6tKtPobis11cT78MMP27pqRNQO1I9WamkzJD179sSYMWMQHR2NO3fuaKVMw7pCFdy4cUNh27x5M0QiEYKDg9u6akTUDhhjn8Pp06fx6KOPYsWKFejWrRsmT56MH3/8sVUzChhdcpBKpQrb3r178fjjj8uGcRERtYYxJochQ4bgww8/RE5ODg4ePAgHBwfMnTsXXbt2xZw5czQq0+iSg7zCwkLs378foaGhzZ5XVVWF0tJShY2IqDHGmBzqiUQiPP7444iOjsbhw4fh5uaGb7/9VqOyjDo5fPvttxCLxZg6dWqz50VGRkIikci2htPsEhG1B3l5eVi3bh08PT3x6KOPwsbGpsnV51pi0MkhNjYWNjY2si05OVnh9c2bN2PmzJktPmIeERGBkpIS2Zabm6vLahOREROJRC12Rhtay+Grr77CmDFj0LNnT2zduhXTpk1DVlYWkpOT8fLLL2tUpkEPZQ0KCoK3t7dsv3v37rKfk5OTkZGRgR07drRYjiarMxERGYtVq1ZhxowZ+O9//4vBgwdrpUyDTg5isRhisbjR1zZt2gQvLy+t/SKIiADjfM6hfnZXbTLo20pNKS0txa5du/DCCy+0dVWIqJ0xxg5pXdTHoFsOTdm+fTsEQcCMGTPauipE1M4YY8tBF4yy5TB37lxUVFRAIpG0dVWIqJ0xxpaDLhhly4GISFfYcqjD5EBE1I7cunULJ0+eRE1NDYYNG4Zu3bppVA6TAxGRHGNuOfz0008IDQ1Fnz59cO/ePWRkZCAqKgqzZ89Wuyyj7HMgItIVY+pzKCsrU9hfvnw5Tp06hVOnTiEtLQ27du3Cu+++q1HZTA5ERHKMKTl4eXlh7969sv0OHTrg5s2bsv3CwkKYm5trVDZvKxERyTGm20qHDh1CWFgYYmJiEBUVhc8//xzTpk1DTU0N7t+/DxMTE8TExGhUNpMDEZEcY0oOPXv2xP79+/HDDz9gzJgxWLBgATIzM5GZmYmamhr07du3xbnnmvJAJwd9regkCILOY8jT5ypVpqameosF6P93qU/379/XW6z/+7//01us9PR0vcWqqKhodRm6SA7Lli3D8uXLFY55eHjg8uXLatevMTNmzEBgYCDefPNNjB07Fl9//TU8PT1bVeYDnRyIiPRlwIABOHz4sGy/QwftfPweOHAAly5dwuDBg/HNN98gKSkJM2fORGBgIFasWAErKyuNymWHNBGRHF2tId2hQweFVSy7dOnS6rq+8cYbmD17Nv7880+89NJLWLlyJcaMGYMzZ87A0tISQ4YMwcGDBzUqm8mBiEgPrl69CicnJ/Tq1QszZ85ETk5Oq8uMiYnBgQMHsH37dvz555/47rvvAADm5uZYuXIl4uLi8MEHH2hUNpMDEZEcdYayNlx+uKqqqtEyvb29ERMTg4SEBGzYsAHZ2dkYNWoU7t6926q6duzYEdnZ2QCA3Nxcpc7n/v37Ky2SpiomByIiOeokB2dnZ4UliCMjIxstMzAwEE899RQGDRoEf39/HDhwAMXFxdi5c2er6hoZGYnnn38eTk5OGDNmDFauXNmq8uSxQ5qISI46o5Vyc3Nha2srO67qipOdOnVCnz59kJmZqXlFAcycORMBAQG4du0aHnroIXTq1KlV5cljciAiakDVoaq2trYKyUFVZWVlyMrKwnPPPaf2exuyt7eHvb19q8tpiLeViIjk6GL6jDfffBNJSUm4fv06jh8/jilTpsDU1NSgFywzyuRw6dIlBAUFQSKRoGPHjhg2bJhWev6JiHQhLy8PM2bMgIeHB55++mnY29sjJSUFDg4ObV21JhndbaWsrCyMHDkSoaGhWL58OWxtbfHXX39p/Ig4EZE8XTwhvX379tZUqU0YXXJ49913MWHCBKxbt052rHfv3m1YIyJqT4xpbiVdMqrbSrW1tdi/fz/69OkDf39/ODo6wtvbG3v27GnrqhFRO2FMU3brklElh5s3b6KsrAxr1qxBQEAAfvnlF0yZMgVTp05FUlJSk++rqqpSeliFiIiaZtDJITY2FjY2NrItIyMDADBp0iQsXLgQnp6eWLx4MZ588kls3LixyXIiIyMVHlRxdnbW1yUQkZFhy6GOQfc5BAUFwdvbW7bv4OCADh06oH///grn9evXD8eOHWuynIiICISHh8v2S0tLmSCIiJph0MlBLBZDLBYrHBs2bJisBVHvypUrcHV1bbIcCwsLlZ9cJKIHGzuk6xh0cmjMokWLMG3aNIwePRqPP/44EhISsG/fPiQmJrZ11YioHWByqGPQfQ6NmTJlCjZu3Ih169Zh4MCB+Oabb/DTTz9h5MiRbV01IqJ2w+haDgAwZ84czJkzp62rQUTtEFsOdYyu5UBERLpnlC0HIiJdYcuhDpMDEZEcJoc6vK1ERERKmByIiEgJbysREcnhbaU6bDkQEZESthyIiOSw5VDngU4OJiYmMDHRfeNJEASdx5Cnz/+4tbW1eosFAGZmZnqLVVNTo7dYAPTyf7GePv/dPD099RZLG9PxMznU4W0lIiJS8kC3HIiIGmLLoQ5bDkREpIQtByIiOWw51GFyICKSw+RQh7eViIj0ICoqCj179oSlpSW8vb1x6tSptq5Ss5gciIjk1LccWtrUsWPHDoSHh2Pp0qU4c+YMBg8eDH9/f9y8eVNHV9F6TA5ERDr2ySef4MUXX8Ts2bPRv39/bNy4EdbW1ti8eXNbV61JRpccQkJClDJ4QEBAW1eLiNqJu3fvqrQBdQ/dyW9VVVVK5VVXVyM1NRW+vr6yYyYmJvD19cWJEyf0dl3qMsoO6YCAAGzZskW2b2Fh0Ya1IaL2wNzcHFKpFM7Oziqdb2Njo3Tu0qVLsWzZMoVjt27dQk1NDbp27apwvGvXrrh8+XKr6qxLRpkcLCwsIJVK27oaRNSOWFpaIjs7G9XV1SqdLwiCUt9De/qiapTJITExEY6OjujcuTPGjRuHVatWwd7evsnzq6qqFJp72ph/hYjaH0tLS1haWmq1zC5dusDU1BSFhYUKxwsLCw36S67R9TkEBARg69atOHLkCNauXYukpCQEBgY2O0laZGQkJBKJbFO12UhE1Frm5ubw8vLCkSNHZMdqa2tx5MgR+Pj4tGHNmicS9D1lqBpiY2Px0ksvyfYPHjyIUaNGKZxz7do19O7dG4cPH8b48eMbLaexloOzszOKi4tha2urm8rL4ays2qPPa9P3rKympqZ6jacv+pxttrS0FBKJBCUlJXr521bVjh07MGvWLHz11Vd49NFH8dlnn2Hnzp24fPmyUl+EoTDo20pBQUHw9vaW7Xfv3l3pnF69eqFLly7IzMxsMjlYWFi0q3uBRGRcpk2bhn/++Qfvv/8+CgoK4OnpiYSEBINNDICBJwexWAyxWNzsOXl5ebh9+za6deump1oREalv3rx5mDdvXltXQ2VG1edQVlaGRYsWISUlBdevX8eRI0cwadIkuLu7w9/fv62rR0TUbhhVcjA1NcW5c+cQFBSEPn36IDQ0FF5eXkhOTuZtIyIiLTLo20oNWVlZ4dChQ21dDSKids+oWg5ERKQfTA5ERKSEyYGIiJQwORARkRImByIiUsLkQERESoxqKKu21M91pK/ZWTm3kvZwbiXjo++5lQD9/821Rw9kcqhfxcnFxaWNa0JEunD37l1IJJK2roZRM+hZWXWltrYW+fn5EIvFan0TrZ/NNTc3V+czPrbXWPqOx2szzniaxhIEAXfv3oWTk5NeWyzt0QPZcjAxMUGPHj00fr+tra3epgNur7H0HY/XZpzxNInFFoN2MLUSEZESJgciIlLC5KAGCwsLLF26VC8zwLbXWPqOx2szznj6vjZS9kB2SBMRUfPYciAiIiVMDkREpITJgYiIlDA5NCEuLg5+fn6wt7eHSCRCenq60jkvvfQSevfuDSsrKzg4OGDSpEm4fPmy1mMVFRVh/vz58PDwgJWVFVxcXLBgwQKUlJTo7Nq+/vprjB07Fra2thCJRCguLtZZrMrKSoSFhcHe3h42NjYIDg5GYWGhRvEaKiwsREhICJycnGBtbY2AgABcvXpVK2U3VFZWhnnz5qFHjx6wsrJC//79sXHjRp3EAuqmEmls+/DDD3US79KlSwgKCoJEIkHHjh0xbNgw5OTk6CRWSEiI0nUFBAToJBY1jsmhCeXl5Rg5ciTWrl3b5DleXl7YsmULLl26hEOHDkEQBPj5+ak9J09LsfLz85Gfn4+PPvoIFy5cQExMDBISEhAaGqpWHFXjAUBFRQUCAgLwzjvvaBRDnVgLFy7Evn37sGvXLiQlJSE/Px9Tp05tVVyg7mnZyZMn49q1a9i7dy/S0tLg6uoKX19flJeXt7r8hsLDw5GQkIDvv/8ely5dwuuvv4558+YhPj5e67EA4MaNGwrb5s2bIRKJEBwcrPVYWVlZGDlyJPr27YvExEScO3cO7733HiwtLbUeq15AQIDC9f3www86i0WNEKhZ2dnZAgAhLS2txXPPnj0rABAyMzN1Hmvnzp2Cubm5cO/ePY1iqRrv999/FwAId+7c0ThOc7GKi4sFMzMzYdeuXbJjly5dEgAIJ06caFXMjIwMAYBw4cIF2bGamhrBwcFBiI6OblXZjRkwYICwYsUKhWOPPPKI8O6772o9VmMmTZokjBs3TidlT5s2TXj22Wd1UnZjZs2aJUyaNElv8UgZWw5aUl5eji1btsDNzQ3Ozs46j1dSUgJbW1t06GDcM6Ckpqbi3r178PX1lR3r27cvXFxccOLEiVaVXVVVBQAK325NTExgYWGBY8eOtarsxjz22GOIj4/H33//DUEQ8Pvvv+PKlSvw8/PTeqyGCgsLsX//fo1bk82pra3F/v370adPH/j7+8PR0RHe3t7Ys2eP1mPJS0xMhKOjIzw8PPDKK6/g9u3bOo1HipgcWunLL7+EjY0NbGxscPDgQfz6668wNzfXacxbt25h5cqVmDt3rk7j6ENBQQHMzc3RqVMnheNdu3ZFQUFBq8quTzIRERG4c+cOqqursXbtWuTl5eHGjRutKrsx69evR//+/dGjRw+Ym5sjICAAUVFRGD16tNZjNfTtt99CLBZr5XZcQzdv3kRZWRnWrFmDgIAA/PLLL5gyZQqmTp2KpKQkrccD6m4pbd26FUeOHMHatWuRlJSEwMBAvU+j/iBjcgAQGxsr+4C3sbFBcnKyyu+dOXMm0tLSkJSUhD59+uDpp59GZWWlTmIBdbNVPvHEE+jfvz+WLVvW4vmtjacOfcZSJX5KSgri4uJw5coV2NnZwdraGr///jsCAwNbPWNnY9e6fv16pKSkID4+Hqmpqfj4448RFhaGw4cPa/3aGv5uN2/ejJkzZ2qlD6BhrIyMDADApEmTsHDhQnh6emLx4sV48skntdLh3ti1TZ8+HUFBQRg4cCAmT56Mn3/+GX/++ScSExNbHY9UY9z3JLQkKCgI3t7esv3u3bur/F6JRAKJRIKHHnoIw4cPR+fOnbF7927MmDFD67Hu3r2LgIAAiMVi7N69G2ZmZi2+pzXx1KVJLKlUiurqahQXFyu0HgoLCyGVSlsd38rKCunp6SgpKUF1dTUcHBzg7e2NoUOHqlW2KrHGjx+P3bt344knngAADBo0COnp6fjoo48UbptpK1695ORkZGRkYMeOHa2K0VQsBwcHdOjQAf3791c4r1+/flq5PafK/5tevXqhS5cuyMzMxPjx41sdk1rG5ABALBZDLBa3uhxBECAIguxetzZjlZaWwt/fHxYWFoiPj1f5G6K2rk1Xsby8vGBmZoYjR47IRtlkZGQgJycHPj4+WotfP43z1atXcfr0aaxcuVKtsluKVVpainv37im1SExNTbWyWl5z17Zp0yZ4eXlh8ODBrY7TVKxhw4bJWhD1rly5AldXV53EaygvLw+3b99Gt27dWh2PVMPk0ISioiLk5OQgPz8fAGR/GFKpFFKpFNeuXcOOHTvg5+cHBwcH5OXlYc2aNbCyssKECRO0Gqu0tBR+fn6oqKjA999/j9LSUtlyiA4ODmovL9lSPKCuL6CgoACZmZkAgPPnz0MsFsPFxQV2dnZaiyWRSBAaGorw8HDY2dnB1tYW8+fPh4+PD4YPH67WdTVm165dcHBwgIuLC86fP4/XXnsNkydP1nonsa2tLcaMGYNFixbBysoKrq6uSEpKwtatW/HJJ59oNZa80tJS7Nq1Cx9//LHOYgDAokWLMG3aNIwePRqPP/44EhISsG/fPp3c5ikrK8Py5csRHBwMqVSKrKwsvPXWW3B3d4e/v7/W41ET2ni0lMHasmWLAEBpW7p0qSAIgvD3338LgYGBgqOjo2BmZib06NFDeOaZZ4TLly9rPVb9cNLGtuzsbK3HEwRBWLp0aaPnbNmyReux/v33X+HVV18VOnfuLFhbWwtTpkwRbty4ofZ1Nebzzz8XevToIZiZmQkuLi7CkiVLhKqqKq2U3dCNGzeEkJAQwcnJSbC0tBQ8PDyEjz/+WKitrdVJPEEQhK+++kqwsrISiouLdRaj3qZNmwR3d3fB0tJSGDx4sLBnzx6dxKmoqBD8/PwEBwcHwczMTHB1dRVefPFFoaCgQCfxqHGclZWIiJRwtBIRESlhciAiIiVMDkREpITJgYiIlDA5EBGREiYHIiJSwuRARERKmByIiEgJkwMZLZFIpPM1BYgeVEwOZJD++ecfvPLKK3BxcYGFhQWkUin8/f3xxx9/yM65ceMGAgMDAQDXr19vco1qdR09ehQTJ06Ek5MTExA9sDjxHhmk4OBgVFdX49tvv0WvXr1QWFiII0eOKKwGpu6U3qoqLy/H4MGDMWfOHJ0snkNkFNp6cieihu7cuSMAEBITE5s9D4Cwe/du2c/y25gxY2TnRUdHC3379hUsLCwEDw8PISoqSuW6yMcgepCw5UAGp35FsD179mD48OGwsLBo8T2nTp3Co48+isOHD2PAgAGypVpjY2Px/vvv44svvsCQIUOQlpaGF198ER07dsSsWbN0fSlERot9DmRwOnTogJiYGHz77bfo1KkTRowYgXfeeQfnzp1r8j0ODg4AAHt7e0ilUtmaE0uXLsXHH3+MqVOnws3NDVOnTsXChQvx1Vdf6eVaiIwVkwMZpODgYOTn5yM+Ph4BAQFITEzEI488gpiYGJXLKC8vR1ZWFkJDQxXWKF61ahWysrJ0V3midoC3lchgWVpa4j//+Q/+85//4L333sMLL7yApUuXIiQkRKX3l5WVAQCio6MV1igGoPbqeUQPGiYHMhr9+/dvclhpfR9DTU2N7FjXrl3h5OSEa9euYebMmfqoIlG7weRABuf27dt46qmnMGfOHAwaNAhisRinT5/GunXrMGnSpEbf4+joCCsrKyQkJKBHjx6wtLSERCLB8uXLsWDBAkgkEgQEBKCqqgqnT5/GnTt3EB4e3mhZZWVlsrWzASA7Oxvp6emws7ODi4uLTq6ZyOC09XApooYqKyuFxYsXC4888oggkUgEa2trwcPDQ1iyZIlQUVEhOw8NhplGR0cLzs7OgomJicJQ1tjYWMHT01MwNzcXOnfuLIwePVqIi4trMn5Ta3bPmjVLB1dLZJi4hjQRESnhaCUiIlLC5EBEREqYHIiISAmTAxERKWFyICIiJUwORESkhMmBiIiUMDkQEZESJgciIlLC5EBEREqYHIiISAmTAxERKfl/TQkP99cJRgIAAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 1, figsize=(4, 4))\n",
"im = axes.imshow(\n",
" m,\n",
" cmap=\"Greys\",\n",
" vmin=0,\n",
")\n",
"axes.set(\n",
" xticks=ticks,\n",
" xlabel=\"Site 1\",\n",
" yticks=ticks,\n",
" xticklabels=positions,\n",
" yticklabels=positions,\n",
" ylabel=\"Site 2\",\n",
" aspect=\"equal\",\n",
")\n",
"plt.colorbar(im, shrink=0.6, label=\"% variance explained\")\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "9fd01c5a",
"metadata": {},
"source": [
"Overall, this summary plots show that positions in the Shine-Dalgarno sequence tend to interact mostly with adjacent positions. Moreover, higher order interactions can involve sites that are further appart from each other in the primary sequence."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "gpmap",
"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.8.20"
}
},
"nbformat": 4,
"nbformat_minor": 5
}