{
"cells": [
{
"cell_type": "markdown",
"id": "7c057692-50ec-4b18-b387-7fabf3c36d23",
"metadata": {
"tags": []
},
"source": [
"Tutorial 3: Computing Morphological Distances in Very Large Datasets\n",
"====================================================================\n",
"The Gromov-Wasserstein distance between two cells with 100 points takes about 9 ms to compute on a standard desktop computer. The number of pairs grows quadratically with the number of cells, and so the total runtime can become large in datasets with several thousands of cells.\n",
"\n",
"For very large datasets we provide two tools to reduce the necessary computation, as well as a hybrid of these.\n",
"\n",
"Reference [[1]](#cite_note-1) established several lower bounds for the Gromov-Wasserstein (GW) distance. CAJAL implements one of the fastest bounds, the second lower bound (SLB) [[2]](#cite_note-2). For many downstream analyses, such as clustering and dimensional reduction, it is not crucial to know the exact values between cells with very disparate morphologies, and it is enough to know the precise Gromov-Wasserstein distance only for cells that are close to each other in the morphology space. Since the SLB is a fast lower bound to the GW distance, it can be used to quickly identify pairs of cells that are located far apart in the morphology space so that their precise GW distance does not need to be computed.\n",
"\n",
"Let us illustrate how the computation of the SLB using CAJAL works on the same neuronal dataset as in Tutorial 1. We start with the file of intracellular distances computed in Tutorial 1:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "99f0eefc-3bcf-499e-be2e-1bf3ce01031c",
"metadata": {},
"outputs": [],
"source": [
"# bd = \"/home/jovyan/\" # Base directory\n",
"bd = \"/home/patn/dropbox/Data/AllenInstitute\""
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "c852ce97-8851-4083-a5a6-3ed563dda25a",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "b9b8cd4bea8b4e16b8552d79f524bc0c",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/129286 [00:00\n",
" window.PlotlyConfig = {MathJaxConfig: 'local'};\n",
" if (window.MathJax && window.MathJax.Hub && window.MathJax.Hub.Config) {window.MathJax.Hub.Config({SVG: {font: \"STIX-Web\"}});}\n",
" if (typeof require !== 'undefined') {\n",
" require.undef(\"plotly\");\n",
" requirejs.config({\n",
" paths: {\n",
" 'plotly': ['https://cdn.plot.ly/plotly-2.24.1.min']\n",
" }\n",
" });\n",
" require(['plotly'], function(Plotly) {\n",
" window._Plotly = Plotly;\n",
" });\n",
" }\n",
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import plotly.io as pio\n",
"\n",
"# Choose the adequate plotly renderer for visualizing plotly graphs in your system\n",
"pio.renderers.default = 'notebook_connected'\n",
"# pio.renderers.default = 'iframe'\n",
"\n",
"from cajal.utilities import read_gw_dists, dist_mat_of_dict\n",
"from cajal.run_gw import cell_iterator_csv\n",
"import plotly.express\n",
"\n",
"names, _ = zip(*cell_iterator_csv(join(bd,\"swc_bdad_100pts_euclidean_icdm.csv\")))\n",
"names=list(names)\n",
"\n",
"_, gw_dist_dict = read_gw_dists(join(bd,\"swc_bdad_100pts_euclidean_GW_dmat.csv\"), True)\n",
"gw100_dist_table = dist_mat_of_dict(gw_dist_dict, names, as_squareform=False)\n",
"\n",
"_, slb_dist_dict = read_gw_dists(join(bd,\"slb_dists.csv\"), True)\n",
"slb_dist_table = dist_mat_of_dict(slb_dist_dict, names, as_squareform=False)\n",
"\n",
"fig = plotly.express.scatter(x=slb_dist_table,\n",
" y=gw100_dist_table,\n",
" template=\"simple_white\",\n",
" labels={\n",
" \"x\" : \"SLB\",\n",
" \"y\" : \"GW distance\"})\n",
"fig.update_traces(marker={'size': 1})\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"id": "d7b35df0-2130-492b-91cd-2184f9fb628f",
"metadata": {},
"source": [
"The mean and stadard deviation of the relative difference between the GW distance and the SLB in this example are:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3b11e2e6-53f4-419a-8a60-b741563e3baf",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.22240966232782416\n",
"0.1966012201739646\n"
]
}
],
"source": [
"import numpy\n",
"\n",
"print(numpy.mean((gw100_dist_table-slb_dist_table)/gw100_dist_table))\n",
"print(numpy.std((gw100_dist_table-slb_dist_table)/gw100_dist_table))"
]
},