domain colocalization code tidied up

Author: Junbum Kim

documentation/IMC Lung Infection.ipynb

@@ -45,12 +45,11 @@
     "import seaborn as sns\n",
     "from utag import utag\n",
     "\n",
-    "from utag.utils import celltype_connectivity, domain_connectivity\n",
+    "from utag.utils import celltype_connectivity\n",
     "from utag.visualize import (\n",
     "    add_scale_box_to_fig,\n",
     "    add_spatial_image,\n",
     "    adj2chord,\n",
-    "    draw_network,\n",
     ")\n",
     "\n",
     "sc.settings.set_figure_params(dpi=200, dpi_save=300, fontsize=6)"
@@ -2715,72 +2714,6 @@
     "adata"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 59,
-   "id": "2b8e1ba2-df46-466d-b1c7-8c43b4b72722",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import matplotlib.pyplot as plt\n",
-    "import networkx as nx\n",
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "import squidpy as sq\n",
-    "from anndata import AnnData\n",
-    "from tqdm import tqdm\n",
-    "\n",
-    "FIG_KWS = dict(bbox_inches=\"tight\", dpi=300)\n",
-    "MAX_BETWEEN_CELL_DIST = 4\n",
-    "\n",
-    "\n",
-    "def measure_cell_type_adjacency(\n",
-    "    adata: AnnData,\n",
-    "    utag_key: str = \"UTAG Label\",\n",
-    "    max_dist: int = 40,\n",
-    "    n_iterations: int = 100,\n",
-    "):\n",
-    "    a_ = adata.copy()\n",
-    "    sq.gr.spatial_neighbors(a_, radius=max_dist, coord_type=\"generic\")\n",
-    "\n",
-    "    G = nx.from_scipy_sparse_matrix(a_.obsp[\"spatial_connectivities\"])\n",
-    "\n",
-    "    utag_map = {i: x for i, x in enumerate(adata.obs[utag_key])}\n",
-    "    nx.set_node_attributes(G, utag_map, name=utag_key)\n",
-    "\n",
-    "    adj, order = nx.linalg.attrmatrix.attr_matrix(G, node_attr=utag_key)\n",
-    "    order = pd.Series(order).astype(adata.obs[utag_key].dtype)\n",
-    "    freqs = pd.DataFrame(adj, order, order).fillna(0) + 1\n",
-    "\n",
-    "    norm_freqs = correct_interaction_background_random(G, freqs, utag_key, n_iterations)\n",
-    "    return norm_freqs\n",
-    "\n",
-    "\n",
-    "def correct_interaction_background_random(\n",
-    "    graph: nx.Graph, freqs: pd.DataFrame, attribute: str, n_iterations: int = 100\n",
-    "):\n",
-    "    values = {x: graph.nodes[x][attribute] for x in graph.nodes}\n",
-    "    shuffled_freqs = list()\n",
-    "    for _ in range(n_iterations):\n",
-    "        g2 = graph.copy()\n",
-    "        shuffled_attr = pd.Series(values).sample(frac=1)\n",
-    "        shuffled_attr.index = values\n",
-    "        nx.set_node_attributes(g2, shuffled_attr.to_dict(), name=attribute)\n",
-    "        rf, rl = nx.linalg.attrmatrix.attr_matrix(g2, node_attr=attribute)\n",
-    "        rl = pd.Series(rl, dtype=freqs.index.dtype)\n",
-    "        shuffled_freqs.append(pd.DataFrame(rf, index=rl, columns=rl))#.fillna(0) + 1)\n",
-    "    shuffled_freq = pd.concat(shuffled_freqs)\n",
-    "    shuffled_freq = shuffled_freq.groupby(level=0).sum()\n",
-    "    shuffled_freq = shuffled_freq.fillna(0) + 1\n",
-    "    \n",
-    "    fl = np.log((freqs / freqs.values.sum()))\n",
-    "    sl = np.log((shuffled_freq / shuffled_freq.values.sum()))\n",
-    "    # make sure both contain all edges/nodes\n",
-    "    fl = fl.reindex(sl.index, axis=0).reindex(sl.index, axis=1)\n",
-    "    sl = sl.reindex(fl.index, axis=0).reindex(fl.index, axis=1)\n",
-    "    return fl - sl"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 107,
@@ -2796,6 +2729,8 @@
     }
    ],
    "source": [
+    "from utag.utils import measure_per_domain_cell_type_colocalization\n",
+    "\n",
     "image_key = \"roi\"\n",
     "utag_key = \"UTAG Label\"\n",
     "\n",
@@ -2834,107 +2769,6 @@
     "adata.obs[\"phenotypes\"].cat.categories"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 341,
-   "id": "d55fb487-8ab9-45e7-9bb3-306404b3dc9e",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def draw_network(\n",
-    "    adata: AnnData,\n",
-    "    node_key: str = \"UTAG Label\",\n",
-    "    adjacency_matrix_key: str = \"UTAG Label_domain_adjacency_matrix\",\n",
-    "    figsize: tuple = (11, 11),\n",
-    "    dpi: int = 200,\n",
-    "    font_size: int = 12,\n",
-    "    node_size_min: int = 1000,\n",
-    "    node_size_max: int = 3000,\n",
-    "    edge_weight: float = 5,\n",
-    "    edge_weight_baseline: float = 1,\n",
-    "    log_transform: bool = True,\n",
-    "    ax=None,\n",
-    "):\n",
-    "    import networkx as nx\n",
-    "\n",
-    "    s1 = adata.obs.groupby(node_key).count()\n",
-    "    s1 = s1[s1.columns[0]]\n",
-    "    node_size = s1.values\n",
-    "    node_size = (node_size - node_size.min()) / (node_size.max() - node_size.min()) * (\n",
-    "        node_size_max - node_size_min\n",
-    "    ) + node_size_min\n",
-    "\n",
-    "    if ax == None:\n",
-    "        fig = plt.figure(figsize=figsize, dpi=dpi)\n",
-    "    G = nx.from_numpy_matrix(\n",
-    "        np.matrix(adata.uns[adjacency_matrix_key]), create_using=nx.Graph\n",
-    "    )\n",
-    "    G = nx.relabel.relabel_nodes(\n",
-    "        G, {i: label for i, label in enumerate(adata.uns[adjacency_matrix_key].index)}\n",
-    "    )\n",
-    "\n",
-    "    edges, weights = zip(*nx.get_edge_attributes(G, \"weight\").items())\n",
-    "\n",
-    "    weights = np.array(list(weights))\n",
-    "    weights = (weights - weights.min()) / (\n",
-    "        weights.max() - weights.min()\n",
-    "    ) * edge_weight + edge_weight_baseline\n",
-    "\n",
-    "    if log_transform:\n",
-    "        weights = np.log(np.array(list(weights)) + 1)\n",
-    "    else:\n",
-    "        weights = np.array(list(weights))\n",
-    "\n",
-    "    weights = tuple(weights.tolist())\n",
-    "\n",
-    "    #pos = nx.spectral_layout(G, weight = 'weight')\n",
-    "    pos = nx.spring_layout(G, weight=\"weight\", seed=42, k=1)\n",
-    "\n",
-    "    if ax:\n",
-    "        nx.draw(\n",
-    "            G,\n",
-    "            pos,\n",
-    "            node_color=\"w\",\n",
-    "            edgelist=edges,\n",
-    "            edge_color=weights,\n",
-    "            width=weights,\n",
-    "            edge_cmap=plt.cm.coolwarm,\n",
-    "            with_labels=True,\n",
-    "            font_size=font_size,\n",
-    "            node_size=node_size,\n",
-    "            ax=ax,\n",
-    "        )\n",
-    "    else:\n",
-    "        nx.draw(\n",
-    "            G,\n",
-    "            pos,\n",
-    "            node_color=\"w\",\n",
-    "            edgelist=edges,\n",
-    "            edge_color=weights,\n",
-    "            width=weights,\n",
-    "            edge_cmap=plt.cm.coolwarm,\n",
-    "            with_labels=True,\n",
-    "            font_size=font_size,\n",
-    "            node_size=node_size,\n",
-    "        )\n",
-    "\n",
-    "    if ax == None:\n",
-    "        ax = plt.gca()\n",
-    "\n",
-    "    color_key = node_key + \"_colors\"\n",
-    "    if color_key in adata.uns:\n",
-    "        ax.collections[0].set_edgecolor(adata.uns[color_key])\n",
-    "        ax.collections[0].set_facecolor(adata.uns[color_key])\n",
-    "    else:\n",
-    "        ax.collections[0].set_edgecolor(\"lightgray\")\n",
-    "    ax.collections[0].set_linewidth(3)\n",
-    "    ax.set_xlim([1.3 * x for x in ax.get_xlim()])\n",
-    "    ax.set_ylim([1 * y for y in ax.get_ylim()])\n",
-    "\n",
-    "    if ax == None:\n",
-    "        return fig"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 154,
@@ -3016,26 +2850,6 @@
     "plt.show()"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 122,
-   "id": "14a6556e-e597-459d-ac3f-5d5c432e2db2",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "array(['Epithelial (Airway / AT)', 'Fibroblasts (Airway Wall)',\n",
-       "       'Alveolar', 'Immune', 'Vessel', 'Neutrophils'], dtype=object)"
-      ]
-     },
-     "execution_count": 122,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": []
-  },
   {
    "cell_type": "code",
    "execution_count": 174,
@@ -3754,133 +3568,6 @@
     "pval_df.sort_values('adj p-val').head(50)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 67,
-   "id": "9aa8c11f-ef17-4332-ba16-ae4c662cb20b",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import networkx as nx\n",
-    "from anndata import AnnData\n",
-    "\n",
-    "\n",
-    "def domain_connectivity(\n",
-    "    adata: AnnData,\n",
-    "    max_dist: int = 40,\n",
-    "    slide_key: str = \"Slide\",\n",
-    "    domain_key: str = \"UTAG Label\",\n",
-    ") -> AnnData:\n",
-    "    import numpy as np\n",
-    "    import squidpy as sq\n",
-    "    from tqdm import tqdm\n",
-    "\n",
-    "    order = pd.Series(adata.obs[domain_key].unique()).astype(\n",
-    "        adata.obs[domain_key].dtype\n",
-    "    )\n",
-    "    global_pairwise_connection = (\n",
-    "        pd.DataFrame(\n",
-    "            np.zeros(shape=(len(order), len(order))), index=order, columns=order\n",
-    "        )\n",
-    "        .astype(int)\n",
-    "        .fillna(0)\n",
-    "    )\n",
-    "\n",
-    "    for slide in tqdm(adata.obs[slide_key].unique()):\n",
-    "        adata_batch = adata[adata.obs[slide_key] == slide].copy()\n",
-    "\n",
-    "        sq.gr.spatial_neighbors(adata_batch, radius=max_dist, coord_type=\"generic\")\n",
-    "        G = nx.from_scipy_sparse_matrix(adata_batch.obsp[\"spatial_connectivities\"])\n",
-    "\n",
-    "        utag_map = {i: x for i, x in enumerate(adata.obs[domain_key])}\n",
-    "        nx.set_node_attributes(G, utag_map, name=domain_key)\n",
-    "\n",
-    "        adj, order = nx.linalg.attrmatrix.attr_matrix(G, node_attr=domain_key)\n",
-    "        order = pd.Series(order).astype(adata.obs[domain_key].dtype)\n",
-    "        freqs = pd.DataFrame(adj, order, order).astype(int)\n",
-    "\n",
-    "        freqs = freqs.fillna(0)\n",
-    "        global_pairwise_connection = global_pairwise_connection + freqs\n",
-    "    adata.uns[f\"{domain_key}_domain_adjacency_matrix\"] = global_pairwise_connection\n",
-    "    return adata"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 103,
-   "id": "c46298e1-07e3-457a-8428-4457ff449498",
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "100%|██████████| 31/31 [00:13<00:00,  2.27it/s]\n",
-      "100%|██████████| 45/45 [00:13<00:00,  3.32it/s]\n",
-      "100%|██████████| 75/75 [00:31<00:00,  2.39it/s]\n",
-      "100%|██████████| 16/16 [00:05<00:00,  2.94it/s]\n",
-      "100%|██████████| 37/37 [00:13<00:00,  2.81it/s]\n",
-      "100%|██████████| 33/33 [00:15<00:00,  2.18it/s]\n"
-     ]
-    }
-   ],
-   "source": [
-    "disease = dict()\n",
-    "for disease_ in results.obs[\"phenotypes\"].unique():\n",
-    "    result_ = results[results.obs[\"phenotypes\"] == disease_].copy()\n",
-    "    disease[disease_] = domain_connectivity(\n",
-    "        result_, slide_key=\"roi\", domain_key=\"UTAG Label\"\n",
-    "    )"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 65,
-   "id": "458de81f-5322-42d3-8abd-b5916b9bb635",
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "'COVID19_late'"
-      ]
-     },
-     "execution_count": 65,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "d_"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "94bc31a4-a0f6-4fc7-a3e3-f20fcb2b6b68",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plt.subplots(1, 6, figsize=(50, 8))\n",
-    "for i, d_ in enumerate(\n",
-    "    [\"Healthy\", \"Flu\", \"ARDS\", \"Pneumonia\", \"COVID19_early\", \"COVID19_late\"]\n",
-    "):\n",
-    "\n",
-    "    fig = draw_network(\n",
-    "        adata=disease[d_],\n",
-    "        font_size=30,\n",
-    "        edge_weight=2,\n",
-    "        edge_weight_baseline=1,\n",
-    "        dpi=50,\n",
-    "        node_size_max=5000,\n",
-    "        node_size_min=500,\n",
-    "        ax=ax[i],\n",
-    "    )\n",
-    "    ax[i].set_title(d_, fontsize=50)\n",
-    "plt.tight_layout()\n",
-    "# plt.savefig(\"figures/infected_lung_Domain_network.pdf\")"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "c2457178-d8f8-4663-8a3c-5766994d0931",
@@ -4870,7 +4557,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.0"
+   "version": "3.7.11"
   },
   "toc-autonumbering": true
  },