{ "cells": [ { "cell_type": "markdown", "id": "c000", "metadata": {}, "source": [ "# What the well median hides\n", "\n", "A well median summarizes a population that is rarely uniform. If a fifth of the cells respond strongly and the rest\n", "do not, the median reports a mild effect, and the fact that there are two populations is lost before analysis\n", "starts.\n", "\n", "This page works on the 13,578 cells of {func}`~mantispy.ds.jump_cells`, 24 wells of the JUMP plate `BR00121438`,\n", "to ask what the median misses." ] }, { "cell_type": "code", "execution_count": 1, "id": "c001", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:02.643313Z", "iopub.status.busy": "2026-10-21T01:15:02.642919Z", "iopub.status.idle": "2026-10-21T01:15:12.622261Z", "shell.execute_reply": "2026-10-21T01:15:12.620042Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":10: UserWarning: 285 of 5839 features have no spread in at least one group of 'Metadata_Plate' among the rows selected by reference='negcon'. Adding epsilon=1e-18 to their scale multiplies them by up to 1e18, so they dominate every distance downstream. They are flagged in var['degenerate_scale'], and pp.feature_select drops them by default; to drop them now, adata = adata[:, ~adata.var['degenerate_scale']].copy(). A variance or outlier cut does not catch a feature that varies across a plate but is constant among its control wells; the flag does.\n", " mt.pp.normalize(cells, by=\"Metadata_Plate\", reference=\"negcon\")\n" ] } ], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import scanpy as sc\n", "\n", "import mantispy as mt\n", "\n", "raw = mt.ds.jump_cells()\n", "cells = raw.copy()\n", "mt.pp.normalize(cells, by=\"Metadata_Plate\", reference=\"negcon\")" ] }, { "cell_type": "markdown", "id": "c002", "metadata": {}, "source": [ "At this resolution a feature that barely varies among the control cells can come back enormous from normalization\n", "without being flagged as degenerate, so `drop_outliers` joins the default selection." ] }, { "cell_type": "code", "execution_count": 2, "id": "c003", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:12.626125Z", "iopub.status.busy": "2026-10-21T01:15:12.625917Z", "iopub.status.idle": "2026-10-21T01:15:20.971132Z", "shell.execute_reply": "2026-10-21T01:15:20.970433Z" } }, "outputs": [ { "data": { "text/plain": [ "AnnData object with n_obs × n_vars = 13578 × 1564\n", " obs: 'Metadata_ImageNumber', 'Metadata_ObjectNumber', 'Metadata_CellsObjectNumber', 'Metadata_NucleiObjectNumber', 'Metadata_AbsPositionZ', 'Metadata_AbsTime', 'Metadata_BinningX', 'Metadata_BinningY', 'Metadata_ChannelID', 'Metadata_ChannelName', 'Metadata_Col', 'Metadata_ExposureTime', 'Metadata_FieldID', 'Metadata_ImageResolutionX', 'Metadata_ImageResolutionY', 'Metadata_ImageSizeX', 'Metadata_ImageSizeY', 'Metadata_MainEmissionWavelength', 'Metadata_MainExcitationWavelength', 'Metadata_MaxIntensity', 'Metadata_ObjectiveMagnification', 'Metadata_ObjectiveNA', 'Metadata_PlaneID', 'Metadata_Plate', 'Metadata_PositionX', 'Metadata_PositionY', 'Metadata_PositionZ', 'Metadata_Row', 'Metadata_Site', 'Metadata_Well', 'Metadata_Center_X', 'Metadata_Center_Y', 'Metadata_ObjectType', 'Metadata_ImageID', 'Metadata_Source', 'Metadata_JCP2022', 'Metadata_InChIKey', 'Metadata_Perturbation', 'Metadata_Perturbation_Type', 'Metadata_Control_Type', 'Metadata_Control'\n", " var: 'object', 'feature_group', 'feature', 'channel', 'scale', 'angle', 'gray_levels', 'radial_bin', 'params', 'is_feature', 'feature_kind', 'selected', 'degenerate_scale'\n", " uns: 'mantispy', 'pca'\n", " obsm: 'X_pca'\n", " varm: 'PCs'\n", " layers: None (.X)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.pp.feature_select(\n", " cells,\n", " operations=(\n", " \"drop_degenerate\",\n", " \"variance_threshold\",\n", " \"correlation_threshold\",\n", " \"drop_na_columns\",\n", " \"blocklist\",\n", " \"drop_outliers\",\n", " ),\n", " na_cutoff=0.0,\n", ")\n", "cells = mt.pp.subset_features(cells)\n", "sc.pp.pca(cells, n_comps=20)\n", "cells" ] }, { "cell_type": "markdown", "id": "c004", "metadata": {}, "source": [ "## The typical cell, or some cells?\n", "\n", "`method=\"mahalanobis\"` asks whether the group's typical cell is far from the controls. `method=\"ks\"` compares the\n", "whole distribution of cell distances against the controls', so it answers a different question: whether *some*\n", "cells moved, even when the middle of the distribution did not." ] }, { "cell_type": "code", "execution_count": 3, "id": "c005", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:20.972853Z", "iopub.status.busy": "2026-10-21T01:15:20.972698Z", "iopub.status.idle": "2026-10-21T01:15:22.215372Z", "shell.execute_reply": "2026-10-21T01:15:22.213822Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ ":2: UserWarning: the well-block null draws from 4 block(s) in the null half, too few to calibrate finely, so the test there is conservative rather than anti-conservative. For an exact well-level test, aggregate to wells with mt.tl.aggregate before testing.\n", " mt.tl.hit_calling(\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ ":2: UserWarning: the well-block null draws from 4 block(s) in the null half, too few to calibrate finely, so the test there is conservative rather than anti-conservative. For an exact well-level test, aggregate to wells with mt.tl.aggregate before testing.\n", " mt.tl.hit_calling(\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
median distanceKS statistic
group
JCP2022_0986886.6951460.534674
JCP2022_0104044.7166320.457240
JCP2022_1159634.7635410.375978
JCP2022_1165604.5274590.345986
JCP2022_0674324.4932210.333079
JCP2022_0614213.8199950.228362
JCP2022_0797153.8857440.206773
JCP2022_0010363.7092990.206682
\n", "
" ], "text/plain": [ " median distance KS statistic\n", "group \n", "JCP2022_098688 6.695146 0.534674\n", "JCP2022_010404 4.716632 0.457240\n", "JCP2022_115963 4.763541 0.375978\n", "JCP2022_116560 4.527459 0.345986\n", "JCP2022_067432 4.493221 0.333079\n", "JCP2022_061421 3.819995 0.228362\n", "JCP2022_079715 3.885744 0.206773\n", "JCP2022_001036 3.709299 0.206682" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "for method in (\"mahalanobis\", \"ks\"):\n", " mt.tl.hit_calling(\n", " cells,\n", " groupby=\"Metadata_Perturbation\",\n", " reference=\"negcon\",\n", " use_rep=\"X_pca\",\n", " method=method,\n", " n_permutations=500,\n", " key_added=method,\n", " )\n", "\n", "comparison = pd.DataFrame(\n", " {\n", " \"median distance\": cells.uns[\"mantispy\"][\"mahalanobis\"].set_index(\"group\")[\"distance\"],\n", " \"KS statistic\": cells.uns[\"mantispy\"][\"ks\"].set_index(\"group\")[\"distance\"],\n", " }\n", ")\n", "comparison.sort_values(\"KS statistic\", ascending=False).head(8)" ] }, { "cell_type": "markdown", "id": "c006", "metadata": {}, "source": [ "One caution applies to every cell-level test on this page. Cells from one well share its plate position, its\n", "seeding density and its focus, so they are not independent draws; the unit the experiment randomized is the well.\n", "A null that shuffles cells assumes more evidence than the design supports, and p-values come out small for reasons\n", "that have nothing to do with the treatment. Read a cell-level call as a ranking, and put the well-level call beside\n", "it before believing it.\n", "\n", "## What the median misses\n", "\n", "Each cell's own distance from the control centroid is on `obs`. For each compound, how far its typical cell moved,\n", "and what share of its cells sit beyond the distance that only one control cell in twenty reaches:" ] }, { "cell_type": "code", "execution_count": 4, "id": "c007", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:22.218356Z", "iopub.status.busy": "2026-10-21T01:15:22.217955Z", "iopub.status.idle": "2026-10-21T01:15:22.236784Z", "shell.execute_reply": "2026-10-21T01:15:22.235719Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
median distanceshare beyond the controls' 95th percentile
Metadata_Perturbation
JCP2022_0986886.700.36
JCP2022_1159634.760.17
JCP2022_0104044.720.16
JCP2022_1165604.530.15
JCP2022_0674324.490.15
JCP2022_0350953.810.12
JCP2022_0797153.890.12
JCP2022_0018903.640.10
JCP2022_0614213.820.07
JCP2022_0507973.200.07
JCP2022_0010363.710.07
JCP2022_0339243.070.05
JCP2022_0258483.530.03
\n", "
" ], "text/plain": [ " median distance \\\n", "Metadata_Perturbation \n", "JCP2022_098688 6.70 \n", "JCP2022_115963 4.76 \n", "JCP2022_010404 4.72 \n", "JCP2022_116560 4.53 \n", "JCP2022_067432 4.49 \n", "JCP2022_035095 3.81 \n", "JCP2022_079715 3.89 \n", "JCP2022_001890 3.64 \n", "JCP2022_061421 3.82 \n", "JCP2022_050797 3.20 \n", "JCP2022_001036 3.71 \n", "JCP2022_033924 3.07 \n", "JCP2022_025848 3.53 \n", "\n", " share beyond the controls' 95th percentile \n", "Metadata_Perturbation \n", "JCP2022_098688 0.36 \n", "JCP2022_115963 0.17 \n", "JCP2022_010404 0.16 \n", "JCP2022_116560 0.15 \n", "JCP2022_067432 0.15 \n", "JCP2022_035095 0.12 \n", "JCP2022_079715 0.12 \n", "JCP2022_001890 0.10 \n", "JCP2022_061421 0.07 \n", "JCP2022_050797 0.07 \n", "JCP2022_001036 0.07 \n", "JCP2022_033924 0.05 \n", "JCP2022_025848 0.03 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "distance = cells.obs[\"mahalanobis_row_distance\"].astype(float)\n", "cutoff = distance[cells.obs[\"Metadata_Control\"].to_numpy()].quantile(0.95)\n", "group = cells.obs[\"Metadata_Perturbation\"]\n", "spread = pd.DataFrame(\n", " {\n", " \"median distance\": distance.groupby(group, observed=True).median(),\n", " \"share beyond the controls' 95th percentile\": (distance > cutoff).groupby(group, observed=True).mean(),\n", " }\n", ")\n", "spread.sort_values(\"share beyond the controls' 95th percentile\", ascending=False).round(2)" ] }, { "cell_type": "markdown", "id": "c008", "metadata": {}, "source": [ "The DMSO wells, JCP2022_033924, sit at 5% by construction. The compound whose typical cell moved furthest, JCP2022_098688, has 36% of\n", "its cells beyond the controls' 95th percentile, and the next four have 15 to 17%. Even the strongest response on\n", "this plate leaves most of a well's cells where the controls are, and a well median summarizes a minority that\n", "moved together with a majority that did not." ] }, { "cell_type": "markdown", "id": "c009", "metadata": {}, "source": [ "## A change in shape\n", "\n", "The same idea applies feature by feature. {func}`~mantispy.tl.effect_size` compares locations;\n", "{func}`~mantispy.tl.wasserstein_features` compares whole distributions, so it sees a treatment that widens a\n", "feature without moving its centre." ] }, { "cell_type": "code", "execution_count": 5, "id": "c010", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:22.238994Z", "iopub.status.busy": "2026-10-21T01:15:22.238856Z", "iopub.status.idle": "2026-10-21T01:15:27.323439Z", "shell.execute_reply": "2026-10-21T01:15:27.322650Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAGuCAYAAAADLgZkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYvhJREFUeJzt3Xl4k1XaBvA7SZM0aZuutBRaWiibICCVfS+LggrqiArMjIiIy4johzjCiDPiDCPO6OACDDAguDMsVpBViwUKsi9alKWUvS3d2yRN0zTJ+f4ozRDaQvo26Zb7d129aN6cnDx5G/LknPcsMiGEABERkZeSN3QAREREDYmJkIiIvBoTIREReTUmQiIi8mpMhERE5NWYCImIyKsxERIRkVdjIiQiIq/GREhERF7Np6EDIKqr1NRU/PLLL7j//vsREBBQ5f6ysjIcP34c165dQ0BAAGJiYhAXFweZTFaljkpKpRKRkZG4++67oVarneqz2+04ffo00tPTERYWhq5du0Kn09UYn8lkwv79+2E0GnHnnXciLi6u2nK1qbe2MbjC1TgBICMjA8ePH4dMJkN8fDwiIyNrLJuVlYVTp06hpKQE0dHR6NGjh9O5v1lOTg6OHTsGHx8f9O7dG4GBgXWut7YxkJcRRE3c66+/LgCIU6dOOR0vLS0Vr776qvDz8xMdOnQQ48aNE6NGjRI6nU506NBBfPzxx1XquO+++8Tjjz8uHnjgAREeHi5CQ0PF2rVrHeX+9a9/idjYWNGiRQvxwAMPiB49eoiAgAAxf/78amP75ptvRHBwsOjZs6d44IEHhFarFZMmTRJlZWVO5WpTb21jcIWrcVosFvH0008LpVIpBg8eLIYMGSJUKpX4v//7P2G3253KmkwmMXnyZKFUKsXAgQPFAw88IMLCwkTnzp3F3r17q8RQUlIinnzySeHn5ydGjBghxo0bJ2JjY8XixYsl11vbGMg7MRFSk1ddIjSbzWLgwIEiODhYbN261al8SUmJmDNnjhg8ePAt6yguLhZ33323UCqVIj09XQghRMeOHcXcuXOF2Wx2lPvoo48EAPHRRx85Pc/p06eFSqUSM2bMcBw7efKk0Gq1YubMmU5la1Nvbcq6ojZxvvrqqwKA2LJli+PYjh07hEwmE++8845T2blz5woAYsOGDY5jxcXFonPnziI8PNwpydrtdjFq1CjRrVs3kZGR4ThuNBrFunXrJNdbm7LkvZgIqcmrLon95S9/EQDE5s2ba3zcgQMHblmHEEKsXLlSABBLliwRQgiRkpJSbV3h4eGiX79+TseeffZZoVKpRGFhodPxZ555Rvj6+gq9Xu84Vpt6a1PWFa7GabVaRWBgoBg2bFiVOu69914REhIiLBaL49jo0aOFr69vlbJ//OMfBQBx7tw5x7EvvvhCABBHjhy5bby1qbc2Zcl7cbAMNTtCCPz73/9G586dcf/999dYrm/fvrety9/fH7h+/QwABg0aVKWM3W6H3W53lKmUlJSEHj16ICgoyOl4QkICzGYz9uzZ4zhWm3prU9YVrsaZn5+P4uJidOzYsUodnTp1QkFBAY4cOeI41rNnT5jNZqSnpzuVPXnyJMLCwhAdHe049tlnn6Fdu3aIj4/H0aNH8fXXX2Pv3r2wWCxVnqs29damLHkvDpahZuf8+fPIycnBfffdV+e6tm3bBgC4++67ayzz3//+F3l5efj973/vOGa323HhwgXcddddVcrHxMQAAM6dO3fL566uXneUvVFt4gwICIBcLkd+fn6Vsnl5eQCAtLQ09O/fHwDwpz/9Cenp6Rg9ejSeffZZBAcH47vvvsPJkyexbt06qFQqx+OPHTuGyMhI9OrVC4WFhbjjjjtw9OhR+Pj44KuvvsLgwYMdZWtTb23KkvdiIqRmp/JDOTQ0tNaP3bp1K06cOAGTyYSdO3fiyy+/xO9+9zsMGzas2vIXLlzA9OnTERERgdmzZzuOm0wm2O12+Pn5VXlM5TGDwVBjHDXVW9eyN6tNnBqNBsOGDcOOHTuQk5OD8PBw4Pr53rJlS5XXpNVq0atXL+zduxdff/01goODcfjwYQwePBht2rRxeq6CggLk5OTgkUcewYEDB6BUKlFcXIzhw4dj3LhxSEtLQ1hYWK3rrU1Z8l5MhNTsVH6Al5aW1vqxO3fuREBAAHx8fBAZGYkdO3bgnnvuqbZsdnY27r33XpSVlWHTpk2OxADAMeWiuq69ymO+vr61rrcuZatT2zgXLVqEhIQE9O/fH9OnT4dMJsOSJUsQFRWF4uJip4Q6e/ZsvPvuu9i8ebOjda7X6zF69GgMHDgQp0+fdkyNUKlUsFqtmD9/PpRKJQAgMDAQb7zxBh5++GGsW7cOzz//fK3rrU1Z8mINfZGSqK5uHuhSWloqfH19xdChQyXXcTt5eXnizjvvFBqNRvzwww/VlmnRooVISEiocvzbb78VAMSXX34pqV4pZW+ltnHm5eWJ9957Tzz55JPimWeeEevWrROffvqpACB+/PFHR7mgoCAxcODAKvUmJiYKAOKTTz5xHOvUqZOQyWTCarU6lf31118FAPHqq69Kqrc2Zcl7sUVIzY6vry8efPBBrF+/HhcuXEDbtm2rLVdcXCypNVBUVIRRo0YhPT0d3377LRISEqot169fP+zduxc2mw0KhcJxvHJASb9+/STVW9uyt1PbOENDQzFz5kynY5MnT0bLli3Rq1cv4Pq1R4PBUO0CB5XHioqKHMcGDBiAM2fOoLCw0NEFimq6uWtTb21jIC/W0JmYqK6qa82lp6eL4OBgMWDAAFFcXFzlMT/++KP4zW9+c8s6qqPX60WfPn2ERqMR33///S3LfvPNN1VaVCaTScTGxooRI0ZIrrc2ZV1Rmzhzc3OdpkgIIcThw4eFSqUSS5cudTo+YMAAERgYKK5du+Z0/A9/+IMAIA4fPuw4tn//fgFAvPfee05ln3rqKSGXy8Uvv/wiqd7alCXvxRYhNRs3LpnVrl07JCcn4/HHH0enTp3wxBNP4I477kBpaSl++OEHJCYm1nqEJQD85je/waFDh/D0008jLy8Pa9ascdynVqvx8MMPO24/+OCDmDp1KqZNm4a0tDS0bt0aK1asgMViwbJlyyTXW5uyrqhNnGfOnMELL7yA8ePHIzo6Gj/99BNWrFiBmTNn4tlnn3Uq++GHH+K+++5D//798fzzzyM4OBhJSUlYt24dZs2a5Wg94nqrc968eXjttddw8eJFdO/eHT/88APWrFmDf/3rX+jSpYukemtTlryXTAghGjoIorqYMWMGPvroI1y+fLnKvDCbzYbt27djz549uHbtGnQ6HWJiYjBmzBh07drVUW7t2rX4+uuv8a9//QutWrWq8bmmTp2KkpKSau/T6XRYvnx5leNbt27F1q1bHWt4PvXUUwgJCZFcr5QYXOFKnLg+PeWzzz7DhQsXEB0djfHjx6NHjx7V1llSUoLPPvsMp06dgtFoRExMDB588MEay+/fvx/r169Hbm4uWrVqhccffxw9e/asU721jYG8DxMhNXmDBg3Cr7/+iry8PMjlXCOCiGqHXaPUJJWXl+Orr77C7t27sW/fPrzzzjtMgkQkCRMhNUnl5eXYvn07QkNDsXXrVowZM6ahQ2pULl26hP3799+yTPfu3Z2uvRF5KyZCapK0Wi2+/PLLhg6j0bpy5Qq++eabW5ZRq9VMhES8RkhERN6OF1WIiMirMRESEZFX84prhHa7HZmZmQgICHCadE1ERM2TEAIGgwGtWrW67Yhyr0iEmZmZ3ICTiMgLXblyBVFRUbcs4xWJsHKB3StXrkCn0zV0OERE5GF6vR7R0dHVLrp+M69IhJXdoTqdjomQiMiLuHI5jINliIjIqzEREhGRV2MiJCIir8ZESEREXo2JkIiIvBoTIREReTUmQiIi8mpMhERE5NWYCImIyKsxERIRkVdjIiQiokbBarNDX2qB1Wav1+f1irVGiYioccsoMGLLsXTkG0oRGqDB/fFxaB3iXy/PzRYhERE1KKvNji3H0nEpVw9flQKXcvXYciy93lqGbBESEVGDMlmsyDeUIkznC51GDQDIN5TCZLFCp1F5/PnZIiQiogalVfkgNECDPL0Z+tIy5OrNCA3QQKuqn7YaEyERETUoH4Uc98fHIaaFDmaLDbEtdLg/Pg4+ivpJUewaJSKiBtc6xB9PJXSDyWKFVuVTb0kQTIRERNRY+Cjk9XJN8GZ1SrlGoxGpqanui4aIiKieSUqEer0eEyZMgE6nQ/fu3R3Hx48fj8OHD7szPiIiIo+SlAhnz56N3NxcnDhxwun41KlT8dZbb7krNiIiIo+TCSFEbR/UunVr7Nu3D7GxsZDJZKisorCwENHR0TAajZ6IVTK9Xo/AwEAUFxdDp9M1dDhERORhtfncl9QizM/PR1hYGABAJpM5jptMJqfbREREjZ2kRBgfH49vv/0WuCkRvvfee+jbt6/7oiMiIvIwSdMn5s+fj3HjxmHfvn0QQmD+/PnYsWMHDh48iOTkZPdHSURE5CGSEmFCQgJ27tyJd955B9HR0Vi6dCni4+ORkpKCPn361Lq+rKws+Pj4oEWLFlXuy8zMREFBgdMxX19ftG/fXkroRERETiQNlnGXxYsX491330VZWRnKysrQokULLF26FMOGDXOUee655/DVV18hOjracaxDhw5ITEx0+Xk4WIaIyLvU5nO/TivLHD58GKdOnQIAdOnSBb169XL5sTabDadOncKuXbsQExMDm82GP/7xj3jwwQdx7tw5p9bhqFGjsH79+rqESkREVC1Jg2WuXr2KQYMGoU+fPnjppZfw0ksvoXfv3hgyZAgyMzNdqkOhUGDRokWIiYlx3J41axb0ej2OHj3qVNZut+P8+fPIz8+XEi4REVGNJCXCp59+GkqlEmlpaSgsLERhYSHS0tKgUCjw9NNPSw7ml19+AQCnblAASExMxLBhwxAdHY0uXbpgz549kp+DiIjoRpKuEWo0Gpw5cwZt2rRxOn758mV06tQJpaWltQ6kuLgYffv2Rfv27bF582bH8S+//BL9+/dH27ZtUVZWhpdffhlffvklTp48WSVhVqq85lhJr9cjOjqa1wiJiLyExyfUt2rVCnJ51YfK5XK0bt261vWZTCaMGzcOarUan376qdN9kyZNQtu2bQEAarUaH3zwAYQQ+Prrr2us7+2330ZgYKDjp6aESUREJCkRTpkyBdOmTcPVq1cdx65evYpp06ZhypQptaqrtLQUDzzwAAoKCpCUlISQkJBbllepVAgPD8eVK1dqLDNnzhwUFxc7fm5VloiIvJukUaOffvop0tLSEBMTg5YtW0IIgezsbNjtdqSnp+Ozzz5zlD19+nSN9VQmwdzcXPzwww9V5hEKIWC1WqFUKh3HLl26hEuXLqFTp0411qtWq6FWq6W8NCIi8jKSEuGMGTPq/MRWqxUPPvggTp48iTVr1iA7OxvZ2dnA9UW9g4ODUV5ejj59+mDGjBno2rUrLl++jDfffBMdO3bEb3/72zrHQERE1GAT6ouKijBo0KBq73vrrbfwm9/8BgBw7tw5/Otf/8Lx48cRHByMwYMH46WXXoJWq3X5uTihnojIu9Tmc79BV5apL0yERETepV5WlklNTUViYiIuX74Mq9XqdN/q1aulVktERFSvJI0aTUxMRJ8+fXDo0CGsXLkSRqMRu3fvxieffFJlgWwiIqLGTFKLcN68eVi1ahUmTJgAmUyG9evXw2q1Yvr06U4T2YmIiBo7ySvL5OfnQ6vVQqlUQq/XQ6PRICcnB126dEFeXp5nopWI1wiJiLyLx1eWMZvNjlGbrVq1wtmzZwEA5eXlMJvNUqokIiJqEHXahgkAxo0bhylTpmDSpElITEzE0KFD3RMZERFRPZDUIty2bZvj9wULFmDgwIFYs2YN2rZtixUrVrgzPiIiIo/iPEIiImp2PH6NkIiIqLmQlAgvXryIsWPHIiQkBDKZrMoPERFRUyFpsMzUqVNht9uxbNkyBAcHuz8qIiKieiIpER48eBDp6emIiIhwf0RERET1SFLXaGRkpPsjISIiagCSEuG0adMwZ84cLqdGRERNnstdo507d3b8brfbkZaWhv/+97+IioqqMkDmVrvSExERNSYuJ8Lp06d7NhIiIqIGwERIREReTdI1wpycHCxatMhxe/ny5Wjfvj3uvfdeZGVluTM+IiIij5KUCGfNmoXQ0FAAQFZWFl566SVMmTIFKpUKM2fOdHeMREREHiNpHuH27dvx4YcfOn4fMWIEXn/9dWRmZuKuu+5yd4xEREQeI6lFaLFYHFMnkpKSMHLkSACAn58fLBaLeyMkIiLyIEmJcPDgwXj++efx4YcfIjExEWPHjgUAHDhwAP3793d3jERERB4jKREuXrwY5eXlWLFiBRYuXIi4uDgAwJIlS/DGG2+4O0YiIiKP4X6ERETU7HA/QiIiIhcxERIRkVdjIiQiIq/GREhERF6NiZCIiLyapER48eJFjB07FiEhIZDJZFV+iIiImgpJS6xNnToVdrsdy5YtQ3BwsPujIiIiqieSEuHBgweRnp6OiIgI90dERERUjyR1jUZGRro/EiIiogYgKRFOmzYNc+bMcSy8TURE1FRJWmKtY8eOSEtLg1arRVRUVJUBMqdPn3ZnjHXGJdaIiLxLbT73JV0jnDFjhtTYiIiIGhVJiXD69Onuj4SIiKgBcEI9ERF5NZdbhFFRUQCAq1evOn6vydWrV+seGRERUT1wORH+7W9/q/Z3IiKipowb8xIRUbPDjXmJiIhcxERIRERejYmQiIi8GhMhERF5tTolQqPRiNTUVPdFQ0REVM8kJUK9Xo8JEyZAp9Ohe/fujuPjx4/H4cOH3RkfERGRR0lKhLNnz0Zubi5OnDjhdHzq1Kl466233BUbERGRx0lKhBs3bsTKlSudWoMA0K9fPyQnJ7tcz8WLF/F///d/6Nu3LwYOHIjXXnsNBQUF1T7fiBEjcOedd2LixIk4d+6clLCJiIiqkJQI8/PzERYWBgBOWzCZTKYqWzLVxGaz4d5770VsbCw++ugjzJ8/H7t27cKIESOc9jncvHkzxo8fj7Fjx2LlypUQQmDw4MHVJkwiIqLakrSyzIABA/Diiy9i4sSJUCgUsNlsAICZM2fi559/RlJSkkv1WK1W+Pj8b5W3M2fOoHPnzti9ezeGDBkCAOjVqxe6deuGVatWAQAsFgsiIyMxc+ZMvP766y49D1eWISLyLh7fj3D+/PkYN24c9u3bByEE5s+fjx07duDgwYO16hq9MQkCQGVOrmxVGgwGHDt2DLNmzXKUUalUGDlyJHbt2uVyIiQiIqqJpK7RhIQE7Ny5E1lZWYiOjsbSpUsRHByMlJQUDBgwQHIwb7zxBmJjY9GnTx/g+i4WQgi0bNnSqVzLli1vucNFWVkZ9Hq90w8REVF1JLUIAaBPnz7YsGGD2wJ58803sXXrVvzwww9Qq9UAALvdDgBQKpVOZVUqlaM7tjpvv/025s2b57bYiIio+WoUK8ssWLAA//znP7F582b07dvXcbxyQE5+fr5T+by8PMd91ZkzZw6Ki4sdP1euXPFg9ERE1JRJahHeKgmp1Wq0a9cOU6ZMwVNPPXXbut555x289dZb2Lx5MxISEpzui4iIQHR0NPbv349x48Y5ju/btw+jR4++ZQyVrUoiIqJbkdQifPHFF1FaWorHH38c77zzDv7xj3/gscceg8lkwhNPPIF+/fphxowZWLFixS3reffddx1JcPjw4dWWee6557BixQqcPn0aALB8+XKcP38e06ZNkxI6ERGRMyHBsGHDxDfffFPl+IYNG8SwYcOEEEKsXbtWdO3atcY6CgoKBADh7+8v4uLinH7Wrl3rKGe1WsXzzz8vVCqVCA4OFmFhYeKrr76qVbzFxcUCgCguLq7V44iIqGmqzee+pHmEOp0OGRkZCAgIcDqu1+sRHR2N4uJi6PV6tGrVCkajsdo67HY7zp8/X+194eHhVeZ9lJSUoKCgAJGRkVWmXdwO5xESEXkXj88j1Ol02Lx5MyZOnOh0fNOmTY4nzMjIQFxcXI11yOVytG/f3uXn9PPzg5+fn5RwiYiIaiQpEf7pT3/C5MmTsXnzZvTq1QtCCBw9ehTr1q3D+++/DwB47733nCbCExERNUaSukYBICkpCR988AFOnToFmUyGzp0746WXXsLIkSPdH2UdsWuUiMi71OZzX3IibEqYCImIvEttPvcbxYR6IiKihiIpEV68eBFjx45FSEgIZDJZlR8iIqKmQtJgmalTp8Jut2PZsmUIDg52f1RERET1RFIiPHjwINLT0xEREeH+iIiIiOqRpK7RyMhI90dCRETUACQlwmnTpmHOnDkoKytzf0RERET1SNL0iY4dOyItLQ1arRZRUVFVBshULpDdWHD6BBGRd/H4EmszZsyQGhsREVGjIikRTp8+3f2REBERNQBOqCciIq8mqUUIAKmpqUhMTMTly5dhtVqd7lu9erU7YiMiIvI4SS3CxMRE9OnTB4cOHcLKlSthNBqxe/dufPLJJygoKHB/lERERB4iqUU4b948rFq1ChMmTIBMJsP69ethtVoxffp0TqkgIqImRdL0CY1Gg/z8fGi1WiiVSuj1emg0GuTk5KBLly7Iy8vzTLQScfoEEZF38fjuE2azGVqtFgDQqlUrnD17FgBQXl4Os9kspUoiIqIGIXmwTKVx48ZhypQpmDRpEhITEzF06FD3REZERFQPJLUIt23b5vh9wYIFGDhwINasWYO2bdtixYoV7oyPiIjIo7hDPRERNTvcoZ6IiMhF3KGeiIi8GneoJyIir8Yd6omIyKtxh3oiIvJq3KGeiIi8mstdo507d3b8brfbkZaWhv/+979NYod6IiKimricCLkZLxERNUdMhERE5NUkXSMUQuD8+fNVjp8/fx5esFANERE1I5IS4cKFC7F48eIqxxctWoQPP/zQHXERERHVC0lrjcbGxmLPnj1o06aN0/FLly5hxIgROHfunDtjrDOuNUpE5F08vtZoTk4OfH19qxz39fVFRkaGlCqJiIgahKREGB8fj+XLl1c5vmzZMvTs2dMdcREREdULSUus/fWvf8U999yDH3/8EUOGDIEQAnv27EFSUhJ27Njh/iiJiLyc1WaHyWKFVuUDHwU3DnInSYkwISEBe/bswYIFC7B48WLIZDL07NkTe/bsQf/+/d0fJRGRF8soMGLLsXTkG0oRGqDB/fFxaB3i39BhNRvcmJeIqBGz2uz4ODkVl3L1CNP5Ik9vRkwLHZ5K6MaW4S3U5nNfUouQiIjqh8liRb6hFGE6X+g0agBAvqEUJosVOo2qocNrFvh1goioEdOqfBAaoEGe3gx9aRly9WaEBmigVbEd4y5MhEREjZiPQo774+MQ00IHs8WG2BY63B8fx25RN5L0laKkpAR+fn7uj4aIiKpoHeKPpxK6cdSoh0jemPfpp5/Gjz/+6P6IiIioCh+FHDqNiknQAySd0WXLluHSpUsYPHgw7rjjDvzzn/9Edna2+6MjIiLyMEmJcOLEifj+++9x/vx5TJgwAUuWLEFUVBQefPBBbNy4EVar1f2REhEReYBb5hEKIbBo0SLMmjULFosFLVu2xPTp0/HKK69UuyZpfeM8QiIi71Jv8wj1ej3WrFmDjz/+GIcPH8aoUaPw9NNPIycnBwsXLsTBgwexadOmujwFERGRR0lKhLt27cLHH3+MDRs2IDQ0FFOmTMHatWudtmV69NFHERUVddu6DAYDvvjiC5w+fRrTp09H+/btne5ft24d9u3b53QsMjISr732mpTQiYiInEi6RnjPPffAaDRi/fr1uHjxIubNm1dlb8IWLVpg2rRpt6xnxYoV6NSpE5KTk/HBBx/g6tWrVcrs3LkTu3btQmxsrOOndevWUsImIiKqQlKL8OrVqwgPD79tuUWLFt3y/t69e+PMmTMoLi7G2rVrayzXvn17vPzyy1JCJSIiuiVJidCVJOiKHj16AACKi4tvWS49PR1/+tOfEBgYiMGDB2PAgAFueX4iIiJJXaMXL17E2LFjERISAplMVuXHnWQyGUJCQuDr64vz589jxIgRePHFF2/5mLKyMuj1eqcfIiKi6khqEU6dOhV2ux3Lli1DcHCw+6O6wezZsxETE+O4PWHCBAwfPhwPPvggRo4cWe1j3n77bcybN8+jcRERUfMgKREePHgQ6enpiIiIcH9EN7kxCeL6psDR0dHYu3dvjYlwzpw5mDlzpuO2Xq9HdHS0x2MlIqKmR1IijIyMdH8ktWCxWFBeXl7j/Wq1Gmq1ul5jIiKipknSNcJp06Zhzpw5KCsrc39EN7Bardi9e7fTsS+//BLZ2dm45557PPrcRETkHSQtsdaxY0ekpaVBq9UiKiqqygCZ06dPu1TPoUOH8OWXX8JoNGLlypV45JFHEBUVhdGjR2P06NGw2Wy47777YDAY0LVrV1y+fBl79uzB66+/jrlz57ocL5dYIyLyLh5fYm3GjBlSY3Pi7++P2NhYAMDChQsdx4OCggAACoUCO3bswLFjx3D8+HGMGTMGK1eudGnFGiIiIle4ZdHtxo4tQiIi71Kbz33u8EhERF5N8u4TqampSExMxOXLl6vsP7h69Wp3xEZERORxklqEiYmJ6NOnDw4dOoSVK1fCaDRi9+7d+OSTT1BQUOD+KImIiDxEUotw3rx5WLVqFSZMmACZTIb169fDarVi+vTpHp9SQURE5E6SBstoNBrk5+dDq9VCqVRCr9dDo9EgJycHXbp0QV5enmeilYiDZYiIvIvHB8uYzWZotVoAQKtWrXD27FkAQHl5Ocxms5QqiYiIGoTkwTKVxo0bhylTpmDSpElITEzE0KFD3RMZERFRPZDUIty2bZvj9wULFmDgwIFYs2YN2rZtixUrVrgzPiIiIo/ihHoiImp2OKGeiIjIRUyERETk1ZgIiYjIqzEREhGRV3NrIlyxYgXOnDnjziqJiIg8yq2JcNq0aejWrRtefvlld1ZLRETkMW5NhAaDAUePHkXLli3dWS0REZHH1HllmRv5+/ujW7du6NatmzurJSIi8pg6tQiNRiNSU1PdFw0REVE9k5QI9Xo9JkyYAJ1Oh+7duzuOjx8/HocPH3ZnfERERB4lKRHOnj0bubm5OHHihNPxqVOn4q233nJXbERERB4naa3R1q1bY9++fYiNjYVMJkNlFYWFhYiOjobRaPRErJJxrVEiIu/i8bVG8/PzERYWBgCQyWSO4yaTyek2ERFRYycpEcbHx+Pbb78FbkqE7733Hvr27eu+6IiIiDxM0vSJ+fPnY9y4cdi3bx+EEJg/fz527NiBgwcPIjk52f1REhEReYikFmFCQgJ27tyJrKwsREdHY+nSpQgODkZKSgoGDBjg/iiJiIg8hBvzEhFRs+PxwTJCCJw/f77K8fPnz8ML8ioRETUjkhLhwoULsXjx4irHFy1ahA8//NAdcREREdULSV2jsbGx2LNnD9q0aeN0/NKlSxgxYgTOnTvnzhjrjF2jRETexeNdozk5OfD19a1y3NfXFxkZGVKqJCIiahCS5xEuX768yvFly5ahZ8+e7oiLiIioXkiaR/jXv/4V99xzD3788UcMGTIEQgjs2bMHSUlJ2LFjh/ujJCIi8hDJ8wj37NkDpVKJxYsX49///jdUKhX27NmD4cOHuz9KIiIiD+E8QiIianY8PliGiIiouZCUCC9evIixY8ciJCQEMpmsyg8REVFTIWmwzNSpU2G327Fs2TIEBwe7PyoiIqJ6IikRHjx4EOnp6YiIiHB/RERERPVIUtdoZGSk+yMhIiJqAJIS4bRp0zBnzhyUlZW5PyIiIqJ6JGn6RMeOHZGWlgatVouoqKgqA2ROnz7tzhjrjNMniIi8S20+9yVdI5wxY4bU2IiIiBoVSYlw+vTp7o+EiIioAXBCPREReTVJLUIASE1NRWJiIi5fvgyr1ep03+rVq90RGxERkcdJahEmJiaiT58+OHToEFauXAmj0Yjdu3fjk08+QUFBgfujJCIi8hBJLcJ58+Zh1apVmDBhAmQyGdavXw+r1Yrp06dzSgURETUpkqZPaDQa5OfnQ6vVQqlUQq/XQ6PRICcnB126dEFeXl6t6rt69SrOnTuHu+66C0FBQdWWuXDhArKzs9GpU6daL+vG6RNERN7F47tPmM1maLVaAECrVq1w9uxZAEB5eTnMZrPL9Rw6dAgPPfQQevXqhYSEBJw4caJKmdLSUjz44IPo1q0bpk2bhlatWmHhwoVSwiYiIqqizqNGx40bhylTpuDdd9/FY489hqFDh7r82DNnzmDy5Mk4ePBgjWXefPNNHD9+HOnp6UhNTcXatWsxc+ZMHDhwoK6hExERSUuE27Ztc/y+YMECDBw4EGvWrEHbtm2xYsUKl+v5/e9/j4cffhgKhaLGMqtXr8bTTz/tWOB77Nix6N69O1atWiUldCIiIieSEuGNXZh+fn746KOPcOTIEXz++ef45JNP3BZcRkYGcnJycPfddzsdv/vuu3H8+PEaH1dWVga9Xu/0Q0REVB1JiXDOnDmS7qutwsJCAEBISIjT8bCwMMd91Xn77bcRGBjo+ImOjnZbTERE1Ly4dWWZtLQ0hIWFua0+lUoFXB8wcyOTyeS4rzpz5sxBcXGx4+fKlStui4mIiJqXWs0jjI2NrfZ3ALDb7cjOzsbkyZPdFlxUVBTkcjkyMjKcjmdkZCAmJqbGx6nVaqjVarfFQUREzVetEuHcuXOB6/sRVv5eSalUIjY2FkOGDHFbcFqtFgMHDsSmTZvw+9//HgBgNBqRlJSEt956y23PQ0RE3qtWifDpp58Grl+je+ihh+r85NeuXcPp06eRm5sL3DAIJzY21tHinD9/PkaMGIFXX30V/fv3x6JFi9CyZUtMmzatzs9PREQk6RrhgAEDsGjRIsft5cuXo3379rj33nuRlZXlcj0nTpzAm2++icWLF2Po0KH45ptv8Oabb2LXrl2OMoMHD8aePXuQnZ2NZcuWIT4+Hvv27YO/v7+U0ImIiJxIWmLtiSeewJgxYzBx4kRkZWWhXbt2mDt3Lg4cOAB/f3989dVXnolWIi6xRkTkXTy+Q/327dvx4YcfOn4fMWIEXn/9dWRmZuKuu+6SFjUREVEDkNQ1arFYHLtMJCUlYeTIkcD1yfUWi8W9ERIREXmQpEQ4ePBgPP/88/jwww+RmJiIsWPHAgAOHDiA/v37uztGIiIij5GUCBcvXozy8nKsWLECCxcuRFxcHABgyZIleOONN9wdIxERkcdIGizT1HCwDBGRd/H4foRERESeYrXZoS+1wGqz18vzuTxqNCoqCri+m3zl7zW5evVq3SMjIiKvk1FgxJZj6cg3lCI0QIP74+PQOsSz88ZdToR/+9vfqv2diIjIHaw2O7YcS8elXD3CdL64lKvHlmPpeCqhG3wUnuvAdDkRPvnkk9X+TkRE5A4mixX5hlKE6Xyh01RsnJBvKIXJYoVOU/OOQ3XFa4RERNQoaFU+CA3QIE9vhr60DLl6M0IDNNCqJK394jKXa6/NPoN5eXlS4yEiIi/lo5Dj/vg4xzXC2BY63B8f59FuUdQmES5dutTx+y+//IIFCxZg8uTJ6N27NwDg8OHD+OSTTzB79mzPREpERM1e6xB/PJXQDSaLFVqVj8eTIKTOI0xISMD06dPxyCOPOB1fv349Fi9ejOTkZHfGWGecR0hE5F1q87kvKRHqdDpcvXq1SuV6vR7R0dEoLi6ufdQexERIRORdPD6hXqfTYdOmTVWOb9q0CYGBgVKqJCIiahCShuK88cYbeOqpp7B161b07t0bQggcOXIE69evx5IlS9wfJRERkYdISoTPPvssOnbsiPfff9+R+Lp06YLvv/8eQ4cOdXeMREREHsNFt4mIqNnhottEREQuYiIkIiKvxkRIRERejYmQiIi8GhMhERF5NUmJ8OLFixg7dixCQkIgk8mq/BARETUVkuYRTp06FXa7HcuWLUNwcLD7oyIiIqonkhLhwYMHkZ6ejoiICPdHREREVI8kdY1GRka6PxIiIqIGICkRTps2DXPmzEFZWZn7IyIiIqpHkpZY69ixI9LS0qDVahEVFVVlgMzp06fdGWOdcYk1IiLvUpvPfUnXCGfMmCE1NiIi8hCrzV6vO7s3F5IS4fTp090fCRERSZZRYMSWY+nIN5QiNECD++Pj0DrEv6HDahL4lYGIqImz2uzYciwdl3L18FUpcClXjy3H0mG12Rs6tCbB5RZhVFQUAODq1auO32ty9erVukdGREQuMVmsyDeUIkznC51GDQDIN5TCZLFCp1E1dHiNnsuJ8G9/+1u1vxMRUcPSqnwQGqDBpVw9ACBXb0ZsCx20KklXv7wON+YlImoGeI3QmcdHjVYyGo24cOECunXrVpdqiIiojlqH+OOphG4cNSqBpDOl1+sxYcIE6HQ6dO/e3XF8/PjxOHz4sDvjIyIiF/ko5NBpVEyCtSTpbM2ePRu5ubk4ceKE0/GpU6firbfecldsREREHifpGmHr1q2xb98+xMbGQiaTobKKwsJCREdHw2g0eiJWyXiNkIjIu9Tmc19SizA/Px9hYWEA4LS8mslk4n6ERETUpEhKhPHx8fj222+BmxLhe++9h759+7ovOiIiIg+TNGp0/vz5GDduHPbt2wchBObPn48dO3bg4MGDSE5Odn+URETUYJr7GqaSXlFCQgJ27tyJrKwsREdHY+nSpQgODkZKSgoGDBjg/iiJiKhBZBQY8XFyKhZvP4aPk1ORUdC4xoC4g6REuGDBAvTp0wcbNmzApUuXcOXKFWzcuBF9+vTBggUL3B8lERHVO29Zw1RSIpwzZ46k+4iIqOm4eQ3TMJ2vYw3T5sStnb1paWmO0aRERNS0Va5hmqc3Q19ahly9GaEBmma3hmmtXk1sbGy1vwOA3W5HdnY2Jk+e7L7oiKjRaO4DJqgqH4Uc98fHOdYwjW2hw/3xcc3u71+rCfUrVqwAAEybNg3/+c9/nO5TKpWIjY3FkCFDGt1cQk6oJ6obLujs3ZrilyCPLbr99NNPAwDCwsLw0EMP1S1KF7388suOBFypS5cuOHToUL08P5G3u3HARJjO1zFg4qmEbk3mQ5HqpnIN0+ZK0rt4wIABWLRokeP28uXL0b59e9x7773IyspyZ3wwm80YNWoUrl275vjZvXu3W5+DiGrmLQMmyHtJSoSzZs1CaGgoACArKwsvvfQSpkyZApVKhZkzZ7o7RigUCvj7+zt+NBqN25+DiKrnLQMmyHtJSoTbt2/HmDFjHL+PGDECr7/+OpYtW4adO3e6O0YkJSUhPDwcHTp0wFNPPeX2VicR1axywERMCx3MFluzHTBB3kvSVzqLxYKysjLgepIaOXIkAMDPzw8Wi8WtAUZHR+Pf//43hg4dioyMDMyaNQuDBg3CTz/9BH//6i/Wl5WVOeLD9YumRCQdN32l5kzSNkxjx46FUqnEsGHDMHv2bKSmpiIuLg47duzA+++/j23btnkmWgB5eXlo3bo1Fi9e7Bi8c7M333wT8+bNq3Kco0aJqDFriqMzGyuPb8O0ePFilJeXY8WKFVi4cCHi4uIAAEuWLMEbb7whLWoXhYWFoU2bNjh79myNZebMmYPi4mLHz5UrVzwaExFRXTX3NT2tNjv0pZZGuTybpK5RnU7n2IbpRhs3bnRHTLdkNBqRkZGB8PDwGsuo1Wqo1WqPx0JE5A7NfYpKY5+HKukMh4SEuD+SapSVleH3v/89fv31V9hsNly5cgW/+93voFarMWnSpHqJgYjI05rzFJWmsHC3pETYsmXLehm5qVar8cADD+B3v/sd/Pz80KNHDwghsH//frRq1crjz09EVKm6rj13dfc15ykqTSHJSxos8/bbbyM1NRVLly6tt8EnNpsNCoVC0mO5xBoR1UV1XXsA3Nrd19i7D6Wy2uz4ODnV0e2bqzcjtoXO492+tfncl5QIO3bsiLS0NCiVSkRFRUGlcl565/Tp07WP2oOYCIlIqps/yPP0ZkSHBUAI4Gq+wXEsxg0f7s111GhDJHmPrTVaacaMGVJjIyJqUm7u2gOA7KISyGQyp2OV3X11WZOzua7p2djnoUpKhNOnT3d/JEREjVDl9btLuRULc+TqzWhzQ4uw8lhsC12zuKbnKY05ydcpLRuNRqSmprovGiKiRqa6JebG3t0e43q157JzzYSka4R6vR7PPPMM1q5dCyEEKqsYP348XnvtNfTu3dsTsUrGa4REVFfVXb9rrtf0mgOPrywze/Zs5Obm4sSJE07Hp06dirfeektKlUREjVpl196NCa+6Y9T0SGoRtm7dGvv27UNsbCxkMpmjRVhYWIjo6GgYjY1raSC2CIkaF7akyNM8Pmo0Pz8fYWFhAACZTOY4bjKZnG4TEd2suc6Xo6ZL0lex+Ph4x1qjNya+9957D3379nVfdETUrDSF5bbI+0hqEc6fPx/jxo3Dvn37IITA/PnzsWPHDhw8eBDJycnuj5KImoXq5uS5Y/4dUV1IahEmJCRg586dyMrKQnR0NJYuXYrg4GCkpKRgwIAB7o+SiJqF5rymJjVdkgbLNDUcLEPuxsEe0jWla4T8OzddHh8sk5OTg7Vr1zpWmFm+fDn+8Y9/IC4uDqtXr0ZkZKS0yImagKb0Qd4YNfbltirx7+w9JL0DZ82ahdDQUABAVlYWXnrpJUyZMgUqlQozZ850d4xEjUZjHuzRmHcAv1l18+8aU/yN+e9M7iepRbh9+3Z8+OGHjt9HjBiB119/HZmZmbjrrrvcHSNRo9FYB3s09dZLQ8R/q27Pxvp3Js+Q1CK0WCwoKysDACQlJWHkyJEAAD8/P1gsFvdGSNSINMbBHk299WK2WPH1wbO4kFNcb/FnFBjxcXIqFm8/ho+TU5FR4LwISGP8O5PnSEqEgwcPxvPPP48PP/wQiYmJGDt2LADgwIED6N+/v7tjJGo0qluAuaEXW24KO4DXJKPAiP/s/Bn7zlxFTrEJgMzj8bvyxaEx/p3JcyR9vVm8eDFeeOEFrFixAgsXLkRcXMVuzUuWLMEbb7zh7hiJGpXGNtijum2CGuOWQDd3RVYmpMxCI/zUSuQZTDh5OQ/hgVq0DQ/0WPyudns2tr8zeY6kd1qbNm0cK8vcaOPGje6IiajRa0x7q1W2XiqvsTXG1kt11wADNCrkG0oRHqhBeKAWJy/noqTMglbB4R6NvzZfHDzxd+aUjMZH0l8hJycHixYtctxevnw52rdvj3vvvRdZWVnujI+IXFDZenlhdDyeSuhWZaCJlBGZ7hrFWVNXpEohd1yHAwTCA/0wsFMUpo3o7tGBMg3Z7Xm7a5PUMCRNqH/iiScwZswYTJw4EVlZWWjXrh3mzp2LAwcOwN/fH1999ZVnopWIE+ql47fXpq+yNZanL0WARoUHe7VHmxa3/n/gyihOV98b+lILFm8/Bl+VAjqNGvrSMpgtNrwwOh6GUkuDjXat7/e21WbHx8mpuJSrR5jOF3l6M2Ja6PBUQjf+3/IAj0+o5/QJ79DUh+Q3Ve78gK5sjZ3OKECxqQxFJWb8cjUPf3q4H9qEVf/hcGMLLkzn62jB3fiBXZv3xq26InUaFZ4Y0hX5RjNC/X3hW4/XNeu7e5tTMhovTp+gajX1IflNlbu7zkwWK/L0pSgsMaOkrBxatRK5xSZsPHyuxr/l7UahVvfe+PboORQazdXWeauuyIwCIz7d8ws+3X0Sn+75pc6vtzFNyr8Zp2Q0XpL+ApXTJ4YNG4bExETHrvScPtF88Ntr/XOlJVZbWpUPZLKKBCsDoJDLEBrgC31pWY1/y9sNJrn5vWE0l2P/2UzkFJvQMsgPo3u0RZC/r1OLtroRmO5+vY29B6MpDGryVpw+QdVqKkPy3a0hr4nemGD81SpYNHbk6av/8lGbOH0UcijkgNUmAAGYy20I8ddApZBDX2qpUsftPrBvfG/YhcCvV/MBAH5qJU5nFODYhWxEhwaghU7rlIwquyIrW202m91tX7asNjs2HTmHCznFiAjUuuVLhCdwSkbjxOkTVC1v/Pba0C2KygRz4/W8FoFaFJWYnRJDbeI0Wayw2uzo17EVMvINMJrLofRR4M42Yfh0zy811nGrD+wb3xuZBQb4yOWIaxkEtVKBwhIz9KYytA3X3fbaYrC/L1Q+iuujRuv2ZSs9uwgH0jJRbrXBaLYgMtiv0fZgNKapN1Sh+X6qUZ3dbkh+c9IYron6KOQY3aMtTJZyFJWYEeSnhlalxPYTFxxx1DbOyuRqKbejU+sQtA7VoW+HSJy8nHfbOqpbGLtS6xB/jO7RFsH+GljtdhxMy8KuXy4jo8AIrdoHIX6a215bvJJnAABEhQbUaRqD1WZH8i+XYbHaYLUJ5OpL8evVfAT5+Tb7HgxyD8nvktTUVCQmJuLy5cuwWp2XQlq9erU7YqNGwFu+vTaWa6JB/r6IDg1A23AdQvw0MJZZnOKobZw3t+w7tAzC0C7R+Prg2Tq9VqvNju0/XcC1ohL4yOUot9pQAkAur+h61ZdakG+89bVFADBbbHhqeDco5HKnlmdtun5NFisKjWZ0iQpFVoERxaUW+CjkGH5nm2bdg0HuIykRJiYmYtKkSRgxYgS2bNmCRx55BEePHsXFixcd644S3aixz0dsLNdEtSoftNBVXONS+SiqxCElzpu7OQHU+bVWJrUgrRpZhUZEBvsBAGJa6HAl34iSsvJbXlu88XkDfJ1bnbXtor6x3rjIIFwrKkVcRCDiIoJcfj3k3SR9Is2bNw+rVq3C5s2bAQDr169HWloann32Wcc+hUSVmsJqGg212sjNw/1vjiM6NABDu0TXGOfN99/q9VV2c7rjtVYmnyJTGRTyikE3Kh8FbHagf8dWeLGa7nRXnldKF/WN9VrK7ejQMghj727fKL9wUeMkaWUZjUaD/Px8aLVaKJVK6PV6aDQa5OTkoEuXLsjLy/NMtBJxZZmG09RW06jPluutWj5Wmx3p2UVI/uUyCo1ml++PCNS6HL/VZoeh1ALIUKVVVpv4r+QZkGswISxAgzZhutu24G51jm+1Cs3tum0be68D1S+PryxjNpuh1WoBAK1atcLZs2fRo0cPlJeXw2w2S4uabqmp/idvLNfeXFWba6J1+Zu4Modu969XcCXP4NL9F3KK8enukwjUqlFsKrttl2JlIv3h5GUUlVRNtK64sctVpZDDYrO7PJ2jpnNcly5qb7meTe5X5wsg48aNw5QpUzBp0iQkJiZi6NCh7omsGavtB6iUYf2NJXE2lmtv7pZRYMSmI+eQqzehhU6Lcb3a1yqJ3O4LQm3uB2TILjLi+MUSBGrUuLNN2C3n0WUUGPHt0XPYfzYTFqsNXaJCJc+7uzH5+Lr8qFvX523TdqjhSfo02rZtm+P3BQsWYPbs2VizZg06d+6Mf/7zn+6Mr9lxNalVJjKVQl7r1Tcaej7cjer7g81qs8NgtgACCKhh6L87nmPNvlM4cv4aZABOZeTDVFaOl+672+Xnq+4LQpuwANhsdlivt6xu9QWi8v4LOcXIKS5BrqEUcshgEwJZhSWIiwiqtuVd2RI9n10Mq80Oq+3W5RsCJ51TfatVInznnXcwfPhwjBo1ynHMz88PH330kSdia3ZcXVLqxkTm76tCdnEJQvx9ofbxQYi/7y0/sNy1bJU7W5T19cGWUWDEmn2n8PPlXABA9zYtMGHgHbfsHnQlppvLGUot+PlyLqw2gUCtCsWmituGUguC/V1rF/ko5Bh9V1tsPHwO+tIyhPr7osRcjqXfn3B8ebk/Pg7fHj2H7KIStAkLcPoCUfkF4+uDZ3EhpwgtArSQASgutaCoxIzsYhPatwyq0vKubElGBGlgLC1DnsF8y/INhd2cVJ9qnQhnz56NwMBADBkyBMOHD8fw4cPRrVs3yGQyz0XZTLhyvezmRJZZaER2UQnOXSuEXCaDANCrXcsaP7DccU3OEy3K232w1TXxWm12fHv0HI6cvwabTUAAOHL+GrRqJaaN6F5t96Arr7G6clp1xbmvfMfLAAgBGMwWl1uhGQVGbD6WjuwiI0ICNLDaBIpNZU5fXkb3aAshAJlMhuqGtLUO8ce0Ed0BAOdzipCnL4XeVAaVUoEWuopYcX0ASuV5vbGlGRnij0JTGZQ+CrQND2QXJHmtWr3r8/LycOzYMfz5z3+GEAJvvvkmevTogfDwcDz66KNYsmQJTp8+7blom7jKD6Gc4lLk6UuRU1zqWH2+chi9wWxxWm9Sp1Gh1GJ1fCC6+hw1rXB/u9X5G2KFlcrpFR9tO4p/f3cCl/P0Lj3uxtdisliRXVQCuUwGnVaFQK0KMgA5xSXILjY5xW+2WK+3pIqrvMYb66w8FxdyiiGXyXA+uwhfHzwLlUKOrlFhkMkAU1k5rHY7hBBYlnQC/9p8GKmXcqqcL7PFiowCI8zXlzxbufNnfHvkHA6nX8N3P13E4fQshASoHbs95OpN2HjkHK7mGxCgUeJqvqHav4OvygcP9mqPUosVprJyRIX6o114EAI0Ktjs9irTVm6cauAjl2Nol2j86eF+Ht8Ml6gxq1WLUC6Xo2fPnujZsydmzpwJm82GY8eOYdeuXUhOTsarr74Kk8kECTMyvIKPQo5e7Vri6PlrSL9WiGB/XzzarxOyi02OVkeQX8X6i1fyjCg2laHAWIpymx0924UjxE8Dq92GsnIbsotNiAjUVpmDZbJYMbpHW2z/6UKVa3I3r/OY0LUN4iKCnOqoS4tSSquuMtmcupqPwpJSFJeU4efLOXhtXB+EB/lXW1fliMedqZeQZyitWNy5ZztEBPnhdGbFOp12AchkQGZhCVbvSnUsAG2z2bHuwBkcOpcFrVqJID+1o7v5bGYBdp68DEOpBWE6DYZ2icaVPANyik24nKeHvrQMadcKYTRbUG61QyYDCq6fmxwBXMgphgCw7cQFDL0jGtNG9kDrEH8cO5+N5Tt/ur5smi8e698Rh9KzYLHaoVEpUGqxocBgRlZBRSLP1ZsRGeR3PQ7nv4PBbKmyCkt1q9EUGMzYePgcMguNCPJT4UJOsaOLnNfgiJzV6YKA0WhEdna248dsNqNTp07ui66ZqVwTsdhUBgGgwGjGjp8vQOWjwNU8A1qF+ONqvgFBWjVKyiqu9QT7+6LQaMbhtCwEaisSYaBW7fTh3jrEv0oX3ui72jrWWvRRyFFitmDNvlPI1Zvg56vE7l+vYP/ZTPTr0MppxKPUUZ616U69cSBQvtGM9OxCpF0rhLncCrtdoNBUhlc+3Y1ubcLQpoXzvLTLuXokHk7DgbRMFJWUIUCjglwmQ1GJGRMGdkZOkQmH0rNgtdmhVCggl8ExoCXPYML57CLk60thsdmRZyjFtaISBGhU6NwqBO9vPYpcQykCNSpkFalgNFtwrdiIbH0J7DY7Si02+PsCJy7mQCYDFDLAWGbFzV/7LFY79qdlokWgFk8O7YrlO39CVqERgddXYflszymUW22Qy2RQyOWQy+yQyYGIID/HRPPRd7XF9hMXnP4Oof6+WLPvdJXpDtWtRhMZ5IdrRUbkFJcgq9AIhVwGhVzm+ELDa3BE/1OrCfVGoxF79+5FcnIykpOTcfz4cXTq1AlDhw51/ERERHg2YgncOaH+5lZPbVpBeXoTpn+chHyDGVa7HTabACAqujxlQKBGjS5RIVAqfCCTCfj5KhGk9UXyL5eRVWSCSiGH1W5HkFaNQXe0Rp7ejFbB/ngq4U58sfdUjZPWj53Pxr+/O4GLucXwVfkgSKtCWbkdPgoZokJ1aN8ySNLu4zcmtE/3/OLSpPmbJ2H7+iiQejkXZuv/uvzk11vPMS10iAz2R9vwQDyV0A2ZhUb8PfEArhVWrCcJAFqVsqK7UCZwb4+2sNoEcopL4KtU4OfLubBY7dcXY7bDDsBHLoOfWomychvKr3czKq+vtlJus0Muk8FitQEQUPrIIZfJIbue3NRKBYL81JABsAtAbyqD3lxe7d9aq/JBvw6R+N2QrvjL2r1Q+SgQoFHBUGqBudwGlY8c14pKoJDJYRN2tAnV4aOnRsAm4Hgv3dyCLzGXo8BorvYc3/w3G9UtBu9sOuRIwMWmMkQG++P9ycPrdRd4oobisQn1wcHBiIyMxMMPP4zZs2djyJAhCAsLq2u8TcbNHza92rXEkfPXakwYNybJ7GIT1u4/jWtFJSi3iYoBFpUFRcXtQqMZB9KuISrED1q1Cr5KH1jtdhQYzPD1USAiSIvMQiPMVhvKrXbkFJfgQk4RLFYb8o2lCNSq4K9WAbr/dWeqFHIsSzqBrGIjlAoZDKVlMJotCPHzRZjODxGBWqeuT6vNjgCNCk8M6eo0QfrmaQk3duf6+6qQVWR0ev5cvcnRfWu12ZFvNCNIq3Jcc8spNuFaUcWalOU25+9i4noitAuBIK26okuw1IKNh88hp9gEtdIHdlMZIIASs8XRDbrnlysQMhkCfH1gtwsYSy0ou6nucptASZn1+nW9itaij1zAaC6HXFYx6KUyJdvK7ZBDwFepgPL6OdAoFSgtt0MIO8qt1V83lQFQKGRoGeyPqBB/BPn5IquwYlm5yoT05NCuWLXrFxQaK1r9z47sAT9f5xbajV2YNpsdS78/UWOX9c3dnSaLFS0CtCi32mGz2xEWoEWLAC0sNrtb5vsRNSe1SoR33303jh49iu+++w4WiwXl5eUYNmxYo2wFutvNozkv5BTj2IVsBPv5IjxQU2Wawo1JM8jPF6ayclwrKoFKoUC5rWp3GmQVrYzKLriSslKU22zIN5igUMihUfrARy6HHDKUl1txKiMfRSVlCAvQ4nKeHhdz9ZDJgBD/ig/KuIggqBRy/Ho1H5fy9LDZBOz26z9CwGYXFXu23bBDQEaBERsPpyGrsASRwX544O44AL4oKjFj3f4zjmkJd0aFATIgV18xDD89uwgXcysGlAT5qeEjlwMyYMXOnyAEcClPD6PZgkCtGkFaX7TQVZwvUzVJENcToclixeXcYqh8FLi7XUvYhB2X8/QoK7dCb7IBArDahdODco0Vqxpd392nRuU3DDiRyQCrqBhlenMo8uv3+/jIoVUpYS63IjzIH+E6raObG1Zblfr9fZUY0KE1HurdAX6+KjwzoofjGmFksD+eGdED8e0icFdsBPKNFV2eNbXSbtzM9nZd1jd2d2pVPogOC3D0IBSZyhAdFtBopkcQNSa1+l9x4MABGAwGpKSkIDk5Gf/85z/x29/+Fh06dMCwYcMcP80xMd48iMRSbseFnCK0DddBp1HDLgSyCo0wmC3QKH3w9cGzyCw0IjxQg/PZRcgsNKJNmD/kNQz8rOyglssAQ2kZZDI5QgN80UKnheL6Nb7MQgNsQkDt4wODyYJgP190jgrBzxdzUFJWDh+5rGK0p1IBjUqBZUk/wWQph7DDsSccrn/AW2026E0WxLUMwui72sJgtmBZ0gkcSsuC3V7RXbv/TAY6tgrBpTw9ikrKAAB2u8D+sxkos9qg06hQaDSh1FJxbU+tUuBSTjFsAhXX//SlyDOUQgjAV6lAgdEMjaoEdnsw9KXmapPgjcrtwKmMAmiUCnxUZMSR9Guw2OzOrWmJfOQyyGQVr0cuq77GymkLYQEaRAT6oW14ICYOugMBvhVdnIUlZlzILkKO3gSbTcBXpcCz99yFoXdEI8jP19EtHN8uAu9HDa+S9HxVPi6P1KztwgQ3l+f0CKKaSVp0+0bFxcXYs2cPvv/+e6xevRoGg6HRjRp1xzVCq82O/+z8GeezixDkp0aB0VzRxeivgVbtg1+u5EEhl+OOqFBAACev5MJPrUJsuA5nMwpwpcBwvTXm2oe4QgZEhwVArfRBsakMZosVZosVCrkMvkqfioSo9IG/2gcZhRWjDZUKuePaWYCvEkofOWx2AbVCjlyD2fG8clS0QP19VejVLgJatRL5xlIcSrsGuxBQKmSwXk+GkUF+yCk2wS7skMvlsNvtsNkr4lP6KGCx2mATgFohgx2oaHlW83puTDUtA7XQl1pgsVVspOrKuVD5yFFabq9SlxRyGRCoVUGnUSO72AQfhQwmsxVyOSCTVXSBCgBqpRx+KhV6xLaosph05WLiF3KKEaBRIldfig4tg/HMyB4eTTa1HZnbWJbaI6pvtfncl5wI7XY7jh8/juTkZOzatQt79uyBwWBAREQErl27JjV2j3BHIswoMGLlzp/x49kMlFqsUPvI0SU6DCH+Gvx6NQ8FBjNkMqC03AaVQoZgP1+YLFbIZDKYyqwot9lwY0+e/HpXqFxW8aFe3V/BX+0DlY8CfmolrhUZUW7/XxJQyGSICNKi0FCKsuuDPIQQTt17PrKK7r6a/sBKeUUMCjkAO2Cpp+8vOrUPbAIosVhdKF2VQg64Oq3R53oTvLIbVSEDgv180TLYHy10GuTqS6HTKHE2qxA2u6j4UmAohVqpwMjuMXi4VweE6rTVJpLGtJQdETnzWCI8ceKEY87gnj17UFRUhLCwMAwdOhQJCQlISEhAly5d3PEa3KquidBqs2NZ0k9I+ukCcgzOu2soZYC4ntRkAqh6xahuFLKKC1U2e9U/k1algN0uKgZEuPl5PU0pr+j6lEqtqHh8NacFMlR0xSp9FPBV+kCtVCDPYILNLhARqEXvuEg82r9TxZzN67smnMnIx8rkVBSVmBGoVeP3g7uiX8dWLm1lxBYXUePjsUQok8kQFBSEIUOGOBJf9+7dG/3yanVNhGcyCzD7811VkiBJp5RXtFYrW7h2UXOXp0IGyGUyR6vOV6lAwp1t0DY8EF+k/AqjuRwKuQzd2oThwd7t0SkyBGqVD3KKTNhz+ioKDGb4+yqR0LUNosICatx7z2yx3nbwChE1DR6bPnHkyBH07NkTcrn3fPO12uz4dM8vTIJupvSRw0fIKga/CAGVQo7wQA00Kh/4qVWICNLip4u5MJotAGQot9vgp/SBTqvCwE6t8dyouxzLi13O0yMkwBeh/s4r7YQFaNGxVYjLLbbaDF4houaj1tMnGkJ6ejo+/vhjZGdno1u3bnjmmWeg0Wjq5bmLTGYcOZdVL8/VHMgA+PkqUVpW7rheqVTIYbNVXOBU+yjwYO/2MJfbcPxCNopNFgRolIiPjcCDfTpAo/RxLA83oFNrlJjLcSozH4ZSC3QaNe5uF4FH+3d2tNj8fFW4I6rmuaxcQYWIbqfR9//8/PPPGDRoEO6//37069cPH3/8Mb744gvs3bsXKpXnP+AMJguMZdIGdTRXmuu7G9hsAmXl5fBRyBERqEFIoD/kkOFSnh7CLhAbHojRPWLRsXUIFDLgQq4B7SMCodP6VkzQL7XAJuxQyOROuzbcODEcqNjVwWq1w8dHXmO3JhGRVHWePuFp9913H2w2G3bs2AEAyMnJQUxMDD788ENMmzbNpTrqco0ws1CPR/+1WVLsTUGovy8CfJVQ+QClFjsidBr4aX1ht9vw69UClJjLYRMVE/Fb6LTo36kV7uvZDu3Cg1BqsQIyQKP0caxCAwCGUgsgA5MWETUYj10jrG8WiwVJSUn497//7TgWHh6O4cOHY/PmzS4nwro4eTnP48/hSUqFDK2D/VBsKodN2GAstUKhkMHXR4He7SvWwowJ08FisztGUFYuCbdm3yn8dCkHdiHQrU0LPNS7g9NuFTcOKLlx2S5XN6clImoMGnUivHz5MsrLyxETE+N0PCYmBikpKTU+rqysDGVlZY7ber1r+9tV54efL0l+bH0L9PXBmPh2GNYlCkF+vjCYrYgK8QdkcihkQJHJggBfHxSZKla/Cfb/3+onlamr8t/WIf546b67ndYXZeuOiJqjRp0IS0tLAQD+/s4j+QICAhz3Veftt9/GvHnz3BJD96gwpJxtHINlAtQydI2JQMIdUdBp1WgRqIGfSoUiUxn81Eq0DtHdcth/5aLOOq1rLTYfhRzBfmzdEVHz1qgTYWBgIACgsLDQ6XhBQYHjvurMmTMHM2fOdNzW6/WIjo6WFMOY3nFY/EOqpMfeTpgKuLNdK1hQjnKzQGyLQLQODoQVNtzVNhQadcWE71KLFQq5HC2D/KtNdG08Eh0RkXdo1IkwOjoaQUFBOHnyJMaMGeM4npqaim7dutX4OLVaDbVa7ZYYgv21mDykKz7Z84vLjwmUA3YfwF8O2H0UiAzSIjBAh64tQ9A+Ohix4UFV5rwREVHDaPSjRp9//nkkJyfj8OHDCAgIwI8//ohBgwZh+/btuOeee1yqwx1rjRYaTdh+9BwOpF+Dr1IOuVyOYI0GHSJD0DOuBQL9tBwlSUTUSNTLotv1pbCwEPfccw9ycnLQtWtXpKSk4JlnnsF7773nch3u3KGeiIgav2aVCAHAZrNh3759jpVlOnfuXKvHMxESEXmXZjOPsJJCocCQIUMaOgwiImqGeEGLiIi8GhMhERF5NSZCIiLyakyERETk1ZgIiYjIqzEREhGRV2MiJCIir9Yk5hHWVeWaAXXZjomIiJqOys97V9aM8YpEaDAYgOuLeBMRkfcwGAy33K0ITWWJtbqy2+3IzMxEQEAAZDKZpDoqt3K6cuWK1y/TxnPhjOfjf3gunPF8/E99nwshBAwGA1q1agW5/NZXAb2iRSiXyxEVFeWWunQ6nde/oSvxXDjj+fgfngtnPB//U5/n4nYtwUocLENERF6NiZCIiLwaE6GL1Go1/vKXv0Dtpp3vmzKeC2c8H//Dc+GM5+N/GvO58IrBMkRERDVhi5CIiLwaEyEREXk1JkIiIvJqXjGP0FXXrl3DlStX0K5dO4SGhnrsMU1BeXk5fvnlF6jVanTu3PmWCxEIIbBv374qxzt06ICIiAgPR1o/srOzkZaWhjvvvBNBQUEuPebq1au4du0a2rdv7/JjmgK73Y6jR49Co9HgzjvvvGVZi8WCQ4cOVTnepUsXhISEeDDK+mG1WnH27Fmo1WrExsZCoVC49LiLFy8iLy8PnTt3hr+/v8fjrC+5ubnIyMhAbGzsbd/zJpMJx44dq3K8W7duLs//cxtBwmq1iqlTpwpfX1/RpUsXoVarxZ///Ge3P6ap2L17t2jZsqWIiYkRYWFh4s477xTnz5+vsXxpaakAILp16yYGDhzo+Pn222/rNW5POH78uHjsscdEeHi4ACC2bdt228eUlZWJxx9/XGg0GnHHHXcIX19f8c9//rNe4vWksrIy8fe//120bdtWBAYGihEjRtz2MVeuXBEARHx8vNN7Y8+ePfUSs6dYrVbx5z//WbRo0UJ06dJFREdHi9jYWPH999/f8nFGo1GMGTNG+Pn5iU6dOgmtVitWrFhRb3F7yo8//igGDRokwsPDxV133SU0Go14+umnRXl5eY2PSU1NFQBEnz59nN4bx44dq9fYhRCCiVAI8f7774ugoCBx5swZIYQQKSkpwsfHR2zcuNGtj2kKDAaDaNGihZg5c6YQQojy8nIxcuRI0b9//xofU5kIU1JS6jHS+vH555+Lr776SmRkZLicCN98803RsmVLcfnyZSGEEFu3bhUymUzs3r27HiL2nIKCAjF79mxx4cIFMXny5FolwlOnTtVLjPXFYDCIefPmiaKiIiGEEHa7XcyaNUvodDpRWFhY4+OmT58u2rVrJ3Jzc4UQQnz22WdCLpeL1NTUeovdEz755BOxd+9ex+0zZ86I4OBgMX/+/BofU5kIs7Ky6inKmjERCiG6d+8unnvuOadjI0eOFA8++KBbH9MUfPnll0KhUIi8vDzHsaSkpFt+mFUmwjVr1ogjR46IgoKCeoy4fuTm5rqcCNu0aSNmz57tdKxXr15i8uTJHoywftU2EW7fvl0cPXpUFBcX10t8DeHcuXMCQI1feMrLy4VOpxPvvvuu0/G2bduKV155pZ6irD/jx48XY8aMqfH+ykS4b98+cfz4cWEwGOo1vht5/WCZymthd999t9PxPn364Pjx4257TFNx/PhxxMbGOl3v7NOnj+O+W5k+fTqmTJmCli1b4tFHH0VBQYHH421sCgoKcPny5Wb53qiLJ598Ek888QTCwsIwZcoUlJSUNHRIbnf48GHIZDK0a9eu2vvPnz8PvV5f5b3Ru3fvZvfesFqtOHHiBNq3b3/bso899hgmTZqEkJAQvPDCC7BYLPUS4428PhEWFxfDZrNVGegSGhpa4we5lMc0FQUFBVVeV0BAAJRKZY2vTS6XY9WqVcjJycHPP/+M06dP4/jx43j++efrKerGo/IcNcf3hhS+vr5Yv349srKycPLkSfz000/Ytm0b/vjHPzZ0aG6VmZmJmTNnYsqUKTUu8O9N7425c+ciOzsbL7/8co1ldDodtm3bhqtXr+LXX3/FwYMH8eWXX2LevHn1GiuYCAGlUgkAMJvNTsdLS0uhUqnc9pimQqlUVnldVqsVVqu1xtemUqnw5JNPOkaWtm3bFrNnz8bXX3/dIN/uGlJzfm9IERYWhkceecRx+4477sDLL7+MNWvWNGhc7pSXl4d7770XnTp1wqJFi2os5y3vjQ8++ADvv/8+1q5dW2PrGADatGmD0aNHO2737NkTzz77bIO8N7w+EQYGBiIoKAgZGRlOxzMyMtCmTRu3PaapiImJQWZmptOuzpW3a/PaIiIiYLVakZOT46FIG6dWrVpBqVQ2y/eGu0RERKCgoKBKQmiK8vPzMXLkSISEhGDz5s3QaDQ1lo2JiQGuvxdu1JzeG4sWLcJrr72GDRs2OCU5V0VERFQ5P/XB6xMhAIwcORLffvut47bVasWWLVswatQox7HLly/j4MGDtXpMUzRq1Cjk5eVh//79jmMbN26EVqvFwIEDgevzyPbu3Yvs7GwAqPZ6z3fffYeQkBBERkbWY/QN4/z58zh69Chw/Vv/0KFDsWnTJsf9ZrMZO3bsaPLvDVeUlZVh7969jq6+mt4bcXFx8PX1bYAI3aegoAAjR45EYGAgtm7dCj8/vyplzpw5g59//hm43jq+6667nN4bhYWF2LNnT7N4byxevBizZs3C+vXrcf/991e5v6SkBHv37kVxcbHj9s2+//77285N9YgGG6bTiPz8889Cq9WK5557TmzatEk88sgjIjw83GlY7+uvvy5CQ0Nr9Zim6tFHHxXt2rUTX331lVi6dKnw9/cXf//73x33GwwGAUD85z//EUIIsXjxYvHoo4+Kzz//XGzdulW8+OKLwsfHx3F/U5abmytSUlLE5s2bBQDx7rvvipSUFHHx4kVHmRdeeEHExcU5bu/fv1+oVCrxyiuviI0bN4oxY8aINm3a3HJYfVNx6NAhkZKSIsaMGSN69eolUlJSnIbNX7hwQQAQiYmJQggh3nrrLfHEE0+Ir776Snz77bdi6tSpwsfHR2zYsKEBX0XdmUwmER8fL1q3bi22b98uUlJSHD+VUyOEEOLxxx8Xffv2ddzeunWrUCgU4i9/+Yv45ptvxODBg0WXLl1EaWlpA70S91i9erWQyWRi9uzZTufi+PHjjjLHjx8XAERycrIQQohZs2aJZ555Rqxdu1Zs3LhRTJw4UahUKvHdd9/Ve/xcWeb6SgY//vgjFi5ciIULF6JTp044cOAAWrZs6SgTExODfv361eoxTdXnn3+ODz74AKtWrYJarcaSJUvw+9//3nG/QqHAwIEDHa/1D3/4A2JiYrB27VpkZ2ejXbt2OHToEHr27NmAr8I9fvrpJ/zlL38BAAwcOBCJiYlITEzE7373Ozz33HMAgLi4OPTq1cvxmH79+mHPnj346KOP8MEHH6Br1674z3/+0yxWl3n77bcd3d1qtRqzZ8+Gj48Pdu3aBVwfHDNw4EDHgJC5c+di/fr1+Oabb1BYWIgOHTogNTUVnTt3btDXUVfFxcXQaDSIjY3FX//6V6f75s2bhxEjRgAAOnfu7LRKypgxY/D9999j6dKl2Lt3L/r27YvXXnutybeOL1y4gAEDBiAlJQUpKSmO4x06dMCqVasAAP7+/hg4cKDjfLzzzjv48ssvsWHDBhgMBnTq1AmnTp265XVFT+E2TERE5NV4jZCIiLwaEyEREXk1JkIiIvJqTIREROTVmAiJiMirMRESEZFXYyIkIiKvxkRIbpGUlIQLFy54pO7MzEysW7fOI3V7WmFhIdasWQOr1eo45slz5QpXzmdOTg7WrFmDm6cZHz9+HOvXr8eBAwdceq6TJ09i586ddYrXHXW4S2OKhdyHiZDcYu7cudi9e3e192VmZmLbtm3YsmULzpw5U+u6jx075rSyjSfl5+djw4YNbqsvPT0dEydOdFpg+lbnqj64cj5//fVXTJw4ETabzXHslVdewf3334///ve/OHr0KFJTU/HDDz/csp7169c7VuaRyh11SFHd62uoWMizuMQaeUxubi6effZZ7NixA3379oW/vz/OnTsHrVaLpUuXOi1L1likpaVh8uTJTlsHeaPw8HA8/vjjkMv/9135P//5D7744guMHTsWuJ7Q9+7di+HDh9dYT4MsoOwm//3vf6u8vqb8eqhmTITkEWVlZRg1ahTUajXOnz+PiIgIx30///wzioqKnMqfO3cOqampCAkJQf/+/Wvcn+3kyZO4ePEi7rjjDsTFxVW5Pzs7G4cOHYK/vz/i4+Od1nnMzMzEvn378Oijj962nhsZDAYcPnwY5eXl6NWrV5WNVW927NgxZGZmokuXLrcs5yqDwYAtW7YA1/d+jIuLQ48ePVx63O3iruk8hIWF4aGHHoJMJkNeXh62bdsGg8GAY8eOoaSkBK1bt8avv/7q6EIFgAEDBlTZTqhz584IDg523D5x4gSKi4vRr18/nDhxAgUFBejduzfCwsJqdU6EEDh8+DCuXr2K2NhYxMfHV1suLS0Nv/76K2JiYnDXXXc5jqenp+Pw4cPA9Y2nu3btitjYWMf9Z8+erfb13fx6XInFXa+ZPIeJkDzis88+Q2pqKn766SenJAgA3bt3d7r94osvYtWqVRg4cCDOnz8Pu92OHTt2oH379o4ydrsdY8aMQWFhIYKCgpCcnIz3338fzz//vKPMP/7xD/z9739Hv379YDabcfLkSXz66ae47777gOsJ6re//S0+/vjjW9Zzo4MHD+K+++5Dx44dERISgl9//RXz58/HpEmTqi3/5JNP4uuvv8bAgQNx+vRp3HHHHXU6j7ie0L755hvg+pZOhw4dQufOnbF169YaF2u+Xdy3O5+VXaPjx49Hfn6+Y8uxlJQUnDp1CoMGDcKZM2eQl5fniC0mJqZKIly/fj2SkpIci1B//vnn2LJlC3x8fBAZGYmioiKcPXsW3333Hfr06ePS+SguLsb999+P9PR09OzZE4cPH0Z8fDy++eYbx36AZrMZkydPxrZt29CvXz/k5eUhKioKGzduhEwmw4ULFxxx6/V6pKSk4IUXXsCCBQuA64myute3Y8cOp9fjSizueM3kYfW+3wU1S3379hWrVq1y3B4/frzo3LnzbR+3efNmoVQqHdu1WCwWcc8994hRo0Y5ynz77beOLZAqLVq0SAQGBgq73S6EEOK7774TISEhIj093VHm888/F2FhYcJgMLhcT1pampg8ebLT63jmmWcct00mk9i2bVu1r2Xjxo1CrVaLU6dOOcr27t1bAHDEUN25qi2TySS6d+/u9Dpudqu4XTkPycnJAoAoLy8XQghRWloqAIiUlBTHY15//XUxdOjQW8b6l7/8RQwcONBx+5VXXhEKhUIcOHDAKdZx48a5XMesWbNEx44dRUFBgRBCiKysLBEZGSnmz5/vKPPqq6+K1q1bO22XtWnTJmGz2ap9jrS0NKHVasWRI0du+fqkxCLlNVP94mAZ8ohr1665tOv2mjVrMGbMGEe3lVKpxB//+Ed8//33js1dAUAmkzm2PQKAYcOGobi42LE58KpVq9ClSxccO3YM69atw9q1a4Hrg19SU1Ndrqd9+/ZYvXq1436NRoOLFy86NhPVaDQ17ry9du1ajB071rHFkEajwYsvvujyObsVIQSOHj2Kb775Bhs3bkR0dDQOHTpUY/nbxX278+BJ8fHx6Nu3r9Nz12YQ1Zo1a/Dcc885uihbtmyJqVOnOrowAeCTTz7Biy++6NgVHgDGjh3rdM3TZDJh7969WL9+PY4cOYKIiIhbnlOpsbjjNZNnsWuUPCIgIMCxb92tXLp0CXfffbfTscprVZcuXUJISAhw/drYjTuAq9Vq4HoXGABcvHgRhYWFWL9+vVNdjz32GJRKpeP27eq52V//+ldMnToVkZGR6N27N+677z784Q9/QEBAQJWyly9fdvqwA4C2bdve9hzcTmZmJkaOHAmj0Yhu3bohICAAFy5cQHh4eI2PuV3ctT0P7lT5N73xuV19XqvVioyMjCp71sXFxeHSpUvA9evTOTk56NixY431JCUlYcKECWjVqhXatWsHX19fGAwGl96ztYmlUl1eM3keEyF5RJ8+ffDuu++iqKjolhvShoWFObX8ADhu12YwgU6nQ1RUVJVv4nUVExODpKQkFBQUYNeuXXj77bexceNG/Pjjj1XKhoaGorCw0OnYzbelWLBgASIiInDy5ElHi+a5557D6dOn3RJ3U+Lj44OgoKBq3zOV7xe1Wg2NRoP8/Pwa65k5cyZeeOEFzJs3z3Gsc+fOVeZN1jUWahrYNUoeMW3aNMhkMsydO7fKfTabDZmZmQCAQYMGYfv27TCZTI77161bh9jYWERFRbn8fKNHj8aWLVtw5coVp+NZWVm1+nC7WUZGBnD9G/1vfvMb/PnPf8bhw4dht9urlB00aBC2bduGsrIyx7Gvv/7apedJTEx06sK90bVr19ChQwdHEiwtLcXWrVvdFrdU/v7+DdKqGTRokNN5FUJgw4YNGDRokOPYqFGj8Nlnnzn97fPy8hy3r127hk6dOjnu+/nnn5GWlub0PK68PldiocaPLULyiNatW2PDhg14/PHHcerUKTz88MPw8/PDuXPn8M0332Du3LmYOHEi/vCHP2DlypVISEjAlClTcObMGSxevBhr166FTCZz+fmef/55bNq0Cf369cMLL7yAsLAwnDhxAklJSTh9+nSt6rrRjBkzoFAoMHToUMjlcixZsgS/+c1vnK41VXruueewZMkSjBw5Er/73e9w7NgxbN682aXneeqppzBz5kx069atyn0PPfQQpk6diujoaISFhWHlypUwGAxui1uqXr164c0338QHH3yAiIiIaqdPeMLbb7+N/v37Y+LEiRg+fDg2b96M06dP44svvnCUeffddzF48GCMHDkS48ePR25uLtatW4effvoJMpkMDz30EP70pz+hqKgIJSUlWLhwIfz9/W/7+qTEQo0fEyF5zD333IP09HR89dVX+Omnn2C329GhQwds2rTJcR3Q19cX+/fvx/Lly7F//34EBwdj7969TsPKW7dujccee8yp7oCAADz++OOO61xqtRrfffcd1q1bh927d+PSpUuIj4/He++95/jwd6Wem61fvx4bNmzArl27YLfb8eqrr2LChAnVlvXz88OBAwfw4Ycf4tChQ+jatSt++OEHzJ8/3+k65c3OnDmD8vJyTJ06tdr7J02aBD8/P2zduhVXrlzBK6+8AplMhrNnz9ZY563iduU83DyhXqFQ4PHHH0eLFi0cjxk5ciRWr16N5ORkHDhwoNrpEzdPQO/Zs2eVbsO4uDjHJP3q3FxH165dceLECaxYsQIpKSno0aMHFi9e7NSD0KFDB6SmpmLFihU4cOAA4uLikJSU5Hg9ixYtwvLly3H48GHodDqsW7cO33//vdNzVff6pMQi5TVT/ZKJuvQbEV3Xr18/PPfcc3jyyScbOpRG7+ZzlZycjPPnz9eYCInIs9giJGpgCQkJSEhIaOgwiLwWB8uQW4waNcotUwW8Ac8VUePCrlEiIvJqbBESEZFXYyIkIiKvxkRIRERejYmQiIi8GhMhERF5NSZCIiLyakyERETk1ZgIiYjIqzEREhGRV/t/WcEkS1RB/SUAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.tl.effect_size(cells, groupby=\"Metadata_Perturbation\", reference=\"negcon\")\n", "mt.tl.wasserstein_features(cells, groupby=\"Metadata_Perturbation\", reference=\"negcon\")\n", "\n", "treated = comparison[\"KS statistic\"].drop(index=\"JCP2022_033924\").idxmax()\n", "effect = cells.uns[\"mantispy\"][\"effect\"]\n", "shift = cells.uns[\"mantispy\"][\"wasserstein\"]\n", "paired = pd.DataFrame(\n", " {\n", " \"location\": effect.loc[effect[\"group\"] == treated].set_index(\"feature\")[\"effect\"].abs(),\n", " \"shape\": shift.loc[shift[\"group\"] == treated].set_index(\"feature\")[\"distance\"],\n", " }\n", ").dropna()\n", "\n", "fig, ax = plt.subplots(figsize=(4.6, 4.4))\n", "ax.scatter(paired[\"location\"], paired[\"shape\"], s=8, alpha=0.5, color=\"#3a7ca5\")\n", "ax.set(xlabel=\"|Cohen's d|, a shift in location\", ylabel=\"Wasserstein distance, any change in shape\", title=treated)\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "id": "c011", "metadata": {}, "source": [ "Most features lie on one rising curve, where the shift in location accounts for the distance. A few sit well\n", "above it: their distributions changed shape by more than the shift of their centre explains, which a test on means\n", "alone does not report." ] }, { "cell_type": "markdown", "id": "c012", "metadata": {}, "source": [ "## Working at this scale\n", "\n", "A whole screen has a million cells by four thousand features, which is 16 GB of float32 before\n", "any transform allocates its output. `pp.downsample` caps the number of cells per group (by\n", "default per well), so the densest wells do not dominate, and `stratify` keeps rare\n", "perturbations represented.\n", "\n", "Explore on a sample, then confirm on the full data." ] }, { "cell_type": "code", "execution_count": 6, "id": "c013", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:27.326843Z", "iopub.status.busy": "2026-10-21T01:15:27.326680Z", "iopub.status.idle": "2026-10-21T01:15:27.390306Z", "shell.execute_reply": "2026-10-21T01:15:27.389275Z" } }, "outputs": [ { "data": { "text/plain": [ "{'before': (13578, 1564),\n", " 'after': (2400, 1564),\n", " 'megabytes': (84.9, 15.0),\n", " 'wells still represented': 24}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "small = mt.pp.downsample(cells, n_per_group=100)\n", "{\n", " \"before\": (cells.n_obs, cells.n_vars),\n", " \"after\": (small.n_obs, small.n_vars),\n", " \"megabytes\": (round(cells.n_obs * cells.n_vars * 4 / 1e6, 1), round(small.n_obs * small.n_vars * 4 / 1e6, 1)),\n", " \"wells still represented\": int(small.obs.groupby([\"Metadata_Plate\", \"Metadata_Well\"], observed=True).ngroups),\n", "}" ] }, { "cell_type": "markdown", "id": "c014", "metadata": {}, "source": [ "## Cell states, and what changed about them\n", "\n", "Clustering uses scanpy. `tl.cluster_composition` then turns the cluster labels into an object whose rows are\n", "wells, whose columns are cell states, and whose values are the fraction of that well's cells in each state. It is\n", "an ordinary mantispy object, so `tl.map`, `pp.normalize` and the plotting functions work on cell-state composition\n", "as they do on features." ] }, { "cell_type": "code", "execution_count": 7, "id": "c015", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:15:27.392527Z", "iopub.status.busy": "2026-10-21T01:15:27.392376Z", "iopub.status.idle": "2026-10-21T01:16:03.374006Z", "shell.execute_reply": "2026-10-21T01:16:03.372924Z" } }, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sc.pp.neighbors(cells)\n", "sc.tl.leiden(cells, key_added=\"leiden\", flavor=\"igraph\", n_iterations=2)\n", "\n", "composition = mt.tl.cluster_composition(cells)\n", "mt.pl.cluster_composition(composition);" ] }, { "cell_type": "code", "execution_count": 8, "id": "c016", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:16:03.377725Z", "iopub.status.busy": "2026-10-21T01:16:03.377567Z", "iopub.status.idle": "2026-10-21T01:16:03.391979Z", "shell.execute_reply": "2026-10-21T01:16:03.390916Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
compoundclusterfractionDMSO
0JCP2022_116560150.5480.010
1JCP2022_098688130.5500.024
2JCP2022_06743290.2840.074
3JCP2022_11596390.2630.074
4JCP2022_010404110.2390.075
\n", "
" ], "text/plain": [ " compound cluster fraction DMSO\n", "0 JCP2022_116560 15 0.548 0.010\n", "1 JCP2022_098688 13 0.550 0.024\n", "2 JCP2022_067432 9 0.284 0.074\n", "3 JCP2022_115963 9 0.263 0.074\n", "4 JCP2022_010404 11 0.239 0.075" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "by_group = pd.DataFrame(np.asarray(composition.X), columns=list(composition.var_names))\n", "by_group = by_group.groupby(composition.obs[\"Metadata_Perturbation\"].astype(str).to_numpy()).mean()\n", "control = str(cells.obs.loc[cells.obs[\"Metadata_Control\"].to_numpy(), \"Metadata_Perturbation\"].iloc[0])\n", "shift = by_group - by_group.loc[control]\n", "largest = shift.stack().abs().nlargest(5).index\n", "pd.DataFrame(\n", " [\n", " {\n", " \"compound\": compound,\n", " \"cluster\": state,\n", " \"fraction\": by_group.loc[compound, state],\n", " \"DMSO\": by_group.loc[control, state],\n", " }\n", " for compound, state in largest\n", " ]\n", ").round(3)" ] }, { "cell_type": "markdown", "id": "c017", "metadata": {}, "source": [ "Two compounds put more than half of their cells into one state: JCP2022_116560 into cluster 15 and\n", "JCP2022_098688 into cluster 13, which hold 1% and 2% of the DMSO cells. Clustering always returns clusters,\n", "most of them shared by every compound. What distinguishes a cell state from a partition of the data is whether\n", "the composition differs between conditions, which is what this object measures." ] }, { "cell_type": "markdown", "id": "c018", "metadata": {}, "source": [ "## Two different questions\n", "\n", "`cluster_composition` asks whether the mix of cell states changed. `tl.subpopulation_hits` asks whether, within a\n", "state that treated and control cells share, the treated cells differ from the controls: a shift small enough that\n", "the cells stay in their states. The statistic is a KS test on each cell's distance to the control center of its own\n", "cluster, so each cell is compared with controls in the same state." ] }, { "cell_type": "code", "execution_count": 9, "id": "c019", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:16:03.394832Z", "iopub.status.busy": "2026-10-21T01:16:03.394692Z", "iopub.status.idle": "2026-10-21T01:16:12.405974Z", "shell.execute_reply": "2026-10-21T01:16:12.404851Z" } }, "outputs": [ { "data": { "text/plain": [ "{'tests': 215, 'significant': 0, 'significant per compound': {}}" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.subpopulation_hits(cells, use_rep=\"X_pca\")\n", "table = cells.uns[\"mantispy\"][\"subpopulation_hits\"]\n", "{\n", " \"tests\": len(table),\n", " \"significant\": int((table[\"qvalue\"] < 0.05).sum()),\n", " \"significant per compound\": table.loc[table[\"qvalue\"] < 0.05, \"group\"].astype(str).value_counts().to_dict(),\n", "}" ] }, { "cell_type": "code", "execution_count": 10, "id": "c020", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:16:12.408942Z", "iopub.status.busy": "2026-10-21T01:16:12.408789Z", "iopub.status.idle": "2026-10-21T01:16:12.489972Z", "shell.execute_reply": "2026-10-21T01:16:12.488862Z" } }, "outputs": [ { "data": { "text/html": [ " \n", " \n", " " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.pl.subpopulation_hits(cells);" ] }, { "cell_type": "markdown", "id": "c021", "metadata": {}, "source": [ "None of the 215 tests is significant at q < 0.05, so on this plate no compound shows a shift within a state that the\n", "controls share. Cells are not independent draws, so read the p-values as a ranking." ] }, { "cell_type": "markdown", "id": "c022", "metadata": {}, "source": [ "## Cell cycle\n", "\n", "Much of the heterogeneity in a Cell Painting image comes from the cell cycle, so label it before attributing\n", "variation to a treatment. `tl.cell_cycle_phase` fits a two-component mixture to log integrated DNA intensity per\n", "plate. The lower component is G1, the upper is G2M, and cells that neither component claims with confidence are\n", "called S. It needs raw intensities, so it runs on the cells as loaded." ] }, { "cell_type": "code", "execution_count": 11, "id": "c023", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:16:12.492949Z", "iopub.status.busy": "2026-10-21T01:16:12.492794Z", "iopub.status.idle": "2026-10-21T01:16:14.478863Z", "shell.execute_reply": "2026-10-21T01:16:14.478132Z" } }, "outputs": [ { "data": { "text/plain": [ "{'G1': 0.549, 'S': 0.297, 'G2M': 0.155}" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.cell_cycle_phase(raw, dna_feature=\"Nuclei_Intensity_IntegratedIntensity_DNA\")\n", "raw.obs[\"Metadata_CellCyclePhase\"].value_counts(normalize=True).round(3).to_dict()" ] }, { "cell_type": "code", "execution_count": 12, "id": "c024", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:16:14.481069Z", "iopub.status.busy": "2026-10-21T01:16:14.480862Z", "iopub.status.idle": "2026-10-21T01:16:14.634366Z", "shell.execute_reply": "2026-10-21T01:16:14.633707Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEiCAYAAABz+/4kAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjIsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvgI3uAAAAAAlwSFlzAAAPYQAAD2EBqD+naQAALn9JREFUeJzt3XlYVGXfB/DvgMMgMKyKoAwgi4JCiriyhHuIkWnKoqk9iLhcGVGa8miWrxnuYqU9jyEkjwaIS4krmvq6gbtg6ouQSggoKDiDKCPL/f4RnMeRRQYGBw6/z3XNdXXOuc+Z36Hm233WW8AYYyCEEB7RUHcBhBCiahRshBDeoWAjhPAOBRshhHco2AghvEPBRgjhHQo2QgjvULARQniHgo0Qwjsd1F0AaZueP3+O8PBwAICGhga6dOmCiRMnwtbWttZykUgEBwcHTJ48GSKRiNuGXC5HTEwM7t27h6FDh8Lb27tRy8LDw/H8+XMAwIoVK6Crq8sty8/PR3x8PJ4+fQofHx+4uro2ar0aJ06cwG+//Ybly5dDLBYDAJ4+fYr//Oc/yMvLw6BBgzB27FgIBAIAwIMHD/DLL7+gsLAQjo6O8Pf3V9hHoh7UYyNNIpfLsXHjRlhYWEAikeDu3btwdXVFfn5+reUmJibYvHkzxo8fr7CNwMBAxMbGQiQSYcaMGYiNjW3UMktLS1hbW2Pjxo1cUAFAeno6RowYgZycHDx9+hTDhg3D3r17X7tejcePH2PJkiX4/vvvUVpaClQH9Ntvv41bt25BKBQiNDQUX3/9NQDg2bNncHV1xdWrV2FoaIhNmzZh+vTpKv07kyZihDRBcXExA8CeP3/OzbO1tWUnT56sc3l2djYDwJ4+fcoYY+z27dtMLBazkpISxhhjhw8fZo6Ojq9d9jIArLCwkJvOy8tjUqmUm165ciXz8/N77Xo1AgMD2cmTJ5mmpibLz89njDEml8tZXl4e1yY2NpYNHz6cMcZYZmYmMzQ0ZFVVVYwxxs6dO8dsbW2V/EuSlkCHoqRZFixYAA0NDWRkZMDFxQVubm51tistLUWHDh24Q7i0tDS4urpCT08PAODl5YVbt25BLpc3uKyhwzxzc3OF6Tt37qBnz56N2o/du3ejU6dO8PLyUpivpaUFc3NzLFmyBA8fPsSNGzewdu1aAICdnR0+/vhjvPPOO7C1tUV6ejo2bNjQqO8jLYuCjTSLlZUVOnToAMYYjh8/jocPH8LCwoJbvmDBAlRWViI5ORlffvkldHR0AADFxcUwMDDg2mlra0NLSwvFxcUNLjMzM2tUXYmJiUhJScHZs2df27awsBBr167FsWPH6m0jkUigqamJkydP4ubNm3Bzc0N5eTkuXLgAAwMDSCQSXLp0Cenp6fD19W1UjaTlULCRZvn444+hra0NAJg5cyY2bdqEiIgIbrmVlRVevHgBLS0tWFpacvM7deqEx48fc9OlpaUoLy+HiYlJg8saIz4+HitXrsSxY8e4CwANWblyJbS1tbF48WIAQFVVFZYuXYrw8HB0794dADBr1iwAwJgxYzBu3DgEBwfj4MGDKCwsxJEjRwAAM2bMgIWFBT799NM6L0yQN4eCjajMs2fPUFFRoTCvJvgmT56MAQMGYMSIEZBIJOjfvz+uXLmCgoICmJqaYv/+/XB1dYVQKGxw2evExMQgMjISR48ehampaaPqHjNmDCQSCTctEAggkUggEolw584dGBoawtjYGKjuadaEloaGBsrKylBVVQUNDQ2UlpaCMYYOHehnpXbqPslH2qaaiwMff/wxCw0NZePGjWOGhobsxo0bCstfvrjw5ZdfssDAQG565syZzM7Ojk2bNo0ZGhqyffv2NWrZhg0bWGhoKAPAZs2axb799lvGGGOnT59mAoGATZ48mYWGhrLQ0FD2/fffv3a9V7188eDmzZusT58+7MMPP2SBgYHMyMiIxcbGMsYYe/bsGXNycmJDhgxhISEhTCKRsHnz5qnwr0yaSsDoDbqkCeRyOX788UeguofTuXNnjB49Gp06dVJYPm/ePGhqagLV94NFR0dj1qxZEIlEYIzh8OHDuHfvHjw9PeHk5MRtv6FlCQkJ3G0lAGBsbIxp06YhKysL+/fvV6jTwsICEydObHC9V3333XcIDg5WOB946NAhVFVVwdPTE1ZWVgp/h0OHDqGgoAC9evWCh4eHCv66pLko2AghvEM36BJCeIeCjRDCOxRshBDeoWAjhPAOBRshhHco2AghvMPrW6SrqqqQl5cHsVjMPXxNCGmbGGMoKSlB165doaHRcJ+M18GWl5en8KgMIaTty8nJUXjRQl14HWw1D0Dn5ORAX19f3eUQQppBJpNBIpE06sUGvA62msNPfX19CjZCeKIxp5V4HWyEtHeMMcjlcnWX0SRCoZB7zlhZFGyE8FheXh5kMpm6y2gSgUAACwsL7k3KyqBgI4SnKioqIJPJYGJi0uZOxTDGUFhYiPv378Pe3l7pnhsFGyE8VfPST7FYzL3luC3p3Lkz7t27h/LycqWDjW7QJYTn2uo9nM2pm3pshLQT4XuuN3sbEROc611WVVWFrVu34vjx4zA1NcWCBQtgYWGBdevW4ejRoxg/fjw3dkRLo2AjhKhEeHg4bty4gWnTpqGgoAAzZszAkSNHMHr0aGhqauLWrVtvrBYKNsI7r/ZMGuplENXZt28f0tLSIBAI4Ovri9LSUpSVlcHZ2RmlpaW4d+/eG6uFzrERQpqtsrISIpEIWlpa0NTUxKefforS0tJao5a9KRRshJBm09TUhFAoRFpaGjQ0NODt7d3ocWBbAh2KEkJUYv369fDx8YGVlRXKy8thbW2Njh07Yvv27di0aRPy8/Mhk8kQHR3d4rVQsBHSTrT0uUZPT0/cvn0b169fh4mJCezt7QEAgwcP5oZlrBlsuqVRsBFCVEZXVxeDBw9WmGdnZwc7O7s3WgedYyOE8A4FGyGEdyjYCCG8Q+fYCO+9fMMu3azbPlCPjRDCO0oH28WLFxEUFIS+ffvi4MGDtZZHRUWhb9++Ch9PT89a7S5cuAA/Pz8MGjQIH374YZ3PkTWmDSGkdXjw4AEiIiIQFhaGyMhItdai1KHo999/j23btmH27NmIiYlBUVFRrTYPHjxAVVUVYmNjuXmvvkspLS0NXl5emDNnDubMmYNt27bB3d0dV69ehZWVVaPbEEIaLykpqdnb8PX1rXfZ+PHj4ebmBisrK5ibmzf7u5pDqWALCgrCvHnzAAAzZ86st52Ojg769u1b7/JvvvkGgwYNwvr16wEAXl5ecHBwwLp16/Ddd981ug0hpPUoKSlBSEgIevbsqe5SlDsUbexdw5mZmfD09MSoUaOwePFiPHnyRGH58ePH4ePj898iNDTg4+OD33//Xak2hJDW49///jcWLVoEZ2dnrFq1Sq21qPyqqJaWFoKDg+Ht7Y0nT57g22+/RVxcHNLS0iAWi1FaWoqioqJaXdWuXbsiJycHABrVpi5yuVxhRJ62OogFIW2Ru7s73N3dUVVVhSFDhiAoKAidO3dWSy0qD7awsDAIhUJueujQobCxscHmzZuxcOFClJeXAwBEIpHCeiKRiFvWmDZ1iYiIwLJly1S6P4SQxlm4cCHkcjkKCwtRWFgIQ0NDtdWi8mB7OdQAwMjICH369MH163/fSyQWi6GlpYXHjx8rtHv8+DH3mpPGtKlLeHg4PvvsM266ZuRoQkjDJ/5VoeatHv369cPGjRtrZcGb9EZu0M3Ly4OjoyNQfYW0T58+uHDhAubMmcO1SUlJgaura6Pb1EUkEtXq5RFC3oy5c+equwSOym/QXbZsGYqLi4HqwR2++eYbZGVlYfLkyVybkJAQ7Nq1C9euXQOqLxScOHFC4UprY9oQQkhdlOqxnT9/XmGUmcWLF2Pt2rWYMGECli5dCgAwMzODk5MTtLW1UVxcDH19fSQmJircpDtjxgzcunULgwYNgpmZGQoKCrB8+XK8++67SrUhhJC6CBhjrLGNnz59iqysrFrzTUxMFM5lMcaQnZ0NHR0dmJqa1ru9J0+eID8/HxKJpN5h7BvTpj4ymQwGBgaQSqVtbiRs0nQNDTPXnp4VLSsrw927d9G9e/c2OWDyq/Ur83tWqsemp6fX4I23NQQCAaytrV/bztDQ8LVXThrThhBCXkYPwRNCeIdeW0QIUZnc3FykpqbC1NQUHh4eEAgESEhIwNOnT9G9e3cMGzYMAoEAR44cgbGxMQYMGAAAOH36NBhjePvtt1VSBwUbIe1E/tKvmr0N8/+p/wb4/fv3Y/bs2XBzc0NBQQG6deuGHTt24OrVq3j06BG+++47jB07Ft9++y02bNgATU1NHDhwAAAwf/58ODo6UrARQlqXxYsX4+zZs7CwsEBMTAzOnj2LiooKrFy5EqjuzX388cdce0tLS9y5cwcFBQXo168fnj9/rrJa6BwbIaTZKisrIRAIYGVlBcYYUlNTkZycjLKyMq7Nxo0b8dVX/+01Tp48Gb/88gvi4+MRGBio0noo2AghzaapqQm5XI7S0lJ06NABUVFR3NNG5eXlmDdvHsaNG6dwV4WLiwuuXr2KgoICdO3aVaX10KEoIUQlPvroIwwbNgwBAQEoLy9HYWEhhEIhPvjgA5ibm+PWrVt49OgRxo0bx60TEBAAAwMDlddCwUYIUYmFCxfCyckJJ06cgImJCY4cOQKRSAR7e3tIpVKkpqbC2toa48aNg7e3N4RCISZNmgQAKCwshJeXl8pqUerJg7aGnjxon+jJg7+15ycP6BwbIYR3KNgIIbxDwUYIz7XVs03NqZsuHhDCUx06/P3zLikpgUAgUHc5SmGMobCwEAKBoElv4qVgI4SnOnToAH19fTx+/LjWa/bbAoFAAAsLi1rjEjcGBRvhhYauhLZnXbt2RadOndrk4ahQKGxSqIGCjRB+EwgE7XIcELp4QAjhHQo2QgjvNOlQVCaT4fbt27CxsYGxsXGdbe7fv4/y8nJYWlrWOk7Oz89Hbm6uYiEdOtT52vGioiLk5ubCysqKnh4ghDSKUj22P//8E7NmzULPnj0xYMAAHDx4sFabmJgY2NnZYfDgwRg6dCi6deuG+Ph4hTZbt27F0KFDMXv2bO4TFham0IYxhtDQUJibm+P999+Hqakpli9f3tT9JIS0I0r12K5evQoXFxesW7cOYrG4zjZZWVk4evQounfvDgD44YcfMHXqVDg7O6N3795cOycnJ6Smptb7XT/99BNiYmJw6dIlODs743//938xcuRI9OnTB++9954yZRNC2hmlemwTJ07E7NmzGxwGb8WKFVyooXp0aA0NDZw7d06hXVVVFTIyMpCdnV3npeiffvoJEydOhLPz3w8te3l5Yfjw4YiKilKmZEJIO9TiFw9u3ryJFy9e1BqO7+LFi/Dx8UH//v1hbm6OnTt3cssqKyuRlpbGDfRQY9CgQbhy5UpLl0wIaeNa9D62srIyzJgxA4MGDcKIESO4+S4uLsjIyECPHj3AGMOqVaswefJk2NjYoH///igpKUF5eTlMTEwUttepU6cG76CWy+WQy+XctEwma6E9I4S0Zi3WYysvL4efnx8ePXqE3bt3Q0Pjv181duxY9OjRA6i+gXDRokWwtbXlLjLUPBv2ckgBwPPnzxt8biwiIgIGBgbc5+XR6Qkh7UeLBFt5eTn8/f3xxx9/4MSJE+jWrdtr1zE3N8f9+/cBALq6ujA2NkZeXp5Cm7y8PFhaWta7jfDwcEilUu6Tk5Ojgr0hhLQ1Kg+2iooKBAQE4Nq1azh58mSdQfTyyDWovlctLS2NG/wBAEaMGMGNOYjqiw0HDx5UOKR9lUgkgr6+vsKHENL+KHWOraSkBBkZGdz03bt3cenSJXTu3BlWVlYAgClTpiA5ORnbtm1DQUEBCgoKgOqHcWtGohk6dCj8/Pzg4uKCR48eYeXKlTAwMMDcuXO5bS9ZsgSDBw/GJ598Al9fX8TGxqKoqAiff/65qvadEMJTSo15cOHCBYXwqfH+++9jyZIlAAB3d/da58YAICQkBCEhIUD1wA2RkZE4f/48dHR0MHDgQISGhta6N+7y5ctYu3Yt/vrrL9jb2yM8PBw9e/Zs9M7RmAftR2Pf7tGexjzgG2V+zzSYC+EFCjb+o8FcCCHtGgUbIYR3KNgIIbxDwUYI4R0KNkII71CwEUJ4h4KNEMI7FGyEEN6hYCOE8A4FGyGEdyjYCCG8Q8FGCOEdCjZCCO9QsBFCeIeCjRDCOxRshBDeoWAjhPAOBRshhHco2AghvKP0SPAymQzbt2/HlStXEBwcjMGDB9dq8+jRI0RHRyM7Oxv29vaYMWNGrYFaVNWGEEJepVSPbdeuXXBwcEB6ejq2bt2KrKysWm3y8vLQt29fHDlyBNbW1oiLi8PAgQMhk8lU3oYQQuqiVLC5uLjg9u3b+Ne//lVvm+XLl8PQ0BCHDx/GggULcOzYMRQXFyMyMlLlbQghpC5KBZutrS309PQabJOUlIQPPvgAQqEQACAWi+Hr64ukpCSVtyGEkLqo9OJBWVkZcnNzYW1trTC/e/fu3GGrqtrURS6XQyaTKXwIIe2PSoPt+fPnQHXv6mVisZhbpqo2dYmIiICBgQH3kUgkKtkvQkjbotJg09PTg4aGBoqLixXmFxUVcSM3q6pNXcLDwyGVSrlPTk6OCveO8EH4nuvch/CX0rd7NEQoFKJHjx64efOmwvwbN27AyclJpW3qIhKJIBKJVLhHhJC2SOU36AYGBiIhIQGFhYUAgD///BMHDx5EYGCgytsQQkhdBIwx1tjGt27dwrp16wAAW7duxdChQ2FrawtPT09Mnz4dqD4/5uPjg8zMTAwaNAinT5+Gp6cndu7cCU1NTZW2eR2ZTAYDAwNIpdIGD2FJ29eUQ8uICc4tUgtpGcr8npUKtry8PBw8eLDWfAcHB3h4eHDTVVVVOH36NP766y/Y29vX+XSCqto0hIKt/aBg478WC7a2hoKt/aBg4z9lfs/0EDwhhHco2AghvEPBRgjhHQo2QgjvULARQniHgo0QwjsUbIQQ3qFgI4TwDgUbIYR3KNgIIbxDwUYI4R0KNkII71CwEUJ4h4KNEMI7FGyEEN6hYCOE8A4FGyGEdyjYCCG8o9Lh9wCgsrISdb1tXENDAxoaf+coYwyVlZUKywUCQaMHaSGEkIaovMfWp08faGtrK3yEQiEmTpzItVmxYgWEQqFCmy5dutTaVkxMDGxsbNChQwc4ODjg119/VXW5hBAeUnmw/fHHH6ioqOA+N27cAAD4+fkptBs0aJBCu0ePHiksP3jwIEJCQvDNN9+guLgYc+fOxaRJk3Dx4kVVl0wI4ZkWP8e2detWmJiYYPz48Uqtt27dOrz33nuYPHkyxGIxPvnkE/Tr1w8bN25ssVoJIfzQosFWUVGB2NhYTJs2DSKRSGHZtWvX0LFjRxgZGeGdd97BtWvXuGWMMaSmpsLLy0thneHDh+PcuXMtWTIhhAdaNNj279+Phw8fIjg4WGF+165d8Z///AcFBQVIT0+Hqakp3n77bWRnZwMASkpK8OzZM3Tu3FlhPVNTUzx8+LDe75PL5ZDJZAofQkj706LBFhUVBXd3d/Tq1UthflBQECZNmgSxWAyJRILo6Gjo6+sjKiqqwe1VVVVBIBDUuzwiIgIGBgbcRyKRqGxfCCFtR4sFW15eHg4fPoyZM2e+tq1QKESPHj2QlZUFABCLxdDV1UVBQYFCu8LCQpiZmdW7nfDwcEilUu6Tk5Ojgj0hhLQ1LRZsP//8M/T09DBp0qTXtn3x4gX+7//+D926dQOq72lzc3PDiRMnFNr9/vvvcHNzq3c7IpEI+vr6Ch9CSPvTIsHGGEN0dDSmTJkCHR2dWssnTJiA48ePo7i4GJmZmZg6dSpkMhlCQkK4NvPnz8eBAwcQHR2NwsJCrFmzBmlpaQgLC2uJkgkhPNIiwXb69Gncu3ev3sPQhQsXYt26dbC3t8eIESNQWVmJCxcuoEePHlyb0aNHY9u2bVizZg2sra2xY8cO/Prrr3BxcWmJkgkhPCJgdT3/xBMymQwGBgaQSqV0WMpz4XuuK71OxATnFqmFtAxlfs/0EDwhhHco2AghvEPBRgjhHQo2QgjvULARQniHgo0QwjsUbIQQ3lH5q8EJeVOacu8aaR+ox0YI4R0KNkII71CwEUJ4h4KNEMI7FGyEEN6hYCOE8A4FGyGEdyjYCCG8Qzfoknbr5Rt8W+qlk/lLvwIAmP/PshbZPqkb9dgIIbxDwUYI4R2VH4peunQJqampCvNEIlGtgV0qKyvx+++/Izs7G/b29vDy8qo1GHJj2hBCyKtUHmyHDx/GDz/8gIkTJ3LzXh2Cr7S0FKNHj0ZeXh7c3Nzw1VdfoV+/fti7dy+EQmGj2xBCSF1a5OKBtbU1fvjhh3qXr1q1CtnZ2UhPT4exsTGys7Ph5OSEqKgozJkzp9FtCCGkLi1yjk0qlWLbtm1ISEhAVlZWreUJCQnw9/eHsbExAMDKygrvvvsuEhISlGpDCCF1aZFgk8lkSE5ORmxsLHr16oUvvviCW1ZeXo7MzEw4OjoqrOPo6IibN282uk1d5HI5ZDKZwocQ0v6oPNjeffdd3LlzBzt27MCBAwdw8OBBrF27Fr/99htQfe6MMQZDQ0OF9YyMjLggakybukRERMDAwID7SCQSVe8eIU1Scz8beTNUHmx9+/aFSCTipkeOHAkXFxccOXIEAKCtrQ0AKCkpUVhPJpNxFxka06Yu4eHhkEql3CcnJ0eFe0YIaSveyJMHIpEIT548AapDSyKR4O7duwpt7ty5A3t7+0a3qe97Xg5VQkj7pPIe2+3btxWmb926hcuXL8PDw4Ob99577yExMRFyuRyovtiQlJSE9957T6k2hBBSFwFjjKlyg8OGDUPnzp3h4uKCR48eITo6Gm5ubtizZw/Xm3r48CEGDRoEiUSC0aNHY+/evSgvL8e5c+cgFosb3eZ1ZDIZDAwMIJVKoa+vr8rdJK2AKgdzaelnRUHPizabMr9nlffYjh8/jqlTp6K0tBQmJiZITEzEgQMHFA4Ru3TpgqtXr2LixIl48uQJQkJCcP78eYXAakwbQgipi8rPsQkEAvj6+sLX17fBdkZGRggNDW12G0IIeRU9BE8I4R0KNkII71CwEUJ4h4KNEMI7FGyEEN6hYCOE8A4FGyGEdyjYCCG8Q8FGCOEdCjZCCO9QsBHyhuQv/YpeOPmGULARQniHgo0QwjsUbIQQ3qFgI4TwDgUbIYR3KNgIIbxDwUYI4R0KNkII76g82IqKirB8+XKMHDkSXl5emD9/Ph4+fKjQZsuWLXByclL4uLm51dpWSkoKJkyYAFdXVwQEBOCPP/5QdbmEEB5SebCNHj0aFRUVWLRoEb7++mtcuXIFbm5u3IDJAFBQUABNTU3Ex8dzn61btyps5+rVqxg2bBjs7OwQGRkJPT09eHh44N69e6oumRDCMyofpers2bMKQ+25uLjAxMQEhw8fRkBAADe/Y8eOcHJyqnc733zzDYYMGYLVq1cDADw8PHD69GmsXbsWP/zwg6rLJoTwiMp7bC+HGgBoampCIBCgqqpKYX5mZiaGDBmCYcOGYeHChSgqKlJYfuLECfj4+HDTAoEAPj4+OHHihKpLJoTwjMp7bK9avnw5dHV1MWrUKG6eSCTC7Nmz4e3tjSdPnmDFihVISEhAeno69PX1UVpaiuLiYpibmytsy9zcHDk5OfV+l1wuh1wu56ZlMlkL7RUhpDVr0WCLjY3F+vXrkZiYiM6dO3Pzw8LC0KHDf7/a09MTNjY22Lx5MxYtWoSKigoAgJaWlsL2RCIRysvL6/2+iIgILFu2rEX2hRDSdrTY7R7x8fEIDg5GdHQ0xo8fr7Ds5VADAENDQ7z11lvcVU+xWAwtLS08fvxYod3jx4/RqVOner8zPDwcUqmU+zTUuyOE8FeL9Nh27tyJ6dOn46effsK0adMatU5ubi569eoFANDQ0EDfvn1x/vx5zJkzh2tz7tw5uLq61rsNkUhU6xwf4ZfwPdfVXQJpA1TeY9u1axemTp2KLVu2YPr06XW2Wbp0KR49egQAqKysxLJly3Dnzh1MmTKFazNr1izs2rULly9fBgAkJyfj5MmTCAkJUXXJhBCeUXmP7aOPPoKmpibWrFmDNWvWcPPnzp2LuXPnAgAsLS3h4uICDQ0NFBcXw9TUFLt374a7uzvXPigoCBkZGXB3d4eJiQmKi4sRERGhcKWUkLao5i265v9D54NbisqD7cKFC7Vu7QAAU1NT7p+Dg4MRHByM+/fvQ0dHB8bGxnVua9WqVViyZAkePHiAbt26QUdHR9XlEkJ4SOXBVnOerDEsLCxe20YsFkMsFjezKkJIe0IPwRNCeIeCjRDCOy3+5AEh7VFSUhIAoL+6C2mnqMdGCOEdCjZCVKymt0bUh4KNEBVSJtRoZPiWQ8FGCOEdCjZCCO/QVVFCVIDOq7Uu1GMjpJko1FofCjZCCO9QsBFSjXpe/EHBRkgDKOzaJrp4QFq9N/HW3IYCLCkpCb6+vrWmKfRaLwo2Ql7ycli9HGYNtSOtDx2KElKPl8Or5p8p0NoGCjZCGolCre2gQ1HS6tBIVKS5Wn2wPXjwADk5ObCxsYGJiYm6yyEq9HKARUxwVmst6kSDu6heqw22qqoqzJkzB9u2bYOtrS3+/PNPfP7551ixYoW6SyPNUF9vjHppRJVabbD9+OOPSEhIwLVr1+Dg4ICUlBR4eXnB1dUVEyZMUHd5vNfc3hQFlfKo56Y6AsYYU3cRdXF1dUW/fv3w008/cfPGjBkDTU1N7N+/v1HbkMlkMDAwgFQqhb6+fgtW2/a19yByE95rke32v3ipSeu9Gm5NCb2X1+FDaCrze26VPbbKykpcv34dM2fOVJg/cOBAREVF1bueXC6HXC7npqVSKVD9ByF/+3rfDXWX0Co9Ez5rke2WvPTfo1LrLVyk1HwAMFv8TwDAgxXf1ruObhv+LdT8jhvVF2OtUHFxMQPAdu7cqTB/48aNTFtbu971vvrqKwaAPvShD48/OTk5r82QVtljEwqFAICysjKF+c+fP4eWlla964WHh+Ozzz7jpquqqlBUVAQTExMIBAKlapDJZJBIJMjJyWnzh7F82Re+7AdoX5qEMYaSkhJ07dr1tW1bZbDp6urCxMQEubm5CvNzc3NhaWlZ73oikQgikUhhnqGhYbNq0dfXb/P/4dXgy77wZT9A+6I0AwODRrVrtU8ejBo1SuFO78rKShw4cACjRo1Sa12EkNav1Qbbl19+ibS0NMyaNQtJSUkIDAxESUkJPv/8c3WXRghp5VptsPXq1QspKSl48eIFIiMjYWxsjNTUVHTr1u2NfL9IJMJXX31V69C2LeLLvvBlP0D70uJa7X1shBDSVK22x0YIIU1FwUYI4R0KNkII77TK+9jUraioCNnZ2bCysoKxsbG6y2mWBw8eID8/H9bW1jAyMlJ3Oc324MEDZGVlwc7ODmZmZuouRykVFRVITU2tNd/BwQGdOnVSS02qkJmZiRcvXsDR0REaGq2kr6SKR6D44tq1a2z06NHM2NiY9e3bl+no6LCAgABWWlqq7tKUdv78eTZkyBBmYWHB+vTpw7S1tdm0adOYXC5Xd2lN9uzZM+bk5MQEAgH78ccf1V2O0goLCxkA1rdvX+bu7s59jh07pu7SmuTixYusV69ezMzMjPXr14+99dZb7MaNG+ouizHGGAXbSxITE1lycjI3nZOTw7p168bCwsLUWldTJCQksLS0NG46KyuLicVitmHDBrXW1RwzZ85koaGhTCQStelgu3r1qrpLabacnBxmaGjIPv30U1ZZWckYYywjI4MdPXpU3aUx1lqfFVWXiRMnKkxbWFjA19cXZ86cUVtNTeXn56cwbWtrC4lEUusxtbYiMTERKSkpuHjxIv71r3+pu5xmycnJQVVVFWxtbRv9iFBrs379eujq6mL16tXc4WePHj3Qo0cPdZcG0Dm2hjHGcPnyZdjZ2am7lCZ58eIFLly4gOfPn+PAgQN4+vQpZs+ere6ylHbv3j3MmzcPycnJ0NbWVnc5zRYcHAxTU1NkZGTAz88PmzdvbnPPi/7+++/w9vYGAFy5cgXGxsawtLRsNefYKNgasG7dOqSlpWHLli3qLqVJpFIpFi1aBJlMhj///BNffPEFrK2t1V2WUioqKhAYGIhFixbhrbfeUnc5zaKlpYW4uDgEBAQAAG7fvo0RI0YgLCwMW7duVXd5SsnLy4NUKoWjoyP09PSQl5cHMzMzxMXFoXfv3uoujy4e1Cc2NpYJhUIWFxen7lJUIjMzk3Xp0oUtWbJE3aUoZdWqVczOzo6dOnWKnT59mp0+fZppaWmx+fPns0uXLqm7vGbbsGED09XVZVVVVeouRSnm5uasY8eO3HncsrIyNmbMGNa3b191l8YYnWOr244dOzBjxgzExMRw/3dt6+zs7DB+/HgcOnQIy5cvV3c5jSYUCtGlSxeEh4dz88rLy7F3717k5OQgPj5erfU1V5cuXVBaWgqpVNrsV2y9SdbW1ujYsSPXixaJRPjHP/4BPz8/lJSUQCwWq7U+CrZXxMXFISgoCNHR0ZgyZYq6y2my0tJS6OrqKszLyspqc0MYhoWFISwsTGGetrY25s+f3+bOF9b17yQ5ORkWFhZtKtQA4J133kF8fDwYY9xLXO/fvw8dHZ1a+6gOFGwv2bdvH6ZOnYqQkBBYW1tzV0NFIhEGDBig7vKUMnbsWIwYMQL9+vVDRUUFdu3ahVOnTiE5OVndpbVbmzZtQnp6OsaOHQt9fX0kJSVh+/btiI2NVXdpSvvkk08QExODoKAgBAQE4M6dO1i+fDkWLFjQKi4g0Ns9XrJhwwbs3r271nxTU1Ps2bNHLTU1VUlJCTZv3oyzZ89CIBDAwcEBc+bMaXMXD+oyYsQIfPLJJxg3bpy6S1Ha3r17sWfPHjx69Ah2dnaYNWsWnJyc1F1Wk+Tn52P16tW4fv06TE1NMWHChFq3TKkLBRshhHfU32ckhBAVo2AjhPAOBRshhHco2AghvEPBRgjhHQo2QgjvULARQniHgo2n7t+/3+ZuKm6utrTPbanWtoiCjadSU1MRFBTU4t/z8OFDxMfHIz4+Hjt37sTBgwdx69YtVFVV1du2rhd3JiUlISMjo9b8U6dOYefOnWjMfeRN2ee//voLe/fuVWodVXi1VnXVwVcUbKRZrl+/jsDAQMTFxWHPnj3YvHkzRo8eDSsrK/z73/+us623tzcePHigsGzevHlISkpSmFdWVoZx48bB398fx48ff20tEokEH3zwgVL1nzt3DjNnzlRqHVV4tVZ11cFX9BB8O3Pz5k1kZGTA1NQUgwcPhqampsJyxhhSU1NRXFwMZ2dnCIVCnD59GpMmTWpwu5s2bYKFhQW3je3btyMoKAhlZWUIDQ1VaNulSxcsW7YMP/74Y4Pb3LVrF3R0dODj44OoqCiMGDGiwfbdunXD2LFjuem7d+/i+vXrePfdd5Geno779+/jrbfegqWlJVDdg0xJScGLFy+41x/17t0bzs7OQPXh4pUrV2BoaIh+/fpBT0+v0duuUVxcjEuXLkEgEKB///7cWzxerrWuOqytrXHv3j2MGTNG4fXhUqkUhw4dgo+PT5t76+4bpe4XwpGWkZiYyAwMDLjpyspKNmXKFKavr8+8vb2ZpaUlc3JyYrm5uVybFy9eMG9vb2ZkZMTGjBnDJBIJe+edd5hIJKr3e44ePcoAsJycnFrLwsLCmIGBATcyVk3bHTt2sA4dOrCMjAyurZWVFVuzZo3C+l5eXmzJkiXszJkzTCQSscePHyu1zzExMczIyIh5enoyDw8PNnz4cCYSiVh8fDxjjLHr16+zIUOGMC0tLebv78/8/f1ZQkICY4yxRYsWMUNDQzZmzBjm5ubGzMzM2KlTpxq9bcYYO3LkCDMwMGCenp7M29ub2djYsP3799eqta464uLimKWlJdu4caPCPkZGRrKuXbuyioqKBv8W7R0FG0+9+iP/+eefmZ6eHsvMzGSseii7gQMHssDAQK7Nli1bmJGREfvrr78YY4xJpVLm6OjY5GA7ceIEA8AuXLig0LawsJCNHDmSffDBB1zbV4MtMzOTaWpqsrt37zLGGOvdu3etH/nr9jkmJoYBYD///DM37+uvv2Y2NjbcdFxcHDMxMVHYzi+//MIsLCxYfn4+Ny8yMpJZWVlxgdKYbQ8dOpT985//5KalUik31N6rtdZVx9dff8369OmjMK9Pnz4sPDy8wb8DYYzOsbUT8fHx8PPz4wam6dixIz777DPs2rULlZWVQPWhn7+/PyQSCQBAX18fs2bNavJ31gxoXFhYWGvZqlWrsHfv3joHEAaArVu3YuTIkdxrlmbOnNmkcQF0dXUxbdo0bnro0KG4e/cuXrx4Ue86MTExcHZ2xpkzZ5CYmIidO3dCW1sb2dnZuHPnTqO33bFjR2RlZaG0tBSo/nu+7nD6ZUFBQbh+/TquXLkCALh8+TLS0tLwj3/8Q8m/QvtDwdZOZGdnw8bGRmGera0tysvLkZeXB1QPC/fq+9qa8/62kpISoDoAXtWvXz/4+fnhiy++qLWssrIS27Ztg6WlJXfFVUtLC+np6bh48aJSNRgaGnJveEX1S0MZYw0G271795Cbm4tdu3Zh9+7d2LNnD06cOAF/f3+Fq72v2/a6detw//59mJqaYuTIkYiMjMTz588bXbtEIsHo0aMRHR0NAIiOjoanpyfs7e2V+hu0R3TxoJ3o1KkTioqKFOYVFRVBIBBwrws3NjbGkydPFNoUFxc3+TtTUlKgoaFR7+hSK1asgIODQ62roQcOHEBZWRlkMhl+/fVXbr6TkxOioqJa/G3G+vr6GDhwIDZv3tys7Tg6OuLs2bMoKCjA8ePHsXz5chw/fhz79u1r9DZmzpyJ4OBgrFixAnFxcVi/fn2zamovqMfWTnh4eGDfvn2oqKjg5iUmJsLFxQU6OjoAAHd3dyQlJSn0Sn777bcmfV9ubi7WrFmDgIAAGBkZ1dnGxsYGs2bNwqJFi7jDYQCIioqCn58f11ur+axevRrx8fHcoZ0q6OnpoaysTGGet7c3du7cWet/BMoONl3T3tTUFAEBAZg/f369h9511QEAvr6+0NLSwkcffYSKiorXXp0mf6MeWzvxxRdfYMeOHRg1ahQCAwNx+fJlxMbGKoyB8PnnnyMmJgZjx47F+PHjcfbsWZw5c0bhcKs++/btg7GxMUpLS5GWlobt27ejf//+r72lY+nSpdi2bRt32Jqfn49Dhw7V6sWh+pXgjDEkJibio48+atLf4VV9+vTBixcvsGTJEjg5OaF3795YuHAhDh8+jAEDBmD27NnQ19fHxYsXce3aNVy6dKnR2/7www9hYWEBd3d3lJeXY+PGjfUGU1111NxuM336dKxevRrBwcGtYqCUtoB6bDz16g2ghoaGuHLlCkaNGoUzZ85ALBbj0qVL8PLy4tqYmpri0qVLcHZ2xsWLFzFgwACsXr26waHUzMzM4O/vj1OnTmHfvn1ITU2FmZkZDh06hOTkZIV7rWraikQibl7nzp0RGRkJf39/ODg44ObNm/Dz88Pw4cNrfZeWlhYWL15c58WIuva5e/futcZF6NSpE/z9/SEUCrl1jhw5gidPnuC3337DrVu3IBaLcfbsWXz55ZfIzMxEWloaPDw8cO7cOaW2fezYMXh7eyMtLQ0ZGRn49ttv8d1339VZa1111Hj//feB6osJpHFozAOioKioCMbGxtz0hx9+iMePH+PQoUNqras9W7p0Kfbs2YM//vhD3aW0GXQoShT4+vpi5MiRkEgk3Inuw4cPq7usdik9PR1nzpxBZGQktmzZou5y2hQ6FCUK9u7dC21tbaSkpKBnz55IT0+Hh4eHustql27cuIGUlBRs2LABAQEB6i6nTaFDUUII71CPjRDCOxRshBDeoWAjhPAOBRshhHco2AghvEPBRgjhHQo2QgjvULARQniHgo0Qwjv/D0xW+O3Oq5r/AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mt.pl.cell_cycle(raw, dna_feature=\"Nuclei_Intensity_IntegratedIntensity_DNA\");" ] }, { "cell_type": "markdown", "id": "c025", "metadata": {}, "source": [ "Two peaks, G1 and G2M, with S between them: 55% of the cells are called G1, 30% S and 16% G2M. Check the plot before\n", "using the labels: two separate peaks with the colours split between them means the fit worked, and a single smeared\n", "distribution means it did not, whatever the counts say. This is a heuristic, and it cannot tell G0 from G1." ] }, { "cell_type": "markdown", "id": "c026", "metadata": {}, "source": [ "## Crowding\n", "\n", "A cell in a dense field looks different from an isolated one regardless of treatment.\n", "`tl.neighbors_local_density` measures crowding as the mean distance to the *k* nearest cells in the same field of\n", "view, never across fields, from the centroids `read_profiles` keeps on `obs`. `pp.regress_out` then removes it." ] }, { "cell_type": "code", "execution_count": 13, "id": "c027", "metadata": { "execution": { "iopub.execute_input": "2026-10-21T01:16:14.635996Z", "iopub.status.busy": "2026-10-21T01:16:14.635813Z", "iopub.status.idle": "2026-10-21T01:16:17.201458Z", "shell.execute_reply": "2026-10-21T01:16:17.200509Z" } }, "outputs": [ { "data": { "text/plain": [ "{'cell area against local density, before': 0.297, 'after': -0.096}" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mt.tl.neighbors_local_density(raw, k=10)\n", "area = \"Cells_AreaShape_Area\"\n", "\n", "\n", "def with_density(adata):\n", " \"\"\"How strongly cell area follows crowding, as a correlation.\"\"\"\n", " return round(float(np.corrcoef(adata.obs[\"Metadata_LocalDensity\"], np.asarray(adata[:, area].X).ravel())[0, 1]), 3)\n", "\n", "\n", "before = with_density(raw)\n", "mt.pp.regress_out(raw, keys=[\"Metadata_LocalDensity\"], by=\"Metadata_Plate\")\n", "{\"cell area against local density, before\": before, \"after\": with_density(raw)}" ] }, { "cell_type": "markdown", "id": "c028", "metadata": {}, "source": [ "Isolated cells are larger: cell area follows the mean distance to the ten nearest cells at r = 0.30, and not at all\n", "once the density is regressed out." ] }, { "cell_type": "markdown", "id": "c029", "metadata": {}, "source": [ "## Summary\n", "\n", "- Aggregate late. Even the strongest compound here moves about a third of its cells, and a median dilutes that\n", " by the cells that did not move.\n", "- `method=\"ks\"` and {func}`~mantispy.tl.wasserstein_features` compare whole distributions rather than their\n", " centres. Cells of one well are not independent, so read a cell-level call as a ranking and put the well-level\n", " call beside it.\n", "- Composition and within-state testing answer different questions, and an experiment usually needs both.\n", "- Label the cell cycle and measure crowding before attributing variation to a treatment.\n", "\n", "Next: [Which measurements moved](../phenotypes/which_features_moved.ipynb)." ] } ], "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.13.15" } }, "nbformat": 4, "nbformat_minor": 5 }