
{
"cell_type": "markdown",
"id": "0682e9b1-6154-4d19-bbfc-c9899b480853",
"metadata": {},
"source": [
"Although the SLB is only a lower bound for the Gromov-Wasserstein distance, using it alone is already fairly accurate as a classifier. Let us repeat the same analysis presented in Tutorial 1 for predicting the molecular type of neurons but using the SLB instead of the GW distance:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "9aaa3495-d339-4125-89dd-b1928b7a7496",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.24036084479003786\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/patn/PGC012_Gromov_Wasserstein/venv/lib/python3.12/site-packages/sklearn/model_selection/_split.py:776: UserWarning:\n",
"\n",
"The least populated class in y has only 1 members, which is less than n_splits=7.\n",
"\n"
]
}
],
"source": [
"import pandas\n",
"from sklearn.neighbors import KNeighborsClassifier\n",
"from sklearn.model_selection import StratifiedKFold, cross_val_predict\n",
"from sklearn.metrics import matthews_corrcoef\n",
"from scipy.spatial.distance import squareform\n",
"\n",
"metadata = pandas.read_csv(join(bd,\"cell_types_specimen_details.csv\"))\n",
"metadata.index = [str(m) for m in metadata[\"specimen__id\"]]\n",
"metadata = metadata.loc[names]\n",
"\n",
"cre_lines = numpy.array(metadata[\"line_name\"])\n",
"\n",
"clf = KNeighborsClassifier(metric=\"precomputed\", n_neighbors=10, weights=\"distance\")\n",
"cv = StratifiedKFold(n_splits=7, shuffle=True)\n",
"\n",
"cvp = cross_val_predict(clf, X=squareform(slb_dist_table), y=cre_lines, cv=cv)\n",
"\n",
"print(matthews_corrcoef(cvp, cre_lines))"
]
},
{
"cell_type": "markdown",
"id": "fd203ccb-816e-45c7-bd12-e00403f126b9",
"metadata": {},
"source": [
"The Matthews correlation coefficient is only slightly lower than the one we derived in the Tutorial 1 using the GW distance."
]
},
{
"cell_type": "markdown",
"id": "aa53fb36-531a-4e03-877a-1a1f9e2db3b7",
"metadata": {},
"source": [
"If the SLB distances between all cells are known, one may envision a simple algorithm[[3]](#cite_note-3) to compute the $k$ nearest neighbors of any given cell under the Gromov-Wasserstein metric. If $k$ is chosen sufficiently large, this is enough to understand the local structure of the morphology space and cluster it."
]
},
{
"cell_type": "markdown",
"id": "bca2d192-a1d0-4f66-9f1d-0ef0d04fb48f",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"The second tool we provide is an implementation of the quantized Gromov-Wasserstein distance proposed by Chowdhury, Miller and Needham[[4]](#cite_note-4). Given cells $X$ and $Y$, the quantized Gromov-Wasserstein distance is given as follows:\n",
"1. Partition the points of $X$ and $Y$ into $n$ clusters. Let $X^n, Y^n$ be the set of medoids of these clusters; $X^n$ can be thought of as the best possible approximation to $X$ in the GW morphology space by a set with at most $n$ points.\n",
"2. Compute the optimal Gromov-Wasserstein transport plan between the subspaces $X^n$ and $Y^n$ formed by the medoids of each cluster.\n",
"3. Extend this to a global transport plan between $X$ and $Y$ by pairing points within paired clusters by their distance from the medoid, and compute the distortion associated to this transport plan.\n",
"This approximation gives an acceptable tradeoff between precision and computation time.\n",
"\n",
"Below, for each cell we cluster the 100 sampled points into 25 clusters."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "104ca14e-3a13-4b24-ab35-ae679ff8465b",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "382ab1e7a60b4d1bb34ee9de88ee2de7",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/509 [00:00 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"_, qgw_dist_dict = read_gw_dists(join(bd,\"quantized_gw.csv\"), header=True)\n",
"qgw_dmat = dist_mat_of_dict(qgw_dist_dict, cell_names=names, as_squareform=False)\n",
"\n",
"fig = plotly.express.scatter(x=qgw_dmat,\n",
" y=gw100_dist_table,\n",
" template=\"simple_white\",\n",
" labels={\n",
" \"x\" : \"Quantized GW distance\",\n",
" \"y\" : \"GW distance\"})\n",
"fig.update_traces(marker={'size': 1})\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"id": "12e12d9b-74b5-4153-8aea-d13fb44baaeb",
"metadata": {},
"source": [
"We can see that the mean and stadard deviation of the relative difference between the GW distance and the quantized GW distance are smaller than those for the SLB:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1174c1bb-dc61-4f62-982f-81602f4a122d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-0.10245317780140377\n",
"0.1328531881517174\n"
]
}
],
"source": [
"print(numpy.mean((gw100_dist_table-qgw_dmat)/gw100_dist_table))\n",
"print(numpy.std((gw100_dist_table-qgw_dmat)/gw100_dist_table))"
]
},