weighted_Adj_mat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ea7c84c9",
   "metadata": {},
   "source": [
    "# Following [SpaGE Tutorial](https://github.com/tabdelaal/SpaGE/blob/master/SpaGE_Tutorial.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de4cc13",
   "metadata": {},
   "source": [
    "Integration of *osmFISH* spatial data with the *AllenSSp* scRNA-seq data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "276e8209",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import loompy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import scipy\n",
    "import scipy.stats as st\n",
    "from scipy import linalg\n",
    "from scipy import sparse as sp\n",
    "import pickle\n",
    "from sklearn.metrics.pairwise import euclidean_distances\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "86e3ff8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = '/Volumes/LaCie/school/combine_lab/SpaGE/SpaGE_Datasets/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b17cbb3e",
   "metadata": {},
   "source": [
    "# import spatial data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9cd87735",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = loompy.connect(path + 'Spatial/osmFISH/osmFISH_SScortex_mouse_all_cells.loom')\n",
    "FISH_Genes = ds.ra['Gene']   \n",
    "colAtr = ds.ca.keys()\n",
    "\n",
    "df = pd.DataFrame()\n",
    "for i in colAtr:\n",
    "    df[i] = ds.ca[i]\n",
    "\n",
    "osmFISH_meta = df.iloc[np.where(df.Valid == 1)[0], :]\n",
    "osmFISH_data = ds[:,:]\n",
    "osmFISH_data = osmFISH_data[:,np.where(df.Valid == 1)[0]]\n",
    "osmFISH_data = pd.DataFrame(data= osmFISH_data, index= FISH_Genes)\n",
    "\n",
    "del ds, colAtr, i, df, FISH_Genes\n",
    "\n",
    "# Select cortical regions only to match the AllenSSp dataset\n",
    "Cortex_Regions = ['Layer 2-3 lateral', 'Layer 2-3 medial', 'Layer 3-4', \n",
    "                  'Layer 4','Layer 5', 'Layer 6', 'Pia Layer 1']\n",
    "Cortical = np.stack(i in Cortex_Regions for i in osmFISH_meta.Region)\n",
    "\n",
    "osmFISH_meta = osmFISH_meta.iloc[Cortical,:]\n",
    "osmFISH_data = osmFISH_data.iloc[:,Cortical]\n",
    "del Cortex_Regions,Cortical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fbd16aa3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CellID</th>\n",
       "      <th>ClusterID</th>\n",
       "      <th>ClusterName</th>\n",
       "      <th>Region</th>\n",
       "      <th>Total_molecules</th>\n",
       "      <th>Valid</th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>_tSNE_1</th>\n",
       "      <th>_tSNE_2</th>\n",
       "      <th>size_pix</th>\n",
       "      <th>size_um2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>778</td>\n",
       "      <td>18</td>\n",
       "      <td>Inhibitory CP</td>\n",
       "      <td>Layer 6</td>\n",
       "      <td>390</td>\n",
       "      <td>1</td>\n",
       "      <td>18171.230942</td>\n",
       "      <td>24590.795275</td>\n",
       "      <td>-58.132385</td>\n",
       "      <td>5.219181</td>\n",
       "      <td>60911.0</td>\n",
       "      <td>257.348975</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3642</td>\n",
       "      <td>18</td>\n",
       "      <td>Inhibitory CP</td>\n",
       "      <td>Layer 4</td>\n",
       "      <td>405</td>\n",
       "      <td>1</td>\n",
       "      <td>11247.433715</td>\n",
       "      <td>36626.892415</td>\n",
       "      <td>-61.572853</td>\n",
       "      <td>-1.185364</td>\n",
       "      <td>44143.0</td>\n",
       "      <td>186.504175</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>769</td>\n",
       "      <td>18</td>\n",
       "      <td>Inhibitory CP</td>\n",
       "      <td>Layer 4</td>\n",
       "      <td>491</td>\n",
       "      <td>1</td>\n",
       "      <td>18374.295796</td>\n",
       "      <td>26286.856654</td>\n",
       "      <td>-66.862244</td>\n",
       "      <td>-1.766500</td>\n",
       "      <td>41014.0</td>\n",
       "      <td>173.284150</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>5870</td>\n",
       "      <td>18</td>\n",
       "      <td>Inhibitory CP</td>\n",
       "      <td>Layer 4</td>\n",
       "      <td>518</td>\n",
       "      <td>1</td>\n",
       "      <td>4197.789280</td>\n",
       "      <td>39465.355903</td>\n",
       "      <td>-59.923431</td>\n",
       "      <td>-3.739808</td>\n",
       "      <td>42376.0</td>\n",
       "      <td>179.038600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5682</td>\n",
       "      <td>18</td>\n",
       "      <td>Inhibitory CP</td>\n",
       "      <td>Layer 2-3 medial</td>\n",
       "      <td>210</td>\n",
       "      <td>1</td>\n",
       "      <td>3757.866349</td>\n",
       "      <td>46281.167166</td>\n",
       "      <td>-64.769470</td>\n",
       "      <td>5.166343</td>\n",
       "      <td>23922.0</td>\n",
       "      <td>101.070450</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4832</th>\n",
       "      <td>3635</td>\n",
       "      <td>25</td>\n",
       "      <td>Vascular Smooth Muscle</td>\n",
       "      <td>Layer 3-4</td>\n",
       "      <td>59</td>\n",
       "      <td>1</td>\n",
       "      <td>11457.806019</td>\n",
       "      <td>38036.523162</td>\n",
       "      <td>-12.763592</td>\n",
       "      <td>7.068469</td>\n",
       "      <td>7993.0</td>\n",
       "      <td>33.770425</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4833</th>\n",
       "      <td>3336</td>\n",
       "      <td>25</td>\n",
       "      <td>Vascular Smooth Muscle</td>\n",
       "      <td>Pia Layer 1</td>\n",
       "      <td>144</td>\n",
       "      <td>1</td>\n",
       "      <td>9214.047907</td>\n",
       "      <td>45314.585476</td>\n",
       "      <td>0.060634</td>\n",
       "      <td>30.242371</td>\n",
       "      <td>14865.0</td>\n",
       "      <td>62.804625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4834</th>\n",
       "      <td>6462</td>\n",
       "      <td>25</td>\n",
       "      <td>Vascular Smooth Muscle</td>\n",
       "      <td>Layer 6</td>\n",
       "      <td>236</td>\n",
       "      <td>1</td>\n",
       "      <td>4476.301656</td>\n",
       "      <td>31097.786260</td>\n",
       "      <td>3.578825</td>\n",
       "      <td>28.862700</td>\n",
       "      <td>24835.0</td>\n",
       "      <td>104.927875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4835</th>\n",
       "      <td>260</td>\n",
       "      <td>25</td>\n",
       "      <td>Vascular Smooth Muscle</td>\n",
       "      <td>Layer 2-3 lateral</td>\n",
       "      <td>67</td>\n",
       "      <td>1</td>\n",
       "      <td>18796.948664</td>\n",
       "      <td>36603.475262</td>\n",
       "      <td>-11.011883</td>\n",
       "      <td>5.821001</td>\n",
       "      <td>8873.0</td>\n",
       "      <td>37.488425</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4837</th>\n",
       "      <td>6431</td>\n",
       "      <td>25</td>\n",
       "      <td>Vascular Smooth Muscle</td>\n",
       "      <td>Layer 6</td>\n",
       "      <td>147</td>\n",
       "      <td>1</td>\n",
       "      <td>884.041752</td>\n",
       "      <td>33406.078418</td>\n",
       "      <td>-20.118292</td>\n",
       "      <td>38.138481</td>\n",
       "      <td>19604.0</td>\n",
       "      <td>82.826900</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3405 rows × 12 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     CellID  ClusterID             ClusterName             Region  \\\n",
       "0       778         18           Inhibitory CP            Layer 6   \n",
       "2      3642         18           Inhibitory CP            Layer 4   \n",
       "5       769         18           Inhibitory CP            Layer 4   \n",
       "7      5870         18           Inhibitory CP            Layer 4   \n",
       "9      5682         18           Inhibitory CP   Layer 2-3 medial   \n",
       "...     ...        ...                     ...                ...   \n",
       "4832   3635         25  Vascular Smooth Muscle          Layer 3-4   \n",
       "4833   3336         25  Vascular Smooth Muscle        Pia Layer 1   \n",
       "4834   6462         25  Vascular Smooth Muscle            Layer 6   \n",
       "4835    260         25  Vascular Smooth Muscle  Layer 2-3 lateral   \n",
       "4837   6431         25  Vascular Smooth Muscle            Layer 6   \n",
       "\n",
       "      Total_molecules  Valid             X             Y    _tSNE_1  \\\n",
       "0                 390      1  18171.230942  24590.795275 -58.132385   \n",
       "2                 405      1  11247.433715  36626.892415 -61.572853   \n",
       "5                 491      1  18374.295796  26286.856654 -66.862244   \n",
       "7                 518      1   4197.789280  39465.355903 -59.923431   \n",
       "9                 210      1   3757.866349  46281.167166 -64.769470   \n",
       "...               ...    ...           ...           ...        ...   \n",
       "4832               59      1  11457.806019  38036.523162 -12.763592   \n",
       "4833              144      1   9214.047907  45314.585476   0.060634   \n",
       "4834              236      1   4476.301656  31097.786260   3.578825   \n",
       "4835               67      1  18796.948664  36603.475262 -11.011883   \n",
       "4837              147      1    884.041752  33406.078418 -20.118292   \n",
       "\n",
       "        _tSNE_2  size_pix    size_um2  \n",
       "0      5.219181   60911.0  257.348975  \n",
       "2     -1.185364   44143.0  186.504175  \n",
       "5     -1.766500   41014.0  173.284150  \n",
       "7     -3.739808   42376.0  179.038600  \n",
       "9      5.166343   23922.0  101.070450  \n",
       "...         ...       ...         ...  \n",
       "4832   7.068469    7993.0   33.770425  \n",
       "4833  30.242371   14865.0   62.804625  \n",
       "4834  28.862700   24835.0  104.927875  \n",
       "4835   5.821001    8873.0   37.488425  \n",
       "4837  38.138481   19604.0   82.826900  \n",
       "\n",
       "[3405 rows x 12 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "osmFISH_meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0cc300f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalization of ST gene expression\n",
    "cell_count = np.sum(osmFISH_data,axis=0)\n",
    "def Log_Norm_spatial(x):\n",
    "    return np.log(((x/np.sum(x))*np.median(cell_count)) + 1)\n",
    "\n",
    "osmFISH_data_norm = osmFISH_data.apply(Log_Norm_spatial,axis=0)\n",
    "\n",
    "# scale\n",
    "genes = st.zscore(osmFISH_data_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "854e4720",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>CellID</th>\n",
       "      <th>X</th>\n",
       "      <th>Y</th>\n",
       "      <th>Gad2</th>\n",
       "      <th>Slc32a1</th>\n",
       "      <th>Crhbp</th>\n",
       "      <th>Cnr1</th>\n",
       "      <th>Vip</th>\n",
       "      <th>Cpne5</th>\n",
       "      <th>Pthlh</th>\n",
       "      <th>...</th>\n",
       "      <th>Ctps</th>\n",
       "      <th>Anln</th>\n",
       "      <th>Mrc1</th>\n",
       "      <th>Hexb</th>\n",
       "      <th>Ttr</th>\n",
       "      <th>Foxj1</th>\n",
       "      <th>Vtn</th>\n",
       "      <th>Flt1</th>\n",
       "      <th>Apln</th>\n",
       "      <th>Acta2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>778</td>\n",
       "      <td>18171.230942</td>\n",
       "      <td>24590.795275</td>\n",
       "      <td>1.621371</td>\n",
       "      <td>1.777401</td>\n",
       "      <td>-1.285498</td>\n",
       "      <td>0.169445</td>\n",
       "      <td>-0.417264</td>\n",
       "      <td>0.932325</td>\n",
       "      <td>-1.285498</td>\n",
       "      <td>...</td>\n",
       "      <td>0.817961</td>\n",
       "      <td>0.050672</td>\n",
       "      <td>-0.237032</td>\n",
       "      <td>0.876704</td>\n",
       "      <td>-0.237032</td>\n",
       "      <td>-0.417264</td>\n",
       "      <td>-0.635124</td>\n",
       "      <td>-1.285498</td>\n",
       "      <td>-0.237032</td>\n",
       "      <td>-0.417264</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3642</td>\n",
       "      <td>11247.433715</td>\n",
       "      <td>36626.892415</td>\n",
       "      <td>2.723553</td>\n",
       "      <td>2.088070</td>\n",
       "      <td>0.162081</td>\n",
       "      <td>-0.609545</td>\n",
       "      <td>-0.353881</td>\n",
       "      <td>-0.609545</td>\n",
       "      <td>-0.609545</td>\n",
       "      <td>...</td>\n",
       "      <td>0.162081</td>\n",
       "      <td>-0.150690</td>\n",
       "      <td>-0.609545</td>\n",
       "      <td>-0.150690</td>\n",
       "      <td>-0.609545</td>\n",
       "      <td>-0.609545</td>\n",
       "      <td>-0.353881</td>\n",
       "      <td>-0.954585</td>\n",
       "      <td>0.017943</td>\n",
       "      <td>-0.609545</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>769</td>\n",
       "      <td>18374.295796</td>\n",
       "      <td>26286.856654</td>\n",
       "      <td>2.459380</td>\n",
       "      <td>2.289477</td>\n",
       "      <td>1.116087</td>\n",
       "      <td>-0.466125</td>\n",
       "      <td>-0.148632</td>\n",
       "      <td>-0.951611</td>\n",
       "      <td>0.186618</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.022988</td>\n",
       "      <td>-0.148632</td>\n",
       "      <td>-0.293912</td>\n",
       "      <td>1.150048</td>\n",
       "      <td>-0.677581</td>\n",
       "      <td>-0.148632</td>\n",
       "      <td>-0.022988</td>\n",
       "      <td>-0.951611</td>\n",
       "      <td>-0.466125</td>\n",
       "      <td>-0.293912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5870</td>\n",
       "      <td>4197.789280</td>\n",
       "      <td>39465.355903</td>\n",
       "      <td>2.253989</td>\n",
       "      <td>1.549556</td>\n",
       "      <td>-0.197524</td>\n",
       "      <td>0.225092</td>\n",
       "      <td>-0.510714</td>\n",
       "      <td>-0.982660</td>\n",
       "      <td>0.135831</td>\n",
       "      <td>...</td>\n",
       "      <td>0.225092</td>\n",
       "      <td>-0.982660</td>\n",
       "      <td>-0.072863</td>\n",
       "      <td>0.515976</td>\n",
       "      <td>-0.072863</td>\n",
       "      <td>-0.717495</td>\n",
       "      <td>-0.717495</td>\n",
       "      <td>-0.982660</td>\n",
       "      <td>-0.717495</td>\n",
       "      <td>-0.072863</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5682</td>\n",
       "      <td>3757.866349</td>\n",
       "      <td>46281.167166</td>\n",
       "      <td>1.496840</td>\n",
       "      <td>2.316222</td>\n",
       "      <td>-0.944582</td>\n",
       "      <td>1.345124</td>\n",
       "      <td>0.013995</td>\n",
       "      <td>-0.944582</td>\n",
       "      <td>-0.944582</td>\n",
       "      <td>...</td>\n",
       "      <td>0.495107</td>\n",
       "      <td>-0.354909</td>\n",
       "      <td>0.013995</td>\n",
       "      <td>0.283120</td>\n",
       "      <td>-0.354909</td>\n",
       "      <td>-0.354909</td>\n",
       "      <td>-0.354909</td>\n",
       "      <td>-0.354909</td>\n",
       "      <td>0.283120</td>\n",
       "      <td>-0.944582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3400</th>\n",
       "      <td>3635</td>\n",
       "      <td>11457.806019</td>\n",
       "      <td>38036.523162</td>\n",
       "      <td>0.309345</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>0.309345</td>\n",
       "      <td>0.309345</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>...</td>\n",
       "      <td>1.309683</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>1.086314</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>-0.843207</td>\n",
       "      <td>0.309345</td>\n",
       "      <td>1.086314</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3401</th>\n",
       "      <td>3336</td>\n",
       "      <td>9214.047907</td>\n",
       "      <td>45314.585476</td>\n",
       "      <td>1.396286</td>\n",
       "      <td>-0.977202</td>\n",
       "      <td>-0.977202</td>\n",
       "      <td>1.229384</td>\n",
       "      <td>1.229384</td>\n",
       "      <td>-0.977202</td>\n",
       "      <td>-0.977202</td>\n",
       "      <td>...</td>\n",
       "      <td>0.381686</td>\n",
       "      <td>0.759097</td>\n",
       "      <td>0.110551</td>\n",
       "      <td>-0.278649</td>\n",
       "      <td>-0.278649</td>\n",
       "      <td>-0.977202</td>\n",
       "      <td>1.024502</td>\n",
       "      <td>0.381686</td>\n",
       "      <td>-0.278649</td>\n",
       "      <td>2.598400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3402</th>\n",
       "      <td>6462</td>\n",
       "      <td>4476.301656</td>\n",
       "      <td>31097.786260</td>\n",
       "      <td>0.232061</td>\n",
       "      <td>-0.423666</td>\n",
       "      <td>-1.342403</td>\n",
       "      <td>-0.423666</td>\n",
       "      <td>-0.156431</td>\n",
       "      <td>-1.342403</td>\n",
       "      <td>-0.423666</td>\n",
       "      <td>...</td>\n",
       "      <td>0.055889</td>\n",
       "      <td>0.232061</td>\n",
       "      <td>2.117512</td>\n",
       "      <td>-0.784520</td>\n",
       "      <td>0.830838</td>\n",
       "      <td>0.232061</td>\n",
       "      <td>-0.423666</td>\n",
       "      <td>-1.342403</td>\n",
       "      <td>0.232061</td>\n",
       "      <td>1.848433</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3403</th>\n",
       "      <td>260</td>\n",
       "      <td>18796.948664</td>\n",
       "      <td>36603.475262</td>\n",
       "      <td>0.632525</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>1.198118</td>\n",
       "      <td>0.632525</td>\n",
       "      <td>0.129437</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>...</td>\n",
       "      <td>0.129437</td>\n",
       "      <td>0.129437</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>0.129437</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>-1.041085</td>\n",
       "      <td>0.632525</td>\n",
       "      <td>1.389190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3404</th>\n",
       "      <td>6431</td>\n",
       "      <td>884.041752</td>\n",
       "      <td>33406.078418</td>\n",
       "      <td>0.647033</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>-0.430347</td>\n",
       "      <td>-0.027124</td>\n",
       "      <td>-0.027124</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>...</td>\n",
       "      <td>0.647033</td>\n",
       "      <td>1.036189</td>\n",
       "      <td>0.471020</td>\n",
       "      <td>-0.027124</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>-0.027124</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>-1.149300</td>\n",
       "      <td>1.584745</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>3405 rows × 36 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     CellID             X             Y      Gad2   Slc32a1     Crhbp  \\\n",
       "0       778  18171.230942  24590.795275  1.621371  1.777401 -1.285498   \n",
       "1      3642  11247.433715  36626.892415  2.723553  2.088070  0.162081   \n",
       "2       769  18374.295796  26286.856654  2.459380  2.289477  1.116087   \n",
       "3      5870   4197.789280  39465.355903  2.253989  1.549556 -0.197524   \n",
       "4      5682   3757.866349  46281.167166  1.496840  2.316222 -0.944582   \n",
       "...     ...           ...           ...       ...       ...       ...   \n",
       "3400   3635  11457.806019  38036.523162  0.309345 -0.843207 -0.843207   \n",
       "3401   3336   9214.047907  45314.585476  1.396286 -0.977202 -0.977202   \n",
       "3402   6462   4476.301656  31097.786260  0.232061 -0.423666 -1.342403   \n",
       "3403    260  18796.948664  36603.475262  0.632525 -1.041085 -1.041085   \n",
       "3404   6431    884.041752  33406.078418  0.647033 -1.149300 -1.149300   \n",
       "\n",
       "          Cnr1       Vip     Cpne5     Pthlh  ...      Ctps      Anln  \\\n",
       "0     0.169445 -0.417264  0.932325 -1.285498  ...  0.817961  0.050672   \n",
       "1    -0.609545 -0.353881 -0.609545 -0.609545  ...  0.162081 -0.150690   \n",
       "2    -0.466125 -0.148632 -0.951611  0.186618  ... -0.022988 -0.148632   \n",
       "3     0.225092 -0.510714 -0.982660  0.135831  ...  0.225092 -0.982660   \n",
       "4     1.345124  0.013995 -0.944582 -0.944582  ...  0.495107 -0.354909   \n",
       "...        ...       ...       ...       ...  ...       ...       ...   \n",
       "3400 -0.843207  0.309345  0.309345 -0.843207  ...  1.309683 -0.843207   \n",
       "3401  1.229384  1.229384 -0.977202 -0.977202  ...  0.381686  0.759097   \n",
       "3402 -0.423666 -0.156431 -1.342403 -0.423666  ...  0.055889  0.232061   \n",
       "3403  1.198118  0.632525  0.129437 -1.041085  ...  0.129437  0.129437   \n",
       "3404 -0.430347 -0.027124 -0.027124 -1.149300  ...  0.647033  1.036189   \n",
       "\n",
       "          Mrc1      Hexb       Ttr     Foxj1       Vtn      Flt1      Apln  \\\n",
       "0    -0.237032  0.876704 -0.237032 -0.417264 -0.635124 -1.285498 -0.237032   \n",
       "1    -0.609545 -0.150690 -0.609545 -0.609545 -0.353881 -0.954585  0.017943   \n",
       "2    -0.293912  1.150048 -0.677581 -0.148632 -0.022988 -0.951611 -0.466125   \n",
       "3    -0.072863  0.515976 -0.072863 -0.717495 -0.717495 -0.982660 -0.717495   \n",
       "4     0.013995  0.283120 -0.354909 -0.354909 -0.354909 -0.354909  0.283120   \n",
       "...        ...       ...       ...       ...       ...       ...       ...   \n",
       "3400  1.086314 -0.843207 -0.843207 -0.843207 -0.843207 -0.843207  0.309345   \n",
       "3401  0.110551 -0.278649 -0.278649 -0.977202  1.024502  0.381686 -0.278649   \n",
       "3402  2.117512 -0.784520  0.830838  0.232061 -0.423666 -1.342403  0.232061   \n",
       "3403 -1.041085 -1.041085  0.129437 -1.041085 -1.041085 -1.041085  0.632525   \n",
       "3404  0.471020 -0.027124 -1.149300 -1.149300 -0.027124 -1.149300 -1.149300   \n",
       "\n",
       "         Acta2  \n",
       "0    -0.417264  \n",
       "1    -0.609545  \n",
       "2    -0.293912  \n",
       "3    -0.072863  \n",
       "4    -0.944582  \n",
       "...        ...  \n",
       "3400  1.086314  \n",
       "3401  2.598400  \n",
       "3402  1.848433  \n",
       "3403  1.389190  \n",
       "3404  1.584745  \n",
       "\n",
       "[3405 rows x 36 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cell_st = osmFISH_meta.loc[:,['CellID', 'X', 'Y']]\n",
    "\n",
    "# add gene expression info for each cell at spot i\n",
    "cell_spot = cell_st.merge(genes.T, left_index=True, right_index=True, how='left')\n",
    "\n",
    "cell_spt = cell_spot.reset_index(drop=True) # reset index\n",
    "cell_spt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fbdfc36a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_express = cell_spt.set_index('CellID', drop=True)\n",
    "gene_express = gene_express.drop(['X', 'Y'], axis=1)\n",
    "gene_express\n",
    "\n",
    "# cell_spot.to_csv(\"ST_gene_expression.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a45ec4cd",
   "metadata": {},
   "source": [
    "# Get ddf(x) for all training samples\n",
    "\n",
    "## Create adj mat (inspired by [GCNG](https://github.com/xiaoyeye/GCNG/blob/master/data_generation_interaction_ten_fold.py) and [SpaGCN](https://github.com/jianhuupenn/SpaGCN))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fd7ce39",
   "metadata": {},
   "source": [
    "From [GCNG](https://github.com/xiaoyeye/GCNG/blob/2fba06bac64e97744a9a8d2a1505aa5c2117730c/data_generation_interaction_ten_fold.py#L19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a224863",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### from GCNG ####\n",
    "\n",
    "############ the distribution of distance\n",
    "distance_list = []\n",
    "print ('calculating distance matrix, it takes a while')\n",
    "\n",
    "# compute euclidean distance, considering both physical location and gene expression\n",
    "for j in range(cell_st.shape[0]):\n",
    "    for i in range(cell_st.shape[0]):\n",
    "        if i!=j:\n",
    "            distance_list.append(np.linalg.norm(cell_spot.iloc[j, 1:]-cell_spot.iloc[i, 1:]))\n",
    "\n",
    "# save out file\n",
    "np.save(path+'dist_array_eucl_xyz.npy',np.array(distance_list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5979f779",
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_array = np.array(distance_list)\n",
    "distance_array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a594ba3b",
   "metadata": {},
   "source": [
    "import saved distance array (instead of re-running since it takes a long time to compute)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d82e1305",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(11590620,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([13885.48235918,  1708.18015734, 20408.56783459, ...,\n",
       "       14532.77639917,  4269.96297519, 18196.03282577])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distance_array = np.load(path+'dist_array_eucl_xyz.npy')\n",
    "print(distance_array.shape)\n",
    "distance_array"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4aad6f8",
   "metadata": {},
   "source": [
    "from [spektral](https://github.com/danielegrattarola/spektral/blob/1dbdf9e919ff4606ee5375567a9b993687c95bcd/spektral/utils/convolution.py#L47) -- unable to import the package (received 'illegal hardware' error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6a5e2fe9",
   "metadata": {},
   "outputs": [],
   "source": [
    "### from spektral ###\n",
    "def degree_power(A, k):\n",
    "    r\"\"\"\n",
    "    Computes \\(\\D^{k}\\) from the given adjacency matrix. Useful for computing\n",
    "    normalised Laplacian.\n",
    "    :param A: rank 2 array or sparse matrix.\n",
    "    :param k: exponent to which elevate the degree matrix.\n",
    "    :return: if A is a dense array, a dense array; if A is sparse, a sparse\n",
    "    matrix in DIA format.\n",
    "    \"\"\"\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")\n",
    "        degrees = np.power(np.array(A.sum(1)), k).ravel()\n",
    "    degrees[np.isinf(degrees)] = 0.0\n",
    "    if sp.issparse(A):\n",
    "        D = sp.diags(degrees)\n",
    "    else:\n",
    "        D = np.diag(degrees)\n",
    "    return D\n",
    "\n",
    "def normalized_adjacency(A, symmetric=True):\n",
    "    r\"\"\"\n",
    "    Normalizes the given adjacency matrix using the degree matrix as either\n",
    "    \\(\\D^{-1}\\A\\) (row normalization) or \\(\\D^{-1/2}\\A\\D^{-1/2}\\) (symmetric normalization).\n",
    "    :param A: rank 2 array or sparse matrix;\n",
    "    :param symmetric: boolean, compute symmetric normalization;\n",
    "    :return: the normalized adjacency matrix.\n",
    "    \n",
    "    motivation for symmetric norm: takes into account both # of neighbours and their own # of neighbours\n",
    "    \"\"\"\n",
    "    if symmetric:\n",
    "        normalized_D = degree_power(A, -0.5)\n",
    "        return normalized_D.dot(A).dot(normalized_D)\n",
    "    else:\n",
    "        normalized_D = degree_power(A, -1.0) # row normalized: consider only nieghbour connected to node v\n",
    "        return normalized_D.dot(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e0356e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_matrix = euclidean_distances(cell_spot.iloc[:, 1:], cell_spot.iloc[:, 1:])\n",
    "\n",
    "############### get normalized sparse adjacent matrix\n",
    "dist_symnorm = normalized_adjacency(distance_matrix, symmetric=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a5b42b7b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[    0.         13885.48235918  1708.18015734 ... 15162.19227364\n",
      "  12028.96694124 19405.05509302]\n",
      " [13885.48235918     0.         12558.20499189 ...  8741.81396431\n",
      "   7549.55471867 10852.35417449]\n",
      " [ 1708.18015734 12558.20499189     0.         ... 14707.1188864\n",
      "  10325.27560572 18883.65330954]\n",
      " ...\n",
      " [15162.19227364  8741.81396431 14707.1188864  ...     0.\n",
      "  15342.54223372  4269.96297519]\n",
      " [12028.96694124  7549.55471867 10325.27560572 ... 15342.54223372\n",
      "      0.         18196.03282577]\n",
      " [19405.05509302 10852.35417449 18883.65330954 ...  4269.96297519\n",
      "  18196.03282577     0.        ]]\n"
     ]
    }
   ],
   "source": [
    "print(distance_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "631425e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00000000e+00 3.85134799e-04 3.76706298e-05 ... 3.66614173e-04\n",
      "  2.89538098e-04 4.37072874e-04]\n",
      " [3.85134799e-04 0.00000000e+00 3.61531199e-04 ... 2.75929298e-04\n",
      "  2.37218049e-04 3.19089077e-04]\n",
      " [3.76706298e-05 3.61531199e-04 0.00000000e+00 ... 3.69097647e-04\n",
      "  2.57955898e-04 4.41460057e-04]\n",
      " ...\n",
      " [3.66614173e-04 2.75929298e-04 3.69097647e-04 ... 0.00000000e+00\n",
      "  4.20261086e-04 1.09447929e-04]\n",
      " [2.89538098e-04 2.37218049e-04 2.57955898e-04 ... 4.20261086e-04\n",
      "  0.00000000e+00 4.64291023e-04]\n",
      " [4.37072874e-04 3.19089077e-04 4.41460057e-04 ... 1.09447929e-04\n",
      "  4.64291023e-04 0.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "print(dist_symnorm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6f7413c9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1., 1., ..., 1., 1., 1.],\n",
       "       [1., 0., 1., ..., 1., 1., 1.],\n",
       "       [1., 1., 0., ..., 1., 1., 1.],\n",
       "       ...,\n",
       "       [1., 1., 1., ..., 0., 1., 1.],\n",
       "       [1., 1., 1., ..., 1., 0., 1.],\n",
       "       [1., 1., 1., ..., 1., 1., 0.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unweighted_Adj = np.zeros(distance_matrix.shape) # initialize A w/ zeros\n",
    "\n",
    "for i in range(unweighted_Adj.shape[0]):\n",
    "    for j in range(unweighted_Adj.shape[1]):\n",
    "        if dist_symnorm[i,j] > 0:\n",
    "            unweighted_Adj[i,j] = 1\n",
    "unweighted_Adj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4ca7c8b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00000000e+00, 2.59327195e-09, 1.71357319e-07, ...,\n",
       "        2.17493361e-09, 3.45551945e-09, 1.32782317e-09],\n",
       "       [2.59327195e-09, 0.00000000e+00, 3.17040591e-09, ...,\n",
       "        6.54284880e-09, 8.77257985e-09, 4.24543393e-09],\n",
       "       [1.71357319e-07, 3.17040591e-09, 0.00000000e+00, ...,\n",
       "        2.31161124e-09, 4.68993365e-09, 1.40216127e-09],\n",
       "       ...,\n",
       "       [2.17493361e-09, 6.54284880e-09, 2.31161124e-09, ...,\n",
       "        0.00000000e+00, 2.12410191e-09, 2.74234308e-08],\n",
       "       [3.45551945e-09, 8.77257985e-09, 4.68993365e-09, ...,\n",
       "        2.12410191e-09, 0.00000000e+00, 1.51013781e-09],\n",
       "       [1.32782317e-09, 4.24543393e-09, 1.40216127e-09, ...,\n",
       "        2.74234308e-08, 1.51013781e-09, 0.00000000e+00]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_weighted = np.zeros(distance_matrix.shape) # initialize A with zeros\n",
    "l = 0.5\n",
    "\n",
    "for i in range(A_weighted.shape[0]):\n",
    "    for j in range(A_weighted.shape[1]):\n",
    "        if i!=j:\n",
    "            A_weighted[i,j] = unweighted_Adj[i,j] / (( (distance_matrix[i,j]**2) / (2*l**2) ) ) \n",
    "A_weighted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fe0704b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# distance_matrix_threshold_I_N = np.float32(whole_distance_matrix_threshold_I) ## do not normalize adjcent matrix\n",
    "distance_matrix_threshold_I_N = np.float32(A_weighted)\n",
    "distance_matrix_threshold_I_N_crs = sp.csr_matrix(distance_matrix_threshold_I_N)\n",
    "with open(path+'symmnorm_weighted_adj_mat_crs.pickle', 'wb') as fp:\n",
    "    pickle.dump(distance_matrix_threshold_I_N_crs, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d1b827eb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.0000000e+00, 2.5932720e-09, 1.7135731e-07, ..., 2.1749336e-09,\n",
       "        3.4555194e-09, 1.3278232e-09],\n",
       "       [2.5932720e-09, 0.0000000e+00, 3.1704059e-09, ..., 6.5428489e-09,\n",
       "        8.7725800e-09, 4.2454338e-09],\n",
       "       [1.7135731e-07, 3.1704059e-09, 0.0000000e+00, ..., 2.3116113e-09,\n",
       "        4.6899338e-09, 1.4021613e-09],\n",
       "       ...,\n",
       "       [2.1749336e-09, 6.5428489e-09, 2.3116113e-09, ..., 0.0000000e+00,\n",
       "        2.1241020e-09, 2.7423431e-08],\n",
       "       [3.4555194e-09, 8.7725800e-09, 4.6899338e-09, ..., 2.1241020e-09,\n",
       "        0.0000000e+00, 1.5101378e-09],\n",
       "       [1.3278232e-09, 4.2454338e-09, 1.4021613e-09, ..., 2.7423431e-08,\n",
       "        1.5101378e-09, 0.0000000e+00]], dtype=float32)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distance_matrix_threshold_I_N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "ab7aabb8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1337.4694\n",
      "110.0\n",
      "3.4643227e-10\n",
      "5796586.5\n"
     ]
    }
   ],
   "source": [
    "# explore some distance metrics\n",
    "print(max(map(max, distance_matrix_threshold_I_N)))\n",
    "print(len(distance_matrix_threshold_I_N[np.where(distance_matrix_threshold_I_N > 0.0001)])/2)\n",
    "print(min(distance_matrix_threshold_I_N[np.where(distance_matrix_threshold_I_N != 0)]))\n",
    "print(len(distance_matrix_threshold_I_N[np.where(distance_matrix_threshold_I_N < 0.0001)])/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41f661a4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e337b1e4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b38f5460",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a6f55dc4",
   "metadata": {},
   "source": [
    "## Alternative calculation of A using threshold "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a8375cc5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 1138 0.1671071953010279\n",
      "140 1812 0.26607929515418505\n",
      "180 2948 0.4328928046989721\n",
      "220 4456 0.6543318649045521\n",
      "260 6594 0.9682819383259912\n",
      "300 8862 1.3013215859030838\n",
      "340 11546 1.6954478707782672\n",
      "380 14618 2.14654919236417\n"
     ]
    }
   ],
   "source": [
    "###try different distance threshold, so that on average, each cell has x neighbor cells, see Tab. S1 for results\n",
    "distance_array = np.load(path+'dist_array_eucl_xyz.npy')\n",
    "\n",
    "for threshold in range(100,400,40): #[100,140,180,210,220,260]:#range (210,211):#(100,400,40):\n",
    "    num_big = np.where(distance_array<threshold)[0].shape[0]\n",
    "    print (threshold,num_big,str(num_big/(len(cell_st)*2)))\n",
    "    distance_matrix_threshold_I_list = []\n",
    "#     distance_matrix = euclidean_distances(cell_spot[['X','Y']], cell_spot[['X','Y']])\n",
    "    distance_matrix = euclidean_distances(cell_spot.iloc[:, 1:], cell_spot.iloc[:, 1:])\n",
    "    distance_matrix_threshold_I = np.zeros(distance_matrix.shape)\n",
    "    for i in range(distance_matrix_threshold_I.shape[0]):\n",
    "        for j in range(distance_matrix_threshold_I.shape[1]):\n",
    "            if distance_matrix[i,j] <= threshold and distance_matrix[i,j] > 0:\n",
    "                distance_matrix_threshold_I[i,j] = 1\n",
    "    distance_matrix_threshold_I_list.append(distance_matrix_threshold_I)\n",
    "\n",
    "    ############### get normalized sparse adjacent matrix\n",
    "    distance_matrix_threshold_I_N = normalized_adjacency(distance_matrix_threshold_I, symmetric=False)\n",
    "    # distance_matrix_threshold_I_N = np.float32(whole_distance_matrix_threshold_I) ## do not normalize adjcent matrix\n",
    "    distance_matrix_threshold_I_N = np.float32(distance_matrix_threshold_I)\n",
    "    distance_matrix_threshold_I_N_crs = sp.csr_matrix(distance_matrix_threshold_I_N)\n",
    "    with open(path+'adj_mat_crs_'+str(threshold)+'.pickle', 'wb') as fp:\n",
    "        pickle.dump(distance_matrix_threshold_I_N_crs, fp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14e505cd",
   "metadata": {},
   "source": [
    "Import saved Adjacency matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0f455790",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 adj_mat_crs_100.pickle\n",
      "1 adj_mat_crs_140.pickle\n",
      "2 adj_mat_crs_180.pickle\n",
      "3 adj_mat_crs_220.pickle\n",
      "4 adj_mat_crs_260.pickle\n",
      "5 adj_mat_crs_300.pickle\n",
      "6 adj_mat_crs_340.pickle\n",
      "7 adj_mat_crs_380.pickle\n",
      "8 symmnorm_weighted_adj_mat_crs.pickle\n",
      "average cell neighborhood for file 0 : 0.1486049926578561\n",
      "average cell neighborhood for file 1 : 0.3465491923641703\n",
      "average cell neighborhood for file 2 : 0.6801762114537445\n",
      "average cell neighborhood for file 3 : 1.1230543318649044\n",
      "average cell neighborhood for file 4 : 1.7509544787077826\n",
      "average cell neighborhood for file 5 : 2.417033773861968\n",
      "average cell neighborhood for file 6 : 3.2052863436123347\n",
      "average cell neighborhood for file 7 : 4.107488986784141\n",
      "average cell neighborhood for file 8 : 3404.0\n"
     ]
    }
   ],
   "source": [
    "import os \n",
    "\n",
    "A = []\n",
    "# read in all pickle files\n",
    "for i, file in enumerate([path for path in os.listdir(path) if path.endswith('pickle')]):\n",
    "    print(i, file)\n",
    "    for i in [path+file]:\n",
    "        with open(i, 'rb') as filepath:\n",
    "            while True:\n",
    "                try:\n",
    "                    A.append(pickle.load(filepath))\n",
    "                except EOFError:\n",
    "                    break    \n",
    "                    \n",
    "# uncompress Compressed Sparse Row matrix\n",
    "adj_mat = []\n",
    "for i in range(len(A)):\n",
    "    adj_mat.append(A[i].todense())\n",
    "\n",
    "# get number of counts for all connected cells of cell i\n",
    "sum_neighb_cell = []\n",
    "for file in range(len(adj_mat)):\n",
    "     sum_neighb_cell.append([(adj_mat[file][i][np.where(adj_mat[file][i] != 0)].size) for i in range(len(adj_mat[file]))])\n",
    "\n",
    "# get average cell neighbourhood by threshold\n",
    "for i in range(len(adj_mat)):\n",
    "    avg_cell_neib = sum(sum_neighb_cell[i])/len(adj_mat[0])\n",
    "    print('average cell neighborhood for file', i,':', avg_cell_neib)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a9efda3",
   "metadata": {},
   "source": [
    "GNCG [supp](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02214-w#Sec16)\n",
    "\n",
    "seqFISH+ avg cell neighbourhood for different distance threshold:\n",
    "- 100: 1.2\n",
    "- 140: 2.5 <-- threshold with the best validation performace\n",
    "- 180: 4.0\n",
    "- 220: 5.8\n",
    "\n",
    "MERFISH:\n",
    "- 130: 1.03 <-- threshold with best validation performance\n",
    "- 150: 1.52\n",
    "- 170: 2.05\n",
    "- 200: 2.95\n",
    "- 240: 4.25\n",
    "\n",
    "Try thresholds:\n",
    "- 220: 1.12\n",
    "- 300: 2.42"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8d0f8c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "A300 = adj_mat[5]\n",
    "A220 = adj_mat[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5ad99940",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8230\n",
      "3824\n"
     ]
    }
   ],
   "source": [
    "# how many connections (ie. count number of non-zeros)\n",
    "print(A220[np.where(A220 > 0)].size)\n",
    "print(A300[np.where(A300 > 0)].size) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b642308",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}