"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"sns.heatmap(assignments)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where each row corresponds to a cell, and each column to a cell type, with the entry being the probability of that cell belonging to a particular cell type.\n",
"\n",
"To fetch an array corresponding to the most likely cell type assignments, call"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
cell_type
\n",
"
\n",
" \n",
" \n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_1
\n",
"
Unknown
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_2
\n",
"
Unknown
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_3
\n",
"
Unknown
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4
\n",
"
Unknown
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_5
\n",
"
Unknown
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4927
\n",
"
epithelial(basal)
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4928
\n",
"
Unknown
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4929
\n",
"
epithelial(basal)
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4930
\n",
"
epithelial(luminal)
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4931
\n",
"
epithelial(luminal)
\n",
"
\n",
" \n",
"
\n",
"
4931 rows × 1 columns
\n",
"
"
],
"text/plain": [
" cell_type\n",
"BaselTMA_SP43_115_X4Y8_1 Unknown\n",
"BaselTMA_SP43_115_X4Y8_2 Unknown\n",
"BaselTMA_SP43_115_X4Y8_3 Unknown\n",
"BaselTMA_SP43_115_X4Y8_4 Unknown\n",
"BaselTMA_SP43_115_X4Y8_5 Unknown\n",
"... ...\n",
"BaselTMA_SP43_115_X4Y8_4927 epithelial(basal)\n",
"BaselTMA_SP43_115_X4Y8_4928 Unknown\n",
"BaselTMA_SP43_115_X4Y8_4929 epithelial(basal)\n",
"BaselTMA_SP43_115_X4Y8_4930 epithelial(luminal)\n",
"BaselTMA_SP43_115_X4Y8_4931 epithelial(luminal)\n",
"\n",
"[4931 rows x 1 columns]"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ast.get_celltypes()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cell type diagnostics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is important to run diagnostics to ensure that cell types express their markers at higher levels than other cell types. To do this, run the `diagnostics_celltype()` function, which will alert to any issues if a cell type doesn't express its marker signficantly higher than an alternative cell type (for which that protein isn't a marker):"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
feature
\n",
"
should be expressed higher in
\n",
"
than
\n",
"
mean cell type 1
\n",
"
mean cell type 2
\n",
"
p-value
\n",
"
note
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
Fibronectin
\n",
"
stromal
\n",
"
B cells
\n",
"
2.045195
\n",
"
1.565420
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
1
\n",
"
Vimentin
\n",
"
stromal
\n",
"
B cells
\n",
"
2.994335
\n",
"
1.095477
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
2
\n",
"
CD20
\n",
"
B cells
\n",
"
stromal
\n",
"
0.403677
\n",
"
0.081944
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
3
\n",
"
CD20
\n",
"
B cells
\n",
"
T cells
\n",
"
0.403677
\n",
"
0.123195
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
4
\n",
"
CD20
\n",
"
B cells
\n",
"
macrophage
\n",
"
0.403677
\n",
"
0.212008
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
5
\n",
"
CD20
\n",
"
B cells
\n",
"
epithelial(basal)
\n",
"
0.403677
\n",
"
0.118228
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
6
\n",
"
CD20
\n",
"
B cells
\n",
"
epithelial(luminal)
\n",
"
0.403677
\n",
"
0.131991
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
7
\n",
"
CD20
\n",
"
B cells
\n",
"
Other
\n",
"
0.403677
\n",
"
0.017356
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
8
\n",
"
CD45
\n",
"
B cells
\n",
"
stromal
\n",
"
0.077596
\n",
"
0.227257
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
"
\n",
"
9
\n",
"
CD45
\n",
"
B cells
\n",
"
epithelial(basal)
\n",
"
0.077596
\n",
"
0.182078
\n",
"
inf
\n",
"
Only 1 cell in a type: comparison not possible
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" feature should be expressed higher in than \\\n",
"0 Fibronectin stromal B cells \n",
"1 Vimentin stromal B cells \n",
"2 CD20 B cells stromal \n",
"3 CD20 B cells T cells \n",
"4 CD20 B cells macrophage \n",
"5 CD20 B cells epithelial(basal) \n",
"6 CD20 B cells epithelial(luminal) \n",
"7 CD20 B cells Other \n",
"8 CD45 B cells stromal \n",
"9 CD45 B cells epithelial(basal) \n",
"\n",
" mean cell type 1 mean cell type 2 p-value \\\n",
"0 2.045195 1.565420 inf \n",
"1 2.994335 1.095477 inf \n",
"2 0.403677 0.081944 inf \n",
"3 0.403677 0.123195 inf \n",
"4 0.403677 0.212008 inf \n",
"5 0.403677 0.118228 inf \n",
"6 0.403677 0.131991 inf \n",
"7 0.403677 0.017356 inf \n",
"8 0.077596 0.227257 inf \n",
"9 0.077596 0.182078 inf \n",
"\n",
" note \n",
"0 Only 1 cell in a type: comparison not possible \n",
"1 Only 1 cell in a type: comparison not possible \n",
"2 Only 1 cell in a type: comparison not possible \n",
"3 Only 1 cell in a type: comparison not possible \n",
"4 Only 1 cell in a type: comparison not possible \n",
"5 Only 1 cell in a type: comparison not possible \n",
"6 Only 1 cell in a type: comparison not possible \n",
"7 Only 1 cell in a type: comparison not possible \n",
"8 Only 1 cell in a type: comparison not possible \n",
"9 Only 1 cell in a type: comparison not possible "
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ast.diagnostics_celltype().head(n=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
".. note:: \n",
" In this tutorial, we end up with many \"Only 1 cell in a type: comparison not possible\" notes - this is simply because the small dataset size results in only a single cell assigned to many types, making statistical testing infeasible."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calling `ast.diagnostics_celltype()` returns a `pd.DataFrame`, where each column corresponds to a particular protein and two cell types, with a warning if the protein is not expressed at higher levels in the cell type for which it is a marker than the cell type for which it is not."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The diagnostics:\n",
"\n",
"1. Iterates through every cell type and every marker for that cell type\n",
"\n",
"2. Given a cell type *c* and marker *g*, find the set of cell types *D* that don't have *g* as a marker\n",
"\n",
"3. For each cell type *d* in *D*, perform a t-test between the expression of marker *g* in *c* vs *d*\n",
"\n",
"4. If *g* is not expressed significantly higher (at significance *alpha*), output a diagnostic explaining this for further investigation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If multiple issues are found, the markers and cell types may need refined."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Fitting cell state "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly as before, to fit cell state, call"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jinelles.h/Documents/Camlab/astir-top-level/astir/astir/models/cellstate.py:176: UserWarning: Delta loss batch size is greater than the number of epochs\n",
" warnings.warn(\"Delta loss batch size is greater than the number of epochs\")\n",
"training restart 1/5: 100%|██████████| 5/5 [ 3.41s/epochs, current loss: 59.8]\n",
"training restart 2/5: 100%|██████████| 5/5 [ 3.40s/epochs, current loss: 98.0] \n",
"training restart 3/5: 100%|██████████| 5/5 [ 3.51s/epochs, current loss: 78.2]\n",
"training restart 4/5: 100%|██████████| 5/5 [ 3.48s/epochs, current loss: 88.0] \n",
"training restart 5/5: 100%|██████████| 5/5 [ 3.37s/epochs, current loss: 60.0]\n",
"training restart (final): 100%|██████████| 50/50 [ 3.60s/epochs, current loss: 14.5]\n",
"/Users/jinelles.h/Documents/Camlab/astir-top-level/astir/astir/astir.py:268: UserWarning: Maximum epochs reached. More iteration may be needed to complete the training.\n",
" warnings.warn(msg)\n"
]
}
],
"source": [
"ast.fit_state(batch_size = 1024, learning_rate=1e-3, max_epochs=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and similary plot the losses via"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 0, 'Epoch')"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAEGCAYAAAD2TVeiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU9b3H8fd3JitkIxBCQgj7FpFFAiLgAq7VulStgkup15ZWbbW91tZu2v1qe+vWqlXrQuuCK1dcqkUUBRf2fY3sECAsSSBkT373jwwWdYAAmZxZPq/nmWdmzkyST3iGT845v3N+x5xziIjI5/m8DiAiEo5UjiIiQagcRUSCUDmKiAShchQRCSLO6wDN0aFDB9etWzevY4hIlJk/f/4u51xWsNciohy7devGvHnzvI4hIlHGzDYe6jVtVouIBKFyFBEJQuUoIhKEylFEJAiVo4hIECpHEZEgVI4iIkFExHGOR2NbeRXPfLKJnIwkcjOSyU1PJicjibSkeK+jiUgEibpyXL9rPw+/v5aGxs/PUzkwL51XbxqFmXmUTEQiSdSV48ieHVj92/PYWVFDcVkVxWXV/HvFDl5bXEzJvhqy05K8jigiESDqyhEgzu8jJz2ZnPRkhnaFjqmJvLa4mBXb9qocRaRZYmJApl9OGgArivd6nEREIkVMlGN6cjx57ZJZuU3lKCLNExPlCFCQk8YKlaOINFPMlGP/nDTW79pPZW2911FEJALETDkW5KbhHKzevs/rKCISAWKnHAODMiu3qRxF5Mhiphzz2iWTmhjHim3lXkcRkQgQM+VoZvTPTdPhPCLSLDFTjtC0ab1q+z4av3BqoYjIF8VcOVbWNrBxT6XXUUQkzMVUOfb/bFBGm9YicngxVY69s1Pw+0z7HUXkiGKqHJPi/fTKSgm65li6v5ZrH5/N8mKNZotIjJUjQP+c1KCnET710QZmFu3irn+t8iCViISbmCvHgtw0tpVXU7q/9rNllbX1TPp4A6mJccws2sXcDXu8CygiYSHmyjHYoMwLczdTVlnHw9cMpUNKIvdOW+NVPBEJEzFbjgc2resaGnls5noKu7ZjdO8OfPf0Hny0djez1+32MqaIeCzmyrFDSiLZaYmfleMbS7axtayK757eE4BrRnQlKzWRe9/R2qNILIu5coSmtccVxXtxzvG399fSu2MKY/t1BJpGtG88oyefrNvDR2t3eZxURLwS0nI0swwze8nMVpnZSjM7xcwyzWyamRUF7tuFMkMwBTlprN1ZwbQVO1i1fR8TT+uBz/efqxKOH55Pdloi900rwjmdaigSi0K95ng/8JZzrh8wCFgJ3A5Md871BqYHnreq/jlp1DU4fjV1OTnpSVw8uPPnXm9ae+zFnA17+Git9j2KxKKQlaOZpQGnAY8DOOdqnXNlwMXApMDbJgGXhCrDoRTkNg3KFJdXc/3o7iTEffmf4cphXchJT+KeaWu09igSg0K55tgD2Ak8aWYLzezvZtYWyHbObQMI3HcM9sVmNtHM5pnZvJ07d7ZosG7t25Ic7yctKY5xw/ODvicp3s9NY3oxf2MpM1a37M8XkfAXynKMA04CHnbODQH2cxSb0M65R51zhc65wqysrBYN5vcZ3zq1Oz+/oD8piYe+dPcVhV3Iz2zDH99erWnORGJMKMtxC7DFOTc78Pwlmspyh5nlAATuS0KY4ZBuPacvVw4LvtZ4QEKcj1vP6cPKbXt5bUlxKyUTkXAQsnJ0zm0HNptZ38CiM4EVwFRgQmDZBODVUGVoCRcOzKVfp1TumbaGuoZGr+OISCsJ9Wj194FnzGwJMBj4A3AXcLaZFQFnB56HLZ/P+PF5fdm4u5IX5m3+0uvlVXX84c2VrN1Z4UE6EQmVQ+9wawHOuUVAYZCXzgzlz21pY/p2pLBrO+5/p4hLh+SRnOAHYPOeSq57ai6fllSwtbSKB68+yeOkItJSYvIMmaNlZvzkK/0o2VfDpI83ALBgUymXPPghJXurOb1PFm8v307JvmpPc4pIy1E5NtOwbpmM6ZvFwzPWMnnOJsY/+gkpSXFMuWkUd1xYQH2j48V5W7yOKSItROV4FG47tx/lVXXc/spSTuyczpQbR9EzK4WeWSmM6JHJ5LmbdMiPSJRQOR6Fgtw0vj+2F9eO6MrT3zqZzLYJn7121cld2bynipmfarIKkWgQ0gGZaHTrOX2DLj/3hGzat03g2dkbOb1Pyx60LiKtT2uOLSQxzs/lhXm8s7KEHXs1MCMS6VSOLWj8sHwaGh3Pz/3y8ZAiEllUji2oW4e2jO7VgclzNtGggRmRiKZybGFXnZxPcXk176/x5JRxEWkhKscWdnZBNlmpiTw7e5PXUUTkOKgcW1i838cVhXm8u6qE4rIqr+OIyDFSOYbAuGH5ONDAjEgEUzmGQJfMNpzWO4vn526mXtOciUQklWOIXHVyPtv3VvOeLrEgEpFUjiFyZr+OdExN5NnZG72OIiLHQOUYInF+H1cO68KMNTvZUlrpdRwROUoqxxC6clgXQAMzIpFI5RhCee3acEafpoEZXX9GJLKoHEPsqpO7UrKvhukrdcaMSCRROYbYmL5ZdEpL4tk5nz9jZldFDQ/N+JTFm8s8SiYih6P5HEPswMDMA+8WsXlPJX6f8egH65g8dxPVdY0kxhVx/7jBnDcgx+uoInIQrTm2gnHDu2DA9ZPmcvqf3uPpTzby1YG5vHzDKRTkpnHDMwt4YtZ6r2OKyEG05tgKctKTOW9AJ6avLOGq4fl8+7Qe5LVrA8Bz3x7BLZMX8pvXV7CltIqfX9Afv888TiwiKsdWcs8Vg6lraCQ1Kf5zy5Pi/Tx09VB+/8ZKnvhwPcVlVdw/fjCJcX6PkooIaLO61STF+79UjAf4fcYdFxbwy68W8Nby7fxq6opWTiciX6Q1xzBy/eju7Kqo4eEZaxmUl8644fleRxKJWVpzDDM/Oqcvp/buwB2vLmeRDvMR8YzKMcz4fcYD44bQMS2RG56ez66KGq8jicQklWMYatc2gb9dM5Q9+2u56ZkFmhNSxAMhLUcz22BmS81skZnNCyzLNLNpZlYUuG8XygyRakDndO667ERmr9/DL19dxv6aeq8jicSU1lhzHOOcG+ycKww8vx2Y7pzrDUwPPJcgvjYkj4mn9eC5OZsZdfe73PfOGsoqa72OJRITzLnQXV/ZzDYAhc65XQctWw2c4ZzbZmY5wAznXN/DfZ/CwkI3b968kOUMdws3lfLQjLVMW7GDNgl+rhnRletGdSMnPdnraCIRzczmH7Ti9vnXQlyO64FSwAGPOOceNbMy51zGQe8pdc59adPazCYCEwHy8/OHbtyoGbVXb9/HwzM+ZeriYsyMM/t15OoRXTm1Vwd8OqtG5Kh5WY65zrliM+sITAO+D0xtTjkeLNbXHL9o855Knpm9iRfnbWb3/lryM9tw1cn5fGt0d+L8GmMTaa7DlWNI/yc554oD9yXAFGA4sCOwOU3gXhMdHqUumW24/Sv9+OinY3lg/BA6pSdx179W8fQnWrsWaSkhK0cza2tmqQceA+cAy4CpwITA2yYAr4YqQ7RLjPNz0aBcXvjOKQzMS+e5OZsJ5ZaASCwJ5ZpjNjDLzBYDc4A3nHNvAXcBZ5tZEXB24Lkcp/HD81m9Yx8LNpV6HUUkKoTs3Grn3DpgUJDlu4EzQ/VzY9VFg3L53esreHb2ZoZ2zfQ6jkjE0977KNE2MY6Lh3Tm9SXFlFfWeR1HJOKpHKPIVcPzqalvZMrCLV5HEYl4KscoMqBzOgPz0nl2ziYNzIgcJ5VjlLlqeD5rdlRoYEbkOKkco8yFg3Jpm+Dn2dmbvY4iEtFUjlFGAzMiLUPlGIUOHphxzrGltJI3l27j7rdW8e6qHV7HE4kIuoZMFDowMHPvO0X85d1P2b3/P9OcvdA2gQ9vH0tSvK5uKHI4KscodfPY3tw/vYi+nVIZlJfOwLwMSitr+eaTc5mycCvjdfEukcNSOUapswqyOasg+3PLnHMM6JzGYzPXcWVhF01zJnIY2ucYQ8yMb5/ag3U79/Peak2GJHI4KscYc/6JOeSkJ/HYzHVeRxEJayrHGBPv9/Ffo7rzybo9LN1S7nUckbClcoxBVw7vQkpinNYeRQ5D5RiD0pLiGTesC28s3cbWsiqv44iEJZVjjLpudHcAnpy13uMkIuFJ5RijOmckc8GJOUyeu5nlxeVs3L2fbeVV7NlfS219o9fxRDyn4xxj2LdP7cHUxcVc8MCszy3PbJvAKzeMpFuHth4lE/GeyjGGnZiXzkvfPYXte6upqWuktqGR6roG/vzvNfzqteU8+c1hmOlAcYlNKscYV9jty9ebaXTw29dX8O8VOzj3hE4epBLxnvY5ypdMOKUr/Tql8pvXVlBV2+B1HBFPqBzlS+L8Pn5z8QC2llXx1/eKvI4j4olmlaOZ9TSzxMDjM8zsZjPLCG008dLw7plcelJnHv1gHet2VngdR6TVNXfN8WWgwcx6AY8D3YFnQ5ZKwsJPv9KfpDg/d05drgt2Scxpbjk2Oufqga8B9znnfgjkhC6WhIOs1ERuPacPM4t28day7V7HEWlVzS3HOjMbD0wAXg8siw9NJAkn14zoSkFOGr98dRnFOtVQYkhzy/E64BTg98659WbWHXg6dLEkXMT5fTwwfjDVdY1M/Oc8jV5LzGhWOTrnVjjnbnbOPWdm7YBU59xdIc4mYaJXx1QeGD+Y5cV7+fHLS7T/UWJCc0erZ5hZmpllAouBJ83sntBGk3Aytl82t53bl9cWF/Pw+2u9jiMScs3drE53zu0FLgWedM4NBc5qzheamd/MFprZ64Hn3c1stpkVmdnzZpZwbNGltd1wek8uHJTLn95erUu8StRrbjnGmVkOcAX/GZBprluAlQc9vxu41znXGygFrj/K7yceMTP+eNlATshN45bnFjF1cTG7K2q8jiUSEs0tx98AbwNrnXNzzawHcMRTJ8wsD7gA+HvguQFjgZcCb5kEXHK0ocU7yQl+Hr22kNSkOG5+biFDf/cO5977Ab+aupzpK3dof6REDQvlh9nMXgL+B0gFfgR8E/jEOdcr8HoX4F/OuQFBvnYiMBEgPz9/6MaNG0OWU45eXUMjS7aU88m63XyybjdzN+yhuq6RC07M4U9fH0ibBM1pIuHPzOY75wqDvdbcAZk8M5tiZiVmtsPMXg6sFR7ua74KlDjn5h+8OMhbg7azc+5R51yhc64wKyurOTGlFcX7fQzt2o6bxvTin9efzJI7z+X2r/TjzWXbuPShj9i8p9LriCLHpbmb1U8CU4FcoDPwWmDZ4YwCLjKzDcBkmjan7wMyzOzAakUeUHyUmSUMJcT5+O7pPXnym8PYWlbFRX+dxUdrd3kdS+SYNbccs5xzTzrn6gO3p4DDrs45537qnMtzznUDxgHvOueuBt4DLg+8bQLw6rFFl3B0Rt+OTP3eaNqnJHLt43N4fNZ67YeUiNTcctxlZtcEDsvxm9k1wO5j/Jk/Af7bzD4F2tM0kYVEke4d2jLlxpGM6duR376+gm88MYft5dVexxI5Ks0akDGzfOCvNJ1C6ICPgJudc5tCG69JYWGhmzdvXmv8KGlBzjmenr2JP7yxkoQ4H7+9ZAAXDcr1OpbIZ457QMY5t8k5d5FzLss519E5dwlNB4SLHJKZce2Irrx5y6n0yGrLzc8t5HvPLqC8ss7raCJHdDwzgf93i6WQqNa9Q1te/M4p3HZuX95atp1rn5itCSwk7B1POeqydNJscX4fN43pxcPXDGXp1nJ+9NJiDdRIWDuectQnW47a2QXZ3H5eP95Yso373tH1aSR8HfY0BjPbR/ASNCA5JIkk6k08rQdFJRXcP72Inh1TNEgjYemw5eicS22tIBI7zIzff20Am3ZXctuLi8nPbMPgLrpem4QXXZpVPJEY5+fha06iY1oi3/7HPNbv2u91JJHPUTmKZ9qnJPL4hGE0NDouefBDZhXpdEMJHypH8VSf7FRevWkUndKSmPDkHP7x8QaNYktYUDmK57pktuHlG0cypm8Wd7y6nF/83zLqGhq9jiUxTuUoYSElMY5Hri3khjN68szsTUx4Yg6VtfVex5IYpnKUsOH3GT85rx//+/VBfLJuN9/553xq6nUmjXhD5Shh5/Khedx16UBmFu3ilucWUa9NbPGAylHC0hXDuvDLrxbw1vLt/OTlpTQ2apBGWpcu9CFh6/rR3amorufed9aQmhTHnRcW0HSNNpHQUzlKWLv5zF5U1NTx2Mz1lFXW8t0zetKvU5rXsSQGqBwlrJkZPzu/P/F+H098uJ7/W1TMyd0zmTCyG+cUZBPn154hCY2QXpq1pWgmcAEoq6zl+bmb+ecnG9lSWkVOehJ/u2Yog3Rethyj454JXCQcZLRJ4Dun9+T928bw2DcKMeCHzy+iuk6H+0jLUzlKxPH7jLMLsvnfrw9i3a79/Pnfq72OJFFI5SgRa2SvDlx9cj5/n7We+Rv3eB1HoozKUSLaT8/vT256Mre9tESb19KiVI4S0VIS47j7soGs27mfe6et8TqORBGVo0S80b07MH54Po/NXMeCTaVex5EooXKUqPCz8/vRKS2JW19YzKrte72OI1FA5ShRITUpnnuuHMye/bWcf/9Mfj5lKbsraryOJRFM5ShRY0SP9rx/2xl845RuTJ67mTP+dwZ/n7mO2nrN6iNHT+UoUSWjTQK/uugE3rrlVIbkt+N3b6zkv56aq0svyFFTOUpU6p2dyqTrhvGLC/oz69NdvLl0u9eRJMKErBzNLMnM5pjZYjNbbma/DizvbmazzazIzJ43s4RQZZDYZmZcN6o7/Tql8oc3V+o4SDkqoVxzrAHGOucGAYOB88xsBHA3cK9zrjdQClwfwgwS4/w+444LC9haVsXjs9Z7HUciSMjK0TWpCDyND9wcMBZ4KbB8EnBJqDKIAIzs2YFzT8jmwfc+Zcfeaq/jSIQI6T5HM/Ob2SKgBJgGrAXKnHMHLiu3Beh8iK+daGbzzGzezp07QxlTYsDPzu9PfYPjj29pkgppnpCWo3OuwTk3GMgDhgP9g73tEF/7qHOu0DlXmJWVFcqYEgO6tm/LdaO78fKCLSzeXOZ1HIkArTJa7ZwrA2YAI4AMMzswA3keUNwaGUS+N6YXHVIS+c3rK3RojxxRKEers8wsI/A4GTgLWAm8B1weeNsE4NVQZRA5WGpSPLed24f5G0v5x8cbvY4jYS6U15DJASaZmZ+mEn7BOfe6ma0AJpvZ74CFwOMhzCDyOZcP7cJby7Zz59Tl1DU08q1Te3gdScJUyMrRObcEGBJk+Tqa9j+KtDq/z3jk2kJ+8PxCfvfGSipq6rnlzN665Kt8ia4+KDEnIc7HA+OG0DZhKfe9U8S+6np+cUF/FaR8jspRYlKc38fdlw2kbWIcj89aT0V1Pf9z6Yn4fCpIaaJylJjl8xl3XlhAalIcf3n3U3pnp2gfpHxGE09ITDMz/vvsPpxTkM0f31rNimJNlCtNVI4S88yMuy4bSEabeH7w/EJNUCGAylEEgMy2Cfzp64NYs6OCu/61yus4EgZUjiIBp/fJ4psju/HURxuYsbrE6zjiMZWjyEFu/0o/+mSn8KMXl+gaNDFO5ShykKR4P/ePG8Leqjp++MJi7X+MYSpHkS/on5PGry8+gQ/W7OQbj8+hrLLW60jiAZWjSBDjh+fzl/FDWLS5jMv/9jFbSiu9jiStTOUocggXDspl0n8NZ8feai596COWF5d7HUlakcpR5DBO6dmel28Yid9nXPnIJ8xet9vrSNJKVI4iR9AnO5UpN44iOy2RG55ZwLbyKq8jSStQOYo0Q6f0JB79RiE1dQ3c9MwC6hoavY4kIaZyFGmmnlkp3H35QBZsKuN/3tRZNNFO5ShyFL46MJdvjuzGEx+u582l27yOIyGkchQ5Sj87vz9D8jP48UtLWLez4shfIBFJ5ShylBLifDx41UnE+40bn1nAvuo6ryNJCKgcRY5BbkYy948bQlFJBV976CPW79rvdSRpYSpHkWN0Wp8snr7+ZHZX1HDxX2fx/pqdXkeSFqRyFDkOp/Rsz9TvjSY3I5nrnpzDox+sxTnndSxpASpHkePUJbMNr9w4kvMGdOIPb67i1hcXU6/jICOeylGkBbRJiOPBq07ih2f14ZUFW7n1xcU0NGoNMpLp6oMiLcTMuOWs3sT5jT+9vZp4v48/XjZQl3uNUCpHkRZ205he1NY3cv/0IhLifPz+kgGYqSAjjcpRJAR+cFZvahsaeXjGWhL8Pu68sEAFGWFUjiIhYGb8+Ny+1NY38vis9dQ3NnLnhScQ79du/kihchQJETPjFxf0J85nPPLBOtZsr+CvVw+hY2qS19GkGUL2Z8zMupjZe2a20syWm9ktgeWZZjbNzIoC9+1ClUHEa2bGT8/vz/3jBrNkaxkX/mUW8zeWeh1LmiGU6/j1wK3Ouf7ACOAmMysAbgemO+d6A9MDz0Wi2sWDOzPlxlEkxvkZ9+jHPDN7ow4WD3MhK0fn3Dbn3ILA433ASqAzcDEwKfC2ScAlocogEk7656Tx2vdGM6pXB34+ZRnXPTVXs/qEMWuNv15m1g34ABgAbHLOZRz0Wqlz7kub1mY2EZgIkJ+fP3Tjxo0hzynSGhoaHU9+uJ773imipr6B60f34Ptje9E2UUMArc3M5jvnCoO+FupyNLMU4H3g9865V8ysrDnleLDCwkI3b968kOYUaW0l+6q5+1+reXnBFrLTEvnZ+f25aFCuDvlpRYcrx5AeV2Bm8cDLwDPOuVcCi3eYWU7g9RygJJQZRMJVx9Qk/nzFIF6+YSQdU5O4ZfIivvfcQsqrND9kOAjlaLUBjwMrnXP3HPTSVGBC4PEE4NVQZRCJBEO7tuP/bhrFj8/ry9vLtnP+/TOZu2GP17FiXijXHEcB1wJjzWxR4HY+cBdwtpkVAWcHnovENL/PuPGMXrx0w0ji/MaVj3zMPdPWaHYfD7XKgMzx0j5HiSUVNfXc8eoyXlmwlRE9Mnnk2kLSk+O9jhWVPNvnKCJHLyUxjnuuGMyfvz6I+RtLueJvH7O9vNrrWDFH5SgSpi4bmsdT1w1na1kVlz70IUU79nkdKaaoHEXC2KheHZg8cQS1DY7L//axBmpakcpRJMwN6JzOlBtHktk2gWv+Ppvfvr6CuRv2aKbxENOAjEiE2LO/lp+9spR3V5VQ29BIh5REzi7I5vQ+HYjz+ahvbKS2wVHf0EjX9m0Y2jXT68hhz9MzZFqCylHkP/ZV1/He6p28vXw7M1aVsL+2Iej7RvVqz63n9OWkfE18dSgqR5EoVV3XwOrt+/CZEec34v2G3+dj+sodPDxjLbv313Jmv47cek5fCnLTvI4bdlSOIjFof009T364nkc+WMe+6noGdclgRPdMTu6RSWG3TNKSdOykylEkhpVX1vGPjzfw/pqdLN5SRl2Dw2fQr1MafTul0jOrLT2zUujZMYX8zDYkxfu9jtxqVI4iAkBVbQMLN5cye90eFmwqZW1JBcVfOMC8Q0oCOenJ5KQnkZuRzPDumZzVP5uEuOg7uEXlKCKHtL+mnvW79rN2ZwUbd1eyrbyK4rJqtpVXsbW0iv21DbRvm8DlQ/O4clgXemSlHPb7NTY6SvbVkBTvI6NNQiv9FsdG5Sgix6Sh0fFB0U4mz9nEOytLaGh0FHZtR5fMNiTG+UiK95MY56OuwbFpTyWb9uxn4+5KauqbJszonJHMgM5pnJCbzgm5aWS2TSA5wU9yfOCW4CclMc6zOSxVjiJy3Er2VfPS/C38a+l2yqvqqKlvoLqukZr6BgyjS2YyXdu3pWtmG7q2b0NlbQPLiveyfGs563btP+T39fuMjOR4MtrEk9Emga7t23Ba7yxG9+5Ah5TEkP5OKkcR8dS+6jrW7NjH3up6qmsbqKprulXWNFBWVUtpZR3llXWUVtayctteSiubJvw9ITeN0/pk0Sc7hey0JDqlJdEpPYk2CS1zSYnDlaMuWiEiIZeaFN/sM3YaGx3LisuZWbSL99fs5LEP1lH/hVMlUxLjSIzzBY7t9BHv95Gdlsjkiae0WGaVo4iEFZ/PGJiXwcC8DG4a04uq2gaKy6vYUV7NtvJqtu+tZldFDbX1jdQ3OOoaGqlrdKS08AXKVI4iEtaSE/xNx2EeYZS8pUXfgUsiIi1A5SgiEoTKUUQkCJWjiEgQKkcRkSBUjiIiQagcRUSCUDmKiAQREedWm9lOYONRflkHYFcI4oRapOYGZfeKsh+7rs65rGAvREQ5Hgszm3eoE8rDWaTmBmX3irKHhjarRUSCUDmKiAQRzeX4qNcBjlGk5gZl94qyh0DU7nMUETke0bzmKCJyzFSOIiJBRF05mtl5ZrbazD41s9u9znM4ZvaEmZWY2bKDlmWa2TQzKwrct/My46GYWRcze8/MVprZcjO7JbA87PObWZKZzTGzxYHsvw4s725mswPZnzezsL2uqJn5zWyhmb0eeB4R2c1sg5ktNbNFZjYvsCwsPzNRVY5m5gceBL4CFADjzazA21SH9RRw3heW3Q5Md871BqYHnoejeuBW51x/YARwU+DfOhLy1wBjnXODgMHAeWY2ArgbuDeQvRS43sOMR3ILsPKg55GUfYxzbvBBxzeG5WcmqsoRGA586pxb55yrBSYDF3uc6ZCccx8Ae76w+GJgUuDxJOCSVg3VTM65bc65BYHH+2j6j9qZCMjvmlQEnsYHbg4YC7wUWB6W2QHMLA+4APh74LkRIdkPISw/M9FWjp2BzQc93xJYFkmynXPboKmAgI4e5zkiM+sGDAFmEyH5A5uli4ASYBqwFihzztUH3hLOn537gB8DjYHn7Ymc7A74t5nNN7OJgWVh+ZmJtgtsWZBlOlYphMwsBXgZ+IFzbm/TSkz4c841AIPNLAOYAvQP9rbWTXVkZvZVoMQ5N9/MzjiwOMhbwy57wCjnXLGZdQSmmdkqrwMdSrStOW4Buhz0PA8o9ijLsdphZjkAgfsSj/MckpnF01SMzzjnXgksjpj8AM65MmAGTftNM8zswApDuH52RgEXmdkGmnYbjaVpTTISsuOcKw7cl9D0R2k4YdPJMigAAAJ1SURBVPqZibZynAv0DozcJQDjgKkeZzpaU4EJgccTgFc9zHJIgf1cjwMrnXP3HPRS2Oc3s6zAGiNmlgycRdM+0/eAywNvC8vszrmfOufynHPdaPp8v+ucu5oIyG5mbc0s9cBj4BxgGeH6mXHORdUNOB9YQ9M+pJ97necIWZ8DtgF1NK31Xk/T/qPpQFHgPtPrnIfIPpqmTbclwKLA7fxIyA8MBBYGsi8D7ggs7wHMAT4FXgQSvc56hN/jDOD1SMkeyLg4cFt+4P9nuH5mdPqgiEgQ0bZZLSLSIlSOIiJBqBxFRIJQOYqIBKFyFBEJQuUoYcnMGgIztxy4tdhkBGbW7eCZkESCibbTByV6VDnnBnsdQmKX1hwlogTmA7w7MB/jHDPrFVje1cymm9mSwH1+YHm2mU0JzN242MxGBr6V38weC8zn+O/AmTIin1E5SrhK/sJm9ZUHvbbXOTcc+CtN5xUTePwP59xA4BnggcDyB4D3XdPcjSfRdGYGQG/gQefcCUAZcFmIfx+JMDpDRsKSmVU451KCLN9A00S16wITX2x3zrU3s11AjnOuLrB8m3Oug5ntBPKcczUHfY9uwDTXNLkqZvYTIN4597vQ/2YSKbTmKJHIHeLxod4TTM1BjxvQ/nf5ApWjRKIrD7r/OPD4I5pmqQG4GpgVeDwduAE+m+A2rbVCSmTTX0sJV8mBmboPeMs5d+BwnkQzm03TH/fxgWU3A0+Y2W3ATuC6wPJbgEfN7Hqa1hBvoGkmJJHD0j5HiSiBfY6FzrldXmeR6KbNahGRILTmKCIShNYcRUSCUDmKiAShchQRCULlKCIShMpRRCSI/wekwukX8akBIgAAAABJRU5ErkJggg==\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.scatter(\n",
" states['RTK_signalling'],\n",
" ast.get_state_dataset().get_exprs_df()['Her2']\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cell state diagnostics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is important to run diagnostics on cell states model for the same reasons\n",
"stated for the cell type model. `Astir.diagnostics_cellstate()` spots any non\n",
" marker protein and pathway pairs whose expressions are higher than those of\n",
" the marker proteins of the pathway."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
pathway
\n",
"
protein A
\n",
"
correlation of protein A
\n",
"
protein B
\n",
"
correlation of protein B
\n",
"
note
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
RTK_signalling
\n",
"
EGFR
\n",
"
0.823298
\n",
"
Cleaved Caspase3
\n",
"
0.831128
\n",
"
EGFR is marker for RTK_signalling but Cleaved ...
\n",
"
\n",
"
\n",
"
1
\n",
"
RTK_signalling
\n",
"
EGFR
\n",
"
0.823298
\n",
"
cleaved PARP
\n",
"
0.831128
\n",
"
EGFR is marker for RTK_signalling but cleaved ...
\n",
"
\n",
"
\n",
"
2
\n",
"
cell_growth
\n",
"
Ki-67
\n",
"
0.058879
\n",
"
Cleaved Caspase3
\n",
"
0.862459
\n",
"
Ki-67 is marker for cell_growth but Cleaved Ca...
\n",
"
\n",
"
\n",
"
3
\n",
"
cell_growth
\n",
"
Ki-67
\n",
"
0.058879
\n",
"
EGFR
\n",
"
0.897650
\n",
"
Ki-67 is marker for cell_growth but EGFR isn't
\n",
"
\n",
"
\n",
"
4
\n",
"
cell_growth
\n",
"
Ki-67
\n",
"
0.058879
\n",
"
Her2
\n",
"
0.750985
\n",
"
Ki-67 is marker for cell_growth but Her2 isn't
\n",
"
\n",
"
\n",
"
5
\n",
"
cell_growth
\n",
"
Ki-67
\n",
"
0.058879
\n",
"
cleaved PARP
\n",
"
0.862459
\n",
"
Ki-67 is marker for cell_growth but cleaved PA...
\n",
"
\n",
"
\n",
"
6
\n",
"
cell_growth
\n",
"
Ki-67
\n",
"
0.058879
\n",
"
phospho S6
\n",
"
0.272989
\n",
"
Ki-67 is marker for cell_growth but phospho S6...
\n",
"
\n",
"
\n",
"
7
\n",
"
mTOR_signalling
\n",
"
phospho S6
\n",
"
0.699387
\n",
"
EGFR
\n",
"
0.720013
\n",
"
phospho S6 is marker for mTOR_signalling but E...
\n",
"
\n",
"
\n",
"
8
\n",
"
mTOR_signalling
\n",
"
phospho S6
\n",
"
0.699387
\n",
"
Her2
\n",
"
0.757867
\n",
"
phospho S6 is marker for mTOR_signalling but H...
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" pathway protein A correlation of protein A protein B \\\n",
"0 RTK_signalling EGFR 0.823298 Cleaved Caspase3 \n",
"1 RTK_signalling EGFR 0.823298 cleaved PARP \n",
"2 cell_growth Ki-67 0.058879 Cleaved Caspase3 \n",
"3 cell_growth Ki-67 0.058879 EGFR \n",
"4 cell_growth Ki-67 0.058879 Her2 \n",
"5 cell_growth Ki-67 0.058879 cleaved PARP \n",
"6 cell_growth Ki-67 0.058879 phospho S6 \n",
"7 mTOR_signalling phospho S6 0.699387 EGFR \n",
"8 mTOR_signalling phospho S6 0.699387 Her2 \n",
"\n",
" correlation of protein B note \n",
"0 0.831128 EGFR is marker for RTK_signalling but Cleaved ... \n",
"1 0.831128 EGFR is marker for RTK_signalling but cleaved ... \n",
"2 0.862459 Ki-67 is marker for cell_growth but Cleaved Ca... \n",
"3 0.897650 Ki-67 is marker for cell_growth but EGFR isn't \n",
"4 0.750985 Ki-67 is marker for cell_growth but Her2 isn't \n",
"5 0.862459 Ki-67 is marker for cell_growth but cleaved PA... \n",
"6 0.272989 Ki-67 is marker for cell_growth but phospho S6... \n",
"7 0.720013 phospho S6 is marker for mTOR_signalling but E... \n",
"8 0.757867 phospho S6 is marker for mTOR_signalling but H... "
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ast.diagnostics_cellstate().head(n=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calling `ast.diagnostics_cellstate()` returns a `pd.DataFrame`, where each\n",
"column corresponds to a particular protein and two cell types, with a warning\n",
" if the protein is not expressed at higher levels in the cell state for which\n",
" it is a marker than the cell state for which it is not."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The diagnostics:\n",
"\n",
"1. Get correlations between all cell states and proteins\n",
"\n",
"2. For each cell state *c*, get the smallest correlation with marker *g*\n",
"\n",
"3. For each cell state *c* and its non marker *g*, find any correlation that is\n",
"bigger than those smallest correlation for *c*.\n",
"\n",
"4. Any *c* and *g* pairs found in step 3 will be included in the output of\n",
"`Astir.diagnostics_cellstate()`, including an explanation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If multiple issues are found, the markers and cell states may need refined.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 4. Saving results "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Both cell type and cell state information can easily be saved to disk via"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"ast.type_to_csv(\"data/cell-types.csv\")\n",
"ast.state_to_csv(\"data/cell-states.csv\")"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
",stromal,B cells,T cells,macrophage,epithelial(basal),epithelial(luminal),Other\r\n",
"BaselTMA_SP43_115_X4Y8_1,5.9812641356773485e-05,0.43514272442399615,4.4133780759206555e-06,4.778704083811438e-05,0.00590960852651017,1.6351770806453748e-06,0.5588340188121422\r\n",
"BaselTMA_SP43_115_X4Y8_2,7.080264289768121e-05,0.38496458175145926,5.889272783675721e-06,4.3142984449199766e-05,0.04471036731616506,7.2442437475747445e-06,0.5701979717884975\r\n"
]
}
],
"source": [
"!head -n 3 data/cell-types.csv"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
",RTK_signalling,cell_growth,mTOR_signalling,apoptosis\r\n",
"BaselTMA_SP43_115_X4Y8_1,0.2191776340692214,0.17884886731101723,0.13343558969907734,0.1441052515306596\r\n",
"BaselTMA_SP43_115_X4Y8_2,0.29646752007327026,0.11967817733248733,0.20014708539822942,0.1579368118796618\r\n"
]
}
],
"source": [
"!head -n 3 data/cell-states.csv"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the first (unnamed) column always corresponds to the cell name/ID."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 5. Accessing internal functions and data "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Data stored in `astir` objects is in the form of an `SCDataSet`. These can be retrieved via"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"celltype_data = ast.get_type_dataset()\n",
"celltype_data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and similarly for cell state via `ast.get_state_dataset()`.\n",
"\n",
"These have several helper functions to retrieve relevant information to the dataset:"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"['BaselTMA_SP43_115_X4Y8_1',\n",
" 'BaselTMA_SP43_115_X4Y8_2',\n",
" 'BaselTMA_SP43_115_X4Y8_3',\n",
" 'BaselTMA_SP43_115_X4Y8_4']"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"celltype_data.get_cell_names()[0:4] # cell names"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"['stromal',\n",
" 'B cells',\n",
" 'T cells',\n",
" 'macrophage',\n",
" 'epithelial(basal)',\n",
" 'epithelial(luminal)']"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"celltype_data.get_classes() # cell type names"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6\n",
"14\n"
]
}
],
"source": [
"print(celltype_data.get_n_classes()) # number of cell types\n",
"print(celltype_data.get_n_features()) # number of features / proteins"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"tensor([[4.2121e-02, 5.2044e-02, 1.7348e-03, ..., 6.1314e-01, 4.4827e-02,\n",
" 3.8841e-01],\n",
" [0.0000e+00, 3.4770e-02, 0.0000e+00, ..., 9.4025e-01, 1.9424e-02,\n",
" 5.6716e-01],\n",
" [1.8200e-02, 7.8596e-02, 0.0000e+00, ..., 6.2852e-01, 4.0905e-02,\n",
" 8.4946e-01],\n",
" ...,\n",
" [1.4403e-01, 7.0877e-02, 2.0325e-01, ..., 1.7303e+00, 2.1434e-01,\n",
" 9.4889e-01],\n",
" [3.6400e-02, 9.9307e-02, 1.1815e-01, ..., 1.6467e+00, 1.3463e-01,\n",
" 1.9238e+00],\n",
" [5.4069e-02, 1.1861e-01, 5.3894e-02, ..., 1.4265e+00, 4.9443e-01,\n",
" 1.6126e+00]], dtype=torch.float64)"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"celltype_data.get_exprs() # Return a torch tensor corresponding to the expression data used"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" CD20 CD3 CD45 CD68 \\\n",
"BaselTMA_SP43_115_X4Y8_1 0.008424 0.010409 0.000347 0.008586 \n",
"BaselTMA_SP43_115_X4Y8_2 0.000000 0.006954 0.000000 0.016915 \n",
"BaselTMA_SP43_115_X4Y8_3 0.003640 0.015719 0.000000 0.008586 \n",
"BaselTMA_SP43_115_X4Y8_4 0.000000 0.003740 0.003136 0.021657 \n",
"BaselTMA_SP43_115_X4Y8_5 0.019956 0.014758 0.018645 0.024490 \n",
"... ... ... ... ... \n",
"BaselTMA_SP43_115_X4Y8_4927 0.000000 0.007295 0.028565 0.076878 \n",
"BaselTMA_SP43_115_X4Y8_4928 0.000000 0.020123 0.017928 0.018262 \n",
"BaselTMA_SP43_115_X4Y8_4929 0.028802 0.014175 0.040638 0.100518 \n",
"BaselTMA_SP43_115_X4Y8_4930 0.007280 0.019860 0.023629 0.089328 \n",
"BaselTMA_SP43_115_X4Y8_4931 0.010813 0.023720 0.010779 0.096799 \n",
"\n",
" Cytokeratin 14 Cytokeratin 19 Cytokeratin 5 \\\n",
"BaselTMA_SP43_115_X4Y8_1 0.009215 0.015349 0.022756 \n",
"BaselTMA_SP43_115_X4Y8_2 0.011813 0.011297 0.019773 \n",
"BaselTMA_SP43_115_X4Y8_3 0.005251 0.025360 0.033187 \n",
"BaselTMA_SP43_115_X4Y8_4 0.014276 0.010551 0.004497 \n",
"BaselTMA_SP43_115_X4Y8_5 0.000000 0.047767 0.038716 \n",
"... ... ... ... \n",
"BaselTMA_SP43_115_X4Y8_4927 0.009306 0.000123 0.027927 \n",
"BaselTMA_SP43_115_X4Y8_4928 0.000000 0.019686 0.000000 \n",
"BaselTMA_SP43_115_X4Y8_4929 0.010249 0.040175 0.027628 \n",
"BaselTMA_SP43_115_X4Y8_4930 0.057028 0.054079 0.021476 \n",
"BaselTMA_SP43_115_X4Y8_4931 0.068105 0.026760 0.030130 \n",
"\n",
" Cytokeratin 7 Cytokeratin 8/18 E-Cadherin \\\n",
"BaselTMA_SP43_115_X4Y8_1 0.014714 0.022104 0.159329 \n",
"BaselTMA_SP43_115_X4Y8_2 0.003848 0.030560 0.193268 \n",
"BaselTMA_SP43_115_X4Y8_3 0.062041 0.050745 0.220020 \n",
"BaselTMA_SP43_115_X4Y8_4 0.018720 0.012455 0.238691 \n",
"BaselTMA_SP43_115_X4Y8_5 0.000000 0.025077 0.242644 \n",
"... ... ... ... \n",
"BaselTMA_SP43_115_X4Y8_4927 0.000000 0.000000 0.339496 \n",
"BaselTMA_SP43_115_X4Y8_4928 0.000000 0.031084 0.180707 \n",
"BaselTMA_SP43_115_X4Y8_4929 0.009111 0.011427 0.415244 \n",
"BaselTMA_SP43_115_X4Y8_4930 0.084552 0.046873 0.368578 \n",
"BaselTMA_SP43_115_X4Y8_4931 0.073124 0.032700 0.360796 \n",
"\n",
" Fibronectin Her2 Vimentin pan Cytokeratin \n",
"BaselTMA_SP43_115_X4Y8_1 0.064861 0.122322 0.008965 0.077604 \n",
"BaselTMA_SP43_115_X4Y8_2 0.076731 0.186959 0.003885 0.113190 \n",
"BaselTMA_SP43_115_X4Y8_3 0.051344 0.125376 0.008181 0.169086 \n",
"BaselTMA_SP43_115_X4Y8_4 0.132006 0.105818 0.018769 0.083547 \n",
"BaselTMA_SP43_115_X4Y8_5 0.124412 0.098638 0.000000 0.059350 \n",
"... ... ... ... ... \n",
"BaselTMA_SP43_115_X4Y8_4927 0.046690 0.252541 0.034240 0.079310 \n",
"BaselTMA_SP43_115_X4Y8_4928 0.082254 0.073101 0.022108 0.157243 \n",
"BaselTMA_SP43_115_X4Y8_4929 0.167446 0.339493 0.042855 0.188656 \n",
"BaselTMA_SP43_115_X4Y8_4930 0.152866 0.323661 0.026924 0.375847 \n",
"BaselTMA_SP43_115_X4Y8_4931 0.230238 0.281564 0.098726 0.317179 \n",
"\n",
"[4931 rows x 14 columns]"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ast.get_type_dataset().get_exprs_df()"
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"source": [
"## 6. Saving models "
]
},
{
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"source": [
"After fixing the models, we can save the cell type/state assignment, the losses, the parameters (e.g. `mu`, `rho`, `log_sigma`, etc) and the run informations (e.g. `batch_size`, `learning_rate`, `delta_loss`, etc) to an hdf5 file."
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"ast.save_models(\"data/astir_summary.hdf5\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The hierarchy of the hdf5 file would be:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Only the model that is trained will be saved (`CellTypeModel` or `CellStateModel` or both). If the functioned is called before any model is trained, exception will be raised. Data saved in the file is either `int` or `np.array`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 7. Plot clustermap of expression data "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"After fixing the cell type model, we can also plot a heatmap of protein expression of cells clustered by type. The heatmap will be saved at the location `plot_name`, which is default to `\"./celltype_protein_cluster.png\"`"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAFzCAYAAAAjVEDpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3xUVfo/8M8zfSa995BCJp0QCEGaKMoCqyAsICiuigtYvhZ0bavufl2+uj/XtWJZZVdRbLgLCAgKFqpGgQSBkE5CGul9JpmZTDm/PzLBGEMSIAnJ8Lxfr3ll5t5z7z33zmSeOefeex4SQoAxxhhjjkNyqSvAGGOMsYHFwZ0xxhhzMBzcGWOMMQfDwZ0xxhhzMBzcGWOMMQfDwZ0xxhhzMLJLXQHGBkJGRoavTCb7N4AE8I9WxpjjswE4abFYVowfP76m+0wO7swhyGSyf/v7+8f6+Pg0SiQSHryBMebQbDYb1dbWxlVVVf0bwLzu87mFwxxFgo+PTwsHdsbY5UAikQgfH59mdPRW/nr+ENeHscEi4cDOGLuc2L/zeozjHNwZYw6ttLRUdv3110eEhIQkREZGxk+fPn30iRMnlL0to9FokgEgLy9PERUVFd913uHDh9UxMTFxMTExcW5ubmODgoISY2Ji4iZPnqwdzP1g7HzwOXfGmMOy2WyYN2/e6Jtvvrl+x44dRQCQlpamrqiokI8ZM8Z0IetMTU015ObmZgPAwoULw66//vrm5cuXNw5kvRm7WNxyZ4w5rB07drjIZDLx6KOP1nZOmzx5smH27Nl6APjzn//sl5CQEKvVauMefPDBwIvd3ty5c8M3btzo1vn6uuuui/j000/dXnrpJe+ZM2dGTp06NSosLCzh0UcfDegs89prr3klJibGxsTExN1yyy2hVqsVZrMZ8+fPD9dqtXFRUVHxzzzzjO/F1o1dXji4M8Yc1okTJ9RJSUltPc3bsmWL66lTp1QnTpzIycnJyT527Jjmyy+/dL6Y7a1YsaLuvffe8wKA2tpa6fHjx50WLlzYDADHjx93+uyzz4pOnDiRvXnzZs+0tDT1kSNHVNu2bXM/evRoTm5ubrbVaqV//etfngcPHnRqaGiQ5efnZxcUFGTddddd9RdTL3b54W55xthladeuXa4HDhxwjYuLiwOAtrY2SW5urmrOnDn6C13n3LlzdQ899FBoVVWV9N133/WcN29eg0zW8TU7ffr0Fh8fHysAzJkzp2nfvn3OFouFTpw44ZSYmBgHAEajURIcHNw+f/785qKiItXy5ctDrr/++uYFCxa0DMAus8sIB3fGmMNKTEw0bN261aOneUIIrF69uvKRRx6pG6jtSSQSLFq0qP6dd97x+vjjj70+/vjjos55RPSLuzmICEII3HTTTXWvvvpqRfd1ZWVlZW3evNnttdde8920aZPHJ598UjJQ9WSOj7vlGWMOa+7cubr29nZ68cUXvTun7d+/X7Nz507nOXPmtHzwwQfezc3NEgA4ffq0/MyZMxfd4Ln77rvrX3/9dT+lUimSkpLOXrR34MABt7q6OqlOp5Ps2rXLffr06fo5c+botm3b5llZWSkDgKqqKmlBQYGioqJCZrPZcMcddzSuWbOmIjMzU3Ox9WKXF265M8YclkQiwfbt2wvvueeekFdeecVfqVSK4OBg02uvvVaWmJhoysrKUk2YMCEGADQaje2jjz46HRQUZLmYbYaFhZnDw8NNixYtaug6fcKECbqFCxeGl5SUqG688cb6yZMnGwDg8ccfr7j66qu1NpsNcrlcvPnmmyVSqRQrV64ME0KAiPDss8+WX0yd2OWHhOBxP9jId/z48eKkpKQB615l7EK1tLRI4uLi4jMzM7M8PDxsAPDSSy95nzx5Uv3uu++WXer6Mcdy/Phx76SkpLDu07lbnjHGBsjmzZtdo6Oj4+++++7qzsDO2KXALXfmELjlzhi7HJ2r5d7rOff/XDF10CP/jT9+R4O9DcYYY+xy0mtwlyrlQ1UPxhhjjA2Q3oO7QjFU9WCMMcbYAOkjuHPLnTHGGBtper1aXqpUDPqDMUchlUrHx8TExEVHR8fFxcXFfv31104Xu861a9d63XrrraEA8NBDDwX+5S9/8bv4mg5fPaVY7Wu/ux4jxi6l4fQ/2mvLXcItdzZCrV3iOn4g13f/py0ZfZVRKpW2zlSgmzdvdn3iiSeCZ86cmTeQ9RhKhsb8AT2Gag9tn8eQjUwzP1g3oJ+Vr3+/6pJ/ViwWCzrzAoxEvbfcFYpBfzDmiJqbm6Vubm49jnT2+uuve2m12rjo6Oi4+fPnhwNARUWFbNasWZEJCQmxCQkJsV999VWvrf5nnnnGNzIyMl6r1cZdf/31EYOxD8NNampq9N133x2UmJgYGxYWlrBr165fZXDbuHGj29ixY2MqKytlCxcuDLv99ttDkpOTY4KDgxPXr1/vAXTkeL/zzjuDo6Ki4rVabdy//vUvDwC45ZZbQj/66CM3AJg5c2bk4sWLwwDg5Zdf9r7//vsD8/LyFBEREfFLly4dNXr06PgpU6ZE6fV6vtvnEsnLy1OEh4fHL1myZFRUVFT8vHnzwrdu3eoybty4mFGjRiXs3btXs3fvXk1ycnJMbGxsXHJycszx48eVQEfgXrVqVbBWq43TarVxzz77rC8ABAUFJT788MMB48ePj3733Xc90tLS1ElJSTFarTZu5syZkbW1tVKg47N4xx13hCQnJ8dERUXF79279+zwwDk5OerU1NTo4ODgxK6peq+99trI+Pj42NGjR8e/8MILZ4dDfvnll73DwsISUlNTo5cuXTqqsxfqfL8TuuNz7owNEJPJJImJiYkzmUxUV1cn/+KLL/K7l0lPT1e98MILAT/88ENuQECApbq6WgoAd955Z8hDDz1UPWvWLH1BQYFi1qxZUUVFRVnn2tbatWv9S0pKMtVqtairq5MO5n4NJxaLhTIzM3M+/fRTtzVr1gTOnj377DHesGGD+6uvvur39ddfF3RmX6uurpanp6fnHjt2TLVgwYLRy5cvb9ywYYN7ZmamOicnJ6uyslKWmpoa+5vf/EZ/5ZVX6g4cOOCybNmy5qqqKkVNTY0AgO+//975pptuagCA0tJS1Ycfflg0efLkkt/+9rcRGzZs8Ljnnnsaeq4tG2xlZWWqTz/9tGj8+PElY8aMif3oo4+80tPTcz/++GP3Z599NuA///nP6cOHD+fK5XJs3brV5dFHHw3evXt34YsvvuhTUlKizMrKypbL5ej8PwQAlUply8jIyAMArVYb9/LLL5ded911+tWrVwc+9thjgZ2jDLa1tUl++umn3C+//NJ51apV4QUFBVkAcOrUKVVaWlpeU1OTNDY2NuGRRx6pVSqV4qOPPir28/Oz6vV6Sk5OjrvlllsajUaj5IUXXgg4evRotru7u23y5Mna+Ph4A3D+3wnd8dXyjA2Qrt3y33zzjdPy5cvD8/PzsySSnzvIdu/e7Tp37tzGgIAACwD4+flZAeD77793LSgoUHeW0+v10sbGxnP2rEVHRxsWLFgQPm/evKZly5Y1DdpODTGinhvCndMXL17cCACTJ09ufeSRR85+QaWlpbkcP35cs3fv3nxPT8+zI8PNmzevSSqVYvz48cb6+no5ABw8eNDlxhtvbJDJZAgJCbFMnDhR/91332lmzpypf+ONN/wyMjJUWq3W0NTUJC0pKZFnZGQ4/etf/yqtqamRBQUFmTrHhE9OTm4rLi5WDt7RYH0JCgoypaamGgBAq9UaZsyY0SKRSDBu3Li2Z555JrChoUG6ZMmS8OLiYhURCbPZTACwZ88e17vuuqtWLu9owHb+HwLArbfe2ggA9fX1Up1OJ73uuuv0ALBy5cr6xYsXn+0lu/nmmxsAYM6cOXq9Xi/p/JH9m9/8pkmtVgu1Wm3x9PQ0l5eXyyIjI81///vf/Xbu3OkOAFVVVfKsrCxVRUWFfOLEibrO7S9YsKAxPz9fBZz7O6G/Ix/2cc6dgztjF+Laa69tbWxslFVWVsq6JiKxJwL51eBQQgikp6fnODs792vgqL179xZ8+eWXLlu3bnV//vnnAwsKCk52flGNZH5+fpbm5uZf9EQ0NDRIw8PDTQCgUqkEAMhkMlit1rO/BEJDQ02lpaXKkydPqq688sq2zumd5YGOY9z1b3fh4eHm5uZm2eeff+42bdo0XUNDg2zDhg0eTk5ONg8PD1tNTQ0UCsXZhaVSqTAYDDyE9yXU9f2QSCRn32+pVAqr1UqPPfZY0PTp03Vff/11YV5enmLGjBnRwLn/DwHAxcWlX8Gz+w/RztdKpbLrZwQWi4V27Njhsn//fpf09PRcFxcXW2pqarTBYJD0NkLs+X4ndNfHOXf5oD8Yc0Q//fSTymazwc/P7xfn3WfPnt2yfft2z6qqKimAs92BU6dObfn73/9+9vxcWlqaGudgtVpRWFiomDt3ru7NN98s1+l00u4BcaRyc3Oz+fr6mrdt2+YCdByfffv2uc2YMUPf23LBwcHtmzdvPrV8+fLw9PR0VW9lp0+frtu0aZOnxWJBRUWF7PDhw87Tpk1rBYDx48fr3377bd9rr71Wf9VVV+nfeOMN/4kTJ/a6bTZ8tbS0SIODg9sB4O233z57nvvaa69teeutt3zMZjMA/KJbvpOXl5fV1dXV2nltxzvvvOM1adKks5+FTz75xAMAdu/e7ezi4mL18vKydl9Hp6amJqmbm5vVxcXF9tNPP6mOHz/uBADTpk1rPXTokEttba3UbDZj27ZtHp3LnM93Qk+45c7YAOk85w50/Or+5z//Wdz9atuUlBTjH//4x8pp06bFSCQSkZCQ0LZ58+bidevWla1YsSJUq9XGWa1Wmjhxom7y5MmlPW3HYrHQzTffHK7T6aRCCLrzzjurvb29z/nFMtK8//77p++5557Qxx57LAQAHnvssYr4+HhTX8slJSWZNmzYULRkyZLI7du3nzpXud///vdNaWlpzrGxsfFEJP7617+Wh4aGWgBg6tSp+oMHD7omJCSYTCZTe3Nzs/TKK6/UDdzesaH02GOPVa1YsSJ87dq1/tOmTWvpnP7ggw/W5ufnK2NiYuJlMpm47bbbap944ona7suvX7/+9N133z3q/vvvl4SGhpo++eST4s55Hh4e1uTk5Bi9Xi9dt27d6d7qsXDhwuZ169b5aLXauMjISGNSUlIr0NFb9OCDD1ZOmDAh1tfX16zVag1ubm5WADif74Se9Jo45sia/x30seUn/OWvfLUpu2icOIYxNlRSU1OjX3jhhbKup4AuVHNzs8TNzc1mNpsxa9as0bfffnvdrbfe2u/raC4ocQy33BljjLHB88gjjwQeOHDA1WQy0fTp01tuueWWAblAtver5eUc3BljjLGuDh8+PGCDU61bt658oNbVFbfcGWOMMQfD97kzxhhjDqbX4E7ySz+uLhGFANgAwB+ADcA6IcSr3cpcBWAbgM4rFrcIIdYMZT0ZY4yx4WIktNwtAP4ohDhKRC4AMojoayFEdrdyB4UQ11+C+jHGGGPDSq+D2EgUykF/9EUIUSmEOGp/rgOQAyBoQPaesQFWWFgov+aaayJHjRqVEBISkrB8+fIQo9FIaWlp6k8//dSts9xwSg05Uq1evTpw69atLgCwZs0aX51Od/b7TKPRJJ/PurqmjX3++ed9Xn/9da/eynd//+64446QL7/80hnoSD5SWVk5aN2enftWUVEhmzZtWtRgbWck6v45GCrn+3kbCiOh5X4WEYUBSAZwqIfZk4joOIAKAA8LIfo9wD5zPGkvXzmgKSgnP3igzxSUNpsN8+fPH71ixYqaBx54oNBiseDmm28e9cADDwTFx8cb0tPTnZYsWdI8EPUZinSUp5rqB/QYjnb3GtA0nq+88kpF5/O3337bb+XKlQ39HTq0N48++uivBjPpTXV1tTQjI8OpM6HIUAkMDLT4+fmZv/rqK6ff/OY3rUO57e6U964Y0M+K6fV/X9BnpbfPwUhP4Xq++mi5Kwb9QUSriCi9y2NVT3UhImcAmwGsFkK0dJt9FMAoIUQSgNcAbB2Ig8PY+fj8889dlEql7YEHHqgHOsY/f+utt8o++eQT76effjr4888/94iJiTmbYvRcqSHffPNNz8TExNiYmJi4m2++eZTF0jGCrUajSV69enXgmDFjYr799ttfpTt1BD3tu0ajSV65cmVwXFxc7KRJk7QVFRUyAFi4cGHY+vXrPZ555hnfmpoa+fTp07UTJ07Udq7rvvvuC4qOjo5LSkqKKSsrkwH9S6PZtVX+4osveickJMRGR0fHzZo1K7KnVuEHH3zgcc011/ziO2nNmjV+iYmJsYmJibEnT55UAsDHH3/sNmbMmJjY2Ni4yZMnazvrtHPnTueYmJi4mJiYuNjY2LjGxkZJc3OzZNKkSdq4uLhYrVYb9+GHH7r3dLzmz5/ftGHDhl57GRxVS0uL5KqrrhodHR0dFxUVFf/HP/4xoPvnoPv/zLZt21xiY2PjtFpt3OLFi8MMBgMBHb0t9957b9DYsWNjEhISYr/77jvN1KlTo0JCQhKef/55H6BjsJn+vCfDRe9jyysVg/4QQqwTQqR0eazrXg8ikqMjsH8khNjSfb4QokUIobc//wKAnIi8u5djbDBlZmaqk5KSfjFilaenpy0oKKj9oYceqpw7d25jbm5u9sqVKxuBjtSQ+/fvzz9y5EjOCy+8EGgymejo0aOqTZs2eaanp+fm5uZmSyQS8dZbb3kBgMFgkCQkJBhOnDiRO2vWLIcb7/xc+24wGCTjxo1ry87OzpkyZYru8ccfD+y63FNPPVXj6+tr3r9/f/6hQ4fygY5jNWnSJH1eXl72pEmT9K+99poP8HMazZMnT+Z89tlnhXfddVdYb3VatmxZ48mTJ3Py8vKyo6OjDWvXrv3V90paWppzSkrKL1rOrq6u1szMzJw777yz5r777gsBgJkzZ+qPHTuWm5OTk71o0aKGNWvW+APAiy++6L927dqS3Nzc7B9//DHX2dnZptFobDt37jyVnZ2ds3///vwnnngi2Gb7dafElClTWg8fPuyQP/T6smXLFld/f39zXl5edkFBQdYTTzzR4+eg839m2rRprXfeeWf4p59+Wpifn59tsVjwj3/8w6dzfSEhIe3Hjh3LnThxov6OO+4I+/zzzwsPHTqU+9xzzwUCQH/fk+Gij275S5/NkDpS7bwDIEcI8dI5yvgDqBZCCCJKRcePlvohrCZjvWZ86ymVaU+pIXft2uVy8uRJTVJSUiwAGI1Gia+vrwXoyDB1++23Nw76jlwi59p3iUSCFStWNADAHXfcUf+73/1udF/rksvlYunSpc0AMH78+NZvvvnGFTj/1LoZGRnqv/zlL0E6nU7a2toqnT59+q9Oq1RXV8u7Jwi67bbbGgBg5cqVDU899VQIAJw+fVoxf/784NraWnl7e7skJCTEBABXXHGF/uGHHw658cYbG2666abGyMhIm8lkotWrVwf/+OOPzhKJBDU1NYry8nJZ5xj4nQIDAy01NTXD6/zpEBk3bpzhySefDLn77ruDbrjhhubZs2f/6gdv1/+Z48ePq4KDg01jxowxAcDtt99e/8Ybb/gCqAGAG2+8sQkAEhMT21pbWyUeHh42Dw8Pm1KptNXV1UldXFxs/XlPhouRMIjNFAC/B5BJRMfs054AEAoAQoi3ACwCcDcRWQAYACwVvQ2az9ggSExMNHTN6gQADQ0NkqqqKoVUKv3V57Gn1JBCCFq8eHH9G2+8caZ7eYVCYXPkc4bn2vdXX301oOvrc+V870omkwmJRNL5HBaLhezbOK80mqtWrQrftGnTqUmTJhnWrl3rtX//fpfuZVQqla176tfObdvrKwDg3nvvDX3ggQeqli1b1rxjxw6XNWvWBALA3/72t6r58+c3b9u2zW3y5Mmxu3btyj948KBTfX29LDMzM0epVIqgoKDEntLLtrW1kVKpHL7Nx0E0ZswY09GjR7M3b97s9uSTTwZ988033U/X/uJ/pq+Q0JkuViKR/CqVrNlsprffftuzP+/JcNFHt7xy0B99EUJ8J4QgIcQYIcRY++MLIcRb9sAOIcTrQoh4IUSSEOIKIUTaAB0fxvpt3rx5OqPRKOm80tpiseCee+4JWbx4cZ2/v79Zr9f3+UUwe/bslh07dnicOXNGBnRcrJWfnz8sfmUPtnPtu81mw/r16z0A4L333vNKTU39VZY2Jycna3Nzc5/H93zTaLa1tUlCQ0PNJpOJNm7c6NlTmejoaGN+fv4vvsw2bNjgCQDvvPOOR3JycisA6HQ6aWhoqLlzPzrLZmVlKVNTUw3PPvtsVWJiYuvJkydVzc3NUm9vb7NSqRSff/65S0VFRY+fgZMnT6q0Wq2hr/12RMXFxXIXFxfbPffc07B69erqY8eOaXr7HIwdO9Z45swZRec1EBs2bPCaNm1avzP+9fc9GS5G1NXyjA1nEokEW7duPbVq1apR//jHPwJsNhtmzJjRvHbt2jMtLS2SF154ISAmJibuj3/8Y+W51jF+/HjjU089deaaa67R2mw2yOVysXbt2lKtVts+lPtyKZxr39VqtS0rK0sdHx/v7+LiYt2yZUtR92Vvu+22ujlz5kT5+vqaO8+39uR802g+/vjjFampqbFBQUHtsbGxbXq9/ld5v+fNm9f8z3/+0+ehhx46m5XQZDLRmDFjYmw2G23cuLEIAJ588smKm266KdLPz689JSWltbS0VAkAzz//vG9aWpqrRCIRWq3WsGjRouampibpnDlzRickJMTGx8e3hYeHG3uq39dff+0ye/bsAbkDY6TJyMhQ/+lPfwqWSCSQyWTizTffLDl48KDzuT4HGo1GvPXWW8WLFy+OtFqtSEpKanv44Yf7fWfEihUrGvrzngwXvaZ8bT5VMOhd226jozjlK7tonPLVcWk0muS2trafLnU9ejN+/Pjo3bt3n/L29rYO5XZTUlKiv/zyy1M+Pj5Dul02fFxQylepklvujDHWl3/84x/lhYWFCm9v7yHrIq+oqJA98MAD1RzYWU+G/dXyjLHL23BvtQPAjBkzhnwQmcDAQMvvf//7Acn9zRwPn3NnjDHGHEwfWeHkQ1UPxhhjjA2Q3m+alQzbW/gYY4wxdg69BncBvpCdMcYYG2l6bZpbrGLQH4w5gqqqKmln8g9vb+8kX1/fMZ2vjUbjBf1K3rFjh8vVV189GvhlSlLWYbimfE1NTY0+cOCA5ny2fy7Tp08fXVdX96t76/ujax1XrVoVvH379l+NrsccV68t93bL4N9hMSD/AYx1U3XsjQFNQek/9n96TUHp7+9vzc3NzQY6vlSdnZ2ta9asqR7IOgy1t3OPDugxvDNmHKd8PU/79+8/NRDrefjhh2uWL18+at68ef0ekY2NbL223NvNlkF/MHY52bRpk2tcXFxsdHR03KRJk7RAR+rKxYsXhyUkJMTGxsb2mUry3Xff9YiKioqPjo6OS0lJiR6amg8NR0n52qlrD8L69es9Fi5cGNZZ92XLloVOnDhRGxwcnLhz507nxYsXh0VERMR3lgE6UpFWVlbK8vLyFBEREfFLly4dNXr06PgpU6ZE6fV66m8dtVpte1NTk6y0tNRxkxOwX+g1uJsttkF/MHa5qKiokN17771hW7ZsKczLy8veunVrIQA88cQTAVdffXXLyZMncw4ePJj31FNPBbe0tJzzf/O5554L+Oqrr/Lz8vKyd+3aNSAtu+HAkVK+9kdzc7Pshx9+yH/uuefKlixZEvXII49UFxQUZOXm5qp7GvO+tLRUdf/999ecOnUqy83NzbphwwaP/tYR6Mh2tmfPnssyPezlqNdfcSZuWTM2YPbt2+eUmpqqi4mJaQcAPz8/q3266+7du93Xrl3rD3SMS37q1KlzDjKRkpKiX7ZsWdjChQsbly1b5jApYB0p5Wt/XHfddU0SiQTjxo1r8/LyMqemphoAQKvVGgoLC5WTJ0/+xWh3QUFBps5pycnJbcXFxcr+1hEAfHx8LGfOnOHBSy4TvZ9zb+dRDRkbKOfK6y6EwKZNm04lJSWZuk6vqKjocaCJjz/+uHTPnj1O27dvdxs7dmz8sWPHsvz9/Uf8P6sjpXztqa4Gg+EXFe9MMSqVSn+VYrSzvl11LSOVSkXnNvtTRwAwGo2kVqu5u/Qy0Wu3vMlsGfQHY5eLq6++uvXQoUMuubm5CqDjQiz79JYXX3zRz2br+N79/vvve01DmpWVpZwxY0brK6+8UuHh4WEpKipyiNaYI6V87eTl5WU+evSoymq1Ytu2bR591e9C9KeOAFBYWKhKSkq6LNPDXo56b7mbR3xjgLFhIzAw0LJ27driBQsWjLbZbPDy8jKnpaUVPPfccxWrVq0KjYmJiRNCUHBwsGnv3r3nPJf+4IMPBhcXFyuFEDR16tSWK664wiG+sB0p5Wunv/71r2duuOGG0QEBAeaYmBhDa2vrgI8M1p86mkwmKi4uVl555ZVDPgY+uzR6Tfn6Xfrgp3ydmtJ3ylcimg3gVQBSAP8WQjzXbb4SwAYA4wHUA1gihCge+Nqy4YpTvjouTvl68TZs2OCekZGhefXVVyv6Ls1GkgtK+TocWu5EJAXwBoCZAMoBHCGi7UKI7C7F/gCgUQgxmoiWAvg7gCVDX1vG2OXoUqR8PR8Wi4X+/Oc/j+hxF9j56WMQm2Fx7UUqgFNCiCIAIKKNAG4A0DW43wDgafvzTQBeJyISvXVLMMZGhOHeagcuTcrX83HHHXc4zF0VrH/6uBXu0rfcAQQB6DryUzmAiecqI4SwEFEzAC8A3E3LGGPsstNHt/zgt9yJaBWAVV0mrRNCrOtapIfFurfI+1OGOTabzWYjiUTC7ztj7LJgs9kIQI+BuveW+xCMLW8P5Ot6KVIOIKTL62AA3S8K6SxTTkQyAG4AGgaynmzYO1lbWxvn4+PTzAGeMebobDYb1dbWugE42dP8S95y74cjAKKIKBzAGQBLAdzcrcx2ALcB+AHAIgB7+Hz75cVisayoqqr6d1VVVQL6GL+BMcYcgA3ASYvFsqKnmcP+gjr7OfR7AexGx61w7wohsohoDYB0IfwsknEAACAASURBVMR2AO8A+ICITqGjxb700tWYXQrjx4+vATDvUteDMcaGgz4uqLv0wR0AhBBfAPii27S/dHluBLB4qOvFGGOMDUfDvuXOGGOMsfPTR8udT1szxhhjIw233BljjDEH08etcNxyZ4wxxkaaPoI7t9wZY4yxkaaP+9y55c4YY4yNNL0G9+ExtDxjjDHGzkevwd3I59wZY4yxEYcvqGOMMcYcTB/BfaiqwRhjjLGB0nu3vHmoqsEYY4yxgdJrcDdwy50xxhgbcTi4M8YYYw6m1+DeysGdMcYYG3FICL4injHGGHMkkktdAcYYY4wNLA7ujDHGmIPh4M4YY4w5GA7ujDHGmIPh4M4YY4w5GA7ujDHGmIPh4M4YY4w5mF4HsWEXp+rYG+c1iIDbqJmDVZVh4Qx5XeoqsH76tqrkUldh0G06kt7j9AOHDg9xTUYu0+v/pktdB9YzbrkzxhhjDoaDO2OMsUuOiG4notftz58moof7KBs4dLUbeTi4M8YYG2luB8DBvRcc3BljjA0aIrqViE4Q0XEi+oCIfIhoMxEdsT+mnOf6FgFIAfARER0jouuI6LMu82cS0Rb7cz0RvUhER4noWyLysU+PJKJdRJRBRAeJKGYg93k44ODOGGPsghDRKiJK7/JY1W1+PIAnAcwQQiQBeADAqwBeFkJMALAQwL/PZ5tCiE0A0gEsE0KMBfAFgNjOwA1gOYD19udOAI4KIcYB2A/gf+3T1wG4TwgxHsDDAN48rx0fAfhqecYYYxdECLEOHYHyXGYA2CSEqLOXbyCiawHEEZ290N6ViFwuog6CiD4AcAsRrQcwCcCt9tk2AJ/an38IYAsROQOYDOC/XeqgvNDtD1cc3BljjA0WAtD9lmAJgElCCMMvCtJF3VW3HsDnAIwA/iuEsJyjnLBvv8ne6ndY3C3PGGNssHwL4EaijkEuiMgTwFcA7u0sQEQXEmR1AM629oUQFQAqADwF4L0u5SQAFtmf3wzgOyFEC4DTRLTYvn0ioqQLqMOwxi13xhhjg0IIkUVEzwLYT0RWAD8BuB/AG0R0Ah0x6ACAu85z1e8BeIuIDPi5F+AjAD5CiOwu5VoBxBNRBoBmAEvs05cB+CcRPQVADmAjgOMXso/DFQlxXoOoMcYYY8OO/R75n4QQ73SZphdCOF/Cal0yHNwZY4yNaPaWeSuAmUIIU5fpHNwZY4yx4YSI3gDQ/T74V4UQ63sqz37GwZ0xxhhzMHxB3SBqb2sTusICuERGob2+Fpa2NjiFjEJxVROsVoGwQHc0thhhsdqw+UAhVl0VgncPluOO6SFY9s8czB/vDjcnOdRKGdrNNqiVMuzLrMWt10Zi/VenEOihxNHiVkyPdcOyqyJR02rBE++fQItRQCkD/nZrPDYfKMKeHAPiA2X43eQgqBQyvLGzCL+/Kgh6gxlymQT55ToAQICnCvMmhODvW7Kx7NrRKChrQHG1HlszdLgmToMdJ9owyoPQ0CYwJ9EZWeUGBHnIIZcRgrzUyC3Xo6jWjMWTfJBT1oL9eUYsnOCGk2V6FNRY8cjcYLy+qxxLrvDG9ox6PHB9ODYeKEVMkBMCPdUorNRjT44eEyNU+P6UEXdf6489mXXwcZVjYrQX2owWnKlvg4+bEk4qOY7kN8DTRY5wf2dIJRI8uekMlk92QYCnGsXVrThS1IbfJntAZzAju7wNo7yVGBPugcr6Nvh7qpFZ3ARtkAuqGo2ICnLFmbpWtLSZUVpnRHK4G1w0cqTl1CMqwAnuzgpkl7XATSODxSrQbrEhyEsNpVyCBl07XNRyvHewDi5KIHmUGgfyDZgZ5wRntQyHCnSYGOWCTw81I9CNUFgvcEOSGifLjbhhgjea9O3YdrQZM2Kd0G6x4XipEfWtAitn+OBfe2oR4S1FqLcCH6YbcFOyCgcLjPB2JuhNAoFuUtS3WmETwD2zgqFrM2NnRi3CfBSobDKjodWGP1wTBJlUApVChtd2nsakKGeE+TlD12aGNsQdG/eX4Oarw/D/Np3Coiu84apRICHCB3uOlp59T/3dpPjd5CBIJRJ8n1UNLxcFfshvQXK4E6KD3XDfR+VYNl6FsREeKDijw7o0PRaNUcJksUEb4IS0/BZE+CpQVm/GK6vG47ZXjmBcmBrHSw2I8JEjr6od7VYgwluGPacseGy2F/69rx4po+Q4ccaMW6d6wdNFibLaVrS0WSCTEpIiPFHd2Aaj2YqqRhPKG0wI81ZBKiGE+zuhpc2M7RkNSB6lhkYhhatGBgkRfshvxq0zQvHPXSUIdJdBISNMivHC0VONCPRUYdk1WrQ3NUEn02DLgVPwdVchKtgdcWE+EABaWk0gIrhoFIDJiHapAkopUNnQBte2eqi8vCGsVuzPb8CUMSE4mleJ2iYDrp8yGgVl9dAo5Ug7WYmlV2vxwTe5uGG0EhKZFHUKdwTIzZCpVMitMWJUgDuU7a1oLiqC99hxqKzXI8BDjcY2MzQmHdp1OrTV1CDD5oNrolzxTX4zZkZ7oPFUAQDAEhEPAPDzdEZ5TTNqm9oQEegBJ7UCZ2paEBnkcfYetZZWE1ydlEjLLEN1owG/jXWHTKWCzMUVpvpamJqaYGpqgkQuh8rLC04BgahqMiIk2Iezwg1TfCscY4wx5mA4uDPGGGMOhoM7Y4wx5mA4uDPGGGMOhoM7Y4yxQUFEVnta1s5HGBGlENFa+/yniejhYVBPdyK6p8vrQCLadCnrdLH4annGGGODxdBDgpZidKRs7RcikgohrANaq19zB3AP7Klf7WPVL+p1iWGOW+6MMcaGDBFdRUQ7ukxKIqI9RFRARCu7lNlLRB8DyLRPe4iITtofq+3Twogoh4j+RURZRPQVEant8yKJaBcRZRDRQSKKsU/3I6LPiOi4/TEZwHMAIu29C/+wr/ekvfztRLTFvq4CInp+6I7WheOWO2OMscGiJqJj9uenhRALeigzBsAVAJwA/EREO+3TUwEkCCFOE9F4AMsBTERHGtlDRLQfQCOAKAA3CSFWEtF/ACxER+72dQDuEkIUENFEdLTKZwBYC2C/EGIBEUkBOAN43L6tsUDHj4ZudRwLIBmACUAeEb0mhCi7qCMzyEZUcE9JSXkcgOpS16M/0tPTn77UdWCMscFERKsArOoyaZ0QYl2X1z11y3e3zZ7VzUBEe9ER1JsAHBZCnLaXmQrgMyFEq327WwBMA7AdHT8aOn9AZAAIIyJnAJMB/LdLnnil/e8MALcCgL27v5mIPPqo47dCiGb7trMBjALAwX0AqThoMsbY8GAP5Ov6LNjHas7xurXLtN5GwjN1eW4FoEbHKeemfvyw6K/u2xj2sZPPuTPGGLuUbiAiFRF5AbgKwJEeyhwAMJ+INETkBGABgIPnWqEQogXAaSJaDADUIck++1sAd9unS4nIFYAOgMtA7dBwwMGdMcbYpXQYwE4APwL4P/uV6r8ghDgK4D172UMA/i2E+KmP9S4D8AciOg4gC8AN9ukPALiaiDLR0Y0fL4SoB/C9/WK9fwzAPl1yw75rgTHG2MjUUy51IcQ+APvsz58+x3Jny3SZ9hKAl7pNKwaQ0OX1C12enwYwu4d1V+PnQN91+s3dJiXYp7+Hjh8WneWu76nOww233BljjDEHw8GdMcYYczAc3BljjDEHI3366acvdR36bd26dVetWrVq36WuR3+lHSt+2ickEMcLqqB0cUFunQkhHiq4qaTwdHMCSSSoaWxFfXMb5k2JwI8FdThR3IKJSaG4JtEbyVo/OKvk+HBfKRZMCYOHqwr7T9bgRFED/u/2CRg72hdHcqtgMttQWquDXEpoM5hwy1UhKK/VIyncAwez6rD4Ch80t5pxRawfiAgnSxqREuWFrYcqEeihxJVJQUiM8MZHe0tx7YRQ7DxyBlV1emSW6FCnM+Oe2aGYP200nCUG3DApCKXVLXhg4Ri4KAWKa9rg76FERpEOFhtwz28jsOWHCsxNDUJVgw7xoS4orDTA11UKg9EMuRQoqjHAWSVBTlkzBAB/dyWOnGpCQ6sF0f4KWKxAfJAKWWU6SCWEhlYLKusNOF3TBqVcgsnxAThR1ICSOhNqWswYN9oDVqsNwmyEq1qKjKIWNLZa4OkkRUmtEVargNFig85oQ2ltG0rqjFDKCIXVRlwR443ECG98vLcYTa1mnGlsh1xKaG41o7LBCJsQsNoEIIBT1Qa0Gq245ZooZJ5uwOkaAwoq26AzWGCx2jA91g2NehPMFgEvJwlqWyxQyQlGsxWTY71Q29gKX1cZPNRATLATHvxdIvYdq0BpnRGuagl0BitaTTZolBJ4OhHqWkxwUxMSQjQoqTMhzIMwPcEbkb5ymNotUMkILmoJ3NRSqOWEwspW3HVdLH7MqUaLwQYPjQxOSkJmqQ6V9QacLGmGTArcPTce735diJrmdhwtbEK7RSAlyhPCYsKxYj1M7RZUNejh5qRAUbUBSjmQOtoNx0834WRJC0xmG+p1ZjirOtoGOWV6LEpxQ3VTOyL8nfBDfhNifKVwd5Kh1WRDkJcKdbp2OCklkEsJO46cgZOSYDILKGSAn5sC7WYbgjxkcFHJ4KkSONNoRKCbFCaLgI+LFGX1RjTq2lFUY0RTmwXNbVZklbaguskEvcGKmhYzwn1UqGxqh1ohQUFlG6w2W8fDKtDQakFNczuIAKkEaNKbYBMdnwujWaCstg0miwAgkDjaFwqJgNJJgwlaH8SM8kKr0YyqBj1cnVRQyGUQNWegcHNDWX0bvFxVsILgqlHAqnaBXKWEVC5HoK8rMgurMT46EC2tRoS4KeAKE1w9PaBWSODUrsP4hFC0KZzgrFFCZdTjQEkbosL94eGqATXVorRNAjeZgEyphMKkh7DZYJUqUNJogquvL+TWdkSHeUPh4YWoEE+UN7fDfdQoyHwC4O6shrNGAaFvhrypFr6hQXCTCZiqKqAjJRpaDFA118BqaEOzTQZUliIi1Bfx4T6QOTtDKJQwtFuxP7sW0WG+ONokg0HthlHhQTDbCG5OCkjl8r9e0i9Zdk5DekHdAAxCEzZAVRkSceE+sNpsZ19PSghGc6sJrvajbrH+fHtnZX0rymp+vq3zzc9z8ODCs9eJYO22HNw7LwYA8MzyCaiq18Pb3QkAMHNcAHJLG8+W/eRAKWSSjttC3dRS7M+qh7umY6MtrV1v1wRCfH6+3kUpJ5TXNJ99bbLYzq6nrLoZx063IC7Ms2NeuwXT4/3xfXYdAOCFFSl47N0MAMCq2ZEI8nUDDnaM8RAbpEZepeHsel1UUrS1dxyXm68aBavVhuxyPQCgQW8BEcHNXl9zl2NksghI6Je3uyplP7+22n4u66SUoNVk61JO8ot15Z75+Viv+yIPv5scjC+OdFyk224RUMs75hnabdAZ2+HprAAA6I02vPRZNoI8Ol6r5D93fn2X24iu1ZNKASF+3qa8S11vSPDE3z75CX5uHesxWwQWTQnE54crAeAX6ymp7fgx0tbesa6jRTpIznHX738PnEJbuw1Se4HO49yp6/Ezmm2YkeiFPZn1AID8KsPZ+c4qGcrr2s6WzSnXw1kl/cWyGoXEfrxsKKn9+f3VKCS/2G52mf7sc73JhkhfNQprOsoLAVQ3m8/ON1ttXZ6Ls/sBdBw/qYTOvs9d5wnxy9dd32sACPRQoKKxHUDH52jV3ET85f1fD29e2WTGsfxKXBHhDgBoKSqES3DI2flKGcFsH+W8oLwRUX5drhez2aAmKwQ6jlPXz0ZqpAeE+HnfRmlskHv6A/Yh02tNErja59W3GOCpJHQOpn5cJ8e0CA3Mel1HHRQ/f22XkBtG1dZC7tpRX09XzS/2R9JuOLsetdWIk9UGhNsXjwzyhKFcd7asxtfv7PPs4jpEaKyAmzcAoLz157qbLDZU1+sQ6uf2q+PHho+hvlr+ogahSUlJueBlGWOMscsFn3NnjDHGHAwHd8YYY8zBcHBnjDHGHAwHd8YYYwOOiPYR0axu01YT0btEtGmQtx1GRDd3eZ1CRGsHc5vDDQd3xhhjg+ETAEu7TVsKYL0QYtEgbzsMwNngLoRIF0LcP8jbHFY4uDPGGBsMmwBcT0RKoKM1DSAQQDkRnbRPu52IthLR50R0mojuJaKHiOgnIvqRiDzt5SKJaBcRZRDRQSKKsU9/j4jWElEaERURUeePhucATCOiY0T0IBFdRUQ77Ms8be892GdfxiGDPgd3xhhjA86eae0wfk7eshTAp/h1/vYEdLSyUwE8C6BNCJEM4AcAt9rLrANwnxBiPICHAbzZZfkAAFMBXI+OoA4AjwM4KIQYK4R4uYfqxQCYZd/m/xKR/EL3c7gaaVnhjCPlXveLuZ+fMcZGAiJaBWBVl0nrhBDrurzu7JrfZv97Rw+r2SuE0AHQEVEzgM/t0zMBjCEiZwCTAfyXfh6ESdll+a2iY4SgbCLyQ//sFEKYAJiIqAaAH4Dyfi47Ioyo4J6env5c36UYY4wNBXsgX9dLka0AXiKicQDUQoij9u75rroOm2nr8tqGjhglAdAkhBh7jm10Xf4cYzf2uowVIywW9gd3yzPGGBsUQgg9OvKyv4uOVvyFrKMFwGkiWgwA1CGpj8V0AFwuZHuOgoM7Y4yxwfQJgCQAGy9iHcsA/IGIjgPIAnBDH+VPALAQ0XEievAitjtiOVxXBGOMseFDCPEZunSXCyGK0XERHYQQ7wF4r8u8sC7Pz84TQpzGzxfmdV337d1eO9v/mgFc0634Pvu8p7stkwAHxC13xhhjzMFwcGeMMcYcDAd3xhhjzMFwcGeMMcYczHlfUJeSkvI4ANUFbi/sApcbkXSHv0NtRjosZyrgvHo1Dv75YcidndFaVY2r1r+P/XfcBu/EeNRlZqEOwOTnX4OrUx1Mp3JQ02KBU3srVB7uIAL+3+3jsDu9BKO8lThV3oCC8kYUV7fC102BfSeqUF7fjkMFLfBzkyPYUwmZhPDhniK4qKW4VuuJ/6ZV4W//yUFiiAZquQSlNTrEBGrw1q5ihHop4OMqR9IoF/h5usDXVYZQXw1unxOHjd/m4sv0ClQ2WTAu3AkffHsa4T4qfPRNHuakhmJKnDdOV+mQVVyP5jYrXt9RhD/MDIWGrAAApVwKF7UciaFStBotEELgf+bF4z/7CiAhwsHMKkQGOGN8pBuOnW6BRiFBcV07lLKO629MZoHHlySg1dAOr7Za/DfPBHn2EUQFxUIpk6C62YS4cF9s2pcPqYQwbUwANMpaLJsZg7teOwQACPaUwGi2wdhuQ0qkO4ztVjipZFAppAj0cUVbuw1p2TW4KtEHB7PqcNUYX3i7qSGXSbH7SBm8XBU4kN0EbYAKTiop3DRyFFW1QSWXYNb4AAT5uOK1bdmYmeSNuHBvbP3uNHzdlXDTKLAzoxZmq0BtkxERvir8VNyKUC8FcptsaGy14k+zvLEurRZ/8CnFt2UapES6IrtMD09nOWpa2qFRSBAb5IzaFhOqm81w0Sjg6STDBK0HjuQ34nStCR5OUuiNNiy8wg+Jkb44U29AcW07JBpAISXEBTujqsmIGyaNwu6MM0jLLEdymAsmxQfgUHYVFk0OQk5JA2KDNNAoZAjy1mDjd5WI8FVAZ7DCZBGYpHVHVLA7fsyuhlwmQZifM44WNqLNZMMDC+JRUFYPXZsZOoMZ0+K84OOuRnGVDrEhUnx5tA7J4c6w2myQSgjJoz2hM9Tgt6nB2PRdKVbM0eJQdiXMVhv0BgueWhaPjd/mIeN0C7QBauScaUOghwJRgc5o0DdhtJ8Kni4KBHg54XBeHYztNgR7KtDUakFSmAuKqtqwfGYkdqeXYbSfGjXN7citMMLXVQZPZzla2iz45JscrF4Qj5ziWrQaLZg9MQKAQHFVM4oqmtBokeGpNw7hlSVh2PJjCXQGC2amhOKu1w7hrfuuwM7idjS1VsKpqBbPZ2swO9kbPu5q7DtRhapmM3QGG1bNCsOEmEA8uf4IfN0UiAtxRU2TEQdzCwEACybqcDC7AXU6C57/wzh8mamHk0oOuUyKu9cdBQBoFBIkhGpgtZUiLMANxiYLXnn/MP7nugg8+8lxRPiqsFNvhldOJm66NgYni2oQ6O2M7zMrER7gghNFjbhu0ig8+0U9rk6wIrNEj+hADdydLSiraUFLqxlVjS1YOTcA//P6IVyT6IGyOgP83ZXwSw3D8x/+BABoNZqRXabHuEg3RNgakKdTo6HFgAlJEUP7pcr67UKulldd6OhrI2V0ucHklZCA1qrqfpV95vPT+NNNyYNcowv3xo58zJ3g36+yeoMF1M/hJa5OcEf+mdZ+lfVzU/ZdaICdKNHBWSntV9mX7kzFfW8e6lfZw6f0mDnGo19lCyt0/So3kPYcq4RG0b/OvvXflODqBM9+lc0tqe93HeamBuBEUUO/yw+Er0804oqo/t0yfUWMN7YequxX2c8O1cDbpX9fwW0mS7/KMdaJu+UZY4wxB8PBnTHGGHMwHNwZY4wxB8PBnTHG2KAhIn8i2khEhUSUTURfEJGWiAz2vO05RHSYiG7rsswNRHTCno89nYimXsp9GIl4+FnGGGODgjpytH4G4H0hxFL7tLHoSLFaaM/bDiKKALCFiCRCiPUAvgWwXQghiGgMgP+gIwc76yduuTPGGBssVwMwCyHe6pwghDgGoKxrISFEEYCHANxvf60XQgj7bCcAAuy8cHBnjDE2WBIAZPSz7FF0aZ0T0QIiygWwE8Adg1A3hzYsu+UvcqCcYeFCxwJgjLGRgohWAVjVZdI6IcS6C11d1xf2bHKfEdGVAP4PwLUXuN7L0rAM7riIgXIYY4wNDXsg7y2YZwFY1M/VJQPI6WEbB4gokoi8hRB1F1DNyxJ3yzPGGBssewAoiWhl5wQimgBgVNdCRBQG4AUAr9lfj7ZfjAciGgdAAaD/QxmyYdtyZ4wxNsLZr3ZfAOAVInocgBFAMYDVACKJ6Cd0nILVAXjNfqU8ACwEcCsRmQEYACzpcoEd6wcO7owxxgaNEKICwI09zFL3sszfAfx90Cp1GeBuecYYY8zBcHBnjDHGHAwHd8YYY8zBcHBnjDHGHMxQX1BnTElJebof5cIGuR5Dwlhfh4ApU6Eek4KaXdvgFBQIfXk5ZGo1Mv73SXiPSYDP+BQ0nSqCZ4wWJVXNKKttw6zEMHg552LNbSl46r10WK3AY+uPIjpQjVAfDbJO10Mhl8BqE6htMcPHVQ4/NznaLTb4uSuQX2mAySwQE6iG0WxFekEjwn2UuH5iMD5LK4XVJtDSasZvJ4ahue0UjO02ZJUb8JuxGlBTLRr0FhRVtaL9x0LMSA7Ge1+fgoeTFEXVBqgVEjTozQjyUuHqBE/kljVDKZdg03elkBDhmkR3bD90Bv5ezgCAqLpsZLcFoanVAqmEMC7CFVmna9GgN2NClAdyyloQeGg7SlOuh5NSirZ2K5r0ZhjNNlwZ74XCSj3e/TIXM8cFYE+hCc1tFlhsBoyPCcSR/ExEB7ngv3vzEB/mgUZ9O2w2gWMlOpx4Jx1XRLnAbLEhMtAFmcXNGOWjQqvRgqpGE+JHueJIQRMSI32QEumKmBAPKBVSSAjIKKjHBK03fNw1+O3EUHz+QwmuindH3plWhPposPDKKHyfWY4zdW34+mglpsSZMSXWEydON2NPZj18XTv+rYqq9PDIq0C4jxJEQFVTOwLc5fB0loOIYBMAefkhwr8N37WOhpezHj+d1sFotqGt3YbEUGfkVbQi43QL1HIJhBCoqGtFdbMZEf4alNWbICFCTJATWo0W/JjXgEP5DfB2UWD+RD98c7wWja025Fe2wtdVgT3HKlBYbcLd80ah6EwDTp6ug9FsxbZDFZg9zg8Zhc0I8lRC2mhATKAKNgH4uMrg46qAVEII9Ow4fgCQW9aC2hYzvF1kaNIZ8M2xGowNd4UQgFIuxeQYN+Sd0aPdYoO3ixRF1W34TbIf9p6oga+HE6pbzNC1tcNqEzBbbMgq7dhvIsJf3s/AuAg3KGUSnK4xQiohRAc5Y5zWH1sOFsHXXQmzRcDYboGXiwL+HmqM8nfF+98UoazOAD93Bfy8nFHbYoabRorxo91x5dhQ7D5UBIVcCv92K1yd5PjLByegDVCiQW9B3hkdgr1UkEokMFttOHC8DF7OUhToJVgw0ReNuna4aJSI8lfipU3HMDbCHSaLFcqxE/HKlSq8vzsHZ+rbEOSlwfhwV+gMZlhtNmw9WIAQLxVc1DJkl7XA01mOCZHOqG5qx4OLkrD+yyzIpQRpYzWqGk1ICFPjn9uzEOmrhK+7EhNjffHPnYUI9NQgyMcNb247gbGjnFBU0YzfjvfDtkNV8HCSwmi24bMDBXDVyPHtsWpcmeADY7sFYyM94e3mhKVTA7A3sxYuKimsNgEnlQzuziqU1bbhilgftBrMeHlxKCqsauSWF6G4xoDvTpRjYpQbvjnRiLI6A8ZFuuGalHAcyT4DicSGD/eVYkJSxCX7fmW9G9Lgnp6e/lx/yvXzB8CI5D12LFpKSvtV9u8bfxrk2lyc6iZjv8sS9V2mk0za/8Jq5dDf8PHlj6fh6iTvV9maJkO/1xvqrUJZvalfZWub+1duIBVUtCLQU9mvsj/ktSDAvX/H6EhBY7/rkBDmfl7HdCDIpP3v4Jw9IQwb9+b3q+yBYyX9Xu/M1HDYbP0uzhh3yzPGGGOOhoM7Y4wx5mA4uDPGGGMOhoM7Y4yxQUNE/kS0kYgKiSibiL4gIi0RGYjoJyLKIaLDRHRbD8tOICIrES3qMs1KRMfsj+1DuzcjBw8/yxhjbFDYk798BuB9IcRS+7SxAPwAFAohku3TIgBsISJJ5/jyRCRFxxC0u7ut1iCEGDtU+zBSccudMcbYYLkagFkI8VbnBCHEMQBlXQsJIYoAPATg/i6T7wOwGUDNENTT4XBwZ4wxNlgSAGT0s+xRADEAQERBABYAcT5Q9QAAIABJREFUeKuHcioiSieiH4lo/sBU0/EM1275/g52M2ylp6c/fanrwBhjg4mIVgFY1WXSOiHEuv/P3n2Hx3XVif9/f6Y3adS7Jdly707kOHF6SEIIkEIIoSwQYAk8+4VdwvLlyxbAwC678NuF3ZBls14eCLCUbEgzJCSk9zi2Y8e9yE299xlNvZ/fHzOyR4rGlhLLdpTzeh49uufezz3n3JFmPnPrebPVZUz/G/D/VDUpb3xIRrWqtqYP5T8lIjtU9eCbbHPGOiuT+2QfdmMYhmGcOelEfqJkvgv44AmWZ1oF7ElP1wO/TSf2IuBaEUmo6oPpIWRR1UMi8kx6PZPcxzGH5Q3DMIzp8hTgFpHPjs4QkdVATWaQiNQC/wL8CEBVZ6tqrarWAr8D/kJVHxSRfBFxp9cpAi4Edp+G7XjbOSv33A3DMIy3P1VVEbkR+DcR+RoQAY4AXwLqRGQr4AGGgB+NXil/AouA/xIRi9TO6T+rqknuEzDJ3TAMw5g26cPoH5pgkXeS69+aMf0SsOzU9GxmM4flDcMwDGOGMcndMAzDMGYYk9wNwzAMY4Yxyd0wDMMwZhj7unXrprTC+vXrL7vtttuemY7OzDQPHkqsa0r6eGxLC23eUhqLFpAsqWR70XI+8H8+SsXaC/nyQ53EV5zPQOUCioIegn4nv3u5lT3tMa5ZVYiVtIglLdoGktSVejh3fimgNHaGaOgYYV6Zl5GYxabDEY70JCnNsbF2YQGWleSqcytJJJLsbwvTOZjggU091BQ6mFvm59k9A9iJseXQMNefX862w4P81Xvn8a+PHOLrNy+mMM/HirmlPLbpKDubR8j32znQGScWt4gmlHnlXp7fM0hlvotfbxzitisrWVabw/6WIcryPSSTSTYfHMJZUc7qBSU8tq2Ha84p5rWDA6yYU4CVTPLivj5GYhZ5K1dit9koDrr4w7YhCn1CTbEXh81Ga2+Uj181jz9uaqKuPEBloZdWVyEuhw2vS7DbBKfDxrM7uinOdZFIJPnMtYtZMiuH7oEw19U5Ka2ppLVriOGRJJ96z1LK8900dgyzuDqXjt4Qzd0j/GlbN4V+G6FIgubeGNddUIPTaWfP0R5K8z088XoPq+bkcOeTPcwvcbCvJQQCS2pyGQjFeGpHL0W5TmIJi8uXl7KopoClsws5Z3Y+//1kM+9aUcw58woZicQIeB1ctKKaulIP7T3DXLKimu0Hu+gLxVk1J8iRzhFyvXaurq8iPBJlR1OEy5bk0dAe4c/eNY99jb1cfW4lB1oG2dyU4NKFOURiSUqCbv7imrk8t7uLK1bNYmh4hP5wnHWfqOeRTc187cOriEdHCEeivNbQzc3nz6KpJ0xzT4QPXlLHz55q4XufPofob/6L2ZdfQsnv/4vc8y+kNN+LpcryOfmsXVpJW/cwkViS962pIhKLU1Oay+YDvVy/tpbkz3/IC4461i7MI5G0mFsR4NLl5QyGIly9eg7nzCskJ9TNE/tD2CRJYY4Lv8fOc/tH+FTLfdwXmc+tV1TSOxylssDNLbOibOkWFlblcrh9AJfDRmWRn4tXVtPVF8ImsKyumO0Hu7i5KkJBdSXRWJKNe9q5+bI5fP/3rdx+/Tzae4a4ZMUscn1u9jT2YxchnkyS63UCcO68AoYjCQ60DtM1GGPNwmKaOkNcEujn3l1hBsJxDrb0MxKzWLuoCI/LwX2vdFJT6OIXTzRwy+VzeXRzKyV5Lp7b1YXLYeOV/f08vWeY5r4YQ+E4eX4HdeU5PLy1l9tvXEhZ5wE85RXct6mP+bV5eFw2Ht7cRUHAgc9tZ3F1kF881UTnsMXLO9u4dnUlz77eRl8ojt9tZ3ldEb96oZNVNV7y/E6Kcj3sPDpIbYmPTQ39qFpcd2Edv3+xgVcP9HPxkiKOdIZ535pq7nzkKIPDI1QWehkaiTMcjtI0mGT5vFK8TqV7IEI0nqSlJ8qOljjFARvP7h2iobGbjv4RXA4bH/A2UrB46bfO7KeskY3Zcz/Nrqmvzrps0/6eMeXvP9yRNdZue8NTm7La3JiYdOxQJDnp2PesKpx0rKXZlxXnjb1otqEzljX23AXlY8prl8/KGnvO/NLJdQ74xBXZ6xmvIMc16dj2V1/JuuyHG8Y+d8OysteTk05Co+5/sTFrrMM++bd1VXFg0rHL6ib/ep7o790Q9Ywpl+V7skTCg21jly2sKcoa++87Jv+e+PC7Fkw6dkGFf9KxeTnZt2W8h3rysi7rGBj7Hmjoyv6+vPU9Sybd5sDI5N/fTodJD29n5q9nGIZhGDOMSe6GYRiGMcOY5G4YhmEYM4xJ7oZhGIYxw5jkbhiGYUwLESkTkd+KyEER2S0ij4jI/BPErxSRaydR7zoR+cqp7e0b2xeR69LPxJ9KHT8VkU4R2Zll+VdERNMD30wbk9wNwzCMU05S47U+ADyjqnWquhj4W+BEt12sBE6a3N9iv040psqY9lV1g6pOdQjyu4FrsrQ9C7gKyH67yylikrthGIYxHS4H4qp61+gMVd2mqs+LyC9F5PrR+SLyKxG5Dvg2cIuIbBORW0SkQEQeFJHtIvKKiCwf34iIfFZE/igiXhGpE5FHRWSLiDwvIgvTMXeLyA9E5GngeyJynoi8JCJb078XiIhrgvZvFZE7M+q4Ix1/SEQmHKdeVZ8DerO8Jj8Evgqc4GbRU8OMCmcYhmFMh6XAlizLfgLcDjwkIkFgLfBJ4BtAvap+AUBEfgRsVdUbROQK4Bek9q5JL/8CcDVwg6pGRWQ98HlVPSAia4AfA1ekw+cDV6pqUkRygUtUNSEiVwLfVdWbRGR8+7eO63c5cBGwENhAaqz5SUl/eWlR1ddTBzWm1xlN7vX19V8jNZbvjLN58+Z1Z7oPhmEY00lEbgNuy5i1XlXXn2w9VX1WRP5DREqADwD3pRPt+NCLgJvS6zwlIoXpLwMAHweaSSX2uIgESH1JuDejHndGXfeq6uhTfILAz0VkHqm96LFPicruQVW1gN0iMumnOomID/g7Ul9EToszvefuMUnQMAzj7SmdyLMl813AhIeu034JfAz4MPDpLDET7eKOHtLeSWovvgo4TOo0c7+qrpxgHYBQxvR3gKdV9UYRqQWeOUE/M0VP0rds6oDZwOheexXwmoicp6rtU6hn0sw5d8MwDGM6PAW4ReSzozNEZLWIXJou3g18CUBVd6XnDQE5GXU8R+oLACJyGdCtqoPpZVuBzwEbRKQiPf+wiNycjhcRWZGlb0GgJT19a8b88e2fEqq6Q1VLVLVWVWtJHXE4Z7oSO5jkbhiGYUwDVVXgRuCq9K1wu4B1QGt6eQewB/hZxmpPA4tHL2hLx9eLyHbgn0mdl89s4wXgK8DD6VvLPgZ8RkReJ3Xk4Hom9n3gn0TkRcB+gvanTER+A7wMLBCRZhH5zJup560604flDcMwjBlKVVuBD020LH0eeh7wm4z4XmD1uNA3JGhVXZcx/RjwWLrYzQS3oanqrePKL5O6wG7U10/Q/t1Z6phw1CVV/chE88fF1J4s5q0ye+6GYRjGaZW+Qn0v8CNVHTjT/ZmJzJ67YRiGcVqp6hNA9vGvjbfM7LkbhmEYxgwjqWseJq++vn7dqbp97VTWdTb6Pz94RmcVOjnUGcPrEuIJJeCxEYoqy6p99IfiNPXEsQm8b3UJy+tKKcj1cvtdG1m7MMjBtjA5XjvNvTGSlhJPKN/55Cqef72Rp3f2Ylnw1ZuX8P17d1EQsNPWn+DmtWU88Eo7X715MUfbB3jq9Q7a+hPMKnDS1h8nnoSiHDvhmMXHL6/hlT2dHOqMMByxuOMvzuMnf9jJx65aSCSa4F/v28l3bq3n9rtexVIIeGzk++wUB13sagoTTyrf/sQK7nhgFwUBB5GYRdJSOgcTRBPKl2+YT3WOELZ7ufuxvaxdXMy2hh4uXFrG/qZ+RKCtN8L719YSDHj44e+2c8nSYo50DFNV5GPjvl4uXlJEWWGAXYd72Hp4gAsWFGC3CaFIgotXVLFxdytzK/No6hyio2+ES1dWsXF3Gxctr+JAUw8VRTl4XA42vHSYOWUB8nLc9A5GmD+rgIHhCACPv9ZGdYmXZbML2d/UT47PSVmBH7/XSXtPiNcP9eJx2YnFk8yvCtI9MILHZaepa4RowqKq0MP8WfnYRGjvDTE8EmdWSQ4Br5OO3hCLaov51ZP7WTE7D6/bQSSWxOmw0dQ1zHsvqKOxvZ/BUIx4Ikln/whHOsNcU19FImlhswl7G/twOWzUluXS3hvC6bDh96Ruy31yWwcel40rVpTh97ro7AvR3hsmYSnLZheyt7EPW/qe3/dcUMfDLzXg9zho74tgS9/Is2x2AQOhGN0DI3jdDroGoty0yIcnvwAA9eXQ3R8iJ9RNIhwmp7qGzq2v4V+1hoDXRfdAmNiWFyleuZKBQ4fwFhcTbm8nHgohK86nMuhGVRlKCC6nnb6NL2Bfuhqnw0Zj+wD5OR6GwjHKCgMU5HqPvX96BkYo9AgdoSTSsIOS885Hkgm6d2wn2teHbfkaygv9JBUirc1Y8TgA3urZOGxC85OPU/muq7DCIew+PwCRznbILyZycB/3H7Zo6o5y6dIiXtnbQ1m+m/a+KH5P6j0XiVksqPAypyzAw1u6sAnEk+C0g90mJC1FRAj6bHQPJXHYIGFBvt9ONGERi6eWx5OKz516sZNJyPHaKMpxcrQ7SjSuWAqOcbtZloLHJYQiitMOIoKqUhp00jmYIGkpTruk1rVDIgkluQ46BxM47OC0C+GoRcKCa88pZF/LMIc7o9htgk3AbgeHTXA5hKDPweHOKHWlbpp7Y1gWxJOK3SZ4XcJA2OKC+Tm82jCEpeB2pLfFUn7815dP/9NYjDfF7LlPs7b++JhyLHn8y1RTz9hl3/rV9jHl3lDi2HQiOWYRs4uzP/vnrz+wcEz54oXBLJHQn9EGwKHOyLHpr35o2ZhlVQWuMWW77fj7ujBn7LITqSoZ25/2nuEx5c7+47eShqNj+1ec5x1T9riOX+i6pLZgzLKdR/qz9uFo+9jTfAtn5R+bHh4Z+3epnz92fAcr4/twbYlvzLKS/LHlpXUlx6brqvLHLAv6xr5mB1oHj003tIztu8Oe/TP02vOqxpRry3KzxlYW+ceUXc7sHwFDTWMff+30H1+37MKLxyzLv+jyMeXC8y44Nt383LNjlpWsOmdMubTw+J1HOw91jlkmzrHPFolnXNjsOLp3zLKRvDKySYaGxpTz5mUdu4SG9siY8qYDfcem7Sf5xBxNfKnYsX+z+Nh/5THv6fi497fPPbahzKoc9rGxJ3rYWa53bD3xZPadud7hsR3MbGfzwbGvn9tpcvrZziR3wzAMw5hhTHI3DMMwjBnGJHfDMAzDmGFMcjcMwzCMGcYkd8MwDONtRUTmi8gjItIgIntE5H9PNEqbiFwmImsnUe/d2cZpfyvGty8inxeRT5zqdjKZh9gYhmEYbxsi4gEeBr6sqr9Pz7scKAY6sqx2GTAMvDSN/XKoaiLL4jHtq+pd09WPUWbP3TAMwzjlRKRWRPaKyM9FZLuI/C79PHlE5BsisklEdorIekmPgyoiz4jI90TkVRHZLyIXT1D1R4GXRxM7gKo+rao7ReR5ETk25KuIvCgiy4HPA7enB4S5WERqROTJdL+eFJE3PC1PRL6T3pO3ici5IvKsiGwRkcdEpDyjv98VkWeBvxKR94vIRhHZKiJPiEhpekjZ8e2vE5GvTGGbp8wkd8MwDGO6LADWq+pyYBD4i/T8O1V1taouBbzA+zLWcajqeaSGg/3mBHUuBbZkae8npIdwFZH5gFtVtwN3AT9U1ZWq+jxwJ/CLdL9+BdyRWYmIfB8oAT5FatS4HwEfVNVzgZ8C/5gRnqeql6rqvwIvAOer6irgt8BXVfXIBO2Pd7JtnrIZd1i+vr7+a0D2J7ycJjP5yXuGYRgAInIbcFvGrPWquj6j3KSqL6an/wf4S+BfgMtF5KuADyggNTzr6J74/enfW4DaKXbpXuDrIvJ/gU+THtFtAhcAH0hP/5LUELCjvg5sVNXbAERkAakvFI+nDzDYgbaM+HsypquAe9J79i7g8CT7/Va2eUIzLrkDHpNYDcMwpl86ka8/Ucj4cvqc+Y+BelVtEpF1jN0hG31EZZKJc9Qu4NIs/QmLyOOkhon9EFB/0o14Yz83AeeKSEF6CFgBdqnqBROvSihj+kfAD1R1g4hcRmo8+sk42TZPmTksbxiGYUyXahEZTYofIXXYejSRd4tIAJjq1em/BtaKyHtHZ4jINSIy+rzsn5A6zL4pnZwBhoCcjDpeAj6cnv5Yul+jHgX+GXhYRHKAfUDx6HaIiFNElmTpWxBoSU9/MmP++PannUnuhmEYxnTZA3xSRLaTOvz+n6raD/w3sAN4kNSe8qSp6gipc/RfFJEDIrKb1Hn2zvTyLaTO7/8sY7XfAzeOXtBG6vTAp9L9+jjwV+PauDfdxw2kDsN/EPieiLwObAOy3Va3DrhXRJ4Huk/Q/rSbiYflDcMwjLODpaqfHz9TVf8e+PsJ5l+WMd1NlvPPqroXuGaiZSJSQWrH9U8Z8fuB5eNCr5ig3lszpn9K6uI5SCX0S07U33T5IeChCeLGt/98xrLLMqazbvNUmT13wzAMY0ZIPxhmI/B3qmqd6f6cSWbP3TAMwzjl0reALT3Nbf4C+MXpbPNsZfbcDcMwDGOGMcn9NHPZZdKxBf6z+8BK0hp/l8upUZLnnnRsJJaclj6ciG3yf0J2HuycdOy8itw30ZvTJx4KnTworefVlycd29EzNOlYJ6f/7716Xv6kY6OJyb8nHPbJ92Ga3mrGDGaS+zTqCSVp6U/yeptFf9iiO2RRV+YlllCqinwsr/HRF7YYjlqsmFtK26DFSPNRDnRZdA/EaOyJsaQmyIuHE2xuThJPQmgkxryqfOw2oX/E4khbP5bC9uYYvWGLSCzJvsZeBoYj7G8eoKrQw77WEMVBF4d7LCxVKgvcxBOKqtLaF8NSpW3Qom8wQvtAko6eYZwOG/uO9hCOxNnZbhGOKbk+Jx6XDbtNiMSVUEzZ3tDJcDRJcdCD3S60DyTY25FkXpmbnYe7sdntRONJbAIuh51cvxO/18n5SypYXFvItWuqKSsMAHDzJTWsnFfKstkFNHeH+eDFNRTn+SgvyqGuMsifv2c+K+aVsXb5LK46bzbO7lauWVNHcb6fORV5tPeG8LgcXLJiFm6nnWvPr8HrdlJkj1M/vwhLlfqFFeT4XMypyCNpKWWFOVyzupL8HA9zKgtYNqeI6y+YTVVJLkV5flYvrmBOZT6ffPdi6ipyqCgK8O41c1gyu4h3nVPJmsUVfOD8GpwOG36vk2TSYtW8EsoLA7icdi5eVkk0nuDa82bR1RfG6bAT8DopLwyQtCwCHifzCpwA5OV4WFhdwEVLSoglLHL9bv7n6aNUFvmZV5XPnqN9FAW9BP1uygsD5Od4WLOggIHhKPGkxbyKPOoXVeJxOZhVHGDXkV7mlOeybE4RN6ydjfZ143E5WF5XwkevXMQN8z28Z80ccnxuaspyKS/0c8W5tVxxThX2slkMFc4iGYvS3R8iEktgFVWgVXNoD1s8u/UoTR0DWKrEE0m6+oZJBPKx6paSUzuH9qI6ButWUpIfYGegjlcburHZhEMtvdj8Odzx4G7yA262H+rB53HxuxcaufOh3Ty5tZWu/jChSJyfPbqXxzYf5Z5nGiisXwPA4NGjtOVWkVNTS7J2EXsbe4jFLZxeH/c8fYBA7Ry27G1lIBRlYM5KfvrITmIODz1xO8RjhDxBDjb38tKedkrzPDT1xnl4Swc7WqK82jDM/o4YAPGE0j5kkbSUfS1DtA1atA1aRBNKy4BFjtfGvk6LCxbkEopaVOY7yfPbWVzppXMwQVHAid0OhQH7sYSftJR4UtnREqcs30M0riyd5cPrEgYiFiNxpanfoqnfIui143WmPp6b+i38bsHlFHpDCYYjFjaBLc1JmvuTJJMQiVvEk8rAiMXgiIWlqS8aFfkuWntGUFV8bhsdQxZBn52+kMWRngSRuBKJWQxGlGhCWVWbw3DUomPIIpGEg10JtrenPrsO91gc7LZo6EqSsNR84TjLmeRuGIZhGDPMtB/3PcnjYGunu33DMAzDeKc5HSd1sz4Otr6+fsL5hmEYhmG8eeawvGEYhmHMMGfN5dincDS32lNQh2EYhvEWiUiS1GNmR/1WVf95grj5wL8B84F4ep0vqmrHJNu5G/iDqv5ORI6QGpSm+8RrZa2rArhDVaf6zPuzylmT3DlFo7mZQ/2GYRhnjRFVXXmigPQocQ8DX1bV36fnXQ4UA5NK7qeKiDhUtZWpD2Zz1jGH5Q3DMIwz6aPAy6OJHUBVn1bVnSJSKyLPi8hr6Z+1AJJyp4jsFpGHgZJxdX4xHb9DRBam1/GLyE9FZJOIbBWR69PzbxWRe0Xk98Cf0m3uzFh2v4g8mh6k5vu8TZxNe+6GYRjGzOIVkW0Z5X9S1XvGxSwFtmRZvxO4SlUjIjIP+A2pMdpvBBYAy4BSYDfHB3kB6FbVc0TkL4CvAH8O/B3wlKp+WkTygFdF5Il0/AXAclXtFZHacX1YCawiNeb6PhH5kao2TXL7z5h3XHI/hef2T+hUnGIwDMM4m4nIbcBtGbPWq+r6jPJJD8ufhBO4U0RWAklS5+QhNULbb1Q1CbSKyFPj1rs//XsL8IH09NXAdSLylXTZA1Snpx/PGPt9vCdVdQAgPbxsDWCS+1nolJzbNwzDeKdLJ/L1Jw3MICJrgP9KF78B7AIuzRJ+O6nz7itInUaOZDZ/gmai6d9Jjuc5AW5S1X0T9OdEz1aOZkxn1ndWM+fcDcMwjNNGVTeq6sr0zwbg18BaEXnvaIyIXCMiy4Ag0JYevvXjwOgT+Z8DPiwidhEpBy6fRNOPkToXL+k2Vp3CzTrrmORuGIZhTBeviGzL+HnDbXCqOgK8j1TiPZA+9H0rqfPtPwY+KSKvkDokP7qH/QBwgNQtc/8JPDuJvnyH1GH+7ekL5r7z1jbt7Pa2OLxgGIZhvP2o6qTGvlPVvcA1EyzqAJZnlP8mHa/AF7LUVZsxvRm4LD09Anxugvi7gbszykdIj0M/wbL3nXRjzhJmz90wDMMwZhiT3A3DMAxjhjHJ3TAMwzBmGJPcDcMwDGOGmYnJPVJfX78u2w+ncWAZt0MIeoQSPxQG7PhdQt9wnNJg6jrGvuE4Jbk2QjHF7XJQV+zAVVHN4nI78ytzyPHaeHxbF4tLhMUlgscplBf62bS3E0uVHLfwyt4u4knFbgOPQ9jc0MeRjmE6+8J0DcQIRZIc6Y7jsNvI8ZBqy2mjO2SxfG4ZlkIsrsQS4Pe68LuEHL+beMKioihAR+8wxX7wOoVQJEE8oXicdopzU9sTiyc52JWkdzjKliMjVOY7OafaSSRmAWB3uWnuHKSy0MtINMHciiDxhEVoJMY9zx4mEksAMByO8fz2Npo7B7HbhKDPSd9QhF1Hemlo6iHod/PSzlb6hyKERmI47DbanXnEkxa5tiS15Xn0h2IUeGzY7ULQJWzZ20EBEcjJ40DzAG6nneauQVbNLyfH7+ac2QUUOuIcbBmgOD9Ae88QSctCk0mCfhd+m0WktZl8n5NoIrU9cwrcbHihAVXo7Aux92g3iWiEQ60DOB12+oej7DnaA8DAcJS+cJwntjRRkOsl6HcTTyTp6h9hdnmQaNxib2M3jtw82nvD7G3sY29jL26nnYpCH6/t7+AL1y3gcNsQvYMRhkbiVBbnUhj00js4gt/rIhpP4nE5iCcs7nv+AG3dg7T3hth1pI/FNQXsa+4nL8fDb59poM/mw+2yIyKEIzG63QV09A7z4EtHAGjqGiYWT+K027EN9TKrNEjcn8fOQ12UFATI9bnweVx4XA5EBL/XyXPbjlJemIPdbiORSFLoiBNqaaK6LEhNeT53bdjOZefUpOoN9ROLJxHgfavLSYbD1JYGiCeSvH9NJbdcWsNly8uIxhMEPA78Hjt5ATfLaoMkerro6g+zL+alsWOQdmcepQV+CnJ9uJ12joaFc+YWEG46iipsePEgS2oLuWR5Ba54mKKgjxG14xno5LUDXaxeVImlkLSgqsBNXbGTulIXXqdwsDv1t7YJnLegiIKAk/JcG+EYxJNgt0FJrovqfGF/S4jSXCc7WqI4bMJLDamLuY92x8j3Odh4OIalMDRi0Tmcqtfvgn0tw9TX5dA5EGNnaxJLU58RQY8QiUNvKEHPcJLSoIOhKDT3JQl67cQTit0mDEaUc6vsDERgMGJRGnQQ9NkpyrHROay83pzAbhN2NY9QUeilutjH0IiVes8mLLpDypIKF0lLSVqK2wEep7C/LUTQZyNhgdspnFPtxm2HoQiU5QpVeUKhX2jtT1KUY67HPpvNuL/O5s2b33CrRaYzPbBMY3eMgsDkXvaeYYvS3EldbHpGlORO/t/n9SODVBT6J1dvvp/DbYOTirWlblk9rT5w6XzauocmFZvvd0263vecP4fHNh6aVOyR9iGK8yf3ep4qoZYmHKWVk4rdur+dNUuqJhXb0NyL3zu51ynod9PZF55U7KnS2hs9eVBaXbGT9oHEpGI3HhiiPG9y76GA+/T/nxtvbzNxz90wDMMw3tFMcjcMwzCMGcYkd8MwDMOYYUxyNwzDMKaFiAyPK98qIndOsY6rRGRLemz2LSJyxant5cw04y6oMwzDMGYGEXEA3cD7VbVVRJaSGgBmcld2voOZ5G4YhmGcdiJSDNzF8THVv6SqL4rIOqCC1G3L3ar60YzVdgEeEXGr6uRvY3gHMsndMAzDmC5eEdmWUS4ANqSn/x34oaq+ICLVpPbIF6WXnQtclB7sJdNNwFaT2E/unZjcI6fjXvfNmzdPexuGYRhnkojcBtyWMWu9qq7PKI+o6sqM+FuBUEwwAAAgAElEQVSB+nTxSmCxHH9WRa6I5KSnN4xP7CKyBPgecPWp24KZ6x2X3E/2kBvDMAxjctKJfP1JAydmAy6YIInD8XHbR+dVkRrD/ROqevBNtveOYq6WNwzDMM6EP5ExJruIrJwoSETygIeBv1HVF09T3972THI3DMMwzoS/BOpFZLuI7AY+nyXuC8Bc4Osisi39U3Laevk29Y47LG8YhmGcHqoaGFe+G7g7Pd0N3DLBOuvGlf8B+Ifp6uNMdaaTe+bFbbVnsB+GYRiGMWOc0eSeeXHbmR6tzTAMwzBmCnPO3TAMwzBmGJPcDcMwDGOGsa9bt25KK6xfv/6y22677ZlTHT/Vet8OjjZ1rCvOdfBnl1bR0RtmWXWAgXCcBRV+Xm0YYPW8fGLxBLG4RWN7H3k+B+cvLuPxrS209UW4cEGQfa0jeJw28nw2LlyYTzJpcf8r7cwv9+GwwycX2WkYhNVzclheE2BoJE5Tb4ziHCd5fgevHw3xycur2HaonxU1Pgr9duJJZXGVDyuZAE3w/vOrqMoT2rqHSFpJBobCvLq3kwdeaccpSRxiMafEQzRhcagzRlmei7b+1O+BcJxFFW5ae6JcubwAUHI8dnpDCS5ZVsZrB7tZW5fH4a4QXQMRNu7vo3dgmIrCAI9v62LHkT7C4Qgv7OxkJGaxr3mA/c39lBd4ePS1TvxuOxcsqeDVPW3sbhrmlb2d9PSHeOTVJg62DvLb5xqpKffz7V/v5OpVZTz0ShMbNrbQ1DXIi/sGqa7I4ed/OkDPcIxwJMET2zp4dnsbHT1DPLerk+b+CD2DMXoHQ0SiCR58uYUHNrbR0tHHuYvK+YcHDtDZP8yGlxvpGYyys2mA3Ue62bCxnXg8QTSe5IFNHexuDnPh4mIef62FwXACp12JJy1+8/RhioIuNu3ror03TDxhsf3IIPULSnlsUzMBj51XdrfROxyjvMDLi7u66B+OUpznJZaweHFXB92DcTr6Rzh/UQk7D3ezuLaI1w50UlUU4BdPN3Hx0mKe2tbG1aur+cPLR/F57AyPJFhQnc+jW9rpGwhxdcEwe4fsvLy7k9ULS/nVE/t4dV83+5oGuP2Dq3jqtUZ2NQ3T1TvEho2tXL12Hs2dAxTneqhwRNh0ZICCsmJyXDb+tLmRNUsq+PWTB1hUnU9bzxDNX/8ym4JLKDy6jUTtIjwuO8OPPUh1/Sr67/sf7Ls2sdFVy4q5pfQNRTjSPoDX56GjN8ySZAcvt1vUluby4s52Ll9RxSMbD3PpyioW1hTR1DGIw++j2BYjjJ3yogCxxx4gb9lKNu1ppTjPR+89P2OLllEb62DOisXMtQ+x5Z++S+3yhewftlOe58FhJbDlFbKs1MPe1kEWVBfikhg2gb98Tx2VZUG8DotrVhVSGHAwPBLn3fVV9A+NMLfcR12Jk1yvDbcdBKgpdmOzCTabcNXyQrwuGyVBN/Vzg+xrDTESs/j2RxdCMkpRjpP5ZW6qizzMKnRhtwkrZ+cwFE5w5dJc2vui1NflkOezUZZr4+pVJVy+opSDrUM4xMJhgwUVPopyHIRjSSyFklwHq6o9XL6sEFWLSNzCLiAoRX7B77Fx+ZICOvoi7GwK8YkrqllY4eHGi+awZV8nc0s9vKe+ArdD8DihKNdFPKF8/c/OZdOeNq49t4RXGwZZVe3mvDo/Qa+Ndy0vIjQS5Ws3L8bzn99g3oc+9K0z/TlrTMzsuZ9mK2fnTjr29SND09iTt87vsU86tqM/NunYzoHJx9YWuSYde6p0DcQnHXtNfcWkY/+0uWnSsTUlgZMHnWKeoe5JxzZ3Tf5/1+HzTTp2X2PvpGNPlY6ByT/ptL0/MunYkWhy0rFLZ03+NTIMMMndMAzDMGYck9wNwzAMY4Yxyd0wDMMwZhiT3A3DMAxjhjHJ3TAMw5gWIlImIr8VkYMisltEHhGR+SeIXyki106i3nUi8pVT29s3ti8i14nI16ZYx+0isktEdorIb0TEk7HsIyLydyKyUEReFpHo+O040fpTYZK7YRiGccpJauzWB4BnVLVOVRcDfwuUnmC1lcBJk/tb7NeJnsw6pn1V3aCqkx4mXEQqSQ+Io6pLATvw4YyQa4BHgd503L9Mcf1JM8ndMAzDmA6XA3FVvWt0hqpuU9XnReSXInL96HwR+ZWIXAd8G7glPfLbLSJSICIPpkeOe0VElo9vREQ+KyJ/FBGviNSJyKMiskVEnheRhemYu0XkByLyNPA9ETlPRF4Ska3p3wtExDVB+7eKyJ0ZddyRjj8kIh/Mst0OwJv+EuEDWtPrC6kvD6+paqeqbgImuq92wvWnyiR3wzAMYzosBbZkWfYT4FMAIhIE1gKPAN8A7lHVlap6D/AtYKuqLie11/+LzEpE5AvA+4EbVHUEWA98UVXPBb4C/DgjfD5wpar+NbAXuERVV6Xb/K6qxiZof7xy4CLgfcAb9uhVtYXU3ngj0AYMqOqf0otXAa+rqmZ5TU62/pSc6VHhTqq+vv5rwJs653Ambd68ed2Z7oNhGMZ0EpHbgNsyZq1X1fUnW09VnxWR/0iPy/4B4D5VTaR2bse4CLgpvc5TIlKY/jIA8HGgmVRij4tIgNSXhHsz6nFn1HWvqo4+OSgI/FxE5gEKOCe5yQ+qqgXsFpE3nF4QkXzgemA20J/uy5+p6v+QOiT/xxNVfpL1p+SsT+6AxyRKwzCMs086kWdL5ruAbIeuAX4JfIzUOeVPZ4l5Q7YnlYwBdpI6zF0FHCZ1JLpfVVdmqSuUMf0d4GlVvVFEaoFnTtDPTJmPK5yob1cCh1W1C0BE7if1heN/gKtJf1E5gROtPyXmsLxhGIYxHZ4C3CLy2dEZIrJaRC5NF+8GvgSgqrvS84aAnIw6niP1BQARuQzoVtXB9LKtwOeADSJSkZ5/WERuTseLiKzI0rcg0JKevjVj/vj2p6oROF9EfOlz7O8C9qSPNjhUtefNrP9mOmKSu2EYhnHKpc8t3whclb4VbhewjvQFYqraQSpx/SxjtaeBxaMXtKXj60VkO6lz3J8c18YLpM6tPywiRaS+CHxGRF4ndeTgeib2feCfRORFUlekZ2t/qtu8Efgd8Bqwg1SOXQ9cBTwxGpe+RbAZ+DLw9yLSLCK5J1h/yt4Oh+UNwzCMtyFVbQU+NNEyEfEB84DfZMT3AqvHhb4hQavquozpx4DH0sVuUue2x8ffOq78MqkL7EZ9/QTt352ljglHb1LVbwLfzJwnIteQuohwNKad1OmESa3/Zpg9d8MwDOO0EpErSV2x/iNVHTjT/ZluqvrnqvrK6WzT7LkbhmEYp5WqPgFUn+l+zGRmz90wDMMwZhiT3KeRx2XjizcsY39TP93DCV7aP0g4mqS1N8KKmhyi8SS9wwlWzQ5w8dJSqkt8NHYM0NyXxOu0MbfMy9UrCrAJrF2Yj8Mu9A6OcMXSPDoHorT1x/jvHQmKgy4OdYTJC7iZWx6gpsjN6kVl3HTpApbX+HnglVZeOhhj+9Ew25siPLsvzBM7BxgKx8j1OensC/Pc7n4uWjELVeX6i+dz40Vz+MxVNZQX+mgfSNA1GONgR5TVdQF6hmIMjVhsOxIiaSktPVF2tkToGYxSHHQTTypNvQl+/cxRtjT088vnG+keStDaG6Usz0VBwM3zO1q5eGEQh0248ZJ5zC33EU9aFAddiAgPvdrF6rlBHHZh8942QpEEi2cFmFvmY05ZgPevqeRz9QHOmZuPx+Xgs1fXcMfDR6gt9bOkys/CWUE+ekkF2w72sKDST9dgArsNVs8NUlPkIT/gwu+2UxBwk7SU7YcHWDy7mIKAg4KAnSW1+exv6uW9q8vZfmSI96+pJOBxsGZhMfGkxcpaP7PL/OT6nNx+01IuWZzH6w2drJyTTziapKU7TDxhsetIL3l+F7OKfaxdUsbGA/209MXZ8OIB5lUEKMn3cdNlC2jqiTIcSTAUSdLYHeFg6yC9QzGuX1vLRy6fw+uNEQZCUY50htm8t42u/gh3/fEAS6q83Pf8EVbPK+CRV44iAnUVQboG4zz00hFW1Obg9zj45ssWxXleKgo8/PSPu1m7uITSoIurziln054W1iwqI5FULlpWQThqsf1fv0dejpeG392Do7ic4ZEY7twgTU8+wQ3n1xBPWMwtD1BeGKA4z8f8D9/CtUUhaq56N66WBrz/+2MGDx6ivTdE/4ED9L7rw5QV+Gho7iXy7B85PyfMvsZeLFWe7HJRUeijqjTI2iVlPL7lKHkBN+UBOy+83kh1aS7J5x9lR1eMRTWFVOW6CM6bh9OeuhOpdzBM4S2f4iOLPQRnz2HL3laGGo/i8HrZFvKysq6E4bgiThd/evUQD2xpY9msPH73bAM3XDyf2tIAhzY8hPfIbq6Z6+el3d1UFPopDToJR+IU57r4w2u9PLlrkCXVuQxHLDoHE2w7GqY06GZBZQAR6BmKUV7g4Zmd3VyyOEhp0Mk//u9eDneO0DUY48KlZcQSFs09UfweO7l+J2X5bjYfHKS2yIXP7SCeVF5oiNLZP8LvX27isuUlzClxket30tEfA8DnsrGyxkfSgqGRBOFoAq/LjsdpY06pj9veXct7zinC7RBe3NdHWYGX0qCTPUf7eHZnNw3NvdSVuukdTvD8zg76QzEGwgkCXgefumYBL+9owmaDl/b0UFfq5tn9EXxuO4/vCtPRN0JxrovP/ngHvpLiM/nxapyESe7T7K7f7xpTXlyTlzV204H+MeWdTaEskdA1ONFTC1POnVc4phxwZ/8zD4bH1pPrzX6mpqrANaa8avaE15MAcOXS4Jiy1zXRLaEpL+9oHlOeV+7LGutx2ceUn9vVfWz6g+eP/bDZ8Gp71noiseSYssuZ/TW6amXJmPLiWdnvlMkLuMeUL15ROek+JJJZAoGbLxi7bYW5x9ux28a+tg579m0Z378TxXoiQ2PKhwcS2eu99Kox5cZ3Z7ttGWz2sX/DeZXBLJFQWzH2/XK4Lfvp2dy6eVmXed1j/6/jtuzPLHE6xr6ec8pzs7fpH1uPw5b9//wnjx0aU97dOJQlEpq6R8aUD7aFs8buOjqYdVlJztj+uR3Z/96bGsbW8+5zy7PGfvPmmqzLjLODSe6GYRiGMcOY5G4YhmEYM8ybuVo+Ul9fv24K8bVvog3DMAzDMN6kKSf3zZs3T3psW4ApfBHI9qWhdirtGYZhGMY73Vlzn3u2Lw1TPEpgGIZhnCVEpAz4N1JPfYsCR4Avqer+LPErgQpVfeQk9a4DhlX1X05xf8e0nx5jfrGqTnqnVkR+SmpI2E5VXZoxfwVwFxAg9Tp8LOM5+aecOeduGIZhnHLpgU8eAJ5R1TpVXUxqTPY3DJWaYSVw7TT360Q7tWPaV9UNU0nsaXczwSNwST1+9muquozU6/J/p1jvlJjkbhiGYUyHy4G4qt41OkNVt6nq8yLySxE59sx4EflVei/528AtowO3iEiBiDwoIttF5BURWT6+ERH5rIj8UUS8IlInIo+KyBYReV5EFqZj7haRH4jI08D3ROQ8EXlJRLamfy8QEdcE7d8qIndm1HFHOv6QiEw4nK2qPgf0TrBoAalR7gAe5+TDv74lJrkbhmEY02EpsCXLsp8AnwJID4e6FngE+AZwj6quVNV7gG8BW1V1Oam9/l9kViIiXwDeD9ygqiOkRlD7oqqeS2q0uB9nhM8HrlTVvyb1XPtLVHVVus3vqmpsgvbHKwcuInXYfap79DuB69LTNwOzprj+lJw159wno76+/muA50z3YzI2b9687kz3wTAMYzqJyG3AbRmz1qvqSYcoVdVnReQ/RKQE+ABwn6omUkfyx7iI9B6uqj4lIoXpLwMAHweaSSX2uIgESH1JuDejnsynNt2rqqOPigoCPxeReYAC2Z9qNNaDqmoBu0XkRKcXJvJp4A4R+QawAYhNcf0peVsld8BjkqZhGMbZIZ3IsyXzXcCEh67Tfklq/PUPk0p8E5nokX+a/r2T1DnyKuAwqSPR/aq6MktdmY/8/A7wtKreKCK1wDMn6Gem6En6lpWq7gWuBhCR+cB7p7L+VJnD8oZhGMZ0eApwi8hnR2eIyGoRuTRdvBv4EoCqjj6newjIfL7zc6S+ACAilwHdGVeYbwU+B2wQkYr0/MMicnM6XtJXqE8kCLSkp2/NmD++/VMmfZQCEbEBf0/qyvlpY5K7YRiGccqpqgI3AleJyEER2QWsA1rTyzuAPcDPMlZ7Glg8ekFbOr5eRLaTOsf9yXFtvEDq3PrDIlJE6ovAZ0TkdVJHDq5nYt8H/klEXgQyBzsY3/6UichvgJeBBSLSLCKfSS/6iIjsJ3W+v3Xcdp9yb7fD8oZhGMbbhKq2Ah+aaJmI+IB5wG8y4ntJ3ROf6Q0JWlXXZUw/BjyWLnYzwW1oqnrruPLLpC6wG/X1E7R/d5Y6Jhw5S1U/kmX+vwP/PtGy6WD23A3DMIzTSkSuJLUH+yNVzT7Un/GmmT13wzAM47RS1SeA6jPdj5nM7LkbhmEYxgwjqWsepk99ff26t3L7Wub6b7Wu0+3H/7tRF9fkseNwPyKQSCrRhIWqkh9w0jccx24TLIXCHBdl+R46+yP0Dse5dFkpL+3uxOOy0TccP1bn8to8bDZh26E+yvM92EToHY5htwmhaOoWTp/LzuwyPwD7W4aIJSxyvQ4GRxIAeF12RmJJllQHaekJ0x+KU1vi40hnmNoSH7k+F4Ph2LFya2+EeFJxO2yIQJ7fSVtfBIDzFxbxyt5u/G4HqxcUMRCKsuvoALGExQcums3WA51EYkmSllKQ48Jus9HUFSJppeoB8LjsOB3C3uYhAJbWBOkeiNLeH8HjtFNZ6OVoV6ovhztCFOe6yQs42ds8RHm+hxyfk8bOMNGERdDnoD8UpzzfQ38ozkgseWx7RYSLl5aw+2g/kGr/cEcIv9vOstkFWKo4HTZ2HenDSr8vKgt97G0ewibgc9vJD7gYGkmQSFrUzy/mia2tFOa4yPE68LodNHWFGY4k8DjtOOxC0OekYyBKjtdO33CcwpzU+rGERVmeh4DXQUPbMJUFXpwO4WjXCHl+B33DcSoLvAD0h+L4PXZyfU4Gw3E6B6Lkeh3MKc/hUNsQkbiFx2kjHLNYWJVDY2fqjp+r6mfx2KamdN8dWKp0DkTxOO1E4kmqirw0d4+wsCoHj8tB71CUZFJxOoT3XTiP1q5BRn7/Wwo++HEKgz5irY20PPMMR//0BK6An5pv/jMtXUMssvXz7Bf/mosfeIiRjc9SvHIlf/zwJ6i+8mIW3fppjsRcLJ6VT/tAlPaeIQ5/7hNc9P/9I1q3hJJcN4+8eoSLls8i1+vgydcamV9dwKziHBRh6/52ytr3YsXjlJ9/PnaXmwfe/T4ArnvyMe5/dj+LawqwPXk/4YuvY55jGP+cuXT2htj7N7eT+OzfUvDYL+l998fxuOzUlzoI+/Lpuu9XPFO4hlnFXhraQly0pAS/18mTr7WSsJRE0gKgPN+D3+PgcEeIpJX6n3A5bMQSFgGPg0g89X4GSFpKnt/JcCRJImnhctiIJxVVxW4TkpbisNtQVSwFj9PGSCz1nnXYbSSSFoU5LsLR5LH1kpYe+3uNyix7XXY8LhuqEI4miSWsN3wO5QecRGLWsbb87tQB21A09X8K4PfYCUeT+Nx2eoZix+JEIBRN4nbYiMSTBDwOhiMJHHYbOV47n7vpvCndDmacPm+HPfdIfX39uvQAMrVnuC9TMhRJsuNw/7GyCMc+CPqG48feWADdg1F2Hj1+6mlfc/+YN36mbYf6js0ryfcQTx7/guZx2li9oOhYefybPfMhEbl+FwHv8TMzK+fk0zWQemMX5nq4YFHxmHUzYzM57DZC0cSxssthw+s6vm02EeLpD8sDrUPMKT9+p0nvcIxwxroAuxqPj6WQ+aHW0DZ87AN2b/MQ584tOFb2uGwsqc49FjuvKnhs+sIlpce+SAAkLevYBxjA4EiCF3d3AvDcjg5K8o4/J6lnMMr5C4+/npFYkr7h1LpPbG0l4Dn+mrR0h6ko9Ga0o7T0jgCpv/eS6uN9EhG6Bo/fMhuKJmjvS5WTSR3TX6/bTufA8diqIi+xRGq7L15eTl358et6dh4dYH56219v6Dr2+gA47TbyA8frjcQsltfmAdDeO8LQyPEvkYfu/m8Gf/uTY+WB4Qj+ytQDteZc916CX153vH/zFnPd7+87VvaUlFH9nz8HoE197Pv0x44tczntXHj/g8fKcbURT/+PxgcGmF9dcGzZ3p/8FzVlx1+zFw70gf346733ruMPHztw/x9YWJv6f913tJsntjQeW3bOV/8fteXH62lo7sXzntTt10PhBEtrji9bNbfgWGIHaOuLcLD9+O3Rme+f4UiCzJ0jj9POQPj4/3Lmey/z7wBgt8mY921mmyOxVDIdlfkeGF8eiSUZiVrH2svNeI9mvgcBcr2OY3//UDRBeX7G//lQbEx/UzsLqW2xyRv7kNlf4+x01p9zzxwtzowQZxiGYRgn93bYczcMwzAMYwpMcjcMwzCMGcYkd8MwDMOYYUxyNwzDMKaFiJSJyG/Tj5/dLSKPpAdNyRa/UkSunUS960TkK6e2t29sX0SuE5GvTWH9BelH147+DIrIl051PyfjrL+gzjAMw3j7kdStBQ8AP1fVD6fnrQRKgf1ZVlsJ1JMa2326+uVQ1USWxWPaV9UNpIZnnRRV3ZeuAxGxkxqc5oG31OE3yey5G4ZhGNPhciCuqsdGP1PVbar6vIj8UkSOPTNeRH4lItcB3wZuGR24RUQKRORBEdkuIq+IyPLxjYjIZ0XkjyLiFZE6EXlURLaIyPMisjAdc7eI/EBEnga+JyLnichLIrI1/XuBiLgmaP9WEbkzo4470vGHROREw9kCvAs4qKpH3+oL+WaYPXfDMAxjOiwFtmRZ9hPgduAhEQkCa0mN+PYNoF5VvwAgIj8CtqrqDSJyBfAL0nvG6eVfIDVG+g2qGhWR9cDnVfWAiKwBfgxckQ6fD1ypqkkRyQUuUdVE+jn331XVm0Tk/2/vzoPjqu5Ej39/6pZaS2uXtcuSZUu25T2RbewYbNYYwkDMYpxMFioBJ/OKyZCQV4+pmpl4JvNeMjAzjxkCD1wQIGQGhj1MbCAYG2PABrzbsrxKtqzNUmvtRb2f90e3QttIdkMsGTe/TxXFvef+7jmn/Yd+99zlnDPbv/2MfpcAS4BpREb0L5zl968iZlGc8XaxJXfvxfKt+8U0k55SSn0WIrIaWB1TtNYYs/Zc5xljNovIQ9E1zm8CXowm2jNDlwA3R8/ZKCL50YsBgG8DrUQSe0BE7EQuEp6PqccWU9fzxpjh2XiygadEpAYwQDLxecUYEwYOiEjRaEHRuwA3AH8dZ73n3UWV3GMntFFKKXVhRRP5aMm8ATjbreuniay/vgr43igxI01vOzy1334io/hyoJnIY+Z+Y8zcEc4BcMds/xzYZIxZISJVwNtn6WcsX8z22abevRbYGV2z/oLQZ+5KKaXGwkbAJiJ3DheIyHwRWRrdfRK4G8AY0xAtcwKZMXW8Q+QCABFZBjiMMcPzU+8CfgC8KiKl0fJmEbk1Gi8iMmeUvmUTedkN4PaY8jPb/6y+wQW8JQ+a3JVSSo0BE5l4fwVwdfRTuAZgDdAePX4KaASeiDltE1A3/EJbNL5eRPYCvyTyXD62jXeBnwLrRKSAyIXA90VkD5E7BzcysvuAX4jIe0DsJPxntv+piUg6cDXw0mc5/3y5qG7LK6WUungYY9qBlSMdiybBGmJGuMaYXmD+GaGfSNDGmDUx228Ab0R3HcDyEeJvP2N/K5EX7Ib97Vnaf3KUOuyMwBjjAfJHOjaedOSulFJqXEXfUD8IPGiMGThXvPr0dOSulFJqXBljNgATL3Q/EpmO3JVSSqkEo8ldKaWUSjCa3MdQnj0yL8LGRhcfHHMz6Aly62XVZNis3HpZNauqw8ybnEdRto0Ftfm81ehl0BPkld0e9h53csm0Au5/y8W+k0N0DQYREXYc7WPIH6bfHeTRzf3MnlxE92CA7kE/6SkWAiFDRVE2ydYk/MEws6tyeOwDH//6tpumLj/dgwFmVObw6l4vz7zTxvpd/UwpsbP8ksl4fEFae31UFGWSkZZMw/E+Fs8sw+0NMeAJ4vYG2X9yiFDY0NDmo8Xh50jbIDWldpbXlzMp6MAfCPPAZjfOoRCvbj3B4pll5NhTcHtDdPZ5SRLBNRQgFDY0nfLwjStrOdjmJs1mZXeLlwFPCEuSYE+zMr0iiwybla8urGZKSQbXzi9ndlUOHl+Q714zjazXf8sVKZ3MrSlk4oQMDrc5mVRs556bZzOzuojmbi997hCNJ/rodQZYODWfL1VkcqLbyzfLXNRV5XPlvFK2H/eSYbPypakl9LlD3PNCN9uOuLh0VinXTRTaHG6unFdKhs3KB0cGKcy2UZBl46+WT2bJzGKSLUkcaXdTU5bFYrsLjy9M16CfGZU5OJxBBj1B5k3OJT87lXcPOQkbyLBZ6HMHsVqEmZXZiIDVInzUPMShjiE8vhCv7eln9uR81u0eYEqJncNtLiYW2nnxwz6WzSmh4uA7rPj3o8yaPIG9LR6mlGSwaPoEGk70U5qfxuNvOzjc4WXWpDz+8Y0BVl05lbcanNjTrNxzyxxWzon8Wzo9QexpVvLsNjLTrVw9fxJp199G/+Ur2VF1GZv3nOR37x1j/YfN/NvgNN60zWDzvlM0tfdTU5GP0+PjV2800dHjpLd6Huu2HuXEKSfeP7udo219LHrhZVocbu5au5e3drUTCoX59RErrV2DPL5uPxv29rLrcCf/8loz+451I3u28sGBdg7ULOOx9Qe570A6qfVLWFBspaXbyZTHfov9nx5lc+ml3LS0loy0ZIb+5/288PZhsiZW8sSGZhb178H2o5+xYa+D9TtOcqS1j7CBXScIZ/cAABO3SURBVI4wP3u+hbZuJzMqs7ksuZOKwkzkpV+TtGU9kz2tTCrKYHKxncZ2H/fcOpfrFpRTkpeKyxvmL1fMoXswQFqKhdRkC/Nr8/nhDbOoKbVTkpdKjzNAe5+fn9xYR0Obj353kPzMFHpdQaaVZ1Jfk0eLw095fip/ddMcHM4gMyZmU5BlIz8zBa8/TEGWjcY2D7OqcrCnWhnwhHA4gwx4QnylrpBp5Zksm11Mi8PPyZ4A3kCIDQ1OvrawgqyMZG5YVIklScjPSmFBbT5LZ5fxvetmUFOWhdsb4palU3A4g8yqzudbV0+lMMfGsVN+SvJS+dY1dcytzqV7MMCkogySBKaVZ1FVmM6HzV6KclJpcfgpzknllqWjrv+iPgc0uY+zLXvb4459+I3WMezJn24oYM4dFNXQOhR3bF1FVtyx71MWd+z58tL2/rhj19xYFXfsV+dXxB2780h33LHny8K60rhj//P9nrhj77+zPv5O7H4v/tjzZMDlO3dQ1DeXVcUd++DLe+KOzUnX16PUp6PJXSmllEowmtyVUkqpBKPJXSmllEowmtyVUkqNCREpFpFno9PPHhCR9SIy6pt4IjJXRK6Lo941IvLT89vbT7YvIjeIyL2fso7jIrIvOoXt9vPdx3jpWxpKKaXOO4msu/oy8JQxZlW0bC5QBBwe5bS5QD2wfgz7ZTXGBONp3xjzKpF12z+ty40xjs/YxfNCR+5KKaXGwuVAwBjzyHCBMWa3MWaLiDwtIn+cM15E/kNEbgD+AbhteOEWEckTkVdEZK+IbBOR2Wc2IiJ3ishrIpImIpNF5HUR2SEiW0RkWjTmSRH5VxHZBPyTiCwQkfdFZFf0/1Oja7Cf2f7tIvKrmDr+PRrfJCJnW872gtORu1JKqbEwE9gxyrHHgB8DvxORbGAxkRXf/g6oN8bcBSAiDwK7jDFfF5ErgN8QGV0TPX4XcA3wdWOMT0TWAj80xhwRkYXAw8AV0fBa4CpjTEhEsoDLjDHB6Dz3/8cYc7OInNn+7Wf0uwRYAkwjMqJ/YYTfZoA/iIgBHo2ueT/uxjW519fX3wukjmebF8r27dvXXOg+KKXUWBKR1cDqmKK18SQzY8xmEXlIRAqBm4AXo4n2zNAlwM3RczaKSH70YgDg20ArkcQeEBE7kYuE52PqscXU9bwxJhTdzgaeEpEaIsk4Oc6f/IoxJgwcEJGiUWK+Yoxpj/62N0XkoDHmnTjrP2/Ge+SeqklPKaUSQzSRj5bMG4Cz3bp+msj666uA740S84lsTyQZA+wnMoovB5qJPGbuN8bMHeEcAHfM9s+BTcaYFSJSBbx9ln7Gip3RaKS+DS9zizGmS0ReBhYA457c9Zm7UkqpsbARsInIncMFIjJfRJZGd58E7gYwxjREy5xAZkwd7xC5AEBElgEOY8xg9Ngu4AfAqyJSGi1vFpFbo/EiInNG6Vs20Bbdvj2m/Mz2PxURyRCRzOFtIo8M9n/W+v4UmtyVUkqdd8YYA6wAro5+CtcArAGGR7angEbgiZjTNgF1wy+0RePrRWQv8Esiz+Vj23gX+CmwTkQKiFwIfF9E9hC5c3AjI7sP+IWIvAdYztL+p1UEvBtt/0NgnTHm9c9Qz59MX6hTSik1JqK3qFeOdExE0oEa4JmY+F5g/hmhn0jQxpg1MdtvAG9Edx3A8hHibz9jfyuRF+yG/e1Z2n9ylDrsI7TTBIx2t2Bc6chdKaXUuIq+oX4QeNAYM3Ch+5OIdOSulFJqXBljNgATL3Q/EpmO3JVSSqkEo8ldKaWUSjCa3MdQarIFty/ID64qYWqxjcYOP8sXVrFg2gSOtvWxL5RLdVkOc6cU8PK2TuZVWPnu8un8859XkZlq4Sd/VsVPr7CzqMbOFbPycQ6FmF+bx/yaXOoq7KyYm8bjrx0kPSWJI6cCbD3qxusPs3ZdI/uOD7DtUB+BUJgFpfDDxWnYrEJ7f4j3DjioyhX+ZdUUJmQm8dKHPdz7+A7qKvOZNymTDxq7SBJhxaVTcA35cfnCZKVZONrpJc9uYUeTi5JsC//vjln4AmHeP9hHTmYqr55MYsNeB9fWJFGen0p7n5+n3zxMjt2GzSpMn5jD4Y4hKgozOTXg51Cnn/ue20flhFQWz6qgtiiF5V8uZsuBHg62urhuQRUef4j3950kGDLYkq1kpqeQkWolLTWZY4tu4VCbi9KCTOomFdDS4yMQDPOfGw+TJiG8fsNVcybwpdoJhA3MmlyEb2CQiQWp2GbV8+aOVmZV5fM3t0yhrCCNky8+S3VRKncvzeB/+DdTlGcns2Yq82oK6RkY4minm3/4zlye+vF8XENBLDYbG3e3s2JGFhUFaWxpcLBlIAOXL8y9q+bx++2dWFNTOdoVIEmEV7a28Te3TaPXFaC+toAkEcoK7CyaWc7mRhcP/f4wJVkWUpOTGPAEWfujBby5s50J9iQKslM51uXDahG+t6wItzfECzKTBaXw6LrDLJiSSY7dxoeHHBgDTk+Qwkxh5ZJSMtKSue+mAgCuqMtkTnUe+5od+Pr7qSrOYs8JJ4tmltE94KUgO43uPjf+QIjUFAvfWFTO+l39AHj9QaoLU8jPSmHl0sk0tgwSCoVx9LsZ8oVwDwWYUV2IyxvkujmZHGodIN1m5aODHby7r51ffGsq/mCYdoeT8vxUbMkWAqEwj/73ft474GDV5dV09Xs5UVRHeWEWWw8PcN2CMqaXpeHx+ultbKT06mtxDwW4anYJV8wrZfOuExw80cusSfn0OP10DPpZMi2X/9tdzeKphdisQorVwsxJBUwqyWbxzDKWz0wj+PD/5rn3OmjMqKIoz46z5STe7m7aNm3keJcHi0VIsQjPvnWIfpePrn4fTd0B7n7kQ/7+u/VMLrGz76SHpg4nL759mIzUZLIzUlj/YStF2cls2NPO7IpUirJTmFiYwanBEEfaXXT2eki3JTGlLIc7Fk1gyfRcjp9yY7MmMel3D+INhDjeNUS6LYlN+xyEwgZf0GC1wF031uH71c/pHvCx/qN2vnPFRKwWEBG+Xp/L439oYt8JF/+9rQWPL8zOJidef4h2h5Nn3zrI+409tPb6ee7tIwAcaxtg6/42jrR//Pm3NRTgN5s7cfsMx7s8eANhnnmvi/ZeL7PKUgiFw1ROsHG4w81fPrJz/P+oqrhpch9nwQO74o595p2TY9iTP91vt5yIO7Z/yJw7KOrmr1TEHXvLpZVxx54v6bb4X1XZ1+4+d1BUbXFa3LGfnMhr7D3wSmPcsZsbXXHHFufazh0U9fi68f9kOOlT/GNf8+XSuGNfbRw8d1BUfmZK3LFKgSZ3pZRSKuFocldKKaUSjCZ3pZRSKsGMx3fu3vr6+jXR7apxaE8ppZT6Qhvz5L59+/ZfDm/HJHmllFJfACJSDDxAZFpXH3CcyIIxe4jMUpdKZMGWh4wxT0XP+XPgf0WrcAF/YYzZEz22HPg3InPCP2aM+WOOUR/TGeqUUkqNCYksrP4y8JQxZlW0bC6RBVaOGWPmRcuqgZdEJMkY8wSRJVyXGmP6RORaIsvKLhQRC/AQcDWRtdw/EpFXjTEHxv3Hfc7pM3ellFJj5XIgYIx5ZLjAGLMbOO073+iCKz8BfhTdf98Y0xc9vI3Imu0QWRv9qDGmyRjjB55l9JXfvtAu+Mi9vr7+XiK3ZRLK9u3b11zoPiil1FgSkdXA6piitcaYtTH7M4EdcVa3E5g2Qvn3gdei22WcfmHQCiyMs/4vlAue3IFUTYRKKXXxiSbytecMjM8nZgsSkcuJJPclo8UA8c+Q9QWit+WVUkqNlQbgy3HGzgP+OA2iiMwGHgNuNMb0RItbgdgpLMuB9vPQz4SjyV0ppdRY2QjYROTO4QIRmQ+cNm+0iFQB/ww8GN2fCLwEfNsYczgm9COgRkQmiUgKsAp4dSx/wMXq83BbXimlVAIyxhgRWQE8ICL3Al4+/hRusojs4uNP4R6MvikP8HdAPvBw5IV7gsaYemNMUETuAt4g8incr40xDeP6oy4SmtyVUkqNGWNMO7ByhEOjrpRkjLkDuGOUY+uB9eend4lLb8srpZRSCUaTu1JKKZVgNLkrpZRSCcayZs2acWts7dq1y1avXv32ucoSRcOxzjW+QJjD7S7uuG4a2bYQ3ZZMFs8qJ0mEPcd6aDnl5FiHk/opOUyvyGLdBy3UlGXz1t4e/vq2mWw90M2tS6dQPiELp9vDtkN9ZKVZOdDqQoC2vgApVkGAqSU2ZlVmk55iwRcIIyIYE8Y1FMSSBJmpFnIzkhj0hslKS6J6Yh4fHOqltjiFHleIuooMwmFDZVEmmekp7GvqprPHQ3P3EIEQ5GZY6HEFsSQJi6fmcrTDTfdggNx0K+09Tj465uayuhx8/hCleWmc6B6itiSDfpcP51CA4uxkmk556erzIALfXFpOOBTCHzD4/H5mVOaw7qMOOno9NHf5mVZu5/WdXQy4/SyuK6LN4eQ/3mmntjQdR7+HK6rSyMjLY9POE7Q7XHT1+5g9KYfmThfdg0NcMjUPx4CXubXFlOTZeHN7C6nZmdSU5bBpVys5Gckcv+eHtE9bzIpLa+natpVg5VRuWlpLyVeW0NTWR574Odrt4aTDTXFOKiX5djJTkmg82cdQMExWejKV1WV0OAbpcwVoOuWhNDeFN3a0kWwRXB4fbX0BMCFsycK2gz14/GHm1+Tj9HgpyE6l8biDFEsYhzPyb+vxhbFahI7uATy+EF+rL+Hlbe0IsGTGBN5v7MbtDdHrCVJgt2BLFpxDITr7hrBZLRRkp1BVlEmbw417KMD0ibn8ZlMrB5q7cfvCFOXY2Nvcy9QZ1ZzqdXOk3YmFIFsPD3LD4kqe3XSMG+aXc6LbRY8nRO+gh6x0Kysvm0RdZS75WalkZdjw+vzMqSkiNzOdP+xo4/pLKnlv70nKJ9ipK08nxZqEPxjG5w+BwJwphbR3O1k22c6TWzox4SDHTnkZcAfpc4fY3dSLYGhocbKorogdRxyc6HTiCxiumV9JqKCEr9TmktrZjMWWQq9f8HiD+AIh6ktTycjKYMPOVtJtFgY9AXIzk3lzby/pyeD0+Nh+2EFWupXXdnaTvuhS5j7/j8xeUMfJF5+j99ARXO0deDq7OVC+kOIcG9lpSdRV5dLS5SbZmoTLG2B6aTrPbTnBpMJ0irOT6XcHyLEns7d5gPzMZHYe68MYaHUMsXBqPjuaBiAcJjvdQrrNQmefj9I8G/uPDzC9Ips9zf24vCHmTcmj8vrryU6z0tzpIifDSm1pBs1dQ/R7wlw1t5BBt4+GgpnsbHZz4yUlVA2eYGt3EhhocfhIS0ni+1+dQlP7ILl2K1+bX07DiT4ml2VzzewSjnYOkJFi4ZLpE6ibmMVbe7pwe4NMLk5nQW0uyy+p5pmNR8i3W+h3B5lekcG1C6rY2thN2Bh8QUNnv4/sdCurr69jZX0+1gz731/gP7NqFOM9cvfW19evif0PXSnujzbs6Tpt/+k3D48SCXUV9rjr9QXjn+OhOD8j7tiZE7PijvX6R+/DZbOLT9tPT7GMGluWm3zafnN73yiR0NrtirN38NHB+D+V/a93jsYdm5Q00pwbEftbPKftV04YfaLGAU/4tP2uweCose82DsTZO8jKSIk7dvOetrhjf7/91KjHTgVPb7M0d/Q+9LhP/50drtCosc+/e3LUY2eaP7Uw7tj2Xm/csZfOLIo7dvaU+GM7ejyjHps5Mf6/BfW1+XHH3nb5SJPFqYvFuL4tH7tC3DBdKU4ppZQ6v/SZu1JKKZVgNLkrpZRSCUaTu1JKKZVgNLkrpZQaMyJSLCLPisgxETkgIutFpFZEhkRkl4g0isiHIvLdM85bJiK7RaRBRDbHlP84WrZfRJ4RkYRbMvx80OlnlVJKjQmJTAz/MvCUMWZVtGwuUAQcM8bMi5ZVAy+JSJIx5gkRyQEeBpYbY1pEpDAaVwb8CKgzxgyJyHNEFo95crx/2+edjtyVUkqNlcuBgDHmkeECY8xu4LTvFo0xTcBPiCRugG8CLxljWqLHY78TtgJpImIF0tElX0ekyV0ppdRYmQnsiDN2JzD8cX0tkCsib4vIDhH5DoAxpo3I0rAtQAcwYIz5w3nuc0L4PNyW9ybit+7bt29fc6H7oJRSY0lEVgOrY4rWGmPWftbqYratwJeBK4msHrdVRLYB3cCNwCSgH3heRL5ljPntZ2wzYV3w5D7SxDZKKaU+/6KJ/GzJvAG4Jc7q5gGN0e1WwGGMcQNuEXkHmBM91myM6QYQkZeAxYAm9zPobXmllFJjZSNgE5E7hwtEZD5QGRskIlVEbrc/GC36HXCpiFhFJB1YSCTxtwCXiEh69GW9K/n4gkDFuOAjd6WUUonJGGNEZAXwgIjcC3iB48DdwGQR2QWkAk7gQWPME9HzGkXkdWAvEAYeM8bsBxCRF4g8nw8Cuzj7nYMvLE3uSimlxowxph1YOcKhtHOcdz9w/wjlPwN+dn56l7j0trxSSimVYDS5K6WUUglGk7tSSimVYDS5f459++raC92Fs9rfMhh3bGqKnDsoyuMPxR07qTQ37tjz5bbLpsQdGw6buGNPdHvjji3MGv/XZZbOKYs79vr6orhj2/v8cceW2C1xx54vpXnxT12+Zf+puGP3Ho0/tiQ/Pe5YpQDEmPj/+CillFLq809H7koppVSC0eSulFJKJRhN7koppVSC0eSulFJKJRhN7koppVSC0eSulFJKJZj/D1gRFG2801akAAAAAElFTkSuQmCC\n",
"text/plain": [
"
"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ast.type_clustermap(plot_name=\"./img/celltype_protein_cluster.png\", threshold = 0.7, figsize=(7, 5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note: `threshold` is the probability threshold above which a cell is assigned to a cell type, default to 0.7. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 8. Hierarchical model specification "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the marker yaml file, the user can also add a section called `hierarchy`, which specifies the hierarchical structure of cell types. Here's an example:\n",
"```\n",
"hierarchy:\n",
" immune:\n",
" - B cells\n",
" - T cells\n",
" - macrophage\n",
" epithelial:\n",
" - epithelial(basal)\n",
" - epithelial(luminal)\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some notes: \n",
"1. The section would be accessed by key `hierarchy`.\n",
"2. In the section, the higher-levelled cell type names should be the keys.\n",
"3. The values in the section should also exist as the cell type names in the `cell_types` section. (e.g. if we have `\"B cells\"` in `marker[\"hierarchy\"][\"immune\"]`, we should also be able to get `marker[\"cell_types\"][\"B cells\"]`)\n",
"\n",
"This section could be used to summarize the cell types assignment at a higher hierarchical level. (e.g. a cell is predicted as \"immune\" instead of \"B cells\" or \"T cells\")"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
immune
\n",
"
epithelial
\n",
"
\n",
" \n",
" \n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_1
\n",
"
0.435195
\n",
"
0.005911
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_2
\n",
"
0.385014
\n",
"
0.044718
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_3
\n",
"
0.188510
\n",
"
0.330741
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_4
\n",
"
0.444039
\n",
"
0.011956
\n",
"
\n",
"
\n",
"
BaselTMA_SP43_115_X4Y8_5
\n",
"
0.468241
\n",
"
0.017085
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" immune epithelial\n",
"BaselTMA_SP43_115_X4Y8_1 0.435195 0.005911\n",
"BaselTMA_SP43_115_X4Y8_2 0.385014 0.044718\n",
"BaselTMA_SP43_115_X4Y8_3 0.188510 0.330741\n",
"BaselTMA_SP43_115_X4Y8_4 0.444039 0.011956\n",
"BaselTMA_SP43_115_X4Y8_5 0.468241 0.017085"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"hierarchy_probs = ast.assign_celltype_hierarchy()\n",
"hierarchy_probs.head()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To make it more clear, here's a heatmap for the cell assignment in a higher hierarchy:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The way it is calculated is simply summing up the probabilities of the cell type assignments under the same hierarchy."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}