
{
"cell_type": "markdown",
"id": "0583e9bc",
"metadata": {},
"source": [
"And benchmark its ability to classify cells according to their morphology using the Matthews correlation coefficient:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "687f9606",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.22899248962487836\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/patn/PGC012_Gromov_Wasserstein/venv/lib/python3.12/site-packages/sklearn/model_selection/_split.py:776: UserWarning:\n",
"\n",
"The least populated class in y has only 1 members, which is less than n_splits=7.\n",
"\n"
]
}
],
"source": [
"cvp = cross_val_predict(clf, X=squareform(qgw_dmat), y=cre_lines, cv=cv)\n",
"\n",
"print(matthews_corrcoef(cvp, cre_lines))"
]
},
{
"cell_type": "markdown",
"id": "d17432d8-fee4-4b34-af11-3fff525404f8",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"Finally, CAJAL implements an approach that combine these two tools (quantized Gromov-Wasserstein and SLB) in an integrated analysis method which allows the user to reduce computation time in three simultaneous ways:\n",
"\n",
"1. by computing only the $k$ nearest neighbors of each cell and estimating the rest roughly using SLB\n",
"\n",
"2. by accepting a small fraction of errors in the reported nearest neighbors list of each cell (i.e., 98% of the nearest neighbors are correct)\n",
"\n",
"3. by using the quantized GW distance as a proxy for the true GW distance"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "56826b9d-6da7-41a7-af9f-78fc63171e08",
"metadata": {},
"outputs": [],
"source": [
"from cajal.combined_slb_qgw import combined_slb_quantized_gw\n",
"combined_slb_quantized_gw(\n",
" join(bd,\"swc_bdad_100pts_euclidean_icdm.csv\"),\n",
" join(bd,\"swc_bdad_100pts_euclidean_gw_estimator.csv\"),\n",
" num_processes=8,\n",
" num_clusters=25,\n",
" accuracy = 0.95,\n",
" nearest_neighbors = 30)"
]
},
{
"cell_type": "markdown",
"id": "d540212b-6705-48bb-80ca-db19f063c62c",
"metadata": {},
"source": [
"The above command does the following:\n",
"\n",
"1. Compute the pairwise SLB between any two cells in the given list of intracell distance matrices.\n",
"\n",
"2. For each cell $X$, identify the 30 estimated nearest neighbors $X_1,\\dots, X_{30}$ based on the SLB distance matrix, and compute the QGW distance $QGW_{25}(X,X_k)$ for all 30 pairs.\n",
"\n",
"3. For each remaining cell pair $X,Y$, we estimate the probability that $QGW_{25}(X,Y)$ is small enough to \"injure\" the existing purported list of nearest neighbors. We sort the cell pairs in descending order by this probability and compute the $QGW_{25}$ distance between pairs in this list until the expected number of \"injuries\" remaining is less than 5% of the nearest neighbor table, so that of the reported 30 nearest nearest neighbors to each point, 28.5 are expected to be correct.\n",
"\n",
"4. For all remaining cells we estimate the correct distance based on the SLB."
]
},
{
"cell_type": "markdown",
"id": "9425ebd2-3172-4db3-9032-64c057260b24",
"metadata": {},
"source": [
"As expected, a UMAP representation based on the fast QGW/SLB combined estimator is very close to the UMAP that we computed in Tutorial 1 for the ame data using the full GW distance:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "c35925fd-e960-4fc4-828a-32bd9ba91d79",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/patn/PGC012_Gromov_Wasserstein/venv/lib/python3.12/site-packages/umap/umap_.py:1780: UserWarning:\n",
"\n",
"using precomputed metric; inverse_transform will be unavailable\n",
"\n",
"/home/patn/PGC012_Gromov_Wasserstein/venv/lib/python3.12/site-packages/plotly/express/_core.py:1992: FutureWarning:\n",
"\n",
"When grouping with a length-1 list-like, you will need to pass a length-1 tuple to get_group in a future version of pandas. Pass `(name,)` instead of `name` to silence this warning.\n",
"\n"
]
},
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from umap import umap_\n",
"import pandas\n",
"\n",
"_, slbqgw_dist_dict = read_gw_dists(join(bd,\"swc_bdad_100pts_euclidean_gw_estimator.csv\"), header=False)\n",
"slbqgw_sq_dist = dist_mat_of_dict(slbqgw_dist_dict, cell_names=names)\n",
"\n",
"# Load metadata\n",
"metadata = pandas.read_csv(join(bd,\"cell_types_specimen_details.csv\"))\n",
"metadata.index = [str(m) for m in metadata[\"specimen__id\"]]\n",
"metadata = metadata.loc[names]\n",
"\n",
"# Compute UMAP representation\n",
"reducer = umap_.UMAP(metric=\"precomputed\", random_state=1)\n",
"embedding = reducer.fit_transform(slbqgw_sq_dist)\n",
"\n",
"# Visualize UMAP colored by cortical layer\n",
"plotly.express.scatter(x=embedding[:,0],\n",
" y=embedding[:,1],\n",
" template=\"simple_white\",\n",
" hover_name=[m + \".swc\" for m in names],\n",
" color = metadata[\"structure__layer\"])"
]
},
{
"cell_type": "markdown",
"id": "567565d0-2503-49b6-9439-7cbd2bf24c52",
"metadata": {
"tags": []
},
"source": [
"Notes\n",
"-----\n",
"\n",
"1.[^](#cite_ref-1) Mémoli, F. P. [Gromov–Wasserstein Distances and the Metric Approach to Object Matching.](https://dblp.uni-trier.de/db/journals/focm/focm11.html#Memoli11) Found Comput Math (2011) 11, 417–487.\n",
"\n",
"2.[^](#cite_ref-2) Specifically we have implemented the expression which appears on the right hand side of the first inequality of Corollary 6.2 on page 462, for p = 2. This expression is not directly named in the paper. We do not compute the quantity Mémoli calls SLB, as it is too computationally expensive for our purposes.\n",
"\n",
"3.[^](#cite_ref-3) A simple algorithm for computing the nearest neighbors of a cell in the Gromov-Wasserstein morphology space if the SLB distance is known is as follows:\n",
"1. First, sort all other cells by their SLB2 distance from $c_0$, as $c_1, c_2, c_3,\\dots$.\n",
"2. Next, compute the Gromov-Wasserstein distance $GW(c_0,c_j)$, as $j = 1, 2, 3,\\dots$. Write $e^k_j$ for the $k$-th element of the set $GW(c_0,c_1),GW(c_0,c_2),\\dots,GW(c_0,c_j)$ when these are ordered from least to greatest. ($e^k_j$ is only defined when $k \\leq j$). Continue computing $GW(c_0,c_j)$ until $j$ reaches a value $\\ell$ such that for all $i> \\ell$, $SLB(c_0,c_i) > e^k_\\ell$. Because SLB is a lower bound, at this point, the $k$ nearest neighbors of $c_0$ are contained in the set $\\{c_1,\\dots, c_\\ell\\}$.\n",
"\n",
"4.[^](#cite_ref-4) Chowdhury, S., Miller, D., Needham, T. (2021). Quantized Gromov-Wasserstein. In: Oliver, N., Pérez-Cruz, F., Kramer, S., Read, J., Lozano, J.A. (eds) Machine Learning and Knowledge Discovery in Databases. Research Track. ECML PKDD 2021. Lecture Notes in Computer Science(), vol 12977. Springer, Cham. https://doi.org/10.1007/978-3-030-86523-8_49"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "venv",
"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.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}